Fix of match sequentially with past (#170)

* imrpove comparison with BFKit

* allow default option in next id

* Revert "allow default option in next id"

This reverts commit 3becca5.

* critical fix of rematch with past with empty attractors

when eg. all initial conditions diverge

* bump project

* add a todo, I am too tired to do this now...

* add test for interrupted attractors with vanished

* make sure empty dictionaries are handled in the matching maps

* make empty attractors handled ONLY at matching maps

* bump changelog

Author: Datseris

Diff for: CHANGELOG.md

@@ -1,6 +1,7 @@
 # v1.24
 
 - Improved `continuation_series`: now it works for any type of input, including attractor continuation info because the default value (now renamed to fill value) is used to deduce the type of the series vectors.
+- Fixed a critical bug where global continuation would fail if there was a parameter where all initial conditions diverge and the matching was with vanished.
 
 # v1.23
 

Diff for: 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.24.2"
+version = "1.24.3"
 
 
 [deps]

Diff for: src/dict_utils.jl

@@ -104,14 +104,21 @@ end
 
 Return the minimum key of the "new" dictionary
 that doesn't exist in the "old" dictionary.
+If one of the two dictionaries are empty, return its maximum key + 1.
+If both are empty, return 1.
+
+The function assumes tha the dictionary keys are integers.
 """
-function next_free_id(a₊::AbstractDict, a₋::AbstractDict)
-    s = setdiff(keys(a₊), keys(a₋))
-    nextid = isempty(s) ? maximum(keys(a₋)) + 1 : minimum(s)
-    return nextid
-end
+next_free_id(a₊::AbstractDict, a₋::AbstractDict) = next_free_id(keys(a₊), keys(a₋))
 function next_free_id(keys₊, keys₋)
-    s = setdiff(keys₊, keys₋)
-    nextid = isempty(s) ? maximum(keys₋) + 1 : minimum(s)
-    return nextid
+    if length(keys₋) == 0 && length(keys(keys₊)) == 0
+        return 1
+    elseif length(keys₋) == 0
+        return maximum(keys₊) + 1
+    elseif length(keys₊) == 0
+        return maximum(keys₋) + 1
+    else
+        s = setdiff(keys₊, keys₋)
+        return isempty(s) ? maximum(keys₋) + 1 : minimum(s)
+    end
 end

Diff for: src/matching/basin_enclosure.jl

@@ -59,6 +59,11 @@ function matching_map(
         current_attractors, prev_attractors, matcher::MatchByBasinEnclosure;
         ds, p, pprev = nothing, next_id = next_free_id(current_attractors, prev_attractors)
     )
+    # if either dictionary is empty there is no matching to be done
+    if isempty(current_attractors) || isempty(prev_attractors)
+        return Dict{keytype(a₊), keytype(a₋)}()
+    end
+
     if matcher.ε === nothing
         e = ε_from_centroids(current_attractors)
     else

Diff for: src/matching/matching_interface.jl

@@ -44,6 +44,8 @@ The values contained in `a₊, a₋` can be anything supported by `matcher`.
 Within Attractors.jl they are typically `StateSpaceSet`s representing attractors.
 Typically the +,- mean after and before some change of parameter of a dynamical system.
 
+`matching_map` always returns an empty dictionary of either `a₊, a₋` is empty.
+
 ## Keyword arguments
 
 - `ds`: the dynamical system that generated `a₊, a₋`.
@@ -168,15 +170,14 @@ function _rematch_ignored!(attractors_cont, matcher;
     for i in 1:length(attractors_cont)-1
         a₊, a₋ = attractors_cont[i+1], attractors_cont[i]
         p, pprev = pcurve[i+1], pcurve[i]
-        # If there are no attractors, skip the matching
-        if (isempty(a₊) || isempty(a₋))
-            push!(rmaps, Dict{keytype(attractors_cont[1]), keytype(attractors_cont[1])}())
-            continue
+        # If there attractors, update the max id
+        if !(isempty(a₊) || isempty(a₋))
+            # we always compute a next id. In this way, if an attractor disappears
+            # and reappears, it will get a different (incremented) ID as it should!
+            next_id_a = max(maximum(keys(a₊)), maximum(keys(a₋)))
+            next_id = max(next_id, next_id_a) + 1
         end
-        # Here we always compute a next id. In this way, if an attractor disappears
-        # and reappears, it will get a different (incremented) ID as it should!
-        next_id_a = max(maximum(keys(a₊)), maximum(keys(a₋)))
-        next_id = max(next_id, next_id_a) + 1
+        # matching_map returns empty dict if the inputs are empty dicts
         rmap = matching_map!(a₊, a₋, matcher; next_id, ds, p, pprev)
         push!(rmaps, rmap)
     end
@@ -194,11 +195,12 @@ function _rematch_with_past!(attractors_cont, matcher;
     rmaps = Dict{keytype(attractors_cont[1]), keytype(attractors_cont[1])}[]
     for i in 1:length(attractors_cont)-1
         a₊, a₋ = attractors_cont[i+1], attractors_cont[i]
-        p, pprev = pcurve[i+1], pcurve[i]
-        # update ghosts
+        # first update ghosts
         for (k, A) in a₋
             latest_ghosts[k] = A
         end
+        # and then match
+        p, pprev = pcurve[i+1], pcurve[i]
         rmap = matching_map!(a₊, latest_ghosts, matcher; pprev, p, ds)
         push!(rmaps, rmap)
     end

Diff for: src/matching/sssdistance.jl

@@ -54,6 +54,8 @@ _use_vanished(m::MatchBySSSetDistance) = m.use_vanished
 function matching_map(a₊::AbstractDict, a₋::AbstractDict, matcher::MatchBySSSetDistance;
         kw...
     )
+    # if either dictionary is empty there is no matching to be done
+    isempty(a₊) || isempty(a₋) && return Dict{keytype(a₊), keytype(a₋)}()
     distances = setsofsets_distances(a₊, a₋, matcher.distance)
     keys₊, keys₋ = sort.(collect.(keys.((a₊, a₋))))
     _matching_map_distances(keys₊, keys₋, distances::Dict, matcher.threshold; kw...)

Diff for: test/continuation/matching_attractors.jl

@@ -141,6 +141,23 @@ end
                 end
             end
         end
+        @testset "no attractors for a step" begin
+            atts = deepcopy(allatts)
+            atts[5] = Dict()
+            match_sequentially!(atts, MatchBySSSetDistance(use_vanished = true))
+            @test unique_keys(atts) == 1:3
+            for i in eachindex(jrange)
+                if iseven(i)
+                    @test sort!(collect(keys(atts[i]))) == 1:2
+                else
+                    if i == 5
+                        @test isempty(atts[i])
+                    else
+                        @test sort!(collect(keys(atts[i]))) == 1:3
+                    end
+                end
+            end
+        end
     end
 end