Merge pull request #1160 from ParamThakkar123/ParameterEstimation

Updating ParameterEstimation Benchmarks

Author: ChrisRackauckas

benchmarks/ParameterEstimation/FitzHughNagumoParameterEstimation.jmd

@@ -8,6 +8,7 @@ author: Vaibhav Dixit, Chris Rackauckas
 ```julia
 using ParameterizedFunctions, OrdinaryDiffEq, DiffEqParamEstim, Optimization
 using OptimizationBBO, OptimizationNLopt, ForwardDiff, Plots, BenchmarkTools 
+using ModelingToolkit: t_nounits as t, D_nounits as D
 gr(fmt=:png)
 ```
 
@@ -19,10 +20,24 @@ glo_init = [2.5,2.5,2.5,2.5]
 ```
 
 ```julia
-fitz = @ode_def FitzhughNagumo begin
-  dv = v - v^3/3 -w + l
-  dw = τinv*(v +  a - b*w)
-end a b τinv l
+@mtkmodel FitzHughNagumo begin
+    @parameters begin
+        a = 0.7      # Parameter for excitability
+        b = 0.8      # Recovery rate parameter
+        τinv = 0.08  # Inverse of the time constant
+        l = 0.5      # External stimulus
+    end
+    @variables begin
+        v(t) = 1.0  # Membrane potential with initial condition
+        w(t) = 1.0  # Recovery variable with initial condition
+    end
+    @equations begin
+        D(v) ~ v - v^3 / 3 - w + l
+        D(w) ~ τinv * (v + a - b * w)
+    end
+end
+
+@mtkbuild fitz = FitzHughNagumo()
 ```
 
 ```julia
@@ -38,7 +53,7 @@ prob_short = ODEProblem(fitz, r0, tspan2,p)
 dt = 30.0/3000
 tf = 30.0
 tinterval = 0:dt:tf
-t  = collect(tinterval)
+time_points  = collect(tinterval)
 ```
 
 ```julia
@@ -54,7 +69,7 @@ t_short = collect(tinterval_short)
 #Generate Data
 data_sol_short = solve(prob_short,Vern9(),saveat=t_short,reltol=1e-9,abstol=1e-9)
 data_short = convert(Array, data_sol_short) # This operation produces column major dataset obs as columns, equations as rows
-data_sol = solve(prob,Vern9(),saveat=t,reltol=1e-9,abstol=1e-9)
+data_sol = solve(prob,Vern9(),saveat=time_points,reltol=1e-9,abstol=1e-9)
 data = convert(Array, data_sol)
 ```
 
@@ -177,7 +192,7 @@ Vern9 solver with reltol=1e-9 and abstol=1e-9 is used and the dataset is increas
 
 
 ```julia
-obj = build_loss_objective(prob,Vern9(),L2Loss(t,data),Optimization.AutoForwardDiff(),tstops=t,reltol=1e-9,abstol=1e-9)
+obj = build_loss_objective(prob,Vern9(),L2Loss(time_points,data),Optimization.AutoForwardDiff(),tstops=time_points,reltol=1e-9,abstol=1e-9)
 optprob = OptimizationProblem(obj, glo_init, lb = first.(glo_bounds), ub = last.(glo_bounds))
 @btime res1 = solve(optprob, BBO_adaptive_de_rand_1_bin(), maxiters = 4e3)
 ```

benchmarks/ParameterEstimation/LorenzParameterEstimation.jmd

@@ -5,13 +5,12 @@ author: finmod, Chris Rackauckas, Vaibhav Dixit
 
 # Estimate the parameters of the Lorenz system from the dataset
 
-
-
 Note: If data is generated with a fixed time step method and then is tested against with the same time step, there is a biased introduced since it's no longer about hitting the true solution, rather it's just about retrieving the same values that the ODE was first generated by! Thus this version uses adaptive timestepping for all portions so that way tests are against the true solution.
 
 ```julia
 using ParameterizedFunctions, OrdinaryDiffEq, DiffEqParamEstim, Optimization
 using OptimizationBBO, OptimizationNLopt, Plots, ForwardDiff, BenchmarkTools
+using ModelingToolkit: t_nounits as t, D_nounits as D
 gr(fmt=:png)
 ```
 
@@ -25,11 +24,25 @@ LocIniPar = [9.0, 20.0, 2.0] # for local optimization
 ```
 
 ```julia
-g1 = @ode_def LorenzExample begin
-  dx = σ*(y-x)
-  dy = x*(ρ-z) - y
-  dz = x*y - β*z
-end σ ρ β
+@mtkmodel LorenzExample begin
+  @parameters begin
+      σ = 10.0  # Parameter: Prandtl number
+      ρ = 28.0  # Parameter: Rayleigh number
+      β = 8/3   # Parameter: Geometric factor
+  end
+  @variables begin
+      x(t) = 1.0  # Initial condition for x
+      y(t) = 1.0  # Initial condition for y
+      z(t) = 1.0  # Initial condition for z
+  end
+  @equations begin
+      D(x) ~ σ * (y - x)
+      D(y) ~ x * (ρ - z) - y
+      D(z) ~ x * y - β * z
+  end
+end
+  
+@mtkbuild g1 = LorenzExample()
 p = [10.0,28.0,2.66] # Parameters used to construct the dataset
 r0 = [1.0; 0.0; 0.0]                #[-11.8,-5.1,37.5] PODES Initial values of the system in space # [0.1, 0.0, 0.0]
 tspan = (0.0, 30.0)                 # PODES sample of 3000 observations over the (0,30) timespan
@@ -42,7 +55,7 @@ prob_short = ODEProblem(g1, r0, tspan2,p)
 dt = 30.0/3000
 tf = 30.0
 tinterval = 0:dt:tf
-t  = collect(tinterval)
+time_points  = collect(tinterval)
 ```
 
 ```julia
@@ -58,7 +71,7 @@ t_short = collect(tinterval_short)
 # Generate Data
 data_sol_short = solve(prob_short,Vern9(),saveat=t_short,reltol=1e-9,abstol=1e-9)
 data_short = convert(Array, data_sol_short) # This operation produces column major dataset obs as columns, equations as rows
-data_sol = solve(prob,Vern9(),saveat=t,reltol=1e-9,abstol=1e-9)
+data_sol = solve(prob,Vern9(),saveat=time_points,reltol=1e-9,abstol=1e-9)
 data = convert(Array, data_sol)
 ```
 
@@ -204,7 +217,7 @@ Vern9 solver with reltol=1e-9 and abstol=1e-9 has been established to be accurat
 
 ```julia
 # BB with Vern9 converges very slowly. The final values are within the NarrowBounds.
-obj = build_loss_objective(prob,Vern9(),L2Loss(t,data),tstops=t,reltol=1e-9,abstol=1e-9)
+obj = build_loss_objective(prob,Vern9(),L2Loss(time_points,data),tstops=time_points,reltol=1e-9,abstol=1e-9)
 optprob = OptimizationProblem(obj, GloIniPar, lb = first.(LooserBounds), ub = last.(LooserBounds))
 @btime res1 = solve(optprob, BBO_adaptive_de_rand_1_bin(); maxiters = 4e3) # Default adaptive_de_rand_1_bin_radiuslimited 33 sec [10.2183, 24.6711, 2.28969]
 #@btime res1 = bboptimize(obj;SearchRange = LooserBounds, Method = :adaptive_de_rand_1_bin, MaxSteps = 4e3) # Method 32 sec [13.2222, 25.8589, 2.56176]

benchmarks/ParameterEstimation/LotkaVolterraParameterEstimation.jmd

@@ -8,21 +8,37 @@ author: Vaibhav Dixit, Chris Rackauckas
 ```julia
 using ParameterizedFunctions, OrdinaryDiffEq, DiffEqParamEstim, Optimization, ForwardDiff
 using OptimizationBBO, OptimizationNLopt, Plots, RecursiveArrayTools, BenchmarkTools
+using ModelingToolkit: t_nounits as t, D_nounits as D
 gr(fmt=:png)
 ```
 
 ```julia
 loc_bounds = Tuple{Float64, Float64}[(0, 5), (0, 5), (0, 5), (0, 5)]
 glo_bounds = Tuple{Float64, Float64}[(0, 10), (0, 10), (0, 10), (0, 10)]
 loc_init = [1,0.5,3.5,1.5]
-glo_init = [5,5,5,5]
+glo_init = [5.0,5.0,5.0,5.0]
 ```
 
 ```julia
-f = @ode_def LotkaVolterraTest begin
-    dx = a*x - b*x*y
-    dy = -c*y + d*x*y
-end a b c d
+@mtkmodel LotkaVolterraTest begin
+    @parameters begin
+        a = 1.5  # Growth rate of prey
+        b = 1.0  # Predation rate
+        c = 3.0  # Death rate of predators
+        d = 1.0  # Reproduction rate of predators
+        e = 0.0  # Additional parameter (if needed)
+    end
+    @variables begin
+        x(t) = 1.0  # Population of prey with initial condition
+        y(t) = 1.0  # Population of predators with initial condition
+    end
+    @equations begin
+        D(x) ~ a * x - b * x * y
+        D(y) ~ -c * y + d * x * y
+    end
+end
+
+@mtkbuild f = LotkaVolterraTest()
 ```
 
 ```julia
@@ -39,7 +55,7 @@ prob_short = ODEProblem(f, u0, tspan2,p)
 dt = 30.0/3000
 tf = 30.0
 tinterval = 0:dt:tf
-t  = collect(tinterval)
+time_points  = collect(tinterval)
 ```
 
 ```julia
@@ -55,7 +71,7 @@ t_short = collect(tinterval_short)
 #Generate Data
 data_sol_short = solve(prob_short,Tsit5(),saveat=t_short,reltol=1e-9,abstol=1e-9)
 data_short = convert(Array, data_sol_short)
-data_sol = solve(prob,Tsit5(),saveat=t,reltol=1e-9,abstol=1e-9)
+data_sol = solve(prob,Tsit5(),saveat=time_points,reltol=1e-9,abstol=1e-9)
 data = convert(Array, data_sol)
 ```
 
@@ -181,8 +197,9 @@ opt = Opt(:LD_TNEWTON_PRECOND_RESTART, 4)
 Vern9 solver with reltol=1e-9 and abstol=1e-9 is used and the dataset is increased to 3000 observations per variable with the same integration time step of 0.01.
 
 ```julia
-obj = build_loss_objective(prob,Vern9(),L2Loss(t,data),tstops=t,reltol=1e-9,abstol=1e-9)
-optprob = OptimizationProblem(obj, glo_init, lb = first.(glo_bounds), ub = last.(glo_bounds))
+t_concrete = collect(0.0:dt:tf)
+obj = build_loss_objective(prob,Vern9(),L2Loss(t_concrete,data),tstops=t_concrete,reltol=1e-9,abstol=1e-9)
+optprob = OptimizationProblem(obj, glo_init, lb =first.(glo_bounds), ub = last.(glo_bounds))
 ```
 
 ```julia