feat: use attractors.jl solvers

Project.toml

@@ -4,6 +4,7 @@ authors = ["Orjan Ameye <orjan.ameye@hotmail.com>", "Jan Kosata <kosataj@phys.et
 version = "0.2.6"
 
 [deps]
+Attractors = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7"
 BijectiveHilbert = "91e7fc40-53cd-4118-bd19-d7fcd1de2a54"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"

attractors.jl

@@ -0,0 +1,93 @@
+using HarmonicBalance, HarmonicSteadyState
+using Symbolics, RuntimeGeneratedFunctions, QuestBase, FunctionWrappers, SciMLBase
+using ModelingToolkit
+import FunctionWrappers: FunctionWrapper
+using Attractors
+
+@variables γ λ α ω0 ω t
+@variables x(t)
+
+nat_eq = d(d(x, t), t) + γ * d(x, t) + ω0^2 * (1 - λ * cos(2 * ω * t)) * x + α * x^3
+dEOM = DifferentialEquation(nat_eq, x)
+add_harmonic!(dEOM, x, ω)
+harmonic_eq = get_harmonic_equations(dEOM);
+
+struct AttractorsProblem{
+    ParType<:Number,
+    Jac<:JacobianFunction(ComplexF64),
+} <: SteadyStateProblem
+    variables::Vector{Num}
+    parameters::Vector{Num}
+    swept_parameters::OrderedDict{Num,Vector{ParType}}
+    fixed_parameters::OrderedDict{Num,ParType}
+    phase_space::NTuple{N,Vector{ParType}}
+    coupled_ode::DynamicalSystemBase.CoupledODEs
+    jacobian::Jac
+    eom::HarmonicEquation
+end
+
+ODEProblem(
+    harmonic_eq,
+    zeros(2),
+    (0, Inf),
+    Dict(ω0 => 1.0, γ => 0.1, λ => 0.05, α => 0.01),
+)
+
+u_range = range(-0.2, 0.2; length=75)
+u_range |> typeof
+grid = (
+    u_range, # u1
+    u_range, # v1
+    u_range, # u2
+    u_range, # v2
+    u_range, # u3
+    u_range, # v3
+)
+mapper = AttractorsViaRecurrences(
+    ds,
+    grid;
+    consecutive_recurrences=5000,
+    attractor_locate_steps=5000,
+    consecutive_lost_steps=500,
+)
+
+fs = basins_fractions(mapper, sampler)
+attractors = extract_attractors(mapper)
+
+plot_attractors(attractors; access=SVector(1, 2))
+
+continuation_range = range(1e-4, 0.002, 100)
+
+sampler, = statespace_sampler(grid)
+algo = AttractorSeedContinueMatch(mapper)
+
+reinit!(ds)
+fractions_cont, attractors_cont = global_continuation(
+    algo, continuation_range, :F₂, sampler; samples_per_parameter=50
+)
+
+
+# eqs = rearrange_standard(harmonic_eq).equations
+# rhs = Num[eq.lhs for eq in eqs]
+# vars = get_variables(harmonic_eq)
+
+# function compile_function(
+#     rhs::Vector{Num}, variables::Vector{Num}; rules=Dict()
+# )::RuntimeGeneratedFunction
+#     rhss = QuestBase.substitute_all.(rhs, Ref(rules)) # Ref makes sure only mat is broadcasted
+#     rhsf = Symbolics.build_function(rhss, variables; expression=Val(false))
+#     return rhsf isa Tuple ? first(rhsf) : rhsf
+# end
+
+# rgf = compile_function(
+#     rhs, vcat(vars, ω); rules=Dict(ω0 => 1.0, γ => 0.1, λ => 0.05, α => 0.01)
+# ) # RuntimeGeneratedFunction
+
+# function f(u,p,t)
+#     rgf(vcat(u,p))
+# end
+# ff = FunctionWrapper{Vector{Float64},Tuple{Vector{Float64},Vector{Float64},Float64}}(f)
+# ff(rand(2), rand(1), 0.0)
+
+# CoupledODEs(ff, rand(2) , rand(1))
+# SciMLBase.isinplace(ff, 100)

docs/Project.toml

@@ -1,17 +1,26 @@
 [deps]
+Attractors = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
+DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
+FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
+HarmonicBalance = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e"
 HarmonicSteadyState = "1158f75c-a779-4b85-8bfb-8fcf6bf02ced"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
-HarmonicBalance = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+QuestBase = "7e80f742-43d6-403d-a9ea-981410111d43"
+RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
-
-[sources]
-HarmonicBalance = {url = "https://github.com/QuantumEngineeredSystems/HarmonicBalance.jl", rev = "QuestBase"}
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 Documenter = "1"
-julia = "1.11"
+julia = "1.10"
+
+[sources.HarmonicBalance]
+rev = "QuestBase"
+url = "https://github.com/QuantumEngineeredSystems/HarmonicBalance.jl"