sparse-ir 0.8.4

// Comparison tests for FiniteTempBasis against Julia's SparseIR.jl
//
// Reference data generated by SparseIR.jl/test/generate_basis_reference.jl

use num_complex::Complex64;
use sparse_ir::basis::{BosonicBasis, FermionicBasis};
use sparse_ir::kernel::LogisticKernel;

/// Parse reference data file generated by Julia
struct BasisReferenceData {
    lambda: f64,
    beta: f64,
    epsilon: f64,
    size_f: usize,
    size_b: usize,
    accuracy: f64,
    svals: Vec<f64>,
    tau_points: Vec<f64>,
    u_tau_f: Vec<Vec<f64>>, // u_tau_f[tau_idx][l]
    omega_points: Vec<f64>,
    v_omega_f: Vec<Vec<f64>>, // v_omega_f[omega_idx][l]
    wn_f: Vec<i64>,
    uhat_wn_f: Vec<Vec<Complex64>>, // uhat_wn_f[wn_idx][l]
    wn_b: Vec<i64>,
    uhat_wn_b: Vec<Vec<Complex64>>, // uhat_wn_b[wn_idx][l]
}

impl BasisReferenceData {
    fn load(lambda: f64, beta: f64) -> Result<Self, String> {
        // Embed reference data as const arrays
        if (lambda - 10.0).abs() < 1e-10 && (beta - 1.0).abs() < 1e-10 {
            Self::load_lambda_10_beta_1()
        } else if (lambda - 1000.0).abs() < 1e-10 && (beta - 100.0).abs() < 1e-10 {
            Self::load_lambda_1000_beta_100()
        } else {
            Err(format!(
                "Reference data for lambda={}, beta={} not yet embedded",
                lambda, beta
            ))
        }
    }

    fn load_lambda_10_beta_1() -> Result<Self, String> {
        // From Julia generated data for lambda=10, epsilon=1e-6, beta=1.0
        let size_f = 10;
        let size_b = 10;
        let accuracy = 5.767491903071612e-7;

        let svals = vec![
            1.2621489299919288,
            0.8363588547029698,
            0.34626225858303167,
            0.1208262696776911,
            0.03387861935965415,
            0.007961300857785424,
            0.0015966925515801364,
            0.000278230515051532,
            4.276437930624658e-5,
            5.871774564004786e-6,
        ];

        // tau points and u evaluations (first 3 basis functions)
        let tau_points = vec![0.001, 0.01, 0.1, 0.5, 1.0];
        let u_tau_f = vec![
            vec![1.6150510987766558, -2.3413582595599562, 2.5603789134045214],
            vec![1.5652063559508924, -2.225808718979322, 2.3083934891804065],
            vec![1.1914850621030963, -1.3619442766806016, 0.5523820605012136],
            vec![
                0.7328276018851676,
                2.489947588455659e-16,
                -1.0916200577894015,
            ],
            vec![1.620762804489461, 2.3546090586947743, 2.5894771278834905],
        ];

        // omega points and v evaluations (first 3 basis functions)
        let omega_points = vec![-9.0, -5.0, -1.0, 0.0, 1.0, 5.0, 9.0];
        let v_omega_f = vec![
            vec![0.109170234170623, 0.19354184563534282, 0.2698463713103008],
            vec![0.1742626193502517, 0.25643909755939975, 0.1875881684252792],
            vec![0.3582767334465639, 0.15354671326496244, -0.2788092490432718],
            vec![
                0.38415493340064055,
                -7.530568469355681e-16,
                -0.3514988084390997,
            ],
            vec![
                0.3582767334465639,
                -0.15354671326496244,
                -0.2788092490432718,
            ],
            vec![0.1742626193502517, -0.25643909755939975, 0.1875881684252792],
            vec![0.109170234170623, -0.19354184563534282, 0.2698463713103008],
        ];

        // Matsubara frequencies and uhat evaluations (first 3 basis functions)
        let wn_f = vec![1, 3, 5, 11, 21, 41]; // 2n+1 for n=0,1,2,5,10,20
        let uhat_wn_f = vec![
            vec![
                Complex64::new(1.7177310195820194e-17, 0.543551067457796),
                Complex64::new(-0.6671120739157604, -1.5821646943056645e-17),
                Complex64::new(2.880567412410956e-20, -0.4423475462221159),
            ],
            vec![
                Complex64::new(-2.469903871159752e-18, 0.28853626564843604),
                Complex64::new(-0.20486609375458592, 2.7239837749935865e-17),
                Complex64::new(1.9508154554485073e-17, 0.23201183041101947),
            ],
            vec![
                Complex64::new(2.418872417345561e-18, 0.19113017095856463),
                Complex64::new(-0.09139131020081802, 1.3560098009823616e-17),
                Complex64::new(-1.3676429850661748e-19, 0.23699196608598908),
            ],
            vec![
                Complex64::new(4.914849321953196e-18, 0.09215005086046042),
                Complex64::new(-0.02143244711282436, 1.1485250412217592e-19),
                Complex64::new(1.4853741823290678e-17, 0.13946051164223078),
            ],
            vec![
                Complex64::new(-2.467493742296196e-18, 0.048889074583760746),
                Complex64::new(-0.00604338552056134, -2.319202543574785e-18),
                Complex64::new(9.643147758641022e-20, 0.07694677102954463),
            ],
            vec![
                Complex64::new(-6.098764611611493e-20, 0.025132911746251308),
                Complex64::new(-0.0015979446460983396, 7.133563859676406e-18),
                Complex64::new(2.2733461398948136e-19, 0.0399965133423288),
            ],
        ];

        let wn_b = vec![0, 2, 4, 10, 20, 40]; // 2n for n=0,1,2,5,10,20
        let uhat_wn_b = vec![
            vec![
                Complex64::new(0.969721476285501, 0.0),
                Complex64::new(3.434577253370758e-17, 0.0),
                Complex64::new(-0.24342154259873408, 0.0),
            ],
            vec![
                Complex64::new(0.15817285039258278, 1.382311310654452e-18),
                Complex64::new(9.815774307188666e-18, -0.44584111500844953),
                Complex64::new(0.6035472465618628, 4.28830202131642e-18),
            ],
            vec![
                Complex64::new(0.058237074706870486, 1.563714564996459e-17),
                Complex64::new(9.69580319091719e-18, -0.31121020547898665),
                Complex64::new(0.26850802448021305, 6.84217525202855e-18),
            ],
            vec![
                Complex64::new(0.011135832227335999, 3.0762335260910614e-18),
                Complex64::new(-2.6681103083606612e-20, -0.14473679864333844),
                Complex64::new(0.055726719914982535, -1.264477217713443e-18),
            ],
            vec![
                Complex64::new(0.0028715562787866135, 1.217154226790379e-18),
                Complex64::new(-1.149604208734437e-19, -0.07427690022468018),
                Complex64::new(0.014570612954625311, -1.4682589596115627e-17),
            ],
            vec![
                Complex64::new(0.0007236816949663905, 2.0218361972563335e-20),
                Complex64::new(-3.1904750290372817e-20, -0.03738977751379602),
                Complex64::new(0.0036852014767584963, -7.286889522932468e-18),
            ],
        ];

        Ok(BasisReferenceData {
            lambda: 10.0,
            beta: 1.0,
            epsilon: 1e-6,
            size_f,
            size_b,
            accuracy,
            svals,
            tau_points,
            u_tau_f,
            omega_points,
            v_omega_f,
            wn_f,
            uhat_wn_f,
            wn_b,
            uhat_wn_b,
        })
    }

    fn load_lambda_1000_beta_100() -> Result<Self, String> {
        // From Julia generated data for lambda=1000, beta=100, epsilon=1e-6
        // Note: omega_max = lambda/beta = 10.0
        let size_f = 34;
        let size_b = 34;
        let accuracy = 9.516587022699077e-7;

        let svals = vec![
            1.5511081043435417,
            1.4289129641805978,
            1.0588362836617466,
            0.8469455306384281,
            0.6030885452053223,
            0.44256246787047904,
            0.31078628341824566,
            0.21894909354401765,
            0.15151295568369325,
            0.10432666006127535,
        ];

        // tau points (scaled by beta=100)
        let tau_points = vec![0.1, 1.0, 10.0, 50.0, 100.0];
        let u_tau_f = vec![
            vec![0.674804399067457, -0.8743028071547695, 0.8957813311031454],
            vec![
                0.28760450033599283,
                -0.3139521934991705,
                0.07245425088390034,
            ],
            vec![
                0.08271254932070381,
                -0.05357592144738383,
                -0.08868595336857119,
            ],
            vec![
                0.03842202880779884,
                2.1581892743528232e-17,
                -0.06087935939023406,
            ],
            vec![0.8591411662160833, 1.14929985611491, 1.3613030053254767],
        ];

        // omega points (omega_max = 10.0)
        let omega_points = vec![-9.0, -5.0, -1.0, 0.0, 1.0, 5.0, 9.0];
        let v_omega_f = vec![
            vec![0.04940224704932481, 0.06975836605821752, 0.0985619438278353],
            vec![
                0.07946346516860682,
                0.11045323877599142,
                0.14496799491300374,
            ],
            vec![
                0.24852828200268648,
                0.31984343636587476,
                0.27510190333344414,
            ],
            vec![
                2.323910839690569,
                -4.168737253615449e-14,
                -2.8877880823356694,
            ],
            vec![
                0.24852828200268648,
                -0.31984343636587476,
                0.27510190333344414,
            ],
            vec![
                0.07946346516860682,
                -0.11045323877599142,
                0.14496799491300374,
            ],
            vec![
                0.04940224704932481,
                -0.06975836605821752,
                0.0985619438278353,
            ],
        ];

        // Matsubara frequencies
        let wn_f = vec![1, 3, 5, 11, 21, 41];
        let uhat_wn_f = vec![
            vec![
                Complex64::new(-2.6241963311660627e-16, 3.136396927990673),
                Complex64::new(-4.089103284229498, 3.0141079499151143e-17),
                Complex64::new(-3.165732792925321e-16, -4.325665692771434),
            ],
            vec![
                Complex64::new(-5.143494217906491e-17, 2.167092602975325),
                Complex64::new(-2.655603341442736, 1.1538817323013114e-16),
                Complex64::new(1.1754788942676862e-16, -1.9923776960743937),
            ],
            vec![
                Complex64::new(1.0317045396797477e-16, 1.7572943649594648),
                Complex64::new(-2.064452575066733, 2.2813563253687663e-17),
                Complex64::new(-1.3006779284083324e-16, -1.140403056568134),
            ],
            vec![
                Complex64::new(-2.9555663711048974e-17, 1.2263114460798898),
                Complex64::new(-1.3218746351405508, 8.847193257760363e-17),
                Complex64::new(4.670124285968068e-17, -0.23384803706626567),
            ],
            vec![
                Complex64::new(2.5737242692819774e-18, 0.8864843180353076),
                Complex64::new(-0.868606915645137, -6.417572178813686e-17),
                Complex64::new(5.6114138351239e-17, 0.17638150897109875),
            ],
            vec![
                Complex64::new(-3.9423719915883975e-17, 0.616093509254019),
                Complex64::new(-0.5283900588410059, 1.527634547456741e-17),
                Complex64::new(1.477144200184723e-17, 0.36465161306162097),
            ],
        ];

        let wn_b = vec![0, 2, 4, 10, 20, 40];
        let uhat_wn_b = vec![
            vec![
                Complex64::new(7.2092738744317275, 0.0),
                Complex64::new(-3.189250306701418e-17, 0.0),
                Complex64::new(-6.11538960220523, 0.0),
            ],
            vec![
                Complex64::new(2.753370433551026, 1.514906638504372e-16),
                Complex64::new(7.361743927124106e-17, -1.6139218354280587),
                Complex64::new(0.4819142542331083, -1.3841777371174086e-16),
            ],
            vec![
                Complex64::new(1.8800642824818488, 3.8972137871760894e-17),
                Complex64::new(2.502720150965713e-16, -1.5417418751317504),
                Complex64::new(1.099351067367909, -8.24794758448945e-17),
            ],
            vec![
                Complex64::new(1.0607021168670052, -2.587184290267125e-17),
                Complex64::new(-1.9852960036320536e-16, -1.2482907198970428),
                Complex64::new(1.2111396289306668, -6.511452582875415e-17),
            ],
            vec![
                Complex64::new(0.6531367734860059, -1.4399749173858074e-16),
                Complex64::new(1.2296508171726334e-16, -0.9806413148224159),
                Complex64::new(1.00377580855312, -3.157348368783964e-17),
            ],
            vec![
                Complex64::new(0.3796438357226972, -9.311405998530255e-17),
                Complex64::new(-1.9737487362497317e-17, -0.7280160398209289),
                Complex64::new(0.7158525799058675, -9.014761249120269e-18),
            ],
        ];

        Ok(BasisReferenceData {
            lambda: 1000.0,
            beta: 100.0,
            epsilon: 1e-6,
            size_f,
            size_b,
            accuracy,
            svals,
            tau_points,
            u_tau_f,
            omega_points,
            v_omega_f,
            wn_f,
            uhat_wn_f,
            wn_b,
            uhat_wn_b,
        })
    }
}

#[test]
fn test_basis_size_lambda_10_beta_1() {
    let ref_data = BasisReferenceData::load(10.0, 1.0).unwrap();

    let kernel = LogisticKernel::new(ref_data.lambda);
    let basis_f = FermionicBasis::new(kernel, ref_data.beta, Some(ref_data.epsilon), None);
    let basis_b = BosonicBasis::new(kernel, ref_data.beta, Some(ref_data.epsilon), None);

    println!(
        "Fermionic basis size: {} (expected {})",
        basis_f.size(),
        ref_data.size_f
    );
    println!(
        "Bosonic basis size: {} (expected {})",
        basis_b.size(),
        ref_data.size_b
    );

    assert_eq!(
        basis_f.size(),
        ref_data.size_f,
        "Fermionic basis size mismatch"
    );
    assert_eq!(
        basis_b.size(),
        ref_data.size_b,
        "Bosonic basis size mismatch"
    );
}

#[test]
fn test_basis_singular_values_lambda_10_beta_1() {
    let ref_data = BasisReferenceData::load(10.0, 1.0).unwrap();

    let kernel = LogisticKernel::new(ref_data.lambda);
    let basis_f = FermionicBasis::new(kernel, ref_data.beta, Some(ref_data.epsilon), None);

    println!("\nSingular values comparison:");
    println!("  i    Rust              Julia             |diff|");
    println!("  {}", "-".repeat(60));

    let tol = 1e-10;
    for i in 0..ref_data.svals.len() {
        let s_rust = basis_f.s()[i];
        let s_julia = ref_data.svals[i];
        let diff = (s_rust - s_julia).abs();
        println!(
            "  {:2}   {:.15e}  {:.15e}  {:.2e}",
            i, s_rust, s_julia, diff
        );
        assert!(
            diff < tol,
            "Singular value s[{}] mismatch: |{} - {}| = {} >= {}",
            i,
            s_rust,
            s_julia,
            diff,
            tol
        );
    }
}

#[test]
fn test_basis_u_tau_lambda_10_beta_1() {
    let ref_data = BasisReferenceData::load(10.0, 1.0).unwrap();

    let kernel = LogisticKernel::new(ref_data.lambda);
    let basis_f = FermionicBasis::new(kernel, ref_data.beta, Some(ref_data.epsilon), None);

    println!("\nu(tau) comparison (first 3 basis functions):");

    let tol = 1e-10;
    for (tau_idx, &tau) in ref_data.tau_points.iter().enumerate() {
        println!("\n  tau = {}", tau);
        for l in 0..3.min(basis_f.size()) {
            let u_rust = basis_f.u()[l].evaluate(tau);
            let u_julia = ref_data.u_tau_f[tau_idx][l];
            let diff = (u_rust - u_julia).abs();
            println!(
                "    u[{}] = {:.15e} (Julia: {:.15e}, |diff| = {:.2e})",
                l, u_rust, u_julia, diff
            );
            assert!(
                diff < tol,
                "u[{}](tau={}) mismatch: |{} - {}| = {} >= {}",
                l,
                tau,
                u_rust,
                u_julia,
                diff,
                tol
            );
        }
    }
}

#[test]
fn test_basis_v_omega_lambda_10_beta_1() {
    let ref_data = BasisReferenceData::load(10.0, 1.0).unwrap();

    let kernel = LogisticKernel::new(ref_data.lambda);
    let basis_f = FermionicBasis::new(kernel, ref_data.beta, Some(ref_data.epsilon), None);

    println!("\nv(omega) comparison (first 3 basis functions):");

    let tol = 1e-10;
    for (omega_idx, &omega) in ref_data.omega_points.iter().enumerate() {
        println!("\n  omega = {}", omega);
        for l in 0..3.min(basis_f.size()) {
            let v_rust = basis_f.v()[l].evaluate(omega);
            let v_julia = ref_data.v_omega_f[omega_idx][l];
            let diff = (v_rust - v_julia).abs();
            println!(
                "    v[{}] = {:.15e} (Julia: {:.15e}, |diff| = {:.2e})",
                l, v_rust, v_julia, diff
            );
            assert!(
                diff < tol,
                "v[{}](omega={}) mismatch: |{} - {}| = {} >= {}",
                l,
                omega,
                v_rust,
                v_julia,
                diff,
                tol
            );
        }
    }
}

#[test]
fn test_basis_uhat_wn_fermionic_lambda_10_beta_1() {
    let ref_data = BasisReferenceData::load(10.0, 1.0).unwrap();

    let kernel = LogisticKernel::new(ref_data.lambda);
    let basis_f = FermionicBasis::new(kernel, ref_data.beta, Some(ref_data.epsilon), None);

    println!("\nuhat(wn) comparison for Fermionic (first 3 basis functions):");

    let tol = 1e-10;
    for (wn_idx, &wn) in ref_data.wn_f.iter().enumerate() {
        println!("\n  wn = {} (n={})", wn, (wn - 1) / 2);
        for l in 0..3.min(basis_f.size()) {
            let uhat_rust = basis_f.uhat().polyvec[l].evaluate_at_n(wn);
            let uhat_julia = ref_data.uhat_wn_f[wn_idx][l];
            let diff = (uhat_rust - uhat_julia).norm();
            println!(
                "    uhat[{}] = {:.6e} + {:.6e}i (Julia: {:.6e} + {:.6e}i, |diff| = {:.2e})",
                l, uhat_rust.re, uhat_rust.im, uhat_julia.re, uhat_julia.im, diff
            );
            assert!(
                diff < tol,
                "uhat[{}](wn={}) mismatch: |({} + {}i) - ({} + {}i)| = {} >= {}",
                l,
                wn,
                uhat_rust.re,
                uhat_rust.im,
                uhat_julia.re,
                uhat_julia.im,
                diff,
                tol
            );
        }
    }
}

#[test]
fn test_basis_uhat_wn_bosonic_lambda_10_beta_1() {
    let ref_data = BasisReferenceData::load(10.0, 1.0).unwrap();

    let kernel = LogisticKernel::new(ref_data.lambda);
    let basis_b = BosonicBasis::new(kernel, ref_data.beta, Some(ref_data.epsilon), None);

    println!("\nuhat(wn) comparison for Bosonic (first 3 basis functions):");

    let tol = 1e-10;
    for (wn_idx, &wn) in ref_data.wn_b.iter().enumerate() {
        println!("\n  wn = {} (n={})", wn, wn / 2);
        for l in 0..3.min(basis_b.size()) {
            let uhat_rust = basis_b.uhat().polyvec[l].evaluate_at_n(wn);
            let uhat_julia = ref_data.uhat_wn_b[wn_idx][l];
            let diff = (uhat_rust - uhat_julia).norm();
            println!(
                "    uhat[{}] = {:.6e} + {:.6e}i (Julia: {:.6e} + {:.6e}i, |diff| = {:.2e})",
                l, uhat_rust.re, uhat_rust.im, uhat_julia.re, uhat_julia.im, diff
            );
            assert!(
                diff < tol,
                "uhat[{}](wn={}) mismatch: |({} + {}i) - ({} + {}i)| = {} >= {}",
                l,
                wn,
                uhat_rust.re,
                uhat_rust.im,
                uhat_julia.re,
                uhat_julia.im,
                diff,
                tol
            );
        }
    }
}

// ========================================
// Tests for beta=100, lambda=1000 (omega_max=10)
// ========================================

#[test]
fn test_basis_size_lambda_1000_beta_100() {
    let ref_data = BasisReferenceData::load(1000.0, 100.0).unwrap();

    let kernel = LogisticKernel::new(ref_data.lambda);
    let basis_f = FermionicBasis::new(kernel, ref_data.beta, Some(ref_data.epsilon), None);
    let basis_b = BosonicBasis::new(kernel, ref_data.beta, Some(ref_data.epsilon), None);

    println!("\nbeta=100, lambda=1000 (omega_max=10):");
    println!(
        "  Fermionic basis size: {} (expected {})",
        basis_f.size(),
        ref_data.size_f
    );
    println!(
        "  Bosonic basis size: {} (expected {})",
        basis_b.size(),
        ref_data.size_b
    );

    assert_eq!(
        basis_f.size(),
        ref_data.size_f,
        "Fermionic basis size mismatch"
    );
    assert_eq!(
        basis_b.size(),
        ref_data.size_b,
        "Bosonic basis size mismatch"
    );
}

#[test]
fn test_basis_uhat_wn_lambda_1000_beta_100() {
    let ref_data = BasisReferenceData::load(1000.0, 100.0).unwrap();

    let kernel = LogisticKernel::new(ref_data.lambda);
    let basis_f = FermionicBasis::new(kernel, ref_data.beta, Some(ref_data.epsilon), None);
    let basis_b = BosonicBasis::new(kernel, ref_data.beta, Some(ref_data.epsilon), None);

    let tol = 1e-10;

    // Fermionic - check first basis function at wn=1
    let uhat_rust = basis_f.uhat().polyvec[0].evaluate_at_n(ref_data.wn_f[0]);
    let uhat_julia = ref_data.uhat_wn_f[0][0];
    let diff = (uhat_rust - uhat_julia).norm();
    assert!(
        diff < tol,
        "uhat_f[0](wn=1) mismatch: diff={:.2e} >= {}",
        diff,
        tol
    );

    // Bosonic - check first basis function at wn=0
    let uhat_rust = basis_b.uhat().polyvec[0].evaluate_at_n(ref_data.wn_b[0]);
    let uhat_julia = ref_data.uhat_wn_b[0][0];
    let diff = (uhat_rust - uhat_julia).norm();
    assert!(
        diff < tol,
        "uhat_b[0](wn=0) mismatch: diff={:.2e} >= {}",
        diff,
        tol
    );
}