EEMD with parabolic option generate shapes error with divide by zero encountered in true_divide

kimyoungjin06 commented 4 years ago

Problem

I try to apply EEMD to a tidal wave.

So, I follow simplest example of EEMD of documentation. The example code is a success, but I get shape error with divide by zero encountered in true_divide in applying with a tidal wave.

Data

An input data is hourly tidal level in 28 days.

CASE 1: without parabolic option

# %%time
# Assign EEMD to `eemd` variable
eemd = EEMD()
# Say we want detect extrema using parabolic method
emd = eemd.EMD
# emd.extrema_detection="parabol"
# Execute EEMD on S
eIMFs = eemd.eemd(_s, _t)
nIMFs = eIMFs.shape[0]
# Plot results
plt.figure(figsize=(12,9))
plt.subplot(nIMFs+1, 1, 1)
plt.plot(_t, _s, 'r')
for n in range(nIMFs):
    plt.subplot(nIMFs+1, 1, n+2)
    plt.plot(_t, eIMFs[n], 'g')
    plt.ylabel("eIMF %i" %(n+1))
    plt.locator_params(axis='y', nbins=5)
plt.xlabel("Time [s]")
plt.tight_layout()
plt.savefig('eemd_example', dpi=120)
plt.show()

CASE 2: with parabolic option

nIMFs# %%time
# Assign EEMD to `eemd` variable
eemd = EEMD()
# Say we want detect extrema using parabolic method
emd = eemd.EMD
emd.extrema_detection="parabol"
# Execute EEMD on S
eIMFs = eemd.eemd(_s, _t)
nIMFs = eIMFs.shape[0]
# Plot results
plt.figure(figsize=(12,18))
plt.subplot(nIMFs+1, 1, 1)
plt.plot(_t, _s, 'r')
for n in range(nIMFs):
    plt.subplot(nIMFs+1, 1, n+2)
    plt.plot(_t, eIMFs[n], 'g')
    plt.ylabel("eIMF %i" %(n+1))
    plt.locator_params(axis='y', nbins=5)
plt.xlabel("Time [s]")
plt.tight_layout()
plt.savefig('eemd_example', dpi=120)
plt.show()

/usr/local/lib/python3.6/dist-packages/PyEMD/EMD.py:555: RuntimeWarning: divide by zero encountered in true_divide
  a = a/scale
...

"""
Traceback (most recent call last):
  File "/usr/lib/python3.6/multiprocessing/pool.py", line 119, in worker
    result = (True, func(*args, **kwds))
  File "/usr/lib/python3.6/multiprocessing/pool.py", line 44, in mapstar
    return list(map(*args))
  File "/usr/local/lib/python3.6/dist-packages/PyEMD/EEMD.py", line 189, in _trial_update
    return self.emd(self._S+noise, self._T, self.max_imf)
  File "/usr/local/lib/python3.6/dist-packages/PyEMD/EEMD.py", line 197, in emd
    return self.EMD.emd(S, T, max_imf)
  File "/usr/local/lib/python3.6/dist-packages/PyEMD/EMD.py", line 797, in emd
    mean[:] = 0.5*(max_env+min_env)
ValueError: operands could not be broadcast together with shapes (1344,) (1341,) 
"""

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
<ipython-input-84-1cd1e8c94d9b> in <module>
      6 emd.extrema_detection="parabol"
      7 # Execute EEMD on S
----> 8 eIMFs = eemd.eemd(_s, _t)
      9 nIMFs = eIMFs.shape[0]
     10 # Plot results

/usr/local/lib/python3.6/dist-packages/PyEMD/EEMD.py in eemd(self, S, T, max_imf)
    174         # For trial number of iterations perform EMD on a signal
    175         # with added white noise
--> 176         all_IMFs = self.pool.map(self._trial_update, range(self.trials))
    177 
    178         max_imfNo = max([IMFs.shape[0] for IMFs in all_IMFs])

/usr/lib/python3.6/multiprocessing/pool.py in map(self, func, iterable, chunksize)
    264         in a list that is returned.
    265         '''
--> 266         return self._map_async(func, iterable, mapstar, chunksize).get()
    267 
    268     def starmap(self, func, iterable, chunksize=None):

/usr/lib/python3.6/multiprocessing/pool.py in get(self, timeout)
    642             return self._value
    643         else:
--> 644             raise self._value
    645 
    646     def _set(self, i, obj):

ValueError: operands could not be broadcast together with shapes (1344,) (1341,)

rfmiotto commented 4 years ago

Dear Kim, I tried to set the extrema_detection="parabol" with an hypothetical signal, and it worked just fine for me. If you could share the signal that you are using (i.e., the tidal wave), it would be easier for others to test and possibly help you. Best

kimyoungjin06 commented 4 years ago

@rfmiotto Thanks! I upload data at my temporary repository.

,date,tidal_level
0,2010-09-02 00:00:00,317.0
1,2010-09-02 01:00:00,271.0
2,2010-09-02 02:00:00,268.0
3,2010-09-02 03:00:00,310.0
4,2010-09-02 04:00:00,381.0
5,2010-09-02 05:00:00,451.0
6,2010-09-02 06:00:00,500.0
...

data

And the code I used for importing the data.

import pandas as pd

station = stations[0]
df = pd.read_csv('./Data_hour_sorted/%s.csv'%station, index_col=0)

t = pd.to_datetime(df.date).values
s = df.tidal_level.values

win_size = 24*28*2
_s = s[:win_size]
_t = [datetime.utcfromtimestamp(_*1e-9) for _ in t[:win_size].astype(int)] # dt64 to datetime
_t = np.array([time.mktime(x.timetuple()) for x in _t])

After that, same code in original text.

FYI

laszukdawid commented 4 years ago

Hi @kimyoungjin06,

Apologies but it's going to take mea while to test it and get back to you.

The traceback you provided points at scaling line which is based on x-distance between points. My suspicions are your time-array (_t) and your systems data formats.

As for helping you debug, I have a couple of suggestions:

Try first without passing the _t. That will remove suspicion from the _s.
Instead of passing timestamp, pass floats or epoch values. EMD was not created with timestamps in mind (but in the future we should be handling them properly).
Instead of computing EEMD please try EMD. Since EEMD is (EMD + noise) x ManyTimes you should probably get the same error but much quicker. Potentially the message might be more verbose.

Again, apologies for not helping more and quicker. Let me know how your search is going and I'll try to help on the earliest chance.

Dawid

kimyoungjin06 commented 4 years ago

@laszukdawid

Thanks for your kind answer!

I get the result with floats instead of timestamp _t. It is so helpful for me.

Have a nice day without COVID-19!

Young Jin Kim

laszukdawid commented 4 years ago

I understand that the problem is now resolved when using floats instead of timestamps. Great to hear that :)

Sorry that it took me sooo long to get back. Crazy times.

Feel free to reopen if the issue isn't closed.

Hope you're also doing well :)

laszukdawid / PyEMD