Merge branch 'main' into addl

.github/workflows/time.yml

@@ -0,0 +1,20 @@
+name: Time
+
+on:
+  workflow_dispatch:
+
+jobs:
+  cargo-run:
+    runs-on: self-hosted
+    env:
+      NPCORE_CONFIG_TUI: false  # Disable TUI
+    steps:
+    - name: Checkout Repository
+      uses: actions/checkout@v2
+
+    - name: Run
+      run: |
+        /usr/bin/time -v cargo run --release --example bimodal_ke 2>&1 | tee time-output.txt  # Redirect output to a file
+
+    - name: Results
+      run: cat time-output.txt  # Print the content of the file

.gitignore

@@ -1,6 +1,6 @@
 /target
 /Cargo.lock
-/log
+*.log
 theta.csv
 cycles.csv
 pred.csv
@@ -14,10 +14,12 @@ simulation_output.csv
 meta*.csv
 /examples/rosuva/*
 /examples/iohexol/*
+/examples/vori/*
 /examples/data/iohexol*
 /examples/data/rosuva*
+/examples/data/vori*
 /.idea
 stop
 .vscode
-
-*.f90
+*.f90
+settings.json

Cargo.toml

@@ -1,51 +1,53 @@
 [package]
 name = "npcore"
-version = "0.1.1"
+version = "0.1.3"
 edition = "2021"
 authors = [
     "Julián D. Otálvaro <juliandavid347@gmail.com>",
     "Markus Hovd",
     "Michael Neely",
     "Walter Yamada",
 ]
-description = "Rust Library with the building blocks needed to create new Non-Parametric algorithms and its integration with Pmetrics."
+description = "Rust library with the building blocks needed to create new Non-Parametric algorithms and its integration with Pmetrics."
 license = "GPL-3.0"
 documentation = "https://lapkb.github.io/NPcore/npcore/"
 repository = "https://github.com/LAPKB/NPcore"
 exclude = [".github/*", ".vscode/*"]
 
 [dependencies]
 dashmap = "5.5.3"
-#cached = "0.26.2"
 lazy_static = "1.4.0"
 csv = "1.2.1"
 ndarray = { version = "0.15.6", features = ["rayon"] }
 serde = "1.0.188"
 serde_derive = "1.0.188"
+serde_json = "1.0.66"
 sobol_burley = "0.5.0"
 toml = { version = "0.8.1", features = ["preserve_order"] }
 ode_solvers = "0.3.7"
 ndarray-stats = "0.5.1"
 linfa-linalg = "0.1.0"
-log = "0.4.20"
-log4rs = "1.2.0"
 rayon = "1.8.0"
 eyre = "0.6.8"
-ratatui = { version = "0.23.0", features = ["crossterm"] }
+ratatui = { version = "0.25.0", features = ["crossterm"] }
 crossterm = "0.27.0"
 tokio = { version = "1.32.0", features = ["sync", "rt"] }
 ndarray-csv = "0.5.2"
 rawpointer = "0.2.1"
-argmin = "0.8.1"
-itertools = "0.11.0"
-faer-core = { version = "0.11.0", features = [] }
+argmin = { version = "0.9.0", features = [] }
+itertools = "0.12.0"
+faer-core = { version = "0.15.0", features = [] }
 # faer-lu = "0.9"
-faer-qr = "0.11.0"
+faer-qr = "0.16.0"
 # faer-cholesky = "0.9"
 # faer-svd = "0.9"
-dyn-stack = "0.9.0"
-faer = { version = "0.11.0", features = ["nalgebra", "ndarray"] }
-
+argmin-math = { version = "0.3.0", features = ["ndarray_v0_15-nolinalg-serde"] }
+dyn-stack = "0.10.0"
+faer = { version = "0.15.0", features = ["nalgebra", "ndarray"] }
+tracing = "0.1.40"
+tracing-subscriber = { version = "0.3.17", features = ["env-filter", "fmt", "time"] }
+chrono = "0.4"
+config = "0.13"
 
 [profile.release]
 codegen-units = 1

README.md

@@ -1,30 +1,43 @@
 # NPcore
 [![Build](https://github.com/LAPKB/NPcore/actions/workflows/rust.yml/badge.svg)](https://github.com/LAPKB/NPcore/actions/workflows/rust.yml)
 [![Security Audit](https://github.com/LAPKB/NPcore/actions/workflows/security_audit.yml/badge.svg)](https://github.com/LAPKB/NPcore/actions/workflows/security_audit.yml)
+![crates.io](https://img.shields.io/crates/v/npcore.svg)
 
-Rust Library with the building blocks needed to create new Non-Parametric algorithms and its integration with [Pmetrics]([https://link-url-here.org](https://github.com/LAPKB/Pmetrics)).
+Rust library with the building blocks to create and implement new non-parametric algorithms for population pharmacokinetic modelling and their integration with [Pmetrics](https://github.com/LAPKB/Pmetrics).
 
 ## Implemented functionality
 
 * Solver for ODE-based population pharmacokinetic models
 * Supports the Pmetrics data format for seamless integration
-* Basic NPAG implementation for parameter estimation
 * Covariate support, carry-forward or linear interpolation
 * Option to cache results for improved speed
 * Powerful simulation engine 
 * Informative Terminal User Interface (TUI)
 
+## Available algoritms
+
+This project aims to implement several algorithms for non-parametric population pharmacokinetic modelling.
+
+- [x] Non Parametric Adaptive Grid (NPAG)
+    - [Yamada et al (2021)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7823953/)
+    - [Neely et al (2012)](https://pubmed.ncbi.nlm.nih.gov/22722776/)
+- [x] Non Parametric Optimal Design (NPOD)
+  - [Otalvaro et al (2023)](https://pubmed.ncbi.nlm.nih.gov/36478350/)
+  - [Leary et al (2003)](https://www.page-meeting.org/default.asp?abstract=421)
+- [ ] Non Parametric Simulated Annealing (NPSA)
+  - [Chen et al (2023)](https://arxiv.org/abs/2301.12656)
+
+In the future we also aim to support parametric algorithms, such as the Iterative 2-Stage Bayesian (IT2B)
 
 ## Examples
 
 There are two examples using NPAG implemented in this repository, `bimodal_ke` and `two_eq_lag`. 
 
-You may run them with the following command
+You may run them with the following command, e.g.
 ```
-cargo run --example example_name --release
+cargo run --example bimodal_ke --release
 ```
-Be sure to replace `example_name` with the respective example.
-Look at the corresponding `examples/_example_/*.toml`-file to change the configuration for each run.
+Look at the corresponding `examples/.../*.toml`-file to change the configuration for each run.
 
 ## Documentation
 

examples/bimodal_ke/config.toml

@@ -1,16 +1,19 @@
 [paths]
 data = "examples/data/bimodal_ke.csv"
-log_out = "log/bimodal_ke.log"
-#prior_dist = "theta_bimodal_ke.csv"
+log = "bimodal_ke.log"
+#prior = "theta_bimodal_ke.csv"
 
 [config]
 cycles = 1024
 engine = "NPAG"
 init_points = 2129
 seed = 347
 tui = true
-pmetrics_outputs = true
+output = true
 cache = true
+idelta = 0.1
+log_level = "info"
+
 
 [random]
 Ke = [0.001, 3.0]

examples/bimodal_ke/main.rs

@@ -8,94 +8,118 @@ use npcore::prelude::{
 };
 use ode_solvers::*;
 
+// Constants for the absolute and relative tolerance for the dynamic steps used for solving the ODEs
 const ATOL: f64 = 1e-4;
 const RTOL: f64 = 1e-4;
 
+// Define the state vector, which must be equal to the number of compartments in the model
+// These are re-exported from the `nalgebra`-crate by `ode_solvers`, see https://github.com/srenevey/ode-solvers?tab=readme-ov-file#type-alias-definition
+// In brief, for up to 6 compartments, use VectorN<f64>, N being the number of compartments.
+// For more than 6 compartments, use `nalgebra::SVector<f64, N>`, where N is the number of compartments.
+type State = Vector1<f64>;
+
+// Time uses f64 precision
+type Time = f64;
+
+// This is the main structure for the model
+// It holds the parameters (defined in config.toml), the scenarios (i.e. datafile), the (possible) infusions and the covariates if any
 #[derive(Debug, Clone)]
-struct Model<'a> {
-    ke: f64,
-    _v: f64,
-    _scenario: &'a Scenario,
+struct Model {
+    params: HashMap<String, f64>,
+    _scenario: Scenario,
     infusions: Vec<Infusion>,
-    cov: Option<&'a HashMap<String, CovLine>>,
+    cov: Option<HashMap<String, CovLine>>,
 }
 
-type State = Vector1<f64>;
-type Time = f64;
+// This is a helper function to get the parameter value by name
+impl Model {
+    pub fn get_param(&self, str: &str) -> f64 {
+        *self.params.get(str).unwrap()
+    }
+}
 
-impl ode_solvers::System<State> for Model<'_> {
+impl ode_solvers::System<State> for Model {
     fn system(&self, t: Time, y: &State, dy: &mut State) {
-        let ke = self.ke;
-
-        let _lag = 0.0;
+        // Get the parameters from the model
+        let ke = self.get_param("ke");
 
+        // Get the infusions that are active at time `t`
         let mut rateiv = [0.0];
         for infusion in &self.infusions {
             if t >= infusion.time && t <= (infusion.dur + infusion.time) {
                 rateiv[infusion.compartment] += infusion.amount / infusion.dur;
             }
         }
+        // The ordinary differential equations (ODEs) are defined here
+        // This example is a one-compartmental model with first-order elimination, and intravenous infusions
 
-        ///////////////////// USER DEFINED ///////////////
-
+        ////// ODE //////
         dy[0] = -ke * y[0] + rateiv[0];
-
-        //////////////// END USER DEFINED ////////////////
     }
 }
 
 #[derive(Debug, Clone)]
 struct Ode {}
 
-impl Predict for Ode {
-    fn predict(&self, params: Vec<f64>, scenario: &Scenario) -> Vec<f64> {
-        let mut system = Model {
-            ke: params[0],
-            _v: params[1],
-            _scenario: scenario,
-            infusions: vec![],
-            cov: None,
-        };
-        // let scenario = scenario.reorder_with_lag(vec![]);
-        let mut yout = vec![];
-        let mut x = State::new(0.0);
-        let mut index: usize = 0;
-        for block in &scenario.blocks {
-            system.cov = Some(&block.covs);
-            for event in &block.events {
-                if event.evid == 1 {
-                    if event.dur.unwrap_or(0.0) > 0.0 {
-                        //infusion
-                        system.infusions.push(Infusion {
-                            time: event.time,
-                            dur: event.dur.unwrap(),
-                            amount: event.dose.unwrap(),
-                            compartment: event.input.unwrap() - 1,
-                        });
-                    } else {
-                        //dose
-                        x[event.input.unwrap() - 1] += event.dose.unwrap();
-                    }
-                } else if event.evid == 0 {
-                    //obs
-                    yout.push(x[event.outeq.unwrap() - 1] / params[1]);
-                }
-                if let Some(next_time) = scenario.times.get(index + 1) {
-                    //TODO: use the last dx as the initial one for the next simulation.
-                    let mut stepper =
-                        Dopri5::new(system.clone(), event.time, *next_time, 1e-3, x, RTOL, ATOL);
-                    let _res = stepper.integrate();
-                    let y = stepper.y_out();
-                    x = *y.last().unwrap();
-                }
-                index += 1;
-            }
+impl<'a> Predict<'a> for Ode {
+    type Model = Model;
+    type State = State;
+    fn initial_system(&self, params: &Vec<f64>, scenario: Scenario) -> (Self::Model, Scenario) {
+        let params = HashMap::from([("ke".to_string(), params[0]), ("v".to_string(), params[1])]);
+        (
+            Model {
+                params,
+                _scenario: scenario.clone(), //TODO remove
+                infusions: vec![],
+                cov: None,
+            },
+            scenario.reorder_with_lag(vec![(0.0, 1)]),
+        )
+    }
+
+    // This function is used to get the output from the model, defined by the output equations (outeq) supplied by the user
+    fn get_output(&self, time: f64, x: &Self::State, system: &Self::Model, outeq: usize) -> f64 {
+        // Get parameters from the model also used for calculating the output equations
+        let v = system.get_param("v");
+        #[allow(unused_variables)]
+        let t = time;
+        match outeq {
+            1 => x[0] / v, // Concentration of the central compartment defined by the amount of drug, x[0], divided by the volume of the central compartment, v
+            _ => panic!("Invalid output equation"),
+        }
+    }
+
+    // Set the initial state of the compartments
+    fn initial_state(&self) -> State {
+        State::default()
+    }
+    // Add any possible infusions
+    fn add_infusion(&self, system: &mut Self::Model, infusion: Infusion) {
+        system.infusions.push(infusion);
+    }
+    // Add any possible covariates
+    fn add_covs(&self, system: &mut Self::Model, cov: Option<HashMap<String, CovLine>>) {
+        system.cov = cov;
+    }
+    // Add any possible doses
+    fn add_dose(&self, state: &mut Self::State, dose: f64, compartment: usize) {
+        state[compartment] += dose;
+    }
+    // Perform a "step" of the model, i.e. solve the ODEs from the current time to the next time
+    // In the next step, we use this result as the initial state
+    fn state_step(&self, x: &mut Self::State, system: &Self::Model, time: f64, next_time: f64) {
+        if time >= next_time {
+            panic!("time error")
         }
-        yout
+        let mut stepper = Dopri5::new(system.clone(), time, next_time, 1e-3, *x, RTOL, ATOL);
+        let _res = stepper.integrate();
+        let y = stepper.y_out();
+        *x = *y.last().unwrap();
     }
 }
 
 fn main() -> Result<()> {
+    // Main entrypoint, see `entrypoints.rs` for more details
     let _result = start(
         Engine::new(Ode {}),
         "examples/bimodal_ke/config.toml".to_string(),