Merge pull request #86 from JuliaDynamics/internalnorm

New internalnorm based on state only

Author: Datseris

CHANGELOG.md

@@ -1,5 +1,9 @@
 Changelog of `DynamicalSystemsBase`.
 
+# v1.3
+* The specialized integrators (tangent & parallel) now implement an internal norm that only evaluates norm of the main state, instead of using the other parallel states or deviation vectors. (#86)
+* Bunch of bugfixes and performance improvements (see git history)
+
 # v1.2
 In version `1.2` we have moved to `DiffEqBase 5.0+`. In addition to that we have changed the default integrator to `SimpleATsit5` from module `SimpleDiffEq`. This is not a breaking change and you can use any of the previous integrators. **It must be stated though that numeric results you obtained using the default integrator will now be slightly different**.
 

Project.toml

@@ -1,7 +1,7 @@
 name = "DynamicalSystemsBase"
 uuid = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
 repo = "https://github.com/JuliaDynamics/DynamicalSystemsBase.jl.git"
-version = "1.2.5"
+version = "1.3.0"
 
 [deps]
 DelayEmbeddings = "5732040d-69e3-5649-938a-b6b4f237613f"
@@ -15,12 +15,12 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-DiffEqBase = "5.0.0"
+DiffEqBase = "5.10.0"
 DiffEqCallbacks = "≥ 2.1.0"
 ForwardDiff = "≥ 0.8.0"
 OrdinaryDiffEq = "≥ 5.0.0"
 Reexport = "≥ 0.1.0"
-SimpleDiffEq = "≥ 0.3.0"
+SimpleDiffEq = "≥ 0.4.0"
 StaticArrays = "≥ 0.8.0"
 julia = "≥ 1.0.0"
 

src/continuous.jl

@@ -63,9 +63,19 @@ function tangent_integrator(ds::CDS{IIP}, Q0::AbstractMatrix;
     tanprob = ODEProblem{IIP}(tangentf, hcat(u, Q), (t0, typeof(t0)(Inf)), ds.p)
 
     solver = _get_solver(diffeq)
-    return  __init(tanprob, solver; DEFAULT_DIFFEQ_KWARGS...,
-                   save_everystep = false, diffeq...)
+    return __init(tanprob, solver; DEFAULT_DIFFEQ_KWARGS..., internalnorm = _tannorm,
+                  save_everystep = false, diffeq...)
+end
+
+function _tannorm(u::AbstractMatrix, t)
+    s = size(u)[1]
+    x = zero(eltype(u))
+    for i in 1:s
+        @inbounds x += u[i, 1]^2
+    end
+    return sqrt(x/length(x))
 end
+_tannorm(u::Real, t) = abs(u)
 
 # Auto-diffed in-place version
 function tangent_integrator(ds::CDS{true, S, D, F, P, JAC, JM, true},
@@ -86,7 +96,7 @@ function tangent_integrator(ds::CDS{true, S, D, F, P, JAC, JM, true},
 
     solver = _get_solver(diffeq)
     return __init(tanprob, solver; DEFAULT_DIFFEQ_KWARGS..., save_everystep = false,
-                diffeq...)
+                  internalnorm = _tannorm, diffeq...)
 end
 
 ############################### Parallel ############################################
@@ -117,10 +127,18 @@ function parallel_integrator(ds::CDS, states; diffeq...)
     # if typeof(solver) ∈ STIFFSOLVERS
     #     error("Stiff solvers can't support a parallel integrator.")
     # end
-    return __init(pprob, solver; DEFAULT_DIFFEQ_KWARGS..., save_everystep = false,
-                  diffeq...)
+    if !(typeof(ds) <: CDS{true})
+        return __init(pprob, solver; DEFAULT_DIFFEQ_KWARGS..., save_everystep = false,
+                      internalnorm = _parallelnorm, diffeq...)
+    else
+        return __init(pprob, solver; DEFAULT_DIFFEQ_KWARGS..., save_everystep = false,
+                      internalnorm = _tannorm, diffeq...)
+    end
 end
 
+_parallelnorm(u::AbstractVector, t) = @inbounds DiffEqBase.ODE_DEFAULT_NORM(u[1], t)
+_parallelnorm(u::Real, t) = abs(u)
+
 #####################################################################################
 #                                 Trajectory                                        #
 #####################################################################################

test/dynsys_types.jl

@@ -66,7 +66,7 @@ for i in 1:8
 
         # Currently the in-place version does not work from DiffEq's side:
         if i > 2
-            pinteg = parallel_integrator(ds, [INITCOD[sysindx], INITCOD[sysindx]]; alg = alg)
+            pinteg = parallel_integrator(ds, [copy(INITCOD[sysindx]), copy(INITCOD[sysindx])]; alg = alg)
             puprev = deepcopy(get_state(pinteg))
             step!(pinteg)
             @test get_state(pinteg, 1) == get_state(pinteg, 2) == get_state(pinteg)

test/norm_tests.jl

@@ -0,0 +1,44 @@
+using DynamicalSystemsBase, SimpleDiffEq, DiffEqBase
+using Statistics, Test
+
+ds = Systems.lorenz()
+ALG = SimpleATsit5()
+
+# %%
+
+i1 = tangent_integrator(ds, 3; alg = ALG)
+i2 = tangent_integrator(ds, 3; alg = ALG,
+internalnorm = DiffEqBase.ODE_DEFAULT_NORM)
+
+step!(i1)
+step!(i2)
+dts = zeros(2, 500)
+for i in 1:500
+    step!(i1); step!(i2)
+    dts[:, i] .= (i1.dt, i2.dt)
+end
+
+dt1 = mean(dts[1, :])
+dt2 = mean(dts[2, :])
+
+@test dt1 > dt2
+
+# %%
+s = [get_state(ds), get_state(ds) .+ rand(3)]
+i1 = parallel_integrator(ds, deepcopy(s);
+alg = ALG, internalnorm = DiffEqBase.ODE_DEFAULT_NORM)
+i2 = parallel_integrator(ds, deepcopy(s);
+alg = ALG)
+
+step!(i1)
+step!(i2)
+dts = zeros(2, 1000)
+for i in 1:1000
+    step!(i1); step!(i2)
+    dts[:, i] .= (i1.dt, i2.dt)
+end
+
+dt1 = mean(dts[1, :])
+dt2 = mean(dts[2, :])
+
+@test dt2 > dt1

test/runtests.jl

@@ -9,6 +9,7 @@ include("dynsys_inference.jl")
 include("continuous_systems.jl")
 include("discrete_systems.jl")
 include("integrators_with_callbacks.jl")
+include("norm_tests.jl")
 
 ti = time() - ti
 println("\nTest took total time of:")