Conversion between Fourier spaces (#41)

* conversion between Fourier

* Tests

* require one endpoint to match

* Tests for coefficients

* cancel concurrent CI runs

* remove rounding in div

* use union instead of union_by_union_rule

* fix indexing for higher multiples

* Add tests

* remove unused names

* fix sum with scaled/sum domains

* fix bugfix with pts

* clean-up imports

---------

Co-authored-by: Jishnu Bhattacharya <jishnub@users.noreply.github.com>

Author: jishnub

.github/workflows/ci.yml

@@ -12,6 +12,11 @@ on:
       - 'LICENSE'
       - 'README.md'
       - '.github/workflows/TagBot.yml'
+
+concurrency:
+  group: build-${{ github.event.pull_request.number || github.ref }}
+  cancel-in-progress: true
+
 jobs:
   pre_job:
     # continue-on-error: true # Uncomment once integration is finished

src/ApproxFunFourier.jl

@@ -43,9 +43,9 @@ using DomainSets
 import DomainSets: Domain, indomain, UnionDomain, Point, Interval,
             boundary, rightendpoint, leftendpoint, endpoints
 
-import Base: convert, getindex, *, +, -, ==,  /, eltype,
-            show, sum, cumsum, conj, issubset, first, last, rand, setdiff,
-            angle, isempty, zeros, one, promote_rule, real, imag
+import Base: convert, getindex, *, +, -, ==, /, eltype,
+            show, sum, conj, issubset, first, last, rand, setdiff,
+            union, angle, isempty, one, promote_rule, real, imag
 
 import LinearAlgebra: norm, mul!, isdiag
 
@@ -531,6 +531,74 @@ canonicalspace(S::Laurent{DD,RR}) where {DD<:PeriodicSegment,RR} = Fourier(domai
 canonicalspace(S::Fourier{DD,RR}) where {DD<:Circle,RR} = Laurent(domain(S))
 canonicalspace(S::Laurent{DD,RR}) where {DD<:PeriodicLine,RR} = S
 
+# avoid cyclic recursion between CosSpace and Fourier
+# check if one domain is a multiple of another
+# return the multiple
+period(F::Fourier) = period(domain(F))
+function period(A::Domain)
+    la, ra = leftendpoint(A), rightendpoint(A)
+    ra - la
+end
+function domainsmultiple(A::Domain, B::Domain)
+    isapprox(A, B) && return 1
+    la, ra = leftendpoint(A), rightendpoint(A)
+    lb, rb = leftendpoint(B), rightendpoint(B)
+    pa = ra - la
+    pb = rb - lb
+    dmin, dmax = minmax(pa, pb)
+    T = promote_type(eltype(A), eltype(B))
+    if isapprox(mod(dmax, dmin), 0, atol=eps(T)) && isapprox(la, lb, atol=eps(T)) || isapprox(ra, rb, atol=eps(T))
+        return round(Int, dmax/dmin, RoundNearest)
+    end
+    return nothing
+end
+domainsmultiple(A::Space, B::Space) = domainsmultiple(map(domain, (A,B))...)
+# check if one domain may be scaled to obtain the other
+function domainsscaled(A::Domain, B::Domain)
+    la, ra = leftendpoint(A), rightendpoint(A)
+    lb, rb = leftendpoint(B), rightendpoint(B)
+    pa = ra - la
+    pb = rb - lb
+    dmin, dmax = minmax(pa, pb)
+    T = promote_type(eltype(A), eltype(B))
+    if isapprox(la, lb, atol=eps(T)) || isapprox(ra, rb, atol=eps(T))
+        return Rational(dmax/dmin)
+    end
+    return nothing
+end
+domainsscaled(A::Space, B::Space) = domainsscaled(map(domain, (A,B))...)
+function union(A::Fourier, B::Fourier)
+    dA, dB = map(domain, (A,B))
+    AB = domainsmultiple(dA, dB)
+    isnothing(AB) || return Fourier(max(dA, dB))
+    scale = domainsscaled(dA, dB)
+    isnothing(scale) || return Fourier(max(dA, dB) * denominator(scale))
+    SumSpace(A, B)
+end
+_ind(i, n) = 2n * div(i, 2) + isodd(i)
+_invind(i, n) = 2div(i - isodd(i), 2n) + isodd(i)
+function coefficients(v::AbstractVector, A::Fourier, B::Fourier)
+    n = domainsmultiple(A, B)
+    isnothing(n) && error("spaces are not compatible")
+    n == 1 && return copy(v)
+    if period(A) < period(B)
+        v2 = zeros(eltype(v), _ind(length(v), n))
+        for (i, vi) in enumerate(v)
+            v2[_ind(i, n)] = vi
+        end
+    else
+        for (ind, vi) in enumerate(v)
+            if isodd(div(ind, n)) && !isapprox(vi, 0, atol=eps(eltype(vi)))
+                throw(ArgumentError("coefficients incompatible with space conversion"))
+            end
+        end
+        v2 = zeros(eltype(v), _invind(length(v), n))
+        for i in eachindex(v2)
+            v2[i] = v[_ind(i, n)]
+        end
+    end
+    return v2
+end
 
 ## Ones and zeros
 

test/runtests.jl

@@ -189,14 +189,62 @@ end
     @test norm(ApproxFunBase.Reverse(Fourier())*Fun(t->cos(cos(t-0.2)-0.1),Fourier()) - Fun(t->cos(cos(-t-0.2)-0.1),Fourier())) < 10eps()
     @test norm(ApproxFunBase.ReverseOrientation(Fourier())*Fun(t->cos(cos(t-0.2)-0.1),Fourier()) - Fun(t->cos(cos(t-0.2)-0.1),Fourier(PeriodicSegment(2π,0)))) < 10eps()
 
-    @testset "issue #741" begin
+    @testset "ApproxFun issue #741" begin
         c = Fun(cos, Fourier())
         @test roots(c) ≈ (1:2:3) * pi/2
         c = Fun(cos, Fourier(0..4pi))
         @test roots(c) ≈ (1:2:7) * pi/2
         c = Fun(cos, Fourier(-2pi..2pi))
         @test roots(c) ≈ (-3:2:3) * pi/2
     end
+
+    @testset "ApproxFun issue #768" begin
+        s1 = Fourier(0..pi)
+        f1 = Fun(t -> sin(2*t),s1)
+        @testset for d2 in [0..2pi, 0..3pi]
+            s2 = Fourier(d2)
+            f2 = Fun(t -> sin(2*t),s2)
+            g = f1 + f2
+            @test g(pi/3) ≈ f1(pi/3) + f2(pi/3)
+        end
+        f1 = Fun(sin, Fourier(0..2pi));
+        f2 = Fun(x->sin((2/3)x), Fourier(0..3pi));
+        g = f1 + f2
+        pts = [0:6;]*pi
+        @test g.(pts) ≈ f1.(pts) .+ f2.(pts)
+        h = Fun(x-> sin(x)+sin(2*x/3),Fourier(0..6pi))
+        @test g.(pts) ≈ h.(pts)
+
+        f1 = Fun(x->sin(2x), Fourier(pi..2pi))
+        f2 = Fun(x->sin(2x), Fourier(3pi..4pi))
+        g = f1 + f2
+        pts = [0:5;]*pi
+        @test g.(pts) ≈ f1.(pts) .+ f2.(pts)
+    end
+
+    @testset "coeff conversion" begin
+        f = t->1+2sin(t)+3cos(t)
+        f1 = Fun(f, Fourier(0..2pi))
+
+        f2 = Fun(f, Fourier(0..4pi))
+        @test coefficients(coefficients(f1), space(f1), space(f2)) ≈ coefficients(f2)
+
+        f3 = Fun(f, Fourier(0..8pi))
+        @test coefficients(coefficients(f1), space(f1), space(f3)) ≈ coefficients(f3)
+
+        @test coefficients([1,0,0,2,3], Fourier(0..4pi), Fourier(0..2pi)) ≈ [1,2,3]
+        @test coefficients([1,2,3], Fourier(0..2pi), Fourier(0..4pi)) ≈ [1; zeros(2); [2,3]]
+        @test coefficients([1; zeros(2); [2,3]], Fourier(0..4pi), Fourier(0..2pi)) ≈ [1,2,3]
+        @test coefficients([4,1], Fourier(0..2pi), Fourier(0..4pi)) ≈ [4; zeros(2); 1]
+        @test coefficients([4,1], Fourier(0..2pi), Fourier(0..8pi)) ≈ [4; zeros(6); 1]
+        @test coefficients([4,1,2], Fourier(0..2pi), Fourier(0..4pi)) ≈ [4; zeros(2); [1, 2]]
+        @test coefficients([4,1,2], Fourier(0..2pi), Fourier(0..8pi)) ≈ [4; zeros(6); [1, 2]]
+        @test coefficients([4,1,2,3], Fourier(0..2pi), Fourier(0..4pi)) ≈ [4; zeros(2); [1, 2]; zeros(2); 3]
+        @test coefficients([4,1,2,3], Fourier(0..2pi), Fourier(0..6pi)) ≈ [4; zeros(4); [1, 2]; zeros(4); 3]
+        @test coefficients([4,1,2,3], Fourier(0..2pi), Fourier(0..8pi)) ≈ [4; zeros(6); [1, 2]; zeros(6); 3]
+        @test coefficients([4; zeros(4); [1, 2]; zeros(4); 3], Fourier(0..6pi), Fourier(0..2pi)) ≈ [4,1,2,3]
+        @test coefficients([4; zeros(6); [1, 2]; zeros(6); 3], Fourier(0..8pi), Fourier(0..2pi)) ≈ [4,1,2,3]
+    end
 end
 
 @testset "Laurent" begin