Adding dims to DFT/IDFT operator #13

1oly · 2019-05-03T12:37:48Z

This is adding an optional argument dims to DFT/IDFT. Discussed in #12

Follows the definition from AbstractFFTs: https://github.com/JuliaMath/AbstractFFTs.jl/blob/2f11b0e86c7b1c1660f7513849986e85f58a28b7/src/definitions.jl#L198
by adding dims=1:ndims(x) or dims=1:N to functions and constructors:
DFT(x::AbstractArray{D,N},dims=1:ndims(x)) where {N,D<:Real} or
DFT(dim_in::NTuple{N,Int},dims=1:N) where {N} = DFT(zeros(dim_in),dims)

I hope I caught all cases and methods but please review and let me know if I'm missing anything.

1oly · 2019-05-03T17:30:00Z

Oh, I see the problem. Will fix and update...

nantonel · 2019-05-06T10:02:51Z

Hi @1oly , probably you already know but you can test your changes locally by typing:

] test AbstractOperators

Additionally what I typically do is to comment out most of the tests in test/runtests.jl and other files in the same folder, leaving out only the tests related to the PR. Once they pass I uncomment everything. 😉

1oly · 2019-05-06T10:22:13Z

Thanks, all tips are welcome :)

Currently only two of the tests I added are failing and I'm not quite sure why. It is for IDFT when dims is a scalar (AbstractOperators.IDFT(Float64,(n,n)) is passing).

op = AbstractOperators.IDFT(Float64,(n,n),1)
x1 = randn(n,n)
y1 = test_op(op, x1, fft(randn(n,n)), verb)      # failing here
y2 = ifft(x1,1)

@test all(norm.(y1 .- y2) .<= 1e-12)

op = AbstractOperators.IDFT(Complex{Float64},(n,n),2)
x1 = randn(n,n)+im*randn(n,n)
y1 = test_op(op, x1, fft(randn(n,n)), verb)     # failing here
y2 = ifft(x1,2)

@test all(norm.(y1 .- y2) .<= 1e-12)

So going into test_op I see that its failing due to the dot call:

A = AbstractOperators.IDFT(Float64,(n,n),2)
x = randn(n,n)
y = fft(randn(n,n))

Ax = A*x
Ax2 = similar(Ax)
mul!(Ax2, A, x) #verify in-place linear operator works
@test norm(Ax .- Ax2) <= 1e-8
Acy = A'*y
Acy2 = similar(Acy)
At = AdjointOperator(A)
mul!(Acy2, At, y) #verify in-place linear operator works

@test norm(Acy .- Acy2) <= 1e-8     # passing test

s1 = real(dot(Ax2, y))                # THESE ARE NOT EQUAL
s2 = real(dot(x, Acy2))             # THESE ARE NOT EQUAL

@test abs( s1 - s2 ) < 1e-8       # failing

I think it might be due to a simple error in domain transformation but cannot find the error.

Any ideas?

1oly · 2019-05-06T11:45:54Z

OK, its the scaling factor in plan_ifft that causes the error. The ratio s1/s2=n, should the scaling factor be in the tests or should the IDFT constructor be modified?

nantonel · 2019-05-06T13:38:22Z

OK, its the scaling factor in plan_ifft that causes the error. The ratio s1/s2=n, should the scaling factor be in the tests or should the IDFT constructor be modified?

Well the scaling changes depending on the dimension you're applying the dft/idft to. I guess you could add a scaling parameter to the IDFT object and then modify these lines accordingly.

Of course you should also add some tests of DFT and IDFT with different dimensions settings possibly below the ones that are already there.

1oly · 2019-05-09T11:03:59Z

Starting to look better now. Tests are now passing. Some words on my approach:

I've looked closely at RDFT/IRDFT and AbstractFFTs.jl. I adopted the dims argument from AbstractFFTs, which has no type declaration (this can also be seen in the doc string of e.g. FFTW.fft, which states: The optional dims argument specifies an iterable subset of dimensions (e.g. an integer, range, tuple, or array) to transform along.).
Therefore, I added a new parameter E to DFT and IDFT (similar to IRDFT) to hold dims:

struct DFT{N,
	      C<:RealOrComplex,
	      D<:RealOrComplex,
		  E,
	      T1<:AbstractFFTs.Plan,
	      T2<:AbstractFFTs.Plan} <: FourierTransform{N,C,D,E,T1,T2}
	dim_in::NTuple{N,Int}
	A::T1
	At::T2
end

struct IDFT{N,
	       C<:RealOrComplex,
	       D<:RealOrComplex,
		   E,
	       T1<:AbstractFFTs.Plan,
	       T2<:AbstractFFTs.Plan} <: FourierTransform{N,C,D,E,T1,T2}
	dim_in::NTuple{N,Int}
	A::T1
	At::T2
end

and added dims to the inner constructor:

function IDFT(x::AbstractArray{D,N},dims=1:ndims(x)) where {N,D<:Real}
	y2 = zeros(complex(D),size(x))
	A = plan_ifft(x,dims)
	At = plan_fft(y2,dims)
	dim_in = size(x)
	IDFT{N,Complex{D},D,dims,typeof(A),typeof(At)}(dim_in,A,At)
end

It might be overkill to add ´E´ to both DFT and IDFT when DFT is not using it, but since they are both subtypes of FourierTransform and not LinearOperator as is in RDFT/IRDFT, I think it is necessary. Please correct me if I'm wrong.

To account for the normalization in IDFT, I added the line: y ./= prod(size(b,e) for e in E) instead of y ./= size(b). This should account for the fact that E can be any iterable subset of dimensions as used in FFTW.fft. It feels a little clumsy but works.

Let me know what you think, I'll be happy to update again.

codecov-io · 2019-05-09T11:05:43Z

Codecov Report

Merging #13 into master will increase coverage by 0.11%.
The diff coverage is 92.85%.

@@            Coverage Diff             @@
##           master      #13      +/-   ##
==========================================
+ Coverage   83.81%   83.92%   +0.11%     
==========================================
  Files          49       49              
  Lines        1643     1655      +12     
==========================================
+ Hits         1377     1389      +12     
  Misses        266      266

Impacted Files	Coverage Δ
src/linearoperators/DFT.jl	`92.98% <92.85%> (+1.87%)`	⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update e16bbe4...90cdc82. Read the comment docs.

codecov-io · 2019-05-09T11:05:44Z

Codecov Report

Merging #13 into master will not change coverage.
The diff coverage is 100%.

@@           Coverage Diff           @@
##           master      #13   +/-   ##
=======================================
  Coverage   83.81%   83.81%           
=======================================
  Files          49       49           
  Lines        1643     1643           
=======================================
  Hits         1377     1377           
  Misses        266      266

Impacted Files	Coverage Δ
src/linearoperators/DFT.jl	`91.11% <100%> (ø)`	⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update e16bbe4...b819d96. Read the comment docs.

nantonel · 2019-05-09T19:31:56Z

maybe you could put E = prod(size(b,e) for e in dims) in IDFT{N,Complex{D},D,E,typeof(A),typeof(At)}(dim_in,A,At) since it is an integer?
Or add a new field?
Didn't even know that it was possible to put iterators inside types 😄

nantonel · 2019-05-10T06:04:14Z

test/test_linear_operators.jl

+
+@test all(norm.(y1 .- y2) .<= 1e-12)
+
+op = AbstractOperators.IDFT(Float64,(n,n))


perhaps here something like (n,m) where m != n would make the test more robust.

Good point, will add some more robust tests.

nantonel · 2019-05-10T06:08:53Z

src/linearoperators/DFT.jl

 	      T1<:AbstractFFTs.Plan,
-	      T2<:AbstractFFTs.Plan} <: FourierTransform{N,C,D,T1,T2}
+	      T2<:AbstractFFTs.Plan} <: FourierTransform{N,C,D,E,T1,T2}
 	dim_in::NTuple{N,Int}
 	A::T1
 	At::T2


Just to clarify what I said in the other comment. I think something like this could be more approriate:

struct DFT{N, C<:RealOrComplex, D<:RealOrComplex, E<:Int, T1<:AbstractFFTs.Plan, T2<:AbstractFFTs.Plan} <: FourierTransform{N,C,D,T1,T2} T2<:AbstractFFTs.Plan} <: FourierTransform{N,C,D,E,T1,T2} dim_in::NTuple{N,E} A::T1 At::T2 dims::E end

Yes, much better!

nantonel · 2019-05-10T06:09:28Z

src/linearoperators/DFT.jl

 	       T1<:AbstractFFTs.Plan,
-	       T2<:AbstractFFTs.Plan} <: FourierTransform{N,C,D,T1,T2}
+	       T2<:AbstractFFTs.Plan} <: FourierTransform{N,C,D,E,T1,T2}
 	dim_in::NTuple{N,Int}
 	A::T1
 	At::T2
 end



Of course have the same structure here

nantonel · 2019-05-10T06:12:39Z

src/linearoperators/DFT.jl

+	A = plan_fft(x,dims)
+	At = plan_bfft(y2,dims)
+	dim_in = size(x)
+	DFT{N,Complex{D},D,dims,typeof(A),typeof(At)}(dim_in,A,At)


and hence

d = prod(size(x,e) for e in dims) DFT{N,Complex{D},D,typeof(d),dims,typeof(A),typeof(At)}(dim_in,A,At,d)

also are you sure of d = prod(size(x,e) for e in dims)?

Makes sense to do the normalization here. But just to make sure, the same normalization also applies to mul!:
https://github.com/kul-forbes/AbstractOperators.jl/blob/e16bbe48e82c334e9435b6692e621f72a73373cd/src/linearoperators/DFT.jl#L155-L160

Are we certain that x from the constructor and b from mul! are the same dimensions? I guess the input b::AbstractArray{C,N} suggests so...

Quite certain about d = prod(size(x,e) for e in dims) but will try to add a more convincing comment later :)

1oly · 2019-05-10T09:27:06Z

I may have overly complicated things 😂

It turns out that plan_ifft.scale has this parameter and it is simply enough to add

function mul!(y::AbstractArray{C,N},
              L::AdjointOperator{IDFT{N,C,C,T1,T2}},
              b::AbstractArray{C,N}) where {N,C<:Complex,T1,T2}
	mul!(y,L.A.At,b)
	y .*= L.A.A.scale
end

Should I just try to clean up and commit that then (with the new tests and doc string of coarse)?

nantonel · 2019-05-10T11:42:12Z

I may have overly complicated things 😂

It turns out that plan_ifft.scale has this parameter and it is simply enough to add
function mul!(y::AbstractArray{C,N},
              L::AdjointOperator{IDFT{N,C,C,T1,T2}},
              b::AbstractArray{C,N}) where {N,C<:Complex,T1,T2}
	mul!(y,L.A.At,b)
	y .*= L.A.A.scale
end

Ah that's great!

Should I just try to clean up and commit that then (with the new tests and doc string of coarse)?

I guess the fastest thing to do here is to checkout to master the single file?

git checkout origin/master src/linearoperators/DFT.jl

1oly · 2019-05-10T17:31:55Z

OK, should be final now. A little detour but got to study the code in detail which is always nice.

Thanks for reviewing!

WIP: adding dims to DFT operator

3185641

Added normalization when transforming along arbitrary dims with IDFT

90cdc82

nantonel reviewed May 10, 2019

View reviewed changes

1oly added 2 commits May 10, 2019 19:04

Added scale factor from plan_ifft to IDFT

d44be9a

fixed doc string

b819d96

nantonel merged commit 3b68a8d into kul-optec:master May 13, 2019

nantonel mentioned this pull request May 13, 2019

DFT transform along dims... #12

Closed

nantonel referenced this pull request Jul 15, 2019

deprecated Hadamard and NonLinearCompose

5cf184f

1oly mentioned this pull request Aug 29, 2019

Multidimensional arrays JuliaFirstOrder/ProximalAlgorithms.jl#26

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Adding dims to DFT/IDFT operator #13

Adding dims to DFT/IDFT operator #13

1oly commented May 3, 2019

1oly commented May 3, 2019

nantonel commented May 6, 2019

1oly commented May 6, 2019

1oly commented May 6, 2019

nantonel commented May 6, 2019

1oly commented May 9, 2019

codecov-io commented May 9, 2019

codecov-io commented May 9, 2019 •

edited

nantonel commented May 9, 2019

nantonel May 10, 2019

1oly May 10, 2019

nantonel May 10, 2019

1oly May 10, 2019

nantonel May 10, 2019

nantonel May 10, 2019

nantonel May 10, 2019

1oly May 10, 2019

1oly May 10, 2019

1oly commented May 10, 2019

nantonel commented May 10, 2019

1oly commented May 10, 2019


		@test all(norm.(y1 .- y2) .<= 1e-12)

		op = AbstractOperators.IDFT(Float64,(n,n))

Adding dims to DFT/IDFT operator #13

Adding dims to DFT/IDFT operator #13

Conversation

1oly commented May 3, 2019

1oly commented May 3, 2019

nantonel commented May 6, 2019

1oly commented May 6, 2019

1oly commented May 6, 2019

nantonel commented May 6, 2019

1oly commented May 9, 2019

codecov-io commented May 9, 2019

Codecov Report

codecov-io commented May 9, 2019 • edited

Codecov Report

nantonel commented May 9, 2019

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

1oly commented May 10, 2019

nantonel commented May 10, 2019

1oly commented May 10, 2019

codecov-io commented May 9, 2019 •

edited