pmcore/estimation/nonparametric/
cycles.rs

1use std::{fs::File, path::Path};
2
3use anyhow::Result;
4use csv::WriterBuilder;
5use pharmsol::{AssayErrorModel, AssayErrorModels};
6use serde::Serialize;
7
8use crate::{
9    algorithms::Status,
10    estimation::nonparametric::{median, theta::Theta, weights::Weights},
11};
12
13#[derive(Debug, Clone, Serialize)]
14pub struct NPCycle {
15    cycle: usize,
16    objf: f64,
17    error_models: AssayErrorModels,
18    theta: Theta,
19    weights: Weights,
20    nspp: usize,
21    delta_objf: f64,
22    status: Status,
23}
24
25impl NPCycle {
26    pub fn new(
27        cycle: usize,
28        objf: f64,
29        error_models: AssayErrorModels,
30        theta: Theta,
31        weights: Weights,
32        nspp: usize,
33        delta_objf: f64,
34        status: Status,
35    ) -> Self {
36        Self {
37            cycle,
38            objf,
39            error_models,
40            theta,
41            weights,
42            nspp,
43            delta_objf,
44            status,
45        }
46    }
47
48    pub fn cycle(&self) -> usize {
49        self.cycle
50    }
51    pub fn objf(&self) -> f64 {
52        self.objf
53    }
54    pub fn error_models(&self) -> &AssayErrorModels {
55        &self.error_models
56    }
57    pub fn theta(&self) -> &Theta {
58        &self.theta
59    }
60    pub fn weights(&self) -> &Weights {
61        &self.weights
62    }
63    pub fn nspp(&self) -> usize {
64        self.nspp
65    }
66    pub fn delta_objf(&self) -> f64 {
67        self.delta_objf
68    }
69    pub fn status(&self) -> &Status {
70        &self.status
71    }
72
73    pub fn placeholder() -> Self {
74        Self {
75            cycle: 0,
76            objf: 0.0,
77            error_models: AssayErrorModels::default(),
78            theta: Theta::new(),
79            weights: Weights::default(),
80            nspp: 0,
81            delta_objf: 0.0,
82            status: Status::Continue,
83        }
84    }
85}
86
87#[derive(Debug, Clone, Serialize)]
88pub struct CycleLog {
89    cycles: Vec<NPCycle>,
90}
91
92impl CycleLog {
93    pub fn new() -> Self {
94        Self { cycles: Vec::new() }
95    }
96
97    pub fn cycles(&self) -> &[NPCycle] {
98        &self.cycles
99    }
100
101    pub fn push(&mut self, cycle: NPCycle) {
102        self.cycles.push(cycle);
103    }
104
105    pub fn write(&self, path: &Path) -> Result<()> {
106        tracing::debug!("Writing cycles...");
107
108        super::create_parent_dir(path)?;
109        let file = File::create(path)?;
110        let mut writer = WriterBuilder::new().has_headers(true).from_writer(file);
111
112        writer.write_field("cycle")?;
113        writer.write_field("converged")?;
114        writer.write_field("status")?;
115        writer.write_field("neg2ll")?;
116        writer.write_field("nspp")?;
117        if let Some(first_cycle) = self.cycles.first() {
118            first_cycle.error_models.iter().try_for_each(
119                |(outeq, errmod): (usize, &AssayErrorModel)| -> Result<(), csv::Error> {
120                    match errmod {
121                        AssayErrorModel::Additive { .. } => {
122                            writer.write_field(format!("gamlam.{}", outeq))?;
123                        }
124                        AssayErrorModel::Proportional { .. } => {
125                            writer.write_field(format!("gamlam.{}", outeq))?;
126                        }
127                        AssayErrorModel::None => {}
128                    }
129                    Ok(())
130                },
131            )?;
132        }
133
134        let names = self
135            .cycles
136            .first()
137            .map(|cycle| cycle.theta.param_names())
138            .expect("No cycles");
139
140        for param_name in names {
141            writer.write_field(format!("{}.mean", param_name))?;
142            writer.write_field(format!("{}.median", param_name))?;
143            writer.write_field(format!("{}.sd", param_name))?;
144        }
145
146        writer.write_record(None::<&[u8]>)?;
147
148        for cycle in &self.cycles {
149            writer.write_field(format!("{}", cycle.cycle))?;
150            writer.write_field(format!("{}", cycle.status.converged()))?;
151            writer.write_field(format!("{}", cycle.status))?;
152            writer.write_field(format!("{}", cycle.objf))?;
153            writer
154                .write_field(format!("{}", cycle.theta.nspp()))
155                .unwrap();
156
157            cycle.error_models.iter().try_for_each(
158                |(_, errmod): (usize, &AssayErrorModel)| -> Result<()> {
159                    match errmod {
160                        AssayErrorModel::Additive { lambda: _, poly: _ } => {
161                            writer.write_field(format!("{:.5}", errmod.factor()?))?;
162                        }
163                        AssayErrorModel::Proportional { gamma: _, poly: _ } => {
164                            writer.write_field(format!("{:.5}", errmod.factor()?))?;
165                        }
166                        AssayErrorModel::None => {}
167                    }
168                    Ok(())
169                },
170            )?;
171
172            for param in cycle.theta.matrix().col_iter() {
173                let param_values: Vec<f64> = param.iter().cloned().collect();
174
175                let mean: f64 = param_values.iter().sum::<f64>() / param_values.len() as f64;
176                let median = median(&param_values);
177                let std = param_values.iter().map(|x| (x - mean).powi(2)).sum::<f64>()
178                    / (param_values.len() as f64 - 1.0);
179
180                writer.write_field(format!("{}", mean))?;
181                writer.write_field(format!("{}", median))?;
182                writer.write_field(format!("{}", std))?;
183            }
184            writer.write_record(None::<&[u8]>)?;
185        }
186        writer.flush()?;
187
188        tracing::debug!("Cycles written to {:?}", path);
189        Ok(())
190    }
191}
192
193impl Default for CycleLog {
194    fn default() -> Self {
195        Self::new()
196    }
197}
pmcore/estimation/nonparametric/cycles.rs

pmcore/estimation/nonparametric/
cycles.rs
