pmcore/algorithms/nonparametric/
error_optim.rs

1//! Error-model factor optimization used by non-parametric algorithms.
2
3use anyhow::Result;
4use pharmsol::prelude::{
5    data::{AssayErrorModels, Data},
6    simulator::Equation,
7};
8use serde::{Deserialize, Serialize};
9
10use crate::estimation::nonparametric::{calculate_psi, ipm::burke, Psi, Theta, Weights};
11
12/// Configuration for the error-model factor (gamma/lambda) optimization.
13#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
14pub struct ErrorOptimConfig {
15    /// Initial and reset perturbation step applied to each error-model factor.
16    pub step: f64,
17    /// Minimum perturbation step; once the adaptive step reaches this value it is reset to `step`.
18    pub min_step: f64,
19    /// Factor by which the step grows after a successful improvement.
20    pub growth: f64,
21    /// Factor by which the step shrinks after each optimization pass.
22    pub shrink: f64,
23}
24
25impl Default for ErrorOptimConfig {
26    fn default() -> Self {
27        Self {
28            step: 0.1,
29            min_step: 0.01,
30            growth: 4.0,
31            shrink: 0.5,
32        }
33    }
34}
35
36impl ErrorOptimConfig {
37    pub fn new() -> Self {
38        Self::default()
39    }
40
41    pub fn step(mut self, step: f64) -> Self {
42        self.step = step;
43        self
44    }
45
46    pub fn min_step(mut self, min_step: f64) -> Self {
47        self.min_step = min_step;
48        self
49    }
50
51    pub fn growth(mut self, growth: f64) -> Self {
52        self.growth = growth;
53        self
54    }
55
56    pub fn shrink(mut self, shrink: f64) -> Self {
57        self.shrink = shrink;
58        self
59    }
60}
61
62/// Perform one pass of the error-model factor (gamma/lambda) optimization.
63///
64/// For each optimizable error model, the factor is perturbed up and down,
65/// `psi` is recomputed for each perturbation, and the IPM ([`burke`]) is run.
66/// The direction with the highest objective function that also improves on the
67/// current `objf` is adopted, updating `error_models`, `objf`, `lambda`, and
68/// `psi` in place. The per-output adaptive step in `gamma_delta` is grown on
69/// improvement and shrunk every pass, resetting once it falls to `config.min_step`.
70///
71/// A failure of the IPM for one perturbation direction is treated as a warning
72/// (that direction is skipped); if both directions fail the factor is left
73/// unchanged for this pass.
74#[allow(clippy::too_many_arguments)]
75pub(crate) fn optimize_error_models<E: Equation + Send + 'static>(
76    equation: &E,
77    data: &Data,
78    theta: &Theta,
79    error_models: &mut AssayErrorModels,
80    gamma_delta: &mut [f64],
81    objf: &mut f64,
82    lambda: &mut Weights,
83    psi: &mut Psi,
84    config: &ErrorOptimConfig,
85) -> Result<()> {
86    error_models
87        .clone()
88        .iter_mut()
89        .filter_map(|(outeq, em)| {
90            if em.optimize() {
91                Some((outeq, em))
92            } else {
93                None
94            }
95        })
96        .try_for_each(|(outeq, em)| -> Result<()> {
97            // Optimize the best value of gamma/lambda
98
99            let gamma_up = em.factor()? * (1.0 + gamma_delta[outeq]);
100            let gamma_down = em.factor()? / (1.0 + gamma_delta[outeq]);
101
102            let mut error_model_up = error_models.clone();
103            error_model_up.set_factor(outeq, gamma_up)?;
104
105            let mut error_model_down = error_models.clone();
106            error_model_down.set_factor(outeq, gamma_down)?;
107
108            let psi_up = calculate_psi(equation, data, theta, &error_model_up, false)?;
109            let psi_down = calculate_psi(equation, data, theta, &error_model_down, false)?;
110
111            // Treat errors in IPM, such as subjects with zero probability or non-finite
112            // likelihoods as warnings, and continue with the other direction. If both
113            // directions fail, we will not update the error model factor.
114            let up = match burke(&psi_up) {
115                Ok((lambda, objf)) => Some((lambda, objf)),
116                Err(err) => {
117                    tracing::warn!(
118                        "Error in IPM during optim (up) for outeq {}: {:?}",
119                        outeq,
120                        err
121                    );
122                    None
123                }
124            };
125            let down = match burke(&psi_down) {
126                Ok((lambda, objf)) => Some((lambda, objf)),
127                Err(err) => {
128                    tracing::warn!(
129                        "Error in IPM during optim (down) for outeq {}: {:?}",
130                        outeq,
131                        err
132                    );
133                    None
134                }
135            };
136
137            // Select the best improving candidate (if any) over the current
138            // objective. Among the two directions, the one with the higher
139            // objective function wins.
140            let mut best: Option<(f64, Weights, Psi, f64)> = None;
141            if let Some((lambda_up, objf_up)) = up {
142                if objf_up > *objf {
143                    best = Some((objf_up, lambda_up, psi_up, gamma_up));
144                }
145            }
146            if let Some((lambda_down, objf_down)) = down {
147                let threshold = best.as_ref().map_or(*objf, |(o, ..)| *o);
148                if objf_down > threshold {
149                    best = Some((objf_down, lambda_down, psi_down, gamma_down));
150                }
151            }
152            if let Some((new_objf, new_lambda, new_psi, gamma)) = best {
153                error_models.set_factor(outeq, gamma)?;
154                *objf = new_objf;
155                gamma_delta[outeq] *= config.growth;
156                *lambda = new_lambda;
157                *psi = new_psi;
158            }
159            gamma_delta[outeq] *= config.shrink;
160            if gamma_delta[outeq] <= config.min_step {
161                gamma_delta[outeq] = config.step;
162            }
163            Ok(())
164        })?;
165
166    Ok(())
167}
pmcore/algorithms/nonparametric/error_optim.rs

pmcore/algorithms/nonparametric/
error_optim.rs
