search

laszukdawid / PyEMD

Python implementation of Empirical Mode Decompoisition (EMD) method
https://pyemd.readthedocs.io/
Apache License 2.0
867 stars 224 forks source link

operands could not be broadcast together with shapes #79

Closed kimyoungjin06 closed 4 years ago

kimyoungjin06 commented 4 years ago

Problem

data

image

      [[1977.        ,   38.        ],
       [1977.00011416,   48.        ],
       [1977.00022831,   65.        ],
       [1977.00034247,   85.        ],
       [1977.00045662,  105.        ],
       [1977.00057078,  120.        ],
       [1977.00068493,  123.        ],
       [1977.00079909,  120.        ],
       [1977.00091324,  106.        ],
       [1977.0010274 ,   88.        ],
       [1977.00114155,   72.        ],
       [1977.00125571,   59.        ],
       [1977.00136986,   50.        ],
       [1977.00148402,   51.        ],
       [1977.00159817,   63.        ],
       [1977.00171233,   79.        ],
       [1977.00182648,   96.        ],
       [1977.00194064,  109.        ],
       [1977.00205479,  115.        ],
       [1977.00216895,  113.        ],
       [1977.00228311,  102.        ],
       [1977.00239726,   84.        ],
       [1977.00251142,   59.        ],
       [1977.00262557,   38.        ],
       [1977.00273973,   29.        ],
       [1977.00285388,   32.        ],
       [1977.00296804,   49.        ],
       [1977.00308219,   71.        ],
       [1977.00319635,   93.        ],
       [1977.0033105 ,  116.        ],
       [1977.00342466,  130.        ],
       [1977.00353881,  134.        ],
       [1977.00365297,  127.        ],
       [1977.00376712,  111.        ],
       [1977.00388128,   85.        ],
       [1977.00399543,   61.        ],
       [1977.00410959,   44.        ],
       [1977.00422374,   39.        ],
       [1977.0043379 ,   46.        ],
       [1977.00445205,   62.        ],
       [1977.00456621,   86.        ],
       [1977.00468037,  111.        ],
       [1977.00479452,  127.        ],
       [1977.00490868,  133.        ],
       [1977.00502283,  124.        ],
       [1977.00513699,  102.        ],
       [1977.00525114,   68.        ],
       [1977.0053653 ,   37.        ],
       [1977.00547945,   17.        ],
       [1977.00559361,   13.        ],
       [1977.00570776,   27.        ],
       [1977.00582192,   50.        ],
       [1977.00593607,   80.        ],
       [1977.00605023,  106.        ],
       [1977.00616438,  127.        ],
       [1977.00627854,  139.        ],
       [1977.00639269,  138.        ],
       [1977.00650685,  124.        ],
       [1977.006621  ,   92.        ],
       [1977.00673516,   62.        ],
       [1977.00684932,   35.        ],
       [1977.00696347,   22.        ],
       [1977.00707763,   24.        ],
       [1977.00719178,   39.        ],
       [1977.00730594,   60.        ],
       [1977.00742009,   85.        ],
       [1977.00753425,  110.        ],
       [1977.0076484 ,  124.        ],
       [1977.00776256,  124.        ],
       [1977.00787671,  113.        ],
       [1977.00799087,   87.        ],
       [1977.00810502,   57.        ],
       [1977.00821918,   26.        ],
       [1977.00833333,   11.        ],
       [1977.00844749,   12.        ],
       [1977.00856164,   27.        ],
       [1977.0086758 ,   57.        ],
       [1977.00878995,   88.        ],
       [1977.00890411,  120.        ],
       [1977.00901826,  143.        ],
       [1977.00913242,  156.        ],
       [1977.00924658,  148.        ],
       [1977.00936073,  121.        ],
       [1977.00947489,   88.        ],
       [1977.00958904,   52.        ],
       [1977.0097032 ,   33.        ],
       [1977.00981735,   23.        ],
       [1977.00993151,   26.        ],
       [1977.01004566,   46.        ],
       [1977.01015982,   73.        ],
       [1977.01027397,  105.        ],
       [1977.01038813,  130.        ],
       [1977.01050228,  143.        ],
       [1977.01061644,  135.        ],
       [1977.01073059,  108.        ],
       [1977.01084475,   76.        ],
       [1977.0109589 ,   40.        ],
       [1977.01107306,   15.        ],
       [1977.01118721,    5.        ],
       [1977.01130137,   12.        ],
       [1977.01141553,   37.        ]]

example code

t, s = df.date2.values[:100], df.tidal_level.values[:100]+2000.
eemd = EEMD(trials=1000)
emd = eemd.EMD
emd.extrema_detection="parabol"
eIMFs = eemd.eemd(s, t)

Error log

/usr/local/lib/python3.6/dist-packages/PyEMD/EMD.py:555: RuntimeWarning: divide by zero encountered in true_divide
  a = a/scale
/usr/local/lib/python3.6/dist-packages/PyEMD/EMD.py:556: RuntimeWarning: divide by zero encountered in true_divide
  b = b/scale
...
/usr/local/lib/python3.6/dist-packages/PyEMD/EMD.py:559: RuntimeWarning: invalid value encountered in true_divide
  tVertex = -0.5*b/a
...
/usr/local/lib/python3.6/dist-packages/PyEMD/EMD.py:560: RuntimeWarning: invalid value encountered in less
  idx = np.r_[tVertex<T0+0.5*(Tn-T0)] & np.r_[tVertex>=T0-0.5*(T0-Tp)]
"""
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 (100,) (82,) 
"""
The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
<ipython-input-27-2d62e661cc39> in <module>
     10 emd.extrema_detection="parabol"
     11 # Execute EEMD on S
---> 12 eIMFs = eemd.eemd(s, t)

/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 (100,) (82,)
laszukdawid commented 4 years ago

Interesting. Thanks for letting me know. I'll try to take at look at this in the coming days and hopefully get some more details.

Directly from the error I can tell is that for a single proto-IMF at least one spanned envelope has smaller length (82) than expected (100). A reason for this might be that there aren't enough extrema to span it (but that's a bug that should have been reported), or that some extrema were identified incorrectly (either bug or numerical instability). Given that you have big temporal resolution (small difference in time values) and integers as values, may I suggest to try something like:

The reason why these might help is that in parabol mode, we're estimating extrema with a parabol interpolation. These require differences in temporal axis and plenty a couple multiplication (code is here). The x-difference is tiny and then dividing something big by that will explode.

Still that's likely a bug.

laszukdawid commented 4 years ago

Just confirmed. If you don't provide time array t, or scale it appropriately, then all is fine. The exceptions I'm getting are indeed in small (zero) number divisions.

laszukdawid commented 4 years ago

In addition I added time array normalization in 664aa28 so scaling isn't necessary (though still suggested).

Hope this helps 👍