From 2f155584aede60c68db0e4a931c604a1951cdbaf Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Sun, 12 Mar 2023 22:07:29 -0400 Subject: [PATCH 1/4] fix nonlinear composition in composeoperator --- src/basic.jl | 22 ++++++++++++++++++++++ test/basic.jl | 18 ++++++++++++++++++ 2 files changed, 40 insertions(+) diff --git a/src/basic.jl b/src/basic.jl index cc93b92e..65f1751c 100644 --- a/src/basic.jl +++ b/src/basic.jl @@ -538,6 +538,28 @@ end #Base.:*(L::ComposedOperator, u::AbstractVecOrMat) = foldl((acc, op) -> op * acc, reverse(L.ops); init=u) #Base.:\(L::ComposedOperator, u::AbstractVecOrMat) = foldl((acc, op) -> op \ acc, L.ops; init=u) +function (L::ComposedOperator)(u, p, t) + v = u + for op in reverse(L.ops) + update_coefficients!(op, v, p, t) + v = op * v + end + + v +end + +function (L::ComposedOperator)(v, u, p, t) + @assert iscached(L) "cache needs to be set up for operator of type $(typeof(L)). + set up cache by calling cache_operator(L::AbstractSciMLOperator, u::AbstractArray)" + + vecs = (v, L.cache[1:end-1]..., u) + for i in reverse(1:length(L.ops)) + update_coefficients!(L.ops[i], vecs[i+1], p, t) + mul!(vecs[i], L.ops[i], vecs[i+1]) + end + v +end + function Base.:\(L::ComposedOperator, u::AbstractVecOrMat) v = u for op in L.ops diff --git a/test/basic.jl b/test/basic.jl index 5b29eb53..b2d74d6e 100644 --- a/test/basic.jl +++ b/test/basic.jl @@ -235,6 +235,24 @@ end @test ldiv!(rand(N), op, u) ≈ op \ u end +@testset "ComposedOperator nonlinear operator composition test" begin + u = rand(N) + p = nothing + t = 0.0 + + A = DiagonalOperator(zeros(N); update_func = (d, u, p, t) -> copy!(d, u)) # ret u .^2 + B = DiagonalOperator(zeros(N); update_func = (d, u, p, t) -> copy!(d, u)) + C = DiagonalOperator(zeros(N); update_func = (d, u, p, t) -> copy!(d, u)) + L = A ∘ B ∘ C + + ans = A(B(C(u, p, t), p, t), p, t) + + @test L(u, p, t) ≈ ans + + L = cache_operator(L, u) + v = rand(N); @test L(v, u, p, t) ≈ ans +end + @testset "Adjoint, Transpose" begin for (op, From b39765d96374f064db48795546171fb99dcbf3ec Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Mon, 13 Mar 2023 09:16:09 -0400 Subject: [PATCH 2/4] inv(::ComposeOp) not doing the correct thing. gotta think about how to propagate. maybe need more L(u, p, t) defs --- test/basic.jl | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/test/basic.jl b/test/basic.jl index b2d74d6e..cb1ccd66 100644 --- a/test/basic.jl +++ b/test/basic.jl @@ -251,6 +251,19 @@ end L = cache_operator(L, u) v = rand(N); @test L(v, u, p, t) ≈ ans + + Ai = inv(A) + Bi = inv(B) + Ci = inv(C) + + Li = inv(L) + for op in (Ai, Bi, Ci, Li) + @test op isa SciMLOperators.InvertedOperator + end + + ans = Ai(Bi(Ci(u, p, t), p, t), p, t) + @test_broken Li(u, p, t) ≈ ans + v = rand(N); @test_broken Li(v, u, p, t) ≈ ans end @testset "Adjoint, Transpose" begin From 616ccbf353bacca106e42cca24705d8b729269ce Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Fri, 5 May 2023 10:31:31 -0400 Subject: [PATCH 3/4] reorg --- src/basic.jl | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/basic.jl b/src/basic.jl index 65f1751c..88540765 100644 --- a/src/basic.jl +++ b/src/basic.jl @@ -560,25 +560,28 @@ function (L::ComposedOperator)(v, u, p, t) v end -function Base.:\(L::ComposedOperator, u::AbstractVecOrMat) +function Base.:*(L::ComposedOperator, u::AbstractVecOrMat) v = u - for op in L.ops - v = op \ v + for op in reverse(L.ops) + v = op * v end v end -function Base.:*(L::ComposedOperator, u::AbstractVecOrMat) +function Base.:\(L::ComposedOperator, u::AbstractVecOrMat) v = u - for op in reverse(L.ops) - v = op * v + for op in L.ops + v = op \ v end v end function cache_self(L::ComposedOperator, u::AbstractVecOrMat) + # TODO - use similar instead of *, \ + # similar(array, [element_type=eltype(array)], [dims=size(array)]) + # similar(u, Float32, (2,3)) if has_mul(L) vec = zero(u) cache = (vec,) From 5affa9e63ea129af1a946560457223fce7257da1 Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Fri, 5 May 2023 16:19:56 -0400 Subject: [PATCH 4/4] more todos --- src/basic.jl | 4 +--- test/basic.jl | 39 +++++++++++++++++++++++++++++++-------- 2 files changed, 32 insertions(+), 11 deletions(-) diff --git a/src/basic.jl b/src/basic.jl index 88540765..961f1d08 100644 --- a/src/basic.jl +++ b/src/basic.jl @@ -176,6 +176,7 @@ end ScaledOperator (λ L)*(u) = λ * L(u) + """ struct ScaledOperator{T, λType, @@ -579,9 +580,6 @@ function Base.:\(L::ComposedOperator, u::AbstractVecOrMat) end function cache_self(L::ComposedOperator, u::AbstractVecOrMat) - # TODO - use similar instead of *, \ - # similar(array, [element_type=eltype(array)], [dims=size(array)]) - # similar(u, Float32, (2,3)) if has_mul(L) vec = zero(u) cache = (vec,) diff --git a/test/basic.jl b/test/basic.jl index cb1ccd66..7b0d1052 100644 --- a/test/basic.jl +++ b/test/basic.jl @@ -240,30 +240,53 @@ end p = nothing t = 0.0 - A = DiagonalOperator(zeros(N); update_func = (d, u, p, t) -> copy!(d, u)) # ret u .^2 + square(u) = u .^ 2 + square(u, p, t) = u .^ 2 + square(v, u, p, t) = v .= u .* u + + root(u) = u .^ 2 + root(u, p, t) = u .^ 2 + root(v, u, p, t) = v .= u .* u + + F = FunctionOperator(square, u; islinear = false, op_inverse = root) + + A = DiagonalOperator(zeros(N); update_func = (d, u, p, t) -> copy!(d, u)) # u .^2 B = DiagonalOperator(zeros(N); update_func = (d, u, p, t) -> copy!(d, u)) C = DiagonalOperator(zeros(N); update_func = (d, u, p, t) -> copy!(d, u)) - L = A ∘ B ∘ C - ans = A(B(C(u, p, t), p, t), p, t) + L = A ∘ B ∘ C + F3 = F ∘ F ∘ F - @test L(u, p, t) ≈ ans + sq = u |> square |> square |> square + + @test A(B(C(u, p, t), p, t), p, t) ≈ sq + @test L(u, p, t) ≈ sq + @test F3(u, p, t) ≈ sq L = cache_operator(L, u) - v = rand(N); @test L(v, u, p, t) ≈ ans + v = rand(N); @test L(v, u, p, t) ≈ sq + + Fi = inv(F) + F3i = inv(F3) + + rt = u |> root |> root |> root + @test F3i(u, p, t) ≈ rt Ai = inv(A) Bi = inv(B) Ci = inv(C) Li = inv(L) + Fi = inv(F) for op in (Ai, Bi, Ci, Li) @test op isa SciMLOperators.InvertedOperator end - ans = Ai(Bi(Ci(u, p, t), p, t), p, t) - @test_broken Li(u, p, t) ≈ ans - v = rand(N); @test_broken Li(v, u, p, t) ≈ ans + rt = Ai(Bi(Ci(u, p, t), p, t), p, t) + @test Ai(u, p, t) ≈ ones(N) + # TODO - overwrite L(u, p, t) for InvertedOperator + @test_broken Li(u, p, t) ≈ ones(N) + v = rand(N); @test_broken Li(v, u, p, t) ≈ ones(N) end @testset "Adjoint, Transpose" begin