#[derive(Debug, Clone, Copy)]
pub struct DfaResult {
pub alpha: f64,
pub r_squared: f64,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum LawQuality {
Exact,
Strong,
Good,
Approx,
Abstain,
Insufficient,
}
#[derive(Debug, Clone, Copy)]
pub struct StructuralLaw {
pub hurst: f64,
pub dfa: DfaResult,
pub acr: DfaResult,
pub mean: f64,
pub std_dev: f64,
pub kurtosis: f64,
pub p99: f64,
pub max: f64,
pub n: usize,
pub quality: LawQuality,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum HealthVerdict {
Healthy,
Watch,
Warning,
Critical,
}
impl HealthVerdict {
pub fn from_shift(shift: f64) -> Self {
let s = if shift < 0.0 { -shift } else { shift };
if s < 0.03 {
HealthVerdict::Healthy
} else if s < 0.08 {
HealthVerdict::Watch
} else if s < 0.15 {
HealthVerdict::Warning
} else {
HealthVerdict::Critical
}
}
}
pub fn dfa(values: &[f64]) -> DfaResult {
let n = values.len();
if n < 64 {
return DfaResult { alpha: 0.5, r_squared: 0.0 };
}
let mean = values.iter().sum::<f64>() / n as f64;
let mut y = Vec::with_capacity(n);
let mut cum = 0.0;
for &v in values {
cum += v - mean;
y.push(cum);
}
const BOXES: [usize; 8] = [16, 24, 36, 54, 81, 121, 181, 271];
let mut log_s = [0.0f64; 8];
let mut log_f = [0.0f64; 8];
let mut pts = 0usize;
for &s in &BOXES {
if s > n / 4 { break; }
let num_segs = n / s;
if num_segs == 0 { continue; }
let mut f2_sum = 0.0;
for seg in 0..num_segs {
let start = seg * s;
let (mut sx, mut sy, mut sxy, mut sx2) = (0.0, 0.0, 0.0, 0.0);
for i in 0..s {
let xi = i as f64;
sx += xi;
sy += y[start + i];
sxy += xi * y[start + i];
sx2 += xi * xi;
}
let k = s as f64;
let det = k * sx2 - sx * sx;
if det.abs() < 1e-15 { continue; }
let a0 = (sx2 * sy - sx * sxy) / det;
let a1 = (k * sxy - sx * sy) / det;
let mut resid = 0.0;
for i in 0..s {
let d = y[start + i] - (a0 + a1 * i as f64);
resid += d * d;
}
f2_sum += resid / k;
}
let f = (f2_sum / num_segs as f64).sqrt();
if f > 0.0 {
log_s[pts] = (s as f64).ln();
log_f[pts] = f.ln();
pts += 1;
}
}
if pts < 3 {
return DfaResult { alpha: 0.5, r_squared: 0.0 };
}
linreg(&log_s[..pts], &log_f[..pts])
}
pub fn acr(values: &[f64]) -> DfaResult {
let n = values.len();
if n < 20 {
return DfaResult { alpha: 0.0, r_squared: 0.0 };
}
let mean = values.iter().sum::<f64>() / n as f64;
let var: f64 = values.iter().map(|&x| (x - mean) * (x - mean)).sum();
if var < 1e-15 {
return DfaResult { alpha: 0.0, r_squared: 0.0 };
}
const LAGS: [usize; 10] = [1, 2, 3, 5, 8, 13, 21, 34, 55, 89];
let mut log_lag = [0.0f64; 10];
let mut log_r = [0.0f64; 10];
let mut pts = 0usize;
for &lag in &LAGS {
if lag >= n / 2 { break; }
let mut num = 0.0;
for i in 0..n - lag {
num += (values[i] - mean) * (values[i + lag] - mean);
}
let r = num / var;
if r > 0.001 {
log_lag[pts] = (lag as f64).ln();
log_r[pts] = r.ln();
pts += 1;
}
}
if pts < 3 {
return DfaResult { alpha: 0.0, r_squared: 0.0 };
}
linreg(&log_lag[..pts], &log_r[..pts])
}
pub fn analyze(values: &[f64]) -> StructuralLaw {
let n = values.len();
if n < 20 {
return StructuralLaw {
hurst: 0.5, dfa: DfaResult { alpha: 0.5, r_squared: 0.0 },
acr: DfaResult { alpha: 0.0, r_squared: 0.0 },
mean: 0.0, std_dev: 0.0, kurtosis: 0.0, p99: 0.0, max: 0.0,
n, quality: LawQuality::Insufficient,
};
}
let mean = values.iter().sum::<f64>() / n as f64;
let var: f64 = values.iter().map(|&x| (x - mean) * (x - mean)).sum::<f64>() / n as f64;
let std_dev = var.sqrt();
let sd = if std_dev > 1e-15 { std_dev } else { 1e-15 };
let kurtosis = values.iter().map(|&v| {
let z = (v - mean) / sd;
z * z * z * z
}).sum::<f64>() / n as f64;
let max = values.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
let mut sorted = values.to_vec();
sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(core::cmp::Ordering::Equal));
let p99 = sorted[((n as f64 * 0.99) as usize).min(n - 1)];
let dfa_result = dfa(values);
let acr_result = acr(values);
let hurst = clamp(1.0 + acr_result.alpha / 2.0, 0.0, 1.0);
let best_r2 = if dfa_result.r_squared > acr_result.r_squared { dfa_result.r_squared } else { acr_result.r_squared };
let quality = if best_r2 > 0.95 { LawQuality::Exact }
else if best_r2 > 0.85 { LawQuality::Strong }
else if best_r2 > 0.7 { LawQuality::Good }
else if best_r2 > 0.3 { LawQuality::Approx }
else { LawQuality::Abstain };
StructuralLaw { hurst, dfa: dfa_result, acr: acr_result, mean, std_dev, kurtosis, p99, max, n, quality }
}
pub fn health_check(current: &StructuralLaw, baseline_alpha: f64) -> HealthVerdict {
HealthVerdict::from_shift(current.dfa.alpha - baseline_alpha)
}
fn clamp(v: f64, lo: f64, hi: f64) -> f64 {
if v < lo { lo } else if v > hi { hi } else { v }
}
fn linreg(x: &[f64], y: &[f64]) -> DfaResult {
let k = x.len() as f64;
let (mut sx, mut sy, mut sxy, mut sx2) = (0.0, 0.0, 0.0, 0.0);
for i in 0..x.len() {
sx += x[i]; sy += y[i]; sxy += x[i] * y[i]; sx2 += x[i] * x[i];
}
let slope = (k * sxy - sx * sy) / (k * sx2 - sx * sx);
let ic = (sy - slope * sx) / k;
let ym = sy / k;
let mut sst = 0.0;
let mut ssr = 0.0;
for i in 0..x.len() {
sst += (y[i] - ym) * (y[i] - ym);
ssr += (y[i] - slope * x[i] - ic) * (y[i] - slope * x[i] - ic);
}
let r2 = 1.0 - ssr / if sst > 1e-15 { sst } else { 1e-15 };
DfaResult { alpha: slope, r_squared: r2 }
}
#[cfg(test)]
mod tests {
use super::*;
fn white_noise(n: usize, seed: u64) -> Vec<f64> {
let mut state = seed;
(0..n).map(|_| {
state = state.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
(state >> 33) as f64 / (1u64 << 31) as f64 - 0.5
}).collect()
}
fn brownian(n: usize, seed: u64) -> Vec<f64> {
let noise = white_noise(n, seed);
let mut walk = Vec::with_capacity(n);
let mut sum = 0.0;
for v in noise {
sum += v;
walk.push(sum);
}
walk
}
#[test]
fn white_noise_alpha_near_half() {
let data = white_noise(4096, 42);
let result = dfa(&data);
assert!(result.alpha > 0.35 && result.alpha < 0.65,
"white noise DFA alpha should be near 0.5, got {}", result.alpha);
assert!(result.r_squared > 0.8, "R2 should be high, got {}", result.r_squared);
}
#[test]
fn brownian_alpha_above_one() {
let data = brownian(4096, 42);
let result = dfa(&data);
assert!(result.alpha > 1.2 && result.alpha < 1.8,
"brownian DFA alpha should be near 1.5, got {}", result.alpha);
}
#[test]
fn deterministic() {
let data = white_noise(1024, 7);
let r1 = dfa(&data);
let r2 = dfa(&data);
assert!((r1.alpha - r2.alpha).abs() < 1e-10);
}
#[test]
fn too_short_returns_half() {
let data = [1.0; 10];
let result = dfa(&data);
assert_eq!(result.alpha, 0.5);
assert_eq!(result.r_squared, 0.0);
}
#[test]
fn analyze_produces_quality() {
let data = white_noise(2048, 7);
let law = analyze(&data);
assert_eq!(law.n, 2048);
assert!(law.quality != LawQuality::Insufficient);
}
#[test]
fn health_verdict_thresholds() {
assert_eq!(HealthVerdict::from_shift(0.01), HealthVerdict::Healthy);
assert_eq!(HealthVerdict::from_shift(0.05), HealthVerdict::Watch);
assert_eq!(HealthVerdict::from_shift(0.10), HealthVerdict::Warning);
assert_eq!(HealthVerdict::from_shift(0.20), HealthVerdict::Critical);
assert_eq!(HealthVerdict::from_shift(-0.20), HealthVerdict::Critical);
}
#[test]
fn acr_detects_correlation() {
let data = brownian(2048, 99);
let result = acr(&data);
assert!(result.alpha < -0.05, "brownian ACR exponent should be negative, got {}", result.alpha);
}
}