cmtool/
sanitizer.rs

1// SPDX-License-Identifier: GPL-3.0-or-later
2
3use std::fmt::Write;
4
5///Symmetric mean absolute percentage error
6///
7/// Details : https://en.wikipedia.org/wiki/Symmetric_mean_absolute_percentage_error
8fn smape(a: f64, b: f64) -> f64 {
9    if a == 0. && b == 0. {
10        return 0.;
11    }
12    (a - b).abs() / (a.abs() + b.abs())
13}
14
15pub fn check_flows(
16    _grid: &dyn cmtool_core::grid::CompartmentMeshManip,
17    raw_flows: &cmtool_data::RawDataFlux,
18) -> Option<String> {
19    let mut mass_balance: Vec<(f64, f64)> = vec![(0.0, 0.0); raw_flows.header.n_zone as usize];
20
21    // let boundary = grid.get_boundary();
22
23    for flow in raw_flows.fluxes.iter() {
24        //out
25        mass_balance[flow.id_source as usize].0 += flow.flux_target_source;
26        //in
27        mass_balance[flow.id_source as usize].1 += flow.flux_source_target;
28
29        //in
30        mass_balance[flow.id_target as usize].0 += flow.flux_source_target;
31        //out
32        mass_balance[flow.id_target as usize].1 += flow.flux_target_source;
33    }
34
35    let mut zone_relative_errors = Vec::with_capacity(raw_flows.header.n_zone as usize);
36    let mut total_relative_error: f64 = 0.0;
37    let mut max_relative_error: f64 = 0.0;
38
39    let mut total_inflow = 0.0;
40    let mut total_outflow = 0.0;
41    let mut f = String::new();
42    let mut n_zone = 0;
43
44    writeln!(&mut f, "zone_id,int,out,relative_error_percent").unwrap();
45    for (i, (inflow, outflow)) in mass_balance.iter().enumerate() {
46        // let is_boundary = boundary.iter().find(|&&ci| ci == i);
47        // if is_boundary.is_some() {
48        //     continue;
49        // }
50        n_zone += 1;
51        let relative_error = smape(*inflow, *outflow);
52
53        zone_relative_errors.push(relative_error);
54        total_relative_error += relative_error;
55        max_relative_error = max_relative_error.max(relative_error);
56
57        total_inflow += inflow;
58        total_outflow += outflow;
59
60        // writeln!(
61        //     &mut f,
62        //     "{},{},{},{},{}",
63        //     i + 1,
64        //     inflow,
65        //     outflow,
66        //     relative_error * 100.0,
67        //     is_boundary.is_some()
68        // )
69        // .unwrap();
70        writeln!(
71            &mut f,
72            "{},{},{},{}",
73            i + 1,
74            inflow,
75            outflow,
76            relative_error * 100.0,
77        )
78        .unwrap();
79    }
80    writeln!(&mut f).unwrap();
81
82    let mean_relative_error = total_relative_error / n_zone as f64;
83
84    writeln!(&mut f, "metric,value").unwrap();
85    writeln!(&mut f, "total_inflow,{:.8}", total_inflow).unwrap();
86    writeln!(&mut f, "total_outflow,{:.8}", total_outflow).unwrap();
87    writeln!(
88        &mut f,
89        "global_net_error_percent,{:.8}",
90        smape(total_inflow, total_outflow) * 100.0
91    )
92    .unwrap();
93    writeln!(
94        &mut f,
95        "mean_relative_error_percent,{:.6}",
96        mean_relative_error * 100.0
97    )
98    .unwrap();
99    writeln!(
100        &mut f,
101        "max_relative_error_percent,{:.6}",
102        max_relative_error * 100.0
103    )
104    .unwrap();
105
106    Some(f)
107}
108
109pub fn divergence_free(raw_flows: &cmtool_data::RawDataFlux) -> cmtool_data::RawDataFlux {
110    use nalgebra::{DMatrix, DVector};
111    let n_zone = raw_flows.header.n_zone as usize;
112    let n_flux = raw_flows.fluxes.len() * 2;
113    let mut new_flows = raw_flows.clone();
114
115    let mut f_vec = DVector::from_element(n_flux, 0.0);
116    for (i, flow) in new_flows.fluxes.iter().enumerate() {
117        f_vec[i * 2] = flow.flux_source_target;
118        f_vec[i * 2 + 1] = flow.flux_target_source;
119    }
120
121    let mut a_mat = DMatrix::zeros(n_zone, n_flux);
122    for (k, flow) in new_flows.fluxes.iter().enumerate() {
123        let s = flow.id_source as usize;
124        let t = flow.id_target as usize;
125
126        a_mat[(s, k * 2)] = 1.0;
127        a_mat[(t, k * 2)] = -1.0;
128
129        a_mat[(t, k * 2 + 1)] = 1.0;
130        a_mat[(s, k * 2 + 1)] = -1.0;
131    }
132
133    let div = &a_mat * &f_vec;
134
135    let a_at = a_mat.transpose();
136    let aat_inv = match (a_mat.clone() * a_at.clone()).try_inverse() {
137        Some(inv) => inv,
138        None => panic!("Cannot invert A*A^T"),
139    };
140    let delta_f = a_at * (aat_inv * (-div));
141
142    for (i, flow) in new_flows.fluxes.iter_mut().enumerate() {
143        flow.flux_source_target += delta_f[i * 2];
144        flow.flux_target_source += delta_f[i * 2 + 1];
145
146        if flow.flux_source_target < 0.0 || flow.flux_target_source < 0.0 {
147            panic!("Negative flux encountered");
148        }
149    }
150
151    new_flows
152}
cmtool/sanitizer.rs

cmtool/
sanitizer.rs
