diff --git a/ext/DiffEqBaseForwardDiffExt.jl b/ext/DiffEqBaseForwardDiffExt.jl index adfb044de..3f3fbc2ee 100644 --- a/ext/DiffEqBaseForwardDiffExt.jl +++ b/ext/DiffEqBaseForwardDiffExt.jl @@ -2,8 +2,9 @@ module DiffEqBaseForwardDiffExt using DiffEqBase, ForwardDiff using DiffEqBase.ArrayInterface -using DiffEqBase: Void, FunctionWrappersWrappers, OrdinaryDiffEqTag, AbstractTimeseriesSolution, - RecursiveArrayTools, reduce_tup, _promote_tspan, has_continuous_callback +using DiffEqBase: Void, FunctionWrappersWrappers, OrdinaryDiffEqTag, + AbstractTimeseriesSolution, + RecursiveArrayTools, reduce_tup, _promote_tspan, has_continuous_callback import DiffEqBase: hasdualpromote, wrapfun_oop, wrapfun_iip, prob2dtmin, promote_tspan, anyeltypedual, isdualtype, value, ODE_DEFAULT_NORM, InternalITP, nextfloat_tdir, DualEltypeChecker, sse @@ -502,7 +503,8 @@ unitfulvalue(x::ForwardDiff.Dual) = unitfulvalue(ForwardDiff.unitfulvalue(x)) sse(x::ForwardDiff.Dual) = sse(ForwardDiff.value(x)) + sum(sse, ForwardDiff.partials(x)) function DiffEqBase.totallength(x::ForwardDiff.Dual) - return DiffEqBase.totallength(ForwardDiff.value(x)) + sum(DiffEqBase.totallength, ForwardDiff.partials(x)) + return DiffEqBase.totallength(ForwardDiff.value(x)) + + sum(DiffEqBase.totallength, ForwardDiff.partials(x)) end @inline ODE_DEFAULT_NORM(u::ForwardDiff.Dual, ::Any) = sqrt(sse(u)) diff --git a/ext/DiffEqBaseGTPSAExt.jl b/ext/DiffEqBaseGTPSAExt.jl index 8214a48e9..65c383f73 100644 --- a/ext/DiffEqBaseGTPSAExt.jl +++ b/ext/DiffEqBaseGTPSAExt.jl @@ -35,4 +35,4 @@ function ODE_DEFAULT_NORM(f::F, u::AbstractArray{<:TPS}, t) where {F} Base.FastMath.sqrt_fast(x / max(length(u), 1)) end -end \ No newline at end of file +end diff --git a/ext/DiffEqBaseMonteCarloMeasurementsExt.jl b/ext/DiffEqBaseMonteCarloMeasurementsExt.jl index 2ddb724d8..2a5d0ed84 100644 --- a/ext/DiffEqBaseMonteCarloMeasurementsExt.jl +++ b/ext/DiffEqBaseMonteCarloMeasurementsExt.jl @@ -26,7 +26,10 @@ end DiffEqBase.value(x::Type{MonteCarloMeasurements.AbstractParticles{T, N}}) where {T, N} = T DiffEqBase.value(x::MonteCarloMeasurements.AbstractParticles) = mean(x.particles) -DiffEqBase.unitfulvalue(x::Type{MonteCarloMeasurements.AbstractParticles{T, N}}) where {T, N} = T +function DiffEqBase.unitfulvalue(x::Type{MonteCarloMeasurements.AbstractParticles{ + T, N}}) where {T, N} + T +end DiffEqBase.unitfulvalue(x::MonteCarloMeasurements.AbstractParticles) = mean(x.particles) # Support adaptive steps should be errorless diff --git a/ext/DiffEqBaseReverseDiffExt.jl b/ext/DiffEqBaseReverseDiffExt.jl index 872d519a1..879989d30 100644 --- a/ext/DiffEqBaseReverseDiffExt.jl +++ b/ext/DiffEqBaseReverseDiffExt.jl @@ -22,7 +22,6 @@ end DiffEqBase.value(x::ReverseDiff.TrackedReal) = x.value DiffEqBase.value(x::ReverseDiff.TrackedArray) = x.value - DiffEqBase.unitfulvalue(x::Type{ReverseDiff.TrackedReal{V, D, O}}) where {V, D, O} = V function DiffEqBase.unitfulvalue(x::Type{ ReverseDiff.TrackedArray{V, D, N, VA, DA}, diff --git a/ext/DiffEqBaseTrackerExt.jl b/ext/DiffEqBaseTrackerExt.jl index 406927ded..6d77714ba 100644 --- a/ext/DiffEqBaseTrackerExt.jl +++ b/ext/DiffEqBaseTrackerExt.jl @@ -16,7 +16,9 @@ DiffEqBase.value(x::Tracker.TrackedReal) = x.data DiffEqBase.value(x::Tracker.TrackedArray) = x.data DiffEqBase.unitfulvalue(x::Type{Tracker.TrackedReal{T}}) where {T} = T -DiffEqBase.unitfulvalue(x::Type{Tracker.TrackedArray{T, N, A}}) where {T, N, A} = Array{T, N} +function DiffEqBase.unitfulvalue(x::Type{Tracker.TrackedArray{T, N, A}}) where {T, N, A} + Array{T, N} +end DiffEqBase.unitfulvalue(x::Tracker.TrackedReal) = x.data DiffEqBase.unitfulvalue(x::Tracker.TrackedArray) = x.data diff --git a/src/norecompile.jl b/src/norecompile.jl index 273fb14ce..d20d22a76 100644 --- a/src/norecompile.jl +++ b/src/norecompile.jl @@ -23,9 +23,11 @@ end # Default dispatch assumes no ForwardDiff, gets added in the new dispatch function wrapfun_iip(ff, inputs) - FunctionWrappersWrappers.FunctionWrappersWrapper(Void(ff), (typeof(inputs),), (Nothing,)) + FunctionWrappersWrappers.FunctionWrappersWrapper( + Void(ff), (typeof(inputs),), (Nothing,)) end function wrapfun_oop(ff, inputs) - FunctionWrappersWrappers.FunctionWrappersWrapper(ff, (typeof(inputs),), (typeof(inputs[1]),)) -end \ No newline at end of file + FunctionWrappersWrappers.FunctionWrappersWrapper( + ff, (typeof(inputs),), (typeof(inputs[1]),)) +end diff --git a/test/downstream/gtpsa.jl b/test/downstream/gtpsa.jl index 496ebc10d..6b506a0e7 100644 --- a/test/downstream/gtpsa.jl +++ b/test/downstream/gtpsa.jl @@ -9,68 +9,69 @@ x = [1.0, 2.0, 3.0] p = [4.0, 5.0, 6.0] prob = ODEProblem(f!, x, (0.0, 1.0), p) -sol = solve(prob, Tsit5(), reltol=1e-16, abstol=1e-16) +sol = solve(prob, Tsit5(), reltol = 1e-16, abstol = 1e-16) # Parametric GTPSA map desc = Descriptor(3, 2, 3, 2) # 3 variables 3 parameters, both to 2nd order dx = @vars(desc) dp = @params(desc) prob_GTPSA = ODEProblem(f!, x .+ dx, (0.0, 1.0), p .+ dp) -sol_GTPSA = solve(prob_GTPSA, Tsit5(), reltol=1e-16, abstol=1e-16) +sol_GTPSA = solve(prob_GTPSA, Tsit5(), reltol = 1e-16, abstol = 1e-16) @test sol.u[end] ≈ scalar.(sol_GTPSA.u[end]) # scalar gets 0th order part # Compare Jacobian against ForwardDiff J_FD = ForwardDiff.jacobian([x..., p...]) do t prob = ODEProblem(f!, t[1:3], (0.0, 1.0), t[4:6]) - sol = solve(prob, Tsit5(), reltol=1e-16, abstol=1e-16) + sol = solve(prob, Tsit5(), reltol = 1e-16, abstol = 1e-16) sol.u[end] end -@test J_FD ≈ GTPSA.jacobian(sol_GTPSA.u[end], include_params=true) +@test J_FD ≈ GTPSA.jacobian(sol_GTPSA.u[end], include_params = true) # Compare Hessians against ForwardDiff for i in 1:3 Hi_FD = ForwardDiff.hessian([x..., p...]) do t prob = ODEProblem(f!, t[1:3], (0.0, 1.0), t[4:6]) - sol = solve(prob, Tsit5(), reltol=1e-16, abstol=1e-16) + sol = solve(prob, Tsit5(), reltol = 1e-16, abstol = 1e-16) sol.u[end][i] end - @test Hi_FD ≈ GTPSA.hessian(sol_GTPSA.u[end][i], include_params=true) + @test Hi_FD ≈ GTPSA.hessian(sol_GTPSA.u[end][i], include_params = true) end - # ODEProblem 2 ======================= -pdot!(dq, p, q, params, t) = dq .= [0.0, 0.0, 0.0] -qdot!(dp, p, q, params, t) = dp .= [p[1] / sqrt((1 + p[3])^2 - p[1]^2 - p[2]^2), - p[2] / sqrt((1 + p[3])^2 - p[1]^2 - p[2]^2), - p[3] / sqrt(1 + p[3]^2) - (p[3] + 1)/sqrt((1 + p[3])^2 - p[1]^2 - p[2]^2)] +pdot!(dq, p, q, params, t) = dq .= [0.0, 0.0, 0.0] +function qdot!(dp, p, q, params, t) + dp .= [p[1] / sqrt((1 + p[3])^2 - p[1]^2 - p[2]^2), + p[2] / sqrt((1 + p[3])^2 - p[1]^2 - p[2]^2), + p[3] / sqrt(1 + p[3]^2) - (p[3] + 1) / sqrt((1 + p[3])^2 - p[1]^2 - p[2]^2)] +end prob = DynamicalODEProblem(pdot!, qdot!, [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], (0.0, 25.0)) -sol = solve(prob, Yoshida6(), dt = 1.0, reltol=1e-16, abstol=1e-16) +sol = solve(prob, Yoshida6(), dt = 1.0, reltol = 1e-16, abstol = 1e-16) desc = Descriptor(6, 2) # 6 variables to 2nd order -dx = @vars(desc) # identity map +dx = @vars(desc) # identity map prob_GTPSA = DynamicalODEProblem(pdot!, qdot!, dx[1:3], dx[4:6], (0.0, 25.0)) -sol_GTPSA = solve(prob_GTPSA, Yoshida6(), dt = 1.0, reltol=1e-16, abstol=1e-16) +sol_GTPSA = solve(prob_GTPSA, Yoshida6(), dt = 1.0, reltol = 1e-16, abstol = 1e-16) @test sol.u[end] ≈ scalar.(sol_GTPSA.u[end]) # scalar gets 0th order part # Compare Jacobian against ForwardDiff J_FD = ForwardDiff.jacobian(zeros(6)) do t prob = DynamicalODEProblem(pdot!, qdot!, t[1:3], t[4:6], (0.0, 25.0)) - sol = solve(prob, Yoshida6(), dt = 1.0, reltol=1e-16, abstol=1e-16) + sol = solve(prob, Yoshida6(), dt = 1.0, reltol = 1e-16, abstol = 1e-16) sol.u[end] end -@test J_FD ≈ GTPSA.jacobian(sol_GTPSA.u[end], include_params=true) +@test J_FD ≈ GTPSA.jacobian(sol_GTPSA.u[end], include_params = true) # Compare Hessians against ForwardDiff for i in 1:6 Hi_FD = ForwardDiff.hessian(zeros(6)) do t - prob = DynamicalODEProblem(pdot!, qdot!, t[1:3], t[4:6], (0.0, 25.0)) - sol = solve(prob, Yoshida6(), dt = 1.0, reltol=1e-16, abstol=1e-16) + prob = DynamicalODEProblem(pdot!, qdot!, t[1:3], t[4:6], (0.0, 25.0)) + sol = solve(prob, Yoshida6(), dt = 1.0, reltol = 1e-16, abstol = 1e-16) sol.u[end][i] end - @test Hi_FD ≈ GTPSA.hessian(sol_GTPSA.u[end][i], include_params=true) + @test Hi_FD ≈ GTPSA.hessian(sol_GTPSA.u[end][i], include_params = true) end diff --git a/test/forwarddiff_dual_detection.jl b/test/forwarddiff_dual_detection.jl index 3685164ea..1fc50fefd 100644 --- a/test/forwarddiff_dual_detection.jl +++ b/test/forwarddiff_dual_detection.jl @@ -387,5 +387,6 @@ u = Dual.(val, par) @test DiffEqBase.totallength(val[1]) == 1 @test DiffEqBase.totallength(val) == length(val) @test DiffEqBase.totallength(par) == length(par) -@test DiffEqBase.totallength(u[1]) == DiffEqBase.totallength(val[1]) + DiffEqBase.totallength(par[1]) +@test DiffEqBase.totallength(u[1]) == + DiffEqBase.totallength(val[1]) + DiffEqBase.totallength(par[1]) @test DiffEqBase.totallength(u) == sum(DiffEqBase.totallength, u)