//! Comprehensive TDD tests for real scientific computing algorithms
#![cfg(feature = "disabled_tests")]

use crate::ScienceError;
use crate::scientific_computing::*;
use nalgebra as na;

#[cfg(test)]
mod scientific_tensor_tests {
    use super::*;

    #[test]
    fn test_scientific_tensor_creation() {
        let data = ndarray::Array2::from_elem((3, 3), 1.0);
        let tensor = ScientificTensor::from_array(data.into_dyn());

        assert_eq!(tensor.shape(), &[3, 3]);
        assert_eq!(tensor.data().sum(), 9.0);
    }

    #[test]
    fn test_tensor_with_units() {
        let data = ndarray::Array1::from_vec(vec![1.0, 2.0, 3.0]);
        let tensor = ScientificTensor::with_units(data.into_dyn(), "meters");

        assert_eq!(tensor.units(), Some("meters"));
        assert_eq!(tensor.shape(), &[3]);
    }

    #[test]
    fn test_tensor_zeros_ones() {
        let zeros = ScientificTensor::zeros(&[2, 3]);
        assert_eq!(zeros.shape(), &[2, 3]);
        assert_eq!(zeros.data().sum(), 0.0);

        let ones = ScientificTensor::ones(&[2, 3]);
        assert_eq!(ones.data().sum(), 6.0);
    }

    #[test]
    fn test_tensor_randn() {
        let tensor = ScientificTensor::randn(&[1000]);
        assert_eq!(tensor.shape(), &[1000]);

        let data = tensor.data().as_slice().unwrap();
        let mean: f64 = data.iter().sum::<f64>() / data.len() as f64;
        let variance: f64 =
            data.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / data.len() as f64;

        // Should approximate N(0,1)
        assert!(mean.abs() < 0.2);
        assert!(variance > 0.5 && variance < 2.0);
    }

    #[test]
    fn test_tensor_arithmetic() {
        let a =
            ScientificTensor::from_array(ndarray::Array1::from_vec(vec![1.0, 2.0, 3.0]).into_dyn());
        let b =
            ScientificTensor::from_array(ndarray::Array1::from_vec(vec![4.0, 5.0, 6.0]).into_dyn());

        let sum = a.add(&b).unwrap();
        assert_eq!(sum.data().as_slice().unwrap(), &[5.0, 7.0, 9.0]);

        let product = a.mul(&b).unwrap();
        assert_eq!(product.data().as_slice().unwrap(), &[4.0, 10.0, 18.0]);

        let scaled = a.mul_scalar(2.0);
        assert_eq!(scaled.data().as_slice().unwrap(), &[2.0, 4.0, 6.0]);
    }

    #[test]
    fn test_tensor_shape_mismatch() {
        let a = ScientificTensor::zeros(&[2, 3]);
        let b = ScientificTensor::zeros(&[3, 2]);

        let result = a.add(&b);
        assert!(result.is_err());
    }

    #[test]
    fn test_tensor_matmul() {
        let a = ScientificTensor::from_array(
            ndarray::Array2::from_shape_vec((2, 3), vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0])
                .unwrap()
                .into_dyn(),
        );
        let b = ScientificTensor::from_array(
            ndarray::Array2::from_shape_vec((3, 2), vec![7.0, 8.0, 9.0, 10.0, 11.0, 12.0])
                .unwrap()
                .into_dyn(),
        );

        let result = a.matmul(&b).unwrap();
        assert_eq!(result.shape(), &[2, 2]);

        // Verify matrix multiplication result
        let expected = ndarray::arr2(&[[58.0, 64.0], [139.0, 154.0]]);
        let result_array = result
            .data()
            .view()
            .into_dimensionality::<ndarray::Ix2>()
            .unwrap();
        assert_eq!(result_array, expected);
    }

    #[test]
    fn test_tensor_statistics() {
        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0];
        let tensor = ScientificTensor::from_array(ndarray::Array1::from_vec(data).into_dyn());

        assert_eq!(tensor.mean(), 3.0);
        assert!((tensor.std() - 1.58).abs() < 0.1);
        assert!((tensor.var() - 2.5).abs() < 0.1);
    }

    #[test]
    fn test_tensor_fft() {
        let data = vec![1.0, 0.0, 1.0, 0.0];
        let tensor = ScientificTensor::from_array(ndarray::Array1::from_vec(data).into_dyn());

        let fft_result = tensor.fft().unwrap();
        assert_eq!(fft_result.len(), 4);

        // Basic sanity check - FFT of alternating signal
        assert!(fft_result[0].re > 0.0);
    }

    #[test]
    fn test_tensor_gradient() {
        let data = vec![1.0, 4.0, 9.0, 16.0]; // x^2 for x = 1,2,3,4
        let tensor = ScientificTensor::from_array(ndarray::Array1::from_vec(data).into_dyn());

        let grad = tensor.gradient(0).unwrap();
        assert_eq!(grad.shape(), tensor.shape());
        // Gradient of x^2 should be approximately 2x
    }

    #[test]
    fn test_tensor_metadata() {
        let mut tensor = ScientificTensor::zeros(&[2, 2]);
        tensor.add_metadata("experiment", "test_001");
        tensor.add_metadata("date", "2023-01-01");

        assert_eq!(
            tensor.metadata.get("experiment"),
            Some(&"test_001".to_string())
        );
        assert_eq!(tensor.metadata.get("date"), Some(&"2023-01-01".to_string()));
    }
}

#[cfg(test)]
mod physics_simulation_tests {
    use super::*;

    #[test]
    fn test_physics_simulation_creation() {
        let sim = PhysicsSimulation::new(0.01);
        assert_eq!(sim.time_step, 0.01);
        assert_eq!(sim.current_time, 0.0);
    }

    #[test]
    fn test_harmonic_oscillator() {
        let mut sim = PhysicsSimulation::new(0.001);

        let results = sim.simulate_harmonic_oscillator(1.0, 1.0, 1.0).unwrap();
        assert!(!results.is_empty());

        // Check that we have time, position, velocity data
        let (time, pos, vel) = results[0];
        assert_eq!(time, 0.0);
        assert_eq!(pos, 1.0); // Initial position
        assert_eq!(vel, 0.0); // Initial velocity

        // Check energy conservation (approximately)
        let last_result = results.last().unwrap();
        let initial_energy = 0.5 * 1.0 * 1.0; // 1/2 * k * x^2 at t=0
        let final_energy = 0.5 * last_result.1.powi(2) + 0.5 * last_result.2.powi(2);
        assert!((initial_energy - final_energy).abs() < 0.1);
    }

    #[test]
    fn test_wave_equation() {
        let mut sim = PhysicsSimulation::new(0.001);

        let result = sim.simulate_wave_equation(1.0, 1.0, 0.1, 50);
        assert!(result.is_ok());

        let solution = result.unwrap();
        assert_eq!(solution.shape(), &[100, 50]); // nt x nx
    }

    #[test]
    fn test_wave_equation_cfl_violation() {
        let mut sim = PhysicsSimulation::new(0.1); // Large time step

        let result = sim.simulate_wave_equation(1.0, 1.0, 0.1, 10);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            ScienceError::Numerical { .. }
        ));
    }

    #[test]
    fn test_heat_equation() {
        let mut sim = PhysicsSimulation::new(0.001);

        let result = sim.simulate_heat_equation(0.1, 1.0, 0.1, 50);
        assert!(result.is_ok());

        let solution = result.unwrap();
        assert_eq!(solution.shape(), &[100, 50]); // nt x nx

        // Check that heat diffuses (temperature smooths out)
        let initial = solution.row(0);
        let final_row = solution.row(99);

        let initial_variance: f64 =
            initial.iter().map(|x| x.powi(2)).sum::<f64>() / initial.len() as f64;
        let final_variance: f64 =
            final_row.iter().map(|x| x.powi(2)).sum::<f64>() / final_row.len() as f64;

        assert!(final_variance < initial_variance); // Heat should diffuse
    }

    #[test]
    fn test_heat_equation_stability_violation() {
        let mut sim = PhysicsSimulation::new(0.1); // Large time step

        let result = sim.simulate_heat_equation(1.0, 1.0, 0.1, 10);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            ScienceError::Numerical { .. }
        ));
    }
}

#[cfg(test)]
mod chemistry_simulation_tests {
    use super::*;

    #[test]
    fn test_chemistry_simulation_creation() {
        let sim = ChemistrySimulation::new(298.15, 101325.0);
        assert_eq!(sim.temperature, 298.15);
        assert_eq!(sim.pressure, 101325.0);
    }

    #[test]
    fn test_water_molecule() {
        let mut sim = ChemistrySimulation::new(298.15, 101325.0);

        let atoms = vec![
            Atom {
                element: "O".to_string(),
                atomic_number: 8,
                position: na::Vector3::new(0.0, 0.0, 0.0),
                charge: -0.8,
            },
            Atom {
                element: "H".to_string(),
                atomic_number: 1,
                position: na::Vector3::new(0.96, 0.0, 0.0),
                charge: 0.4,
            },
            Atom {
                element: "H".to_string(),
                atomic_number: 1,
                position: na::Vector3::new(-0.24, 0.93, 0.0),
                charge: 0.4,
            },
        ];

        let water = Molecule {
            formula: "H2O".to_string(),
            molecular_weight: 18.015,
            atoms,
            bonds: vec![],
            geometry: ndarray::Array2::zeros((3, 3)),
        };

        sim.add_molecule("water", water);

        let properties = sim.calculate_molecular_properties("water").unwrap();
        assert!(properties.contains_key("molecular_weight"));
        assert!(properties.contains_key("center_of_mass_x"));
        assert!(properties.contains_key("moment_of_inertia"));

        assert_eq!(properties["molecular_weight"], 18.015);
    }

    #[test]
    fn test_reaction_kinetics() {
        let mut sim = ChemistrySimulation::new(298.15, 101325.0);

        // Add a simple reaction: A -> B
        let reaction = ChemicalReaction {
            reactants: vec!["A".to_string()],
            products: vec!["B".to_string()],
            rate_constant: 0.1,
            activation_energy: 50000.0, // J/mol
        };
        sim.reactions.push(reaction);

        let mut initial_concentrations = std::collections::HashMap::new();
        initial_concentrations.insert("A".to_string(), 1.0);
        initial_concentrations.insert("B".to_string(), 0.0);

        let results = sim
            .simulate_reaction_kinetics(initial_concentrations, 1.0)
            .unwrap();
        assert!(!results.is_empty());

        // Check that A decreases and B increases
        let initial = &results[0];
        let final_step = results.last().unwrap();

        assert!(final_step["A"] < initial["A"]);
        assert!(final_step["B"] > initial["B"]);
    }

    #[test]
    fn test_atomic_mass_lookup() {
        let sim = ChemistrySimulation::new(298.15, 101325.0);

        assert_eq!(sim.get_atomic_mass("H").unwrap(), 1.008);
        assert_eq!(sim.get_atomic_mass("C").unwrap(), 12.011);
        assert_eq!(sim.get_atomic_mass("O").unwrap(), 15.999);

        assert!(sim.get_atomic_mass("Xx").is_err());
    }
}

#[cfg(test)]
mod materials_simulation_tests {
    use super::*;

    #[test]
    fn test_materials_simulation_creation() {
        let crystal = CrystalStructure {
            lattice_parameters: [4.0, 4.0, 4.0, 90.0, 90.0, 90.0],
            space_group: "Pm-3m".to_string(),
            atoms: vec![],
        };

        let sim = MaterialsSimulation::new(crystal, 300.0, 101325.0);
        assert_eq!(sim.temperature, 300.0);
        assert_eq!(sim.pressure, 101325.0);
    }

    #[test]
    fn test_elastic_properties() {
        let crystal = CrystalStructure {
            lattice_parameters: [4.0, 4.0, 4.0, 90.0, 90.0, 90.0],
            space_group: "Pm-3m".to_string(),
            atoms: vec![],
        };

        let sim = MaterialsSimulation::new(crystal, 300.0, 101325.0);
        let properties = sim.calculate_elastic_properties().unwrap();

        assert!(properties.contains_key("bulk_modulus"));
        assert!(properties.contains_key("youngs_modulus"));
        assert!(properties.contains_key("poissons_ratio"));
        assert!(properties.contains_key("shear_modulus"));

        // Basic sanity checks
        let bulk_modulus = properties["bulk_modulus"];
        let youngs_modulus = properties["youngs_modulus"];
        let poissons_ratio = properties["poissons_ratio"];

        assert!(bulk_modulus > 0.0);
        assert!(youngs_modulus > 0.0);
        assert!(poissons_ratio > 0.0 && poissons_ratio < 0.5);
    }

    #[test]
    fn test_thermal_properties() {
        let crystal = CrystalStructure {
            lattice_parameters: [3.5, 3.5, 3.5, 90.0, 90.0, 90.0],
            space_group: "Fd-3m".to_string(),
            atoms: vec![],
        };

        let sim = MaterialsSimulation::new(crystal, 300.0, 101325.0);
        let properties = sim.calculate_thermal_properties().unwrap();

        assert!(properties.contains_key("debye_temperature"));
        assert!(properties.contains_key("heat_capacity"));
        assert!(properties.contains_key("thermal_expansion"));
        assert!(properties.contains_key("thermal_conductivity"));

        // Basic sanity checks
        assert!(properties["debye_temperature"] > 0.0);
        assert!(properties["heat_capacity"] > 0.0);
        assert!(properties["thermal_expansion"] > 0.0);
        assert!(properties["thermal_conductivity"] > 0.0);
    }

    #[test]
    fn test_temperature_dependence() {
        let crystal = CrystalStructure {
            lattice_parameters: [4.0, 4.0, 4.0, 90.0, 90.0, 90.0],
            space_group: "Pm-3m".to_string(),
            atoms: vec![],
        };

        let sim_low_t = MaterialsSimulation::new(crystal.clone(), 100.0, 101325.0);
        let sim_high_t = MaterialsSimulation::new(crystal, 1000.0, 101325.0);

        let props_low = sim_low_t.calculate_thermal_properties().unwrap();
        let props_high = sim_high_t.calculate_thermal_properties().unwrap();

        // Thermal conductivity should decrease with temperature
        assert!(props_high["thermal_conductivity"] < props_low["thermal_conductivity"]);
    }
}

#[cfg(test)]
mod numerical_methods_tests {
    use super::*;

    #[test]
    fn test_linear_system_solve() {
        // Solve: 2x + 3y = 7, x + y = 3
        let a = ndarray::arr2(&[[2.0, 3.0], [1.0, 1.0]]);
        let b = ndarray::arr1(&[7.0, 3.0]);

        let x = NumericalMethods::solve_linear_system(&a, &b).unwrap();

        // Solution should be x=2, y=1
        assert!((x[0] - 2.0).abs() < 1e-10);
        assert!((x[1] - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_eigenvalues() {
        // Simple 2x2 symmetric matrix
        let matrix = ndarray::arr2(&[[3.0, 1.0], [1.0, 3.0]]);

        let (eigenvals, _eigenvecs) = NumericalMethods::eigenvalues(&matrix).unwrap();

        // Eigenvalues should be 2 and 4
        let mut sorted_vals = eigenvals.to_vec();
        sorted_vals.sort_by(|a, b| a.total_cmp(b));

        assert!((sorted_vals[0] - 2.0).abs() < 1e-10);
        assert!((sorted_vals[1] - 4.0).abs() < 1e-10);
    }

    #[test]
    fn test_trapezoidal_integration() {
        // Integrate x^2 from 0 to 2
        let x = ndarray::Array1::linspace(0.0, 2.0, 1000);
        let y = x.mapv(|val| val * val);

        let integral = NumericalMethods::integrate_trapezoidal(&x, &y).unwrap();

        // Analytical result is 8/3 ≈ 2.667
        assert!((integral - 8.0 / 3.0).abs() < 0.01);
    }

    #[test]
    fn test_runge_kutta_ode() {
        // Solve dy/dt = -y with y(0) = 1
        // Analytical solution: y(t) = exp(-t)
        let f = |_t: f64, y: f64| -y;

        let (t, y) = NumericalMethods::runge_kutta_4(f, 1.0, (0.0, 2.0), 1000).unwrap();

        // Check solution at t=1: should be e^(-1) ≈ 0.368
        let idx = t.iter().position(|&x| (x - 1.0).abs() < 0.01).unwrap();
        let analytical = (-1.0f64).exp();

        assert!((y[idx] - analytical).abs() < 0.01);
    }

    #[test]
    fn test_newton_raphson() {
        // Find root of f(x) = x^2 - 4, which is x = 2
        let f = |x: f64| x * x - 4.0;
        let df = |x: f64| 2.0 * x;

        let root = NumericalMethods::newton_raphson(f, df, 3.0, 1e-10, 100).unwrap();

        assert!((root - 2.0).abs() < 1e-10);
    }

    #[test]
    fn test_newton_raphson_convergence_failure() {
        // Function with zero derivative
        let f = |x: f64| x * x;
        let df = |_x: f64| 0.0; // Always zero derivative

        let result = NumericalMethods::newton_raphson(f, df, 1.0, 1e-10, 100);
        assert!(result.is_err());
    }

    #[test]
    fn test_integration_error_cases() {
        let x = ndarray::arr1(&[1.0, 2.0, 3.0]);
        let y = ndarray::arr1(&[1.0, 2.0]); // Mismatched length

        let result = NumericalMethods::integrate_trapezoidal(&x, &y);
        assert!(result.is_err());

        let x_short = ndarray::arr1(&[1.0]);
        let y_short = ndarray::arr1(&[1.0]);

        let result = NumericalMethods::integrate_trapezoidal(&x_short, &y_short);
        assert!(result.is_err());
    }
}

#[cfg(test)]
mod integration_tests {
    use super::*;

    #[test]
    fn test_multiphysics_workflow() {
        // Test a combined physics-chemistry-materials workflow

        // 1. Physics: Simulate temperature distribution
        let mut physics_sim = PhysicsSimulation::new(0.001);
        let heat_solution = physics_sim
            .simulate_heat_equation(0.1, 1.0, 0.1, 50)
            .unwrap();

        assert_eq!(heat_solution.shape(), &[100, 50]);

        // 2. Chemistry: Use temperature for reaction kinetics
        let temp = 350.0; // K
        let mut chem_sim = ChemistrySimulation::new(temp, 101325.0);

        let reaction = ChemicalReaction {
            reactants: vec!["A".to_string()],
            products: vec!["B".to_string()],
            rate_constant: 1e6,         // Pre-exponential factor
            activation_energy: 50000.0, // J/mol
        };
        chem_sim.reactions.push(reaction);

        let mut concentrations = std::collections::HashMap::new();
        concentrations.insert("A".to_string(), 1.0);
        concentrations.insert("B".to_string(), 0.0);

        let kinetics_result = chem_sim
            .simulate_reaction_kinetics(concentrations, 0.1)
            .unwrap();
        assert!(!kinetics_result.is_empty());

        // 3. Materials: Calculate properties at this temperature
        let crystal = CrystalStructure {
            lattice_parameters: [4.0, 4.0, 4.0, 90.0, 90.0, 90.0],
            space_group: "Pm-3m".to_string(),
            atoms: vec![],
        };

        let materials_sim = MaterialsSimulation::new(crystal, temp, 101325.0);
        let thermal_props = materials_sim.calculate_thermal_properties().unwrap();

        assert!(thermal_props.contains_key("heat_capacity"));
        assert!(thermal_props["heat_capacity"] > 0.0);
    }

    #[test]
    fn test_scientific_tensor_physics_integration() {
        // Use scientific tensors in physics simulation

        // Create initial temperature distribution
        let temp_data = ndarray::Array1::linspace(100.0, 200.0, 50);
        let temp_tensor = ScientificTensor::with_units(temp_data.into_dyn(), "Kelvin");

        assert_eq!(temp_tensor.units(), Some("Kelvin"));
        assert_eq!(temp_tensor.shape(), &[50]);

        // Apply some transformations
        let temp_celsius = temp_tensor.mul_scalar(1.0); // Simplified conversion
        let temp_squared = temp_celsius.mul(&temp_celsius).unwrap();

        assert_eq!(temp_squared.shape(), &[50]);

        // Statistical analysis
        let mean_temp = temp_tensor.mean();
        let std_temp = temp_tensor.std();

        assert!(mean_temp > 100.0 && mean_temp < 200.0);
        assert!(std_temp > 0.0);
    }

    #[test]
    fn test_numerical_methods_physics_integration() {
        // Solve physics ODE using numerical methods

        // Simple harmonic oscillator: d²x/dt² + ω²x = 0
        // Convert to system: dx/dt = v, dv/dt = -ω²x
        let omega = 1.0;
        let f = |_t: f64, y: f64| {
            // This is simplified - in reality would need system of ODEs
            -omega * omega * y
        };

        let (t, x) =
            NumericalMethods::runge_kutta_4(f, 1.0, (0.0, 2.0 * std::f64::consts::PI), 1000)
                .unwrap();

        // Should complete one full oscillation
        assert_eq!(t.len(), 1001);
        assert!((t[0] - 0.0).abs() < 1e-10);
        assert!((t[1000] - 2.0 * std::f64::consts::PI).abs() < 0.01);

        // Check that we return close to initial condition after one period
        assert!((x[1000] - x[0]).abs() < 0.1);
    }

    #[test]
    fn test_performance_benchmark() {
        // Performance test for large-scale scientific computing

        let size = 100;

        let start = std::time::Instant::now();

        // Large tensor operations
        let a = ScientificTensor::randn(&[size, size]);
        let b = ScientificTensor::randn(&[size, size]);
        let c = a.add(&b).unwrap();
        let d = a.matmul(&b).unwrap();

        let tensor_time = start.elapsed();

        let start = std::time::Instant::now();

        // Physics simulation
        let mut physics_sim = PhysicsSimulation::new(0.01);
        let _oscillator = physics_sim
            .simulate_harmonic_oscillator(1.0, 1.0, 1.0)
            .unwrap();

        let physics_time = start.elapsed();

        let start = std::time::Instant::now();

        // Numerical methods
        use rand::Rng;
        use rand_distr::{Distribution, StandardNormal};
        let mut rng = rand::thread_rng();
        let data: Vec<f64> = (0..2500).map(|_| StandardNormal.sample(&mut rng)).collect();
        let matrix = ndarray::Array2::from_shape_vec((50, 50), data).unwrap();
        let _eigenvals = NumericalMethods::eigenvalues(&matrix);

        let numerical_time = start.elapsed();

        println!("Performance benchmark:");
        println!("Tensor operations ({}x{}): {:?}", size, size, tensor_time);
        println!("Physics simulation: {:?}", physics_time);
        println!("Numerical methods (50x50 eigenvals): {:?}", numerical_time);

        // Basic performance expectations
        assert!(tensor_time.as_millis() < 1000);
        assert!(physics_time.as_millis() < 1000);
        assert!(numerical_time.as_millis() < 5000);
    }
}