feat: Improve error handling (#58)

* chore: Implement serialie/deserialize for ErrorModel

* refactor: Creating predictions from observations

Also adding tests to ErrorModel

* WIP

* WIP

* WIP, seems to be working

* Update BKE example results print

* Update bke.rs

* Handle ShapeError

* Handle CSV errors

* More errors

* Fix tests for exa

* Handle progress bar errors

* Remove unwrap

* Do not cache errors

* Cache results, impl From instead of Into for Prediction to Observation

* Convenience function for variance

* Provide more methods for calculating sigma and variance

* Handle errors in psi

* Error handling for likelihood

* Provide to_observation for Prediction

* ErrorPoly as a struct

Author: mhovd

benches/ode.rs

@@ -32,7 +32,10 @@ fn one_compartment() {
         },
         (2, 1),
     );
-    black_box(ode.estimate_predictions(&subject, &vec![0.3, 0.5, 0.1, 70.0]));
+    black_box(
+        ode.estimate_predictions(&subject, &vec![0.3, 0.5, 0.1, 70.0])
+            .unwrap(),
+    );
 }
 
 fn two_compartment() {
@@ -55,7 +58,10 @@ fn two_compartment() {
         },
         (2, 1),
     );
-    black_box(ode.estimate_predictions(&subject, &vec![0.3, 0.5, 0.1, 0.04, 0.08, 70.0]));
+    black_box(
+        ode.estimate_predictions(&subject, &vec![0.3, 0.5, 0.1, 0.04, 0.08, 70.0])
+            .unwrap(),
+    );
 }
 
 fn criterion_benchmark(c: &mut Criterion) {

benches/performance.rs

@@ -39,7 +39,9 @@ fn readme(n: usize) {
         (2, 1),
     );
     for _ in 0..n {
-        let op = ode.estimate_predictions(&subject, &vec![0.3, 0.5, 0.1, 70.0]);
+        let op = ode
+            .estimate_predictions(&subject, &vec![0.3, 0.5, 0.1, 70.0])
+            .unwrap();
         black_box(op);
     }
 }

examples/bke.rs

@@ -41,26 +41,45 @@ fn main() {
         (1, 1),
     );
 
-    let em = ErrorModel::additive((0.0, 0.05, 0.0, 0.0), 0.0);
-    let ll = an.estimate_likelihood(
-        &subject,
-        &vec![1.02282724609375, 194.51904296875],
-        &em,
-        false,
-    );
-    let op = an.estimate_predictions(&subject, &vec![1.02282724609375, 194.51904296875]);
-    println!(
-        "Analytical: \n-2ll:{:#?}\n{:#?}",
-        -2.0 * ll,
-        op.flat_predictions()
-    );
+    let em = ErrorModel::additive(ErrorPoly::new(0.0, 0.05, 0.0, 0.0), 0.0);
 
-    let ll = ode.estimate_likelihood(
-        &subject,
-        &vec![1.02282724609375, 194.51904296875],
-        &em,
-        false,
-    );
-    let op = ode.estimate_predictions(&subject, &vec![1.02282724609375, 194.51904296875]);
-    println!("ODE: \n-2ll:{:#?}\n{:#?}", -2.0 * ll, op.flat_predictions());
+    // Compute likelihoods and predictions for both models
+    let ll_an = an
+        .estimate_likelihood(
+            &subject,
+            &vec![1.02282724609375, 194.51904296875],
+            &em,
+            false,
+        )
+        .unwrap();
+    let op_an = an
+        .estimate_predictions(&subject, &vec![1.02282724609375, 194.51904296875])
+        .unwrap();
+
+    let ll_ode = ode
+        .estimate_likelihood(
+            &subject,
+            &vec![1.02282724609375, 194.51904296875],
+            &em,
+            false,
+        )
+        .unwrap();
+    let op_ode = ode
+        .estimate_predictions(&subject, &vec![1.02282724609375, 194.51904296875])
+        .unwrap();
+
+    // Display likelihoods side by side
+    println!("Likelihoods:");
+    println!("Analytical\tODE");
+    println!("{:.6}\t:{:.6}", -2.0 * ll_an, -2.0 * ll_ode);
+    println!();
+
+    // Display predictions
+    println!("Predictions:");
+    println!("Analytical\tODE\tDifference");
+    op_an
+        .flat_predictions()
+        .iter()
+        .zip(op_ode.flat_predictions())
+        .for_each(|(a, b)| println!("{:.9}\t{:.9}\t{:.9}", a, b, a - b));
 }

examples/error.rs

@@ -39,6 +39,8 @@ fn main() {
         .repeat(120, 0.1)
         .build();
 
-    let op = ode.estimate_predictions(&subject, &vec![0.1, 0.1, 0.1, 0.1, 70.0]);
+    let op = ode
+        .estimate_predictions(&subject, &vec![0.1, 0.1, 0.1, 0.1, 70.0])
+        .unwrap();
     dbg!(op);
 }

examples/exa.rs

@@ -68,8 +68,8 @@ fn main() {
     let params = vec![1.02282724609375, 194.51904296875];
 
     // Get predictions from both models
-    let dyn_predictions = dyn_ode.estimate_predictions(&subject, &params);
-    let ode_predictions = ode.estimate_predictions(&subject, &params);
+    let dyn_predictions = dyn_ode.estimate_predictions(&subject, &params).unwrap();
+    let ode_predictions = ode.estimate_predictions(&subject, &params).unwrap();
 
     // Check that predictions are the same
     let dyn_flat = dyn_predictions.flat_predictions();

examples/gendata.rs

@@ -62,7 +62,7 @@ fn main() {
 
     let mut data = vec![];
     for (i, spp) in support_points.iter().enumerate() {
-        let trajectories = sde.estimate_predictions(&subject, spp);
+        let trajectories = sde.estimate_predictions(&subject, spp).unwrap();
         let trajectory = trajectories.row(0);
         // dbg!(&trajectory);
         let mut sb = data::Subject::builder(format!("id{}", i)).bolus(0.0, 20.0, 0);
@@ -72,7 +72,8 @@ fn main() {
         data.push(sb.build());
     }
     let data = data::Data::new(data);
-    data.write_pmetrics(&File::create(Path::new("test.csv")).unwrap());
+    data.write_pmetrics(&File::create(Path::new("test.csv")).unwrap())
+        .unwrap();
     // let mut theta = Array2::zeros((1, 3));
     // theta[[0, 0]] = 0.7;
     // theta[[0, 1]] = 0.1;

examples/ode_readme.rs

@@ -35,6 +35,8 @@ fn main() {
         (2, 1),
     );
 
-    let op = ode.estimate_predictions(&subject, &vec![0.3, 0.5, 0.1, 70.0]);
+    let op = ode
+        .estimate_predictions(&subject, &vec![0.3, 0.5, 0.1, 70.0])
+        .unwrap();
     println!("{:#?}", op.flat_predictions());
 }

examples/pf.rs

@@ -28,10 +28,14 @@ fn main() {
         (2, 1),
         10000,
     );
-    let em = ErrorModel::additive((0.0, 0.05, 0.0, 0.0), 0.0);
+    let em = ErrorModel::additive(ErrorPoly::new(0.0, 0.05, 0.0, 0.0), 0.0);
 
-    let ll = sde.estimate_likelihood(&subject, &vec![1.0], &em, false);
+    let ll = sde
+        .estimate_likelihood(&subject, &vec![1.0], &em, false)
+        .unwrap();
 
-    dbg!(sde.estimate_likelihood(&subject, &vec![1.0], &em, false));
+    dbg!(sde
+        .estimate_likelihood(&subject, &vec![1.0], &em, false)
+        .unwrap());
     println!("{ll:#?}");
 }

examples/sde.rs

@@ -106,8 +106,8 @@ fn main() {
         .repeat(5, 0.2)
         .build();
 
-    let ode_predictions = ode.estimate_predictions(&subject, &vec![1.0]);
-    let sde_predictions = sde.estimate_predictions(&subject, &vec![1.0, 0.0]);
+    let ode_predictions = ode.estimate_predictions(&subject, &vec![1.0]).unwrap();
+    let sde_predictions = sde.estimate_predictions(&subject, &vec![1.0, 0.0]).unwrap();
 
     let mut sde_flat_predictions: Vec<Vec<f64>> = Vec::new();
     for trajectory in sde_predictions.rows() {
@@ -141,8 +141,8 @@ fn main() {
         0.0,
     ];
 
-    let ode_predictions = ode.estimate_predictions(&subject, &spp_ode);
-    let sde_predictions = sde.estimate_predictions(&subject, &spp_sde);
+    let ode_predictions = ode.estimate_predictions(&subject, &spp_ode).unwrap();
+    let sde_predictions = sde.estimate_predictions(&subject, &spp_sde).unwrap();
 
     let mut sde_flat_predictions: Vec<Vec<f64>> = Vec::new();
     for trajectory in sde_predictions.rows() {
@@ -164,8 +164,8 @@ fn main() {
     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);
-    let sde_predictions = sde.estimate_predictions(subject, &spp_sde);
+    let ode_predictions = ode.estimate_predictions(subject, &spp_ode).unwrap();
+    let sde_predictions = sde.estimate_predictions(subject, &spp_sde).unwrap();
 
     let mut sde_flat_predictions: Vec<Vec<f64>> = Vec::new();
     for trajectory in sde_predictions.rows() {

src/data/builder.rs

@@ -115,7 +115,7 @@ impl SubjectBuilder {
         time: f64,
         value: f64,
         outeq: usize,
-        errorpoly: Option<(f64, f64, f64, f64)>,
+        errorpoly: Option<ErrorPoly>,
         ignore: bool,
     ) -> Self {
         let observation = Observation::new(time, value, outeq, errorpoly, ignore);
@@ -299,7 +299,13 @@ mod tests {
             .bolus(0.0, 50.0, 0)
             .observation(1.0, 45.3, 0)
             .observation(2.0, 40.1, 0)
-            .observation_with_error(3.0, 36.5, 0, Some((0.1, 0.05, 0.0, 0.0)), false)
+            .observation_with_error(
+                3.0,
+                36.5,
+                0,
+                Some(ErrorPoly::new(0.1, 0.05, 0.0, 0.0)),
+                false,
+            )
             .bolus(4.0, 50.0, 0)
             .repeat(1, 12.0) // Repeat bolus at 16.0
             .reset()