use crate::pod::problem::parameter::ParameterOrdering;
#[cfg(feature = "serde")]
use serde::{Deserialize, Serialize};
#[derive(Debug, Clone, PartialEq)]
#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
pub struct ParameterCovariance {
pub params: ParameterOrdering,
pub data: Vec<f64>,
}
#[derive(Debug, thiserror::Error)]
pub enum CovarianceError {
#[error("dimension mismatch: expected {expected}x{expected}, got {actual} entries")]
Dimension {
expected: usize,
actual: usize,
},
#[error("not symmetric: |C[{i},{j}] - C[{j},{i}]| = {delta} > tol {tol}")]
NotSymmetric {
i: usize,
j: usize,
delta: f64,
tol: f64,
},
#[error("not positive semi-definite: diagonal entry {i} = {value} < 0")]
NotPsdDiagonal {
i: usize,
value: f64,
},
}
impl ParameterCovariance {
pub fn n(&self) -> usize {
self.params.len()
}
pub fn validate_dim(&self) -> Result<(), CovarianceError> {
let n = self.n();
let expected = n * n;
if self.data.len() != expected {
return Err(CovarianceError::Dimension {
expected: n,
actual: self.data.len(),
});
}
Ok(())
}
pub fn validate_symmetric(&self, tol: f64) -> Result<(), CovarianceError> {
self.validate_dim()?;
let n = self.n();
for i in 0..n {
for j in (i + 1)..n {
let a = self.data[i * n + j];
let b = self.data[j * n + i];
let delta = (a - b).abs();
if delta > tol {
return Err(CovarianceError::NotSymmetric { i, j, delta, tol });
}
}
}
Ok(())
}
pub fn validate_diagonal_nonneg(&self) -> Result<(), CovarianceError> {
self.validate_dim()?;
let n = self.n();
for i in 0..n {
let v = self.data[i * n + i];
if v < 0.0 {
return Err(CovarianceError::NotPsdDiagonal { i, value: v });
}
}
Ok(())
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::pod::problem::parameter::{Parameter, ParameterKind};
fn ord(n: usize) -> ParameterOrdering {
ParameterOrdering {
params: (0..n)
.map(|_| Parameter {
kind: ParameterKind::DragScale,
initial_value: 1.0,
apriori_sigma: None,
})
.collect(),
}
}
#[test]
fn validates_symmetry() {
let cov = ParameterCovariance {
params: ord(2),
data: vec![1.0, 0.5, 0.5, 2.0],
};
cov.validate_symmetric(1e-12).unwrap();
cov.validate_diagonal_nonneg().unwrap();
}
#[test]
fn detects_asymmetry() {
let cov = ParameterCovariance {
params: ord(2),
data: vec![1.0, 0.5, 0.6, 2.0],
};
assert!(cov.validate_symmetric(1e-12).is_err());
}
#[test]
fn detects_negative_diagonal() {
let cov = ParameterCovariance {
params: ord(1),
data: vec![-1.0],
};
assert!(cov.validate_diagonal_nonneg().is_err());
}
}