diff --git a/README.md b/README.md
index be562a9..e384c20 100644
--- a/README.md
+++ b/README.md
@@ -1,132 +1,54 @@
-# IFDIFF - A Matlab Toolkit for ODEs with State˗Dependent Switches
+# IFDIFF - A MATLAB Toolkit for ODEs with State˗Dependent Switches
[](https://matlab.mathworks.com/open/github/v1?repo=andreassommer/ifdiff&file=toolbox/doc/GettingStarted.mlx)
-The software package IFDIFF deals with the solution and algorithmic generation of sensitivities in
-ordinary differential equations with implicit (state-dependent) non-differentiabilites ("switches")
-in the right-hand side that is given as Matlab program code with non-differentiable operators such as
-`min`, `max`, `abs`, `sign`, as well as `if`-branching. IFDIFF automatically generates only necessary
-switching functions, outputs them as Matlab code, and detects switching points accurately up to
-machine precision.
+The software package IFDIFF comprises:
+- **Automatic detection and processing** of state-dependent switching events in ODE IVPs
+- **Automatic generation** of only the necessary **switching functions** (exported as MATLAB code)
+- Accurate switching point detection up to machine precision
+- Detection and handling of **Filippov sliding mode**
+- ODEs with state-jumps and model switches
+- **Sensitivity generation** for switched ODEs, ODEs with state-jumps and Filippov ODEs
Sensitivities can be generated w.r.t. the initial values and w.r.t. a given parameter set.
-See the [IFDIFF page](https://andreassommer.github.io/ifdiff/) for a mathematical introduction with example.
+For a mathematical introduction and illustrative examples, see the
+[IFDIFF project page](https://andreassommer.github.io/ifdiff/).
-The file [Readme_Example.m](./toolbox/examples/Readme_Example.m) contains a self-explaining Matlab script similar to the
- contents below.
+A compact, self-explanatory MATLAB example is provided in this file. ([`Readme_Example.m`](./toolbox/examples/Readme_Example.m))
-
-
-
-
-# First Run Prerequisites
-
-Before using IFDIFF, it is mandatory to run the `make_mtreeplus` script once.
-After starting Matlab, change to the IFDIFF directory and type `make_mtreeplus`.
-This scripts generates a modified copy of Matlab's own parser class `mtree`, on which IFDIFF heavily relies.
-
-It is advisable to initialize the paths needed for IFDIFF by invoking
-`initIFDIFF();` once on every new matlab session.
+# Installation
-# Usage
-
-Make sure that you've followed the [First Run Prerequisites](#first-run-prerequisites).
-
-The following commands exemplarily show the usage on the canonical example.
-
-## 1. Preprocessing of Right Hand Side
-
-The preprocessing analyses and transforms the right hand side function `canonicalExampleRHS.m`.
-Preprocessing generates the `datahandle`, the central structure for switching detection and handling.
-We set ODE solver and its options as usual. If not set, default values are used.
-
- ```matlab
- initIFDIFF(); % Initialise the paths for ifdiff -- needed only once per Matlab session
- integrator = @ode45;
- odeoptions = odeset('AbsTol', 1e-14, 'RelTol', 1e-12);
- datahandle = prepareDatahandleForIntegration('canonicalExampleRHS', 'integrator', integrator, 'options', odeoptions);
- ```
+There are two ways to install IFDIFF.
+## Installation Method (Recommended)
-## 2. Integration (Forward Solution)
+1. Download the current release file `IFDIFF_Toolbox.mltbx`.
+2. Open MATLAB and navigate to the directory containing the file.
+3. Right-click the file and select **Install**.
+4. For a detailed introduction, open `GettingStarted.mlx`, which is included with the toolbox.
-Define initial values, parameter values, and integration horizon, and call `solveODE` to start the integration.
-`solveODE` returns a Matlab `sol` structure, that can be evalated using `deval`.
-It is an augmented version of the solution structures returned by Matlab's very own integrator suite
-(see https://de.mathworks.com/help/matlab/ref/deval.html#bu7iw_j-sol)
+Note: You can also open the Getting Started guide directly in
+[MATLAB Online](https://matlab.mathworks.com/open/github/v1?repo=andreassommer/ifdiff&file=toolbox/doc/GettingStarted.mlx) and then continue with the [First Run Prerequisites](#first-run-prerequisites).
- ```matlab
- tspan = [0 20];
- initialvalues = [1;0];
- parameters = 5.437;
- sol = solveODE(datahandle, tspan, initialvalues, parameters);
- T = 0:0.1:20;
- X = deval(sol,T);
- plot(T,X)
-```
+## Installation Method (Alternative)
+1. Navigate to a location of your choice on your PC (e.g. Desktop).
+2. Clone the repository `git clone https://github.com/andreassommer/ifdiff.git`
+3. Open MATLAB and navigate to the cloned `ifdiff` directory.
-## 3. Sensitivity Generation
-
-The sensitivity generation currently supports 1st order forward sensitivities w.r.t. initial state and parameters.
-It is possible to generate sensitivities using external numerical differentiation (flags `END_full`, `END_piecewise`)
-or using the variational differential equations (flag `VDE`).
-
-Usually, method `VDE` delivers vest results in terms of accuracy, as it calculates error-controlled sensitivities.
-It uses the variational differential equations on each interval and performes updates at the switching points to
-ensure accurate sensitivities at the precalculated forward solution. However, the occuring augmented differential
-system might be large size and thus slow to compute.
-Required state derivatives of the right hand side are approximated using automated finite differencing.
-
-The method `END_piecewise` computes the interval sensitivities using external numerical differentiation (finite differencing)
-and connects these using the same updates as used in the VDE method. This requires multiple evaluations of interval solutions,
-but might be faster on larger systems.
-
-When using `END_full`, external numerical differention is used on multiple full horizon solutions. The individual
-trajectories are calculated with switching point detection. Thus, no updates between switches are required.
-In general, this approach is less accurate and slower than `END_piecewise`, but might be the a good choice
-for highly instable ODEs.
-
-1. Choose step sizes for finite differencing (also used in method `VDE` for generating state derivatives of the rhs).
-
- ```matlab
- dim_y = size(sol.y,1); // state dimension
- dim_p = length(parameters); // number of parameters
- FDstep = generateFDstep(dim_y,dim_p);
- ```
-
- The `generateDFstep` function accepts several options influencing e.g. step length.
- See the documentation or the file for more information
-
-
-2. Build the sensitivity function. In this example, the `END_piecewise` method is chosen.
-
- ```matlab
- sensitivity_function = generateSensitivityFunction(datahandle, sol, FDstep, 'method', 'END_piecewise');
- ```
+
- The following table lists several name-value pairs that can be used to configure `generateSensitivityFunction`
+# Usage
- | Parameters | Possible values | Defaults |
- | ----------------------- | ----------------------------------------------------------------------------------- | ------------------------------------- |
- | calcGy | true/false - flag indicating to calculate state sensitivities | true |
- | calcGp | true/false - flag indicating to calculate parameter sensitivities | true |
- | Gmatrices_intermediate | true/false - flag indicating to store update matrics | false |
- | save_intermediates | true/false - flag indicating to store intermediate calculations | true |
- | integrator | Function handle for ODE solver in Matlab (e.g. ode45) | Integrator used by ifdiff |
- | integrator_options | Options struct generated for ODE solver | Integrator options used by ifdiff |
- | method | String with VDE/END_piecewise/END_full | VDE |
- | directions_y | Matrix containing directions for directional derivatives w.r.t initial values. | Identity matrix with dimension n_y |
- | directions_p | Matrix containing directions for directional derivatives w.r.t parameters. | Identity matrix with dimension n_p |
+Once every new MATLAB session run
+ `initIFDIFF();`
-3. Evaluate the sensitivity function at specific times.
+to initialize the required paths. Then IFDIFF is ready to use.
- ```matlab
- t = 0:0.1:20;
- sensitivities = sensitivity_function(t);
- ```
+The first time the command is executed, a copy of MATLAB's internal parser class `mtree` is created on which IFDIFF heavily relies.
+For a step-by-step usage example and background explanation, see the [IFDIFF project page](https://andreassommer.github.io/ifdiff/).
diff --git a/docs/index.md b/docs/index.md
index 56de951..13da4fb 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -1,38 +1,33 @@
-# IFDIFF - A Matlab Toolkit for ODEs with State˗Dependent Switches
+# IFDIFF - A MATLAB Toolkit for ODEs with State˗Dependent Switches
-The software package IFDIFF deals with the solution and algorithmic generation of sensitivities
-in ordinary differential equations with implicit (state-dependent) non-differentiabilites ("switches")
-in the right-hand side that is given as Matlab program code with non-differentiable operators such as
-`min`, `max`, `abs`, `sign`, as well as `if`-branching. IFDIFF automatically generates only necessary
-switching functions, outputs them as Matlab code, and detects switching points accurately up to machine
-precision.
-IFDIFF handles multidimensional state and parameter vectors and can be transparently used in existing code,
-as it generates `sol` solution structures compatible to the Matlab ode solvers. No modifications to the
-right hand side or manual processing are necessary. A single preparation call is sufficient.
+The software package **IFDIFF** provides tools for the accurate numerical solution and sensitivity analysis of ordinary differential equations (ODEs) with state-dependent (implicit) non-differentiabilities ("switches") in the right-hand side as well as Filippov Sliding mode and ODEs with state-jumps.
-In the PErFDiff project, the existing implementation is to be further developed and extended for Filippov systems.
+The right-hand side is supplied as standard MATLAB program code and may contain non-differentiable operators such as `min`, `max`, `abs`, `sign` or `if`-branching.
+IFDIFF automatically detects and processes these switches, generates only the necessary switching functions (exported as MATLAB code), and determines switching points accurately up to machine precision.
-
-
+IFDIFF handles multidimensional state and parameter vectors and produces solution structures compatible with MATLAB's built-in ODE solvers, allowing transparent use of functions such as `deval` without modification of existing code.
+A single preparation call is sufficient to enable switching-aware integration.
+
+
# Canonical Example
Consider the following "canonical example" for a switched ODE system:
-$$
+$${
\dot x = f(t,x,p) = \binom{f_1(t,x,p)}{f_2(t,x,p)}
-$$
+}$$
with
-$$
+$${
f_1(t,x,p) = 0.01 \cdot t^2 + x_2^3
\qquad
f_2(t,x,p) = \begin{cases} 0 ~~~if~~ x_1 < p \\ 5 ~~~if~~ p \leq x_1 < p+0.5 \\ 0 ~~~if~~ x_1 \geq p+0.5 \end{cases}
-$$
+}$$
-and with initial value $ x(0) = (1,0)^T $, parameter $ p = 5.437 $, over time span $ t \in [0,20] $.
+and with initial value ${ x(0) = (1,0)^T }$, parameter ${ p = 5.437 }$, over time span ${ t \in [0,20] }$.
This switched ODE system translates straightforward into the following matlab program:
@@ -58,46 +53,48 @@ and a second "kink" once the first component reaches `p+0.5`.
Let's see what happens if we do not consider appropriate switching handling.
-
-
+
# Naive Integration with ode45 fails unnoticed
-We initialize the variables and start the integration using Matlab's default integrator `ode45` (explicit Runge-Kutta 4(5)-solver),
+We initialize the variables and start the integration using MATLAB's default integrator `ode45` (explicit Runge-Kutta 4(5)-solver),
without caring for the non-differentiable `if` statements.
+
```matlab
- tspan = [0 20]; % time horizon
- x0 = [1;0]; % initial values
- p = 5.437; % parameter values
- solX = ode45(@(t,x) canonicalExampleRHS(t,x,p), tspan, x0)
- T = 0:0.1:20; X = deval(solX2,T); plot(T,X); legend('x_1','x_2');
+tspan = [0 20]; % time horizon
+x0 = [1;0]; % initial values
+p = 5.437; % parameter values
+solX = ode45(@(t,x) canonicalExampleRHS(t,x,p), tspan, x0)
+T = 0:0.1:20; X = deval(solX2,T); plot(T,X); legend('x_1','x_2');
```
Contrary to wide-spread beliefs, tightening the integration tolerances is not a remedy!
```matlab
- odeopts = odeset('AbsTol', 1e-20, 'RelTol', 1e-14);
- solX2 = ode45(@(t,x) canonicalExampleRHS(t,x,p), tspan, x0)
+odeopts = odeset('AbsTol', 1e-20, 'RelTol', 1e-14);
+solX2 = ode45(@(t,x) canonicalExampleRHS(t,x,p), tspan, x0)
```
Without any warning or error, the __naive approach with tight tolerances leads to the same (wrong) result__. Again, no kinks visible.
-
-
+
# Reliable integration with IFDIFF
With the switching point detection in IFDIFF, after a single call to a preparation routine,
integration is just as simple as before:
```matlab
- initIFDIFF(); % initialise IFDIFF (only once)
- tspan = [0 20]; x0 = [1;0]; p = 5.437; % set time horizon, initial value, parameter
- integrator = @ode45; % choose integrator
- odeoptions = odeset('AbsTol', 1e-5, 'RelTol', 1e-3); % set integrator options, here: low accuracy
- datahandle = prepareDatahandleForIntegration('canonicalExampleRHS', 'integrator', integrator, 'options', odeoptions);
- sol = solveODE(datahandle, tspan, x0, p);
- T = 0:0.1:20; X = deval(solX2,T); plot(T,X); legend('x_1','x_2');
+initIFDIFF(); % initialise IFDIFF (only once)
+tspan = [0 20]; x0 = [1;0]; p = 5.437; % set time horizon, initial value, parameter
+integrator = @ode45; % choose integrator
+odeoptions = odeset('AbsTol', 1e-5, 'RelTol', 1e-3); % set integrator options, here: low accuracy
+
+datahandle = prepareDatahandleForIntegration('canonicalExampleRHS', ...
+ 'integrator', integrator, ...
+ 'options', odeoptions);
+sol = solveODE(datahandle, tspan, x0, p);
+T = 0:0.1:20; X = deval(solX2,T); plot(T,X); legend('x_1','x_2');
```
and __IFDIFF delivers the correct result__:
@@ -107,36 +104,73 @@ are exact up to integration tolerance. Using, e.g., absolute and relative tolera
first switching times is computed an error less than $10^{-14}$, and the second with an error less than $10^{-11}$.
Note that the `sol` structure returned by `solveODE` is an augmented version of the solution structures returned
-by Matlab's very own integrators (see [Matlab documentation](https://de.mathworks.com/help/matlab/ref/deval.html#bu7iw_j-sol)),
+by MATLAB's very own integrators (see [MATLAB documentation](https://de.mathworks.com/help/matlab/ref/deval.html#bu7iw_j-sol)),
and can thus be evaluated using `deval`. That means, IFDIFF can be used transparently within existing code!
-
+
## Remarks
-While in many (not all) cases, forcing a small integration step size (e.g. using `odeset('MaxStep',1e-2)` in the above example),
-may lead to approximative solutions, this is not a remedy at all.
+While in some cases forcing a small integration step size (for example by using `odeset('MaxStep',1e-2)` in the example above) may produce a visually plausible solution, this approach is not a reliable remedy.
+
+Firstly, restricting the maximum step size contradicts the fundamental idea of adaptive step-size integrators. It globally limits the step size, even in regions where no switching occurs.
+In the example above, choosing a parameter value of `p = 100` would eliminate all switches within the integration horizon, yet the step size would still be unnecessarily constrained.
+
+Secondly, when non-differentiabilities are ignored, there is no principled way to choose an appropriate maximum step size. The integrator provides neither warnings nor error estimates indicating that switching events have been missed.
+
+Thirdly, and as demonstrated in the following section, this approach completely fails to produce meaningful sensitivities.
+
+
+
+# Sensitivity Generation
+
+Correct state and parameter sensitivities are crucial for derivative-based optimization and control.
+While restricting the maximum integrator step size (e.g. to `0.1`) may allow the naive approach to reproduce the forward trajectory, the resulting sensitivities are completely unreliable and generally unusable.
+
+IFDIFF currently supports first-order forward sensitivities w.r.t. both the initial state and model parameters. Sensitivities can be computed using either external numerical differentiation or variational differential equations:
-First, restricting the maximum step size is completely against the idea of adaptive step size integrators, as it reduces the maximum
-step size everywhere - even if there is no switch at all (consider, in the above example, a parameter value of `p=100`; then there is
-no switch in the specified integration horizon).
+- **Variational Differential Equations** (`VDE`)
+This method typically provides the highest accuracy, as it computes error-controlled sensitivities alongside the forward solution. The variational equations are integrated on each interval between switching events, with appropriate updates applied at the detected switching points to ensure consistency with the forward trajectory.
+The required derivatives of the right-hand side are approximated using automated finite differencing.
+Due to the augmented system size, this method may become computationally expensive for large-scale problems.
-Second, using the naive approach (i.e. ignoring non-differentiabilities completely), we don't know how to choose ths maximum step size
-appropriately, as we don't get any error messages or warnings at all.
+- **External Numerical Differentiation, Piecewise** (`END_piecewise`)
+In this approach, sensitivities are computed via finite differencing on each interval between switching events and then connected using the same update rules as in the VDE method.
+Although this requires multiple forward integrations per interval, it can be more efficient than VDE for larger systems while still providing accurate sensitivities.
-Third, and this we will show in the next section, sensitivites cannot be computed at all by this approach.
+- **External Numerical Differentiation, Full Horizon** (`END_full`) **NOT Recommended**
+Here, finite differencing is applied to multiple full-horizon solutions, each computed with switching point detection. Since sensitivities are not propagated across switching points, no intermediate updates are required.
+This method is generally not accurate and we do not recommend it, but may be faster for a certain class of problems.
+
+
-
-
+## Example sensitivity generation for the Canonical Example
-# Computing Sensitivities
+1. Choose step sizes for finite differencing (also used in method `VDE` for generating state derivatives of the rhs).
-Having correct state and parameters sensitivities, i.e. the derivative of the trajectory w.r.t. changes in the initial values or parameters,
-is cruicial for derivative-based optimization.
-By restricting the maximum integrator step size to $0.1$, also the naive approach is capable of reproducing the forward trajectory.
-However, the sensitivities generated by the naive approach are completely useless.
+ ```matlab
+ dim_y = size(sol.y,1);
+ dim_p = length(parameters);
+ FDstep = generateFDstep(dim_y,dim_p);
+ ```
-The following picture shows the trajectories of the state sensitivities, $G_{y,ij}(t,t_0) := \frac{d y_i}{d x_{0,j}}(t)$, i.e.
+ The `generateDFstep` function accepts several options influencing e.g. step length.
+ See the documentation below for more information.
+
+
+2. Build the sensitivity function. In this example, the `END_piecewise` method is chosen.
+
+ ```matlab
+ sensitivity_function = generateSensitivityFunction(datahandle, sol, FDstep, 'method', 'END_piecewise');
+ ```
+3. Evaluate the sensitivity function at specific times.
+
+ ```matlab
+ t = 0:0.1:20;
+ sensitivities = sensitivity_function(t);
+ ``
+
+The following picture shows the trajectories of the state sensitivities for the example aboves, $G_{y,ij}(t,t_0) := \frac{d y_i}{d y_{0,j}}(t)$, i.e.
$G_{y,ij}$ denotes the sensitivity of the $i$-th solution component w.r.t. the $j$-th component of the initial value.
| Comparison: Sensitivities w.r.t. initial state | |
@@ -144,30 +178,34 @@ $G_{y,ij}$ denotes the sensitivity of the $i$-th solution component w.r.t. the $
|  |  |
| __Naive approach__: Sensitivities are useless | __IFDIFF__: Correct and accurate sensitivities |
-
-
While IFDIFF produces correct and accurate sensitivities, the sensitivities using the naive approach, again without any warning or hint,
are not only inaccurate (notice the scales!), but plain wrong.
-Details of how to compute sensitivities with IFDIFF are given in the project's [README.md](https://github.com/andreassommer/ifdiff/blob/public/README.md).
+
+
+The following table lists several name-value pairs that can be used to configure `generateSensitivityFunction`
+
+ | Parameters | Possible values | Defaults |
+ | ----------------------- | ----------------------------------------------------------------------------------- | ------------------------------------- |
+ | calcGy | true/false - flag indicating to calculate state sensitivities | true |
+ | calcGp | true/false - flag indicating to calculate parameter sensitivities | true |
+ | Gmatrices_intermediate | true/false - flag indicating to store update matrics | false |
+ | save_intermediates | true/false - flag indicating to store intermediate calculations | true |
+ | integrator | Function handle for ODE solver in MATLAB (e.g. ode45) | Integrator used by ifdiff |
+ | integrator_options | Options struct generated for ODE solver | Integrator options used by ifdiff |
+ | method | String with VDE/END_piecewise/END_full | VDE |
+ | directions_y | Matrix containing directions for directional derivatives w.r.t initial values. | Identity matrix with dimension n_y |
+ | directions_p | Matrix containing directions for directional derivatives w.r.t parameters. | Identity matrix with dimension n_p |
+
+
-
-
# IFDIFF is simple to use!
-Correct treatment of switched systems requires elaborate formulation of switching functions and tailored integrators,
-placing high mathematical demands on modelers. Even small model changes often imply considerable reformulation effort.
-Furthermore: $n$ switches generate up to $2^n$ possible program flows and switching functions, rendering a-priori
-formulations not feasible already in medium-sized models.
+Correct treatment of switched systems requires elaborate formulation of switching functions and tailored integrators, placing high mathematical demands on modelers.
+Even small model changes often imply considerable reformulation effort. Furthermore: ${n}$ switches generate up to ${ 2^n}$ possible program flows and switching functions, rendering a-priori formulations not feasible already in medium-sized models.
-IFDIFF programmatically handles switching events, auto-generating only required switching functions.
-It determines switching times up to machine precision, and ensures accurate simulation and sensitivity results.
-Transparently extending the Matlab integrators (ode45, ode15s, etc.), IFDIFF is applicable to existing code with
-state- and parameter-dependent conditionals, thus enabling fast prototyping and relieving modelers of
-mathematical-technical effort.
+IFDIFF programmatically handles switching events, auto-generating only required switching functions. It determines switching times up to machine precision, and ensures accurate simulation and sensitivity results. Transparently extending the Matlab integrators (ode45, ode15s, etc.), IFDIFF is applicable to existing code with state- and parameter-dependent conditionals, thus enabling fast prototyping and relieving modelers of mathematical-technical effort.
Calculation times of the naive approach are a little lower than the ones of IFDIFF.
-But it can generate arbitrarily wrong results without any notice (see the example above!).
-IFDIFF provides correct integration results, correct first-order sensitivities and information
-about the switching structure of your model.
+But it can generate arbitrarily wrong results without any notice (see the example above!). IFDIFF provides correct integration results, correct first-order sensitivities and information about the switching structure of your model.
\ No newline at end of file