Commit df123c7a authored by Gregor Moenke's avatar Gregor Moenke

fixed thresholded ridge plotting bug

parent 7cf15404
### pyBOAT 0.8.23
- Fixed readout plotting bug for thresholded ridges
- added confidence regions from empirical background spectra to API
### pyBOAT 0.8.22
- Fixed export paths for Windows platforms
......
......@@ -3,7 +3,7 @@
import sys,os
import argparse
__version__ = '0.8.22'
__version__ = '0.8.23'
# the object oriented API
from .api import WAnalyzer
......
......@@ -203,24 +203,26 @@ class WAnalyzer:
if smoothing_wsize and smoothing_wsize % 2 == 0:
smoothing_wsize = smoothing_wsize + 1
ridge_y = core.get_maxRidge_ys(modulus)
ridge_power = modulus[ridge_y, np.arange(Nt)]
ridge_ys, ridge_powers = core.get_maxRidge_ys(modulus)
rd = core.eval_ridge(
ridge_y,
ridge_ys,
self.wlet,
self.ana_signal,
sel.periods,
self.periods,
tvec=tvec,
smoothing_wsize=smoothing_wsize,
)
# now threshold
inds = (
ridge_power > power_thresh
) # boolean array of positions exceeding power threshold
rd = rd[inds]
# boolean array of positions exceeding power threshold
b_inds = ridge_powers > power_thresh
# threshold the ridge data
rd = rd[b_inds]
# hitchhike a bit more data
rd.dt = self.dt
rd.Nt = Nt
self.ridge_data = rd
self._has_ridge = True
......@@ -243,6 +245,9 @@ class WAnalyzer:
""" not implemented yet """
print("Not properly implemented yet..")
return
if not self._has_spec:
print("Need to compute a wavelet spectrum first!")
return
......
......@@ -108,13 +108,16 @@ def get_maxRidge_ys(modulus):
Returns
-------
ridge_y : the y-coordinates of the ridge
ridge_ys : the y-coordinates of the ridge
ridge_powers : the power values along the ridge
"""
ridge_y = np.argmax(modulus, axis=0)
ridge_ys = np.argmax(modulus, axis=0)
Nt = modulus.shape[1]
ridge_powers = modulus[ridge_ys, np.arange(Nt)]
return ridge_y
return ridge_ys, ridge_powers
def eval_ridge(
......@@ -279,19 +282,20 @@ def find_COI_crossing(rd):
# first time point outside left COI
coi_inds = left_inds[coi[left_inds] > rd.periods[left_inds]]
# left ridge might be entirely outside COI
i_left = coi_inds[0] if coi_inds.size > 0 else 0
i_left = coi_inds[0] if coi_inds.size > 0 else rd.index[0]
right_inds = np.intersect1d(N2 + np.arange(N2), rd.index)
# last time point outside right COI
coi_inds = right_inds[coi[right_inds] > rd.periods[right_inds]]
# right ridge might be entirely outside COI
i_right = coi_inds[-1] if coi_inds.size > 0 else -1
i_right = coi_inds[-1] if coi_inds.size > 0 else rd.index[-1]
return i_left, i_right
# ============ Snake Annealing =====================================
# ============ Snake Annealing - too slow, needs C-extension ===================
def find_ridge_anneal(landscape, y0, T_ini, Nsteps, mx_jump=2, curve_pen=0):
......
......@@ -335,9 +335,9 @@ def plot_readout(ridge_data, time_unit="a.u.", draw_coi = False, fig=None):
powers = ridge_data["power"]
amplitudes = ridge_data["amplitude"]
ridge_t = ridge_data["time"] # indexable
# check for discontinuous ridge
if np.all(np.diff(ridge_t, n = 2) < 1e-12):
if np.all(np.diff(ridge_t) == ridge_data.dt):
# draw continuous lines
lstyle = '-'
mstyle = ''
......
......@@ -123,6 +123,7 @@ class WaveletAnalyzer(QMainWindow):
self.signal_id = signal_id
self.signal = signal
self.p_max = p_max
self.dt = dt
self.time_unit = time_unit
self.periods = np.linspace(T_min, T_max, step_num)
......@@ -141,8 +142,8 @@ class WaveletAnalyzer(QMainWindow):
self.anneal_pars = None
# =============Compute Spectrum================================================
self.modulus, self.wlet = core.compute_spectrum(self.signal, dt, self.periods)
# =============================================================================
self.modulus, self.wlet = core.compute_spectrum(self.signal, self.dt, self.periods)
# ============================================================================
# Wavelet ridge-readout results
self.ResultWindows = {}
......@@ -369,10 +370,9 @@ class WaveletAnalyzer(QMainWindow):
def do_maxRidge_detection(self):
ridge_ys = core.get_maxRidge_ys(self.modulus)
self.ridge = ridge_y
ridge_ys, ridge_powers = core.get_maxRidge_ys(self.modulus)
if not np.any(ridge_y):
if not np.any(ridge_ys):
self.e = MessageWindow(
"No ridge found..check spectrum!", "Ridge detection error"
)
......@@ -386,6 +386,14 @@ class WaveletAnalyzer(QMainWindow):
self.tvec,
smoothing_wsize=self.rsmoothing,
)
# boolean array of positions exceeding power threshold
b_inds = ridge_powers > self.power_thresh
# threshold the ridge data
ridge_data = ridge_data[b_inds]
# hitchhike a bit more data
ridge_data.dt = self.dt
ridge_data.Nt = len(self.tvec)
self._has_ridge = True
self.ridge_data = ridge_data
......
......@@ -456,19 +456,26 @@ class BatchProcessWindow(QMainWindow):
# compute the spectrum
modulus, wlet = pyboat.compute_spectrum(signal, self.parentDV.dt, periods )
# get maximum ridge
ridge = pyboat.get_maxRidge_ys(modulus)
ridge_ys, ridge_powers = pyboat.get_maxRidge_ys(modulus)
# generate time vector
tvec = np.arange(len(signal)) * self.parentDV.dt
# evaluate along the ridge
ridge_data = pyboat.eval_ridge(
ridge,
ridge_ys,
wlet,
signal,
periods,
tvec,
power_thresh,
smoothing_wsize = rsmooth)
# boolean array of positions exceeding power threshold
b_inds = ridge_powers > power_thresh
# threshold the ridge data
ridge_data = ridge_data[b_inds]
# hitchhike a bit more data
ridge_data.dt = self.parentDV.dt
ridge_data.Nt = len(tvec)
# from ridge thresholding..
if ridge_data.empty:
EmptyRidge += 1
......
......@@ -55,7 +55,7 @@ ppl.tight_layout()
modulus, wlet = pyboat.compute_spectrum(detr_signal, dt, periods)
# get maximum ridge
ridge_ys = pyboat.get_maxRidge_ys(modulus)
ridge_ys, ridge_powers = pyboat.get_maxRidge_ys(modulus)
# evaluate along the ridge, ridge_results is a pandas DataFrame
ridge_results = pyboat.eval_ridge(ridge_ys, wlet, signal, periods, tvec)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment