pmcore/algorithms/nonparametric/
npmap.rs

1use crate::{
2    algorithms::{NonParametricRunner, Status, StopReason},
3    estimation::nonparametric::{
4        calculate_psi, CycleLog, NPCycle, NonParametricResult, Psi, Theta, Weights,
5    },
6};
7
8use anyhow::{Context, Result};
9use pharmsol::prelude::{
10    data::{AssayErrorModels, Data},
11    simulator::Equation,
12};
13
14use crate::estimation::nonparametric::ipm::burke;
15use serde::{Deserialize, Serialize};
16
17/// Configuration options for the non-parametric maximum a posteriori (NPMAP) algorithm
18#[derive(Debug, Clone, PartialEq, Default, Serialize, Deserialize)]
19pub struct NpmapConfig {}
20
21impl NpmapConfig {
22    pub fn new() -> Self {
23        Self::default()
24    }
25}
26
27/// Non-parametric maximum a posteriori (NPMAP) algorithm
28///
29/// This algorithm is a wrapper around the IPM algorithm that calculates the posterior probabilities of the support points
30/// given a prior distribution and the likelihood of the data.
31#[derive(Debug)]
32pub struct NPMAP<E: Equation + Send + 'static> {
33    equation: E,
34    psi: Psi,
35    theta: Theta,
36    w: Weights,
37    objf: f64,
38    cycle: usize,
39    status: Status,
40    data: Data,
41    cyclelog: CycleLog,
42    error_models: AssayErrorModels,
43    prior: Theta,
44}
45
46impl<E: Equation + Send + 'static> NPMAP<E> {
47    pub(crate) fn from_parts(
48        equation: E,
49        data: Data,
50        error_models: AssayErrorModels,
51        theta: Theta,
52        _config: NpmapConfig,
53    ) -> Result<Self> {
54        Ok(Self {
55            equation,
56            psi: Psi::new(),
57            theta: theta.clone(),
58            w: Weights::default(),
59            objf: f64::INFINITY,
60            cycle: 0,
61            status: Status::Continue,
62            data,
63            cyclelog: CycleLog::new(),
64            error_models,
65            prior: theta,
66        })
67    }
68}
69
70impl<E: Equation + Send + 'static> NonParametricRunner<E> for NPMAP<E> {
71    fn into_result(&self) -> Result<NonParametricResult<E>> {
72        NonParametricResult::new(
73            self.equation.clone(),
74            self.data.clone(),
75            self.error_models.clone(),
76            self.prior.clone(),
77            self.theta.clone(),
78            self.psi.clone(),
79            self.w.clone(),
80            self.objf,
81            self.cycle,
82            self.status.clone(),
83            self.cyclelog.clone(),
84        )
85    }
86
87    fn error_models(&self) -> &AssayErrorModels {
88        &self.error_models
89    }
90
91    fn equation(&self) -> &E {
92        &self.equation
93    }
94
95    fn data(&self) -> &Data {
96        &self.data
97    }
98
99    fn likelihood(&self) -> f64 {
100        self.objf
101    }
102
103    fn increment_cycle(&mut self) -> usize {
104        0
105    }
106
107    fn cycle(&self) -> usize {
108        0
109    }
110
111    fn set_theta(&mut self, theta: Theta) {
112        self.theta = theta;
113    }
114
115    fn theta(&self) -> &Theta {
116        &self.theta
117    }
118
119    fn psi(&self) -> &Psi {
120        &self.psi
121    }
122
123    fn set_status(&mut self, status: Status) {
124        self.status = status;
125    }
126
127    fn status(&self) -> &Status {
128        &self.status
129    }
130
131    fn evaluation(&mut self) -> Result<Status> {
132        self.status = Status::Stop(StopReason::Converged);
133        Ok(self.status.clone())
134    }
135
136    fn estimation(&mut self) -> Result<()> {
137        self.psi = calculate_psi(
138            &self.equation,
139            &self.data,
140            &self.theta,
141            &self.error_models,
142            false,
143        )?;
144        (self.w, self.objf) = burke(&self.psi).context("Error in IPM")?;
145        Ok(())
146    }
147
148    fn condensation(&mut self) -> Result<()> {
149        Ok(())
150    }
151
152    fn optimizations(&mut self) -> Result<()> {
153        Ok(())
154    }
155
156    fn expansion(&mut self) -> Result<()> {
157        Ok(())
158    }
159
160    fn log_cycle_state(&mut self) {
161        // NPMAP doesn't track last_objf, so we use 0.0 as the delta
162        let state = NPCycle::new(
163            self.cycle,
164            self.objf,
165            self.error_models.clone(),
166            self.theta.clone(),
167            self.w.clone(),
168            self.theta.nspp(),
169            0.0,
170            self.status.clone(),
171        );
172        self.cyclelog.push(state);
173    }
174
175    /// NPMAP is a single-pass reweighting: it evaluates the likelihood of the
176    /// fixed prior support points once, rather than iterating cycles.
177    fn fit(&mut self) -> Result<NonParametricResult<E>> {
178        self.estimation()?;
179        self.evaluation()?;
180        self.log_cycle_state();
181
182        self.into_result()
183    }
184}
pmcore/algorithms/nonparametric/npmap.rs

pmcore/algorithms/nonparametric/
npmap.rs
