Limit cycle classes (#55)

* fixed loading for version changes

* ordering revamped to have fixed order

* limit cycle unique classification

* 2d result test added

Project.toml

@@ -1,7 +1,7 @@
 name = "HarmonicBalance"
 uuid = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e"
 authors = ["Jan Kosata <kosataj@phys.ethz.ch>", "Javier del Pino <jdelpino@phys.ethz.ch>"]
-version = "0.6.1"
+version = "0.6.2"
 
 [deps]
 BijectiveHilbert = "91e7fc40-53cd-4118-bd19-d7fcd1de2a54"

src/classification.jl

@@ -21,14 +21,19 @@ classify_solutions!(res, "sqrt(u1^2 + v1^2) > 1.0" , "large_amplitude")
 
 """
 function classify_solutions!(res::Result, condition::String, name::String; physical=true)
+    values = classify_solutions(res, condition; physical=physical)
+    res.classes[name] = values
+end
+
+
+function classify_solutions(res::Result, condition::String; physical=true)
     expr = Num(eval(Meta.parse(condition)))
     function cond_func(s::OrderedDict, res)
         physical && !is_physical(s, res) && return false
         s = [key => real(s[key]) for key in keys(s)] # make values real
         Bool(substitute_all(expr, s).val)
     end
-    values = classify_solutions(res, cond_func)
-    res.classes[name] = values
+    classify_solutions(res, cond_func)
 end
 
 
@@ -43,7 +48,7 @@ end
 "Return an array of booleans classifying the solution in `res`
 according to `f` (`f` takes a solution dictionary, return a boolean)"
 function classify_solutions(res::Result, f::Function)
-    values = similar(res.solutions, Vector{Bool})
+    values = similar(res.solutions, BitVector)
     for (idx, soln) in enumerate(res.solutions)
         values[idx] = [f(get_single_solution(res, index=idx, branch=b), res) for b in 1:length(soln)]
     end

src/modules/LimitCycles.jl

@@ -5,5 +5,6 @@ using Symbolics
 using DocStringExtensions
 
 include("LimitCycles/Hopf.jl")
+include("LimitCycles/analysis.jl")
 
 end

src/modules/LimitCycles/Hopf.jl

@@ -4,7 +4,7 @@ export replace_Hopf_variable!
 
 using HarmonicBalance: is_rearranged, rearrange_standard, _remove_brackets
 using HarmonicBalance.LinearResponse: get_implicit_Jacobian, get_Jacobian
-import HarmonicBalance: is_stable, is_physical, is_Hopf_unstable, order_branches!, classify_binaries!
+import HarmonicBalance: is_stable, is_physical, is_Hopf_unstable, order_branches!, classify_binaries!, find_branch_order
 #import HarmonicBalance.get_steady_states; export get_steady_states
 
 
@@ -97,7 +97,13 @@ function _classify_limit_cycles!(res::Result, Δω::Num)
         res.classes[c][idx] .*= abs.(getindex.(res.solutions[idx], Δω_idx)) .> 1E-10
     end
 
-    order_branches!(res, ["physical", "stable", "Hopf"]) # shuffle the branches to have relevant ones first
+    classify_unique!(res, Δω)
+
+    unique_stable = find_branch_order(map(.*, res.classes["stable"], res.classes["unique"]))
+
+    # branches which are unique but never stable
+    unique_unstable = setdiff(find_branch_order(map(.*, res.classes["unique"], res.classes["physical"])), unique_stable)
+    order_branches!(res, vcat(unique_stable, unique_unstable)) # shuffle to have relevant branches first
 end
 
 

src/modules/LimitCycles/analysis.jl

@@ -0,0 +1,15 @@
+
+
+import HarmonicBalance.classify_solutions
+
+function classify_unique!(res::Result, Δω; class_name="unique")
+
+    # 1st degeneracy: arbitrary sign of Δω
+    c1 = classify_solutions(res, string(Δω) * ">= 0")
+
+    # 2nd degeneracy: ambiguity in the fixed phase, manifests as the sign of var
+    var = HarmonicBalance._remove_brackets(get_Hopf_variables(res.problem.eom, Δω)[1])
+    c2 = classify_solutions(res, string(var) * ">= 0") 
+
+    res.classes[class_name] = map(.*, c1, c2)
+end

src/saving.jl

@@ -19,8 +19,7 @@ function save(filename, x::Result)
     x_nofunc = deepcopy(x)
 
     # compiled functions cause problems in saving: ignore J now, compile when loading
-    null_J() = 0
-    x_nofunc.jacobian = null_J
+    x_nofunc.jacobian = 0
     JLD2.save(_jld2_name(filename), Dict("object" => x_nofunc))
 end
 
@@ -37,7 +36,6 @@ reinstated in the `HarmonicBalance` namespace.
 """
 function load(filename)
     loaded = JLD2.load(filename)
-
     if haskey(loaded,"object") #otherwise save data is from a plot
         loaded = loaded["object"]
 

src/solve_homotopy.jl

@@ -120,7 +120,7 @@ function _classify_default!(result)
     classify_solutions!(result, is_physical, "physical")
     classify_solutions!(result, is_stable, "stable")
     classify_solutions!(result, is_Hopf_unstable, "Hopf")
-    order_branches!(result, ["physical", "stable", "Hopf"]) # shuffle the branches to have relevant ones first
+    order_branches!(result, ["physical", "stable"]) # shuffle the branches to have relevant ones first
     classify_binaries!(result) # assign binaries to solutions depending on which branches are stable
 end
 
@@ -160,37 +160,38 @@ function compile_matrix(matrix, variables, fixed_parameters)
 end
 
 
-"Order the solution branches in `res` such that close classified positively by `class` are first."
-function order_branches!(res::Result, class::String)
-    indices = findall( x-> x==1, any.(classify_branch(res, class)))
-    order = cat(indices, setdiff(1:length(res.solutions[1]), indices), dims=1) # permutation of indices which puts stable branches first
-    reorder_solutions!(res, order)
+"Find a branch order according `classification`. Place branches where true occurs earlier first."
+function find_branch_order(classification::Vector{BitVector})
+    branches = [getindex.(classification, k) for k in 1:length(classification[1])] # array of branches
+    indices = replace(findfirst.(branches), nothing => Inf)
+    negative = findall(x -> x == Inf, indices) # branches not true anywhere - leave out
+    order = setdiff(sortperm(indices), negative)
 end
 
-"Order the solution branches in `res` such that close classified positively by `classes` are first.
-The order of classes has descending precedence."
-function order_branches!(s::Result, classes::Vector{String})
-    for class in reverse(classes)
-        order_branches!(s, class)
-    end
+find_branch_order(classification::Array) = collect(1:length(classification[1])) # no ordering for >1D
+
+"Order the solution branches in `res` such that close classified positively by `classes` are first."
+function order_branches!(res::Result, classes::Vector{String})
+    for class in classes
+        order_branches!(res, find_branch_order(res.classes[class]))
+    end 
 end
 
+order_branches!(res::Result, class::String) = order_branches!(res, [class])
+
 "Reorder the solutions in `res` to match the index permutation `order`."
-function reorder_solutions!(res::Result, order::Vector{Int64})
-    res.solutions = reorder_array(res.solutions, order)
+function order_branches!(res::Result, order::Vector{Int64})
+    res.solutions = _reorder_nested(res.solutions, order)
     for key in keys(res.classes)
-        res.classes[key] = reorder_array(res.classes[key], order)
+        res.classes[key] = _reorder_nested(res.classes[key], order)
     end
 end
 
-"Reorder EACH ELEMENT of `a` to match the index permutation `order`."
-function reorder_array(a::Array, order::Vector{Int64})
-    a[1] isa Array || return a
-    new_array = similar(a)
-    for (i,el) in enumerate(a)
-        new_array[i] = el[order]
-    end
-    return new_array
+"Reorder EACH ELEMENT of `a` to match the index permutation `order`. If length(order) < length(array), the remanining positions are kept."
+function _reorder_nested(a::Array, order::Vector{Int64})
+    a[1] isa Union{Array, BitVector} || return a
+    order = length(order) == length(a) ? order : vcat(order, setdiff(1:length(a[1]), order)) # pad if needed
+    new_array = [el[order] for el in a]
 end
 
 

src/types.jl

@@ -223,7 +223,7 @@ mutable struct Result
     If problem.jacobian is a symbolic matrix, this holds a compiled function.
     If problem.jacobian was `false`, this holds a function that rearranges the equations to find J
     only after numerical values are inserted (preferable in cases where the symbolic J would be very large)."
-    jacobian::Function
+    jacobian::Union{Function, Int64}
 
     Result(sol,swept, fixed, problem, classes, J) = new(sol, swept, fixed, problem, classes, J)
     Result(sol,swept, fixed, problem, classes) = new(sol, swept, fixed, problem, classes)

test/parametron.jl

@@ -7,19 +7,24 @@ using Test
 
 natural_equation =  d(d(x,t),t) + γ*d(x,t) + Ω^2*(1-λ*cos(2*ω*t+ψ))*x + α * x^3 +η *d(x,t) * x^2
 forces =  F*cos(ω*t+θ)
-dEOM = HarmonicBalance.DifferentialEquation(natural_equation + forces, x)
-HarmonicBalance.add_harmonic!(dEOM, x, ω)
-harmonic_eq = HarmonicBalance.get_harmonic_equations(dEOM, slow_time=T, fast_time=t);
+dEOM = DifferentialEquation(natural_equation + forces, x)
+add_harmonic!(dEOM, x, ω)
+harmonic_eq = get_harmonic_equations(dEOM, slow_time=T, fast_time=t);
 p = HarmonicBalance.Problem(harmonic_eq);
 
 
-fixed_parameters = (Ω => 1.0,γ => 1E-2, λ => 5E-2, F => 1E-3,  α => 1.,  η=>0.3, θ => 0, ψ => 0)
-sweep = ω => LinRange(0.9, 1.1, 100)
-soln = HarmonicBalance.get_steady_states(p, sweep, fixed_parameters)
+fixed = (Ω => 1.0,γ => 1E-2, λ => 5E-2, F => 1E-3,  α => 1.,  η=>0.3, θ => 0, ψ => 0)
+varied = ω => LinRange(0.9, 1.1, 100)
+res = get_steady_states(p, varied, fixed)
 
 # save the result, try and load in the next step
 #current_path = @__DIR__
-#HarmonicBalance.save(current_path * "/parametron_result.jld2", soln)
+#HarmonicBalance.save(current_path * "/parametron_result.jld2", res)
+
+# try to run a 2D calculation
+fixed = (Ω => 1.0,γ => 1E-2, F => 1E-3,  α => 1.,  η=>0.3, θ => 0, ψ => 0)
+varied = (ω => LinRange(0.9, 1.1, 10), λ => LinRange(0.01, 0.05, 10))
+res = get_steady_states(p, varied, fixed)
 
 
 ###