use crate::changepoint::detector::Cost;
use crate::changepoint::signal::Signal;
use crate::error::{ForecastError, Result};
#[derive(Debug, Default, Clone)]
pub struct CostPoisson {
n: usize,
d: usize,
cumsum: Vec<f64>,
cumsum_log: Vec<f64>,
}
impl CostPoisson {
pub fn new() -> Self {
Self::default()
}
}
impl Cost for CostPoisson {
fn fit(&mut self, signal: &Signal) -> Result<()> {
self.n = signal.n();
self.d = signal.d();
let stride = self.d;
self.cumsum = vec![0.0; (self.n + 1) * stride];
self.cumsum_log = vec![0.0; (self.n + 1) * stride];
for i in 0..self.n {
let row = signal.row(i);
for (j, &v) in row.iter().enumerate() {
let prev = i * stride + j;
let next = (i + 1) * stride + j;
self.cumsum[next] = self.cumsum[prev] + v;
let log_term = if v >= 1.0 { v * v.ln() - v } else { 0.0 };
self.cumsum_log[next] = self.cumsum_log[prev] + log_term;
}
}
Ok(())
}
fn error(&self, start: usize, end: usize) -> Result<f64> {
if end <= start || end > self.n {
return Err(ForecastError::InvalidParameter(format!(
"CostPoisson: invalid segment [{}, {}) for n = {}",
start, end, self.n
)));
}
let stride = self.d;
let len = (end - start) as f64;
let mut total = 0.0;
for j in 0..self.d {
let sum = self.cumsum[end * stride + j] - self.cumsum[start * stride + j];
let mean = sum / len;
if mean <= 1e-12 {
continue;
}
let log_sum_terms =
self.cumsum_log[end * stride + j] - self.cumsum_log[start * stride + j];
total += -2.0 * (sum * mean.ln() - len * mean - log_sum_terms);
}
Ok(total)
}
fn min_size(&self) -> usize {
1
}
fn name(&self) -> &str {
"poisson"
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn count_series_returns_finite_cost() {
let values: Vec<f64> = (1..=10).map(|i| i as f64).collect();
let s = Signal::univariate(&values);
let mut c = CostPoisson::new();
c.fit(&s).unwrap();
let e = c.error(0, 10).unwrap();
assert!(e.is_finite());
}
#[test]
fn zero_segment_zero_cost() {
let values = vec![0.0; 10];
let s = Signal::univariate(&values);
let mut c = CostPoisson::new();
c.fit(&s).unwrap();
let e = c.error(0, 10).unwrap();
assert_eq!(e, 0.0);
}
}