use current parameters in jacobian... (#191)

* use current parameters in jacobian...;

* Add test for parameter usage in jacobian during continuation

* Add test for parameter usage in jacobian during continuation

* bump version

---------

Co-authored-by: andreasmorr <[email protected]>

Author: Datseris

Project.toml

@@ -2,7 +2,7 @@ name = "Attractors"
 uuid = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7"
 authors = ["George Datseris <[email protected]>", "Kalel Rossi", "Alexandre Wagemakers"]
 repo = "https://github.com/JuliaDynamics/Attractors.jl.git"
-version = "1.31"
+version = "1.31.1"
 
 [deps]
 BlackBoxOptim = "a134a8b2-14d6-55f6-9291-3336d3ab0209"

src/mapping/stability_measures_accumulator.jl

@@ -333,10 +333,10 @@ function finalize_accumulator(accumulator::StabilityMeasuresAccumulator)
             if isinplace(ds)
               # For in-place systems, pre-allocate J and then compute it
               J = Array{Float64}(undef, length(A[1]), length(A[1]))
-              jac(J, Array(A[1]), initial_parameters(ds), 0)
+              jac(J, Array(A[1]), current_parameters(ds), 0)
             else
               # For out-of-place systems, compute J directly
-              J = jac(Array(A[1]), initial_parameters(ds), 0)
+              J = jac(Array(A[1]), current_parameters(ds), 0)
             end
 
             λ = min(0, maximum(real.(eigvals(J))))

test/mapping/stability_measures_accumulator.jl

@@ -131,12 +131,16 @@ end
 # Now we will test the local stability measures in `StabilityMeasuresAccumulator` in a
 # linear system.
 function linear_evolution(z, p, n)
-    A = [-0.5 0.0; 0.0 -0.5]  # Linear transformation matrix
+    if p[1] < 0.0
+        A = [1.0 0.0; 0.0 1.0]  # Linear transformation matrix
+    else
+        A = [-0.5 0.0; 0.0 -0.5]  # Linear transformation matrix
+    end
     return SVector(A * [z[1], z[2]]...)  # Convert matrix multiplication result to SVector
 end
 
 # Create the dynamical system
-dynamics = CoupledODEs(linear_evolution, [1.0, 1.0], [0.0])
+dynamics = CoupledODEs(linear_evolution, [1.0, 1.0], [1.0])
 
 # Define a grid for initial conditions
 grid = ([-1.0, -0.1, 0.3, 1.0], [-1.0, -0.3, 0.1, 1.0])
@@ -171,6 +175,41 @@ results_expected = Dict(
     end
 end
 
+# Now we test the continuation of local stability measures using the linear map.
+pcurve_local = [[1 => p] for p in [-1.0, 1.0]]
+attractors_cont_local = [
+    Dict(1 => StateSpaceSet([SVector(0.0, 0.0)])),
+    Dict(1 => StateSpaceSet([SVector(0.0, 0.0)]))
+]
+
+proximity_mapper_options_local = (
+    Ttr=0, stop_at_Δt = false, horizon_limit = 1e2, consecutive_lost_steps = 10000
+)
+measures_cont_local = stability_measures_along_continuation(
+    dynamics, attractors_cont_local, pcurve_local, ics_from_grid(grid), ε=0.1, finite_time=0.5,
+    proximity_mapper_options = proximity_mapper_options_local
+)
+
+measures_cont_local_expected = Dict(
+    "characteristic_return_time" => [Dict(1 => Inf, -1 => NaN), Dict(1 => 2.0, -1 => NaN)],
+    "reactivity"       => [Dict(1 => 1.0, -1 => NaN), Dict(1 => -0.5, -1 => NaN)],
+    "maximal_amplification" => [Dict(1 => Inf, -1 => NaN), Dict(1 => 1.0, -1 => NaN)],
+    "maximal_amplification_time" => [Dict(1 => Inf, -1 => NaN), Dict(1 => 0.0, -1 => NaN)]
+)
+@testset "Local Stability Measures Continuation" begin
+    # Validate the results
+    for (key, value) in measures_cont_local_expected
+        @test key in keys(measures_cont_local)
+        for k in [1, 2]
+            @test isapprox(
+                sort(collect(values(value[k]))),
+                sort(collect(values(measures_cont_local[key][k]))),
+                atol = 1e-5,
+                nans = true,
+            )
+        end
+    end
+end
 
 @testset "Discrete time" begin
     # For these parameters the map has 1 fixed point and one period 3 orbit.