From be847dd469179075550e29bdca374cebee5d70ff Mon Sep 17 00:00:00 2001 From: apkille Date: Wed, 7 Aug 2024 10:52:53 -0400 Subject: [PATCH 1/6] reorganize promote_u0 methods --- src/schroedinger.jl | 29 +++++++++++++++++++---------- 1 file changed, 19 insertions(+), 10 deletions(-) diff --git a/src/schroedinger.jl b/src/schroedinger.jl index 0910066d..2b785ca4 100644 --- a/src/schroedinger.jl +++ b/src/schroedinger.jl @@ -128,16 +128,25 @@ function _promote_time_and_state(u0, H::AbstractOperator, tspan) Ts = eltype(H) Tt = real(Ts) p = Vector{Tt}(undef,0) - u0data_promote = DiffEqBase.promote_u0(u0.data, p, tspan[1]) - tspan_promote = DiffEqBase.promote_tspan(u0data_promote, p, tspan, nothing, Dict{Symbol, Any}()) - if u0data_promote !== u0.data - u0_promote = rebuild(u0, u0data_promote) - return tspan_promote, u0_promote - end - return tspan_promote, u0 + u0_promote = DiffEqBase.promote_u0(u0, p, tspan[1]) + tspan_promote = DiffEqBase.promote_tspan(u0_promote, p, tspan, nothing, Dict{Symbol, Any}()) + return tspan_promote, u0_promote end _promote_time_and_state(u0, f, tspan) = _promote_time_and_state(u0, f(first(tspan), u0), tspan) -rebuild(op::Operator, new_data) = Operator(op.basis_l, op.basis_r, new_data) -rebuild(state::Ket, new_data) = Ket(state.basis, new_data) -rebuild(state::Bra, new_data) = Bra(state.basis, new_data) +@inline function DiffEqBase.promote_u0(u0::T, p, t0) where {T<:Union{Bra,Ket}} + u0data_promote = DiffEqBase.promote_u0(u0.data, p, t0) + if u0data_promote !== u0.data + u0_promote = T(u0.basis, u0data_promote) + return u0_promote + end + return u0 +end +@inline function DiffEqBase.promote_u0(u0::Operator, p, t0) + u0data_promote = DiffEqBase.promote_u0(u0.data, p, t0) + if u0data_promote !== u0.data + u0_promote = Operator(u0.basis_l, u0.basis_r, u0data_promote) + return u0_promote + end + return u0 +end From 247da7366f967c8e258dbae58cb0209b7ee55847 Mon Sep 17 00:00:00 2001 From: apkille Date: Fri, 9 Aug 2024 17:58:13 -0400 Subject: [PATCH 2/6] promote dual in master equations --- Project.toml | 1 + src/master.jl | 4 ++++ src/schroedinger.jl | 28 ------------------------- src/timeevolution_base.jl | 44 +++++++++++++++++++++++++++++++++++++++ test/test_ForwardDiff.jl | 33 +++++++++++++++++++++++++++++ 5 files changed, 82 insertions(+), 28 deletions(-) diff --git a/Project.toml b/Project.toml index c6080e37..c10469c0 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" diff --git a/src/master.jl b/src/master.jl index 0cf63ff6..e402d611 100644 --- a/src/master.jl +++ b/src/master.jl @@ -14,6 +14,7 @@ function master_h(tspan, rho0::Operator, H::AbstractOperator, J; _check_const.(J) _check_const.(Jdagger) check_master(rho0, H, J, Jdagger, rates) + tspan, rho0 = _promote_time_and_state(rho0, H, J, tspan) tmp = copy(rho0) dmaster_(t, rho, drho) = dmaster_h!(drho, H, J, Jdagger, rates, rho, tmp) integrate_master(tspan, dmaster_, rho0, fout; kwargs...) @@ -41,6 +42,7 @@ function master_nh(tspan, rho0::Operator, Hnh::AbstractOperator, J; _check_const.(J) _check_const.(Jdagger) check_master(rho0, Hnh, J, Jdagger, rates) + tspan, rho0 = _promote_time_and_state(rho0, Hnh, J, tspan) tmp = copy(rho0) dmaster_(t, rho, drho) = dmaster_nh!(drho, Hnh, Hnhdagger, J, Jdagger, rates, rho, tmp) integrate_master(tspan, dmaster_, rho0, fout; kwargs...) @@ -86,6 +88,7 @@ function master(tspan, rho0::Operator, H::AbstractOperator, J; _check_const(H) _check_const.(J) _check_const.(Jdagger) + tspan, rho0 = _promote_time_and_state(rho0, H, J, tspan) isreducible = check_master(rho0, H, J, Jdagger, rates) if !isreducible tmp = copy(rho0) @@ -124,6 +127,7 @@ function master(tspan, rho0::Operator, L::SuperOperator; fout=nothing, kwargs... b = GenericBasis(dim) rho_ = Ket(b,reshape(rho0.data, dim)) L_ = Operator(b,b,L.data) + tspan, rho_ = _promote_time_and_state(rho_, L_, tspan) dmaster_(t,rho,drho) = dmaster_liouville!(drho,L_,rho) # Rewrite into density matrix when saving diff --git a/src/schroedinger.jl b/src/schroedinger.jl index 2b785ca4..eddecc6c 100644 --- a/src/schroedinger.jl +++ b/src/schroedinger.jl @@ -122,31 +122,3 @@ function check_schroedinger(psi::Bra, H) check_multiplicable(psi, H) check_samebases(H) end - - -function _promote_time_and_state(u0, H::AbstractOperator, tspan) - Ts = eltype(H) - Tt = real(Ts) - p = Vector{Tt}(undef,0) - u0_promote = DiffEqBase.promote_u0(u0, p, tspan[1]) - tspan_promote = DiffEqBase.promote_tspan(u0_promote, p, tspan, nothing, Dict{Symbol, Any}()) - return tspan_promote, u0_promote -end -_promote_time_and_state(u0, f, tspan) = _promote_time_and_state(u0, f(first(tspan), u0), tspan) - -@inline function DiffEqBase.promote_u0(u0::T, p, t0) where {T<:Union{Bra,Ket}} - u0data_promote = DiffEqBase.promote_u0(u0.data, p, t0) - if u0data_promote !== u0.data - u0_promote = T(u0.basis, u0data_promote) - return u0_promote - end - return u0 -end -@inline function DiffEqBase.promote_u0(u0::Operator, p, t0) - u0data_promote = DiffEqBase.promote_u0(u0.data, p, t0) - if u0data_promote !== u0.data - u0_promote = Operator(u0.basis_l, u0.basis_r, u0data_promote) - return u0_promote - end - return u0 -end diff --git a/src/timeevolution_base.jl b/src/timeevolution_base.jl index 8949df0c..c34bb0d9 100644 --- a/src/timeevolution_base.jl +++ b/src/timeevolution_base.jl @@ -121,3 +121,47 @@ macro skiptimechecks(ex) end Base.@pure pure_inference(fout,T) = Core.Compiler.return_type(fout, T) + +function _promote_time_and_state(u0, H::AbstractOperator, tspan) + Ts = eltype(H) + Tt = real(Ts) + p = Vector{Tt}(undef,0) + u0_promote = DiffEqBase.promote_u0(u0, p, tspan[1]) + tspan_promote = DiffEqBase.promote_tspan(u0_promote.data, p, tspan, nothing, Dict{Symbol, Any}()) + return tspan_promote, u0_promote +end +function _promote_time_and_state(u0, H::AbstractOperator, J, tspan) + Ts = DiffEqBase.promote_dual(eltype(H), DiffEqBase.anyeltypedual(J)) + Tt = real(Ts) + p = Vector{Tt}(undef,0) + u0_promote = DiffEqBase.promote_u0(u0, p, tspan[1]) + tspan_promote = DiffEqBase.promote_tspan(u0_promote.data, p, tspan, nothing, Dict{Symbol, Any}()) + return tspan_promote, u0_promote +end + +_promote_time_and_state(u0, f, tspan) = _promote_time_and_state(u0, f(first(tspan)..., u0), tspan) + +@inline function DiffEqBase.promote_u0(u0::Ket, p, t0) + u0data_promote = DiffEqBase.promote_u0(u0.data, p, t0) + if u0data_promote !== u0.data + u0_promote = Ket(u0.basis, u0data_promote) + return u0_promote + end + return u0 +end +@inline function DiffEqBase.promote_u0(u0::Bra, p, t0) + u0data_promote = DiffEqBase.promote_u0(u0.data, p, t0) + if u0data_promote !== u0.data + u0_promote = Bra(u0.basis, u0data_promote) + return u0_promote + end + return u0 +end +@inline function DiffEqBase.promote_u0(u0::Operator, p, t0) + u0data_promote = DiffEqBase.promote_u0(u0.data, p, t0) + if u0data_promote !== u0.data + u0_promote = Operator(u0.basis_l, u0.basis_r, u0data_promote) + return u0_promote + end + return u0 +end \ No newline at end of file diff --git a/test/test_ForwardDiff.jl b/test/test_ForwardDiff.jl index 91b889fc..365282f2 100644 --- a/test/test_ForwardDiff.jl +++ b/test/test_ForwardDiff.jl @@ -1,6 +1,7 @@ using Test using OrdinaryDiffEq, QuantumOptics import ForwardDiff +import FiniteDiff # for some caese ForwardDiff.jl returns NaN due to issue with DiffEq.jl. see https://github.com/SciML/DiffEqBase.jl/issues/861 # Here we test; @@ -73,3 +74,35 @@ Ftdop(1.0) @test ForwardDiff.derivative(Ftdop, 1.0) isa Any end # testset + + +@testset "ForwardDiff with `master`" begin + +b = SpinBasis(1//2) +psi0 = spindown(b) +rho0 = dm(psi0) +params = [10.0, -3.0] + +# test to see if parameter propagates through Hamiltonian +H(p) = p[1]*sigmax(b) + p[2]*sigmam(b) # Hamiltonian +function cost1(p) # + tf, psif = timeevolution.master((0.0, pi), rho0, H(p), [sigmax(b)]) + return 1 - norm(psif) +end + +forwarddiff1 = ForwardDiff.gradient(cost1, params) +finitediff1 = FiniteDiff.finite_difference_gradient(cost1, params) +@test isapprox(forwarddiff1, finitediff1; atol=1e-5) + +# test to see if parameter propagates through Jump operator +J(p) = p[1]*sigmax(b) + p[2]*sigmam(b) # jump operator +function cost2(p) + tf, psif = timeevolution.master((0.0, pi), rho0, sigmax(b), [J(p)]) + return 1 - norm(psif) +end + +forwarddiff2 = ForwardDiff.gradient(cost2, params) +finitediff2 = FiniteDiff.finite_difference_gradient(cost2, params) +@test isapprox(forwarddiff2, finitediff2; atol=1e-5) + +end From 84855e99bb1b85e38821af4c1cce96d171a8e49d Mon Sep 17 00:00:00 2001 From: apkille Date: Fri, 9 Aug 2024 18:40:19 -0400 Subject: [PATCH 3/6] add more tests --- test/test_ForwardDiff.jl | 42 +++++++++++++++++++++------------------- 1 file changed, 22 insertions(+), 20 deletions(-) diff --git a/test/test_ForwardDiff.jl b/test/test_ForwardDiff.jl index 365282f2..4c43eb70 100644 --- a/test/test_ForwardDiff.jl +++ b/test/test_ForwardDiff.jl @@ -76,33 +76,35 @@ Ftdop(1.0) end # testset -@testset "ForwardDiff with `master`" begin +@testset "ForwardDiff with master" begin b = SpinBasis(1//2) psi0 = spindown(b) rho0 = dm(psi0) -params = [10.0, -3.0] +params = [rand(), rand()] + +for f in (:(timeevolution.master), :(timeevolution.master_h), :(timeevolution.master_nh)) + # test to see if parameter propagates through Hamiltonian + H(p) = p[1]*sigmax(b) + p[2]*sigmam(b) # Hamiltonian + function cost_H(p) # + tf, psif = eval(f)((0.0, pi), rho0, H(p), [sigmax(b)]) + return 1 - norm(psif) + end -# test to see if parameter propagates through Hamiltonian -H(p) = p[1]*sigmax(b) + p[2]*sigmam(b) # Hamiltonian -function cost1(p) # - tf, psif = timeevolution.master((0.0, pi), rho0, H(p), [sigmax(b)]) - return 1 - norm(psif) -end + forwarddiff_H = ForwardDiff.gradient(cost_H, params) + finitediff_H = FiniteDiff.finite_difference_gradient(cost_H, params) + @test isapprox(forwarddiff_H, finitediff_H; atol=1e-4) -forwarddiff1 = ForwardDiff.gradient(cost1, params) -finitediff1 = FiniteDiff.finite_difference_gradient(cost1, params) -@test isapprox(forwarddiff1, finitediff1; atol=1e-5) + # test to see if parameter propagates through Jump operator + J(p) = p[1]*sigmax(b) + p[2]*sigmam(b) # jump operator + function cost_J(p) + tf, psif = eval(f)((0.0, pi), rho0, sigmax(b), [J(p)]) + return 1 - norm(psif) + end -# test to see if parameter propagates through Jump operator -J(p) = p[1]*sigmax(b) + p[2]*sigmam(b) # jump operator -function cost2(p) - tf, psif = timeevolution.master((0.0, pi), rho0, sigmax(b), [J(p)]) - return 1 - norm(psif) + forwarddiff_J = ForwardDiff.gradient(cost_J, params) + finitediff_J = FiniteDiff.finite_difference_gradient(cost_J, params) + @test isapprox(forwarddiff_J, finitediff_J; atol=1e-4) end -forwarddiff2 = ForwardDiff.gradient(cost2, params) -finitediff2 = FiniteDiff.finite_difference_gradient(cost2, params) -@test isapprox(forwarddiff2, finitediff2; atol=1e-5) - end From 03067d9ebc0c9b09aa18829449a986bb3c61ddfa Mon Sep 17 00:00:00 2001 From: apkille Date: Sat, 10 Aug 2024 19:23:03 -0400 Subject: [PATCH 4/6] make FiniteDiff a test dependency --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index c10469c0..164bf2da 100644 --- a/Project.toml +++ b/Project.toml @@ -7,7 +7,6 @@ Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" -FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" @@ -40,10 +39,11 @@ WignerSymbols = "1, 2" julia = "1.3" [extras] +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["LinearAlgebra", "SparseArrays", "Random", "Test"] +test = ["FiniteDiff", "LinearAlgebra", "SparseArrays", "Random", "Test"] From 95cbe6b902eb2cbb0100455911a91dc1bcbf870d Mon Sep 17 00:00:00 2001 From: apkille Date: Sat, 10 Aug 2024 21:45:57 -0400 Subject: [PATCH 5/6] change atol to 1e-2 --- test/test_ForwardDiff.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_ForwardDiff.jl b/test/test_ForwardDiff.jl index 4c43eb70..7c50e63e 100644 --- a/test/test_ForwardDiff.jl +++ b/test/test_ForwardDiff.jl @@ -93,7 +93,7 @@ for f in (:(timeevolution.master), :(timeevolution.master_h), :(timeevolution.ma forwarddiff_H = ForwardDiff.gradient(cost_H, params) finitediff_H = FiniteDiff.finite_difference_gradient(cost_H, params) - @test isapprox(forwarddiff_H, finitediff_H; atol=1e-4) + @test isapprox(forwarddiff_H, finitediff_H; atol=1e-2) # test to see if parameter propagates through Jump operator J(p) = p[1]*sigmax(b) + p[2]*sigmam(b) # jump operator @@ -104,7 +104,7 @@ for f in (:(timeevolution.master), :(timeevolution.master_h), :(timeevolution.ma forwarddiff_J = ForwardDiff.gradient(cost_J, params) finitediff_J = FiniteDiff.finite_difference_gradient(cost_J, params) - @test isapprox(forwarddiff_J, finitediff_J; atol=1e-4) + @test isapprox(forwarddiff_J, finitediff_J; atol=1e-2) end end From 8f8c7fccc40f5728bda6137e828e71f54d852830 Mon Sep 17 00:00:00 2001 From: apkille Date: Thu, 15 Aug 2024 10:03:59 -0400 Subject: [PATCH 6/6] add ODE test for ForwardDiff --- test/test_ForwardDiff.jl | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/test/test_ForwardDiff.jl b/test/test_ForwardDiff.jl index 7c50e63e..e3ccba9c 100644 --- a/test/test_ForwardDiff.jl +++ b/test/test_ForwardDiff.jl @@ -13,6 +13,27 @@ import FiniteDiff # here we partially control the gradient error by limiting step size (dtmax) +@testset "ForwardDiff on ODE Problems" begin + +# schroedinger equation +b = SpinBasis(10//1) +psi0 = spindown(b) +H(p) = p[1]*sigmax(b) + p[2]*sigmam(b) +f_schrod!(dpsi, psi, p, t) = timeevolution.dschroedinger!(dpsi, H(p), psi) +function cost_schrod(p) + prob = ODEProblem(f_schrod!, psi0, (0.0, pi), p) + sol = solve(prob, DP5(); save_everystep=false) + return 1 - norm(sol[end]) +end + +p = [rand(), rand()] +fordiff_schrod = ForwardDiff.gradient(cost_schrod, p) +findiff_schrod = FiniteDiff.finite_difference_gradient(cost_schrod, p) + +@test isapprox(fordiff_schrod, findiff_schrod; atol=1e-2) + +end + @testset "ForwardDiff with schroedinger" begin # system