Merge pull request #141 from vpuri3/isconstant

create and export isconstant

Author: xtalax

Project.toml

@@ -1,7 +1,7 @@
 name = "SciMLOperators"
 uuid = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 authors = ["xtalax <[email protected]>"]
-version = "0.1.17"
+version = "0.1.18"
 
 [deps]
 ArrayInterfaceCore = "30b0a656-2188-435a-8636-2ec0e6a096e2"

src/SciMLOperators.jl

@@ -51,6 +51,7 @@ export ScalarOperator,
 export update_coefficients!,
        update_coefficients,
 
+       isconstant,
        iscached,
        cache_operator,
 

src/interface.jl

@@ -107,9 +107,18 @@ has_ldiv!(L::AbstractSciMLOperator) = false # ldiv!(du, L, u)
 
 ### Extra standard assumptions
 
-isconstant(L) = true
+isconstant(::Union{
+                   # LinearAlgebra
+                   AbstractMatrix,
+                   UniformScaling,
+                   Factorization,
+
+                   # Base
+                   Number,
+
+                  }
+          ) = true
 isconstant(L::AbstractSciMLOperator) = all(isconstant, getops(L))
-#isconstant(L::AbstractSciMLOperator) = L.update_func = DEFAULT_UPDATE_FUNC
 
 #islinear(L) = false
 islinear(::AbstractSciMLOperator) = false

src/matrix.jl

@@ -35,11 +35,13 @@ for op in (
            :adjoint,
            :transpose,
           )
-    @eval function Base.$op(L::MatrixOperator) # TODO - test this thoroughly
-        MatrixOperator(
-                       $op(L.A);
-                       update_func= (A,u,p,t) -> $op(L.update_func($op(L.A),u,p,t)) # TODO - test
-                      )
+    @eval function Base.$op(L::MatrixOperator)
+        if isconstant(L)
+            MatrixOperator($op(L.A))
+        else
+            update_func = (A,u,p,t) -> $op(L.update_func($op(L.A),u,p,t))
+            MatrixOperator($op(L.A); update_func = update_func)
+        end
     end
 end
 Base.conj(L::MatrixOperator) = MatrixOperator(
@@ -50,6 +52,7 @@ Base.conj(L::MatrixOperator) = MatrixOperator(
 has_adjoint(A::MatrixOperator) = has_adjoint(A.A)
 update_coefficients!(L::MatrixOperator,u,p,t) = (L.update_func(L.A,u,p,t); nothing)
 
+getops(L::MatrixOperator) = (L.A)
 isconstant(L::MatrixOperator) = L.update_func == DEFAULT_UPDATE_FUNC
 Base.iszero(L::MatrixOperator) = iszero(L.A)
 
@@ -76,8 +79,6 @@ Base.ndims(::Type{<:MatrixOperator{T,AType}}) where{T,AType} = ndims(AType)
 ArrayInterfaceCore.issingular(L::MatrixOperator) = ArrayInterfaceCore.issingular(L.A)
 Base.copy(L::MatrixOperator) = MatrixOperator(copy(L.A);update_func=L.update_func)
 
-getops(L::MatrixOperator) = (L.A)
-
 # operator application
 Base.:*(L::MatrixOperator, u::AbstractVecOrMat) = L.A * u
 Base.:\(L::MatrixOperator, u::AbstractVecOrMat) = L.A \ u
@@ -104,10 +105,11 @@ an operator of size `(N, N)` where `N = size(diag, 1)` is the leading length of
 `L` then is the elementwise-scaling operation on arrays of `length(u) = length(diag)`
 with leading length `size(u, 1) = N`.
 """
-function DiagonalOperator(diag::AbstractVector; update_func=DEFAULT_UPDATE_FUNC)
-    function diag_update_func(A, u, p, t)
-        update_func(A.diag, u, p, t)
-        A
+function DiagonalOperator(diag::AbstractVector; update_func = DEFAULT_UPDATE_FUNC)
+    diag_update_func = if update_func == DEFAULT_UPDATE_FUNC
+        DEFAULT_UPDATE_FUNC
+    else
+        (A, u, p, t) -> (update_func(A.diag, u, p, t); A)
     end
     MatrixOperator(Diagonal(diag); update_func=diag_update_func)
 end
@@ -214,7 +216,7 @@ struct AffineOperator{T,AType,BType,bType,cType,F} <: AbstractSciMLOperator{T}
     b::bType
 
     cache::cType
-    update_func::F
+    update_func::F # updates b
 
     function AffineOperator(A, B, b, cache, update_func)
         T = promote_type(eltype.((A,B,b))...)
@@ -234,7 +236,7 @@ end
 function AffineOperator(A::Union{AbstractMatrix,AbstractSciMLOperator},
                         B::Union{AbstractMatrix,AbstractSciMLOperator},
                         b::AbstractArray;
-                        update_func=DEFAULT_UPDATE_FUNC,
+                        update_func = DEFAULT_UPDATE_FUNC,
                        )
     @assert size(A, 1) == size(B, 1) "Dimension mismatch: A, B don't output vectors
     of same size"
@@ -250,7 +252,7 @@ end
     L = AddVector(b[; update_func])
     L(u) = u + b
 """
-function AddVector(b::AbstractVecOrMat; update_func=DEFAULT_UPDATE_FUNC)
+function AddVector(b::AbstractVecOrMat; update_func = DEFAULT_UPDATE_FUNC)
     N  = size(b, 1)
     Id = IdentityOperator{N}()
 
@@ -261,19 +263,20 @@ end
     L = AddVector(B, b[; update_func])
     L(u) = u + B*b
 """
-function AddVector(B, b::AbstractVecOrMat; update_func=DEFAULT_UPDATE_FUNC)
+function AddVector(B, b::AbstractVecOrMat; update_func = DEFAULT_UPDATE_FUNC)
     N = size(B, 1)
     Id = IdentityOperator{N}()
 
     AffineOperator(Id, B, b; update_func=update_func)
 end
 
 getops(L::AffineOperator) = (L.A, L.B, L.b)
-Base.size(L::AffineOperator) = size(L.A)
 
 update_coefficients!(L::AffineOperator,u,p,t) = (L.update_func(L.b,u,p,t); nothing)
-
+isconstant(L::AffineOperator) = (L.update_func == DEFAULT_UPDATE_FUNC) & all(isconstant, (L.A, L.B))
 islinear(::AffineOperator) = false
+
+Base.size(L::AffineOperator) = size(L.A)
 Base.iszero(L::AffineOperator) = all(iszero, getops(L))
 has_adjoint(L::AffineOperator) = all(has_adjoint, L.ops)
 has_mul(L::AffineOperator) = has_mul(L.A)

test/basic.jl

@@ -36,6 +36,7 @@ K = 12
     @test iscached(Id)
     @test size(Id) == (N, N)
     @test Id' isa IdentityOperator{N}
+    @test isconstant(Id)
 
     for op in (
                *, \,
@@ -67,6 +68,7 @@ end
     @test issquare(Z)
     @test islinear(Z)
     @test NullOperator(u) isa NullOperator{N}
+    @test isconstant(Z)
     @test zero(A) isa NullOperator{N}
     @test convert(AbstractMatrix, Z) == zeros(size(Z))
 
@@ -105,6 +107,7 @@ end
     op = ScaledOperator(α, MatrixOperator(A))
 
     @test op isa ScaledOperator
+    @test isconstant(op)
     @test iscached(op)
     @test issquare(op)
     @test islinear(op)
@@ -115,6 +118,7 @@ end
     opF = factorize(op)
 
     @test opF isa ScaledOperator
+    @test isconstant(opF)
     @test iscached(opF)
 
     @test α * A  ≈ convert(AbstractMatrix, op) ≈ convert(AbstractMatrix, opF)
@@ -148,6 +152,11 @@ end
         @test op3 isa AddedOperator
         @test op4 isa AddedOperator
 
+        @test isconstant(op1)
+        @test isconstant(op2)
+        @test isconstant(op3)
+        @test isconstant(op4)
+
         @test op1 * u ≈ op(  A*u,   B*u)
         @test op2 * u ≈ op(α*A*u,   B*u)
         @test op3 * u ≈ op(  A*u, β*B*u)
@@ -175,13 +184,16 @@ end
     op = ∘(MatrixOperator.((A, B, C))...)
 
     @test op isa ComposedOperator
+    @test isconstant(op)
+
     @test *(op.ops...) isa ComposedOperator
     @test issquare(op)
     @test islinear(op)
 
     opF = factorize(op)
 
     @test opF isa ComposedOperator
+    @test isconstant(opF)
     @test issquare(opF)
     @test islinear(opF)
 
@@ -256,6 +268,9 @@ end
         AAt = LType(AA)
         DDt = LType(DD)
 
+        @test isconstant(AAt)
+        @test isconstant(DDt)
+
         @test AAt.L === AA
         @test op(u) isa VType
 
@@ -280,6 +295,7 @@ end
 
     @test !iscached(Di)
     Di = cache_operator(Di, u)
+    @test isconstant(Di)
     @test iscached(Di)
 
     @test issquare(Di)

test/matrix.jl

@@ -24,6 +24,9 @@ K = 19
     @test AA  isa MatrixOperator
     @test AAt isa MatrixOperator
 
+    @test isconstant(AA)
+    @test isconstant(AAt)
+
     @test issquare(AA)
     @test islinear(AA)
 
@@ -33,6 +36,9 @@ K = 19
     @test FF  isa InvertibleOperator
     @test FFt isa InvertibleOperator
 
+    @test isconstant(FF)
+    @test isconstant(FFt)
+
     @test eachindex(A)  === eachindex(AA)
     @test eachindex(A') === eachindex(AAt) === eachindex(MatrixOperator(At))
 
@@ -58,9 +64,11 @@ end
     t = rand()
 
     L = MatrixOperator(zeros(N,N);
-                       update_func= (A,u,p,t) -> (A .= p*p'; nothing)
+                       update_func = (A,u,p,t) -> (A .= p*p'; nothing)
                       )
 
+    @test !isconstant(L)
+
     A = p*p'
     ans = A * u
     @test L(u,p,t) ≈ ans
@@ -73,9 +81,10 @@ end
     t = rand()
 
     D = DiagonalOperator(zeros(N);
-                         update_func= (diag,u,p,t) -> (diag .= p*t; nothing)
+                         update_func = (diag,u,p,t) -> (diag .= p*t; nothing)
                         )
 
+    @test !isconstant(D)
     @test issquare(D)
     @test islinear(D)
 
@@ -91,6 +100,7 @@ end
     β = rand()
 
     L = DiagonalOperator(d)
+    @test isconstant(L)
 
     @test issquare(L)
     @test islinear(L)
@@ -114,6 +124,7 @@ end
     β = rand()
 
     L = AffineOperator(MatrixOperator(A), MatrixOperator(B), b)
+    @test isconstant(L)
     @test issquare(L)
     @test !islinear(L)
 
@@ -161,9 +172,11 @@ end
     t = rand()
 
     L = AffineOperator(A, B, b;
-                       update_func= (b,u,p,t) -> (b .= Diagonal(p*t)*b; nothing)
+                       update_func = (b,u,p,t) -> (b .= Diagonal(p*t)*b; nothing)
                       )
 
+    @test !isconstant(L)
+
     b = Diagonal(p*t)*b
     ans = A * u + B * b
     @test L(u,p,t) ≈ ans
@@ -208,6 +221,12 @@ for square in [false, true] #for K in [1, K]
     @test opAB  isa TensorProductOperator
     @test opABC isa TensorProductOperator
 
+    @test isconstant(opAB)
+    @test isconstant(opABC)
+
+    @test islinear(opAB)
+    @test islinear(opABC)
+
     if square
         @test issquare(opAB)
         @test issquare(opABC)
@@ -216,16 +235,16 @@ for square in [false, true] #for K in [1, K]
         @test !issquare(opABC)
     end
 
-    @test islinear(opAB)
-    @test islinear(opABC)
-
     @test AB ≈ convert(AbstractMatrix, opAB)
     @test ABC ≈ convert(AbstractMatrix, opABC)
 
     # factorization tests
     opAB_F = factorize(opAB)
     opABC_F = factorize(opABC)
 
+    @test isconstant(opAB_F)
+    @test isconstant(opABC_F)
+
     @test opAB_F isa TensorProductOperator
     @test opABC_F isa TensorProductOperator
 

test/scalar.jl

@@ -22,6 +22,7 @@ K = 12
     @test convert(ScalarOperator, a) isa ScalarOperator
 
     @test size(α) == ()
+    @test isconstant(α)
 
     v=copy(u); @test lmul!(α, u) ≈ v * x
     v=copy(u); @test rmul!(u, α) ≈ x * v
@@ -68,6 +69,9 @@ end
     α = ScalarOperator(0.0; update_func=(a,u,p,t) -> p)
     β = ScalarOperator(0.0; update_func=(a,u,p,t) -> t)
 
+    @test !isconstant(α)
+    @test !isconstant(β)
+
     @test α(u,p,t)   ≈ p * u
     @test α(v,u,p,t) ≈ p * u