From 50b1068a9b1e8022be24239f11718d332de3b999 Mon Sep 17 00:00:00 2001
From: Arthur Carcano <53921575+krtab@users.noreply.github.com>
Date: Fri, 11 Jun 2021 11:30:04 +0200
Subject: [PATCH] doc

---
 Readme.md  | 98 ++++++++++++++++++++++++++++++++++++++++++++++++++++--
 src/lib.rs | 98 +++++++++++++++++++++++++++++++++++++++++++++++++++++-
 2 files changed, 193 insertions(+), 3 deletions(-)

diff --git a/Readme.md b/Readme.md
index 1bf5cf5..89ee107 100644
--- a/Readme.md
+++ b/Readme.md
@@ -1,5 +1,99 @@
-A wrapper around the sundials ODE solver.
+A wrapper around the cvode(S) ODE solver from sundials.
+
+[ [documentation](https://docs.rs/cvode-wrap) ] [ [lib.rs](https://lib.rs/crates/cvode-wrap) ] [ [git repository](https://gitlab.inria.fr/InBio/Public/cvode-rust-wrap) ]
 
 # Examples
 
-Examples computing the behavior of an oscillatory system defined by `x'' = -k * x` are included in the examples/ directory. In the example computing the sensitivities, sensitivities are computed with respect to `x(0)`, `x'(0)` and `k`.
\ No newline at end of file
+## Oscillator
+
+An oscillatory system defined by `x'' = -k * x`.
+
+### Without sensitivities
+
+```rust
+let y0 = [0., 1.];
+//define the right-hand-side
+fn f(_t: Realtype, y: &[Realtype; 2], ydot: &mut [Realtype; 2], k: &Realtype) -> RhsResult {
+    *ydot = [y[1], -y[0] * k];
+    RhsResult::Ok
+}
+//initialize the solver
+let mut solver = SolverNoSensi::new(
+    LinearMultistepMethod::Adams,
+    f,
+    0.,
+    &y0,
+    1e-4,
+    AbsTolerance::scalar(1e-4),
+    1e-2,
+)
+.unwrap();
+//and solve
+let ts: Vec<_> = (1..100).collect();
+println!("0,{},{}", y0[0], y0[1]);
+for &t in &ts {
+    let (_tret, &[x, xdot]) = solver.step(t as _, StepKind::Normal).unwrap();
+    println!("{},{},{}", t, x, xdot);
+}
+```
+
+### With sensitivities
+
+The sensitivities are computed with respect to `x(0)`, `x'(0)` and `k`.
+
+```rust
+let y0 = [0., 1.];
+//define the right-hand-side
+fn f(_t: Realtype, y: &[Realtype; 2], ydot: &mut [Realtype; 2], k: &Realtype) -> RhsResult {
+    *ydot = [y[1], -y[0] * k];
+    RhsResult::Ok
+}
+//define the sensitivity function for the right hand side
+fn fs(
+    _t: Realtype,
+    y: &[Realtype; 2],
+    _ydot: &[Realtype; 2],
+    ys: [&[Realtype; 2]; N_SENSI],
+    ysdot: [&mut [Realtype; 2]; N_SENSI],
+    k: &Realtype,
+) -> RhsResult {
+    // Mind that when indexing sensitivities, the first index
+    // is the parameter index, and the second the state variable
+    // index
+    *ysdot[0] = [ys[0][1], -ys[0][0] * k];
+    *ysdot[1] = [ys[1][1], -ys[1][0] * k];
+    *ysdot[2] = [ys[2][1], -ys[2][0] * k - y[0]];
+    RhsResult::Ok
+}
+
+const N_SENSI: usize = 3;
+
+// the sensitivities in order are d/dy0[0], d/dy0[1] and d/dk
+let ys0 = [[1., 0.], [0., 1.], [0., 0.]];
+
+//initialize the solver
+let mut solver = SolverSensi::new(
+    LinearMultistepMethod::Adams,
+    f,
+    fs,
+    0.,
+    &y0,
+    &ys0,
+    1e-4,
+    AbsTolerance::scalar(1e-4),
+    SensiAbsTolerance::scalar([1e-4; N_SENSI]),
+    1e-2,
+)
+.unwrap();
+//and solve
+let ts: Vec<_> = (1..100).collect();
+println!("0,{},{}", y0[0], y0[1]);
+for &t in &ts {
+    let (_tret, &[x, xdot], [&[dy0_dy00, dy1_dy00], &[dy0_dy01, dy1_dy01], &[dy0_dk, dy1_dk]]) =
+        solver.step(t as _, StepKind::Normal).unwrap();
+    println!(
+        "{},{},{},{},{},{},{},{},{}",
+        t, x, xdot, dy0_dy00, dy1_dy00, dy0_dy01, dy1_dy01, dy0_dk, dy1_dk
+    );
+}
+```
\ No newline at end of file
diff --git a/src/lib.rs b/src/lib.rs
index ffdd138..95d16b9 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -1,7 +1,103 @@
 //! A wrapper around cvode and cvodes from the sundials tool suite.
 //!
 //! Users should be mostly interested in [`SolverSensi`] and [`SolverNoSensi`].
-
+//! # Examples
+//!
+//! ## Oscillator
+//!
+//! An oscillatory system defined by `x'' = -k * x`.
+//!
+//! ### Without sensitivities
+//!
+//! ```rust
+//! use cvode_wrap::*;
+//! let y0 = [0., 1.];
+//! //define the right-hand-side
+//! fn f(_t: Realtype, y: &[Realtype; 2], ydot: &mut [Realtype; 2], k: &Realtype) -> RhsResult {
+//!     *ydot = [y[1], -y[0] * k];
+//!     RhsResult::Ok
+//! }
+//! //initialize the solver
+//! let mut solver = SolverNoSensi::new(
+//!     LinearMultistepMethod::Adams,
+//!     f,
+//!     0.,
+//!     &y0,
+//!     1e-4,
+//!     AbsTolerance::scalar(1e-4),
+//!     1e-2,
+//! )
+//! .unwrap();
+//! //and solve
+//! let ts: Vec<_> = (1..100).collect();
+//! println!("0,{},{}", y0[0], y0[1]);
+//! for &t in &ts {
+//!     let (_tret, &[x, xdot]) = solver.step(t as _, StepKind::Normal).unwrap();
+//!     println!("{},{},{}", t, x, xdot);
+//! }
+//! ```
+//!
+//! ### With sensitivities
+//!
+//! The sensitivities are computed with respect to `x(0)`, `x'(0)` and `k`.
+//!
+//! ```rust
+//! use cvode_wrap::*;
+//! let y0 = [0., 1.];
+//! //define the right-hand-side
+//! fn f(_t: Realtype, y: &[Realtype; 2], ydot: &mut [Realtype; 2], k: &Realtype) -> RhsResult {
+//!     *ydot = [y[1], -y[0] * k];
+//!     RhsResult::Ok
+//! }
+//! //define the sensitivity function for the right hand side
+//! fn fs(
+//!     _t: Realtype,
+//!     y: &[Realtype; 2],
+//!     _ydot: &[Realtype; 2],
+//!     ys: [&[Realtype; 2]; N_SENSI],
+//!     ysdot: [&mut [Realtype; 2]; N_SENSI],
+//!     k: &Realtype,
+//! ) -> RhsResult {
+//!     // Mind that when indexing sensitivities, the first index
+//!     // is the parameter index, and the second the state variable
+//!     // index
+//!     *ysdot[0] = [ys[0][1], -ys[0][0] * k];
+//!     *ysdot[1] = [ys[1][1], -ys[1][0] * k];
+//!     *ysdot[2] = [ys[2][1], -ys[2][0] * k - y[0]];
+//!     RhsResult::Ok
+//! }
+//!
+//! const N_SENSI: usize = 3;
+//!
+//! // the sensitivities in order are d/dy0[0], d/dy0[1] and d/dk
+//! let ys0 = [[1., 0.], [0., 1.], [0., 0.]];
+//!
+//! //initialize the solver
+//! let mut solver = SolverSensi::new(
+//!     LinearMultistepMethod::Adams,
+//!     f,
+//!     fs,
+//!     0.,
+//!     &y0,
+//!     &ys0,
+//!     1e-4,
+//!     AbsTolerance::scalar(1e-4),
+//!     SensiAbsTolerance::scalar([1e-4; N_SENSI]),
+//!     1e-2,
+//! )
+//! .unwrap();
+//! //and solve
+//! let ts: Vec<_> = (1..100).collect();
+//! println!("0,{},{}", y0[0], y0[1]);
+//! for &t in &ts {
+//!     let (_tret, &[x, xdot], [&[dy0_dy00, dy1_dy00], &[dy0_dy01, dy1_dy01], &[dy0_dk, dy1_dk]]) =
+//!         solver.step(t as _, StepKind::Normal).unwrap();
+//!     println!(
+//!         "{},{},{},{},{},{},{},{},{}",
+//!         t, x, xdot, dy0_dy00, dy1_dy00, dy0_dy01, dy1_dy01, dy0_dk, dy1_dk
+//!     );
+//! }
+//! ```
 use std::{ffi::c_void, os::raw::c_int, ptr::NonNull};
 
 use sundials_sys::realtype;