pharmsol/simulator/equation/analytical/
mod.rs

1pub mod one_compartment_models;
2pub mod three_compartment_models;
3pub mod two_compartment_models;
4
5use diffsol::{NalgebraContext, Vector, VectorHost};
6pub use one_compartment_models::*;
7pub use three_compartment_models::*;
8pub use two_compartment_models::*;
9
10use crate::mapping::Mappings;
11use crate::PharmsolError;
12use crate::{
13    data::Covariates, simulator::*, Equation, EquationPriv, EquationTypes, Observation, Subject,
14};
15use cached::proc_macro::cached;
16use cached::UnboundCache;
17
18/// Model equation using analytical solutions.
19///
20/// This implementation uses closed-form analytical solutions for the model
21/// equations rather than numerical integration.
22#[repr(C)]
23#[derive(Clone, Debug)]
24pub struct Analytical {
25    eq: AnalyticalEq,
26    seq_eq: SecEq,
27    lag: Lag,
28    fa: Fa,
29    init: Init,
30    out: Out,
31    neqs: Neqs,
32    mappings: Option<Mappings>,
33}
34
35impl Analytical {
36    /// Create a new Analytical equation model.
37    ///
38    /// # Parameters
39    /// - `eq`: The analytical equation function
40    /// - `seq_eq`: The secondary equation function
41    /// - `lag`: The lag time function
42    /// - `fa`: The fraction absorbed function
43    /// - `init`: The initial state function
44    /// - `out`: The output equation function
45    /// - `neqs`: The number of states and output equations
46    pub fn new(
47        eq: AnalyticalEq,
48        seq_eq: SecEq,
49        lag: Lag,
50        fa: Fa,
51        init: Init,
52        out: Out,
53        neqs: Neqs,
54    ) -> Self {
55        Self {
56            eq,
57            seq_eq,
58            lag,
59            fa,
60            init,
61            out,
62            neqs,
63            mappings: None,
64        }
65    }
66}
67
68impl EquationTypes for Analytical {
69    type S = V;
70    type P = SubjectPredictions;
71}
72
73impl EquationPriv for Analytical {
74    // #[inline(always)]
75    // fn get_init(&self) -> &Init {
76    //     &self.init
77    // }
78
79    // #[inline(always)]
80    // fn get_out(&self) -> &Out {
81    //     &self.out
82    // }
83
84    // #[inline(always)]
85    // fn get_lag(&self, spp: &[f64]) -> Option<HashMap<usize, f64>> {
86    //     Some((self.lag)(&V::from_vec(spp.to_owned())))
87    // }
88
89    // #[inline(always)]
90    // fn get_fa(&self, spp: &[f64]) -> Option<HashMap<usize, f64>> {
91    //     Some((self.fa)(&V::from_vec(spp.to_owned())))
92    // }
93
94    #[inline(always)]
95    fn lag(&self) -> &Lag {
96        &self.lag
97    }
98
99    #[inline(always)]
100    fn fa(&self) -> &Fa {
101        &self.fa
102    }
103
104    #[inline(always)]
105    fn get_nstates(&self) -> usize {
106        self.neqs.0
107    }
108
109    #[inline(always)]
110    fn get_nouteqs(&self) -> usize {
111        self.neqs.1
112    }
113    #[inline(always)]
114    fn solve(
115        &self,
116        x: &mut Self::S,
117        support_point: &Vec<f64>,
118        covariates: &Covariates,
119        infusions: &Vec<Infusion>,
120        ti: f64,
121        tf: f64,
122    ) -> Result<(), PharmsolError> {
123        if ti == tf {
124            return Ok(());
125        }
126
127        // 1) Build and sort event times
128        let mut ts = Vec::new();
129        ts.push(ti);
130        ts.push(tf);
131        for inf in infusions {
132            let t0 = inf.time();
133            let t1 = t0 + inf.duration();
134            if t0 > ti && t0 < tf {
135                ts.push(t0)
136            }
137            if t1 > ti && t1 < tf {
138                ts.push(t1)
139            }
140        }
141        ts.sort_by(|a, b| a.partial_cmp(b).unwrap());
142        ts.dedup_by(|a, b| (*a - *b).abs() < 1e-12);
143
144        // 2) March over each sub-interval
145        let mut current_t = ts[0];
146        for &next_t in &ts[1..] {
147            // prepare support and infusion rate for [current_t .. next_t]
148            let mut sp = V::from_vec(support_point.to_owned(), NalgebraContext);
149            let mut rateiv = V::from_vec(vec![0.0; 3], NalgebraContext);
150            for inf in infusions {
151                let s = inf.time();
152                let e = s + inf.duration();
153                if current_t >= s && next_t <= e {
154                    rateiv[inf.input()] += inf.amount() / inf.duration();
155                }
156            }
157
158            // advance the support-point to next_t
159            (self.seq_eq)(&mut sp, next_t, covariates);
160
161            // advance state by dt
162            let dt = next_t - current_t;
163            *x = (self.eq)(x, &sp, dt, rateiv, covariates);
164
165            current_t = next_t;
166        }
167
168        Ok(())
169    }
170
171    #[inline(always)]
172    fn process_observation(
173        &self,
174        support_point: &Vec<f64>,
175        observation: &Observation,
176        error_models: Option<&ErrorModels>,
177        _time: f64,
178        covariates: &Covariates,
179        x: &mut Self::S,
180        likelihood: &mut Vec<f64>,
181        output: &mut Self::P,
182    ) -> Result<(), PharmsolError> {
183        let mut y = V::zeros(self.get_nouteqs(), NalgebraContext);
184        let out = &self.out;
185        (out)(
186            x,
187            &V::from_vec(support_point.clone(), NalgebraContext),
188            observation.time(),
189            covariates,
190            &mut y,
191        );
192        let pred = y[observation.outeq()];
193        let pred = observation.to_prediction(pred, x.as_slice().to_vec());
194        if let Some(error_models) = error_models {
195            likelihood.push(pred.likelihood(error_models)?);
196        }
197        output.add_prediction(pred);
198        Ok(())
199    }
200    #[inline(always)]
201    fn initial_state(&self, spp: &Vec<f64>, covariates: &Covariates, occasion_index: usize) -> V {
202        let init = &self.init;
203        let mut x = V::zeros(self.get_nstates(), NalgebraContext);
204        if occasion_index == 0 {
205            (init)(
206                &V::from_vec(spp.to_vec(), NalgebraContext),
207                0.0,
208                covariates,
209                &mut x,
210            );
211        }
212        x
213    }
214}
215
216impl Equation for Analytical {
217    fn estimate_likelihood(
218        &self,
219        subject: &Subject,
220        support_point: &Vec<f64>,
221        error_models: &ErrorModels,
222        cache: bool,
223    ) -> Result<f64, PharmsolError> {
224        _estimate_likelihood(self, subject, support_point, error_models, cache)
225    }
226    fn kind() -> crate::EqnKind {
227        crate::EqnKind::Analytical
228    }
229    fn mappings_ref(&self) -> Option<&Mappings> {
230        self.mappings.as_ref()
231    }
232    fn mappings_mut(&mut self) -> &mut Mappings {
233        self.mappings.get_or_insert_with(Mappings::new)
234    }
235}
236fn spphash(spp: &[f64]) -> u64 {
237    spp.iter().fold(0, |acc, x| acc + x.to_bits())
238}
239
240#[inline(always)]
241#[cached(
242    ty = "UnboundCache<String, SubjectPredictions>",
243    create = "{ UnboundCache::with_capacity(100_000) }",
244    convert = r#"{ format!("{}{}", subject.id(), spphash(support_point)) }"#,
245    result = "true"
246)]
247fn _subject_predictions(
248    ode: &Analytical,
249    subject: &Subject,
250    support_point: &Vec<f64>,
251) -> Result<SubjectPredictions, PharmsolError> {
252    Ok(ode.simulate_subject(subject, support_point, None)?.0)
253}
254
255fn _estimate_likelihood(
256    ode: &Analytical,
257    subject: &Subject,
258    support_point: &Vec<f64>,
259    error_models: &ErrorModels,
260    cache: bool,
261) -> Result<f64, PharmsolError> {
262    let ypred = if cache {
263        _subject_predictions(ode, subject, support_point)
264    } else {
265        _subject_predictions_no_cache(ode, subject, support_point)
266    }?;
267    ypred.likelihood(error_models)
268}
pharmsol/simulator/equation/analytical/mod.rs

pharmsol/simulator/equation/analytical/
mod.rs
