#[cfg(feature = "gpu-wgpu")]
mod gpu_visual {
use either::Either;
use ndarray::Array1;
use stochastic_rs::stats::fd::FractalDim;
use stochastic_rs::stochastic::noise::fgn::Fgn;
use stochastic_rs::stochastic::process::fbm::Fbm;
use stochastic_rs::traits::ProcessExt;
use stochastic_rs::visualization::GridPlotter;
fn gpu_fgn_paths(h: f32, n: usize, m: usize) -> Vec<Vec<f64>> {
let fgn = Fgn::<f32>::new(h, n, Some(1.0));
match fgn.sample_gpu(m).expect("GPU Fgn sampling failed") {
Either::Left(p) => vec![p.iter().map(|&x| x as f64).collect()],
Either::Right(ps) => ps
.outer_iter()
.map(|r| r.iter().map(|&x| x as f64).collect())
.collect(),
}
}
fn cpu_fgn_paths(h: f64, n: usize, m: usize) -> Vec<Vec<f64>> {
let fgn = Fgn::<f64>::new(h, n, Some(1.0));
fgn.sample_par(m).into_iter().map(|p| p.to_vec()).collect()
}
fn gpu_fbm_paths(h: f32, n: usize, m: usize) -> Vec<Vec<f64>> {
let fbm = Fbm::<f32>::new(h, n, Some(1.0));
match fbm.sample_gpu(m).expect("GPU Fbm sampling failed") {
Either::Left(p) => vec![p.iter().map(|&x| x as f64).collect()],
Either::Right(ps) => ps
.outer_iter()
.map(|r| r.iter().map(|&x| x as f64).collect())
.collect(),
}
}
fn empirical_autocovariance(paths: &[Vec<f64>], max_lag: usize) -> Vec<f64> {
let all: Vec<f64> = paths.iter().flatten().copied().collect();
let mean = all.iter().sum::<f64>() / all.len() as f64;
(0..=max_lag)
.map(|lag| {
let mut s = 0.0;
let mut c = 0usize;
for p in paths {
for i in 0..(p.len() - lag) {
s += (p[i] - mean) * (p[i + lag] - mean);
c += 1;
}
}
s / c as f64
})
.collect()
}
fn theoretical_autocovariance(h: f64, n: usize, max_lag: usize) -> Vec<f64> {
let dt = 1.0 / n as f64;
let var = dt.powf(2.0 * h);
(0..=max_lag)
.map(|k| {
if k == 0 {
var
} else {
var
* 0.5
* (((k + 1) as f64).powf(2.0 * h) - 2.0 * (k as f64).powf(2.0 * h)
+ ((k - 1) as f64).powf(2.0 * h))
}
})
.collect()
}
#[test]
#[ignore = "visual: writes HTML, no assertions; run with --ignored"]
fn plot_cpu_vs_gpu_autocovariance() {
let n = 1024;
let m = 4096;
let max_lag = 20;
let hursts = [0.25_f64, 0.5, 0.72, 0.9];
let mut grid = GridPlotter::new()
.title("fGN autocovariance: theory vs CPU vs GPU")
.cols(2)
.line_width(2.0)
.show_legend(true);
for &h in &hursts {
let theory = theoretical_autocovariance(h, n, max_lag);
let cpu_acov = empirical_autocovariance(&cpu_fgn_paths(h, n, m), max_lag);
let gpu_acov = empirical_autocovariance(&gpu_fgn_paths(h as f32, n, m), max_lag);
eprintln!("H={h}: autocovariance (lag 0..{max_lag})");
eprintln!(" lag theory CPU GPU CPU/th GPU/th");
for k in 0..=max_lag.min(10) {
eprintln!(
" {k:<4} {:<15.8} {:<15.8} {:<15.8} {:<7.4} {:<7.4}",
theory[k],
cpu_acov[k],
gpu_acov[k],
cpu_acov[k] / theory[k],
gpu_acov[k] / theory[k],
);
}
grid = grid.register_paths(
vec![theory, cpu_acov, gpu_acov],
&format!("H={h} (blue=theory, red=CPU, green=GPU)"),
);
}
let plot = grid.plot();
plot.write_html("target/gpu_autocovariance.html");
eprintln!("\nWrote target/gpu_autocovariance.html");
}
#[test]
#[ignore = "visual: writes HTML, no assertions; run with --ignored"]
fn plot_gpu_fgn_trajectories() {
let n = 1024;
let traj = 8;
let hursts = [0.2_f32, 0.35, 0.5, 0.72, 0.85, 0.95];
let mut grid = GridPlotter::new()
.title("GPU fGN trajectories (CubeCL)")
.cols(3)
.line_width(1.0)
.show_legend(false);
for &h in &hursts {
let paths = gpu_fgn_paths(h, n, traj);
grid = grid.register_paths(paths, &format!("fGN H={h}"));
}
let plot = grid.plot();
plot.write_html("target/gpu_fgn_trajectories.html");
eprintln!("Wrote target/gpu_fgn_trajectories.html");
}
#[test]
#[ignore = "visual: writes HTML, no assertions; run with --ignored"]
fn plot_gpu_fbm_trajectories() {
let n = 1024;
let traj = 8;
let hursts = [0.2_f32, 0.35, 0.5, 0.72, 0.85, 0.95];
let mut grid = GridPlotter::new()
.title("GPU fBM trajectories (CubeCL)")
.cols(3)
.line_width(1.0)
.show_legend(false);
for &h in &hursts {
let paths = gpu_fbm_paths(h, n, traj);
grid = grid.register_paths(paths, &format!("fBM H={h}"));
}
let plot = grid.plot();
plot.write_html("target/gpu_fbm_trajectories.html");
eprintln!("Wrote target/gpu_fbm_trajectories.html");
}
#[test]
fn gpu_fbm_hurst_estimation() {
let n = 4096;
let m = 64;
let hursts = [0.2_f32, 0.35, 0.5, 0.72, 0.85];
eprintln!("\nGPU fBM Hurst estimation (n={n}, m={m}):");
eprintln!(
"{:<8} {:<12} {:<12} {:<12}",
"H_true", "H_vario", "H_higuchi", "avg_err"
);
for &h_true in &hursts {
let paths = gpu_fbm_paths(h_true, n, m);
let mut h_vario_sum = 0.0;
let mut h_higuchi_sum = 0.0;
for path in &paths {
let arr = Array1::from_vec(path.to_vec());
let fd = FractalDim::new(arr);
h_vario_sum += 2.0 - fd.variogram(Some(2.0)).expect("variogram on fGN path");
h_higuchi_sum += 2.0 - fd.higuchi_fd(32).expect("Higuchi on fGN path");
}
let h_vario = h_vario_sum / m as f64;
let h_higuchi = h_higuchi_sum / m as f64;
let h64 = h_true as f64;
let avg_err = ((h_vario - h64).abs() + (h_higuchi - h64).abs()) / 2.0;
eprintln!(
"{:<8.2} {:<12.4} {:<12.4} {:<12.4}",
h_true, h_vario, h_higuchi, avg_err
);
if h64 <= 0.75 {
assert!(
(h_vario - h64).abs() < 0.08,
"H={h_true}: variogram estimate {h_vario} too far"
);
}
assert!(
(h_higuchi - h64).abs() < 0.08,
"H={h_true}: higuchi estimate {h_higuchi} too far"
);
}
}
}