HarmonicBalance Extentions (#154)

* pffffff

Revert "pffffff"

This reverts commit a27cadd.

* Make TimeEvolution module an extention

* reduce namespace

* add ModelingToolkit Ext

* add SteadyStateDiffEqExt

* Format files

* raname file for extention functions

* support for 1.9

support for 1.9
try compatibility with julia 1.8
update test.yml to v4
support for julia < 1.9
add OrdinaryDiffEq still as [deps] such that it can be used in 1.8

* bump project version to v0.9 (breaking)

.github/dependabot.yml

@@ -0,0 +1,7 @@
+# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates
+version: 2
+updates:
+  - package-ecosystem: "github-actions"
+    directory: "/" # Location of package manifests
+    schedule:
+      interval: "weekly"

.github/workflows/CI.yml

@@ -0,0 +1,54 @@
+name: CI
+on:
+  push:
+    tags:
+      - '*'
+    branches:
+      - 'master'
+  pull_request:
+    branches:
+      - 'master'
+concurrency:
+  # Skip intermediate builds: always.
+  # Cancel intermediate builds: only if it is a pull request build.
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
+jobs:
+  test:
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}
+    runs-on: ${{ matrix.os }}
+    timeout-minutes: 60
+    permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created
+      actions: write
+      contents: read
+    strategy:
+      matrix:
+        version:
+          - '1.10'
+          - '1.9'
+          # - 'nightly'
+        os:
+          - ubuntu-latest
+        arch:
+          - x64
+    steps:
+      - uses: actions/checkout@v4
+        with:
+          persist-credentials: false
+      - uses: julia-actions/setup-julia@v1
+        with:
+          show-versioninfo: true
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: julia-actions/cache@v1
+      - uses: julia-actions/julia-buildpkg@latest
+      - uses: julia-actions/julia-runtest@latest
+      # - uses: julia-actions/julia-processcoverage@v1
+      # - uses: codecov/codecov-action@v4
+      #   with:
+      #     files: lcov.info
+      #     token: ${{ secrets.CODECOV_TOKEN }}
+      #     fail_ci_if_error: false
+      # - uses: julia-actions/julia-uploadcoveralls@v1
+      #   env:
+      #     COVERALLS_TOKEN: ${{ secrets.COVERALLS_TOKEN }}

.github/workflows/CompatHelper.yml

@@ -3,43 +3,14 @@ on:
   schedule:
     - cron: 0 0 * * *
   workflow_dispatch:
-permissions:
-  contents: write
-  pull-requests: write
 jobs:
   CompatHelper:
     runs-on: ubuntu-latest
     steps:
-      - name: Check if Julia is already available in the PATH
-        id: julia_in_path
-        run: which julia
-        continue-on-error: true
-      - name: Install Julia, but only if it is not already available in the PATH
-        uses: julia-actions/setup-julia@v1
-        with:
-          version: '1'
-          arch: ${{ runner.arch }}
-        if: steps.julia_in_path.outcome != 'success'
-      - name: "Add the General registry via Git"
-        run: |
-          import Pkg
-          ENV["JULIA_PKG_SERVER"] = ""
-          Pkg.Registry.add("General")
-        shell: julia --color=yes {0}
-      - name: "Install CompatHelper"
-        run: |
-          import Pkg
-          name = "CompatHelper"
-          uuid = "aa819f21-2bde-4658-8897-bab36330d9b7"
-          version = "3"
-          Pkg.add(; name, uuid, version)
-        shell: julia --color=yes {0}
-      - name: "Run CompatHelper"
-        run: |
-          import CompatHelper
-          CompatHelper.main()
-        shell: julia --color=yes {0}
+      - name: Pkg.add("CompatHelper")
+        run: julia -e 'using Pkg; Pkg.add("CompatHelper")'
+      - name: CompatHelper.main()
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
           COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }}
-          # COMPATHELPER_PRIV: ${{ secrets.COMPATHELPER_PRIV }}
+        run: julia -e 'using CompatHelper; CompatHelper.main()'

.gitignore

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

Project.toml

@@ -1,7 +1,7 @@
 name = "HarmonicBalance"
 uuid = "e13b9ff6-59c3-11ec-14b1-f3d2cc6c135e"
 authors = ["Jan Kosata <[email protected]>", "Javier del Pino <[email protected]>"]
-version = "0.8.0"
+version = "0.9.0"
 
 [deps]
 BijectiveHilbert = "91e7fc40-53cd-4118-bd19-d7fcd1de2a54"
@@ -15,15 +15,25 @@ JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Peaks = "18e31ff7-3703-566c-8e60-38913d67486b"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
+[weakdeps]
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
+
+[extensions]
+ModelingToolkitExt = "ModelingToolkit"
+SteadyStateDiffEqExt = "SteadyStateDiffEq"
+TimeEvolution = "OrdinaryDiffEq"
+
 [compat]
+ModelingToolkit = "< 9"
 BijectiveHilbert = "0.3.0"
 DSP = "0.7.4"
 DelimitedFiles = "1"
@@ -39,11 +49,15 @@ Peaks = "0.4.0"
 Plots = "1.35.0"
 ProgressMeter = "1.7.2"
 Symbolics = "5.0.0"
-julia = "1.8.2"
+julia = "1.9.0"
 
 [extras]
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Pkg", "Test"]
+test = ["Pkg", "Test", "OrdinaryDiffEq", "ModelingToolkit", "SteadyStateDiffEq", "NonlinearSolve"]

ext/ModelingToolkitExt.jl

@@ -0,0 +1,67 @@
+module ModelingToolkitExt
+
+export ODESystem, ODEProblem
+
+using HarmonicBalance: HarmonicEquation, is_rearranged,
+                       rearrange_standard, get_variables, simplify, ParameterList
+using ModelingToolkit
+import ModelingToolkit: ODESystem, ODEProblem, NonlinearProblem, SteadyStateProblem
+
+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
+    return eval(var_sym)
+end
+
+function 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
+
+    slow_time = (@variables T; T)
+    par_names = declare_parameter.(eom.parameters)
+    vars = get_variables(eom)
+
+    eqs = deepcopy(eom.equations)
+    eqs = swapsides.(eqs)
+    eqs = simplify.(expand.(eqs))
+
+    # compute jacobian for performance
+    @named sys = ODESystem(eqs, slow_time, vars, par_names) #mtk v9 need @mtkbuild
+    return sys
+end
+
+function ODEProblem(eom::HarmonicEquation, u0, tspan::Tuple,
+        p::ParameterList; in_place = true, kwargs...)
+    sys = ODESystem(eom)
+    param = ModelingToolkit.varmap_to_vars(p, parameters(sys))
+    if !in_place # out-of-place
+        prob = ODEProblem{false}(sys, u0, tspan, param; jac = true, kwargs...)
+    else # in-place
+        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...)
+    return NonlinearProblem(ss_prob) # gives warning of some internal deprication
+end
+
+function SteadyStateProblem(
+        eom::HarmonicEquation, u0, p::ParameterList; in_place = true, kwargs...)
+    sys = ODESystem(eom)
+    param = ModelingToolkit.varmap_to_vars(p, parameters(sys))
+    if !in_place # out-of-place
+        prob = SteadyStateProblem{false}(sys, u0, param; jac = true, kwargs...)
+    else # in-place
+        prob = SteadyStateProblem{true}(sys, u0, param; jac = true, kwargs...)
+    end
+    return prob
+end
+
+end # module

ext/SteadyStateDiffEqExt.jl

@@ -0,0 +1,59 @@
+module SteadyStateDiffEqExt
+
+export steady_state_sweep
+import HarmonicBalance: steady_state_sweep
+using SteadyStateDiffEq: solve, NonlinearProblem, SteadyStateProblem, DynamicSS, remake
+using LinearAlgebra: norm, eigvals
+
+function 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
+    result = [similar(prob.p, 2) for _ in sweep_range]
+
+    foreach(pairs(sweep_range)) do (i, value)
+        u0 = i == 1 ? [0.0, 0.0] : result[i - 1]
+        # make type-stable: FD.Dual or Float64
+        parameters = eltype(prob_np.p)[i == varied_idx ? value : x
+                                       for (i, x) in enumerate(prob_np.p)]
+        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...)
+    varied_idx, sweep_range = varied
+    # if p is dual number (AD) result is dual number
+    result = [similar(prob_np.p) for _ in sweep_range]
+
+    foreach(pairs(sweep_range)) do (i, value)
+        u0 = i == 1 ? Base.zeros(length(prob_np.u0)) : result[i - 1]
+        # make type-stable: FD.Dual or Float64
+        parameters = eltype(prob_np.p)[i == varied_idx ? value : x
+                                       for (i, x) in enumerate(prob_np.p)]
+        sol_nn = solve(remake(prob_np, p = parameters, u0 = u0), alg_np; kwargs...)
+
+        # last argument is time but does not matter
+        zeros = norm(prob_np.f.f.f.f.f_oop(sol_nn.u, parameters, 0))
+        jac = prob_np.f.jac.f.f.f_oop(sol_nn.u, parameters, 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)
+            result[i] = sol_ss.u
+        else
+            result[i] = sol_nn.u
+        end
+    end
+
+    return result
+end
+
+end # module