pmcore/algorithms/
npag.rs

1use crate::algorithms::{Status, StopReason};
2use crate::prelude::algorithms::Algorithms;
3
4pub use crate::routines::estimation::ipm::burke;
5pub use crate::routines::estimation::qr;
6use crate::routines::settings::Settings;
7
8use crate::routines::output::{cycles::CycleLog, cycles::NPCycle, NPResult};
9use crate::structs::psi::{calculate_psi, Psi};
10use crate::structs::theta::Theta;
11use crate::structs::weights::Weights;
12
13use anyhow::bail;
14use anyhow::Result;
15use pharmsol::prelude::{
16    data::{AssayErrorModels, Data},
17    simulator::Equation,
18};
19
20use pharmsol::prelude::AssayErrorModel;
21
22use crate::routines::initialization;
23
24use crate::routines::expansion::adaptative_grid::adaptative_grid;
25
26const THETA_E: f64 = 1e-4; // Convergence criteria
27const THETA_G: f64 = 1e-4; // Objective function convergence criteria
28const THETA_F: f64 = 1e-2;
29const THETA_D: f64 = 1e-4;
30
31#[derive(Debug)]
32pub struct NPAG<E: Equation + Send + 'static> {
33    equation: E,
34    ranges: Vec<(f64, f64)>,
35    psi: Psi,
36    theta: Theta,
37    lambda: Weights,
38    w: Weights,
39    eps: f64,
40    last_objf: f64,
41    objf: f64,
42    f0: f64,
43    f1: f64,
44    cycle: usize,
45    gamma_delta: Vec<f64>,
46    error_models: AssayErrorModels,
47    status: Status,
48    cycle_log: CycleLog,
49    data: Data,
50    settings: Settings,
51}
52
53impl<E: Equation + Send + 'static> Algorithms<E> for NPAG<E> {
54    fn new(settings: Settings, equation: E, data: Data) -> Result<Box<Self>, anyhow::Error> {
55        Ok(Box::new(Self {
56            equation,
57            ranges: settings.parameters().ranges(),
58            psi: Psi::new(),
59            theta: Theta::new(),
60            lambda: Weights::default(),
61            w: Weights::default(),
62            eps: 0.2,
63            last_objf: -1e30,
64            objf: f64::NEG_INFINITY,
65            f0: -1e30,
66            f1: f64::default(),
67            cycle: 0,
68            gamma_delta: vec![0.1; settings.errormodels().len()],
69            error_models: settings.errormodels().clone(),
70            status: Status::Continue,
71            cycle_log: CycleLog::new(),
72            settings,
73            data,
74        }))
75    }
76
77    fn equation(&self) -> &E {
78        &self.equation
79    }
80    fn into_npresult(&self) -> Result<NPResult<E>> {
81        NPResult::new(
82            self.equation.clone(),
83            self.data.clone(),
84            self.theta.clone(),
85            self.psi.clone(),
86            self.w.clone(),
87            -2. * self.objf,
88            self.cycle,
89            self.status.clone(),
90            self.settings.clone(),
91            self.cycle_log.clone(),
92        )
93    }
94
95    fn settings(&self) -> &Settings {
96        &self.settings
97    }
98
99    fn data(&self) -> &Data {
100        &self.data
101    }
102
103    fn get_prior(&self) -> Theta {
104        initialization::sample_space(&self.settings).unwrap()
105    }
106
107    fn likelihood(&self) -> f64 {
108        self.objf
109    }
110
111    fn increment_cycle(&mut self) -> usize {
112        self.cycle += 1;
113        self.cycle
114    }
115
116    fn cycle(&self) -> usize {
117        self.cycle
118    }
119
120    fn set_theta(&mut self, theta: Theta) {
121        self.theta = theta;
122    }
123
124    fn theta(&self) -> &Theta {
125        &self.theta
126    }
127
128    fn psi(&self) -> &Psi {
129        &self.psi
130    }
131
132    fn evaluation(&mut self) -> Result<Status> {
133        tracing::info!("Objective function = {:.4}", -2.0 * self.objf);
134        tracing::debug!("Support points: {}", self.theta.nspp());
135
136        self.error_models.iter().for_each(|(outeq, em)| {
137            if AssayErrorModel::None == *em {
138                return;
139            }
140            tracing::debug!(
141                "Error model for outeq {}: {:.2}",
142                outeq,
143                em.factor().unwrap_or_default()
144            );
145        });
146
147        tracing::debug!("EPS = {:.4}", self.eps);
148        // Increasing objf signals instability or model misspecification.
149        if self.last_objf > self.objf + 1e-4 {
150            tracing::warn!(
151                "Objective function decreased from {:.4} to {:.4} (delta = {})",
152                -2.0 * self.last_objf,
153                -2.0 * self.objf,
154                -2.0 * self.last_objf - -2.0 * self.objf
155            );
156        }
157
158        let psi = self.psi.matrix();
159        let w = &self.w;
160        if (self.last_objf - self.objf).abs() <= THETA_G && self.eps > THETA_E {
161            self.eps /= 2.;
162            if self.eps <= THETA_E {
163                let pyl = psi * w.weights();
164                self.f1 = pyl.iter().map(|x| x.ln()).sum();
165                if (self.f1 - self.f0).abs() <= THETA_F {
166                    tracing::info!("The model converged after {} cycles", self.cycle,);
167                    self.set_status(Status::Stop(StopReason::Converged));
168                    self.log_cycle_state();
169                    return Ok(self.status().clone());
170                } else {
171                    self.f0 = self.f1;
172                    self.eps = 0.2;
173                }
174            }
175        }
176
177        // Stop if we have reached maximum number of cycles
178        if self.cycle >= self.settings.config().cycles {
179            tracing::warn!("Maximum number of cycles reached");
180            self.set_status(Status::Stop(StopReason::MaxCycles));
181            self.log_cycle_state();
182            return Ok(self.status().clone());
183        }
184
185        // Stop if stopfile exists
186        if std::path::Path::new("stop").exists() {
187            tracing::warn!("Stopfile detected - breaking");
188            self.set_status(Status::Stop(StopReason::Stopped));
189            self.log_cycle_state();
190            return Ok(self.status().clone());
191        }
192
193        // Continue with normal operation
194        self.set_status(Status::Continue);
195        self.log_cycle_state();
196        Ok(self.status().clone())
197    }
198
199    fn estimation(&mut self) -> Result<()> {
200        self.psi = calculate_psi(
201            &self.equation,
202            &self.data,
203            &self.theta,
204            &self.error_models,
205            self.cycle == 1 && self.settings.config().progress,
206        )?;
207
208        if let Err(err) = self.validate_psi() {
209            bail!(err);
210        }
211
212        (self.lambda, _) = match burke(&self.psi) {
213            Ok((lambda, objf)) => (lambda, objf),
214            Err(err) => {
215                bail!("Error in IPM during estimation: {:?}", err);
216            }
217        };
218        Ok(())
219    }
220
221    fn condensation(&mut self) -> Result<()> {
222        // Filter out the support points with lambda < max(lambda)/1000
223
224        let max_lambda = self
225            .lambda
226            .iter()
227            .fold(f64::NEG_INFINITY, |acc, x| x.max(acc));
228
229        let mut keep = Vec::<usize>::new();
230        for (index, lam) in self.lambda.iter().enumerate() {
231            if lam > max_lambda / 1000_f64 {
232                keep.push(index);
233            }
234        }
235        if self.psi.matrix().ncols() != keep.len() {
236            tracing::debug!(
237                "Lambda (max/1000) dropped {} support point(s)",
238                self.psi.matrix().ncols() - keep.len(),
239            );
240        }
241
242        self.theta.filter_indices(keep.as_slice());
243        self.psi.filter_column_indices(keep.as_slice());
244
245        //Rank-Revealing Factorization
246        let (r, perm) = qr::qrd(&self.psi)?;
247
248        let mut keep = Vec::<usize>::new();
249
250        // The minimum between the number of subjects and the actual number of support points
251        let keep_n = self.psi.matrix().ncols().min(self.psi.matrix().nrows());
252        for i in 0..keep_n {
253            let test = r.col(i).norm_l2();
254            let r_diag_val = r.get(i, i);
255            let ratio = r_diag_val / test;
256            if ratio.abs() >= 1e-8 {
257                keep.push(*perm.get(i).unwrap());
258            }
259        }
260
261        // If a support point is dropped, log it as a debug message
262        if self.psi.matrix().ncols() != keep.len() {
263            tracing::debug!(
264                "QR decomposition dropped {} support point(s)",
265                self.psi.matrix().ncols() - keep.len(),
266            );
267        }
268
269        // Filter to keep only the support points (rows) that are in the `keep` vector
270        self.theta.filter_indices(keep.as_slice());
271        // Filter to keep only the support points (columns) that are in the `keep` vector
272        self.psi.filter_column_indices(keep.as_slice());
273
274        self.validate_psi()?;
275        (self.lambda, self.objf) = match burke(&self.psi) {
276            Ok((lambda, objf)) => (lambda, objf),
277            Err(err) => {
278                return Err(anyhow::anyhow!(
279                    "Error in IPM during condensation: {:?}",
280                    err
281                ));
282            }
283        };
284        self.w = self.lambda.clone();
285        Ok(())
286    }
287
288    fn optimizations(&mut self) -> Result<()> {
289        self.error_models
290            .clone()
291            .iter_mut()
292            .filter_map(|(outeq, em)| {
293                if em.optimize() {
294                    Some((outeq, em))
295                } else {
296                    None
297                }
298            })
299            .try_for_each(|(outeq, em)| -> Result<()> {
300                // OPTIMIZATION
301
302                let gamma_up = em.factor()? * (1.0 + self.gamma_delta[outeq]);
303                let gamma_down = em.factor()? / (1.0 + self.gamma_delta[outeq]);
304
305                let mut error_model_up = self.error_models.clone();
306                error_model_up.set_factor(outeq, gamma_up)?;
307
308                let mut error_model_down = self.error_models.clone();
309                error_model_down.set_factor(outeq, gamma_down)?;
310
311                let psi_up = calculate_psi(
312                    &self.equation,
313                    &self.data,
314                    &self.theta,
315                    &error_model_up,
316                    false,
317                )?;
318                let psi_down = calculate_psi(
319                    &self.equation,
320                    &self.data,
321                    &self.theta,
322                    &error_model_down,
323                    false,
324                )?;
325
326                let (lambda_up, objf_up) = match burke(&psi_up) {
327                    Ok((lambda, objf)) => (lambda, objf),
328                    Err(err) => {
329                        bail!("Error in IPM during optim: {:?}", err);
330                    }
331                };
332                let (lambda_down, objf_down) = match burke(&psi_down) {
333                    Ok((lambda, objf)) => (lambda, objf),
334                    Err(err) => {
335                        bail!("Error in IPM during optim: {:?}", err);
336                    }
337                };
338                if objf_up > self.objf {
339                    self.error_models.set_factor(outeq, gamma_up)?;
340                    self.objf = objf_up;
341                    self.gamma_delta[outeq] *= 4.;
342                    self.lambda = lambda_up;
343                    self.psi = psi_up;
344                }
345                if objf_down > self.objf {
346                    self.error_models.set_factor(outeq, gamma_down)?;
347                    self.objf = objf_down;
348                    self.gamma_delta[outeq] *= 4.;
349                    self.lambda = lambda_down;
350                    self.psi = psi_down;
351                }
352                self.gamma_delta[outeq] *= 0.5;
353                if self.gamma_delta[outeq] <= 0.01 {
354                    self.gamma_delta[outeq] = 0.1;
355                }
356                Ok(())
357            })?;
358
359        Ok(())
360    }
361
362    fn expansion(&mut self) -> Result<()> {
363        adaptative_grid(&mut self.theta, self.eps, &self.ranges, THETA_D)?;
364        Ok(())
365    }
366
367    fn set_status(&mut self, status: Status) {
368        self.status = status;
369    }
370
371    fn status(&self) -> &Status {
372        &self.status
373    }
374
375    fn log_cycle_state(&mut self) {
376        let state = NPCycle::new(
377            self.cycle,
378            -2. * self.objf,
379            self.error_models.clone(),
380            self.theta.clone(),
381            self.theta.nspp(),
382            (self.last_objf - self.objf).abs(),
383            self.status.clone(),
384        );
385        self.cycle_log.push(state);
386        self.last_objf = self.objf;
387    }
388}
pmcore/algorithms/npag.rs

pmcore/algorithms/
npag.rs
