zilla-muf 0.1.1

use num_traits::Float;

/// Causal Toeplitz matrix-vector product — the convolutional view of an
/// SSM. Implicitly applies the lower-triangular matrix
///   M[i][j] = kernel[i - j]  for j <= i,  0 otherwise
/// to `x`, i.e. computes the causal convolution:
///   y[i] = sum_{j=0}^{i} kernel[i - j] * x[j]
///
/// `kernel` and `x` must be the same length n; cost is O(n^2). This is
/// the direct reference implementation — `fft_conv` computes the exact
/// same quantity in O(n log n) and is tested against this function.
pub fn toeplitz_matvec<T: Float>(kernel: &[T], x: &[T]) -> Vec<T> {
	let n = x.len();
	assert_eq!(kernel.len(), n, "kernel and x must be the same length");

	// For each output i, sum kernel[i - j] * x[j] over j <= i. The
	// `i - j` lag index is what makes this Toeplitz (entries are constant
	// along each diagonal); the `j <= i` bound is what makes it causal
	// (no future input can influence an earlier output).
	(0..n)
		.map(|i| {
			(0..=i).fold(T::zero(), |acc, j| acc + kernel[i - j] * x[j])
		})
		.collect()
}

#[cfg(test)]
mod tests {
	use super::*;

	/// Builds the dense n x n lower-triangular Toeplitz matrix directly
	/// and does a naive matvec — an oracle independent of the
	/// convolution-loop formulation used above.
	fn dense_reference(kernel: &[f64], x: &[f64]) -> Vec<f64> {
		let n = x.len();
		let mut y = vec![0.0; n];
		for i in 0..n {
			for j in 0..=i {
				y[i] += kernel[i - j] * x[j];
			}
		}
		y
	}

	#[test]
	fn matches_dense_reference() {
		let n = 8;
		let kernel: Vec<f64> = (0..n).map(|i| (i as f64 * 0.41).cos() * 0.8_f64.powi(i as i32)).collect();
		let x: Vec<f64> = (0..n).map(|i| (i as f64 * 0.23).sin()).collect();

		let expected = dense_reference(&kernel, &x);
		let actual = toeplitz_matvec(&kernel, &x);
		for (e, got) in expected.iter().zip(actual.iter()) {
			assert!((e - got).abs() < 1e-9, "expected {e}, got {got}");
		}
	}

	#[test]
	fn impulse_kernel_passes_signal_through() {
		// kernel = [1, 0, 0, ...] is the identity matrix -> y == x
		let n = 5;
		let mut kernel = vec![0.0; n];
		kernel[0] = 1.0;
		let x: Vec<f64> = (0..n).map(|i| i as f64 + 1.0).collect();
		assert_eq!(toeplitz_matvec(&kernel, &x), x);
	}
}