CSchoel / nolds

Nonlinear measures for dynamical systems (based on one-dimensional time series)
MIT License
263 stars 59 forks source link

Error computing lyap_r #53

Open lcoandrade opened 1 month ago

lcoandrade commented 1 month ago

I have a time series of meteorological data. I'm trying to compute lyap_r with this code:

lyap_r = nolds.lyap_r(self.data[self.variable])

But I get this error:

Traceback (most recent call last):
  File "/Users/luiz/Codes/meteorological_analysis/process/process.py", line 23, in multi_run
    processor.main()
  File "/Users/luiz/Codes/meteorological_analysis/process/process.py", line 387, in main
    lyap_r = self.compute_Lyapunov()
  File "/Users/luiz/Codes/meteorological_analysis/process/process.py", line 349, in compute_Lyapunov
    lyap_r = nolds.lyap_r(self.data[self.variable])
  File "/Users/luiz/venvs/geoproc/lib/python3.9/site-packages/nolds/measures.py", line 249, in lyap_r
    min_tsep = int(np.ceil(1.0 / mf))
TypeError: ufunc 'ceil' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

This is the series I'm using:

DATE;T2M_MAX;T2M_MIN;PRECTOTCORR;HR;PA;WS2M;WD2M
5/6/1984;11,76;5,36;4,22;91,19;69,72;4,22;93,44
6/6/1984;10,58;4,79;6,02;92,12;69,8;4,1;112,88
7/6/1984;10,99;4,28;3,31;89,5;69,86;3,38;117,81
8/6/1984;13,03;4,08;3,3;89,75;69,78;3,39;111,06
9/6/1984;14,04;5,24;1,27;89;69,67;3,73;83
10/6/1984;12,55;6,05;1,82;92,12;69,7;4,95;96,12
11/6/1984;11,61;5,3;3,59;91,38;69,71;4,95;104,31
12/6/1984;11,98;4,96;3,92;90,38;69,72;5,15;105,19
13/6/1984;10,34;5,65;2,32;92,5;69,7;4,96;106
14/6/1984;12,69;4,78;0,48;90,81;69,71;4,69;110,69
15/6/1984;10,28;5,39;2,21;92,19;69,76;4,66;101,88
16/6/1984;11,96;5,45;1,15;91,94;69,77;5,28;99,44
17/6/1984;11,88;5,33;1,58;90,81;69,72;5,25;98,31
18/6/1984;10,87;4,91;0,49;91,31;69,72;4,65;105,94
19/6/1984;11,21;5,09;3,42;91,5;69,77;5,27;109,69
20/6/1984;11,05;5,19;1,2;91,5;69,8;4,91;103,69
......

I'm using numpy==1.26.4 and nolds==0.6.

How can I solve this problem?

CSchoel commented 1 month ago

Not sure what is happening here. Since a TypeError is reported, I first thought your input data might have a weird, type, but the first line in lyap_r is data = np.asarray(data, dtype=np.float64), so everything should be a float after that anyway. :thinking:

A few things you might want to try:

lcoandrade commented 1 month ago

Using min_tsep=1 it works.... But I would need a good approach for min_tsep, that would be the mean period, but how to make it work?

CSchoel commented 1 month ago

Can you just manually execute the following code from lyap_r in a script or notebook and tell me what mf is? (Both print(mf) and print(type(mf)).)

data = np.asarray(data, dtype=np.float64)
n = len(data)
max_tsep_factor = 0.25
f = np.fft.rfft(data, n * 2 - 1)
mf = np.fft.rfftfreq(n * 2 - 1) * f**2
mf = np.sum(mf[1:]) / np.sum(f[1:]**2)

Alternatively, you can also just send me a sample dataset and the code that you use to load it and I can try to debug this on my end. This week, I actually have a bit of time to work on nolds, so I would be interested to investigate this. If this happens with sane input data, it might be an actual bug in the code for calculating min_tsep.

lcoandrade commented 1 month ago

Thanks for your help! Let me give one of the datasets, the respective yaml and the code I'm using.

This is one of the datasets: el_carmen.csv

This is the yaml file for the data: el_carmen.yaml.zip

This is the code: el_carmen.yaml.zip

The prints you asked me are these:

MF:  (0.0001414471551393651+1.1866513857669307e-05j)
MF Type: <class 'numpy.complex128'>
CSchoel commented 1 month ago

Nice, thanks. Now it all makes a lot more sense. mf somehow ends up as a complex number (albeit not very complex given the e-5) and the ceiling of a complex number is not defined. I was under the impression that you couldn't get complex numbers if all your inputs are real numbers, but it seems I made a mistake there. :thinking: For now, you can take np.ceil(1/mf.real) as an estimate for min_tsep. That would end up at 7070, which is quite high. Nolds would reduce that down to about half to ensure that min_tsep is at maximum 25% of the number of datapoints in your input, but you might want to experiment with different values and see how the debug plot responds to that. I don't think it has to be that high.

lcoandrade commented 1 month ago

Would you agree with this?

mtsep = np.ceil(1/mf.real)
mtsep = min(mtsep, int(0.25 * n))
pawel-biernacki-roche commented 1 week ago

We are observing something similar in our code with pandas dataframe of floats. Before upgrade from 0.5.2 to 0.6.x this was not a problem.

Can we hope for a proper fix in 0.6.2 soon?