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
112        // 1) Build and sort event times
113        let mut ts = Vec::new();
114        ts.push(ti);
115        ts.push(tf);
116        for inf in infusions {
117            let t0 = inf.time();
118            let t1 = t0 + inf.duration();
119            if t0 > ti && t0 < tf {
120                ts.push(t0)
121            }
122            if t1 > ti && t1 < tf {
123                ts.push(t1)
124            }
125        }
126        ts.sort_by(|a, b| a.partial_cmp(b).unwrap());
127        ts.dedup_by(|a, b| (*a - *b).abs() < 1e-12);
128
129        // 2) March over each sub-interval
130        let mut current_t = ts[0];
131        for &next_t in &ts[1..] {
132            // prepare support and infusion rate for [current_t .. next_t]
133            let mut sp = V::from_vec(support_point.to_owned());
134            let mut rateiv = V::from_vec(vec![0.0; 3]);
135            for inf in infusions {
136                let s = inf.time();
137                let e = s + inf.duration();
138                if current_t >= s && next_t <= e {
139                    rateiv[inf.input()] += inf.amount() / inf.duration();
140                }
141            }
142
143            // advance the support-point to next_t
144            (self.seq_eq)(&mut sp, next_t, covariates);
145
146            // advance state by dt
147            let dt = next_t - current_t;
148            *x = (self.eq)(x, &sp, dt, rateiv, covariates);
149
150            current_t = next_t;
151        }
152
153        Ok(())
154    }
155
156    #[inline(always)]
157    fn process_observation(
158        &self,
159        support_point: &Vec<f64>,
160        observation: &Observation,
161        error_models: Option<&ErrorModels>,
162        _time: f64,
163        covariates: &Covariates,
164        x: &mut Self::S,
165        likelihood: &mut Vec<f64>,
166        output: &mut Self::P,
167    ) -> Result<(), PharmsolError> {
168        let mut y = V::zeros(self.get_nouteqs());
169        let out = &self.out;
170        (out)(
171            x,
172            &V::from_vec(support_point.clone()),
173            observation.time(),
174            covariates,
175            &mut y,
176        );
177        let pred = y[observation.outeq()];
178        let pred = observation.to_prediction(pred, x.as_slice().to_vec());
179        if let Some(error_models) = error_models {
180            likelihood.push(pred.likelihood(error_models)?);
181        }
182        output.add_prediction(pred);
183        Ok(())
184    }
185    #[inline(always)]
186    fn initial_state(&self, spp: &Vec<f64>, covariates: &Covariates, occasion_index: usize) -> V {
187        let init = &self.init;
188        let mut x = V::zeros(self.get_nstates());
189        if occasion_index == 0 {
190            (init)(&V::from_vec(spp.to_vec()), 0.0, covariates, &mut x);
191        }
192        x
193    }
194}
195
196impl Equation for Analytical {
197    fn estimate_likelihood(
198        &self,
199        subject: &Subject,
200        support_point: &Vec<f64>,
201        error_models: &ErrorModels,
202        cache: bool,
203    ) -> Result<f64, PharmsolError> {
204        _estimate_likelihood(self, subject, support_point, error_models, cache)
205    }
206    fn kind() -> crate::EqnKind {
207        crate::EqnKind::Analytical
208    }
209}
210fn spphash(spp: &[f64]) -> u64 {
211    spp.iter().fold(0, |acc, x| acc + x.to_bits())
212}
213
214#[inline(always)]
215#[cached(
216    ty = "UnboundCache<String, SubjectPredictions>",
217    create = "{ UnboundCache::with_capacity(100_000) }",
218    convert = r#"{ format!("{}{}", subject.id(), spphash(support_point)) }"#,
219    result = "true"
220)]
221fn _subject_predictions(
222    ode: &Analytical,
223    subject: &Subject,
224    support_point: &Vec<f64>,
225) -> Result<SubjectPredictions, PharmsolError> {
226    Ok(ode.simulate_subject(subject, support_point, None)?.0)
227}
228
229fn _estimate_likelihood(
230    ode: &Analytical,
231    subject: &Subject,
232    support_point: &Vec<f64>,
233    error_models: &ErrorModels,
234    cache: bool,
235) -> Result<f64, PharmsolError> {
236    let ypred = if cache {
237        _subject_predictions(ode, subject, support_point)
238    } else {
239        _subject_predictions_no_cache(ode, subject, support_point)
240    }?;
241    ypred.likelihood(error_models)
242}