Format, restructure and invalidiations

.JuliaFormatter.toml

@@ -0,0 +1,6 @@
+style = "blue"
+align_assignment = true
+align_struct_field = true
+align_conditional = true
+align_pair_arrow = true
+align_matrix = true

.github/workflows/Format.yml

@@ -0,0 +1,9 @@
+name: Format suggestions
+
+on: [pull_request]
+
+jobs:
+  code-style:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: julia-actions/julia-format@v3

.gitignore

@@ -8,5 +8,4 @@ local_tests
 .DS_Store
 .DS_Store
 docs/build/**
-examples/jdp.mplstyle
-.JuliaFormatter.toml
+examples/jdp.mplstyle

Project.toml

@@ -1,6 +1,6 @@
 name = "HarmonicBalance"
 uuid = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e"
-authors = ["Jan Kosata <kosataj@phys.ethz.ch>", "Javier del Pino <jdelpino@phys.ethz.ch>"]
+authors = ["Jan Kosata <kosataj@phys.ethz.ch>", "Javier del Pino <jdelpino@phys.ethz.ch>", "Orjan Ameye <orjan.ameye@hotmail.com>"]
 version = "0.9.2"
 
 [deps]
@@ -20,6 +20,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [weakdeps]
@@ -33,26 +34,37 @@ SteadyStateDiffEqExt = "SteadyStateDiffEq"
 TimeEvolution = "OrdinaryDiffEq"
 
 [compat]
+Aqua = "0.8"
 BijectiveHilbert = "0.3.0"
 DSP = "0.7.9"
 DelimitedFiles = "1.9"
 Distances = "0.10.11"
+Documenter = "1.4"
 DocStringExtensions = "0.9.3"
+ExplicitImports = "1.6"
 FFTW = "1.8"
 HomotopyContinuation = "2.9"
 JLD2 = "0.4.48"
+JET = "0.9"
 Latexify = "0.16"
 ModelingToolkit = "9.17"
+NonlinearSolve = "3.12"
 OrderedCollections = "1.6"
-OrdinaryDiffEq = "v6.82.0"
+OrdinaryDiffEq = "v6.82"
 Peaks = "0.5"
 Plots = "1.39"
 ProgressMeter = "1.7.2"
-Symbolics = "5.30"
 SteadyStateDiffEq = "2"
+SymbolicUtils = "2.0"
+Symbolics = "5.30"
+
 julia = "1.10.0"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
+JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
@@ -61,4 +73,4 @@ SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Pkg", "Test", "OrdinaryDiffEq", "ModelingToolkit", "SteadyStateDiffEq", "NonlinearSolve"]
+test = ["Pkg", "Test", "OrdinaryDiffEq", "ModelingToolkit", "SteadyStateDiffEq", "NonlinearSolve", "ExplicitImports", "Aqua", "JET", "Documenter"]

README.md

@@ -4,6 +4,11 @@
 [![docs](https://img.shields.io/badge/docs-online-blue.svg)](https://nonlinearoscillations.github.io/HarmonicBalance.jl/stable/)
 [![Downloads](https://img.shields.io/badge/dynamic/json?url=http%3A%2F%2Fjuliapkgstats.com%2Fapi%2Fv1%2Fmonthly_downloads%2FHarmonicBalance&query=total_requests&suffix=%2Fmonth&label=Downloads)](https://juliapkgstats.com/pkg/HarmonicBalance)
 
+[![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle)
+[![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
+[![JET](https://img.shields.io/badge/%E2%9C%88%EF%B8%8F%20tested%20with%20-%20JET.jl%20-%20red)](https://github.com/aviatesk/JET.jl)
+
+
 **HarmonicBalance.jl** is a Julia package for solving nonlinear differential equations using the method of harmonic balance.
 
 ## Installation

docs/make.jl

@@ -6,48 +6,27 @@ using ModelingToolkit
 using OrdinaryDiffEq
 using SteadyStateDiffEq
 
-makedocs(
-	sitename="HarmonicBalance.jl",
-    modules = [
+include("pages.jl")
+
+makedocs(;
+    sitename="HarmonicBalance.jl",
+    authors="Nonlinear Oscillations Group",
+    modules=[
         HarmonicBalance,
         Base.get_extension(HarmonicBalance, :TimeEvolution),
         Base.get_extension(HarmonicBalance, :ModelingToolkitExt),
-        Base.get_extension(HarmonicBalance, :SteadyStateDiffEqExt)
-        ],
-    warnonly = true,
-	format = Documenter.HTML(
-		mathengine=MathJax2(),
+        Base.get_extension(HarmonicBalance, :SteadyStateDiffEqExt),
+    ],
+    warnonly=true,
+    format=Documenter.HTML(;
+        mathengine=MathJax2(),
         canonical="https://nonlinearoscillations.github.io/HarmonicBalance.jl/stable/",
-		assets = ["assets/favicon.ico", "assets/docs.css"]
+        assets=["assets/favicon.ico", "assets/docs.css"],
         # size_threshold = nothing
-	),
-	pages = [
-		"Background" => Any[
-			"background/harmonic_balance.md"
-			"background/stability_response.md"
-			"background/limit_cycles.md"
-		],
-		"Examples" => Any[
-			"examples/simple_Duffing.md"
-			"examples/linear_response.md"
-			"examples/time_dependent.md"
-			"examples/parametron.md"
-			"examples/limit_cycles.md"
-			],
-		"Manual" => Any[
-			"manual/entering_eom.md"
-			"manual/extracting_harmonics.md"
-			"manual/solving_harmonics.md"
-            "manual/Krylov-Bogoliubov_method.md"
-			"manual/plotting.md"
-			"manual/time_dependent.md"
-			"manual/linear_response.md"
-			"manual/saving.md"
-		]
-	]
+    ),
+    pages=pages,
 )
 
-deploydocs(
-    repo = "github.com/NonlinearOscillations/HarmonicBalance.jl.git",
-    push_preview = false
+deploydocs(;
+    repo="github.com/NonlinearOscillations/HarmonicBalance.jl.git", push_preview=false
 )

docs/pages.jl

@@ -0,0 +1,26 @@
+#! format: off
+pages = [
+		"Background" => Any[
+			"background/harmonic_balance.md"
+			"background/stability_response.md"
+			"background/limit_cycles.md"
+		],
+		"Examples" => Any[
+			"examples/simple_Duffing.md"
+			"examples/linear_response.md"
+			"examples/time_dependent.md"
+			"examples/parametron.md"
+			"examples/limit_cycles.md"
+			],
+		"Manual" => Any[
+			"manual/entering_eom.md"
+			"manual/extracting_harmonics.md"
+			"manual/solving_harmonics.md"
+            "manual/Krylov-Bogoliubov_method.md"
+			"manual/plotting.md"
+			"manual/time_dependent.md"
+			"manual/linear_response.md"
+			"manual/saving.md"
+		]
+	]
+#! format: on

ext/ModelingToolkitExt.jl

@@ -1,22 +1,31 @@
 module ModelingToolkitExt
 
-export ODESystem, ODEProblem
+export ODESystem, ODEProblem, SteadyStateProblem, NonlinearProblem
 
-using HarmonicBalance: HarmonicEquation, is_rearranged,
-                       rearrange_standard, get_variables, simplify, ParameterList
-using ModelingToolkit
-import ModelingToolkit: ODESystem, ODEProblem, NonlinearProblem, SteadyStateProblem
+using HarmonicBalance:
+    HarmonicEquation, is_rearranged, rearrange_standard, get_variables, ParameterList
+using Symbolics: simplify, Equation, substitute, Num, @variables, expand
+using ModelingToolkit:
+    ModelingToolkit,
+    ODESystem,
+    ODEProblem,
+    NonlinearProblem,
+    SteadyStateProblem,
+    varmap_to_vars,
+    parameters,
+    @parameters,
+    @mtkbuild
 
 swapsides(eq::Equation) = Equation(eq.rhs, eq.lhs)
 
 function declare_parameter(var::Num)
     var_sym = Symbol(var)
     new_var = @parameters $var_sym
-    @eval($(var_sym)=first($new_var)) # store the variable under "name" in this namespace
+    @eval($(var_sym) = first($new_var)) # store the variable under "name" in this namespace
     return eval(var_sym)
 end
 
-function ODESystem(eom::HarmonicEquation)
+function ModelingToolkit.ODESystem(eom::HarmonicEquation)
     if !is_rearranged(eom) # check if time-derivatives of the variable are on the right hand side
         eom = rearrange_standard(eom)
     end
@@ -35,32 +44,35 @@ function ODESystem(eom::HarmonicEquation)
     return sys
 end
 
-function ODEProblem(eom::HarmonicEquation, u0, tspan::Tuple,
-        p::ParameterList; in_place = true, kwargs...)
+function ModelingToolkit.ODEProblem(
+    eom::HarmonicEquation, u0, tspan::Tuple, p::ParameterList; in_place=true, kwargs...
+)
     sys = ODESystem(eom)
-    param = ModelingToolkit.varmap_to_vars(p, parameters(sys))
+    param = varmap_to_vars(p, parameters(sys))
     if !in_place # out-of-place
-        prob = ODEProblem{false}(sys, u0, tspan, param; jac = true, kwargs...)
+        prob = ODEProblem{false}(sys, u0, tspan, param; jac=true, kwargs...)
     else # in-place
-        prob = ODEProblem{true}(sys, u0, tspan, param; jac = true, kwargs...)
+        prob = ODEProblem{true}(sys, u0, tspan, param; jac=true, kwargs...)
     end
     return prob
 end
 
-function NonlinearProblem(
-        eom::HarmonicEquation, u0, p::ParameterList; in_place = true, kwargs...)
-    ss_prob = SteadyStateProblem(eom, u0, p::ParameterList; in_place = in_place, kwargs...)
+function ModelingToolkit.NonlinearProblem(
+    eom::HarmonicEquation, u0, p::ParameterList; in_place=true, kwargs...
+)
+    ss_prob = SteadyStateProblem(eom, u0, p::ParameterList; in_place=in_place, kwargs...)
     return NonlinearProblem(ss_prob) # gives warning of some internal deprication
 end
 
-function SteadyStateProblem(
-        eom::HarmonicEquation, u0, p::ParameterList; in_place = true, kwargs...)
+function ModelingToolkit.SteadyStateProblem(
+    eom::HarmonicEquation, u0, p::ParameterList; in_place=true, kwargs...
+)
     sys = ODESystem(eom)
-    param = ModelingToolkit.varmap_to_vars(p, parameters(sys))
+    param = varmap_to_vars(p, parameters(sys))
     if !in_place # out-of-place
-        prob = SteadyStateProblem{false}(sys, u0, param; jac = true, kwargs...)
+        prob = SteadyStateProblem{false}(sys, u0, param; jac=true, kwargs...)
     else # in-place
-        prob = SteadyStateProblem{true}(sys, u0, param; jac = true, kwargs...)
+        prob = SteadyStateProblem{true}(sys, u0, param; jac=true, kwargs...)
     end
     return prob
 end

ext/SteadyStateDiffEqExt.jl

@@ -1,14 +1,15 @@
 module SteadyStateDiffEqExt
 
 export steady_state_sweep
-import HarmonicBalance: steady_state_sweep
+using HarmonicBalance: HarmonicBalance, steady_state_sweep
+
 using SteadyStateDiffEq: solve, NonlinearProblem, SteadyStateProblem, DynamicSS, remake
 using LinearAlgebra: norm, eigvals
 using SteadyStateDiffEq.SciMLBase.SciMLStructures: isscimlstructure, Tunable, replace
 
-function steady_state_sweep(
-        prob::SteadyStateProblem, alg::DynamicSS;
-        varied::Pair, kwargs...)
+function HarmonicBalance.steady_state_sweep(
+    prob::SteadyStateProblem, alg::DynamicSS; varied::Pair, kwargs...
+)
     varied_idx, sweep_range = varied
 
     # if p is dual number (AD) result is dual number
@@ -18,16 +19,20 @@ function steady_state_sweep(
         u0 = i == 1 ? [0.0, 0.0] : result[i - 1]
         # make type-stable: FD.Dual or Float64
         parameters = get_new_parameters(prob, varied_idx, value)
-        sol = solve(remake(prob, p = parameters, u0 = u0), alg; kwargs...)
+        sol = solve(remake(prob; p=parameters, u0=u0), alg; kwargs...)
         result[i] = sol.u
     end
     return result
 end
 
-function steady_state_sweep(
-        prob_np::NonlinearProblem, prob_ss::SteadyStateProblem,
-        alg_np, alg_ss::DynamicSS;
-        varied::Pair, kwargs...)
+function HarmonicBalance.steady_state_sweep(
+    prob_np::NonlinearProblem,
+    prob_ss::SteadyStateProblem,
+    alg_np,
+    alg_ss::DynamicSS;
+    varied::Pair,
+    kwargs...,
+)
     varied_idx, sweep_range = varied
     # if p is dual number (AD) result is dual number
     result = [similar(prob_np.u0) for _ in sweep_range]
@@ -36,17 +41,18 @@ function steady_state_sweep(
         u0 = i == 1 ? Base.zeros(length(prob_np.u0)) : result[i - 1]
         # make type-stable: FD.Dual or Float64
         parameters = get_new_parameters(prob_np, varied_idx, value)
-        sol_nn = solve(remake(prob_np, p = parameters, u0 = u0), alg_np; kwargs...)
+        sol_nn = solve(remake(prob_np; p=parameters, u0=u0), alg_np; kwargs...)
 
         # last argument is time but does not matter
         param_val = tunable_parameters(parameters)
         zeros = norm(prob_np.f.f.f.f.f_oop(sol_nn.u, param_val, 0))
         jac = prob_np.f.jac.f.f.f_oop(sol_nn.u, param_val, 0)
         eigval = jac isa Vector ? jac : eigvals(jac) # eigvals favourable supports FD.Dual
 
-        if !isapprox(zeros, 0, atol = 1e-5) || any(λ -> λ > 0, real.(eigval))
-            sol_ss = solve(remake(prob_ss, p = parameters, u0 = u0),
-                alg_ss, abstol = 1e-5, reltol = 1e-5)
+        if !isapprox(zeros, 0; atol=1e-5) || any(λ -> λ > 0, real.(eigval))
+            sol_ss = solve(
+                remake(prob_ss; p=parameters, u0=u0), alg_ss; abstol=1e-5, reltol=1e-5
+            )
             result[i] = sol_ss.u
         else
             result[i] = sol_nn.u
@@ -57,7 +63,7 @@ function steady_state_sweep(
 end
 
 function tunable_parameters(param)
-    hasfield(typeof(param), :tunable) ? param.tunable[1] : param
+    return hasfield(typeof(param), :tunable) ? param.tunable[1] : param
 end
 
 function get_new_parameters(prob, varied_idx, value)
@@ -66,16 +72,17 @@ function get_new_parameters(prob, varied_idx, value)
         rest = [prob.p.discrete, prob.p.nonnumeric, prob.p.dependent, prob.p.constant]
 
         all(isempty.(rest)) || error("Only tunable parameters are supported")
-        length(prob.p.tunable) == 1 ||
-            error("The type of the parameters should be uniform")
+        length(prob.p.tunable) == 1 || error("The type of the parameters should be uniform")
 
         old_parameters_values = prob.p.tunable[1]
-        parameter_values = eltype(old_parameters_values)[i == varied_idx ? value : x
-                                          for (i, x) in enumerate(old_parameters_values)]
+        parameter_values = eltype(old_parameters_values)[
+            i == varied_idx ? value : x for (i, x) in enumerate(old_parameters_values)
+        ]
         parameters = replace(Tunable(), prob.p, parameter_values)
     else
-        parameters = eltype(prob.p)[i == varied_idx ? value : x
-                                    for (i, x) in enumerate(prob.p)]
+        parameters = eltype(prob.p)[
+            i == varied_idx ? value : x for (i, x) in enumerate(prob.p)
+        ]
     end
     return parameters
 end