Merge pull request #184 from vpuri3/fixes

Fixes

Author: xtalax

src/basic.jl

@@ -82,12 +82,12 @@ for op in (
     end
 end
 
-function Base.:\(::IdentityOperator, A::AbstractSciMLOperator)
+function Base.:\(ii::IdentityOperator, A::AbstractSciMLOperator)
     @assert size(A, 1) == ii.len
     A
 end
 
-function Base.:/(A::AbstractSciMLOperator, ::IdentityOperator)
+function Base.:/(A::AbstractSciMLOperator, ii::IdentityOperator)
     @assert size(A, 2) == ii.len
     A
 end
@@ -330,8 +330,9 @@ AddedOperator(L::AbstractSciMLOperator) = L
 # constructors
 Base.:+(A::AbstractSciMLOperator, B::AbstractMatrix) = A + MatrixOperator(B)
 Base.:+(A::AbstractMatrix, B::AbstractSciMLOperator) = MatrixOperator(A) + B
-Base.:+(ops::AbstractSciMLOperator...) = AddedOperator(ops...)
 
+Base.:+(ops::AbstractSciMLOperator...) = reduce(+, ops)
+Base.:+(A::AbstractSciMLOperator, B::AbstractSciMLOperator) = AddedOperator(A, B)
 Base.:+(A::AbstractSciMLOperator, B::AddedOperator) = AddedOperator(A, B.ops...)
 Base.:+(A::AddedOperator, B::AbstractSciMLOperator) = AddedOperator(A.ops..., B)
 Base.:+(A::AddedOperator, B::AddedOperator) = AddedOperator(A.ops..., B.ops...)
@@ -471,16 +472,15 @@ function ComposedOperator(ops::AbstractSciMLOperator...; cache = nothing)
 end
 
 # constructors
-Base.:∘(ops::AbstractSciMLOperator...) = ComposedOperator(ops...)
-Base.:∘(A::ComposedOperator, B::ComposedOperator) = ComposedOperator(A.ops..., B.ops...)
-Base.:∘(A::AbstractSciMLOperator, B::ComposedOperator) = ComposedOperator(A, B.ops...)
-Base.:∘(A::ComposedOperator, B::AbstractSciMLOperator) = ComposedOperator(A.ops..., B)
-
-Base.:*(ops::AbstractSciMLOperator...) = ComposedOperator(ops...)
-Base.:*(A::AbstractSciMLOperator, B::AbstractSciMLOperator) = ∘(A, B)
-Base.:*(A::ComposedOperator, B::AbstractSciMLOperator) = ∘(A.ops[1:end-1]..., A.ops[end] * B)
-Base.:*(A::AbstractSciMLOperator, B::ComposedOperator) = ∘(A * B.ops[1], B.ops[2:end]...)
-Base.:*(A::ComposedOperator, B::ComposedOperator) = ComposedOperator(A.ops..., B.ops...)
+for op in (
+           :*, :∘,
+          )
+    @eval Base.$op(ops::AbstractSciMLOperator...) = reduce($op, ops)
+    @eval Base.$op(A::AbstractSciMLOperator, B::AbstractSciMLOperator) = ComposedOperator(A, B)
+    @eval Base.$op(A::ComposedOperator, B::AbstractSciMLOperator) = ComposedOperator(A.ops..., B)
+    @eval Base.$op(A::AbstractSciMLOperator, B::ComposedOperator) = ComposedOperator(A, B.ops...)
+    @eval Base.$op(A::ComposedOperator, B::ComposedOperator) = ComposedOperator(A.ops..., B.ops...)
+end
 
 for op in (
            :*, :∘,
@@ -606,11 +606,20 @@ function cache_self(L::ComposedOperator, u::AbstractVecOrMat)
     K = size(u, 2)
     cache = (zero(u),)
     for i in reverse(2:length(L.ops))
+        op = L.ops[i]
 
-        M = size(L.ops[i], 1)
-        T = promote_type(eltype.((L.ops[i], cache[1]))...)
+        M = size(op, 1)
         sz = u isa AbstractMatrix ? (M, K) : (M,)
 
+        T = if op isa FunctionOperator # 
+            # FunctionOperator isn't guaranteed to play by the rules of
+            # `promote_type`. For example, an rFFT is a complex operation
+            # that accepts and complex vector and returns a real one.
+            op.traits.eltypes[2]
+        else
+            promote_type(eltype.((op, cache[1]))...)
+        end
+
         cache = (similar(u, T, sz), cache...)
     end
 
@@ -623,12 +632,12 @@ function cache_internals(L::ComposedOperator, u::AbstractVecOrMat)
         L = cache_self(L, u)
     end
 
-    vecs = L.cache
+    ops = ()
     for i in reverse(1:length(L.ops))
-        @set! L.ops[i] = cache_operator(L.ops[i], vecs[i])
+        ops = (cache_operator(L.ops[i], L.cache[i]), ops...)
     end
 
-    L
+    @set! L.ops = ops
 end
 
 function LinearAlgebra.mul!(v::AbstractVecOrMat, L::ComposedOperator, u::AbstractVecOrMat)

src/func.jl

@@ -86,10 +86,11 @@ end
 
 function FunctionOperator(op,
                           input::AbstractVecOrMat,
-                          output::AbstractVecOrMat =  input;
+                          output::AbstractVecOrMat = input;
 
                           isinplace::Union{Nothing,Bool}=nothing,
                           outofplace::Union{Nothing,Bool}=nothing,
+                          isconstant::Bool = false,
                           has_mul5::Union{Nothing,Bool}=nothing,
                           cache::Union{Nothing, NTuple{2}}=nothing,
                           T::Union{Type{<:Number},Nothing}=nothing,
@@ -112,8 +113,10 @@ function FunctionOperator(op,
                           isposdef::Bool = false,
                          )
 
+    # store eltype of input/output for caching with ComposedOperator.
+    eltypes = eltype.((input, output))
     sz = (size(output, 1), size(input, 1))
-    T  = isnothing(T) ? promote_type(eltype.((input, output))...) : T
+    T  = isnothing(T) ? promote_type(eltypes...) : T
     t  = isnothing(t) ? zero(real(T)) : t
 
     isinplace = if isnothing(isinplace)
@@ -164,6 +167,7 @@ function FunctionOperator(op,
 
     traits = (;
               islinear = islinear,
+              isconstant = isconstant,
 
               opnorm = opnorm,
               issymmetric = issymmetric,
@@ -176,6 +180,7 @@ function FunctionOperator(op,
               ifcache = ifcache,
               T = T,
               size = sz,
+              eltypes = eltypes,
              )
 
     L = FunctionOperator(
@@ -197,6 +202,11 @@ function FunctionOperator(op,
 end
 
 function update_coefficients(L::FunctionOperator, u, p, t)
+
+    if isconstant(L)
+        return L
+    end
+
     @set! L.op = update_coefficients(L.op, u, p, t)
     @set! L.op_adjoint = update_coefficients(L.op_adjoint, u, p, t)
     @set! L.op_inverse = update_coefficients(L.op_inverse, u, p, t)
@@ -209,6 +219,11 @@ function update_coefficients(L::FunctionOperator, u, p, t)
 end
 
 function update_coefficients!(L::FunctionOperator, u, p, t)
+
+    if isconstant(L)
+        return L
+    end
+
     for op in getops(L)
         update_coefficients!(op, u, p, t)
     end
@@ -250,6 +265,7 @@ function Base.adjoint(L::FunctionOperator)
 
     traits = L.traits
     @set! traits.size = reverse(size(L))
+    @set! traits.eltypes = reverse(traits.eltypes)
 
     p = L.p
     t = L.t
@@ -284,6 +300,7 @@ function Base.inv(L::FunctionOperator)
 
     traits = L.traits
     @set! traits.size = reverse(size(L))
+    @set! traits.eltypes = reverse(traits.eltypes)
 
     @set! traits.opnorm = if traits.opnorm isa Number
         1 / traits.opnorm
@@ -353,8 +370,8 @@ function getops(L::FunctionOperator)
     ops
 end
 
-#TODO - isconstant(L::FunctionOperator)
 islinear(L::FunctionOperator) = L.traits.islinear
+isconstant(L::FunctionOperator) = L.traits.isconstant
 has_adjoint(L::FunctionOperator) = !(L.op_adjoint isa Nothing)
 has_mul(L::FunctionOperator{iip}) where{iip} = true
 has_mul!(L::FunctionOperator{iip}) where{iip} = iip

src/scalar.jl

@@ -214,10 +214,11 @@ end
 for op in (
            :*, :∘,
           )
-    @eval Base.$op(ops::AbstractSciMLScalarOperator...) = ComposedScalarOperator(ops...)
-    @eval Base.$op(A::ComposedScalarOperator, B::ComposedScalarOperator) = ComposedScalarOperator(A.ops..., B.ops...)
-    @eval Base.$op(A::AbstractSciMLScalarOperator, B::ComposedScalarOperator) = ComposedScalarOperator(A, B.ops...)
+    @eval Base.$op(ops::AbstractSciMLScalarOperator...) = reduce($op, ops)
+    @eval Base.$op(A::AbstractSciMLScalarOperator, B::AbstractSciMLScalarOperator) = ComposedScalarOperator(A, B)
     @eval Base.$op(A::ComposedScalarOperator, B::AbstractSciMLScalarOperator) = ComposedScalarOperator(A.ops..., B)
+    @eval Base.$op(A::AbstractSciMLScalarOperator, B::ComposedScalarOperator) = ComposedScalarOperator(A, B.ops...)
+    @eval Base.$op(A::ComposedScalarOperator, B::ComposedScalarOperator) = ComposedScalarOperator(A.ops..., B.ops...)
 
     for T in SCALINGNUMBERTYPES[2:end]
         @eval Base.$op(α::AbstractSciMLScalarOperator, x::$T) = ComposedScalarOperator(α, ScalarOperator(x))

test/basic.jl

@@ -171,6 +171,15 @@ end
 
     v=rand(N,K); @test mul!(v, op, u) ≈ (A+B) * u
     v=rand(N,K); w=copy(v); @test mul!(v, op, u, α, β) ≈ α*(A+B)*u + β*w
+
+    # ensure AddedOperator doesn't nest
+    A = MatrixOperator(rand(N, N))
+    L = A + (A + A) + A
+    @test L isa AddedOperator
+    for op in L.ops
+        @test !isa(op, AddedOperator)
+    end
+
 end
 
 @testset "ComposedOperator" begin
@@ -221,6 +230,14 @@ end
     v=rand(N,K); @test ldiv!(v, op, u) ≈ (A * B * C) \ u
     v=copy(u);   @test ldiv!(op, u)    ≈ (A * B * C) \ v
 
+    # ensure composedoperators doesn't nest
+    A = MatrixOperator(rand(N, N))
+    L = A * (A * A) * A
+    @test L isa ComposedOperator
+    for op in L.ops
+        @test !isa(op, ComposedOperator)
+    end
+
     # Test caching of composed operator when inner ops do not support Base.:*
     # ComposedOperator caching was modified in PR # 174
     inner_op = qr(MatrixOperator(rand(N, N)))

test/scalar.jl

@@ -1,6 +1,14 @@
 #
-using SciMLOperators, LinearAlgebra
-using Random
+using SciMLOperators
+using SciMLOperators: AbstractSciMLScalarOperator,
+                      ComposedScalarOperator,
+                      AddedScalarOperator,
+                      InvertedScalarOperator,
+                      IdentityOperator,
+                      AddedOperator,
+                      ScaledOperator
+
+using LinearAlgebra, Random
 
 Random.seed!(0)
 N = 8
@@ -38,31 +46,31 @@ K = 12
 
     # Test that ScalarOperator's remain AbstractSciMLScalarOperator's under common ops
     β = α + α
-    @test β isa SciMLOperators.AddedScalarOperator
+    @test β isa AddedScalarOperator
     @test β * u ≈ x * u + x * u
     @inferred convert(Float32, β)
     @test convert(Number, β) ≈ x + x
 
     β = α * α
-    @test β isa SciMLOperators.ComposedScalarOperator
+    @test β isa ComposedScalarOperator
     @test β * u ≈ x * x * u
     @inferred convert(Float32, β)
     @test convert(Number, β) ≈ x * x
 
     β = inv(α)
-    @test β isa SciMLOperators.InvertedScalarOperator
+    @test β isa InvertedScalarOperator
     @test β * u ≈ 1 / x * u
     @inferred convert(Float32, β)
     @test convert(Number, β) ≈ 1 / x
 
     β = α * inv(α)
-    @test β isa SciMLOperators.ComposedScalarOperator
+    @test β isa ComposedScalarOperator
     @test β * u ≈ u
     @inferred convert(Float32, β)
     @test convert(Number, β) ≈ true
 
     β = α / α
-    @test β isa SciMLOperators.ComposedScalarOperator
+    @test β isa ComposedScalarOperator
     @test β * u ≈ u
     @inferred convert(Float32, β)
     @test convert(Number, β) ≈ true
@@ -77,6 +85,14 @@ K = 12
         @test (α / op) * u ≈ (op \ α) * u ≈ α * (op \ u)
         @test (op / α) * u ≈ (α \ op) * u ≈ 1/α * op * u
     end
+
+    # ensure composedscalaroperators doesn't nest
+    α = ScalarOperator(rand())
+    L = α * (α * α) * α
+    @test L isa ComposedScalarOperator
+    for op in L.ops
+        @test !isa(op, ComposedScalarOperator)
+    end
 end
 
 @testset "ScalarOperator update test" begin

test/total.jl

@@ -35,11 +35,17 @@ K = 12
     ik = im * DiagonalOperator(k)
     Dx = ftr \ ik * ftr
     Dx = cache_operator(Dx, x)
+    D2x = cache_operator(Dx * Dx, x)
 
-    u  = @. sin(5x)cos(7x);
-    du = @. 5cos(5x)cos(7x) - 7sin(5x)sin(7x);
+    u   = @. sin(5x)cos(7x);
+    du  = @. 5cos(5x)cos(7x) - 7sin(5x)sin(7x);
+    d2u = @. 5(-5sin(5x)cos(7x) -7cos(5x)sin(7x)) +
+           - 7(5cos(5x)sin(7x) + 7sin(5x)cos(7x))
 
     @test ≈(Dx * u, du; atol=1e-8)
+    @test ≈(D2x * u, d2u; atol=1e-8)
+
+    v = copy(u); @test ≈(mul!(v, D2x, u), d2u; atol=1e-8)
     v = copy(u); @test ≈(mul!(v, Dx, u), du; atol=1e-8)
 
     itr = inv(ftr)