Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bifurcation kit extension #2300

Closed
wants to merge 63 commits into from
Closed

Bifurcation kit extension #2300

wants to merge 63 commits into from

Conversation

TorkelE
Copy link
Member

@TorkelE TorkelE commented Oct 4, 2023

This extension creates dispatches for BifurcationKit.jl's BifurcationProblem when used on NonlinearSystem or ODESystems.

This should work now:

using ModelingToolkit
using BifurcationKit

# Define ODEFunction
begin
    @variables t x(t) y(t)
    @parameters μ α
    eqs = [0 ~ μ*x - x^3 + α*y,
           0 ~ -y]
    @named nsys = NonlinearSystem(eqs, [x, y], [μ, α])
end

# Define other bifurcation inputs
begin
    p_span = (-4.0, 6.0)
    bif_par = μ
    plot_var = x
    p_start ==> -1.0, α => 1.0]
    u0_guess = [x => 1.0, y => 1.0]
end

# Creates a BifurcationProblem.
# This dispatch is what is implemented here, then it is BK stuff all the way.
bprob = BifurcationProblem(nsys, u0_guess, p_start, bif_par; plot_var=plot_var, jac=false)

# Sets bifurcation parameters
opt_newton = NewtonPar(tol = 1e-9, max_iterations = 20)
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = 0.01,
	max_steps = 100, nev = 2, newton_options = opt_newton,
	p_max = p_span[2], p_min = p_span[1],
	detect_bifurcation = 3, n_inversion = 4, tol_bisection_eigenvalue = 1e-8, dsmin_bisection = 1e-9)

# Creates the bifurcation diagram.
bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true)

# You can plot the diagram like
using Plots
plot(bf; putspecialptlegend=false, markersize=2, plotfold=false, title = "#branches = $(size(bf))")

image

Tests and docs not added yet, will do if people think this approach looks good.

# When input is a NonlinearSystem.
function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem, u0_guess, p_start, bif_par, args...; plot_var=nothing, record_from_solution=BifurcationKit.record_sol_default, kwargs...)
# Creates F and J functions.
ofun = NonlinearFunction(nsys; jac=true)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is the jac required? The user should request it if it's not.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I will make it optional.

@ChrisRackauckas
Copy link
Member

This looks like a great idea. Just needs tests and docs.

@TorkelE
Copy link
Member Author

TorkelE commented Oct 6, 2023

Sounds good, there is a bug in this for for systems with more species/parameters that I need to fix. I might also drop the p_span input (will discuss it with Romain). Once I'm sure of this I will make an updated version, and test+docs (and then make an extension in Catalyst that builds on this one).

ChrisRackauckas and others added 6 commits October 6, 2023 12:03
The recursive structure of SystemStructure was giving precompilation an issue so I removed that module. Then I setup precompilation on the basic interface.  Test case:

```julia
@time using ModelingToolkit
@time using DifferentialEquations
@time using Plots
@time begin
    @variables t x(t)
    D = Differential(t)
    @nAmed sys = ODESystem([D(x) ~ -x])
end;
@time prob = ODEProblem(structural_simplify(sys), [x => 30.0], (0, 100), [], jac = true);
@time sol = solve(prob);
@time plot(sol, idxs=[x]);
```

Before:

```
11.082586 seconds (19.32 M allocations: 1.098 GiB, 4.09% gc time, 0.71% compilation time: 87% of which was recompilation)
0.639738 seconds (661.39 k allocations: 101.321 MiB, 4.33% gc time, 6.46% compilation time)
3.703724 seconds (5.71 M allocations: 322.840 MiB, 5.22% gc time, 9.92% compilation time: 86% of which was recompilation)
7.795297 seconds (8.25 M allocations: 483.041 MiB, 2.50% gc time, 99.88% compilation time)
21.719376 seconds (44.11 M allocations: 2.485 GiB, 5.68% gc time, 99.48% compilation time)
2.602250 seconds (4.04 M allocations: 253.058 MiB, 4.60% gc time, 99.90% compilation time)
2.450509 seconds (5.17 M allocations: 332.101 MiB, 5.89% gc time, 99.41% compilation time: 30% of which was recompilation)
```

After:

```
9.129141 seconds (22.77 M allocations: 1.291 GiB, 4.65% gc time, 0.62% compilation time: 87% of which was recompilation)
0.784464 seconds (667.59 k allocations: 101.524 MiB, 3.95% gc time, 4.16% compilation time)
3.111142 seconds (5.42 M allocations: 305.594 MiB, 3.82% gc time, 6.39% compilation time: 82% of which was recompilation)
0.105567 seconds (157.39 k allocations: 10.522 MiB, 8.81% gc time, 95.49% compilation time: 74% of which was recompilation)
1.993642 seconds (4.03 M allocations: 218.310 MiB, 2.69% gc time, 96.95% compilation time: 82% of which was recompilation)
1.806758 seconds (4.06 M allocations: 254.371 MiB, 4.44% gc time, 99.91% compilation time)
1.694666 seconds (5.27 M allocations: 339.088 MiB, 6.18% gc time, 99.39% compilation time: 31% of which was recompilation)
```

And that's on v1.9, so v1.10 should be even better. 20 seconds off hehe.
Fix recursive structure and setup precompilation
- All sections, except building hierarchal models, use `@mtkmodel`.
- "Non-DSL way of defining an ODESystem" section is added.
@TorkelE
Copy link
Member Author

TorkelE commented Oct 14, 2023

Added some docs and a test. If things seems to work out I could add a couple of more tests.

@@ -0,0 +1,40 @@
module MTKBifurcationKitExt
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fix the naming.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what should I name it?

@ChrisRackauckas
Copy link
Member

Tests fail because the Project.toml needs to be updated.

@ChrisRackauckas
Copy link
Member

That rebase did not go well

@TorkelE
Copy link
Member Author

TorkelE commented Oct 17, 2023

Yeah, no idea what happened. Just tried putting it back at master.

@ChrisRackauckas
Copy link
Member

Well I think it may be ready to go if it's cleaned up.

@TorkelE
Copy link
Member Author

TorkelE commented Oct 17, 2023

I will probably just make a new PR, easier that way.

@TorkelE
Copy link
Member Author

TorkelE commented Oct 17, 2023

Replaced by: #2321

@TorkelE TorkelE closed this Oct 17, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants