# Generating surrogates

## Introduction

In order to demonstrate how JuliaSim.jl could be used in conjunction with OrdinaryDiffEq.jl, we construct a concise example, namely, the so-called ROBER problem, which consists of a stiff system of 3 non-linear ordinary differential equations (ODEs).

First, we define the ODE in question:

``````using QuasiMonteCarlo, OrdinaryDiffEq, ModelingToolkit, Surrogates, JuliaSim
function rober(du, u, p, t)
y₁, y₂, y₃ = u # initial vector
k₁, k₂, k₃ = p # rate constants
du = -k₁ * y₁ + k₃ * y₂ * y₃
du =  k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃
du =  k₂ * y₂^2
nothing
end

tstop = 1e4
p = [0.04, 3e7, 1e4]
u0 = [1.0, 0.0, 0.0] # initial condition
tspan = (0.0, tstop)

prob = ODEProblem(rober, u0, tspan, p)
sol = solve(prob, Rosenbrock23())``````

This creates an `ODEProblem`, which can be simulated for a set of parameters. Now we need to specify a space of parameters over which we want to train our surrogate:

``param_space = [(0.036, 0.044), (2.7e7, 3.3e7), (0.9e4, 1.1e4)]``

Next, we specify the surrogate algorithm and cast the `ODEProblem` to `DEProblemSimulation`:

``````surralg = LPCTESN(1000)
sim = DEProblemSimulation(prob)``````

The abbreviation `LPCTESN` stands for 'Linear Projection Continuous-Time Echo State Network' (consult this paper for more information). The integer `1000` specifies reservoir size. It should always be provided by the user.

Finally, we call `surrogatize` and `simulate_ensemble` as follows:

``````odesurrogate = JuliaSim.surrogatize(
sim,
param_space,
surralg,
100; # 'n_sample_pts'
ensemble_kwargs = (;),
component_name = :robertson_surrogate,
verbose=true)
odesimresults = JuliaSim.simulate_ensemble(sim, param_space, 100)``````

Note that the number of sets to be sampled from the parameter space (`n_sample_pts`) has to be provided by the user.

As can be inferred from this example, JuliaSim is a very intuitive tool allowing for high composability within the SciML ecosystem (and beyond). It ought to be emphasized that this simple example does not fully reflect JuliaSim's capabilities (as the paragraphs above clearly demonstrate).

### Plotting functions

Numerical results would be incomplete without proper visualization and summaries. JuliaSim provides the necessary plotting functionalities for an easier analysis of the results. In its simplest form, the user may just call:

``plot(odesimresults; ns = 1:2, output_dim = 2) # ns stands for 'samples to be plotted'``

to obtain a visualization of the ROBER problem (see above).

The user may also wish to inspect the training statistics of the CTESN:

``plot(odetrainstats, odesurrogate; ns = 1:2, output_dim = 2, log_scale = true)``

Below is a sample plot: Finally, a comprehensive diagnostic report can be obtained by calling the `weave_diagnostics` function:

``````JuliaSim.weave_diagnostics(
surrogate,
"DCPM Temperature Diagnostic Report",  # Heading of the Report
"output_fmu.pdf";                      # Path to save the Generated Report
doctype = "md2pdf",                    # Any output format from markdown works.
log_scale = false,
include_errors = [:pointwise_error, :curve_distance],
use_absolute_err = false,
visualize_n_samples = 1
)``````

Here is a part of a sample diagnostic report of the ROBER problem:  