Added DTCWT operator

[dikwickley]

Added Dual-Tree Complex Wavelet Transform Operator. Issues: #46

• added operator
• added tests
• added example

[mrava87]

Well done, this is a very good start :)

I have a few general comments (apart from the specific ones in the code review):

• I always find better to make multiple operators for 1D, 2D, ND. Unless the library that we leverage handles that for us, you will see that we end up having a messy code full of if..else statements if we try to pack all together. Also for 1D you have the special case where the array is 1D, for 2D this doesn't make sense (similar for 3D it doesnt make sense to have 1D and 2D arrays).
• For the 1D case we still want to allow ND arrays as in the philosophy of most pylops operators (see other comments for more on this)

[cako]

@dikwickley let me know when this is ready for review. Thanks!

[mrava87]

@cako i was planning to look at it over the weekend (@dikwickley told me is ready) but feel free to go ahead if you have time

[cako]

The code looks solid, thanks for the feature. One think I think would be nice (but is not necessary) is the addition of a more "intuitive" example. Using random input signals is quite hard to interpret. I would prefer something like two wavelets with different frequencies, and show how they (maybe) will be separated in the DTCWT domain. Just my two cents!

I will let @mrava87 have one final look before we approve.

[cako]

Should be pylops.signalprocessing.DTCWT

[cako]

this needs to be wrapped in a pylops.utils.deps check like here:

https://github.com/PyLops/pylops/blob/dev/pylops/_torchoperator.py#LL5-LL10

[mrava87]

I haven't been able to look into this yet, but I agree with @cako on the example. Tests are meant to check the correctness of the code so should use as random signals as possible (generally this may bring to surface problems that specials signals - eg constant - dont). But the examples should be as communicative as possible. Whilst it is not always needed to have a tutorial, and example of DTCWT like @cako suggests would be much nicer for the users than using random signals.

@dikwickley RTD seems to fail, not sure if this is due to your new code (likely not) or some changes in the RTD process... for now let's not worry, once you push changes we will see if things go back to green light :)

[mrava87]

A couple of additional minor comments to avoid diverging from common practice in matvec/rmatvec

[mrava87]

@this is not good, you should use swapaxes like here

pylops/pylops/basicoperators/firstderivative.py

Line 139 in e4d914a

@reshaped(swapaxis=True)

When something is readily available and used in many other operators we should not deviate from it unless there is a special requirement, I do not see it here :)

[mrava87]

Same

[mrava87]

This needs to go all the way to the end of the title

[mrava87]

@dikwickley I am reviewing your code and I found quite a few unclear points:

• your _nd_to_2d does arr_nd.reshape((self.dims[0], -1)) but what happens if the dims you want to apply the transform is 0? Should it not be -1, self.dims[self.axis]? I think yours would just work correctly as a coincidence for 2 arrays if you want to apply it over axis=1
• I dont understand the role of elements = np.prod(T.shape[1:]) it seems to me that the dcwt transform already handles internally multiple dimensions. See an example:

import dtcwt

transform = dtcwt.Transform1d()

x = np.ones(50)
transform.forward(x).lowpass.shape

x = np.ones((50, 10))
transform.forward(x).lowpass.shape

which will give output (14, 1) and (14, 10). Also note that x is a 2d array when you compute elements, why do you need np.prod as by definition you are getting only one element from shape?

For now, I limited myself to make some cleaning of the docstring and include dtcwt in installation files where missing, but we need to make sure what you do above is correct before moving ahead :)

[mrava87]

PS: I made another PR to fix the doc issue and also done a bit of cleaning: requirements-doc.txt isn't really needed, we can just use requirements-dev.txt as we used to do before... so the reason for the conflict. I suggest to keep working on this PR as a draft and once ready make a new clean one

[dikwickley]

@mrava87

• i created _nd_to_2d assuming we would always need to scale down the dimensions of the array to 2 dimensions. For this a check can be added if the dimensions are greater than two, only then, reduce the input.
• elements = np.prod(T.shape[1:]) is important, it converts array of any shape to 2D, eg shape (m, n, o, p...) is reduced to (m, nop*..). dtcwt does handles multiple dimensions only upto 2D https://dtcwt.readthedocs.io/en/0.12.0/reference.html#dtcwt.Transform1d (see parameter x)

[mrava87]

@mrava87

• i created _nd_to_2d assuming we would always need to scale down the dimensions of the array to 2 dimensions. For this a check can be added if the dimensions are greater than two, only then, reduce the input.

I see the point of _nd_to_2d but I am not sure about its implementation. Since you do swapping of dimensions first, you should make sure the 2d array is consistent with it, whilst doing arr_nd.reshape((self.dims[0], -1)) would always make the first dimension of size self.dims[0] no matter if that is what is there after swapping or not. Try a simple example when you swap and dont swap and then run this function, you will see the output will always have the same size, but this is not correct :)

• elements = np.prod(T.shape[1:]) is important, it converts array of any shape to 2D, eg shape (m, n, o, p...) is reduced to (m, n_o_p*..). dtcwt does handles multiple dimensions only upto 2D https://dtcwt.readthedocs.io/en/0.12.0/reference.html#dtcwt.Transform1d (see parameter x)

Ok, I see, but again you have already _nd_to_2d so why do you need to also have elements. If you do print(pyr.lowpass.shape) inside _interpret_coeffs for an example with ND array (say your test test_dtcwt1D_input3D) I see that this is already a 2D-array where the second dimension is the same of that of the input (all your extra dimensions multiplied), so why is this needed self.lowpass_size *= elements?

[mrava87]

@dikwickley i am checking if you are still interested to finalize this PR?

Since we were close perhaps we can try to finalize it (I remember I only had a few concerns with handling of ndarrays which I didn’t agree with your code… other than that everything seemed good to me)

[mrava87]

@dikwickley I created a new PR as this one is quite hold and has some conflicts. In the new PR you will see that the code is changed quite a bit, because as I mentioned above I think your code is not correct for beyond 1d arrays since it moves axis to the last dimension prior to computing the transform, instead this should move to the first dimension.

Let me know if you are still interested to follow that PR or not?

[mrava87]

move to #538