pharmsol/simulator/equation/analytical/
mod.rs

1pub mod one_compartment_models;
2pub mod three_compartment_models;
3pub mod two_compartment_models;
4
5pub use one_compartment_models::*;
6pub use three_compartment_models::*;
7pub use two_compartment_models::*;
8
9use crate::PharmsolError;
10use crate::{
11    data::Covariates, simulator::*, Equation, EquationPriv, EquationTypes, Observation, Subject,
12};
13use cached::proc_macro::cached;
14use cached::UnboundCache;
15
16/// Model equation using analytical solutions.
17///
18/// This implementation uses closed-form analytical solutions for the model
19/// equations rather than numerical integration.
20#[derive(Clone, Debug)]
21pub struct Analytical {
22    eq: AnalyticalEq,
23    seq_eq: SecEq,
24    lag: Lag,
25    fa: Fa,
26    init: Init,
27    out: Out,
28    neqs: Neqs,
29}
30
31impl Analytical {
32    /// Create a new Analytical equation model.
33    ///
34    /// # Parameters
35    /// - `eq`: The analytical equation function
36    /// - `seq_eq`: The secondary equation function
37    /// - `lag`: The lag time function
38    /// - `fa`: The fraction absorbed function
39    /// - `init`: The initial state function
40    /// - `out`: The output equation function
41    /// - `neqs`: The number of states and output equations
42    pub fn new(
43        eq: AnalyticalEq,
44        seq_eq: SecEq,
45        lag: Lag,
46        fa: Fa,
47        init: Init,
48        out: Out,
49        neqs: Neqs,
50    ) -> Self {
51        Self {
52            eq,
53            seq_eq,
54            lag,
55            fa,
56            init,
57            out,
58            neqs,
59        }
60    }
61}
62
63impl EquationTypes for Analytical {
64    type S = V;
65    type P = SubjectPredictions;
66}
67
68impl EquationPriv for Analytical {
69    // #[inline(always)]
70    // fn get_init(&self) -> &Init {
71    //     &self.init
72    // }
73
74    // #[inline(always)]
75    // fn get_out(&self) -> &Out {
76    //     &self.out
77    // }
78
79    #[inline(always)]
80    fn get_lag(&self, spp: &[f64]) -> Option<HashMap<usize, f64>> {
81        Some((self.lag)(&V::from_vec(spp.to_owned())))
82    }
83
84    #[inline(always)]
85    fn get_fa(&self, spp: &[f64]) -> Option<HashMap<usize, f64>> {
86        Some((self.fa)(&V::from_vec(spp.to_owned())))
87    }
88
89    #[inline(always)]
90    fn get_nstates(&self) -> usize {
91        self.neqs.0
92    }
93
94    #[inline(always)]
95    fn get_nouteqs(&self) -> usize {
96        self.neqs.1
97    }
98    #[inline(always)]
99    fn solve(
100        &self,
101        x: &mut Self::S,
102        support_point: &Vec<f64>,
103        covariates: &Covariates,
104        infusions: &Vec<Infusion>,
105        ti: f64,
106        tf: f64,
107    ) -> Result<(), PharmsolError> {
108        if ti == tf {
109            return Ok(());
110        }
111        let mut support_point = V::from_vec(support_point.to_owned());
112        let mut rateiv = V::from_vec(vec![0.0, 0.0, 0.0]);
113        //TODO: This should be pre-calculated
114        for infusion in infusions {
115            if tf >= infusion.time() && tf <= infusion.duration() + infusion.time() {
116                rateiv[infusion.input()] += infusion.amount() / infusion.duration();
117            }
118        }
119        (self.seq_eq)(&mut support_point, tf, covariates);
120        *x = (self.eq)(x, &support_point, tf - ti, rateiv, covariates);
121        Ok(())
122    }
123    #[inline(always)]
124    fn process_observation(
125        &self,
126        support_point: &Vec<f64>,
127        observation: &Observation,
128        error_model: Option<&ErrorModel>,
129        _time: f64,
130        covariates: &Covariates,
131        x: &mut Self::S,
132        likelihood: &mut Vec<f64>,
133        output: &mut Self::P,
134    ) -> Result<(), PharmsolError> {
135        let mut y = V::zeros(self.get_nouteqs());
136        let out = &self.out;
137        (out)(
138            x,
139            &V::from_vec(support_point.clone()),
140            observation.time(),
141            covariates,
142            &mut y,
143        );
144        let pred = y[observation.outeq()];
145        let pred = observation.to_prediction(pred, x.as_slice().to_vec());
146        if let Some(error_model) = error_model {
147            likelihood.push(pred.likelihood(error_model)?);
148        }
149        output.add_prediction(pred);
150        Ok(())
151    }
152    #[inline(always)]
153    fn initial_state(&self, spp: &Vec<f64>, covariates: &Covariates, occasion_index: usize) -> V {
154        let init = &self.init;
155        let mut x = V::zeros(self.get_nstates());
156        if occasion_index == 0 {
157            (init)(&V::from_vec(spp.to_vec()), 0.0, covariates, &mut x);
158        }
159        x
160    }
161}
162
163impl Equation for Analytical {
164    fn estimate_likelihood(
165        &self,
166        subject: &Subject,
167        support_point: &Vec<f64>,
168        error_model: &ErrorModel,
169        cache: bool,
170    ) -> Result<f64, PharmsolError> {
171        _estimate_likelihood(self, subject, support_point, error_model, cache)
172    }
173}
174fn spphash(spp: &[f64]) -> u64 {
175    spp.iter().fold(0, |acc, x| acc + x.to_bits())
176}
177
178#[inline(always)]
179#[cached(
180    ty = "UnboundCache<String, SubjectPredictions>",
181    create = "{ UnboundCache::with_capacity(100_000) }",
182    convert = r#"{ format!("{}{}", subject.id(), spphash(support_point)) }"#,
183    result = "true"
184)]
185fn _subject_predictions(
186    ode: &Analytical,
187    subject: &Subject,
188    support_point: &Vec<f64>,
189) -> Result<SubjectPredictions, PharmsolError> {
190    Ok(ode.simulate_subject(subject, support_point, None)?.0)
191}
192
193fn _estimate_likelihood(
194    ode: &Analytical,
195    subject: &Subject,
196    support_point: &Vec<f64>,
197    error_model: &ErrorModel,
198    cache: bool,
199) -> Result<f64, PharmsolError> {
200    let ypred = if cache {
201        _subject_predictions(ode, subject, support_point)
202    } else {
203        _subject_predictions_no_cache(ode, subject, support_point)
204    }?;
205    ypred.likelihood(error_model)
206}