PyWavelets / pywt

PyWavelets - Wavelet Transforms in Python
http://pywavelets.readthedocs.org
MIT License
2.04k stars 470 forks source link

Can't remove cA and then reconstruct #302

Closed samyachour closed 7 years ago

samyachour commented 7 years ago

I have some EKG data i'm running wavelet decompisition on (thanks for this library!!). Here is an example of my code running level 6 db4 transformation on a 1 dimensional data set (list of voltage values, i.e. [1, 3, -1, 2...], the time intervals are +0.003s).

I used wavedecn because later I want to use waverecn later (where my problem begins)

wavelet = 'db4'
levels = 6
mode = 'constant'
coeffs = pywt.wavedecn(cA, wavelet, level=levels, mode=mode)

So I rebuilt using all the levels to get the original signal back. That was fine.

rebuilt = pywt.waverecn(coeffs, wavelet, mode=mode)
plotWave(rebuilt, "rebuilt1", "hopefully correct indices")

Then if I try to remove cA and cD6

coeffs.pop(6)
coeffs.pop(0)
rebuilt = pywt.waverecn(coeffs, wavelet, mode=mode)
plotWave(rebuilt, "rebuilt2", "hopefully correct indices")

I get all sorts of errors. If I pop cA coeffs.pop(0) or set it to an empty list coeffs[0] = [] I get, "All coefficients must have a matching number of dimensions", if I set coeffs[0] = None it says "The following detail coefficients were set to None: ['a']."

I fixed my error by forking your library and using my version, changing the code on line 40 in multilevel.py from:

a, ds = coeffs[0], coeffs[1:]

to

# if we pop(0) and remove cA, all the elements will be of type dictionary
if type(coeffs[0]) is dict:
        a = None
        ds = coeffs
    else:
        a, ds = coeffs[0], coeffs[1:]

I'm not 100% what effect I had, my code works now and I can remove cA, c1, and c6 and everything runs fine. I think before _fix_coeffs(coeffs) had a problem with setting cA to none (which I want to do if I want to ignore it in reconstruction), so I set it programmatically in the library. Is this correct? Why did my change work?

grlee77 commented 7 years ago

Will take a closer look at this soon. In the meantime, can you tell me which version of PyWavelets you are working from? One related issue with coefficients as None was fixed recently by #291 (which was included in the 0.5.2 release made last week).

grlee77 commented 7 years ago

Okay, a few more comments:

removing c[0] altogether is not supported. Setting it to c[0] = None is supported and seems to be working as of 0.5.2 (and in current master). Another equivalent alternative is c[0] = np.zeros_like(c[0]) if which may be clearer in intent.

Using pop to remove detail coefficients should not be done (unless it is the last element of the coefficient list). Ideally wavedecn should have also been updated as idwtn was in #291 so that detail coefficients as None is supported in the multilevel transform. A workaround is to set the detail coefficients to an empty dictionary at the desired level (e.g. c[3] = {} to act is if all detail coeffs were zero). Out of curiosity is there a concrete application where you would want to discard specific stages of detail coeffs?

A small concrete example of the above:

import numpy as np
import pywt

cam = pywt.data.camera().astype(np.float32)

c = pywt.wavedecn(cam, 'db4', level=6, mode='constant')

c[0] = None
# alternatively:  c[0] = np.zeros_like(c[0])
cr = pywt.waverecn(c, 'db4', mode='constant')

# remove finest detail coefficients
c = c[:-1]
cr_256 = pywt.waverecn(c, 'db4', mode='constant')

# try various methods of removing an intermediate stage of detail coefficients
if False:
    # Case 1:  pop from list.
    # This will not work and I don't think should be supported
    c.pop(2)
    cr_ = pywt.waverecn(c, 'db4', mode='constant')
elif False:
    # Case2:  set individual detail coeffs to None
    # currently does not work, but perhaps it should have a modification
    # similar to that done in #291 for idwtn so that this would be supported.
    c[2] = {k: None for k, v in c[2].items()}
    cr_ = pywt.waverecn(c, 'db4', mode='constant')
else:
    # Case 3: set to an empty dict.  This currently works
    c[2] = {}
    cr_ = pywt.waverecn(c, 'db4', mode='constant')

I can make a simple PR so that Case 2 from the example above is supported in addition to Case 3.

Also open to other suggestions as to what might be a clearer way to do this.

edit as discussed below, case 3 will also not always get the size correct and should also not be recommended. In general the shape of the detail coeffs needs to be known to get the size right. The only truly safe way is to use something like:

c[2] = {k: np.zeros_like(v) for k, v in c[2].items()}
grlee77 commented 7 years ago

at a minimum, whether we decide to keep cases 2 and/or 3, we should add a test so this will not become broken in the future

grlee77 commented 7 years ago

One other comment. If you remove the last stage detail coeffs the reconstructed volume will be smaller than the input as the 6-level transform will be truncated to 5-levels. If you instead set it to an empty dict, the final stage will be run as if all detail coefficients were zeros, but due to a subtle implementation detail that relies on shapes of both the approx and detail coefficients to determine how to appropriately shape the output, the array may be slightly larger than the original input by some small padding amount. The only workaround I see for this is as in the final case below:

import numpy as np
import pywt

cam = pywt.data.camera().astype(np.float32)  # shape (512, 512)

c = pywt.wavedecn(cam, 'db4', level=6, mode='constant')

if False:
    # remove finest detail coefficients (results in a 5 level inverse transform
    c = c[:-1]
    cr_256 = pywt.waverecn(c, 'db4', mode='constant')
elif False:
   # remains a 6-level transform, but as if all detail coeffs at the finest detail were zeros
    c[-1] = {}
    cr_514 = pywt.waverecn(c, 'db4', mode='constant')
    # shape is (514, 514), not (512, 512) as expected!
    cr_512 = cr_514[:-2, :-2]  # have to manually truncate output to the original shape
else:
    c[-1] = {k: np.zeros_like(v) for k, v in c[-1].items()}
    cr_512 = pywt.waverecn(c, 'db4', mode='constant')
    # now size (512, 512) as expected
samyachour commented 7 years ago

@grlee77 your PR (#303) was very helpful in solving this issue. I am not familiar with wavelet theory, so it's news to me that omitting detail coefficients is different from setting them to zero (though intuitively that makes sense!)

Also, not sure if it matters, but on your PR you have a typo on line 470 in _multilevel.py.

Thanks again, everything is in working order 👍