feat: Make Plots into an Extension (breaking) (#371)

NonlinearOscillations · Jan 29, 2025 · e224717 · e224717
1 parent 3908971
commit e224717
Show file tree

Hide file tree

Showing 29 changed files with 813 additions and 725 deletions.
diff --git a/Project.toml b/Project.toml
@@ -11,25 +11,26 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
-Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [weakdeps]
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 
 [extensions]
 ModelingToolkitExt = "ModelingToolkit"
+PlotsExt = "Plots"
 SteadyStateDiffEqExt = "SteadyStateDiffEq"
 TimeEvolution = "OrdinaryDiffEqTsit5"
 
@@ -41,15 +42,15 @@ DocStringExtensions = "0.9.3"
 FunctionWrappers = "1.1.3"
 HomotopyContinuation = "2.12"
 JLD2 = "0.4.48, 0.5"
-Latexify = "0.16"
 ModelingToolkit = "9.60"
-OrdinaryDiffEqTsit5 = "1.1"
+NaNMath = "1.1"
 OrderedCollections = "1.6"
+OrdinaryDiffEqTsit5 = "1.1"
 Plots = "1.40.9"
 PrecompileTools = "1.2"
 ProgressMeter = "1.7.2"
+SciMLBase = "2.72.0"
 SteadyStateDiffEq = "2.4.1"
 SymbolicUtils = "3.11"
 Symbolics = "6.23"
 julia = "1.10.0"
-NaNMath = "1.1"
diff --git a/README.md b/README.md
@@ -31,7 +31,7 @@ Let's find the steady states of a driven Duffing oscillator with nonlinear dampi
 <img src="/docs/images/github_readme_eq.png" width="450">
 
 ```julia
-using HarmonicBalance
+using HarmonicBalance, Plots
 @variables α ω ω0 F η t x(t) # declare constant variables and a function x(t)
 diff_eq = DifferentialEquation(d(x,t,2) + ω0^2*x + α*x^3 + η*d(x,t)*x^2 ~ F*cos(ω*t), x)
 add_harmonic!(diff_eq, x, ω) # specify the ansatz x = u(T) cos(ωt) + v(T) sin(ωt)

diff --git a/docs/make.jl b/docs/make.jl
@@ -21,6 +21,7 @@ bib = CitationBibliography(
 )
 
 using Plots
+PlotsExt = Base.get_extension(HarmonicBalance, :PlotsExt)
 default(; fmt=:png)
 # Gotta set this environment variable when using the GR run-time on CI machines.
 # This happens as examples will use Plots.jl to make plots and movies.
@@ -40,6 +41,7 @@ makedocs(;
         ModelingToolkitExt,
         SteadyStateDiffEqExt,
         HarmonicBalance.LinearResponse,
+        PlotsExt,
     ],
     format=DocumenterVitepress.MarkdownVitepress(;
         repo="github.com/NonlinearOscillations/HarmonicBalance.jl",

diff --git a/docs/src/examples/parametron.md b/docs/src/examples/parametron.md
@@ -14,7 +14,7 @@ where for completeness we also considered an external drive term $F_\text{d}(t)=
 To implement this system in Harmonic Balance, we first import the library
 
 ````@example parametron
-using HarmonicBalance
+using HarmonicBalance, Plots
 ````
 
 Subsequently, we type define parameters in the problem and the oscillating amplitude function $x(t)$ using the `variables` macro from `Symbolics.jl`

diff --git a/docs/src/examples/state_dependent_perturbation.md b/docs/src/examples/state_dependent_perturbation.md
@@ -228,7 +228,7 @@ harmonic_tmp.equations = HB.Symbolics.substitute(
     HB.rearrange_standard(harmonic_normal).equations[1:2], Dict(u2 => ua, v2 => va)
 )
 harmonic_tmp.parameters = push!(harmonic_tmp.parameters, ua, va)
-prob = HarmonicBalance.Problem(harmonic_tmp)
+prob = HarmonicBalance.Problem(harmonic_tmp, varied, fixed)
 ````
 
 We will sweep over the $\omega-\lambda$ plane and substitute the non-zero amplitude solution of the antisymmetric mode into the coupled equations of thesymmetric mode.
@@ -267,8 +267,9 @@ function solve_perturbed_system(prob, input)
     rounded_solutions = HB.HomotopyContinuation.solutions.(result_full)
     solutions = HB.pad_solutions(rounded_solutions)
 
+    jacobian = HarmonicBalance.get_Jacobian(harmonic_tmp)
     J_variables = cat(prob.variables, collect(keys(varied)), [ua, va]; dims=1)
-    compiled_J = HB.compile_matrix(prob.jacobian, J_variables; rules=fixed)
+    compiled_J = HB.compile_matrix(jacobian, J_variables; rules=fixed)
     compiled_J = HB.JacobianFunction(HB.solution_type(solutions))(compiled_J)
     result = HB.Result(
         solutions,

diff --git a/docs/src/introduction/index.md b/docs/src/introduction/index.md
@@ -26,7 +26,7 @@ Let us find the steady states of an external driven Duffing oscillator with nonl
 ```
 
 ```@example getting_started
-using HarmonicBalance
+using HarmonicBalance, Plots
 @variables α ω ω0 F t η x(t) # declare constant variables and a function x(t)
 eom = d(x,t,2) + ω0^2*x + α*x^3 + η*d(x,t)*x^2 ~ F*cos(ω*t)
 diff_eq = DifferentialEquation(eom, x)

diff --git a/docs/src/manual/linear_response.md b/docs/src/manual/linear_response.md
@@ -14,11 +14,7 @@ HarmonicBalance.get_Jacobian
 
 ## Linear response
 
-The response to white noise can be shown with `plot_linear_response`. Depending on the `order` argument, different methods are used. 
-
-```@docs; canonical=false
-HarmonicBalance.LinearResponse.plot_linear_response
-```
+The response to white noise can be shown with `plot_linear_response`. Depending on the `order` argument, different methods are used.
 
 ### First order
 
@@ -29,6 +25,8 @@ The advantage of this method is that for a given parameter set, only one matrix
 Behind the scenes, the spectra are stored using the dedicated structs `Lorentzian` and `JacobianSpectrum`.
 
 ```@docs; canonical=false
+HarmonicBalance.LinearResponse.get_jacobian_response
+HarmonicBalance.LinearResponse.get_rotframe_jacobian_response
 HarmonicBalance.LinearResponse.JacobianSpectrum
 HarmonicBalance.LinearResponse.Lorentzian
 ```
@@ -38,6 +36,7 @@ HarmonicBalance.LinearResponse.Lorentzian
 Setting `order > 1` increases the accuracy of the response spectra. However, unlike for the Jacobian, here we must perform a matrix inversion for each response frequency.  
 
 ```@docs; canonical=false
+HarmonicBalance.LinearResponse.get_linear_response
 HarmonicBalance.LinearResponse.ResponseMatrix
 HarmonicBalance.LinearResponse.get_response
 HarmonicBalance.LinearResponse.get_response_matrix

diff --git a/docs/src/manual/plotting.md b/docs/src/manual/plotting.md
@@ -12,12 +12,10 @@ The function `plot` is multiple-dispatched to plot 1D and 2D datasets.
 In 1D, the solutions are colour-coded according to the branches obtained by `sort_solutions`. 
 
 ```@docs; canonical=false
-HarmonicBalance.plot(::HarmonicBalance.Result, varags...)
-HarmonicBalance.plot!
+plot(::HarmonicBalance.Result, varags...)
+plot!
 ```
 
-
-
 ## Plotting phase diagrams
 
 In many problems, rather than in any property of the solutions themselves, we are interested in the phase diagrams, encoding the number of (stable) solutions in different regions of the parameter space. `plot_phase_diagram` handles this for 1D and 2D datasets.

diff --git a/docs/src/manual/time_dependent.md b/docs/src/manual/time_dependent.md
@@ -11,7 +11,7 @@ AdiabaticSweep
 ## Plotting
 
 ```@docs; canonical=false
-HarmonicBalance.plot(::OrdinaryDiffEqTsit5.ODESolution, ::Any, ::HarmonicEquation)
+plot(::OrdinaryDiffEqTsit5.ODESolution, ::Any, ::HarmonicEquation)
 ```
 
 ## Miscellaneous

diff --git a/docs/src/tutorials/classification.md b/docs/src/tutorials/classification.md
@@ -3,7 +3,7 @@
 Given that you obtained some steady states for a parameter sweep of a specific model it can be useful to classify these solution. Let us consider a simple pametric oscillator 
 
 ```@example classification
-using HarmonicBalance
+using HarmonicBalance, Plots
 
 @variables ω₀ γ λ α ω t x(t)
 

diff --git a/docs/src/tutorials/limit_cycles.md b/docs/src/tutorials/limit_cycles.md
@@ -6,7 +6,7 @@ In contrast to the previous tutorials, limit cycle problems feature harmonic(s)
 
 Here we solve the equation of motion of the [van der Pol oscillator](https://en.wikipedia.org/wiki/Van_der_Pol_oscillator). This is a single-variable second-order ODE with continuous time-translation symmetry (i.e., no 'clock' imposing a frequency and/or phase), which displays periodic solutions known as _relaxation oscillations_. For more detail, refer also to [arXiv:2308.06092](https://arxiv.org/abs/2308.06092).
 ```@example lc
-using HarmonicBalance
+using HarmonicBalance, Plots
 @variables ω_lc, t, ω0, x(t), μ
 diff_eq = DifferentialEquation(d(d(x,t),t) - μ*(1-x^2) * d(x,t) + x, x)
 ```

diff --git a/docs/src/tutorials/steady_states.md b/docs/src/tutorials/steady_states.md
@@ -14,7 +14,7 @@ which constraints the spectrum of ``x(t)`` to a single harmonic. Fixing the quad
 
 First we need to declare the symbolic variables (the excellent [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl) is used here).
 ```@example steady_state
-using HarmonicBalance
+using HarmonicBalance, Plots
 @variables α ω ω0 F γ t x(t) # declare constant variables and a function x(t)
 ```
 Next, we have to input the equations of motion. This will be stored as a `DifferentialEquation`. The input needs to specify that only `x` is a mathematical variable, the other symbols are parameters:

diff --git a/docs/src/tutorials/time_dependent.md b/docs/src/tutorials/time_dependent.md
@@ -21,7 +21,7 @@ Here we primarily demonstrate on the parametrically driven oscillator.
 
 We start by defining our system.
 ```@example time_dependent
-using HarmonicBalance
+using HarmonicBalance, Plots
 @variables ω0 γ λ F θ η α ω t x(t)
 
 eq =  d(d(x,t),t) + γ*d(x,t) + ω0^2*(1 - λ*cos(2*ω*t))*x + α*x^3 + η*d(x,t)*x^2 ~ F*cos(ω*t + θ)

diff --git a/examples/parametron.jl b/examples/parametron.jl
@@ -9,7 +9,7 @@
 
 # To implement this system in Harmonic Balance, we first import the library
 
-using HarmonicBalance
+using HarmonicBalance, Plots
 
 # Subsequently, we type define parameters in the problem and the oscillating amplitude function $x(t)$ using the `variables` macro from `Symbolics.jl`
 

diff --git a/ext/PlotsExt/PlotsExt.jl b/ext/PlotsExt/PlotsExt.jl
@@ -0,0 +1,36 @@
+module PlotsExt
+
+using DocStringExtensions
+using Plots
+using HarmonicBalance:
+    HarmonicBalance,
+    Result,
+    branch_count,
+    transform_solutions,
+    _apply_mask,
+    _get_mask,
+    _realify,
+    dim,
+    get_class,
+    get_variable_solutions
+using Symbolics: Num
+using LinearAlgebra: LinearAlgebra, eigvals, eigvecs
+
+using Plots: heatmap, theme_palette, scatter, RGB, cgrad
+using Symbolics.Latexify: Latexify, latexify, @L_str
+using Symbolics.Latexify.LaTeXStrings: LaTeXStrings
+using SciMLBase: SciMLBase
+
+const _set_Plots_default = Dict{Symbol,Any}([
+    :fontfamily => "computer modern",
+    :titlefont => "computer modern",
+    :tickfont => "computer modern",
+    :linewidth => 2,
+    :legend_position => :outerright,
+])
+
+include("steady_states.jl")
+include("linear_response.jl")
+include("time_evolution.jl")
+
+end #module