pharmsol/simulator/equation/analytical/
mod.rs1pub 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#[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 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)]
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 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 let mut current_t = ts[0];
131 for &next_t in &ts[1..] {
132 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 (self.seq_eq)(&mut sp, next_t, covariates);
145
146 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}