Sde dynamic stepsize

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "pharmsol"
-version = "0.7.7"
+version = "0.7.6"
 edition = "2021"
 authors = ["Julián D. Otálvaro <juliandavid347@gmail.com>", "Markus Hovd"]
 description = "Rust library for solving analytic and ode-defined pharmacometric models."

examples/new_gd.rs

@@ -0,0 +1,60 @@
+use csv::ReaderBuilder;
+use pharmsol::*;
+use prelude::{data::read_pmetrics, simulator::Prediction};
+
+fn main() {
+    // path to theta
+    let path = "../PMcore/outputs/theta.csv";
+    // read theta into an Array2<f64>
+    let mut rdr = ReaderBuilder::new().from_path(path).unwrap();
+
+    let mut spps = vec![];
+    for result in rdr.records() {
+        let record = result.unwrap();
+        let mut row = vec![];
+        for field in record.iter() {
+            row.push(field.parse::<f64>().unwrap());
+        }
+        spps.push(row);
+    }
+
+    let sde = equation::SDE::new(
+        |x, p, _t, dx, _rateiv, _cov| {
+            // automatically defined
+            fetch_params!(p, ke0, _ske);
+            // let ke0 = 1.2;
+            dx[1] = -x[1] + ke0;
+            let ke = x[1];
+            // user defined
+            dx[0] = -ke * x[0];
+        },
+        |p, d| {
+            fetch_params!(p, _ke0, ske);
+            d[1] = ske;
+        },
+        |_p| lag! {},
+        |_p| fa! {},
+        |p, _t, _cov, x| {
+            fetch_params!(p, ke0, ske);
+            x[1] = ke0;
+        },
+        |x, p, _t, _cov, y| {
+            fetch_params!(p, _ke0, _ske);
+            y[0] = x[0] / 50.0;
+        },
+        (2, 1),
+        1,
+    );
+
+    let data = read_pmetrics("../PMcore/examples/iov/test.csv").unwrap();
+
+    for (i, spp) in spps.iter().enumerate() {
+        for (j, subject) in data.get_subjects().iter().enumerate() {
+            let trajectories: ndarray::Array2<Prediction> =
+                sde.estimate_predictions(&subject, &spp);
+            let trajectory = trajectories.row(0);
+            println!("{}, {}", i, j);
+            dbg!(trajectory);
+        }
+    }
+}

examples/sde.rs

@@ -165,7 +165,7 @@ fn main() {
     let ode = three_c_ode();
     let sde = three_c_sde();
 
-    let data = read_pmetrics("../PMcore/examples/w_vanco_sde/test.csv").unwrap();
+    let data = read_pmetrics("../PMcore/examples/vanco_sde/data.csv").unwrap();
     let subject = data.get_subject("51").unwrap();
 
     let ode_predictions = ode.estimate_predictions(&subject, &spp_ode);

examples/sde_paper/main.rs

@@ -0,0 +1,73 @@
+use pharmsol::equation;
+use pmcore::prelude::*;
+use rand::rngs::StdRng;
+use rand::SeedableRng;
+use rand::weighted::WeightedIndex;
+use rand_distr::(Distribution, Normal);
+
+fn model() -> equation::SDE {
+    let sde = equation::SDE::new(
+        drift:
+        diffeq: |x, p, _t, dx, rateiv, _cov| {
+            // fetch_cov!(cov, t, wt);
+            fetch_params!(p, ke, _v);
+            dx[0] = -ke * x[0] + rateiv[0];
+        },
+        lag: |_p| lag! {},
+        fa: |_p| fa! {},
+        init: |_p, _t, _cov, _x| {},
+        out: |x, p, _t, _cov, y| {
+            fetch_params!(p, _ke, v);
+            y[0] = x[0] / v;
+        },
+        neqs: (1, 1),
+    );
+
+}
+
+let fn sample_k0(rng: &mut StdRng, n1: Normal<f64>,      ) -> f64 {
+
+}
+
+
+fn sample_v()
+
+
+
+
+
+const N_SAMPLES: usize = 100;
+
+fn main() {
+    let m1: f64 = 0.5;
+    let s1: f64 = 0.05;
+    let m2: f64 = 1.5;
+    let s2: f64 = 0.15;
+
+    let n1 = Normal::new(m1,s1).unwrap();
+    let n2 = Normal::new(m2,s2).unwrap();
+
+    // let weights = [0.5, 0.5];
+    // let dist = WeightedIndex::new(&weights).unwrap();
+
+    let mut rng = seed_from_u64(state: 42);
+    let mut_k0_pop: Vec<f64> = Vec::new();
+    let v_pop: Vec<f64> = Vec::new();
+
+    let n3 = Normal::new(mean: 0.0, std_dev: 1.0).unwrap();
+
+    for _ in 0..N_SAMPLES {
+        k0_pop.push(sample_k0(&mut rng, n1, n2));
+        v_pop.push(sample)
+    }
+
+    // plot the distributions
+    let trace
+
+
+    let k0_dist = 0.5 * rand_distr::Normal::new(mean: m1, std_dev: s2)
+        + 0.5 * rand_distr::Normal::new(mean: m2, std_dev: s2).unwrap();
+    let seed: u64 = 42;
+    let rng: StdRng = rand::rngs::StdRng::seed_from_u64(state:: seed);
+
+}

src/simulator/equation/sde/em.rs

@@ -1,6 +1,7 @@
 use crate::{
     data::Covariates,
     simulator::{Diffusion, Drift},
+    Infusion,
 };
 use nalgebra::DVector;
 use rand::rng;
@@ -13,7 +14,7 @@ pub struct EM {
     params: DVector<f64>,
     state: DVector<f64>,
     cov: Covariates,
-    rateiv: DVector<f64>,
+    infusions: Vec<Infusion>,
     rtol: f64,
     atol: f64,
     max_step: f64,
@@ -28,7 +29,7 @@ impl EM {
         params: DVector<f64>,
         initial_state: DVector<f64>,
         cov: Covariates,
-        rateiv: DVector<f64>,
+        infusions: Vec<Infusion>,
         rtol: f64,
         atol: f64,
     ) -> Self {
@@ -38,7 +39,7 @@ impl EM {
             params,
             state: initial_state,
             cov,
-            rateiv,
+            infusions,
             rtol,
             atol,
             max_step: 0.1,  // Can be made configurable
@@ -64,27 +65,35 @@ impl EM {
         new_dt
     }
 
-    fn euler_maruyama_step(&self, time: f64, dt: f64, state: &mut DVector<f64>) {
+    fn euler_maruyama_step(&self, time: f64, dt: f64, state: &mut DVector<f64>, sigma: f64) {
         let n = state.len();
+        let mut rateiv = DVector::from_vec(vec![0.0, 0.0, 0.0]);
+        //TODO: This should be pre-calculated
+        for infusion in &self.infusions {
+            if time >= infusion.time() && time <= infusion.duration() + infusion.time() {
+                rateiv[infusion.input()] += infusion.amount() / infusion.duration();
+            }
+        }
         let mut drift_term = DVector::zeros(n);
         (self.drift)(
             state,
             &self.params,
             time,
             &mut drift_term,
-            self.rateiv.clone(),
+            rateiv,
             &self.cov,
         );
 
         let mut diffusion_term = DVector::zeros(n);
         (self.diffusion)(&self.params, &mut diffusion_term);
 
-        let mut rng = rng();
-        let normal_dist = Normal::new(0.0, 1.0).unwrap();
+        let mut _rng = rng();
+        let _normal_dist = Normal::new(0.0, 1.0).unwrap();
 
         for i in 0..n {
             state[i] +=
-                drift_term[i] * dt + diffusion_term[i] * normal_dist.sample(&mut rng) * dt.sqrt();
+                // drift_term[i] * dt + diffusion_term[i] * normal_dist.sample(&mut rng) * dt.sqrt();
+                drift_term[i] * dt + diffusion_term[i] * sigma * dt.sqrt();
         }
     }
 
@@ -94,17 +103,24 @@ impl EM {
         let safety = 0.9;
         let mut times = vec![t0];
         let mut solution = vec![self.state.clone()];
+        let mut sigma = 0.00000001;
+
 
+        let mut rng = rng();
+        let normal_dist = Normal::new(0.0, 1.0).unwrap();
+        // sigma = normal_dist.sample(&mut rng); // oops!!! this makes one sigma for ENTIRE time
         while t < tf {
             let mut y1 = self.state.clone();
             let mut y2 = self.state.clone();
 
+            sigma = normal_dist.sample(&mut rng); // here is correct, only affects one step of the variable time step
+
             // Single step
-            self.euler_maruyama_step(t, dt, &mut y1);
+            self.euler_maruyama_step(t, dt, &mut y1, sigma);
 
             // Two half steps
-            self.euler_maruyama_step(t, dt / 2.0, &mut y2);
-            self.euler_maruyama_step(t + dt / 2.0, dt / 2.0, &mut y2);
+            self.euler_maruyama_step(t, dt / 2.0, &mut y2, sigma);
+            self.euler_maruyama_step(t + dt / 2.0, dt / 2.0, &mut y2, sigma);
 
             let error = self.calculate_error(&y1, &y2);
 

src/simulator/equation/sde/mod.rs

@@ -34,21 +34,14 @@ pub(crate) fn simulate_sde_event(
     if ti == tf {
         return x;
     }
-    let mut rateiv = V::from_vec(vec![0.0, 0.0, 0.0]);
-    //TODO: This should be pre-calculated
-    for infusion in infusions {
-        if tf >= infusion.time() && tf <= infusion.duration() + infusion.time() {
-            rateiv[infusion.input()] += infusion.amount() / infusion.duration();
-        }
-    }
 
     let mut sde = em::EM::new(
         *drift,
         *difussion,
         DVector::from_column_slice(support_point),
         x,
         cov.clone(),
-        rateiv,
+        infusions.to_vec(),
         1e-2,
         1e-2,
     );
@@ -109,8 +102,29 @@ impl Predictions for Array2<Prediction> {
     fn get_predictions(&self) -> Vec<Prediction> {
         //TODO: This is only returning the first particle, not the best, not the worst, THE FIRST
         // CHANGE THIS
-        let row = self.row(0).to_vec();
-        row
+        // let row = self.row(0).to_vec();
+        // row
+        // Make this return the mean prediction across all particles
+        if self.is_empty() || self.ncols() == 0 {
+            return Vec::new();
+        }
+
+        let mut result = Vec::with_capacity(self.ncols());
+
+        for col in 0..self.ncols() {
+            let column = self.column(col);
+            let mean_prediction: f64 = column
+                .iter()
+                .map(|pred: &Prediction| pred.prediction())
+                .sum::<f64>()
+                / self.nrows() as f64;
+
+                let mut prediction = column.first().unwrap().clone();
+                prediction.set_prediction(mean_prediction);
+                result.push(prediction);
+        }
+
+        result
     }
 }
 
@@ -209,8 +223,16 @@ impl EquationPriv for SDE {
         // q = pdf.(Distributions.Normal(0, 0.5), e)
         if let Some(em) = error_model {
             let mut q: Vec<f64> = Vec::with_capacity(self.nparticles);
-
+            //
+            // wmy centering_function is a running Chi^2 w/exp = support point
+            // let centering_function = p.pred[2]; // move this inside the iteration.
+            //
+            // pred.iter().for_each(|p| q.push(p.state[4] * p.likelihood(em)));
+            //
+            // note: Above doesn't compile b/c pred is private; but also: I'm not sure if pred is the state, and state is x ???
+            //
             pred.iter().for_each(|p| q.push(p.likelihood(em)));
+            //
             let sum_q: f64 = q.iter().sum();
             let w: Vec<f64> = q.iter().map(|qi| qi / sum_q).collect();
             let i = sysresample(&w);

src/simulator/likelihood.rs

@@ -194,6 +194,9 @@ impl Prediction {
     pub fn prediction(&self) -> f64 {
         self.prediction
     }
+    pub fn set_prediction(&mut self, prediction: f64) {
+        self.prediction = prediction;
+    }
     pub fn outeq(&self) -> usize {
         self.outeq
     }