OverLordGoldDragon / ssqueezepy

Synchrosqueezing, wavelet transforms, and time-frequency analysis in Python
MIT License
647 stars 95 forks source link

Auto-scales failure with short inputs, default configs #90

Open rmallof opened 1 year ago

rmallof commented 1 year ago

Description

Hi. Thanks for the library, is being really handy!

I found a way to crash the icwt function in the current version - 0.6.3

Provided are the time series and code required to reproduce the issue

data Error message: --------------------------------------------------------------------------- ValueError Traceback (most recent call last) Input In [2], in () 76 wavelet = 'morlet' 78 _, cwt, _, scales = ssq_cwt(data, wavelet=wavelet) ---> 80 icwt(cwt, wavelet=wavelet, scales=scales) File /opt/conda/envs/gst/lib/python3.8/site-packages/ssqueezepy/_cwt.py:391, in icwt(Wx, wavelet, scales, nv, one_int, x_len, x_mean, padtype, rpadded, l1_norm) 389 idx = logscale_transition_idx(scales) 390 x = icwt(Wx[:idx], scales=scales[:idx], **kw) --> 391 x += icwt(Wx[idx:], scales=scales[idx:], **kw) 392 return x 393 ########################################################################## 394 395 #### Invert ############################################################## File /opt/conda/envs/gst/lib/python3.8/site-packages/ssqueezepy/_cwt.py:378, in icwt(Wx, wavelet, scales, nv, one_int, x_len, x_mean, padtype, rpadded, l1_norm) 376 wavelet = Wavelet._init_if_not_isinstance(wavelet) 377 # will override `nv` to match `scales`'s --> 378 scales, scaletype, _, nv = process_scales(scales, x_len, wavelet, nv=nv, 379 get_params=True) 380 assert (len(scales) == na), "%s != %s" % (len(scales), na) 382 #### Handle piecewise scales case ######################################## 383 # `nv` must be left unspecified so it's inferred automatically from `scales` 384 # in `process_scales` for each piecewise case File /opt/conda/envs/gst/lib/python3.8/site-packages/ssqueezepy/utils/cwt_utils.py:247, in process_scales(scales, N, wavelet, nv, get_params, use_padded_N) 244 nv = int(nv) 245 return scaletype, nv, preset --> 247 scaletype, nv, preset = _process_args(scales, nv, wavelet) 248 if isinstance(scales, (np.ndarray, torch.Tensor)): 249 scales = scales.reshape(-1, 1) File /opt/conda/envs/gst/lib/python3.8/site-packages/ssqueezepy/utils/cwt_utils.py:228, in process_scales.._process_args(scales, nv, wavelet) 225 if scales.squeeze().ndim != 1: 226 raise ValueError("`scales`, if array, must be 1D " 227 "(got shape %s)" % str(scales.shape)) --> 228 scaletype, _nv = infer_scaletype(scales) 229 if scaletype == 'log': 230 if nv is not None and _nv != nv: File /opt/conda/envs/gst/lib/python3.8/site-packages/ssqueezepy/utils/cwt_utils.py:289, in infer_scaletype(scales) 286 # ceil to avoid faulty float-int roundoffs 287 nv = int(np.round(1 / np.diff(np.log2(scales), axis=0)[0])) --> 289 elif logscale_transition_idx(scales) is None: 290 raise ValueError("could not infer `scaletype` from `scales`; " 291 "`scales` array must be linear or exponential. " 292 "(got diff(scales)=%s..." % np.diff(scales, axis=0)[:4]) 294 else: File /opt/conda/envs/gst/lib/python3.8/site-packages/ssqueezepy/utils/cwt_utils.py:380, in logscale_transition_idx(scales) 378 scales = asnumpy(scales) 379 scales_diff2 = np.abs(np.diff(np.log(scales), 2, axis=0)) --> 380 idx = np.argmax(scales_diff2) + 2 381 diff2_max = scales_diff2.max() 382 # every other value must be zero, assert it is so File <__array_function__ internals>:180, in argmax(*args, **kwargs) File /opt/conda/envs/gst/lib/python3.8/site-packages/numpy/core/fromnumeric.py:1216, in argmax(a, axis, out, keepdims) 1129 """ 1130 Returns the indices of the maximum values along an axis. 1131 (...) 1213 (2, 1, 4) 1214 """ 1215 kwds = {'keepdims': keepdims} if keepdims is not np._NoValue else {} -> 1216 return _wrapfunc(a, 'argmax', axis=axis, out=out, **kwds) File /opt/conda/envs/gst/lib/python3.8/site-packages/numpy/core/fromnumeric.py:57, in _wrapfunc(obj, method, *args, **kwds) 54 return _wrapit(obj, method, *args, **kwds) 56 try: ---> 57 return bound(*args, **kwds) 58 except TypeError: 59 # A TypeError occurs if the object does have such a method in its 60 # class, but its signature is not identical to that of NumPy's. This (...) 64 # Call _wrapit from within the except clause to ensure a potential 65 # exception has a traceback chain. 66 return _wrapit(obj, method, *args, **kwds) ValueError: attempt to get argmax of an empty sequence ### Code for reproduction ```python from ssqueezepy import ssq_cwt, icwt import numpy as np data = np.array([ 41135817341, 29083562061, 26175547452, 16588370958, 17264085441, 31947336829, 40770974039, 30242059107, 21692004719, 29867476527, 27246574439, 34163220274, 68204556440, 50913575242, 54912007015, 31183975654, 27132421514, 42009436760, 35329942625, 30818458597, 28970212744, 28574793478, 26188097173, 24957784918, 18372538715, 18027170497, 20965695707, 21381535161, 23552740328, 26267239923, 30767551159, 18100418740, 16390821947, 21594638208, 26715546990, 24598943708, 25814972520, 49899834488, 29641127858, 28688807249, 24150249025, 25810220018, 33042430345, 31158743333, 25905575359, 24302954056, 22927802083, 39974475562, 48765202697, 42932549127, 33631012204, 31421555646, 24021799169, 23565495303, 35574561406, 28624673855, 31758955233, 40212386158, 35887249746, 28148218301, 23553591896, 25849159141, 28389250717, 26288169966, 25120229769, 28881249043, 15978259885, 15886817043, 28575544847, 23555719219, 32837431722, 37127036580, 27265804688, 22987346289, 22994133555, 35123501685, 27753685646, 30931623076, 23747613147, 40509610260, 27595671000, 23102307723, 31666498758, 31878280659, 31962253368, 31028679593, 42326789564, 30116729776, 24366810591, 32637854078, 34483360283, 33225232872, 30182031010, 29123998928, 23613051457, 25245861652, 28813460025, 43403978910, 35239757134, 32194477850, 48469528171, 36913738894, 34493951963, 50212088965, 51091116622, 37872380889, 36389011503, 30123362273, 24957448100, 31254779144, 40177002624, 36791346508, 46363793975, 41135767926, 38896078052, 26149643168, 23359966112, 44148798321, 58571439619, 53071298734, 41037843771, 43975248085, 18719537670, 20765955327, 30484729489, 35887278685, 33223790572, 34711412966, 29227315390, 16437423167, 16837262532, 27425022774, 28711532910, 24950173846, 44219840004, 38452356727, 16192235532, 17988916650, 27472552998, 30580012344, 22425387184, 24493974420, 32459287866, 16104440957, 22128794335, 30202235805, 47761524910, 58895950537, 49625110402, 43994715910, 40369840645, 31486345556, 45668466815, 39819303159, 55552169483, 43228750179, 64072727950, 37846047609, 35082693210, 53510852236, 118992465607, 102905151606, 83202283721, 55871616488, 29717699419, 27209183682, 49630243054, 36599436183, 33925512989, 27868914022, 26862218609, 16106223492, 21313378652, 37429485518, 30726828760, 32958875628, 26129037414, 18678255976, 18000008764, 20443898509, 27743025156, 23581685468, 29523576583, 22895392882, 19539705127, 16217776704, 16824520830, 22209086834, 19889922369, 19675404389, 20496603770, 20328426366, 12706781969, 14122486832, 19617581341, 26634741631, 25534481470, 20964448341, 24031608960, 14463581825, 10924354698, 17221074814, 22722096615, 14882945045, 16441573050, 15329265213, 9744636213, 11656379938, 11886957804, 15748580239, 17005713920, 14472237479, 15929162910, 11239186456, 9244361700, 12097775227, 13903079207, 18421743322, 13692758566, 14413662913, 7714767174, 9768827914, 18624736866, 15808338949, 18372283782, 34971338710, 29225029694, 38967784639, 19298407543, 26792494050, 24999983362, 30005625418, 21152848261, 28799154319, 32442278429, 24746386230, 26518700512, 26405069715, 30685366709, 26357839322, 25383335641, 14712928379, 27423687259, 27205595568, 22837828665, 26683255504, 32066936882, 27083066007, 15639298538, 19564262605, 23825006542, 27187964471, 25371367758, 32572572185, 27078406594, 16356226232, 17821046406, 23918742607, 26792596581, 32483312909, 39316664596, 41358451255, 19625427158, 25555105670, 28987376573, 31252098714, 30199996781, 30476264066, 26811744928, 16100721565, 16644534842, 22660763494, 20535363434, 24662841200, 20386398516, 26062404610, 11166012913, 13488941056]) wavelet = 'morlet' _, cwt, _, scales = ssq_cwt(data, wavelet=wavelet) icwt(cwt, wavelet=wavelet, scales=scales) ```
OverLordGoldDragon commented 1 year ago

Thanks for the report, you've reproduced an edge case with short inputs.

Until it's fixed, try one of: 1) higher nv (e.g. 40); 2) scales='log'; 3) longer data.

rmallof commented 1 year ago

Thank you for the response and the suggestions. It does work with longer data indeed. In this case, setting a higher nv parameter or the scales='log' still gives errors. Adding them here for convenience:

icwt(cwt, wavelet=wavelet, nv=40)

---> 80 icwt(cwt, wavelet=wavelet, nv=40)

File /opt/conda/envs/gst/lib/python3.8/site-packages/ssqueezepy/_cwt.py:380, in icwt(Wx, wavelet, scales, nv, one_int, x_len, x_mean, padtype, rpadded, l1norm) 377 # will override nv to match scales's 378 scales, scaletype, , nv = process_scales(scales, x_len, wavelet, nv=nv, 379 get_params=True) --> 380 assert (len(scales) == na), "%s != %s" % (len(scales), na) 382 #### Handle piecewise scales case ######################################## 383 # nv must be left unspecified so it's inferred automatically from scales 384 # in process_scales for each piecewise case 385 if scaletype == 'log-piecewise':

AssertionError: 276 != 222

icwt(cwt, wavelet=wavelet, scales='log') ---> 80 icwt(cwt, wavelet=wavelet, scales='log')

File /opt/conda/envs/gst/lib/python3.8/site-packages/ssqueezepy/_cwt.py:380, in icwt(Wx, wavelet, scales, nv, one_int, x_len, x_mean, padtype, rpadded, l1norm) 377 # will override nv to match scales's 378 scales, scaletype, , nv = process_scales(scales, x_len, wavelet, nv=nv, 379 get_params=True) --> 380 assert (len(scales) == na), "%s != %s" % (len(scales), na) 382 #### Handle piecewise scales case ######################################## 383 # nv must be left unspecified so it's inferred automatically from scales 384 # in process_scales for each piecewise case 385 if scaletype == 'log-piecewise':

AssertionError: 168 != 222

Please, let me know if I'm missing something. Again, great library :)

OverLordGoldDragon commented 1 year ago

Right, I meant in cwt (and therefore also in icwt).