use crate::helpers::simpsons_weights;
use crate::matrix::FdMatrix;
use super::{cross_distance_matrix, self_distance_matrix};
fn shifted_l2_distance(x: &[f64], y: &[f64], weights: &[f64], shift: i32) -> Option<f64> {
let n = x.len();
let mut sum = 0.0;
let mut valid_points = 0;
for i in 0..n {
let j = i as i32 + shift;
if j >= 0 && (j as usize) < n {
let diff = x[i] - y[j as usize];
sum += weights[i] * diff * diff;
valid_points += 1;
}
}
if valid_points >= n / 2 {
Some(sum.sqrt())
} else {
None
}
}
fn hshift_distance(x: &[f64], y: &[f64], weights: &[f64], max_shift: usize) -> f64 {
if x.is_empty() || y.len() != x.len() || weights.len() != x.len() {
return f64::INFINITY;
}
let mut min_dist = f64::INFINITY;
for shift in -(max_shift as i32)..=(max_shift as i32) {
if let Some(dist) = shifted_l2_distance(x, y, weights, shift) {
if dist < min_dist {
min_dist = dist;
}
}
}
min_dist
}
pub fn hshift_self_1d(data: &FdMatrix, argvals: &[f64], max_shift: usize) -> FdMatrix {
let n = data.nrows();
let m = data.ncols();
if n == 0 || m == 0 || argvals.len() != m {
return FdMatrix::zeros(0, 0);
}
let weights = simpsons_weights(argvals);
let rows: Vec<Vec<f64>> = (0..n).map(|i| data.row(i)).collect();
self_distance_matrix(n, |i, j| {
hshift_distance(&rows[i], &rows[j], &weights, max_shift)
})
}
pub fn hshift_cross_1d(
data1: &FdMatrix,
data2: &FdMatrix,
argvals: &[f64],
max_shift: usize,
) -> FdMatrix {
let n1 = data1.nrows();
let n2 = data2.nrows();
let m = data1.ncols();
if n1 == 0 || n2 == 0 || m == 0 || argvals.len() != m || data2.ncols() != m {
return FdMatrix::zeros(0, 0);
}
let weights = simpsons_weights(argvals);
let rows1: Vec<Vec<f64>> = (0..n1).map(|i| data1.row(i)).collect();
let rows2: Vec<Vec<f64>> = (0..n2).map(|i| data2.row(i)).collect();
cross_distance_matrix(n1, n2, |i, j| {
hshift_distance(&rows1[i], &rows2[j], &weights, max_shift)
})
}