From 75e8ae38f0302f4f1886478e07e72eeb8bcae297 Mon Sep 17 00:00:00 2001 From: iguinn Date: Thu, 21 Mar 2024 08:24:37 -0700 Subject: [PATCH 1/5] Added multi_time_point_thresh --- src/dspeed/processors/__init__.py | 3 +- src/dspeed/processors/time_point_thresh.py | 157 +++++++++++++++++++++ 2 files changed, 159 insertions(+), 1 deletion(-) diff --git a/src/dspeed/processors/__init__.py b/src/dspeed/processors/__init__.py index 12bc5d9..a0e8b12 100644 --- a/src/dspeed/processors/__init__.py +++ b/src/dspeed/processors/__init__.py @@ -92,7 +92,7 @@ from .soft_pileup_corr import soft_pileup_corr, soft_pileup_corr_bl from .svm import svm_predict from .time_over_threshold import time_over_threshold -from .time_point_thresh import interpolated_time_point_thresh, time_point_thresh +from .time_point_thresh import interpolated_time_point_thresh, time_point_thresh, multi_time_point_thresh from .transfer_function_convolver import transfer_function_convolver from .trap_filters import asym_trap_filter, trap_filter, trap_norm, trap_pickoff from .upsampler import interpolating_upsampler, upsampler @@ -143,6 +143,7 @@ "svm_predict", "time_point_thresh", "interpolated_time_point_thresh", + "multi_time_point_thresh", "asym_trap_filter", "trap_filter", "trap_norm", diff --git a/src/dspeed/processors/time_point_thresh.py b/src/dspeed/processors/time_point_thresh.py index f5de8c1..5e303bb 100644 --- a/src/dspeed/processors/time_point_thresh.py +++ b/src/dspeed/processors/time_point_thresh.py @@ -181,3 +181,160 @@ def interpolated_time_point_thresh( ) else: raise DSPFatal("Unrecognized interpolation mode") + + +@guvectorize( + [ + "void(float32[:], float32[:], float32, float32, char, float32[:])", + "void(float64[:], float64[:], float64, float64, char, float64[:])", + ], + "(n),(m),(),(),()->(m)", + **nb_kwargs, +) +def multi_time_point_thresh( + w_in: np.ndarray, + a_threshold: np.ndarray, + t_start: int, + polarity: int, + mode_in: np.int8, + t_out: np.ndarray, +) -> None: + """Find the time where the waveform value crosses the threshold + + Search performed walking either forward or backward from the starting + index. Use interpolation to estimate a time between samples. Interpolation + mode selected with `mode_in`. + + Parameters + ---------- + w_in + the input waveform. + a_threshold + list of threshold values. + t_start + the starting index. + polarity + is the average slope of the waveform up (>0) or down (<0) in the + search region; only sign matters, not value. Raise Exception if 0. + mode_in + Character selecting which interpolation method to use. Note this + must be passed as a ``int8``, e.g. ``ord('i')``. Options: + + * ``i`` -- integer `t_in`; equivalent to + :func:`~.dsp.processors.fixed_sample_pickoff` + * ``b`` -- before; closest integer sample before threshold crossing + * ``a`` -- after; closest integer sample after threshold crossing + * ``r`` -- round; round to nearest integer sample to threshold crossing + * ``l`` -- linear interpolation + The following modes are meant to mirror the options + dspeed.upsampler + * ``f`` -- floor; interpolated values are at previous neighbor. + Equivalent to ``a`` + * ``c`` -- ceiling, interpolated values are at next neighbor. + Equivalent to ``b`` + * ``n`` -- nearest-neighbor interpolation; threshold crossing is + half-way between samples + * ``h`` -- Hermite cubic spline interpolation (*not implemented*) + * ``s`` -- natural cubic spline interpolation (*not implemented*) + t_out + the index where the waveform value crosses the threshold. + + JSON Configuration Example + -------------------------- + + .. code-block :: json + + "tp_0": { + "function": "time_point_thresh", + "module": "dspeed.processors", + "args": ["wf_atrap", "bl_std", "tp_start", 0, "'l'", "tp_0"], + "unit": "ns" + } + """ + t_out[:] = np.nan + + if ( + np.isnan(w_in).any() + or np.isnan(a_threshold).any() + or np.isnan(t_start) + ): + return + + if t_start < 0 or t_start >= len(w_in): + return + + # make polarity +/- 1 + if polarity>0: polarity = 1 + elif polarity<0: polarity = -1 + else: raise DSPFatal("polarity cannot be 0") + + sorted_idx = np.argsort(a_threshold) + + # Get initial values for search + t_start = int(t_start) + a_start = w_in[t_start] + i_start = len(sorted_idx) + for i in range(len(sorted_idx)): + if a_threshold[sorted_idx[i]] >= a_start: + i_start = i + break + + # Search for timepoints at larger values + i_tp = i_start + if i_tp0 else -1, polarity): + if i_tp >= len(sorted_idx): break + while w_in[i_wf] <= a_threshold[idx] < w_in[i_wf + polarity]: + if mode_in == ord("i"): # return index closest to start of search + t_out[idx] = i_wf + elif mode_in in (ord("a"), ord("f")): # return index after crossing + t_out[idx] = i_wf if polarity<0 else i_wf + 1 + elif mode_in in (ord("b"), ord("c")): # return index before crossing + t_out[idx] = i_wf if polarity>0 else i_wf - 1 + elif mode_in == ord("r"): # round; return closest index + if a_threshold[idx] - w_in[i_wf] < w_in[i_wf+polarity] - a_threshold[sorted_idx[i_tp]]: + t_out[idx] = i_wf + else: + t_out[idx] = i_wf + polarity + elif mode_in == ord("n"): # nearest-neighbor; return half-way between samps + t_out[idx] = i_wf + 0.5*polarity + elif mode_in == ord("l"): # linear + t_out[idx] = i_wf + (a_threshold[idx] - w_in[i_wf]) / ( + w_in[i_wf + polarity] - w_in[i_wf] + ) + else: + raise DSPFatal("Unrecognized interpolation mode") + i_tp += 1 + if i_tp >= len(sorted_idx): break + idx = sorted_idx[i_tp] + + # Search for timepoints at smaller values + i_tp = i_start-1 + if i_tp >= 0: + idx = sorted_idx[i_tp] + for i_wf in range(t_start-1, len(w_in)-1 if polarity<0 else -1, -polarity): + if i_tp < 0: break + while w_in[i_wf] <= a_threshold[idx] < w_in[i_wf+polarity]: + if mode_in == ord("i"): # return index closest to start of search + t_out[idx] = i_wf + elif mode_in in (ord("a"), ord("f")): # return index after crossing + t_out[idx] = i_wf if polarity<0 else i_wf + 1 + elif mode_in in (ord("b"), ord("c")): # return index before crossing + t_out[idx] = i_wf if polarity>0 else i_wf - 1 + elif mode_in == ord("r"): # round; return closest index + if a_threshold[idx] - w_in[i_wf] < w_in[i_wf + polarity] - a_threshold[idx]: + t_out[idx] = i_wf + else: + t_out[idx] = i_wf + polarity + elif mode_in == ord("n"): # nearest-neighbor; return half-way between samps + t_out[idx] = i_wf + 0.5*polarity + elif mode_in == ord("l"): # linear + t_out[idx] = i_wf + (a_threshold[idx] - w_in[i_wf]) / ( + w_in[i_wf + polarity] - w_in[i_wf] + ) + else: + raise DSPFatal("Unrecognized interpolation mode") + i_tp -= 1 + if i_tp < 0: break + idx = sorted_idx[i_tp] From 3b54ae7dfd96ca8d39ec45aa7b6bd4dfdf05f40d Mon Sep 17 00:00:00 2001 From: iguinn Date: Thu, 21 Mar 2024 09:30:52 -0700 Subject: [PATCH 2/5] Fixed bug in timepoint threshold that caused it to stop searching one sample early --- src/dspeed/processors/time_point_thresh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dspeed/processors/time_point_thresh.py b/src/dspeed/processors/time_point_thresh.py index 5e303bb..23992fc 100644 --- a/src/dspeed/processors/time_point_thresh.py +++ b/src/dspeed/processors/time_point_thresh.py @@ -71,7 +71,7 @@ def time_point_thresh( t_out[0] = i return else: - for i in range(int(t_start), 1, -1): + for i in range(int(t_start), 0, -1): if w_in[i - 1] < a_threshold <= w_in[i]: t_out[0] = i return From 6538d61acecc53f9c67a4e6171820b344ad317a0 Mon Sep 17 00:00:00 2001 From: iguinn Date: Thu, 21 Mar 2024 09:31:35 -0700 Subject: [PATCH 3/5] Added some more intuitive interpolation options to interpolated_time_point_thresh --- src/dspeed/processors/time_point_thresh.py | 24 +++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/src/dspeed/processors/time_point_thresh.py b/src/dspeed/processors/time_point_thresh.py index 23992fc..1b20c29 100644 --- a/src/dspeed/processors/time_point_thresh.py +++ b/src/dspeed/processors/time_point_thresh.py @@ -115,13 +115,18 @@ def interpolated_time_point_thresh( * ``i`` -- integer `t_in`; equivalent to :func:`~.dsp.processors.fixed_sample_pickoff` - * ``f`` -- floor; interpolated values are at previous neighbor, so - threshold crossing is at next neighbor - * ``c`` -- ceiling, interpolated values are at next neighbor, so - threshold crossing is at previous neighbor + * ``b`` -- before; closest integer sample before threshold crossing + * ``a`` -- after; closest integer sample after threshold crossing + * ``r`` -- round; round to nearest integer sample to threshold crossing + * ``l`` -- linear interpolation + The following modes are meant to mirror the options + dspeed.upsampler + * ``f`` -- floor; interpolated values are at previous neighbor. + Equivalent to ``a`` + * ``c`` -- ceiling, interpolated values are at next neighbor. + Equivalent to ``b`` * ``n`` -- nearest-neighbor interpolation; threshold crossing is half-way between samples - * ``l`` -- linear interpolation * ``h`` -- Hermite cubic spline interpolation (*not implemented*) * ``s`` -- natural cubic spline interpolation (*not implemented*) t_out @@ -169,10 +174,15 @@ def interpolated_time_point_thresh( if mode_in == ord("i"): # return index before crossing t_out[0] = i_cross - elif mode_in == ord("f"): # return index after crossing + elif mode_in in (ord("a"), ord("f")): # return index after crossing t_out[0] = i_cross + 1 - elif mode_in == ord("c"): # return index before crossing + elif mode_in in (ord("b"), ord("c")): # return index before crossing t_out[0] = i_cross + elif mode_in == ord("r"): # return closest index to crossing + if abs(a_threshold - w_in[i_cross]) < abs(a_threshold - w_in[i_cross+1]): + t_out[0] = i_cross + else: + t_out[0] = i_cross + 1 elif mode_in == ord("n"): # nearest-neighbor; return half-way between samps t_out[0] = i_cross + 0.5 elif mode_in == ord("l"): # linear From 676e0deab1d5a31d3bf9b890b4d0e1e358bcd197 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 21 Mar 2024 16:41:59 +0000 Subject: [PATCH 4/5] style: pre-commit fixes --- src/dspeed/processors/__init__.py | 6 +- src/dspeed/processors/time_point_thresh.py | 75 +++++++++++++--------- 2 files changed, 50 insertions(+), 31 deletions(-) diff --git a/src/dspeed/processors/__init__.py b/src/dspeed/processors/__init__.py index a0e8b12..fde47cf 100644 --- a/src/dspeed/processors/__init__.py +++ b/src/dspeed/processors/__init__.py @@ -92,7 +92,11 @@ from .soft_pileup_corr import soft_pileup_corr, soft_pileup_corr_bl from .svm import svm_predict from .time_over_threshold import time_over_threshold -from .time_point_thresh import interpolated_time_point_thresh, time_point_thresh, multi_time_point_thresh +from .time_point_thresh import ( + interpolated_time_point_thresh, + multi_time_point_thresh, + time_point_thresh, +) from .transfer_function_convolver import transfer_function_convolver from .trap_filters import asym_trap_filter, trap_filter, trap_norm, trap_pickoff from .upsampler import interpolating_upsampler, upsampler diff --git a/src/dspeed/processors/time_point_thresh.py b/src/dspeed/processors/time_point_thresh.py index 1b20c29..081adc5 100644 --- a/src/dspeed/processors/time_point_thresh.py +++ b/src/dspeed/processors/time_point_thresh.py @@ -119,7 +119,7 @@ def interpolated_time_point_thresh( * ``a`` -- after; closest integer sample after threshold crossing * ``r`` -- round; round to nearest integer sample to threshold crossing * ``l`` -- linear interpolation - The following modes are meant to mirror the options + The following modes are meant to mirror the options dspeed.upsampler * ``f`` -- floor; interpolated values are at previous neighbor. Equivalent to ``a`` @@ -179,7 +179,7 @@ def interpolated_time_point_thresh( elif mode_in in (ord("b"), ord("c")): # return index before crossing t_out[0] = i_cross elif mode_in == ord("r"): # return closest index to crossing - if abs(a_threshold - w_in[i_cross]) < abs(a_threshold - w_in[i_cross+1]): + if abs(a_threshold - w_in[i_cross]) < abs(a_threshold - w_in[i_cross + 1]): t_out[0] = i_cross else: t_out[0] = i_cross + 1 @@ -236,7 +236,7 @@ def multi_time_point_thresh( * ``a`` -- after; closest integer sample after threshold crossing * ``r`` -- round; round to nearest integer sample to threshold crossing * ``l`` -- linear interpolation - The following modes are meant to mirror the options + The following modes are meant to mirror the options dspeed.upsampler * ``f`` -- floor; interpolated values are at previous neighbor. Equivalent to ``a`` @@ -263,20 +263,19 @@ def multi_time_point_thresh( """ t_out[:] = np.nan - if ( - np.isnan(w_in).any() - or np.isnan(a_threshold).any() - or np.isnan(t_start) - ): + if np.isnan(w_in).any() or np.isnan(a_threshold).any() or np.isnan(t_start): return if t_start < 0 or t_start >= len(w_in): return # make polarity +/- 1 - if polarity>0: polarity = 1 - elif polarity<0: polarity = -1 - else: raise DSPFatal("polarity cannot be 0") + if polarity > 0: + polarity = 1 + elif polarity < 0: + polarity = -1 + else: + raise DSPFatal("polarity cannot be 0") sorted_idx = np.argsort(a_threshold) @@ -291,24 +290,30 @@ def multi_time_point_thresh( # Search for timepoints at larger values i_tp = i_start - if i_tp0 else -1, polarity): - if i_tp >= len(sorted_idx): break + for i_wf in range(t_start, len(w_in) - 1 if polarity > 0 else -1, polarity): + if i_tp >= len(sorted_idx): + break while w_in[i_wf] <= a_threshold[idx] < w_in[i_wf + polarity]: if mode_in == ord("i"): # return index closest to start of search t_out[idx] = i_wf elif mode_in in (ord("a"), ord("f")): # return index after crossing - t_out[idx] = i_wf if polarity<0 else i_wf + 1 + t_out[idx] = i_wf if polarity < 0 else i_wf + 1 elif mode_in in (ord("b"), ord("c")): # return index before crossing - t_out[idx] = i_wf if polarity>0 else i_wf - 1 + t_out[idx] = i_wf if polarity > 0 else i_wf - 1 elif mode_in == ord("r"): # round; return closest index - if a_threshold[idx] - w_in[i_wf] < w_in[i_wf+polarity] - a_threshold[sorted_idx[i_tp]]: + if ( + a_threshold[idx] - w_in[i_wf] + < w_in[i_wf + polarity] - a_threshold[sorted_idx[i_tp]] + ): t_out[idx] = i_wf else: t_out[idx] = i_wf + polarity - elif mode_in == ord("n"): # nearest-neighbor; return half-way between samps - t_out[idx] = i_wf + 0.5*polarity + elif mode_in == ord( + "n" + ): # nearest-neighbor; return half-way between samps + t_out[idx] = i_wf + 0.5 * polarity elif mode_in == ord("l"): # linear t_out[idx] = i_wf + (a_threshold[idx] - w_in[i_wf]) / ( w_in[i_wf + polarity] - w_in[i_wf] @@ -316,29 +321,38 @@ def multi_time_point_thresh( else: raise DSPFatal("Unrecognized interpolation mode") i_tp += 1 - if i_tp >= len(sorted_idx): break + if i_tp >= len(sorted_idx): + break idx = sorted_idx[i_tp] # Search for timepoints at smaller values - i_tp = i_start-1 + i_tp = i_start - 1 if i_tp >= 0: idx = sorted_idx[i_tp] - for i_wf in range(t_start-1, len(w_in)-1 if polarity<0 else -1, -polarity): - if i_tp < 0: break - while w_in[i_wf] <= a_threshold[idx] < w_in[i_wf+polarity]: + for i_wf in range( + t_start - 1, len(w_in) - 1 if polarity < 0 else -1, -polarity + ): + if i_tp < 0: + break + while w_in[i_wf] <= a_threshold[idx] < w_in[i_wf + polarity]: if mode_in == ord("i"): # return index closest to start of search t_out[idx] = i_wf elif mode_in in (ord("a"), ord("f")): # return index after crossing - t_out[idx] = i_wf if polarity<0 else i_wf + 1 + t_out[idx] = i_wf if polarity < 0 else i_wf + 1 elif mode_in in (ord("b"), ord("c")): # return index before crossing - t_out[idx] = i_wf if polarity>0 else i_wf - 1 + t_out[idx] = i_wf if polarity > 0 else i_wf - 1 elif mode_in == ord("r"): # round; return closest index - if a_threshold[idx] - w_in[i_wf] < w_in[i_wf + polarity] - a_threshold[idx]: + if ( + a_threshold[idx] - w_in[i_wf] + < w_in[i_wf + polarity] - a_threshold[idx] + ): t_out[idx] = i_wf else: t_out[idx] = i_wf + polarity - elif mode_in == ord("n"): # nearest-neighbor; return half-way between samps - t_out[idx] = i_wf + 0.5*polarity + elif mode_in == ord( + "n" + ): # nearest-neighbor; return half-way between samps + t_out[idx] = i_wf + 0.5 * polarity elif mode_in == ord("l"): # linear t_out[idx] = i_wf + (a_threshold[idx] - w_in[i_wf]) / ( w_in[i_wf + polarity] - w_in[i_wf] @@ -346,5 +360,6 @@ def multi_time_point_thresh( else: raise DSPFatal("Unrecognized interpolation mode") i_tp -= 1 - if i_tp < 0: break + if i_tp < 0: + break idx = sorted_idx[i_tp] From 758c7c9d2b7099f8cb6ddcf4471bc60ab35c9bd6 Mon Sep 17 00:00:00 2001 From: iguinn Date: Thu, 21 Mar 2024 10:43:51 -0700 Subject: [PATCH 5/5] Fixed failed docs hopefully --- src/dspeed/processors/time_point_thresh.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/dspeed/processors/time_point_thresh.py b/src/dspeed/processors/time_point_thresh.py index 081adc5..30595ac 100644 --- a/src/dspeed/processors/time_point_thresh.py +++ b/src/dspeed/processors/time_point_thresh.py @@ -119,8 +119,10 @@ def interpolated_time_point_thresh( * ``a`` -- after; closest integer sample after threshold crossing * ``r`` -- round; round to nearest integer sample to threshold crossing * ``l`` -- linear interpolation + The following modes are meant to mirror the options dspeed.upsampler + * ``f`` -- floor; interpolated values are at previous neighbor. Equivalent to ``a`` * ``c`` -- ceiling, interpolated values are at next neighbor. @@ -236,8 +238,10 @@ def multi_time_point_thresh( * ``a`` -- after; closest integer sample after threshold crossing * ``r`` -- round; round to nearest integer sample to threshold crossing * ``l`` -- linear interpolation + The following modes are meant to mirror the options dspeed.upsampler + * ``f`` -- floor; interpolated values are at previous neighbor. Equivalent to ``a`` * ``c`` -- ceiling, interpolated values are at next neighbor.