csrk 1.1.4

/// Implementation tests for random number generators
#[test]
#[ignore]
fn next_f64_performance() {
    use csrk::PCG64Stream;
    use std::time::Instant;
    // Get a random number generator
    let mut rs = PCG64Stream::from_OsRng();
    let doubles = 1_000_000;
    let draws = 1_000;
    // Generate some running means
    let mut mean_inline:     f64 = 0.0;
    let mut mean_vec:         f64 = 0.0;
    let mut mean_fill_fast: f64 = 0.0;
    let mut mean_fill_slow: f64 = 0.0;
    let mut time_inline:     u128 = 0;
    let mut time_vec:         u128 = 0;
    let mut time_fill_fast: u128 = 0;
    let mut time_fill_slow: u128 = 0;
    let mut size_inline:     usize = 0;
    let mut size_vec:         usize = 0;
    let mut size_fill_fast: usize = 0;
    let mut size_fill_slow: usize = 0;
    // Inline test
    for _ in 0..draws {
        let t0 = Instant::now();
        let uniform = (0..doubles).map(
                |_| rs.next_f64()
            ).collect::<Vec<f64>>();
        time_inline += t0.elapsed().as_nanos();
        mean_inline += uniform.iter().sum::<f64>();
        size_inline += uniform.len();
    }
    // Vector test
    for _ in 0..draws {
        let t0 = Instant::now();
        let uniform = rs.vec_f64(doubles);
        time_vec += t0.elapsed().as_nanos();
        mean_vec += uniform.iter().sum::<f64>();
        size_vec += uniform.len();
    }
    // fill fast test
    for _ in 0..draws {
        let mut uniform = vec![0.;doubles];
        let t0 = Instant::now();
        rs.fill_f64(&mut uniform);
        time_fill_fast += t0.elapsed().as_nanos();
        mean_fill_fast += uniform.iter().sum::<f64>();
        size_fill_fast += uniform.len();
    }
    // fill slow test
    for _ in 0..draws {
        let t0 = Instant::now();
        let mut uniform = vec![0.;doubles];
        rs.fill_f64(&mut uniform);
        time_fill_slow += t0.elapsed().as_nanos();
        mean_fill_slow += uniform.iter().sum::<f64>();
        size_fill_slow += uniform.len();
    }
    let sec_inline    = (time_inline    as f64) / 1_000_000_000.0;
    let sec_vec       = (time_vec       as f64) / 1_000_000_000.0;
    let sec_fill_fast = (time_fill_fast as f64) / 1_000_000_000.0;
    let sec_fill_slow = (time_fill_slow as f64) / 1_000_000_000.0;
    mean_inline /= size_inline as f64;
    mean_vec /= size_vec as f64;
    mean_fill_fast /= size_fill_fast as f64;
    mean_fill_slow /= size_fill_slow as f64;
    println!(
        "Method inline generated {:.1e} doubles in {} sec (with mean: {})",
        size_inline,
        sec_inline,
        mean_inline,
    );
    println!(
        "Method vec generated {:.1e} doubles in {} sec (with mean: {})",
        size_vec,
        sec_vec,
        mean_vec,
    );
    println!(
        "Method fill_fast generated {:.1e} doubles in {} sec (with mean: {})",
        size_fill_fast,
        sec_fill_fast,
        mean_fill_fast,
    );
    println!(
        "Method fill_slow generated {:.1e} doubles in {} sec (with mean: {})",
        size_fill_slow,
        sec_fill_slow,
        mean_fill_slow,
    );
}
#[test]
#[ignore]
fn gaussian_implementation_test() {
    // Imports at the top
    use csrk::PCG64Stream;
    use std::time::Instant;
    // Parameters
    let n_realizations = 200;   // number of independent runs
    let n_draws = 50_000;       // samples per run
    // Initialize stats
    let mut means   = Vec::with_capacity(n_realizations);
    let mut vars    = Vec::with_capacity(n_realizations);
    // Keep track of time
    let t0 = Instant::now();
    // Loop realizations
    for _ in 0..n_realizations {
        // Define a random number stream
        let mut rs = PCG64Stream::from_OsRng();
        // Initialize normal draws
        let mut norm = vec![0.;n_draws];
        rs.fill_standard_normal(&mut norm);
        // do statistics
        let mut sum = 0.;
        let mut sum2 = 0.;
        for &x in &norm {
            sum += x;
            sum2 += x * x;
        }
        let mean = sum / n_draws as f64;
        let var = sum2 / n_draws as f64 - (mean * mean);
        // Update our statistics vectors for the whole test
        means.push(mean);
        vars.push(var);
    }
    // Keep track of time
    let elapsed = t0.elapsed();
    println!(
        "Generated {} realizations of {} samples in {:?}",
        n_realizations, n_draws, elapsed,
    );
    // empirical stats
    let mean_of_means = means.iter().sum::<f64>() / n_realizations as f64;
    let mean_of_vars  = vars.iter().sum::<f64>()  / n_realizations as f64;
    let var_means = means.iter()
        .map(|m| (m - mean_of_means).powi(2))
        .sum::<f64>() / n_realizations as f64;
    let var_vars = vars.iter()
        .map(|v| (v - mean_of_vars).powi(2))
        .sum::<f64>() / n_realizations as f64;
    let std_mean_emp = var_means.sqrt();
    let std_var_emp  = var_vars.sqrt();
    // Analytic expresions for these quantities
    let std_mean_expected = 1.0 / (n_draws as f64).sqrt();
    let std_var_expected  = (2.0 / n_draws as f64).sqrt();
    println!("Mean of means     = {:.4e}", mean_of_means);
    println!("Mean of variances = {:.6e}", mean_of_vars);
    println!("Std(mean): empirical {:.4e}, theory {:.4e}",
        std_mean_emp, std_mean_expected);
    println!("Std(var):  empirical {:.4e}, theory {:.4e}",
        std_var_emp,  std_var_expected);
    // relative error checks
    let rel_err_mean = (std_mean_emp - std_mean_expected).abs() / 
        std_mean_expected;
    let rel_err_var  = (std_var_emp  - std_var_expected ).abs() / 
        std_var_expected;
    println!("Relative error std(mean): {:.2}%", 100.0 * rel_err_mean);
    println!("Relative error std(var):  {:.2}%", 100.0 * rel_err_var);
    assert!(rel_err_mean < 0.25, "mean spread off by too much");
    assert!(rel_err_var < 0.30, "var spread off by too much");
    // tail frequency check
    // amount things can fail by
    let k = 3.0; // number of sigma
    let mut mean_exceed = 0;
    let mut var_exceed  = 0;
    // loop means
    for &m in &means {
        if m.abs() > k * std_mean_expected {
            mean_exceed += 1;
        }
    }
    // loop vars
    for &v in &vars {
        if (v - 1.0).abs() > k * std_var_expected {
            var_exceed += 1;
        }
    }
    // Identify the fraction of failures
    let frac_mean = mean_exceed as f64 / n_realizations as f64;
    let frac_var  = var_exceed  as f64 / n_realizations as f64;
    // Test against 3 sigma
    println!("3σ exceedance: mean {:.3}%, var {:.3}%",
        100. * frac_mean, 100. * frac_var);
    assert!(frac_mean < 0.02, "too many extreme sample means");
    assert!(frac_var  < 0.05, "too many extreme variances");
    // Identify the chisquare value
    let sigma_m2  = 1.0 / n_draws as f64;
    let chi2_means: f64 = means.iter()
        .map(|&m| m*m / sigma_m2)
        .sum();
    println!(
        "χ^2 means = {:.2} for {} dof (χ^2/dof = {:.3})",
        chi2_means,
        n_realizations,
        chi2_means / n_realizations as f64,
    );
    // Identify the chi-square value for the variances
    let sigma_v2 = 2.0 / n_draws as f64;
    let chi2_vars: f64 = vars.iter()
        .map(|&v| (v - 1.0)*(v - 1.0) / sigma_v2)
        .sum();
    println!(
        "χ^2 vars  = {:.2} for {} dof (chi^2/dof = {:.3})",
        chi2_vars,
        n_realizations,
        chi2_vars / n_realizations as f64,
    );

}
#[test]
#[ignore]
fn gaussian_two_point_correlation_integration() {
    use csrk::PCG64Stream;
    use std::time::Instant;

    let n_realizations = 80;
    let n_draws = 100_000;
    let max_lag = 40;

    let sigma2 = 1.0 / n_draws as f64; // Var(C_k)

    let mut all_corrs = vec![vec![0.0; max_lag]; n_realizations];

    let t0 = Instant::now();

    for r in 0..n_realizations {
        let mut rs = PCG64Stream::from_OsRng();
        let mut x = vec![0.0; n_draws];
        rs.fill_standard_normal(&mut x);

        for lag in 1..=max_lag {
            let mut sum = 0.0;
            let mut i = lag;
            while i < n_draws {
                sum += x[i] * x[i - lag];
                i += 1;
            }
            all_corrs[r][lag - 1] = sum / (n_draws - lag) as f64;
        }
    }

    let elapsed = t0.elapsed();
    println!(
        "Computed {} realizations of {} lags (N={}) in {:?}",
        n_realizations, max_lag, n_draws, elapsed
    );

    // ---- Per-lag statistics ----

    let mut chi2 = 0.0;

    for lag in 0..max_lag {
        let mut mean = 0.0;
        for r in 0..n_realizations {
            mean += all_corrs[r][lag];
        }
        mean /= n_realizations as f64;

        let mut var = 0.0;
        for r in 0..n_realizations {
            let d = all_corrs[r][lag] - mean;
            var += d * d;
        }
        var /= n_realizations as f64;

        let std_emp = var.sqrt();
        let std_th = sigma2.sqrt();

        println!(
            "lag {:2}: mean={:+.3e}, std_emp={:.3e}, std_th={:.3e}, ratio={:.3}",
            lag + 1,
            mean,
            std_emp,
            std_th,
            std_emp / std_th
        );

        // contribute to chi-square (means should be 0)
        for r in 0..n_realizations {
            let c = all_corrs[r][lag];
            chi2 += c * c / sigma2;
        }
    }

    let dof = (n_realizations * max_lag) as f64;
    let chi2_per = chi2 / dof;
    let tol = 5.0 * (2.0 / dof).sqrt();

    println!(
        "Total chi^2 = {:.1} for {} dof (chi^2/dof = {:.4}, tol ≈ {:.4})",
        chi2, dof, chi2_per, tol
    );

    assert!((chi2_per - 1.0).abs() < tol, "global correlation chi^2 too large");
}