mod support;
use gspx::conv::{ChebyConvolver, DynamicConvolver, LowpassOptions, StaticConvolver};
use gspx::fit::{ChebyFitOptions, ChebyModel};
use gspx::kernel::{ChebyKernel, VfKernel, impulse};
use ndarray::{Array1, Array2, Axis, Zip, array};
use num_complex::Complex64;
use sprs::{CsMat, TriMat};
#[test]
fn static_filters_shape_and_scalar_scale() {
let l = support::path_laplacian(80);
let x = impulse(&l, 10, 1).expect("impulse should build");
let mut conv = StaticConvolver::new(l).expect("static convolver should initialize");
let scales = [0.1, 1.0, 10.0];
let lp = conv
.lowpass(x.view(), &scales, 1)
.expect("lowpass should run");
let bp = conv
.bandpass(x.view(), &scales, 1)
.expect("bandpass should run");
let hp = conv
.highpass(x.view(), &scales)
.expect("highpass should run");
assert_eq!(lp.len(), scales.len());
assert_eq!(bp.len(), scales.len());
assert_eq!(hp.len(), scales.len());
let y = conv
.lowpass_one(x.view(), 1.0, 1)
.expect("scalar lowpass should run");
assert_eq!(y.ndim(), 2);
assert_eq!(y.ncols(), 1);
}
#[test]
fn static_filters_match_spectral_definitions() {
let l = support::path_laplacian(18);
let x = support::random_signal(18, 2);
let scale = 2.75;
let order = 3;
let q = 1.0 / scale;
let mut conv = StaticConvolver::new(l.clone()).expect("static convolver should initialize");
let lp = conv
.lowpass_one(x.view(), scale, order)
.expect("lowpass should run");
let bp = conv
.bandpass_one(x.view(), scale, order)
.expect("bandpass should run");
let hp = conv
.highpass_one(x.view(), scale)
.expect("highpass should run");
let lp_ref =
apply_spectral_filter_dense(&l, &x, |lam| (1.0 / (scale * lam + 1.0)).powi(order as i32));
let bp_ref = apply_spectral_filter_dense(&l, &x, |lam| {
((4.0 * q * lam) / (lam + q).powi(2)).powi(order as i32)
});
let hp_ref = apply_spectral_filter_dense(&l, &x, |lam| (scale * lam) / (scale * lam + 1.0));
let err_lp = (&lp - &lp_ref)
.mapv(|v| v.abs())
.fold(0.0_f64, |a, b| a.max(*b));
let err_bp = (&bp - &bp_ref)
.mapv(|v| v.abs())
.fold(0.0_f64, |a, b| a.max(*b));
let err_hp = (&hp - &hp_ref)
.mapv(|v| v.abs())
.fold(0.0_f64, |a, b| a.max(*b));
assert!(err_lp < 1e-9);
assert!(err_bp < 1e-8);
assert!(err_hp < 1e-9);
}
#[test]
fn static_complex_signal_handling() {
let l = support::path_laplacian(60);
let x = impulse(&l, 10, 1).expect("impulse should build");
let xc = x.mapv(|v| Complex64::new(v, v));
let mut conv = StaticConvolver::new(l).expect("static convolver should initialize");
let yc = conv
.lowpass_one(xc.mapv(|z| z.re).view(), 1.0, 1)
.expect("lowpass on real part should run");
let y = conv
.lowpass_one(x.view(), 1.0, 1)
.expect("lowpass on real signal should run");
assert_eq!(yc.shape(), y.shape());
}
#[test]
fn vfkernel_convolution_and_direct_term() {
let l = support::path_laplacian(50);
let x = impulse(&l, 10, 1).expect("impulse should build");
let mut conv = StaticConvolver::new(l).expect("static convolver should initialize");
let k = VfKernel {
poles: array![1.0],
residues: array![[1.0]],
direct: array![5.0],
};
let y = conv
.convolve(x.view(), &k)
.expect("vf convolution should run");
let lp = conv
.lowpass_one(x.view(), 1.0, 1)
.expect("lowpass should run");
let ycol = y.index_axis(Axis(2), 0);
let expected = &lp + &(x.to_owned() * 5.0);
let err = (&ycol - &expected).mapv(|v| v.abs()).sum();
assert!(err < 1e-7);
}
#[test]
fn dynamic_topology_update_changes_response() {
let l = support::path_laplacian(80);
let x = impulse(&l, 10, 1).expect("impulse should build");
let mut conv = DynamicConvolver::new(l).expect("dynamic convolver should initialize");
let y0 = conv
.lowpass_one(x.view(), 1.0, 1)
.expect("baseline lowpass should run");
let ok = conv
.add_branch(10, 20, 1.0)
.expect("add_branch should succeed");
assert!(ok);
let y1 = conv
.lowpass_one(x.view(), 1.0, 1)
.expect("post-update lowpass should run");
let diff = (&y0 - &y1).mapv(|v| v.abs()).sum();
assert!(diff > 1e-8);
let oob = conv
.add_branch(1000, 2000, 1.0)
.expect("out-of-bounds branch add should return Ok(false)");
assert!(!oob);
assert!(conv.add_branch(1, 2, -1.0).is_err());
}
#[test]
fn cheby_convolution_basics() {
let l = support::path_laplacian(30);
let mut conv = ChebyConvolver::new(l).expect("cheby convolver should initialize");
let i_sig = Array2::<f64>::eye(30);
let k = ChebyKernel {
coefficients: array![[3.0]],
spectrum_bound: 4.0,
min_lambda: 0.0,
};
let y = conv
.convolve(i_sig.view(), &k)
.expect("cheby convolution should run");
let y0 = y.index_axis_move(Axis(2), 0);
let err = (&y0 - &(i_sig * 3.0)).mapv(|v| v.abs()).sum();
assert!(err < 1e-8);
}
#[test]
fn bandpass_order_equivalence() {
let l = support::path_laplacian(24);
let mut conv = StaticConvolver::new(l).expect("static convolver should initialize");
let x = support::random_signal(24, 2);
let s = 0.8;
let y_direct = conv
.bandpass_one(x.view(), s, 3)
.expect("direct order-3 bandpass should run");
let mut y_iter = x.clone();
for _ in 0..3 {
y_iter = conv
.bandpass_one(y_iter.view(), s, 1)
.expect("iterated order-1 bandpass should run");
}
let max_abs = (&y_direct - &y_iter)
.mapv(|v| v.abs())
.fold(0.0_f64, |a, b| a.max(*b));
assert!(max_abs < 1e-10);
}
#[test]
fn bandpass_order_zero_returns_input() {
let l = support::path_laplacian(20);
let x = support::random_signal(20, 3);
let mut static_conv = StaticConvolver::new(l.clone()).expect("static convolver should init");
let y_static = static_conv
.bandpass_one(x.view(), 1.0, 0)
.expect("order=0 static bandpass should run");
let err_static = (&y_static - &x).mapv(|v| v.abs()).sum();
assert!(err_static < 1e-12);
let mut dyn_conv = DynamicConvolver::new(l).expect("dynamic convolver should init");
let y_dynamic = dyn_conv
.bandpass_one(x.view(), 1.0, 0)
.expect("order=0 dynamic bandpass should run");
let err_dynamic = (&y_dynamic - &x).mapv(|v| v.abs()).sum();
assert!(err_dynamic < 1e-12);
}
#[test]
fn static_vf_convolution_matches_dense_reference() {
let l = support::path_laplacian(12);
let b = support::random_signal(12, 3);
let k = VfKernel {
poles: array![0.5, 2.0],
residues: array![[1.5, -0.5], [0.25, 0.75]],
direct: array![0.1, -0.2],
};
let mut conv = StaticConvolver::new(l.clone()).expect("static convolver should initialize");
let y = conv
.convolve(b.view(), &k)
.expect("vf convolution should run");
let mut expected = ndarray::Array3::<f64>::zeros((b.nrows(), b.ncols(), k.residues.ncols()));
for pole_idx in 0..k.poles.len() {
let z = solve_shifted_dense(&l, k.poles[pole_idx], &b);
for d in 0..k.residues.ncols() {
expected
.slice_mut(ndarray::s![.., .., d])
.scaled_add(k.residues[[pole_idx, d]], &z.view());
}
}
for d in 0..k.direct.len() {
expected
.slice_mut(ndarray::s![.., .., d])
.scaled_add(k.direct[d], &b.view());
}
let max_abs = (&y - &expected)
.mapv(|v| v.abs())
.fold(0.0_f64, |a, b| a.max(*b));
assert!(max_abs < 1e-10);
}
#[test]
fn dynamic_matches_static_rebuilt_laplacian_after_updates() {
let l0 = support::path_laplacian(40);
let x = support::random_signal(40, 2);
let updates = [(5usize, 30usize, 0.7f64), (10usize, 20usize, 0.25f64)];
let mut dyn_conv =
DynamicConvolver::new(l0.clone()).expect("dynamic convolver should initialize");
for &(i, j, w) in &updates {
assert!(
dyn_conv
.add_branch(i, j, w)
.expect("dynamic update should succeed")
);
}
let mut l_updated = l0.clone();
for &(i, j, w) in &updates {
l_updated = add_branch_to_laplacian(&l_updated, i, j, w);
}
let mut static_conv =
StaticConvolver::new(l_updated).expect("reference static convolver should initialize");
let y_dyn_lp = dyn_conv
.lowpass_one(x.view(), 1.2, 2)
.expect("dynamic lowpass should run");
let y_ref_lp = static_conv
.lowpass_one(x.view(), 1.2, 2)
.expect("reference static lowpass should run");
let y_dyn_bp = dyn_conv
.bandpass_one(x.view(), 0.9, 2)
.expect("dynamic bandpass should run");
let y_ref_bp = static_conv
.bandpass_one(x.view(), 0.9, 2)
.expect("reference static bandpass should run");
let err_lp = (&y_dyn_lp - &y_ref_lp)
.mapv(|v| v.abs())
.fold(0.0_f64, |a, b| a.max(*b));
let err_bp = (&y_dyn_bp - &y_ref_bp)
.mapv(|v| v.abs())
.fold(0.0_f64, |a, b| a.max(*b));
assert!(err_lp < 1e-8);
assert!(err_bp < 1e-8);
}
#[test]
fn cheby_multi_equivalence_and_high_order_stability() {
let l = support::path_laplacian(32);
let mut conv = ChebyConvolver::new(l).expect("cheby convolver should initialize");
let b = Array2::<f64>::eye(32);
let k1 = ChebyKernel {
coefficients: array![[1.0], [0.25], [-0.1], [0.05]],
spectrum_bound: 4.0,
min_lambda: 0.0,
};
let k2 = ChebyKernel {
coefficients: array![[0.8], [-0.2], [0.1], [0.02]],
spectrum_bound: 4.0,
min_lambda: 0.0,
};
let ys = conv
.convolve_multi(b.view(), &[k1.clone(), k2.clone()])
.expect("multi-kernel convolution should run");
assert_eq!(ys.len(), 2);
let y1 = conv
.convolve(b.view(), &k1)
.expect("single kernel 1 convolution should run");
let y2 = conv
.convolve(b.view(), &k2)
.expect("single kernel 2 convolution should run");
let err1 = (&ys[0] - &y1).mapv(|v| v.abs()).sum();
let err2 = (&ys[1] - &y2).mapv(|v| v.abs()).sum();
assert!(err1 < 1e-10);
assert!(err2 < 1e-10);
let empty = conv
.convolve_multi(b.view(), &[])
.expect("empty multi-kernel convolution should run");
assert!(empty.is_empty());
let model = ChebyModel::fit(
|x: &Array1<f64>| x.mapv(|v| (-v).exp()).insert_axis(Axis(1)),
60,
4.0,
Default::default(),
)
.expect("high-order cheby fit should succeed");
let k_stable = ChebyKernel {
coefficients: model.coefficients,
spectrum_bound: model.spectrum_bound,
min_lambda: model.min_lambda,
};
let y_stable = conv
.convolve(b.view(), &k_stable)
.expect("high-order cheby convolution should run");
assert!(y_stable.iter().all(|v| v.is_finite()));
}
#[test]
fn dynamic_vfkernel_direct_term_matches_lowpass_form() {
let l = support::path_laplacian(36);
let x = impulse(&l, 8, 1).expect("impulse should build");
let k = VfKernel {
poles: array![1.0],
residues: array![[1.0]],
direct: array![2.0],
};
let mut conv = DynamicConvolver::new(l).expect("dynamic convolver should initialize");
let y0 = conv
.convolve(x.view(), &k)
.expect("dynamic vf convolution should run");
let lp0 = conv
.lowpass_one(x.view(), 1.0, 1)
.expect("dynamic lowpass should run");
let expected0 = &lp0 + &(x.to_owned() * 2.0);
let err0 = (&y0.index_axis(Axis(2), 0) - &expected0)
.mapv(|v| v.abs())
.sum();
assert!(err0 < 1e-8);
assert!(
conv.add_branch(5, 24, 0.9)
.expect("dynamic branch update should succeed")
);
let y1 = conv
.convolve(x.view(), &k)
.expect("dynamic vf convolution after update should run");
let lp1 = conv
.lowpass_one(x.view(), 1.0, 1)
.expect("dynamic lowpass after update should run");
let expected1 = &lp1 + &(x.to_owned() * 2.0);
let err1 = (&y1.index_axis(Axis(2), 0) - &expected1)
.mapv(|v| v.abs())
.sum();
assert!(err1 < 1e-8);
}
#[test]
fn dynamic_constructor_with_poles_and_kernel_prefactors() {
let l = support::path_laplacian(18);
let x = support::random_signal(18, 1);
let mut conv_poles = DynamicConvolver::with_poles(l.clone(), &[0.5, 1.0, 2.0])
.expect("dynamic constructor with poles should succeed");
let y_poles = conv_poles
.lowpass_one(x.view(), 1.0, 1)
.expect("lowpass should run");
assert!(y_poles.iter().all(|v| v.is_finite()));
let k = VfKernel {
poles: array![0.5, 2.0],
residues: array![[1.0], [0.5]],
direct: array![0.0],
};
let mut conv_kernel = DynamicConvolver::with_kernel(l, &k)
.expect("dynamic constructor with kernel should succeed");
let y_kernel = conv_kernel
.convolve(x.view(), &k)
.expect("dynamic kernel convolution should run");
assert!(y_kernel.iter().all(|v| v.is_finite()));
}
#[test]
fn static_vf_complex_convolution_is_linear() {
let l = support::path_laplacian(30);
let x = impulse(&l, 6, 1).expect("impulse should build");
let xc = x.mapv(|v| Complex64::new(v, -0.4 * v));
let k = VfKernel {
poles: array![1.0, 2.0],
residues: array![[1.0], [0.3]],
direct: array![0.5],
};
let mut conv = StaticConvolver::new(l).expect("static convolver should initialize");
let y_complex = conv
.convolve_complex(xc.view(), &k)
.expect("complex vf convolution should run");
let y_real = conv
.convolve(xc.mapv(|z| z.re).view(), &k)
.expect("real-part vf convolution should run");
let y_imag = conv
.convolve(xc.mapv(|z| z.im).view(), &k)
.expect("imag-part vf convolution should run");
let mut expected = ndarray::Array3::<Complex64>::zeros(y_complex.raw_dim());
Zip::from(&mut expected)
.and(&y_real)
.and(&y_imag)
.for_each(|o, &r, &i| *o = Complex64::new(r, i));
let max_abs = (&y_complex - &expected)
.mapv(|z| z.norm())
.fold(0.0_f64, |a, b| a.max(*b));
assert!(max_abs < 1e-10);
}
#[test]
fn static_analytic_complex_filters_are_linear() {
let l = support::path_laplacian(24);
let x = impulse(&l, 7, 2).expect("impulse should build");
let xc = x.mapv(|v| Complex64::new(v, 0.5 * v));
let mut conv = StaticConvolver::new(l).expect("static convolver should init");
let lp_complex = conv
.lowpass_complex_one(xc.view(), 1.0, 2)
.expect("complex lowpass should run");
let lp_real = conv
.lowpass_one(xc.mapv(|z| z.re).view(), 1.0, 2)
.expect("real lowpass should run");
let lp_imag = conv
.lowpass_one(xc.mapv(|z| z.im).view(), 1.0, 2)
.expect("imag lowpass should run");
let lp_expected =
lp_real.mapv(|v| Complex64::new(v, 0.0)) + lp_imag.mapv(|v| Complex64::new(0.0, v));
let lp_err = (&lp_complex - &lp_expected).mapv(|z| z.norm()).sum();
assert!(lp_err < 1e-10);
let bp_complex = conv
.bandpass_complex_one(xc.view(), 1.0, 1)
.expect("complex bandpass should run");
let bp_real = conv
.bandpass_one(xc.mapv(|z| z.re).view(), 1.0, 1)
.expect("real bandpass should run");
let bp_imag = conv
.bandpass_one(xc.mapv(|z| z.im).view(), 1.0, 1)
.expect("imag bandpass should run");
let bp_expected =
bp_real.mapv(|v| Complex64::new(v, 0.0)) + bp_imag.mapv(|v| Complex64::new(0.0, v));
let bp_err = (&bp_complex - &bp_expected).mapv(|z| z.norm()).sum();
assert!(bp_err < 1e-10);
}
#[test]
fn static_lowpass_with_options_matches_baseline() {
let l = support::path_laplacian(30);
let x = impulse(&l, 5, 1).expect("impulse should build");
let mut conv = StaticConvolver::new(l).expect("static convolver should init");
let y_ref = conv
.lowpass(x.view(), &[1.0, 2.0], 1)
.expect("baseline lowpass should run");
let y_opt = conv
.lowpass_with_options(
x.view(),
&[1.0, 2.0],
LowpassOptions::with_refactor(false),
1,
)
.expect("lowpass with options should run");
assert_eq!(y_ref.len(), y_opt.len());
for (a, b) in y_ref.iter().zip(y_opt.iter()) {
let err = (a - b).mapv(|v| v.abs()).sum();
assert!(err < 1e-10);
}
}
#[test]
fn cheby_complex_convolution_is_linear() {
let l = support::path_laplacian(24);
let mut conv = ChebyConvolver::new(l).expect("cheby convolver should initialize");
let x = Array2::<f64>::eye(24);
let xc = x.mapv(|v| Complex64::new(v, 0.25 * v));
let k = ChebyKernel {
coefficients: array![[1.0], [0.4], [-0.15], [0.05]],
spectrum_bound: 4.0,
min_lambda: 0.0,
};
let y_complex = conv
.convolve_complex(xc.view(), &k)
.expect("complex cheby convolution should run");
let y_real = conv
.convolve(xc.mapv(|z| z.re).view(), &k)
.expect("real cheby convolution should run");
let y_imag = conv
.convolve(xc.mapv(|z| z.im).view(), &k)
.expect("imag cheby convolution should run");
let mut expected = ndarray::Array3::<Complex64>::zeros(y_complex.raw_dim());
Zip::from(&mut expected)
.and(&y_real)
.and(&y_imag)
.for_each(|o, &r, &i| *o = Complex64::new(r, i));
let max_abs = (&y_complex - &expected)
.mapv(|z| z.norm())
.fold(0.0_f64, |a, b| a.max(*b));
assert!(max_abs < 1e-10);
}
#[test]
fn cheby_fitted_analytic_filters_match_dense_reference() {
let l = support::path_laplacian(18);
let x = support::random_signal(18, 1);
let eig = nalgebra::linalg::SymmetricEigen::new(dense_from_sparse(&l));
let spectrum_bound = eig
.eigenvalues
.iter()
.fold(0.0_f64, |acc, value| acc.max(*value));
let scale = 1.0;
let lowpass_order = 2usize;
let bandpass_order = 2usize;
let cheby_degree = 32usize;
let lowpass_kernel = ChebyModel::fit(
|lambda: &Array1<f64>| {
gspx::functions::lowpass(lambda, scale)
.mapv(|v| v.powi(lowpass_order as i32))
.insert_axis(Axis(1))
},
cheby_degree,
spectrum_bound,
ChebyFitOptions {
rtol: 0.0,
..ChebyFitOptions::default()
},
)
.expect("cheby lowpass fit should succeed");
let bandpass_kernel = ChebyModel::fit(
|lambda: &Array1<f64>| {
gspx::functions::bandpass(lambda, scale, bandpass_order).insert_axis(Axis(1))
},
cheby_degree,
spectrum_bound,
ChebyFitOptions {
rtol: 0.0,
..ChebyFitOptions::default()
},
)
.expect("cheby bandpass fit should succeed");
let highpass_kernel = ChebyModel::fit(
|lambda: &Array1<f64>| gspx::functions::highpass(lambda, scale).insert_axis(Axis(1)),
cheby_degree,
spectrum_bound,
ChebyFitOptions {
rtol: 0.0,
..ChebyFitOptions::default()
},
)
.expect("cheby highpass fit should succeed");
let kernels = [
(
"lowpass",
ChebyKernel {
coefficients: lowpass_kernel.coefficients,
spectrum_bound: lowpass_kernel.spectrum_bound,
min_lambda: lowpass_kernel.min_lambda,
},
apply_spectral_filter_dense(&l, &x, |lam| {
gspx::functions::lowpass(&array![lam], scale)[0].powi(lowpass_order as i32)
}),
),
(
"bandpass",
ChebyKernel {
coefficients: bandpass_kernel.coefficients,
spectrum_bound: bandpass_kernel.spectrum_bound,
min_lambda: bandpass_kernel.min_lambda,
},
apply_spectral_filter_dense(&l, &x, |lam| {
gspx::functions::bandpass(&array![lam], scale, bandpass_order)[0]
}),
),
(
"highpass",
ChebyKernel {
coefficients: highpass_kernel.coefficients,
spectrum_bound: highpass_kernel.spectrum_bound,
min_lambda: highpass_kernel.min_lambda,
},
apply_spectral_filter_dense(&l, &x, |lam| {
gspx::functions::highpass(&array![lam], scale)[0]
}),
),
];
let mut conv = ChebyConvolver::new(l).expect("cheby convolver should initialize");
for (label, kernel, expected) in kernels {
let y = conv
.convolve_1d(x.column(0), &kernel)
.expect("cheby convolution should run");
let y_col = y.column(0).to_owned().insert_axis(Axis(1));
let max_abs = (&y_col - &expected)
.mapv(|v| v.abs())
.fold(0.0_f64, |acc, value| acc.max(*value));
assert!(max_abs < 1e-4, "{label} max_abs={max_abs}");
}
}
#[test]
fn cheby_lstsq_fitted_filters_match_dense_reference() {
let l = support::path_laplacian(18);
let x = support::random_signal(18, 1);
let eig = nalgebra::linalg::SymmetricEigen::new(dense_from_sparse(&l));
let spectrum_bound = eig
.eigenvalues
.iter()
.fold(0.0_f64, |acc, value| acc.max(*value));
let scale = 1.0;
let bandpass_order = 2usize;
for sampling in [
gspx::fit::Sampling::Linear,
gspx::fit::Sampling::Logarithmic,
] {
let kernel_model = ChebyModel::fit(
|lambda: &Array1<f64>| {
gspx::functions::bandpass(lambda, scale, bandpass_order).insert_axis(Axis(1))
},
32,
spectrum_bound,
ChebyFitOptions {
sampling,
rtol: 0.0,
..ChebyFitOptions::default()
},
)
.expect("cheby bandpass fit should succeed");
let kernel = ChebyKernel {
coefficients: kernel_model.coefficients,
spectrum_bound: kernel_model.spectrum_bound,
min_lambda: kernel_model.min_lambda,
};
let expected = apply_spectral_filter_dense(&l, &x, |lam| {
gspx::functions::bandpass(&array![lam], scale, bandpass_order)[0]
});
let mut conv = ChebyConvolver::new(l.clone()).expect("cheby convolver should initialize");
let y = conv
.convolve_1d(x.column(0), &kernel)
.expect("cheby convolution should run");
let y_col = y.column(0).to_owned().insert_axis(Axis(1));
let max_abs = (&y_col - &expected)
.mapv(|v| v.abs())
.fold(0.0_f64, |acc, value| acc.max(*value));
assert!(
max_abs < 1e-4,
"bandpass {sampling:?} convolution error too large: {max_abs}"
);
}
}
fn dense_from_sparse(a: &CsMat<f64>) -> nalgebra::DMatrix<f64> {
let n = a.rows();
let mut data = vec![0.0_f64; n * n];
for (col, vec) in a.outer_iterator().enumerate() {
for (row, val) in vec.iter() {
data[row * n + col] = *val;
}
}
nalgebra::DMatrix::<f64>::from_row_slice(n, n, &data)
}
fn apply_spectral_filter_dense<F>(l: &CsMat<f64>, b: &Array2<f64>, f: F) -> Array2<f64>
where
F: Fn(f64) -> f64,
{
let dense = dense_from_sparse(l);
let eig = nalgebra::linalg::SymmetricEigen::new(dense);
let v = eig.eigenvectors;
let vt = v.transpose();
let lam = eig.eigenvalues;
let mut out = Array2::<f64>::zeros((b.nrows(), b.ncols()));
for col in 0..b.ncols() {
let x = nalgebra::DVector::<f64>::from_column_slice(&b.column(col).to_vec());
let mut coeff = &vt * x;
for i in 0..coeff.len() {
coeff[i] *= f(lam[i]);
}
let y = &v * coeff;
for row in 0..b.nrows() {
out[[row, col]] = y[row];
}
}
out
}
fn solve_shifted_dense(l: &CsMat<f64>, q: f64, b: &Array2<f64>) -> Array2<f64> {
let n = l.rows();
let mut a = dense_from_sparse(l);
for i in 0..n {
a[(i, i)] += q;
}
let lu = a.lu();
let mut out = Array2::<f64>::zeros((n, b.ncols()));
for c in 0..b.ncols() {
let rhs = nalgebra::DVector::<f64>::from_column_slice(&b.column(c).to_vec());
let x = lu
.solve(&rhs)
.expect("dense shifted solve should be nonsingular for q > 0");
for r in 0..n {
out[[r, c]] = x[r];
}
}
out
}
fn add_branch_to_laplacian(l: &CsMat<f64>, i: usize, j: usize, w: f64) -> CsMat<f64> {
let n = l.rows();
let mut tri = TriMat::<f64>::with_capacity((n, n), l.nnz() + 4);
for (col, vec) in l.outer_iterator().enumerate() {
for (row, val) in vec.iter() {
tri.add_triplet(row, col, *val);
}
}
tri.add_triplet(i, i, w);
tri.add_triplet(j, j, w);
tri.add_triplet(i, j, -w);
tri.add_triplet(j, i, -w);
tri.to_csc()
}