Merge branch 'main' into NPOD

.gitignore

@@ -19,5 +19,4 @@ meta*.csv
 /.idea
 stop
 .vscode
-
 *.f90

Cargo.toml

@@ -8,15 +8,14 @@ authors = [
     "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/*"]
+exclude = [".github/*", ".vscode/*", "debug.rs", "examples/faer_qr.rs", "examples/drusano/*", "examples/two_eq_lag/*", "examples/vori/*"]
 
 [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"] }
@@ -27,25 +26,26 @@ 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.24.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 = { version = "0.8.1", features = [] }
 itertools = "0.11.0"
-faer-core = { version = "0.11.0", features = [] }
+faer-core = { version = "0.14.0", features = [] }
 # faer-lu = "0.9"
-faer-qr = "0.11.0"
+faer-qr = "0.14.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.14.0", features = ["nalgebra", "ndarray"] }
+tracing = "0.1.40"
+tracing-subscriber = { version = "0.3.17", features = ["env-filter", "fmt", "time"] }
+chrono = "0.4"
 
 
 [profile.release]

README.md

@@ -2,18 +2,29 @@
 [![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)
 
-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 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/)
+- [ ] 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)
 
 ## Examples
 

examples/bimodal_ke/config.toml

@@ -11,6 +11,8 @@ seed = 347
 tui = true
 pmetrics_outputs = true
 cache = true
+idelta = 0.1
+log_level = "info"
 
 [random]
 Ke = [0.001, 3.0]

examples/bimodal_ke/main.rs

@@ -2,32 +2,34 @@ use std::collections::HashMap;
 
 use eyre::Result;
 use npcore::prelude::{
+    datafile,
     datafile::{CovLine, Infusion, Scenario},
     predict::{Engine, Predict},
-    start,
+    settings, start, start_with_data,
 };
 use ode_solvers::*;
 
 const ATOL: f64 = 1e-4;
 const RTOL: f64 = 1e-4;
 
+type State = Vector1<f64>;
+type Time = f64;
 #[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>>,
+}
+impl Model {
+    pub fn get_param(&self, str: &str) -> f64 {
+        *self.params.get(str).unwrap()
+    }
 }
 
-type State = Vector1<f64>;
-type Time = f64;
-
-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;
+        let ke = self.get_param("ke");
 
         let mut rateiv = [0.0];
         for infusion in &self.infusions {
@@ -47,51 +49,58 @@ impl ode_solvers::System<State> for Model<'_> {
 #[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,
+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,
-        };
-        // 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;
-            }
+        },
+        scenario.reorder_with_lag(vec![(0.0, 1)])
+    )
+
+    }
+    fn get_output(&self, x: &Self::State, system: &Self::Model, outeq: usize) -> f64 {
+        let v = system.get_param("v");
+        match outeq {
+            1 => x[0] / v,
+            _ => panic!("Invalid output equation"),
         }
-        yout
+    }
+    fn initial_state(&self) -> State {
+        State::default()
+    }
+    fn add_infusion(&self, mut system: Self::Model, infusion: Infusion) -> Model {
+        system.infusions.push(infusion);
+        system
+    }
+    fn add_covs(&self, mut system: Self::Model, cov: Option<HashMap<String, CovLine>>) -> Model {
+        system.cov = cov;
+        system
+    }
+    fn add_dose(&self, mut state: Self::State, dose: f64, compartment: usize) -> Self::State {
+        state[compartment] += dose;
+        state
+    }
+    fn state_step(
+        &self,
+        mut x: Self::State,
+        system: Self::Model,
+        time: f64,
+        next_time: f64,
+    ) -> State {
+        if time == next_time {
+            return x;
+        }
+        let mut stepper = Dopri5::new(system, time, next_time, 1e-3, x, RTOL, ATOL);
+        let _res = stepper.integrate();
+        let y = stepper.y_out();
+        x = *y.last().unwrap();
+        x
     }
 }
 
@@ -103,3 +112,23 @@ fn main() -> Result<()> {
 
     Ok(())
 }
+
+#[allow(dead_code)]
+fn new_entry_test() -> Result<()> {
+    let settings_path = "examples/bimodal_ke/config.toml".to_string();
+    let settings = settings::run::read(settings_path);
+    let mut scenarios = datafile::parse(&settings.parsed.paths.data).unwrap();
+    if let Some(exclude) = &settings.parsed.config.exclude {
+        for val in exclude {
+            scenarios.remove(val.as_integer().unwrap() as usize);
+        }
+    }
+
+    let _result = start_with_data(
+        Engine::new(Ode {}),
+        "examples/bimodal_ke/config.toml".to_string(),
+        scenarios,
+    )?;
+
+    Ok(())
+}

examples/debug.rs

@@ -92,7 +92,7 @@ impl Predict for Ode {
                                 let y = stepper.y_out();
                                 x = *y.last().unwrap();
                             } else if *next_time < event.time {
-                                log::error!("Panic: Next event's time is in the past!");
+                                tracing::error!("Panic: Next event's time is in the past!");
                                 panic!("Panic: Next event's time is in the past!");
                             }
                         }
@@ -103,7 +103,7 @@ impl Predict for Ode {
                             //     Rk4::new(system.clone(), event.time, x, lag_time, 0.1);
                             if let Some(next_time) = scenario.times.get(index + 1) {
                                 if *next_time < lag_time {
-                                    log::error!("Panic: lag time overpasses next observation, not implemented. Stopping.");
+                                    tracing::error!("Panic: lag time overpasses next observation, not implemented. Stopping.");
                                     panic!("Panic: lag time overpasses next observation, not implemented. Stopping.");
                                 }
                                 let mut stepper = Dopri5::new(
@@ -139,7 +139,7 @@ impl Predict for Ode {
                                 let y = stepper.y_out();
                                 x = *y.last().unwrap();
                             } else if *next_time > event.time {
-                                log::error!("Panic: Next event's time is in the past!");
+                                tracing::error!("Panic: Next event's time is in the past!");
                                 panic!("Panic: Next event's time is in the past!");
                             }
                         }
@@ -164,7 +164,7 @@ impl Predict for Ode {
                             let y = stepper.y_out();
                             x = *y.last().unwrap();
                         } else if *next_time < event.time {
-                            log::error!("Panic: Next event's time is in the past!");
+                            tracing::error!("Panic: Next event's time is in the past!");
                             panic!("Panic: Next event's time is in the past!");
                         }
                     }

examples/vori/main.rs

@@ -127,8 +127,11 @@ impl Predict for Ode {
                         let y = stepper.y_out();
                         x = *y.last().unwrap();
                     } else if event.time > *next_time {
-                        log::error!("next time is in the past!");
-                        log::error!("event_time: {}\nnext_time: {}", event.time, *next_time);
+                        tracing::error!(
+                            "The next event before the current: event_time: {}, next_time: {}",
+                            event.time,
+                            *next_time
+                        );
                     }
                 }
                 index += 1;

src/algorithms.rs

@@ -1,7 +1,7 @@
-use crate::prelude::{self, output::NPCycle, settings::run::Data};
+use crate::prelude::{self, settings::run::Data};
 
 use output::NPResult;
-use prelude::*;
+use prelude::{datafile::Scenario, *};
 use simulation::predict::{Engine, Predict};
 use tokio::sync::mpsc;
 
@@ -23,25 +23,21 @@ pub trait Algorithm {
 pub fn initialize_algorithm<S>(
     engine: Engine<S>,
     settings: Data,
-    tx: mpsc::UnboundedSender<NPCycle>,
+    scenarios: Vec<Scenario>,
+    tx: mpsc::UnboundedSender<Comm>,
 ) -> Box<dyn Algorithm>
 where
-    S: Predict + std::marker::Sync + Clone + 'static,
+    S: Predict<'static> + std::marker::Sync + Clone + 'static,
 {
     if std::path::Path::new("stop").exists() {
         match std::fs::remove_file("stop") {
-            Ok(_) => log::info!("Removed previous stop file"),
+            Ok(_) => tracing::info!("Removed previous stop file"),
             Err(err) => panic!("Unable to remove previous stop file: {}", err),
         }
     }
     let ranges = settings.computed.random.ranges.clone();
     let theta = initialization::sample_space(&settings, &ranges);
-    let mut scenarios = datafile::parse(&settings.parsed.paths.data).unwrap();
-    if let Some(exclude) = &settings.parsed.config.exclude {
-        for val in exclude {
-            scenarios.remove(val.as_integer().unwrap() as usize);
-        }
-    }
+
     //This should be a macro, so it can automatically expands as soon as we add a new option in the Type Enum
     match settings.parsed.config.engine.as_str() {
         "NPAG" => Box::new(npag::NPAG::new(