Skip to content

Commit

Permalink
up
Browse files Browse the repository at this point in the history
  • Loading branch information
TorkelE committed May 29, 2024
1 parent 8fd8cf6 commit 635e824
Show file tree
Hide file tree
Showing 11 changed files with 200 additions and 201 deletions.
22 changes: 11 additions & 11 deletions docs/pages.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,10 +15,10 @@ pages = Any[
#"model_creation/parametric_stoichiometry.md",# Distributed parameters, rates, and initial conditions.
# Loading and writing models to files.
# Model visualisation.
"model_creation/network_analysis.md",
#"model_creation/network_analysis.md",
"model_creation/chemistry_related_functionality.md",
"Model creation examples" => Any[
"model_creation/examples/basic_CRN_library.md",
#"model_creation/examples/basic_CRN_library.md",
"model_creation/examples/programmatic_generative_linear_pathway.md",
#"model_creation/examples/hodgkin_huxley_equation.md",
#"model_creation/examples/smoluchowski_coagulation_equation.md"
Expand All @@ -29,20 +29,20 @@ pages = Any[
# Simulation introduction.
"model_simulation/simulation_plotting.md",
"model_simulation/simulation_structure_interfacing.md",
"model_simulation/ensemble_simulations.md",
#"model_simulation/ensemble_simulations.md",
# Stochastic simulation statistical analysis.
"model_simulation/ode_simulation_performance.md",
# ODE Performance considerations/advice.
# SDE Performance considerations/advice.
# Jump Performance considerations/advice.
# Finite state projection
# "model_simulation/ode_simulation_performance.md",
# # ODE Performance considerations/advice.
# # SDE Performance considerations/advice.
# # Jump Performance considerations/advice.
# # Finite state projection
],
"Steady state analysis" => Any[
"steady_state_functionality/homotopy_continuation.md",
# "steady_state_functionality/homotopy_continuation.md",
"steady_state_functionality/nonlinear_solve.md",
"steady_state_functionality/steady_state_stability_computation.md",
"steady_state_functionality/bifurcation_diagrams.md",
"steady_state_functionality/dynamical_systems.md"
# "steady_state_functionality/bifurcation_diagrams.md",
# "steady_state_functionality/dynamical_systems.md"
],
"Inverse Problems" => Any[
# Inverse problems introduction.
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -80,8 +80,8 @@ The `@reaction_network` command is followed by the `begin` keyword, which is fol
Here, we create a simple [*birth-death* model](@ref basic_CRN_library_bd), where a single species ($X$) is created at rate $b$, and degraded at rate $d$. The model is stored in the variable `rn`.
```@example ex2
rn = @reaction_network begin
b, 0 --> X
d, X --> 0
b, 0 --> X
d, X --> 0
end
```
For more information on how to use the Catalyst model creator (also known as *the Catalyst DSL*), please read [the corresponding documentation](https://docs.sciml.ai/Catalyst/stable/catalyst_functionality/dsl_description/).
Expand Down Expand Up @@ -152,8 +152,8 @@ Each reaction is also associated with a specific rate (corresponding to a parame
We declare the model using the `@reaction_network` macro, and store it in the `sir_model` variable.
```@example ex2
sir_model = @reaction_network begin
b, S + I --> 2I
k, I --> R
b, S + I --> 2I
k, I --> R
end
```
Note that the first reaction contains two different substrates (separated by a `+` sign). While there is only a single product (*I*), two copies of *I* are produced. The *2* in front of the product *I* denotes this.
Expand Down
24 changes: 12 additions & 12 deletions docs/src/inverse_problems/examples/ode_fitting_oscillation.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
In this example we will use [Optimization.jl](https://github.com/SciML/Optimization.jl) to fit the parameters of an oscillatory system (the Brusselator) to data. Here, special consideration is taken to avoid reaching a local minimum. Instead of fitting the entire time series directly, we will start with fitting parameter values for the first period, and then use those as an initial guess for fitting the next (and then these to find the next one, and so on). Using this procedure is advantageous for oscillatory systems, and enables us to reach the global optimum.

First, we fetch the required packages.
```@example pe1
```@example pe_osc_example
using Catalyst
using OrdinaryDiffEq
using Optimization
Expand All @@ -11,7 +11,7 @@ using SciMLSensitivity # Required for `Optimization.AutoZygote()` automatic diff
```

Next, we declare our model, the Brusselator oscillator.
```@example pe1
```@example pe_osc_example
brusselator = @reaction_network begin
A, ∅ --> X
1, 2X + Y --> 3X
Expand All @@ -24,7 +24,7 @@ nothing # hide

We simulate our model, and from the simulation generate sampled data points
(to which we add noise). We will use this data to fit the parameters of our model.
```@example pe1
```@example pe_osc_example
u0 = [:X => 1.0, :Y => 1.0]
tspan = (0.0, 30.0)
Expand All @@ -37,7 +37,7 @@ nothing # hide
```

We can plot the real solution, as well as the noisy samples.
```@example pe1
```@example pe_osc_example
using Plots
default(; lw = 3, framestyle = :box, size = (800, 400))
Expand All @@ -49,7 +49,7 @@ Next, we create a function to fit the parameters using the `ADAM` optimizer. For
a given initial estimate of the parameter values, `pinit`, this function will
fit parameter values, `p`, to our data samples. We use `tend` to indicate the
time interval over which we fit the model.
```@example pe1
```@example pe_osc_example
function optimise_p(pinit, tend)
function loss(p, _)
newtimes = filter(<=(tend), sample_times)
Expand All @@ -71,12 +71,12 @@ nothing # hide
```

Next, we will fit a parameter set to the data on the interval `(0, 10)`.
```@example pe1
```@example pe_osc_example
p_estimate = optimise_p([5.0, 5.0], 10.0)
```

We can compare this to the real solution, as well as the sample data
```@example pe1
```@example pe_osc_example
newprob = remake(prob; tspan = (0., 10.), p = p_estimate)
sol_estimate = solve(newprob, Rosenbrock23())
plot(sol_real; color = [:blue :red], label = ["X real" "Y real"], linealpha = 0.2)
Expand All @@ -88,7 +88,7 @@ plot!(sol_estimate; color = [:darkblue :darkred], linestyle = :dash,

Next, we use this parameter estimate as the input to the next iteration of our
fitting process, this time on the interval `(0, 20)`.
```@example pe1
```@example pe_osc_example
p_estimate = optimise_p(p_estimate, 20.)
newprob = remake(prob; tspan = (0., 20.), p = p_estimate)
sol_estimate = solve(newprob, Rosenbrock23())
Expand All @@ -101,20 +101,20 @@ plot!(sol_estimate; color = [:darkblue :darkred], linestyle = :dash,

Finally, we use this estimate as the input to fit a parameter set on the full
time interval of the sampled data.
```@example pe1
```@example pe_osc_example
p_estimate = optimise_p(p_estimate, 30.0)
newprob = remake(prob; tspan = (0., 30.0), p = p_estimate)
sol_estimate = solve(newprob, Rosenbrock23())
plot(sol_real; color = [:blue :red], label = ["X real" "Y real"], linealpha = 0.2)
scatter!(sample_times, sample_vals'; color = [:blue :red],
label = ["Samples of X" "Samples of Y"], alpha = 0.4)
label = ["Samples of X" "Samples of Y"], alpha = 0.4)
plot!(sol_estimate; color = [:darkblue :darkred], linestyle = :dash,
label = ["X estimated" "Y estimated"], xlimit = tspan)
```

The final parameter estimate is then
```@example pe1
```@example pe_osc_example
p_estimate
```
which is close to the actual parameter set of `[1.0, 2.0]`.
Expand All @@ -125,7 +125,7 @@ then extend the interval, is to avoid getting stuck in a local minimum. Here
specifically, we chose our initial interval to be smaller than a full cycle of
the oscillation. If we had chosen to fit a parameter set on the full interval
immediately we would have obtained poor fit and an inaccurate estimate for the parameters.
```@example pe1
```@example pe_osc_example
p_estimate = optimise_p([5.0,5.0], 30.0)
newprob = remake(prob; tspan = (0.0,30.0), p = p_estimate)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ using DiffEqParamEstim, Optimization
ps_dummy = [:kB => 0.0, :kD => 0.0, :kP => 0.0]
oprob = ODEProblem(rn, u0, (0.0, 10.0), ps_dummy)
loss_function = build_loss_objective(oprob, Tsit5(), L2Loss(data_ts, data_vals), Optimization.AutoForwardDiff();
maxiters = 10000, verbose = false, save_idxs = 4)
maxiters = 10000, verbose = false, save_idxs = 4)
nothing # hide
```
To `build_loss_objective` we provide the following arguments:
Expand Down
Loading

0 comments on commit 635e824

Please sign in to comment.