Skip to content

Commit

Permalink
h cut
Browse files Browse the repository at this point in the history
  • Loading branch information
remnrem committed Nov 17, 2024
1 parent da45078 commit 976f68f
Show file tree
Hide file tree
Showing 3 changed files with 89 additions and 107 deletions.
46 changes: 6 additions & 40 deletions edf/edf.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1644,7 +1644,7 @@ bool edf_t::attach( const std::string & f ,
}
else
{
// TODO... need to check EDFZ file. e.g. try reading the last record?
// TODO: check EDFZ file. e.g. try reading the last record?
}


Expand Down Expand Up @@ -2247,12 +2247,6 @@ bool edf_t::is_actually_standard_edf()
// EDF Annotations (other than time track)?
if ( has_edf_annots ) return false;

// for (int s=0;s<header.ns;s++)
// {
// if ( ! header.is_data_channel(s) )
// if ( s != header.t_track ) return false;
// }

// discontinuous?
if ( is_actually_discontinuous() ) return false;

Expand Down Expand Up @@ -3634,11 +3628,9 @@ signal_list_t edf_header_t::signal_list( const std::string & s ,
{
if ( cmd_t::primary_alias.find( cmd_t::primary_upper2orig[ uc_lb ] ) != cmd_t::primary_alias.end() )
{
//std::cout << "tok2_ " << tok2_[t2] << "and alias ["<<alias<< "]\n";
if ( alias == "" ) alias = cmd_t::primary_upper2orig[ uc_lb ] ; // OLD = tok2_[t2];
if ( alias == "" ) alias = cmd_t::primary_upper2orig[ uc_lb ] ;
else if ( ! Helper::iequals( alias , uc_lb ) )
Helper::halt( "more than one alias implied" );
//std::cout << "tok2_ " << tok2_[t2] << "and alias ["<<alias<< "]\n";
Helper::halt( "more than one alias implied" );
}
}
else if ( cmd_t::label_aliases.find( uc_lb ) != cmd_t::label_aliases.end() )
Expand All @@ -3653,8 +3645,6 @@ signal_list_t edf_header_t::signal_list( const std::string & s ,
//
// update list if needed
//
// std::cout << " alias = [" << alias << "]\n";
// std::cout << " primary alias size = " << cmd_t::primary_alias.size() << "\n";

std::vector<std::string> tok2;
if ( alias != "" )
Expand All @@ -3681,32 +3671,22 @@ signal_list_t edf_header_t::signal_list( const std::string & s ,

//
// proceed as before
//

// std::cout << " tok2 = " << tok2.size() << "\n";
//

for (int t2=0;t2<tok2.size();t2++)
{

// std::cout << "t2 = " << t2 << "\t" << tok2[t2] << "\n";
// std::cout << "label2header.find size " << label2header.size() << "\n";

// add first match found
if ( label2header.find( Helper::toupper( tok2[t2] ) ) != label2header.end() )
{

const int l = label2header[ Helper::toupper( tok2[t2] ) ];

// std::cout << "found match " << l << "\n";

if ( t2 > 0 ) // relabel if wasn't first choice?
{
label2header[ Helper::toupper( tok2[0] ) ] = l;
//label[l] = tok2[0];
label2header[ Helper::toupper( tok2[0] ) ] = l;
}

//std::cout << "adding N " << label2header[ Helper::toupper( tok2[0] ) ] << "\n";

const int l0 = label2header[ Helper::toupper( tok2[0] ) ] ;

if ( added.find( l0 ) == added.end() )
Expand Down Expand Up @@ -4217,9 +4197,6 @@ void edf_t::update_records( int a , int b , int s , const std::vector<double> *
if ( x < pmin ) x = pmin;
else if ( x > pmax ) x = pmax;

// std::cout << "edit\t" << edf_record_t::dig2phys( data[p] , bv , os )
// << "\t"
// << (*d)[cnt] << "\n";
data[p] = edf_record_t::phys2dig( (*d)[cnt] , bv , os );
++cnt;
}
Expand Down Expand Up @@ -4724,10 +4701,6 @@ int edf_t::add_time_track( const std::vector<uint64_t> * tps )
r = timeline.next_record(r);
}

// std::cout << "DET1 " << records.begin()->second.pdata.size() << "\t"
// << records.begin()->second.data.size() << "\n";


return header.time_track();

}
Expand Down Expand Up @@ -4803,15 +4776,11 @@ uint64_t edf_t::timepoint_from_EDF( int r )

if ( div10 != 0 )
{
// std::cout << " fixin " << tp << "\t" << div100 << "\t";
// round down, or up
if ( div10 < 5 ) tp -= div10 ;
else tp += 10LLU - div10;
// std::cout << tp << "\n";
}

//std::cout << " READ [ " << tt.substr(0,e) << " ]\t" << tt_sec << "\t" << tp << "\n";

// cache in case recalled
cached_EDF_timepoints[ r ] = tp;

Expand Down Expand Up @@ -4851,10 +4820,7 @@ void edf_t::reverse( const int s )
const int np = d->size();
std::vector<double> reversed( np );
for (int i=0;i<np;i++)
{
reversed[i] = (*d)[np-i-1];
// std::cout << " reversed[i] = " << i << "\t" << reversed[i] << "\n";
}
reversed[i] = (*d)[np-i-1];
update_signal_retain_range( s , &reversed );
}

Expand Down
30 changes: 0 additions & 30 deletions spindles/spindles.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -844,21 +844,10 @@ annot_t * spindle_wavelet( edf_t & edf , param_t & param )
// redundant code, but fine to leave as is for now...
double t_mean = MiscMath::mean( averaged );
double t_median = MiscMath::median( averaged );
// std::vector<double> mean_adj_vals( averaged.size() );
// std::vector<double> median_adj_vals( averaged.size() );
// for (int i=0;i<adj_vals.size();i++) mean_adj_vals[i] = averaged[i] / t_mean;
// for (int i=0;i<adj_vals.size();i++) median_adj_vals[i] = averaged[i] / t_median;

if ( t_median > 0 )
writer.value( "MEAN_OVER_MEDIAN" , t_mean / t_median );

// writer.value( "ADJMEAN_MEAN" , MiscMath::mean( mean_adj_vals ) );
// writer.value( "ADJMEAN_MEDIAN" , MiscMath::median( mean_adj_vals ) );

// writer.value( "ADJMEDIAN_MEAN" , MiscMath::mean( median_adj_vals ) );
// writer.value( "ADJMEDIAN_MEDIAN" , MiscMath::median( median_adj_vals ) );


std::map<double,double>::const_iterator tt = tvals.begin();
while ( tt != tvals.end() )
{
Expand Down Expand Up @@ -892,8 +881,6 @@ annot_t * spindle_wavelet( edf_t & edf , param_t & param )
// e.g. a local sliding window average instead of the whole-night baseline
const int sz = averaged.size();

// std::cout << " averaged size = " << sz << "\n";

// required core to be a spindle
std::vector<double> threshold( sz , multiplicative_threshold * mean );

Expand Down Expand Up @@ -2103,23 +2090,6 @@ annot_t * spindle_wavelet( edf_t & edf , param_t & param )
std::vector<double> so_angle( nbins );


// // ----------------------------------dumm
// std::vector<double> dumm_x( p_chirp_if->size() , 22 );

// std::vector<int> dumm_n ;
// std::vector<bool> dumm_b( in_spindle.size() , true );
// std::vector<double> dumm_dumm = p_sw->phase_locked_averaging( &dumm_x , nbins ,
// &dumm_b , &dumm_n );
// if ( dumm_dumm.size() == nbins )
// {
// for (int j=0;j<nbins;j++)
// {
// std::cout << " DUM PH " << ph << "\t" << dumm_dumm[j] << "\t" << dumm_n[j] << "\n";
// ph += inc;
// }
// }
// //duumy ---------------------------------

std::vector<int> pl_chirp_n ;
std::vector<double> pl_chirp = p_sw->phase_locked_averaging( p_chirp_if , nbins ,
&in_spindle , &pl_chirp_n );
Expand Down
120 changes: 83 additions & 37 deletions timeline/hypno.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -776,25 +776,31 @@ void hypnogram_t::edit( timeline_t * timeline , param_t & param )
// wwwwwwSwwwwwwwwwwwSSwSSwwwwwwwwSSSSSSSSwSSSwSSSSSSSSSSwww
// XXXX XXXXXXXXX XXXXXXXXXXXXXXXXXXXXXXXXXXX
// XXXX XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX (drop W < 1hr between S)
// - calculate cumulate ratio of S/W per epoch moving forward (then back)
// note::: this won't work well w/ GAPPED EDF+D


const bool trim1b = param.has( "trim" );
const bool trim1b = param.has( "cut" );
const bool verbose = param.has( "verbose" );

if ( trim1b )
{
std::vector<double> p = param.dblvector( "trim" );
if ( p.size() != 3 ) Helper::halt( "expecting trim=gap,flank,th" );
const int egap = p[0];
const int eflnk = p[1];
const double th = p[2];

std::vector<double> p;
if ( ! param.empty( "cut" ) ) p = param.dblvector( "cut" );

if ( p.size() > 4 )
Helper::halt( "expecting 0 to 4 params: cut=th,fac,gap,flank" );

const double th = p.size() > 0 ? p[0] : 50; // diff (W vs weighted-S)
const double fac = p.size() > 1 ? p[1] : 3; // S vs W weighting
const int egap = p.size() > 2 ? p[2] : 10; // 5 mins
const int eflnk = p.size() > 3 ? p[3] : 10; // 5 mins

// major sleep period = msp
std::vector<bool> msp( ne_gaps , false );
for (int e=0; e<ne_gaps; e++ )
msp[e] = is_sleep( stages[e] );

// extend blobs (from second epoch)
// extend blobs (from second epoch) based on egap
for (int e=1; e<ne_gaps; e++ )
{
if ( msp[e-1] && ! msp[e] ) // S->W transition?
Expand All @@ -814,8 +820,11 @@ void hypnogram_t::edit( timeline_t * timeline , param_t & param )
e = q;
}
}

// make a copy so we can ignore flankers when getting the prop W/S
std::vector<bool> msp0 = msp;

// add left flankers?
// add left flankers? (based on eflnk)
for (int e=1; e<ne_gaps; e++ )
{
if ( (!msp[e-1] ) && msp[e] ) // W->S transition?
Expand All @@ -832,54 +841,91 @@ void hypnogram_t::edit( timeline_t * timeline , param_t & param )

// add right flankers?
for (int e=0; e<ne_gaps-1; e++ )
{
std::cerr << "testing " << e << "\n";
{
if ( msp[e] && ! msp[e+1] ) // S->W transition?
{
std::cerr << " at trans\n";
{
int q=e;
while (1) {
++q;
if ( q == ne_gaps ) break;
if ( q - e > eflnk ) break;
std::cerr << " right F ing from " << e << " to " << q << "\n";
if ( q - e > eflnk ) break;
msp[q] = true;
}
e = q; // skip ahead
std::cerr << " done, moving onto next...\n";
e = q; // skip ahead
}
}

// get prob: leading
std::vector<double> clead( ne_gaps , 0 );
int cntS = 0 , cntW = 0;
// forward scoring
std::vector<double> scr( ne_gaps , 0 );
int s = 0;
for (int e=0; e<ne_gaps; e++)
{
if ( msp[e] ) cntS++;
else cntW++;
clead[e] = cntS / (double)(cntS + cntW);
if ( msp0[e] ) s -= fac; else s++;
scr[e] = s;
}

// get prob: trailing
std::vector<double> ctrail( ne_gaps , 0 );
cntS = 0; cntW = 0;
// reverse scoring
std::vector<double> rscr( ne_gaps , 0 );
s = 0;
for (int e=ne_gaps-1; e>=0; e--)
{
if ( msp[e] ) cntS++;
else cntW++;
ctrail[e] = cntS / (double)(cntS + cntW);
if ( msp0[e] ) s -= fac; else s++;
rscr[e] = s;
}

// verbose dump
// get cut-points (not allow in msp[])
int cut = -1 , rcut = -1;
double mx = 0;
double rmx = 0;
double mx0 = 0 , rmx0 = 0;

for (int e=0; e<ne_gaps; e++)
{
std::cout << "tb\t" << e+1 << "\t"
<< is_sleep( stages[e] ) << "\t"
<< msp[e] << "\t"
<< clead[e] << "\t"
<< ctrail[e] << "\n";
if ( msp[e] ) scr[e] = 0;
else if ( scr[e] > mx && scr[e] > 0 && scr[e] >= th ) { mx = scr[e] ; cut = e; }
if ( scr[e] > mx0 ) mx0 = scr[e];
}

for (int e=ne_gaps-1; e>0; e--)
{
if ( msp[e] ) rscr[e] = 0;
if ( rscr[e] > rmx && rscr[e] > 0 && rscr[e] >= th ) { rmx = rscr[e] ; rcut = e; }
if ( rscr[e] > rmx0 ) rmx0 = rscr[e];
}

// verbose dump
if ( verbose )
{
for (int e=0; e<ne_gaps; e++)
{
std::cout << "tb\t" << e+1 << "\t"
<< is_sleep( stages[e] ) << "\t"
<< msp[e] << "\t"
<< ( scr[e] > 0 ? scr[e] : 0 ) << "\t"
<< ( rscr[e] > 0 ? rscr[e] : 0 )<< "\n";
}
}

logger << " maximum left-cut and right-scores are " << mx0 << " and " << rmx0 << " (th = " << th << ")\n";

if ( cut != -1 && rcut != -1 )
logger << " placing left-cut at epoch " << cut + 1 << " and right-cut at " << rcut + 1 << "\n";
else if ( cut != -1 )
logger << " placing left-cut at epoch " << cut + 1 << "\n";
else if ( rcut != - 1 )
logger << " placing right-cut at epoch " << rcut + 1 << "\n";
else
logger << " no cuts placed\n";

// check we haven't spliced out everything
if ( cut != 1 && rcut != -1 )
if ( cut >= rcut ) logger << " *** cut isn't leaving any sleep epochs\n";

if ( cut != -1 )
for (int e=0; e<=cut; e++) stages[e] = UNKNOWN;
if ( rcut != -1 )
for (int e=ne_gaps-1; e>= rcut ; e-- ) stages[e] = UNKNOWN;

}

// trimming method 1a
Expand Down

0 comments on commit 976f68f

Please sign in to comment.