From 153b0dca82df125401aa4aa83b9ea59786428803 Mon Sep 17 00:00:00 2001 From: Jonas Koziorek Date: Sat, 13 Jul 2024 22:59:58 +0200 Subject: [PATCH 1/7] jacobian code WIP --- src/core/jacobian.jl | 56 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) create mode 100644 src/core/jacobian.jl diff --git a/src/core/jacobian.jl b/src/core/jacobian.jl new file mode 100644 index 00000000..5787c09f --- /dev/null +++ b/src/core/jacobian.jl @@ -0,0 +1,56 @@ +# TODO: Where to put this file? + +export jacobian + +import ForwardDiff + +function jacobian(ds::CoreDynamicalSystem{IIP}, J0=nothing) where {IIP} + f = dynamic_rule(ds) + u0 = current_state(ds) + _jacobian(f, J0, Val{IIP}(), u0) + +end + +function _jacobian(f, J0, ::Val{true}, u0) + cfg = ForwardDiff.JacobianConfig( + (du, u) -> f(du, u, p, p), deepcopy(u0), deepcopy(u0) + ) + Jf = (du, u, p, t) -> begin + uv = @view u[:, 1] + ForwardDiff.jacobian!( + J0, (du, u) -> f(du, u, p, t), view(du, :, 1), uv, cfg, Val{false}() + ) + nothing + end + return Jf +end + +function _jacobian(f, J0, ::Val{false}, u0) + # Initial matrix `J0` is ignored + Jf = (u, p, t) -> ForwardDiff.jacobian((x) -> f(x, p, t), u) + return Jf +end + + +function oop(u, p, t) + return p[1] * SVector(u[1], -u[2]) +end + +function iip(du, u, p, t) + du[1] = p[1] * u[1] + du[2] = -p[1] * u[2] + return nothing +end + +#%% +oopds=CoupledODEs(oop, [1.0, 1.0], [1.0]) +oopjac=jacobian(oopds) +oopjac(current_state(oopds), current_parameters(oopds), 0.0) + +iipds=CoupledODEs(iip, [1.0, 1.0], [1.0]) +J0 = zeros(2, 2) +iipjac=jacobian(iipds, J0) +du = zeros(2) +iipjac(du, current_state(iipds), current_parameters(iipds), 0.0) +J0 + From ad27eacf7b895405f1788739fa14730fdbcce949 Mon Sep 17 00:00:00 2001 From: Jonas Koziorek Date: Sat, 20 Jul 2024 12:59:55 +0200 Subject: [PATCH 2/7] Separating Jacobian from tangent_rule --- src/DynamicalSystemsBase.jl | 1 + src/core/jacobian.jl | 56 ------------------------------------ src/core_systems/jacobian.jl | 30 +++++++++++++++++++ 3 files changed, 31 insertions(+), 56 deletions(-) delete mode 100644 src/core/jacobian.jl create mode 100644 src/core_systems/jacobian.jl diff --git a/src/DynamicalSystemsBase.jl b/src/DynamicalSystemsBase.jl index fb215f98..07c4474f 100644 --- a/src/DynamicalSystemsBase.jl +++ b/src/DynamicalSystemsBase.jl @@ -24,6 +24,7 @@ include("core_systems/discrete_time_map.jl") include("core_systems/continuous_time_ode.jl") include("core_systems/arbitrary_steppable.jl") include("core_systems/additional_supertypes.jl") +include("core_systems/jacobian.jl") include("derived_systems/stroboscopic_map.jl") include("derived_systems/parallel_systems.jl") diff --git a/src/core/jacobian.jl b/src/core/jacobian.jl deleted file mode 100644 index 5787c09f..00000000 --- a/src/core/jacobian.jl +++ /dev/null @@ -1,56 +0,0 @@ -# TODO: Where to put this file? - -export jacobian - -import ForwardDiff - -function jacobian(ds::CoreDynamicalSystem{IIP}, J0=nothing) where {IIP} - f = dynamic_rule(ds) - u0 = current_state(ds) - _jacobian(f, J0, Val{IIP}(), u0) - -end - -function _jacobian(f, J0, ::Val{true}, u0) - cfg = ForwardDiff.JacobianConfig( - (du, u) -> f(du, u, p, p), deepcopy(u0), deepcopy(u0) - ) - Jf = (du, u, p, t) -> begin - uv = @view u[:, 1] - ForwardDiff.jacobian!( - J0, (du, u) -> f(du, u, p, t), view(du, :, 1), uv, cfg, Val{false}() - ) - nothing - end - return Jf -end - -function _jacobian(f, J0, ::Val{false}, u0) - # Initial matrix `J0` is ignored - Jf = (u, p, t) -> ForwardDiff.jacobian((x) -> f(x, p, t), u) - return Jf -end - - -function oop(u, p, t) - return p[1] * SVector(u[1], -u[2]) -end - -function iip(du, u, p, t) - du[1] = p[1] * u[1] - du[2] = -p[1] * u[2] - return nothing -end - -#%% -oopds=CoupledODEs(oop, [1.0, 1.0], [1.0]) -oopjac=jacobian(oopds) -oopjac(current_state(oopds), current_parameters(oopds), 0.0) - -iipds=CoupledODEs(iip, [1.0, 1.0], [1.0]) -J0 = zeros(2, 2) -iipjac=jacobian(iipds, J0) -du = zeros(2) -iipjac(du, current_state(iipds), current_parameters(iipds), 0.0) -J0 - diff --git a/src/core_systems/jacobian.jl b/src/core_systems/jacobian.jl new file mode 100644 index 00000000..873a7053 --- /dev/null +++ b/src/core_systems/jacobian.jl @@ -0,0 +1,30 @@ +export jacobian + +import ForwardDiff + +function jacobian(ds::CoreDynamicalSystem{IIP}) where {IIP} + _jacobian(ds, Val{IIP}()) +end + +function _jacobian(ds, ::Val{true}) + f = dynamic_rule(ds) + u0 = current_state(ds) + cfg = ForwardDiff.JacobianConfig( + (du, u) -> f(du, u, p, p), deepcopy(u0), deepcopy(u0) + ) + Jf = (J0, u, p, t) -> begin + uv = @view u[:, 1] + du = copy(u) + ForwardDiff.jacobian!( + J0, (du, u) -> f(du, u, p, t), view(du, :, 1), uv, cfg, Val{false}() + ) + nothing + end + return Jf +end + +function _jacobian(ds, ::Val{false}) + f = dynamic_rule(ds) + Jf = (u, p, t) -> ForwardDiff.jacobian((x) -> f(x, p, t), u) + return Jf +end \ No newline at end of file From 47da9ea329f61078acb5617ed7f0f953a7098cc6 Mon Sep 17 00:00:00 2001 From: Jonas Koziorek Date: Sat, 20 Jul 2024 13:00:20 +0200 Subject: [PATCH 3/7] inserting new jacobian back into tangent_rule --- src/derived_systems/tangent_space.jl | 28 ++-------------------------- 1 file changed, 2 insertions(+), 26 deletions(-) diff --git a/src/derived_systems/tangent_space.jl b/src/derived_systems/tangent_space.jl index 6cdb8d25..116f514e 100644 --- a/src/derived_systems/tangent_space.jl +++ b/src/derived_systems/tangent_space.jl @@ -113,6 +113,7 @@ function TangentDynamicalSystem(ds::CoreDynamicalSystem{IIP}; u0_correct = correct_state(Val{IIP}(), u0) Q0_correct = correct_matrix_type(Val{IIP}(), Q0) newstate = hcat(u0_correct, Q0_correct) + J = isnothing(J) ? jacobian(ds) : J newrule = tangent_rule(f, J, J0, Val{IIP}(), Val{k}(), u0_correct) # Pass everything to analytic system constructors @@ -149,37 +150,12 @@ function tangent_rule(f::F, J::JAC, J0, ::Val{true}, ::Val{k}, u0) where {F, JAC end return tangentf end -# for the case of autodiffed systems, a specialized version is created -# so that f! is not called twice in ForwardDiff -function tangent_rule(f::F, ::Nothing, J0, ::Val{true}, ::Val{k}, u0) where {F, k} - let - cfg = ForwardDiff.JacobianConfig( - (du, u) -> f(du, u, p, p), deepcopy(u0), deepcopy(u0) - ) - tangentf = (du, u, p, t) -> begin - uv = @view u[:, 1] - ForwardDiff.jacobian!( - J0, (du, u) -> f(du, u, p, t), view(du, :, 1), uv, cfg, Val{false}() - ) - mul!((@view du[:, 2:k+1]), J0, (@view u[:, 2:k+1])) - nothing - end - return tangentf - end -end # OOP Tangent space dynamics function tangent_rule(f::F, J::JAC, J0, ::Val{false}, ::Val{k}, u0) where {F, JAC, k} - # out of place - if JAC == Nothing - # There is no config needed here - Jf = (u, p, t) -> ForwardDiff.jacobian((x) -> f(x, p, t), u) - else - Jf = J - end # Initial matrix `J0` is ignored ws_index = SVector{k, Int}(2:(k+1)...) - tangentf = TangentOOP(f, Jf, ws_index) + tangentf = TangentOOP(f, J, ws_index) return tangentf end struct TangentOOP{F, JAC, k} <: Function From 4d0f95e3a3f71fcf3a05e6536aea8c16e9e4afa9 Mon Sep 17 00:00:00 2001 From: Jonas Koziorek Date: Sat, 20 Jul 2024 13:00:30 +0200 Subject: [PATCH 4/7] tests --- test/jacobian.jl | 29 +++++++++++++++++++++++++++++ test/runtests.jl | 1 + 2 files changed, 30 insertions(+) create mode 100644 test/jacobian.jl diff --git a/test/jacobian.jl b/test/jacobian.jl new file mode 100644 index 00000000..a36869d8 --- /dev/null +++ b/test/jacobian.jl @@ -0,0 +1,29 @@ +using DynamicalSystemsBase, Test + +function oop(u, p, t) + return p[1] * SVector(u[1], -u[2]) +end + +function iip(du, u, p, t) + du .= oop(u, p, t) + return nothing +end + +#%% +@testset "IDT=$(IDT), IIP=$(IIP)" for IDT in (true, false), IIP in (false, true) + SystemType = IDT ? DeterministicIteratedMap : CoupledODEs + rule = IIP ? iip : oop + p = 3.0 + u0 = [1.0, 1.0] + result = [p 0.0; 0.0 -p] + + ds = SystemType(rule, u0, p) + J0 = zeros(dimension(ds), dimension(ds)) + J = jacobian(ds) + if IIP + J(J0, current_state(ds), current_parameters(ds), 0.0) + @test J0 == result + else + @test J(current_state(ds), current_parameters(ds), 0.0) == result + end +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index ec852aec..3084f88f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,4 +16,5 @@ testfile(file, testname=defaultname(file)) = @testset "$testname" begin; include testfile("successful_step.jl") testfile("mtk_integration.jl") testfile("trajectory.jl") + testfile("jacobian.jl") end From 2be772ea95dff37931bf5a9692dfcabb2ed7c87c Mon Sep 17 00:00:00 2001 From: Jonas Koziorek Date: Sun, 21 Jul 2024 14:41:56 +0200 Subject: [PATCH 5/7] import not needed --- src/derived_systems/tangent_space.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/derived_systems/tangent_space.jl b/src/derived_systems/tangent_space.jl index 116f514e..5edc6458 100644 --- a/src/derived_systems/tangent_space.jl +++ b/src/derived_systems/tangent_space.jl @@ -138,7 +138,6 @@ correct_matrix_type(::Val{true}, Q::AbstractMatrix) = ismutable(Q) ? Q : Array(Q ########################################################################################### # Creation of tangent rule ########################################################################################### -import ForwardDiff # IIP Tangent space dynamics function tangent_rule(f::F, J::JAC, J0, ::Val{true}, ::Val{k}, u0) where {F, JAC, k} tangentf = (du, u, p, t) -> begin From ac6f54bce3303eda3e8e0a75d349379223162eb4 Mon Sep 17 00:00:00 2001 From: Jonas Koziorek Date: Sun, 21 Jul 2024 15:13:43 +0200 Subject: [PATCH 6/7] docstring --- docs/src/index.md | 1 + src/core_systems/jacobian.jl | 20 ++++++++++++++++++-- 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 57e336df..342a9f67 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -74,6 +74,7 @@ CoreDynamicalSystem TangentDynamicalSystem current_deviations set_deviations! +jacobian orthonormal ``` diff --git a/src/core_systems/jacobian.jl b/src/core_systems/jacobian.jl index 873a7053..4face80e 100644 --- a/src/core_systems/jacobian.jl +++ b/src/core_systems/jacobian.jl @@ -2,6 +2,22 @@ export jacobian import ForwardDiff +""" + jacobian(ds::CoreDynamicalSystem) + +Construct the Jacobian rule for the dynamical system `ds`. +This is done via automatic differentiation using module +[`ForwardDiff`](https://github.com/JuliaDiff/ForwardDiff.jl). + +## Description + +For out-of-place systems, `jacobian` returns the Jacobian rule as a +function `Jf(u, p, t) -> J0::SMatrix`. Calling `Jf(u, p, t)` will compute the Jacobian +at the state `u`, parameters `p` and time `t` and return the result as `J0`. +For in-place systems, `jacobian` returns the Jacobian rule as a function +`Jf!(J0, u, p, t)`. Calling `Jf!(J0, u, p, t)` will compute the Jacobian +at the state `u`, parameters `p` and time `t` and save the result in `J0`. +""" function jacobian(ds::CoreDynamicalSystem{IIP}) where {IIP} _jacobian(ds, Val{IIP}()) end @@ -12,7 +28,7 @@ function _jacobian(ds, ::Val{true}) cfg = ForwardDiff.JacobianConfig( (du, u) -> f(du, u, p, p), deepcopy(u0), deepcopy(u0) ) - Jf = (J0, u, p, t) -> begin + Jf! = (J0, u, p, t) -> begin uv = @view u[:, 1] du = copy(u) ForwardDiff.jacobian!( @@ -20,7 +36,7 @@ function _jacobian(ds, ::Val{true}) ) nothing end - return Jf + return Jf! end function _jacobian(ds, ::Val{false}) From fbe3bd21f1c931a94e6318e7364ea44a791f6c0f Mon Sep 17 00:00:00 2001 From: Jonas Koziorek Date: Sun, 21 Jul 2024 16:09:23 +0200 Subject: [PATCH 7/7] changelog and version change --- CHANGELOG.md | 5 +++++ Project.toml | 2 +- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f7811d27..ec799380 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,8 @@ +# v3.9.0 + +A new function `jacobian` that generates Jacobian rule for any `ds<:CoreDynamicalSystem` via +automatic differentiation with `ForwardDiff.jl`. + # v3.8.0 Even more amazing level of integration with ModelingToolkit.jl! diff --git a/Project.toml b/Project.toml index 72cac47e..6695a2d7 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "DynamicalSystemsBase" uuid = "6e36e845-645a-534a-86f2-f5d4aa5a06b4" repo = "https://github.com/JuliaDynamics/DynamicalSystemsBase.jl.git" -version = "3.8.3" +version = "3.9.0" [deps] ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"