Allow initial conditions to depend on parameters during a global continuation (#183)

* allow parameter-dependent continuation

* versions

* add test for the new feature

Author: Datseris

CHANGELOG.md

@@ -1,3 +1,7 @@
+# v1.27
+
+- It is now possible to specify initial conditions that depend on the parameters during a global continuation. See `PerParameterInitialConditions`.
+
 # v1.26
 
 - **Brand new, ground-breaking functionality** for finding (almost) all stability

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.26.4"
+version = "1.27.0"
 
 
 [deps]

src/continuation/continuation_ascm_generic.jl

@@ -128,7 +128,10 @@ function global_continuation(acam::AttractorSeedContinueMatch, pcurve, ics;
     # first parameter is run in isolation, as it has no prior to seed from
     set_parameters!(referenced_dynamical_system(mapper), pcurve[1])
     if ics isa Function
-        fs = basins_fractions(mapper, ics; show_progress = false, N = samples_per_parameter)
+        fs = basins_fractions(mapper, ics; show_progress = false, N)
+    elseif ics isa PerParameterInitialConditions
+        u0s = ics.generator(pcurve[1], N)
+        fs, = basins_fractions(mapper, u0s; show_progress = false)
     else # we ignore labels in this continuation algorithm
         fs, = basins_fractions(mapper, ics; show_progress = false)
     end
@@ -152,7 +155,7 @@ function global_continuation(acam::AttractorSeedContinueMatch, pcurve, ics;
         fs = if allows_mapper_u0(mapper)
             seed_attractors_to_fractions_individual(mapper, prev_attractors, ics, N, acam.seeding)
         else
-            seed_attractors_to_fractions_grouped(mapper, prev_attractors, ics, N, acam.seeding)
+            seed_attractors_to_fractions_grouped(mapper, prev_attractors, ics, N, acam.seeding, p)
         end
         current_attractors = deepcopy(extract_attractors(mapper))
         # we don't match attractors here, this happens directly at the end.
@@ -196,7 +199,7 @@ function seed_attractors_to_fractions_individual(mapper, prev_attractors, ics, N
     return fs
 end
 
-function seed_attractors_to_fractions_grouped(mapper, prev_attractors, ics, N, seeding)
+function seed_attractors_to_fractions_grouped(mapper, prev_attractors, ics, N, seeding, parameters)
     # what makes this version different is that we can't just use `mapper(u0)`,
     # so we need to store the seeded initial conditions and then combine them with
     # the the ones generated from `ics`.
@@ -212,6 +215,8 @@ function seed_attractors_to_fractions_grouped(mapper, prev_attractors, ics, N, s
         for _ in 1:N
             push!(u0s, copy(ics()))
         end
+    elseif ics isa PerParameterInitialConditions
+        append!(u0s, ics.generator(parameters, N))
     else
         append!(u0s, vec(ics))
     end

src/continuation/continuation_grouping.jl

@@ -87,7 +87,12 @@ function _get_features_pcurve(mapper::AttractorsViaFeaturizing, ics, n, spp, pcu
     # Collect features
     for (i, p) in enumerate(pcurve)
         set_parameters!(mapper.ds, p)
-        current_features = extract_features(mapper, ics; show_progress, N = spp)
+        if ics isa PerParameterInitialConditions
+            u0s = ics.generator(p, spp)
+        else
+            u0s = ics
+        end
+        current_features = extract_features(mapper, u0s; show_progress, N = spp)
         features[((i - 1)*spp + 1):i*spp] .= current_features
         ProgressMeter.next!(progress)
     end

src/continuation/global_continuation_api.jl

@@ -1,4 +1,5 @@
 export global_continuation, GlobalContinuationAlgorithm, continuation_series, stability_measures_along_continuation
+export PerParameterInitialConditions
 
 """
     GlobalContinuationAlgorithm
@@ -17,9 +18,22 @@ abstract type GlobalContinuationAlgorithm end
     global_continuation(gca::GlobalContinuationAlgorithm, pcurve, ics; kwargs...)
 
 Find and continue attractors (or representations of attractors)
-and the fractions of their basins of attraction across a parameter range `pcurve`
+and the fractions of their basins of attraction across a parameter curve `pcurve`
 by sampling given initial conditions `ics` according to algorithm `gca`.
 
+Possible subtypes of a `GlobalContinuationAlgorithm` are:
+
+- [`AttractorSeedContinueMatch`](@ref)
+- [`FeaturizeGroupAcrossParameter`](@ref)
+
+`ics` are the initial conditions to use when sampling the state space.
+They can be specified in one of three ways:
+
+1. A set vector of initial conditions (vector of vectors).
+2. A 0-argument function that generates random initial conditions.
+3. The special type [`PerParameterInitialConditions`](@ref) that allows
+   different initial conditions for different parameter values.
+
 Return:
 
 1. `fractions_cont::Vector{Dict{Int, Float64}}`. The fractions of basins of attraction.
@@ -57,27 +71,31 @@ across the parameter range `prange`, for the parameter of the system with index
 traditional continuation (see online Tutorial for a comparison), global continuation
 can be performed over arbitrary user-defined curves in parameter space.
 The second call signature with `pcurve` allows for this possibility. In this case
-`pcurve` is a vector of iterables, where each itereable maps parameter indices
+`pcurve` is a vector of iterables, where each iterable maps parameter indices
 to parameter values. These iterables can be dictionaries, named tuples, `Vector{Pair}`,
-etc., and the sequence of the iterables defines a curve in parameter space.
+anything that can be given in `set_parameters!`.
+The sequence of the iterables defines a curve in parameter space.
 In fact, the version with `prange, pidx` simply defines
 `pcurve = [[pidx => p] for p in prange]` and calls the second method.
-
-`ics` are the initial conditions to use when globally sampling the state space.
-Like in [`basins_fractions`](@ref) it can be either a set vector of initial conditions,
-or a 0-argument function that generates random initial conditions.
-
-Possible subtypes of `GlobalContinuationAlgorithm` are:
-
-- [`AttractorSeedContinueMatch`](@ref)
-- [`FeaturizeGroupAcrossParameter`](@ref)
 """
 function global_continuation(alg::GlobalContinuationAlgorithm, prange::AbstractVector, pidx, sampler; kw...)
     # everything is propagated to the curve setting
     pcurve = [[pidx => p] for p in prange]
     return global_continuation(alg, pcurve, sampler; kw...)
 end
 
+"""
+    PerParameterInitialConditions(generator)
+
+Wrapper around a function `generator`, to be called as
+`generator(parameter_pairs, N::Int)`.
+It inputs the current parameter(s) of a [`global_continuation`](@ref)
+(elements of `pcurve`), and generates a vector of `N` initial conditions.
+"""
+struct PerParameterInitialConditions{F}
+    generator::F
+end
+
 
 """
     continuation_series(continuation_info, fillval = NaN)

test/continuation/seed_continue_generic.jl

@@ -85,4 +85,53 @@ using Random
         end
     end
 
+end
+
+@testset "parameter dependent ic" begin
+    # fake map with two equlibria at (-r, r) and (r, r)
+    function dumb_map2(dz, z, p, n)
+        x, y = z
+        r = p[1]
+
+        if y < r - 1
+            dz[1] = dz[2] = 0.0
+        else
+            if x < 0.0
+                dz[1] = -1.0
+                dz[2] = r
+            else
+                dz[1] = 1.0
+                dz[2] = r
+            end
+        end
+        return nothing
+    end
+
+    rs = [0, 2.0, 4.0]
+    ds = DeterministicIteratedMap(dumb_map2, [0., 0.], [1.0])
+    N = 5
+
+    function make_ics(parameters, n)
+        I = parameters[1].second
+        @show I
+        N = round(Int, sqrt(n))
+        grid = (range(-1.0, 1.0; length = N), range(-0.5, 0.5; length = N))
+        ics = [SVector(v, w + I) for w in grid[2] for v in grid[1]]
+        return ics
+    end
+    icsgen = PerParameterInitialConditions(make_ics)
+
+    featurizer(A, t) = A[end]
+    gconfig = GroupViaPairwiseComparison(threshold = 0.25, rescale_features = false)
+    mapper = AttractorsViaFeaturizing(ds, featurizer, gconfig; Ttr = 1, T = 1)
+    gca = AttractorSeedContinueMatch(mapper; seeding = A -> [])
+
+    pcurve = [[1 => r] for r in rs]
+    fractions_cont, attractors_cont = global_continuation(gca, pcurve, icsgen; samples_per_parameter = 25)
+
+    # all times 2 attractors, never the 0.0 attractor
+    @test all(A -> length(A) == 2, attractors_cont)
+    for attractors in attractors_cont
+        @test all(A -> A[end][1] != 0, values(attractors))
+    end
 end