From 66f0d5fbf8d27c2c86958594dfc8421625a61ec0 Mon Sep 17 00:00:00 2001 From: Arthur Carcano Date: Thu, 10 Jun 2021 19:34:57 +0200 Subject: [PATCH] Clean up and add doc --- Cargo.toml | 22 +++-- Readme.md | 38 +------- cvode-wrap/Cargo.toml | 12 --- example/Cargo.toml | 10 --- example/src/main.rs | 80 ----------------- examples/oscillator_no_sensi.rs | 28 ++++++ examples/oscillator_sensi.rs | 58 +++++++++++++ .../plot.py => examples/plot_oscillator.py | 0 {cvode-wrap/src => src}/cvode.rs | 62 +++++-------- {cvode-wrap/src => src}/cvode_sens.rs | 86 ++++++------------- {cvode-wrap/src => src}/lib.rs | 63 +++++++++++++- {cvode-wrap/src => src}/nvector.rs | 0 12 files changed, 211 insertions(+), 248 deletions(-) delete mode 100644 cvode-wrap/Cargo.toml delete mode 100644 example/Cargo.toml delete mode 100644 example/src/main.rs create mode 100644 examples/oscillator_no_sensi.rs create mode 100644 examples/oscillator_sensi.rs rename example/plot.py => examples/plot_oscillator.py (100%) rename {cvode-wrap/src => src}/cvode.rs (82%) rename {cvode-wrap/src => src}/cvode_sens.rs (85%) rename {cvode-wrap/src => src}/lib.rs (65%) rename {cvode-wrap/src => src}/nvector.rs (100%) diff --git a/Cargo.toml b/Cargo.toml index 50aeb65..a721599 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,5 +1,17 @@ -[workspace] -members = [ - "cvode-wrap", - "example" -] \ No newline at end of file +[package] +name = "cvode-wrap" +version = "0.1.0" +authors = ["Arthur Carcano "] +edition = "2018" +license = "BSD-3" +description="A wrapper around cvode and cvodeS from sundials, allowing to solve ordinary differential equations (ODEs) with or without their sensitivities." +repository="https://gitlab.inria.fr/InBio/Public/cvode-rust-wrap/" +readme="Readme.md" +keywords=["sundials","cvode","cvodes","ode","sensitivities"] +categories=["science","simulation","api-bindings"] + +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html + +[dependencies] +sundials-sys = {version="0.2.1", default-features=false, features=["cvodes"]} +array-init = "2.0" \ No newline at end of file diff --git a/Readme.md b/Readme.md index 9b67df0..1bf5cf5 100644 --- a/Readme.md +++ b/Readme.md @@ -1,39 +1,5 @@ A wrapper around the sundials ODE solver. -# Example +# Examples -An oscillatory 2-D system. - -```rust -use cvode_wrap::*; - -let y0 = [0., 1.]; -// define the right-hand-side as a rust function of type RhsF -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 = cvode::Solver::new( - LinearMultistepMethod::Adams, - wrapped_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); -} -``` \ No newline at end of file +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 diff --git a/cvode-wrap/Cargo.toml b/cvode-wrap/Cargo.toml deleted file mode 100644 index d4aba2e..0000000 --- a/cvode-wrap/Cargo.toml +++ /dev/null @@ -1,12 +0,0 @@ -[package] -name = "cvode-wrap" -version = "0.1.0" -authors = ["Arthur Carcano "] -edition = "2018" - -# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html - -[dependencies] -#cvode-5-sys = {path = "../cvode-5-sys"} -sundials-sys = {path = "../../sundials-sys", default-features=false, features=["cvodes"]} -array-init = "2.0" \ No newline at end of file diff --git a/example/Cargo.toml b/example/Cargo.toml deleted file mode 100644 index ff70778..0000000 --- a/example/Cargo.toml +++ /dev/null @@ -1,10 +0,0 @@ -[package] -name = "test-solver" -version = "0.1.0" -authors = ["Arthur Carcano "] -edition = "2018" - -# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html - -[dependencies] -cvode-wrap = {path = "../cvode-wrap"} \ No newline at end of file diff --git a/example/src/main.rs b/example/src/main.rs deleted file mode 100644 index bc9099b..0000000 --- a/example/src/main.rs +++ /dev/null @@ -1,80 +0,0 @@ -use std::env::args; - -use cvode_wrap::*; - -fn main() { - 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 - } - // If there is any command line argument compute the sensitivities, else don't. - if false && args().nth(1).is_none() { - //initialize the solver - let mut solver = cvode::Solver::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); - } - } else { - 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.]]; - - fn fs( - _t: Realtype, - y: &[Realtype; 2], - _ydot: &[Realtype; 2], - ys: [&[Realtype; 2]; N_SENSI], - ysdot: [&mut [Realtype; 2]; N_SENSI], - k: &Realtype, - ) -> RhsResult { - *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 - } - - //initialize the solver - let mut solver = cvode_sens::Solver::new( - LinearMultistepMethod::Adams, - f, - fs, - 0., - &y0, - &ys0, - 1e-4, - AbsTolerance::scalar(1e-4), - cvode_sens::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 - ); - } - } -} diff --git a/examples/oscillator_no_sensi.rs b/examples/oscillator_no_sensi.rs new file mode 100644 index 0000000..ec08bd9 --- /dev/null +++ b/examples/oscillator_no_sensi.rs @@ -0,0 +1,28 @@ +use cvode_wrap::*; + +fn main() { + 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); + } +} diff --git a/examples/oscillator_sensi.rs b/examples/oscillator_sensi.rs new file mode 100644 index 0000000..f3e5597 --- /dev/null +++ b/examples/oscillator_sensi.rs @@ -0,0 +1,58 @@ +use cvode_wrap::*; + +fn main() { + 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 + ); + } +} diff --git a/example/plot.py b/examples/plot_oscillator.py similarity index 100% rename from example/plot.py rename to examples/plot_oscillator.py diff --git a/cvode-wrap/src/cvode.rs b/src/cvode.rs similarity index 82% rename from cvode-wrap/src/cvode.rs rename to src/cvode.rs index 571a4fd..527985f 100644 --- a/cvode-wrap/src/cvode.rs +++ b/src/cvode.rs @@ -1,54 +1,30 @@ -use std::{convert::TryInto, ffi::c_void, os::raw::c_int, pin::Pin, ptr::NonNull}; +//! Wrapper around cvode, without sensitivities + +use std::{convert::TryInto, os::raw::c_int, pin::Pin}; use sundials_sys::{SUNLinearSolver, SUNMatrix}; use crate::{ - check_flag_is_succes, check_non_null, AbsTolerance, LinearMultistepMethod, NVectorSerial, - NVectorSerialHeapAllocated, Realtype, Result, RhsResult, StepKind, + check_flag_is_succes, check_non_null, AbsTolerance, CvodeMemoryBlock, + CvodeMemoryBlockNonNullPtr, LinearMultistepMethod, NVectorSerial, NVectorSerialHeapAllocated, + Realtype, Result, RhsResult, StepKind, }; -#[repr(C)] -struct CvodeMemoryBlock { - _private: [u8; 0], -} - -#[repr(transparent)] -#[derive(Debug, Clone, Copy)] -struct CvodeMemoryBlockNonNullPtr { - ptr: NonNull, -} - -impl CvodeMemoryBlockNonNullPtr { - fn new(ptr: NonNull) -> Self { - Self { ptr } - } - - fn as_raw(self) -> *mut c_void { - self.ptr.as_ptr() as *mut c_void - } -} - -impl From> for CvodeMemoryBlockNonNullPtr { - fn from(x: NonNull) -> Self { - Self::new(x) - } -} - struct WrappingUserData { actual_user_data: UserData, f: F, } -/// The main struct of the crate. Wraps a sundials solver. +/// The ODE solver without sensitivities. /// -/// Args -/// ---- -/// `UserData` is the type of the supplementary arguments for the +/// # Type Arguments +/// +/// - `F` is the type of the right-hand side function +/// +/// - `UserData` is the type of the supplementary arguments for the /// right-hand-side. If unused, should be `()`. /// -/// `N` is the "problem size", that is the dimension of the state space. -/// -/// See [crate-level](`crate`) documentation for more. +/// - `N` is the "problem size", that is the dimension of the state space. pub struct Solver { mem: CvodeMemoryBlockNonNullPtr, y0: NVectorSerialHeapAllocated, @@ -58,9 +34,6 @@ pub struct Solver { user_data: Pin>>, } -/// The wrapping function. -/// -/// Internally used in [`wrap`]. extern "C" fn wrap_f( t: Realtype, y: *const NVectorSerial, @@ -88,6 +61,7 @@ impl Solver where F: Fn(Realtype, &[Realtype; N], &mut [Realtype; N], &UserData) -> RhsResult, { + /// Create a new solver. pub fn new( method: LinearMultistepMethod, f: F, @@ -170,6 +144,11 @@ where Ok(res) } + /// Takes a step according to `step_kind` (see [`StepKind`]). + /// + /// Returns a tuple `(t_out,&y(t_out))` where `t_out` is the time + /// reached by the solver as dictated by `step_kind`, and `y(t_out)` is an + /// array of the state variables at that time. pub fn step( &mut self, tout: Realtype, @@ -225,6 +204,7 @@ mod tests { 1e-4, AbsTolerance::Scalar(1e-4), (), - ).unwrap(); + ) + .unwrap(); } } diff --git a/cvode-wrap/src/cvode_sens.rs b/src/cvode_sens.rs similarity index 85% rename from cvode-wrap/src/cvode_sens.rs rename to src/cvode_sens.rs index d8e8838..13a34e9 100644 --- a/cvode-wrap/src/cvode_sens.rs +++ b/src/cvode_sens.rs @@ -1,76 +1,35 @@ -use std::{convert::TryInto, ffi::c_void, os::raw::c_int, pin::Pin, ptr::NonNull}; +//! Wrapper around cvodeS, with sensitivities + +use std::{convert::TryInto, os::raw::c_int, pin::Pin}; use sundials_sys::{SUNLinearSolver, SUNMatrix, CV_STAGGERED}; use crate::{ - check_flag_is_succes, check_non_null, AbsTolerance, LinearMultistepMethod, NVectorSerial, - NVectorSerialHeapAllocated, Realtype, Result, RhsResult, StepKind, + check_flag_is_succes, check_non_null, AbsTolerance, CvodeMemoryBlock, + CvodeMemoryBlockNonNullPtr, LinearMultistepMethod, NVectorSerial, NVectorSerialHeapAllocated, + Realtype, Result, RhsResult, SensiAbsTolerance, StepKind, }; -#[repr(C)] -struct CvodeMemoryBlock { - _private: [u8; 0], -} - -#[repr(transparent)] -#[derive(Debug, Clone, Copy)] -struct CvodeMemoryBlockNonNullPtr { - ptr: NonNull, -} - -impl CvodeMemoryBlockNonNullPtr { - fn new(ptr: NonNull) -> Self { - Self { ptr } - } - - fn as_raw(self) -> *mut c_void { - self.ptr.as_ptr() as *mut c_void - } -} - -pub enum SensiAbsTolerance { - Scalar([Realtype; N_SENSI]), - Vector([NVectorSerialHeapAllocated; N_SENSI]), -} - -impl SensiAbsTolerance { - pub fn scalar(atol: [Realtype; N_SENSI]) -> Self { - SensiAbsTolerance::Scalar(atol) - } - - pub fn vector(atol: &[[Realtype; SIZE]; N_SENSI]) -> Self { - SensiAbsTolerance::Vector( - array_init::from_iter( - atol.iter() - .map(|arr| NVectorSerialHeapAllocated::new_from(arr)), - ) - .unwrap(), - ) - } -} - -impl From> for CvodeMemoryBlockNonNullPtr { - fn from(x: NonNull) -> Self { - Self::new(x) - } -} - struct WrappingUserData { actual_user_data: UserData, f: F, fs: FS, } -/// The main struct of the crate. Wraps a sundials solver. +/// The ODE solver with sensitivities. /// -/// Args -/// ---- -/// `UserData` is the type of the supplementary arguments for the +/// # Type Arguments +/// +/// - `F` is the type of the right-hand side function +/// +/// - `FS` is the type of the sensitivities right-hand side function +/// +/// - `UserData` is the type of the supplementary arguments for the /// right-hand-side. If unused, should be `()`. /// -/// `N` is the "problem size", that is the dimension of the state space. +/// - `N` is the "problem size", that is the dimension of the state space. /// -/// See [crate-level](`crate`) documentation for more. +/// - `N_SENSI` is the number of sensitivities computed pub struct Solver { mem: CvodeMemoryBlockNonNullPtr, y0: NVectorSerialHeapAllocated, @@ -83,9 +42,6 @@ pub struct Solver { sensi_out_buffer: [NVectorSerialHeapAllocated; N_SENSI], } -/// The wrapping function. -/// -/// Internally used in [`wrap`]. extern "C" fn wrap_f( t: Realtype, y: *const NVectorSerial, @@ -168,6 +124,7 @@ where &UserData, ) -> RhsResult, { + /// Creates a new solver. #[allow(clippy::clippy::too_many_arguments)] pub fn new( method: LinearMultistepMethod, @@ -293,6 +250,12 @@ where Ok(res) } + /// Takes a step according to `step_kind` (see [`StepKind`]). + /// + /// Returns a tuple `(t_out,&y(t_out),[&dy_dp(tout)])` where `t_out` is the time + /// reached by the solver as dictated by `step_kind`, `y(t_out)` is an + /// array of the state variables at that time, and the i-th `dy_dp(tout)` is an array + /// of the sensitivities of all variables with respect to parameter i. #[allow(clippy::clippy::type_complexity)] pub fn step( &mut self, @@ -379,6 +342,7 @@ mod tests { AbsTolerance::scalar(1e-4), SensiAbsTolerance::scalar([1e-4; 4]), (), - ).unwrap(); + ) + .unwrap(); } } diff --git a/cvode-wrap/src/lib.rs b/src/lib.rs similarity index 65% rename from cvode-wrap/src/lib.rs rename to src/lib.rs index ab828a9..ffdd138 100644 --- a/cvode-wrap/src/lib.rs +++ b/src/lib.rs @@ -1,12 +1,19 @@ -use std::{os::raw::c_int, ptr::NonNull}; +//! A wrapper around cvode and cvodes from the sundials tool suite. +//! +//! Users should be mostly interested in [`SolverSensi`] and [`SolverNoSensi`]. + +use std::{ffi::c_void, os::raw::c_int, ptr::NonNull}; use sundials_sys::realtype; mod nvector; pub use nvector::{NVectorSerial, NVectorSerialHeapAllocated}; -pub mod cvode; -pub mod cvode_sens; +mod cvode; +mod cvode_sens; + +pub use cvode::Solver as SolverNoSensi; +pub use cvode_sens::Solver as SolverSensi; /// The floatting-point type sundials was compiled with pub type Realtype = realtype; @@ -80,6 +87,29 @@ impl AbsTolerance { } } +/// An enum representing the choice between scalars or vectors absolute tolerances +/// for sensitivities. +pub enum SensiAbsTolerance { + Scalar([Realtype; N_SENSI]), + Vector([NVectorSerialHeapAllocated; N_SENSI]), +} + +impl SensiAbsTolerance { + pub fn scalar(atol: [Realtype; N_SENSI]) -> Self { + SensiAbsTolerance::Scalar(atol) + } + + pub fn vector(atol: &[[Realtype; SIZE]; N_SENSI]) -> Self { + SensiAbsTolerance::Vector( + array_init::from_iter( + atol.iter() + .map(|arr| NVectorSerialHeapAllocated::new_from(arr)), + ) + .unwrap(), + ) + } +} + /// A short-hand for `std::result::Result` pub type Result = std::result::Result; @@ -94,3 +124,30 @@ fn check_flag_is_succes(flag: c_int, func_id: &'static str) -> Result<()> { Err(Error::ErrorCode { flag, func_id }) } } + +#[repr(C)] +struct CvodeMemoryBlock { + _private: [u8; 0], +} + +#[repr(transparent)] +#[derive(Debug, Clone, Copy)] +struct CvodeMemoryBlockNonNullPtr { + ptr: NonNull, +} + +impl CvodeMemoryBlockNonNullPtr { + fn new(ptr: NonNull) -> Self { + Self { ptr } + } + + fn as_raw(self) -> *mut c_void { + self.ptr.as_ptr() as *mut c_void + } +} + +impl From> for CvodeMemoryBlockNonNullPtr { + fn from(x: NonNull) -> Self { + Self::new(x) + } +} diff --git a/cvode-wrap/src/nvector.rs b/src/nvector.rs similarity index 100% rename from cvode-wrap/src/nvector.rs rename to src/nvector.rs