bcore/
impl_unique.rs

1use crate::api::{ModelEstimator, PostProcessReader};
2use crate::datamodel::{f_get_probes, make_histogram, read_spatial_model_properties, Results};
3use crate::datamodel::{
4    get_n_export_real, read_avg_model_properties, read_model_mass, read_model_properties,
5    tallies::Tallies, vec_to_array_view2, vec_to_array_view3, Dim, Weight,
6};
7use crate::process::{
8    spatial_average_concentration, variance_concentration, Histogram,
9};
10use crate::{api::Estimator, api::Phase, error::ApiError};
11use ndarray::{s, Array1, Array2, ArrayView2, ArrayView3, Axis};
12
13/// The `PostProcess` struct handles post-processing of simulation results.
14///
15/// It contains the path to the results folder, the root directory, and the processed results.
16///
17/// The struct uses the `Results` struct internally, which is expected to contain time and other data.
18#[derive(Debug)]
19pub struct PostProcess {
20    results: Results, // The results of the simulation, which will be accessed for time and other data.
21}
22
23impl PostProcess {
24    /// Creates a new instance of `PostProcess`.
25    ///
26    /// # Arguments
27    /// * `folder` - The name of the folder containing the simulation results.
28    /// * `root` - Optional root directory. Defaults to "./results/" if not provided.
29    ///
30    /// # Returns
31    /// * `Result<Self, String>` - Returns the `PostProcess` instance or an error message if initialization fails.
32    pub fn new(folder: &str, root: Option<String>) -> Result<Self, ApiError> {
33        let _root = root.unwrap_or_else(|| "./results/".to_string());
34        let result_path = format!("{}/{}/{}.h5", _root, folder, folder);
35        let main = Results::new(&result_path, &_root, folder)?;
36        Ok(Self { results: main })
37    }
38}
39
40impl PostProcessReader for PostProcess {
41    fn time(&self) -> &[f64] {
42        self.results.main.time()
43    }
44
45    fn weight(&self) -> &Weight {
46        &self.results.main.weight
47    }
48
49    fn get_property_names(&self) -> Vec<String> {
50        self.results.property_name.clone()
51    }
52
53    fn get_spatial_average_mtr(&self, species: usize) -> Result<Array1<f64>, ApiError> {
54        let r = &self.results.main.records;
55        let nt = r.time.len();
56        let dim = &r.dim;
57
58        if let Some(mtr) = &r.mtr {
59            let mtr = vec_to_array_view3(mtr, dim, nt);
60
61            return match mtr.slice(s![.., .., species]).mean_axis(Axis(1)) {
62                Some(avg) => Ok(avg),
63                None => Err(ApiError::ShapeError),
64            };
65        }
66
67        Err(ApiError::RecordsError("mtr".to_string()))
68
69        // let mtr = vec_to_array_view3(self.results., &dim, nt);
70    }
71
72    fn v_liquid(&self) -> ArrayView2<'_, f64> {
73        let nt = self.results.main.records.time.len();
74        let dim = &self.results.main.records.dim;
75        vec_to_array_view2(&self.results.main.records.volume_liquid, nt, dim.0)
76    }
77
78    fn get_variance_concentration(
79        &self,
80        species: usize,
81        phase: Phase,
82    ) -> Result<Array1<f64>, ApiError> {
83        // Helper
84        fn process_phase(
85            concentration: &Vec<f64>,
86            volume: &Vec<f64>,
87            nt: usize,
88            dim: &Dim,
89            species: usize,
90        ) -> Array1<f64> {
91            let c: ndarray::ArrayBase<ndarray::ViewRepr<&f64>, ndarray::Dim<[usize; 3]>> =
92                vec_to_array_view3(concentration, dim, nt);
93            let vol = vec_to_array_view2(volume, nt, dim.0);
94            let c_slice = &c.slice(s![.., .., species]);
95            
96            variance_concentration(c_slice, &vol)
97        }
98
99        let records = &self.results.main.records;
100        let nt = records.time.len();
101        let dim = &records.dim;
102
103        match phase {
104            Phase::Gas => {
105                if let (Some(c), Some(v)) = (&records.concentration_gas, &records.volume_gas) {
106                    return Ok(process_phase(c, v, nt, dim, species));
107                }
108
109                panic!("Gas is not present");
110            }
111            Phase::Liquid => Ok(process_phase(
112                &records.concentration_liquid,
113                &records.volume_liquid,
114                nt,
115                dim,
116                species,
117            )),
118        }
119    }
120
121    /// Returns an `ArrayView1<f64>` representing the time data from the simulation results.
122    ///
123    /// # Returns
124    /// * An `ArrayView1<f64>` containing the time data.
125    fn time_array(&self) -> Array1<f64> {
126        let time = self.results.main.time();
127        Array1::<f64>::from_vec(time.to_vec()) // TODO: safe unwrap
128    }
129
130    /// Retrieves the maximum number of export events specifically related to biological dumps from the simulation results.
131    /// This value represents the real number of biological export events, which may not always align with the total number of exports.
132    ///
133    /// # Returns
134    /// * `usize` - The maximum number of biological export events, or `0` if no events are found or if an error occurs.
135    fn get_max_n_export_bio(&self) -> usize {
136        get_n_export_real(self.results.get_files()).unwrap_or(0)
137    }
138
139    /// Returns the total number of export events from the simulation results.
140    /// This count includes all export actions, regardless of whether they are biological or not.
141    ///
142    /// # Returns
143    /// * `usize` - The total number of export events.
144    fn n_export(&self) -> usize {
145        self.results.main.records.time.len()
146    }
147
148    fn get_concentrations(&self, phase: Phase) -> ArrayView3<f64> {
149        let records = &self.results.main.records;
150        let nt = records.time.len();
151        let dim = &records.dim;
152
153        match phase {
154            Phase::Gas => {
155                if let (Some(c), Some(_v)) = (&records.concentration_gas, &records.volume_gas) {
156                    return vec_to_array_view3(c, dim, nt);
157                }
158
159                panic!("Gas is not present");
160            }
161            Phase::Liquid => vec_to_array_view3(&records.concentration_liquid, dim, nt),
162        }
163    }
164
165    fn get_spatial_average_concentration(&self, species: usize, phase: Phase) -> Array1<f64> {
166        // Helper
167        fn process_phase(
168            concentration: &Vec<f64>,
169            volume: &Vec<f64>,
170            nt: usize,
171            dim: &Dim,
172            species: usize,
173        ) -> Array1<f64> {
174            let cl = vec_to_array_view3(concentration, dim, nt);
175            let vol = vec_to_array_view2(volume, nt, dim.0);
176            let res = spatial_average_concentration(&cl.slice(s![.., .., species]), &vol);
177
178            res
179        }
180
181        let records = &self.results.main.records;
182        let nt = records.time.len();
183        let dim = &records.dim;
184
185        match phase {
186            Phase::Gas => {
187                if let (Some(c), Some(v)) = (&records.concentration_gas, &records.volume_gas) {
188                    return process_phase(c, v, nt, dim, species);
189                }
190
191                panic!("Gas is not present");
192            }
193            Phase::Liquid => process_phase(
194                &records.concentration_liquid,
195                &records.volume_liquid,
196                nt,
197                dim,
198                species,
199            ),
200        }
201    }
202
203    fn get_spatial_average_biomass_concentration(&self) -> Result<Array1<f64>, ApiError> {
204        let num_dimensions = self.results.main.records.dim.0;
205        let nt = self.results.main.records.time.len();
206        let volume =
207            vec_to_array_view2(&self.results.main.records.volume_liquid, nt, num_dimensions);
208        let vtot = volume.sum_axis(Axis(1));
209
210        let mut biomass_matrix = Array1::zeros(nt);
211
212        for i in 0..nt {
213            let m = self.get_properties("mass", i)?;
214            biomass_matrix[i] = m.sum() / vtot[i];
215        }
216
217        Ok(self.results.main.initial.initial_weight * biomass_matrix)
218    }
219
220    fn get_time_average_concentration(
221        &self,
222        species: usize,
223        position: usize,
224        phase: Phase,
225    ) -> Result<Array1<f64>, ApiError> {
226        let r = &self.results.main.records;
227        let nt = r.time.len();
228        let dim = &r.dim;
229
230        let callback = |c: &Vec<f64>| {
231            let cl = vec_to_array_view3(c, dim, nt);
232            cl.slice(s![.., .., species]).mean_axis(Axis(0)).unwrap()
233        };
234
235        match phase {
236            Phase::Liquid => Ok(callback(&r.concentration_liquid)),
237            Phase::Gas => {
238                if let Some(c) = &r.concentration_gas {
239                    return Ok(callback(c));
240                }
241
242                Err(ApiError::RecordsError("Gas".to_string()))
243            }
244        }
245    }
246
247    //Actually compute spatial average property cannot be use to follow distribution
248    fn get_spatial_average_property(&self, key: &str) -> Result<Array2<f64>, ApiError> {
249        let nt: usize = self.results.main.records.time.len(); // Number of time steps
250        let num_dimensions = self.results.main.records.dim.0; // Dimensionality
251
252        // Initialize the biomass matrix
253        let mut biomass_matrix = Array2::zeros((nt, num_dimensions));
254
255        if !self.results.property_name.iter().any(|x| x == key) {
256            return Err(ApiError::KeyError(key.to_string()));
257        }
258
259        read_spatial_model_properties(key, self.results.get_files(), &mut biomass_matrix, nt)?;
260
261        // Calculate biomass concentration
262        biomass_matrix /= self.get_number_particle();
263
264        Ok(biomass_matrix)
265    }
266
267    fn get_probes(&self) -> Result<Array1<f64>, ApiError> {
268        f_get_probes(&self.results.files)
269    }
270
271    /// Calculates the biomass concentration over time for the simulation.
272    ///
273    /// # Returns
274    /// * `Result<Array2<f64>, String>` - A 2D array with biomass concentrations or an error message.
275    fn get_biomass_concentration(&self) -> Result<Array2<f64>, ApiError> {
276        let nt: usize = self.results.main.records.time.len(); // Number of time steps
277        let num_dimensions = self.results.main.records.dim.0; // Dimensionality
278
279        // Initialize the biomass matrix
280        let mut biomass_matrix = Array2::zeros((nt, num_dimensions));
281
282        if !self.results.property_name.iter().any(|x| x == "mass") {
283            return Err(ApiError::KeyError("mass".to_string()));
284        }
285
286        // Attempt to read model mass
287        read_model_mass(self.results.get_files(), &mut biomass_matrix, nt)?;
288
289        // Convert volume to an array view
290        let volume =
291            vec_to_array_view2(&self.results.main.records.volume_liquid, nt, num_dimensions);
292
293        // Calculate biomass concentration
294        biomass_matrix = self.results.main.initial.initial_weight * (biomass_matrix / volume);
295
296        Ok(biomass_matrix)
297    }
298
299    /// Calculates the total growth in number
300    ///
301    /// # Returns
302    /// * `Array1<f64>` - 1D array containing the summed growth numbers.
303    fn get_growth_in_number(&self) -> Array1<f64> {
304        self.results.total_particle_repetition.sum_axis(Axis(1))
305    }
306
307    /// Retrieves the 2D array of particle numbers.
308    ///
309    /// # Returns
310    /// * `&Array2<f64>` - Reference to the 2D array containing particle numbers.
311    fn get_number_particle(&self) -> &Array2<f64> {
312        &self.results.total_particle_repetition
313    }
314
315    /// Fetches specific properties of the model at a given export index.
316    ///
317    /// # Arguments
318    /// * `key` - The property key to fetch.
319    /// * `i_export` - The index of the export to retrieve.
320    ///
321    /// # Returns
322    /// * `Result<Array1<f64>, String>` - A 1D array of property values or an error message.
323    fn get_properties(&self, key: &str, i_export: usize) -> Result<Array1<f64>, ApiError> {
324        if i_export >= self.results.main.records.time.len() {
325            return Err(ApiError::OutOfRange(
326                i_export,
327                self.results.main.records.time.len(),
328            ));
329        }
330
331        if !self.results.property_name.iter().any(|x| x == key) {
332            return Err(ApiError::KeyError(key.to_string()));
333        }
334
335        match read_model_properties(key, self.results.get_files(), i_export) {
336            Ok(res) => Ok(res),
337            // Err(e) => Err(format!("Failed to read model properties: {:?}", e)),
338            Err(e) => Err(ApiError::Io(e)),
339        }
340    }
341
342    /// Calculates the population mean over time for a given property key.
343    ///
344    /// # Arguments
345    /// * `key` - The property key for which to calculate the mean.
346    ///
347    /// # Returns
348    /// * `Result<Array1<f64>, String>` - A 1D array of mean values over time or an error message.
349    fn get_time_population_mean(&self, key: &str) -> Result<Array1<f64>, ApiError> {
350        if !self.results.property_name.iter().any(|x| x == key) {
351            return Err(ApiError::KeyError(key.to_string()));
352        }
353
354        match read_avg_model_properties(key, self.results.get_files(), self.n_export()) {
355            Ok(res) => Ok(res),
356            Err(e) => Err(ApiError::Io(e)),
357        }
358    }
359
360    fn get_histogram_array(
361        &self,
362        n_bins: usize,
363        i_export: usize,
364        key: &str,
365    ) -> Result<(Array1<f64>, Array1<f64>), ApiError> {
366        let (b, c) = self.get_histogram(n_bins, i_export, key)?;
367        let b = Array1::from_vec(b);
368        let c = Array1::from_vec(c);
369        Ok((b, c))
370    }
371
372    fn get_histogram(
373        &self,
374        n_bins: usize,
375        i_export: usize,
376        key: &str,
377    ) -> Result<(Vec<f64>, Vec<f64>), ApiError> {
378        if i_export > self.n_export() {
379            return Err(ApiError::OutOfRange(i_export, self.n_export()));
380        }
381        if !self.results.property_name.iter().any(|x| x == key) {
382            return Err(ApiError::KeyError(key.to_string()));
383        }
384        // let np = n_bins;//*self.results.total_particle_repetition.sum_axis(Axis(1)).last().unwrap() as usize;
385        let mut hist = Histogram::new(n_bins);
386
387        if let Err(e) = make_histogram(self.results.get_files(), i_export, key, &mut hist) {
388            return Err(ApiError::Io(e));
389        }
390
391        let b = hist.get_bins().to_vec();
392        let c = hist.get_counts().to_vec();
393        Ok((b, c))
394    }
395
396    fn get_population_mean(&self, key: &str, i_export: usize) -> Result<f64, ApiError> {
397        // Check if the index is out of range
398        if i_export >= self.results.main.records.time.len() {
399            return Err(ApiError::OutOfRange(
400                i_export,
401                self.results.main.records.time.len(),
402            ));
403        }
404        if !self.results.property_name.iter().any(|x| x == key) {
405            return Err(ApiError::KeyError(key.to_string()));
406        }
407
408        match read_model_properties(key, self.results.get_files(), i_export) {
409            Ok(res) => res
410                .mean()
411                .ok_or(ApiError::Default("get_population_mean".to_string())),
412            Err(e) => Err(ApiError::Io(e)),
413        }
414    }
415    fn tallies(&self) -> Option<&Tallies> {
416        self.results.main.records.tallies.as_ref()
417    }
418}
419
420impl ModelEstimator for PostProcess {
421    fn mu_direct(&self) -> Result<Array1<f64>, ApiError> {
422        match self.weight() {
423            Weight::Single(_) => {
424                let nt = self.results.main.records.time.len();
425                let time = self.time();
426
427                let mu_functor = |i: usize, j: usize, im: usize| -> Result<f64, ApiError> {
428                    let mass_i = self.get_properties("mass", i)?;
429                    let mass_mi = self.get_properties("mass", j)?;
430                    let mass_tot_i = self.get_properties("mass", im)?.sum();
431                    let dm: f64 = mass_i.sum() - mass_mi.sum();
432                    let dt = time[i] - time[j];
433                    Ok(dm / dt / mass_tot_i)
434                };
435
436                let mut mu = Array1::zeros(nt);
437
438                mu[0] = mu_functor(1, 0, 0)?; //Forward
439
440                for i in 1..nt - 1 {
441                    mu[i] = mu_functor(i + 1, i - 1, i)?; //Center
442                }
443                mu[nt - 1] = mu_functor(nt - 1, nt - 2, nt - 1)?; //Backward
444
445                Ok(mu)
446            }
447            Weight::Multiple(_) => todo!("Mu with different weight"),
448        }
449    }
450
451    fn estimate(&self, etype: Estimator, key: &str, i_export: usize) -> Result<f64, ApiError> {
452        crate::process::estimate(etype, self.weight(), &self.get_properties(key, i_export)?)
453    }
454
455    fn estimate_time(&self, etype: Estimator, key: &str) -> Result<Array1<f64>, ApiError> {
456        let nt = self.results.main.records.time.len();
457        let mut estimator = Array1::<f64>::zeros(nt);
458        for i in 0..nt {
459            estimator[i] = self.estimate(etype, key, i)?;
460        }
461        Ok(estimator)
462    }
463}
bcore/impl_unique.rs

bcore/
impl_unique.rs
