pharmsol 0.27.1

Rust library for solving analytic and ode-defined pharmacometric models.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106

fn main() -> Result<(), pharmsol::PharmsolError> {
    use pharmsol::{prelude::*, Parameters};

    let analytical = analytical! {
        name: "one_cmt_iv",
        params: [ke, v],
        states: [central],
        outputs: [cp],
        routes: [
            infusion(iv) -> central,
        ],
        structure: one_compartment,
        out: |x, _p, _t, _cov, y| {
            y[cp] = x[central] / v;
        },
    };

    let ode = ode! {
        name: "one_cmt_iv",
        params: [ke, v],
        states: [central],
        outputs: [cp],
        routes: [
            infusion(iv) -> central,
        ],
        diffeq: |x, _p, _t, dx, _cov| {
            dx[central] = -ke * x[central];
        },
        out: |x, _p, _t, _cov, y| {
            y[cp] = x[central] / v;
        },
    };

    // Create a subject using route and output labels directly.
    let subject = Subject::builder("Nikola Tesla")
        .infusion(0., 500.0, "iv", 0.5)
        .observation(0.5, 1.645, "cp")
        .observation(1., 1.216, "cp")
        .observation(2., 0.462, "cp")
        .observation(3., 0.169, "cp")
        .observation(4., 0.063, "cp")
        .observation(6., 0.009, "cp")
        .observation(8., 0.001, "cp")
        .missing_observation(12.0, "cp")
        .build();

    // Define the assay error models once by label and reuse them across both
    // equations.
    let ems = AssayErrorModels::new().add(
        "cp",
        AssayErrorModel::additive(ErrorPoly::new(0.0, 0.05, 0.0, 0.0), 0.0),
    )?;

    // Define the parameter values for the simulations
    let ke = 1.022; // Elimination rate constant
    let v = 194.0; // Volume of distribution
    let analytical_parameters = Parameters::with_model(&analytical, [("ke", ke), ("v", v)])
        .expect("valid named parameters");
    let ode_parameters =
        Parameters::with_model(&ode, [("ke", ke), ("v", v)]).expect("valid named parameters");

    // Compute likelihoods and predictions for both models
    let analytical_likelihoods =
        analytical.estimate_log_likelihood(&subject, &analytical_parameters, &ems)?;

    let analytical_predictions =
        analytical.estimate_predictions(&subject, &analytical_parameters)?;

    let ode_likelihoods = ode.estimate_log_likelihood(&subject, &ode_parameters, &ems)?;

    let ode_predictions = ode.estimate_predictions(&subject, &ode_parameters)?;

    // Print comparison table
    println!("\n┌───────────┬─────────────────┬─────────────────┬─────────────────────┐");
    println!("│           │   Analytical    │       ODE       │     Difference      │");
    println!("├───────────┼─────────────────┼─────────────────┼─────────────────────┤");
    println!(
        "│ Log-Likeli│ {:>15.6} │ {:>15.6} │ {:>19.2e} │",
        analytical_likelihoods,
        ode_likelihoods,
        analytical_likelihoods - ode_likelihoods
    );
    println!("├───────────┼─────────────────┼─────────────────┼─────────────────────┤");
    println!("│   Time    │   Prediction    │   Prediction    │                     │");
    println!("├───────────┼─────────────────┼─────────────────┼─────────────────────┤");

    let times = analytical_predictions.flat_times();
    let analytical_preds = analytical_predictions.flat_predictions();
    let ode_preds = ode_predictions.flat_predictions();

    for ((t, a), b) in times
        .iter()
        .zip(analytical_preds.iter())
        .zip(ode_preds.iter())
    {
        let diff = a - b;
        println!(
            "│ {:>9.2} │ {:>15.9} │ {:>15.9} │ {:>19.2e} │",
            t, a, b, diff
        );
    }

    println!("└───────────┴─────────────────┴─────────────────┴─────────────────────┘\n");

    Ok(())
}