pmcore/estimation/
problem.rs

1use anyhow::{anyhow, Result};
2use pharmsol::{
3    AssayErrorModel, AssayErrorModels, Data, Equation, Event, ResidualErrorModel,
4    ResidualErrorModels,
5};
6use std::collections::{BTreeSet, HashSet};
7
8use crate::estimation::nonparametric::Theta;
9use crate::model::parameter_space::{BoundedParameter, ParameterSpace, UnboundedParameter};
10use crate::model::{EquationMetadataSource, Model, ModelBuilder};
11
12pub trait Framework {
13    type ErrorModels;
14    /// The prior that seeds the algorithm.
15    ///
16    /// For the non-parametric framework this is a [`Theta`] (a discrete prior
17    /// distribution that also carries the parameter space). For the parametric
18    /// framework it is the [`ParameterSpace`] of unbounded parameters.
19    type Prior;
20}
21
22#[derive(Debug, Clone, Copy)]
23pub struct Parametric;
24
25impl Framework for Parametric {
26    type ErrorModels = ResidualErrorModels;
27    type Prior = ParameterSpace<UnboundedParameter>;
28}
29
30#[derive(Debug, Clone, Copy)]
31pub struct NonParametric;
32
33impl Framework for NonParametric {
34    type ErrorModels = AssayErrorModels;
35    type Prior = Theta;
36}
37
38#[derive(Debug, Clone)]
39pub struct EstimationProblem<E: Equation, F: Framework> {
40    pub(crate) model: Model<E>,
41    pub(crate) data: Data,
42    pub(crate) error_models: F::ErrorModels,
43    /// The prior that seeds the algorithm.
44    ///
45    /// For the non-parametric framework this is the prior [`Theta`], which also
46    /// carries the parameter space. The parameter space is therefore not stored
47    /// separately.
48    pub(crate) prior: F::Prior,
49}
50
51impl<E: Equation + EquationMetadataSource> EstimationProblem<E, NonParametric> {
52    /// Creates a non-parametric estimation problem.
53    ///
54    /// The `prior` is a [`Theta`] holding the prior distribution (the initial
55    /// set of support points) together with the [`ParameterSpace`] it was built
56    /// from. The parameter space is taken directly from the prior, so there is
57    /// no separate parameter-declaration step.
58    pub fn nonparametric(
59        equation: E,
60        data: Data,
61        prior: Theta,
62        error_models: AssayErrorModels,
63    ) -> Result<Self> {
64        let model_builder = Model::builder(equation);
65
66        validate_nonparametric_parameters(&model_builder, prior.parameters())?;
67
68        let model = model_builder.build()?;
69
70        validate_nonparametric_error_models(&model, &data, &error_models)?;
71
72        Ok(EstimationProblem {
73            model,
74            data,
75            error_models,
76            prior,
77        })
78    }
79}
80
81impl<E: Equation> EstimationProblem<E, Parametric> {
82    /// Begins building a parametric estimation problem.
83    pub fn parametric(equation: E, data: Data) -> ParametricBuilder<E> {
84        ParametricBuilder {
85            model: Model::builder(equation),
86            data,
87            parameters: ParameterSpace::<UnboundedParameter>::new(),
88            error_models: Vec::new(),
89        }
90    }
91
92    /// Returns the parameter space defined for this problem.
93    pub fn parameters(&self) -> &ParameterSpace<UnboundedParameter> {
94        &self.prior
95    }
96}
97
98impl<E: Equation> EstimationProblem<E, NonParametric> {
99    /// Returns the parameter space carried by the prior [`Theta`].
100    pub fn parameters(&self) -> &ParameterSpace<BoundedParameter> {
101        self.prior.parameters()
102    }
103}
104
105pub struct ParametricBuilder<E: Equation> {
106    model: ModelBuilder<E>,
107    data: Data,
108    parameters: ParameterSpace<UnboundedParameter>,
109    error_models: Vec<(String, ResidualErrorModel)>,
110}
111
112impl<E: Equation> ParametricBuilder<E> {
113    pub fn parameter(mut self, parameter: impl Into<UnboundedParameter>) -> Self {
114        self.parameters.push(parameter.into());
115        self
116    }
117
118    pub fn parameters<P, I>(mut self, parameters: I) -> Self
119    where
120        P: Into<UnboundedParameter>,
121        I: IntoIterator<Item = P>,
122    {
123        for param in parameters {
124            self.parameters.push(param.into());
125        }
126        self
127    }
128
129    pub fn error_model(mut self, name: impl Into<String>, model: ResidualErrorModel) -> Self {
130        self.error_models.push((name.into(), model));
131        self
132    }
133}
134
135impl<E: Equation + EquationMetadataSource> ParametricBuilder<E> {
136    pub fn build(self) -> Result<EstimationProblem<E, Parametric>> {
137        validate_parametric_parameters(&self.model, &self.parameters)?;
138        validate_parametric_error_models(&self.model, &self.error_models)?;
139
140        let mut all_errors = ResidualErrorModels::new();
141        for (name, error_model) in self.error_models {
142            let outeq = self
143                .model
144                .output_index(&name)
145                .ok_or_else(|| anyhow!("unknown equation output label: {name}"))?;
146
147            all_errors = all_errors.add(outeq, error_model);
148        }
149
150        Ok(EstimationProblem {
151            model: self.model.build()?,
152            data: self.data,
153            error_models: all_errors,
154            prior: self.parameters,
155        })
156    }
157}
158
159fn validate_nonparametric_parameters<E: Equation + EquationMetadataSource>(
160    model: &ModelBuilder<E>,
161    parameters: &ParameterSpace<BoundedParameter>,
162) -> Result<()> {
163    if parameters.is_empty() {
164        anyhow::bail!("at least one parameter is required for non-parametric models");
165    }
166
167    for parameter in parameters.iter() {
168        if !parameter.lower.is_finite() || !parameter.upper.is_finite() {
169            anyhow::bail!(
170                "invalid bounds for parameter '{}': bounds must be finite numbers",
171                parameter.name
172            );
173        }
174
175        if parameter.lower >= parameter.upper {
176            anyhow::bail!(
177                "invalid bounds for parameter '{}': lower bound ({}) must be strictly less than upper bound ({})",
178                parameter.name,
179                parameter.lower,
180                parameter.upper
181            );
182        }
183    }
184
185    let names: Vec<String> = parameters
186        .iter()
187        .map(|parameter| parameter.name.clone())
188        .collect();
189    validate_parameter_declarations(model, &names)
190}
191
192fn validate_parametric_parameters<E: Equation + EquationMetadataSource>(
193    model: &ModelBuilder<E>,
194    parameters: &ParameterSpace<UnboundedParameter>,
195) -> Result<()> {
196    if parameters.is_empty() {
197        anyhow::bail!("at least one parameter is required for parametric models");
198    }
199
200    let names: Vec<String> = parameters
201        .iter()
202        .map(|parameter| parameter.name.clone())
203        .collect();
204    validate_parameter_declarations(model, &names)
205}
206
207fn validate_parameter_declarations<E: Equation + EquationMetadataSource>(
208    model: &ModelBuilder<E>,
209    provided_names: &[String],
210) -> Result<()> {
211    let mut seen: HashSet<&str> = HashSet::new();
212    let mut duplicates: Vec<String> = Vec::new();
213    for name in provided_names {
214        if !seen.insert(name.as_str()) {
215            duplicates.push(name.clone());
216        }
217    }
218
219    if !duplicates.is_empty() {
220        duplicates.sort();
221        duplicates.dedup();
222        anyhow::bail!(
223            "duplicate parameter declarations found: {}",
224            duplicates.join(", ")
225        );
226    }
227
228    let declared = model.parameter_names();
229
230    let unknown: Vec<String> = provided_names
231        .iter()
232        .filter(|name| model.parameter_index(name).is_none())
233        .cloned()
234        .collect();
235    if !unknown.is_empty() {
236        anyhow::bail!(
237            "unknown parameter name(s): {}. Valid parameters are: {}",
238            unknown.join(", "),
239            declared.join(", ")
240        );
241    }
242
243    let provided: HashSet<&str> = provided_names.iter().map(|name| name.as_str()).collect();
244    let missing: Vec<String> = declared
245        .iter()
246        .filter(|name| !provided.contains(name.as_str()))
247        .cloned()
248        .collect();
249
250    if !missing.is_empty() {
251        anyhow::bail!("missing parameter declaration(s): {}", missing.join(", "));
252    }
253
254    Ok(())
255}
256
257fn validate_nonparametric_error_models<E: Equation + EquationMetadataSource>(
258    model: &Model<E>,
259    data: &Data,
260    error_models: &AssayErrorModels,
261) -> Result<()> {
262    // Bind the (label-first) error models to the equation. This resolves and
263    // validates that every declared output label maps to a valid model output.
264    let bound = model
265        .equation
266        .bind_error_models(error_models)
267        .map_err(|e| anyhow!("invalid assay error model output(s): {e}"))?;
268
269    // Collect the set of model output indices that are actually observed in the
270    // data, resolving each observation's output label the same way the simulator
271    // does (exact name, then the `outeq_<N>` numeric alias).
272    let mut observed_outputs: BTreeSet<usize> = BTreeSet::new();
273    let mut unresolved_labels: BTreeSet<String> = BTreeSet::new();
274    for subject in data.subjects() {
275        for occasion in subject.occasions() {
276            for event in occasion.events() {
277                if let Event::Observation(obs) = event {
278                    let label = obs.outeq().to_string();
279                    match resolve_output_index(model, &label) {
280                        Some(outeq) => {
281                            observed_outputs.insert(outeq);
282                        }
283                        None => {
284                            unresolved_labels.insert(label);
285                        }
286                    }
287                }
288            }
289        }
290    }
291
292    if !unresolved_labels.is_empty() {
293        let labels: Vec<String> = unresolved_labels.into_iter().collect();
294        anyhow::bail!(
295            "the data references output label(s) that are not defined by the model: {}",
296            labels.join(", ")
297        );
298    }
299
300    if observed_outputs.is_empty() {
301        anyhow::bail!("the data contains no observations to fit");
302    }
303
304    // Every observed output must have a (non-`None`) assay error model.
305    for &outeq in &observed_outputs {
306        let has_model = matches!(
307            bound.error_model(outeq),
308            Ok(error_model) if *error_model != AssayErrorModel::None
309        );
310
311        if !has_model {
312            let label = model
313                .output_name(outeq)
314                .map(|name| name.to_string())
315                .unwrap_or_else(|| outeq.to_string());
316            anyhow::bail!(
317                "no assay error model defined for output '{}' (index {}), which is observed in the data",
318                label,
319                outeq
320            );
321        }
322    }
323
324    Ok(())
325}
326
327/// Resolves an observation output `label` to a model output index, mirroring the
328/// simulator: first by exact output name, then via the `outeq_<N>` numeric alias.
329fn resolve_output_index<E: Equation + EquationMetadataSource>(
330    model: &Model<E>,
331    label: &str,
332) -> Option<usize> {
333    model.output_index(label).or_else(|| {
334        if !label.is_empty() && label.bytes().all(|b| b.is_ascii_digit()) {
335            model.output_index(&format!("outeq_{label}"))
336        } else {
337            None
338        }
339    })
340}
341
342fn validate_parametric_error_models<E: Equation + EquationMetadataSource>(
343    model: &ModelBuilder<E>,
344    error_models: &[(String, ResidualErrorModel)],
345) -> Result<()> {
346    if error_models.is_empty() {
347        anyhow::bail!("at least one residual error model is required");
348    }
349
350    validate_error_model_labels(model, error_models.iter().map(|(name, _)| name.as_str()))
351}
352
353fn validate_error_model_labels<'a, E, I>(model: &ModelBuilder<E>, labels: I) -> Result<()>
354where
355    E: Equation + EquationMetadataSource,
356    I: IntoIterator<Item = &'a str>,
357{
358    let valid_outputs = model.output_names();
359    let mut seen_output_indexes: HashSet<usize> = HashSet::new();
360
361    for name in labels {
362        let outeq = model.output_index(name).ok_or_else(|| {
363            anyhow!(
364                "unknown equation output label: {}. Valid outputs are: {}",
365                name,
366                valid_outputs.join(", ")
367            )
368        })?;
369
370        if !seen_output_indexes.insert(outeq) {
371            anyhow::bail!(
372                "duplicate error model declaration for output '{}' (index {})",
373                name,
374                outeq
375            );
376        }
377    }
378
379    Ok(())
380}
pmcore/estimation/problem.rs

pmcore/estimation/
problem.rs
