From 976f68f8402150c57a9aee4c732f3e354c36dfcd Mon Sep 17 00:00:00 2001 From: remnrem Date: Sun, 17 Nov 2024 12:39:07 -0500 Subject: [PATCH] h cut --- edf/edf.cpp | 46 +++------------- spindles/spindles.cpp | 30 ----------- timeline/hypno.cpp | 120 +++++++++++++++++++++++++++++------------- 3 files changed, 89 insertions(+), 107 deletions(-) diff --git a/edf/edf.cpp b/edf/edf.cpp index eb11a82..1a6b8f9 100644 --- a/edf/edf.cpp +++ b/edf/edf.cpp @@ -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? } @@ -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 tok2; if ( alias != "" ) @@ -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 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() ) @@ -4217,9 +4197,6 @@ void edf_t::update_records( int a , int b , int s , const std::vector * 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; } @@ -4724,10 +4701,6 @@ int edf_t::add_time_track( const std::vector * 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(); } @@ -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; @@ -4851,10 +4820,7 @@ void edf_t::reverse( const int s ) const int np = d->size(); std::vector reversed( np ); for (int i=0;i mean_adj_vals( averaged.size() ); - // std::vector median_adj_vals( averaged.size() ); - // for (int i=0;i 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::const_iterator tt = tvals.begin(); while ( tt != tvals.end() ) { @@ -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 threshold( sz , multiplicative_threshold * mean ); @@ -2103,23 +2090,6 @@ annot_t * spindle_wavelet( edf_t & edf , param_t & param ) std::vector so_angle( nbins ); - // // ----------------------------------dumm - // std::vector dumm_x( p_chirp_if->size() , 22 ); - - // std::vector dumm_n ; - // std::vector dumm_b( in_spindle.size() , true ); - // std::vector dumm_dumm = p_sw->phase_locked_averaging( &dumm_x , nbins , - // &dumm_b , &dumm_n ); - // if ( dumm_dumm.size() == nbins ) - // { - // for (int j=0;j pl_chirp_n ; std::vector pl_chirp = p_sw->phase_locked_averaging( p_chirp_if , nbins , &in_spindle , &pl_chirp_n ); diff --git a/timeline/hypno.cpp b/timeline/hypno.cpp index 1e11e31..2959858 100644 --- a/timeline/hypno.cpp +++ b/timeline/hypno.cpp @@ -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 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 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 msp( ne_gaps , false ); for (int e=0; eW transition? @@ -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 msp0 = msp; - // add left flankers? + // add left flankers? (based on eflnk) for (int e=1; eS transition? @@ -832,54 +841,91 @@ void hypnogram_t::edit( timeline_t * timeline , param_t & param ) // add right flankers? for (int e=0; eW 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 clead( ne_gaps , 0 ); - int cntS = 0 , cntW = 0; + // forward scoring + std::vector scr( ne_gaps , 0 ); + int s = 0; for (int e=0; e ctrail( ne_gaps , 0 ); - cntS = 0; cntW = 0; + // reverse scoring + std::vector 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 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 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