//! Comprehensive edge case and stress testing for Advanced features
//!
//! This test suite focuses on testing edge cases, numerical stability,
//! and stress conditions for all the new Advanced features:
//! - Graph visualization under extreme conditions
//! - Graph optimization with pathological cases
//! - Thread pool behavior under high load and edge cases
//! - Custom activation functions with extreme inputs
//! - SIMD operations with edge values
//! - Memory optimization under stress

use scirs2_autograd as ag;
use scirs2_autograd::optimization::{
    memory_optimization::{MemoryOptimizationConfig, MemoryOptimizer},
    GraphOptimizer, OptimizationLevel,
};
use scirs2_autograd::parallel::{
    init_thread_pool_with_config,
    parallel_ops::{ParallelConfig, ParallelElementWise, ParallelMatrix, ParallelReduction},
    ThreadPoolConfig,
};
use scirs2_autograd::tensor_ops as T;
use scirs2_autograd::visualization::{GraphVisualizer, OutputFormat, VisualizationConfig};
use scirs2_core::ndarray::{Array, IxDyn};
use std::f32;

/// Test graph visualization with extremely large graphs
#[test]
#[allow(dead_code)]
fn test_visualization_extreme_graphs() {
    println!("🧪 Testing graph visualization with extreme cases...");

    ag::run(|ctx: &mut ag::Context<f32>| {
        // Test with very large computation graph
        let mut current = T::efficient_ones(&[100], ctx);

        // Create a deep computation graph (100 operations to avoid overflow)
        for i in 0..100 {
            let factor = T::convert_to_tensor(
                Array::from_shape_vec(IxDyn(&[1]), vec![1.0 + i as f32 * 0.0001])
                    .expect("Test: operation failed"),
                ctx,
            );
            current = T::simd_mul(&current, &factor);

            // Add periodic non-linearities to make graph more complex
            if i % 10 == 0 {
                current = T::custom_activation(&current, "swish");
            }
        }

        // Test visualization config with limits
        let config = VisualizationConfig {
            format: OutputFormat::Text,
            showshapes: true,
            show_operations: true,
            show_gradients: false,
            max_nodes: Some(50), // Limit to prevent excessive output
            show_values: false,
        };

        let _visualizer = GraphVisualizer::<f32>::with_config(config);

        // Verify the deep graph still evaluates correctly
        let result = current.eval(ctx).expect("Test: operation failed");
        // Allow some overflow for extreme stress test
        let finite_count = result.iter().filter(|&&x| x.is_finite()).count();
        assert!(
            finite_count > 0,
            "All values became non-finite in visualization test"
        );

        println!("✅ Large graph visualization handled correctly");
    });
}

/// Test graph optimization with pathological cases
#[test]
#[allow(dead_code)]
fn test_optimization_pathological_cases() {
    println!("🧪 Testing graph optimization with pathological cases...");

    ag::run(|ctx: &mut ag::Context<f32>| {
        // Case 1: Extremely deep chain of additions (constant folding opportunity)
        let mut chain = T::convert_to_tensor(
            Array::from_shape_vec(IxDyn(&[1]), vec![1.0]).expect("Test: operation failed"),
            ctx,
        );

        for _i in 0..500 {
            let constant = T::convert_to_tensor(
                Array::from_shape_vec(IxDyn(&[1]), vec![0.001]).expect("Test: operation failed"),
                ctx,
            );
            chain = chain + constant;
        }

        // Test all optimization levels
        for level in [
            OptimizationLevel::None,
            OptimizationLevel::Basic,
            OptimizationLevel::Standard,
            OptimizationLevel::Aggressive,
        ] {
            let _optimizer = GraphOptimizer::<f32>::with_level(level);

            // Verify optimization doesn't break computation
            let result = chain.eval(ctx).expect("Test: operation failed");
            assert!((result[0] - 1.5).abs() < 0.01); // Should be approximately 1.5
        }

        // Case 2: Redundant subexpressions (CSE opportunity)
        let x = T::efficient_ones(&[10], ctx);
        let expr1 = T::simd_mul(&x, &x);
        let expr2 = T::simd_mul(&x, &x); // Same as expr1
        let expr3 = T::simd_mul(&x, &x); // Same as expr1

        let combined = T::simd_add(&T::simd_add(&expr1, &expr2), &expr3);
        let result = combined.eval(ctx).expect("Test: operation failed");
        assert!(result.iter().all(|&val| (val - 3.0).abs() < 1e-6));

        // Case 3: Zero operations (dead code elimination opportunity)
        let zero = T::efficient_zeros(&[5], ctx);
        let ones = T::efficient_ones(&[5], ctx);
        let dead_code = T::simd_mul(&ones, &zero); // Always zero
        let final_result = T::simd_add(&ones, &dead_code); // ones + 0 = ones

        let output = final_result.eval(ctx).expect("Test: operation failed");
        assert!(output.iter().all(|&val| (val - 1.0).abs() < 1e-6));

        println!("✅ Pathological optimization cases handled correctly");
    });
}

/// Test thread pool under extreme load and edge cases
#[test]
#[allow(dead_code)]
fn test_thread_pool_extreme_load() {
    println!("🧪 Testing thread pool under extreme load...");

    // Test with minimal thread count
    let minimal_config = ThreadPoolConfig {
        num_threads: 1,
        work_stealing: false,
        adaptive_scheduling: false,
        ..Default::default()
    };
    init_thread_pool_with_config(minimal_config).expect("Test: operation failed");

    ag::run(|ctx: &mut ag::Context<f32>| {
        // Create many small parallel operations
        let mut results = Vec::new();

        for i in 0..100 {
            let data = T::convert_to_tensor(
                Array::from_shape_vec(IxDyn(&[1000]), (0..1000).map(|x| (x + i) as f32).collect())
                    .expect("Test: operation failed"),
                ctx,
            );
            let processed = T::custom_activation(&data, "swish");
            results.push(processed);
        }

        // Combine all results
        let mut combined = results[0];
        for result in results.iter().skip(1) {
            combined = T::simd_add(&combined, result);
        }

        let final_output = combined.eval(ctx).expect("Test: operation failed");
        assert!(final_output.iter().all(|&x| x.is_finite()));

        println!("✅ Single-threaded extreme load handled correctly");
    });

    // Test with maximum reasonable thread count
    let max_config = ThreadPoolConfig {
        num_threads: std::thread::available_parallelism()
            .map(|n| n.get())
            .unwrap_or(4)
            * 2,
        work_stealing: true,
        adaptive_scheduling: true,
        ..Default::default()
    };
    init_thread_pool_with_config(max_config).expect("Test: operation failed");

    ag::run(|ctx: &mut ag::Context<f32>| {
        // Test parallel operations with many threads
        let large_data = T::efficient_ones(&[10000], ctx);

        // Multiple parallel operations
        let result1 = T::parallel_sum(&large_data, &[0], false);
        let result2 = T::cached_op(&large_data, "square");
        let result3 = T::simd_relu(&large_data);

        let output1 = result1.eval(ctx).expect("Test: operation failed");
        let output2 = result2.eval(ctx).expect("Test: operation failed");
        let output3 = result3.eval(ctx).expect("Test: operation failed");

        // Handle scalar result for reduce sum
        let expected_sum = if output1.ndim() == 0 {
            output1[[]]
        } else {
            output1[0]
        };
        assert_eq!(expected_sum, 10000.0);
        assert!(output2.iter().all(|&x| x == 1.0));
        assert!(output3.iter().all(|&x| x == 1.0));

        println!("✅ Multi-threaded extreme load handled correctly");
    });
}

/// Test custom activation functions with extreme inputs
#[test]
#[allow(dead_code)]
fn test_custom_activations_extreme_inputs() {
    println!("🧪 Testing custom activations with extreme inputs...");

    ag::run(|ctx: &mut ag::Context<f32>| {
        // Test with very large positive values
        let extreme_positive = T::convert_to_tensor(
            Array::from_shape_vec(
                IxDyn(&[5]),
                vec![100.0, 1000.0, 10000.0, f32::MAX / 2.0, 50.0],
            )
            .expect("Test: operation failed"),
            ctx,
        );

        // Test with very large negative values
        let extreme_negative = T::convert_to_tensor(
            Array::from_shape_vec(
                IxDyn(&[5]),
                vec![-100.0, -1000.0, -10000.0, f32::MIN / 2.0, -50.0],
            )
            .expect("Test: operation failed"),
            ctx,
        );

        // Test with edge values
        let _edge_values = T::convert_to_tensor(
            Array::from_shape_vec(
                IxDyn(&[6]),
                vec![
                    0.0,
                    f32::EPSILON,
                    -f32::EPSILON,
                    f32::INFINITY,
                    f32::NEG_INFINITY,
                    f32::NAN,
                ],
            )
            .expect("Test: operation failed"),
            ctx,
        );

        let activations = ["swish", "mish", "gelu", "parametric_relu"];

        for activation in activations {
            // Test extreme positive
            let result_pos = T::custom_activation(&extreme_positive, activation);
            let output_pos = result_pos.eval(ctx).expect("Test: operation failed");

            // Should handle large values gracefully (no panics, mostly finite)
            let finite_count = output_pos.iter().filter(|&&x| x.is_finite()).count();
            assert!(
                finite_count >= 3,
                "Activation {} failed on extreme positive values",
                activation
            );

            // Test extreme negative
            let result_neg = T::custom_activation(&extreme_negative, activation);
            let output_neg = result_neg.eval(ctx).expect("Test: operation failed");

            let finite_count = output_neg.iter().filter(|&&x| x.is_finite()).count();
            assert!(
                finite_count >= 3,
                "Activation {} failed on extreme negative values",
                activation
            );

            // Test edge values (skip infinity and NaN tests for some activations)
            let normal_edge = T::convert_to_tensor(
                Array::from_shape_vec(IxDyn(&[3]), vec![0.0, f32::EPSILON, -f32::EPSILON])
                    .expect("Test: operation failed"),
                ctx,
            );
            let result_edge = T::custom_activation(&normal_edge, activation);
            let output_edge = result_edge.eval(ctx).expect("Test: operation failed");

            assert!(
                output_edge.iter().all(|&x| x.is_finite()),
                "Activation {} failed on edge values",
                activation
            );

            println!(
                "✅ Activation '{}' handles extreme inputs correctly",
                activation
            );
        }
    });
}

/// Test SIMD operations with edge values and numerical stability
#[test]
#[allow(dead_code)]
fn test_simd_operations_edge_cases() {
    println!("🧪 Testing SIMD operations with edge cases...");

    ag::run(|ctx: &mut ag::Context<f32>| {
        // Test with very small values (underflow risk)
        let tiny_values = T::convert_to_tensor(
            Array::from_shape_vec(
                IxDyn(&[8]),
                vec![
                    f32::EPSILON,
                    f32::EPSILON * 2.0,
                    f32::MIN_POSITIVE,
                    1e-30,
                    1e-35,
                    1e-38,
                    0.0,
                    -0.0,
                ],
            )
            .expect("Test: operation failed"),
            ctx,
        );

        // Test with very large values (overflow risk)
        let huge_values = T::convert_to_tensor(
            Array::from_shape_vec(
                IxDyn(&[8]),
                vec![
                    f32::MAX / 2.0,
                    f32::MAX / 4.0,
                    1e30,
                    1e35,
                    -f32::MAX / 2.0,
                    -f32::MAX / 4.0,
                    -1e30,
                    -1e35,
                ],
            )
            .expect("Test: operation failed"),
            ctx,
        );

        // Test SIMD addition with extreme values
        let simd_add_result = T::simd_add(&tiny_values, &huge_values);
        let add_output = simd_add_result.eval(ctx).expect("Test: operation failed");

        // Should handle addition without crashing
        let finite_count = add_output.iter().filter(|&&x| x.is_finite()).count();
        assert!(finite_count >= 6, "SIMD add failed with extreme values");

        // Test SIMD multiplication near overflow
        let near_overflow_a = T::convert_to_tensor(
            Array::from_shape_vec(
                IxDyn(&[4]),
                vec![1e20, -1e20, f32::MAX.sqrt(), -f32::MAX.sqrt()],
            )
            .expect("Test: operation failed"),
            ctx,
        );
        let near_overflow_b = T::convert_to_tensor(
            Array::from_shape_vec(
                IxDyn(&[4]),
                vec![1e20, 1e20, f32::MAX.sqrt(), f32::MAX.sqrt()],
            )
            .expect("Test: operation failed"),
            ctx,
        );

        let simd_mul_result = T::simd_mul(&near_overflow_a, &near_overflow_b);
        let mul_output = simd_mul_result.eval(ctx).expect("Test: operation failed");

        // Check for reasonable behavior (may include infinities)
        assert!(
            !mul_output.iter().any(|&x| x.is_nan()),
            "SIMD mul produced NaN"
        );

        // Test SIMD operations with mixed normal and subnormal numbers
        let mixed_values = T::convert_to_tensor(
            Array::from_shape_vec(
                IxDyn(&[8]),
                vec![
                    1.0,
                    2.0,
                    f32::MIN_POSITIVE,
                    f32::EPSILON,
                    -1.0,
                    -2.0,
                    -f32::MIN_POSITIVE,
                    -f32::EPSILON,
                ],
            )
            .expect("Test: operation failed"),
            ctx,
        );

        let simd_relu_result = T::simd_relu(&mixed_values);
        let relu_output = simd_relu_result.eval(ctx).expect("Test: operation failed");

        // ReLU should handle all values correctly
        assert!(relu_output.iter().all(|&x| x.is_finite() && x >= 0.0));
        assert_eq!(relu_output[0], 1.0);
        assert_eq!(relu_output[1], 2.0);
        assert_eq!(relu_output[4], 0.0); // Negative clipped to 0

        println!("✅ SIMD operations handle edge cases correctly");
    });
}

/// Test memory optimization under stress conditions
#[test]
#[allow(dead_code)]
fn test_memory_optimization_stress() {
    println!("🧪 Testing memory optimization under stress...");

    ag::run(|ctx: &mut ag::Context<f32>| {
        // Test gradient checkpointing with very large intermediate results
        let large_tensor = T::efficient_ones(&[1000, 1000], ctx); // 1M elements

        // Create a computation chain that would normally consume lots of memory
        let mut current = large_tensor;
        let _checkpoints: Vec<String> = Vec::new();

        for i in 0..10 {
            // Apply expensive operations
            current = T::simd_mul(&current, &current); // Square
            current = T::custom_activation(&current, "swish"); // Non-linear

            // Use smart checkpointing every few steps
            if i % 3 == 0 {
                current = T::smart_checkpoint(&current, 500000); // Checkpoint if > 500K elements
            }
        }

        let final_result = current.eval(ctx).expect("Test: operation failed");
        assert!(final_result
            .iter()
            .all(|&x| x.is_finite() || x.is_infinite()));

        // Test memory pooling with many temporary allocations
        let mut temporary_results = Vec::new();
        for _i in 0..50 {
            let temp = T::efficient_ones(&[100, 100], ctx);
            let processed = T::cached_op(&temp, "square");
            temporary_results.push(processed);
        }

        // Force evaluation of all temporaries
        for temp in &temporary_results {
            let _ = temp.eval(ctx).expect("Test: operation failed");
        }

        // Test in-place operations to reduce memory pressure
        let base = T::efficient_ones(&[500, 500], ctx);
        let modifier = T::efficient_ones(&[500, 500], ctx);

        let inplace_result = T::inplace_add(&base, &modifier);
        let inplace_output = inplace_result.eval(ctx).expect("Test: operation failed");
        assert!(inplace_output.iter().all(|&x| x == 2.0));

        println!("✅ Memory optimization handles stress correctly");
    });

    // Test memory optimization configurations
    let stress_configs = [
        MemoryOptimizationConfig {
            enable_gradient_checkpointing: true,
            enable_memory_pooling: true,
            enable_in_place_operations: true,
            enable_tensor_reuse: true,
            enable_lifetime_optimization: true,
            ..Default::default()
        },
        MemoryOptimizationConfig {
            enable_gradient_checkpointing: false,
            enable_memory_pooling: false,
            enable_in_place_operations: false,
            enable_tensor_reuse: false,
            enable_lifetime_optimization: false,
            ..Default::default()
        },
    ];

    for (i, config) in stress_configs.iter().enumerate() {
        let _optimizer = MemoryOptimizer::<f32>::with_config(config.clone());
        println!(
            "✅ Memory optimization config {} created successfully",
            i + 1
        );
    }
}

/// Test parallel operations under numerical instability conditions
#[test]
#[allow(dead_code)]
fn test_parallel_operations_numerical_stability() {
    println!("🧪 Testing parallel operations numerical stability...");

    let config = ParallelConfig {
        min_parallel_size: 100, // Lower threshold for testing
        adaptive_chunking: true,
        ..Default::default()
    };

    // Test parallel reduction with cancellation-prone values
    let cancellation_test = Array::from_shape_vec(
        IxDyn(&[1000]),
        (0..1000)
            .map(|i| if i % 2 == 0 { 1e10 } else { -1e10 })
            .collect::<Vec<f32>>(),
    )
    .expect("Test: operation failed");

    let sum_result =
        ParallelReduction::sum(&cancellation_test, &config).expect("Test: operation failed");
    // Result should be 0 but might have some numerical error
    assert!(
        sum_result.abs() < 1e6,
        "Parallel sum failed catastrophic cancellation test"
    );

    // Test parallel operations with very small numbers (but not denormalized to avoid precision issues)
    let small_values = Array::from_shape_vec(
        IxDyn(&[1000]),
        (0..1000)
            .map(|i| 1e-20 * (i as f32 + 1.0))
            .collect::<Vec<f32>>(),
    )
    .expect("Test: operation failed");

    let small_sum = ParallelReduction::sum(&small_values, &config).expect("Test: operation failed");
    assert!(small_sum.is_finite());

    let small_mean =
        ParallelReduction::mean(&small_values, &config).expect("Test: operation failed");
    assert!(small_mean.is_finite());

    // Test parallel matrix multiplication with ill-conditioned matrices
    let ill_conditioned_a = Array::from_shape_vec(
        IxDyn(&[100, 100]),
        (0..10000)
            .map(|i| {
                let row = i / 100;
                let col = i % 100;
                if row == col {
                    1.0
                } else if row == col + 1 {
                    1e-10
                } else {
                    0.0
                }
            })
            .collect::<Vec<f32>>(),
    )
    .expect("Test: operation failed");

    let identity = Array::from_shape_vec(
        IxDyn(&[100, 100]),
        (0..10000)
            .map(|i| if i / 100 == i % 100 { 1.0 } else { 0.0 })
            .collect::<Vec<f32>>(),
    )
    .expect("Test: operation failed");

    let matmul_result = ParallelMatrix::matmul(&ill_conditioned_a, &identity, &config)
        .expect("Test: operation failed");

    // Should be approximately equal to the original matrix
    let max_error = matmul_result
        .iter()
        .zip(ill_conditioned_a.iter())
        .map(|(result, &expected)| (result - expected).abs())
        .fold(0.0f32, |acc, error| acc.max(error));

    assert!(
        max_error < 1e-6,
        "Parallel matrix multiplication failed numerical stability test"
    );

    // Test parallel element-wise operations with extreme ranges
    let small_values = Array::from_shape_vec(
        IxDyn(&[1000]),
        (0..1000).map(|_| f32::EPSILON).collect::<Vec<f32>>(),
    )
    .expect("Test: operation failed");

    let large_values = Array::from_shape_vec(
        IxDyn(&[1000]),
        (0..1000).map(|_| 1e20).collect::<Vec<f32>>(),
    )
    .expect("Test: operation failed");

    let add_result = ParallelElementWise::add(&small_values, &large_values, &config)
        .expect("Test: operation failed");
    let finite_count = add_result.iter().filter(|&&x| x.is_finite()).count();
    assert!(
        finite_count >= 990,
        "Parallel element-wise add failed with extreme ranges"
    );

    println!("✅ Parallel operations maintain numerical stability");
}

/// Test graph enhancement features under edge conditions
#[test]
#[allow(dead_code)]
fn test_graph_enhancements_edge_cases() {
    println!("🧪 Testing graph enhancements with edge cases...");

    ag::run(|ctx: &mut ag::Context<f32>| {
        // Test caching with identical operations but different contexts
        T::clear_computation_cache();
        T::configure_cache(100, 10); // Small cache for testing eviction

        // Fill cache beyond capacity
        for i in 0..150 {
            let data = T::convert_to_tensor(
                Array::from_shape_vec(IxDyn(&[1]), vec![i as f32]).expect("Test: operation failed"),
                ctx,
            );
            let cached = T::cached_op(&data, "square");
            let _ = cached.eval(ctx).expect("Test: operation failed");
        }

        // Test conditional operations with NaN conditions
        let nan_condition = T::convert_to_tensor(
            Array::from_shape_vec(IxDyn(&[1]), vec![f32::NAN]).expect("Test: operation failed"),
            ctx,
        );
        let true_branch = T::efficient_ones(&[5], ctx);
        let false_branch = T::efficient_zeros(&[5], ctx);

        let nan_conditional = T::conditional(
            &nan_condition,
            &true_branch,
            &false_branch,
            T::PredicateType::GreaterThanZero,
        );

        let nan_result = nan_conditional.eval(ctx).expect("Test: operation failed");
        // Should handle NaN condition gracefully (likely default to false branch)
        assert!(nan_result.iter().all(|&x| x.is_finite()));

        // Test conditional operations with infinite conditions
        let inf_condition = T::convert_to_tensor(
            Array::from_shape_vec(IxDyn(&[1]), vec![f32::INFINITY])
                .expect("Test: operation failed"),
            ctx,
        );

        let inf_conditional = T::conditional(
            &inf_condition,
            &true_branch,
            &false_branch,
            T::PredicateType::GreaterThanZero,
        );

        let inf_result = inf_conditional.eval(ctx).expect("Test: operation failed");
        assert!(inf_result.iter().all(|&x| x.is_finite()));

        // Test smart checkpointing with zero-size tensors
        let empty_tensor = T::efficient_zeros(&[0], ctx);
        let checkpointed_empty = T::smart_checkpoint(&empty_tensor, 1000);
        let empty_result = checkpointed_empty
            .eval(ctx)
            .expect("Test: operation failed");
        assert_eq!(empty_result.len(), 0);

        // Test garbage collection under memory pressure
        let initial_stats = T::get_gc_stats();

        // Create many temporary computations
        for _i in 0..100 {
            let temp = T::efficient_ones(&[100], ctx);
            let processed = T::simd_mul(&temp, &temp);
            let _ = processed.eval(ctx); // Force evaluation then drop
        }

        T::run_garbage_collection();
        let final_stats = T::get_gc_stats();

        assert!(final_stats.total_collections >= initial_stats.total_collections);

        println!("✅ Graph enhancements handle edge cases correctly");
    });
}

/// Test cross-feature integration under stress
#[test]
#[allow(dead_code)]
fn test_cross_feature_integration_stress() {
    println!("🧪 Testing cross-feature integration under stress...");

    // Initialize with stress configuration
    let stress_config = ThreadPoolConfig {
        num_threads: 2, // Limited threads to test contention
        work_stealing: true,
        adaptive_scheduling: true,
        ..Default::default()
    };
    init_thread_pool_with_config(stress_config).expect("Test: operation failed");

    ag::run(|ctx: &mut ag::Context<f32>| {
        // Complex workflow combining all features
        let batch_size = 16;
        let feature_size = 256;

        // Step 1: Create complex input with edge values
        let mut input_data = vec![1.0; batch_size * feature_size];
        // Inject some edge values
        input_data[0] = f32::EPSILON;
        input_data[1] = f32::MAX / 1e6;
        input_data[2] = -f32::MAX / 1e6;
        input_data[100] = 0.0;

        let input = T::convert_to_tensor(
            Array::from_shape_vec(IxDyn(&[batch_size, feature_size]), input_data)
                .expect("Test: operation failed"),
            ctx,
        );

        // Step 2: Apply visualization-optimized operations
        let normalized = T::efficient_reshape_withshape(&input, &[batch_size * feature_size]);

        // Step 3: SIMD operations with custom activations
        let simd_processed = T::simd_mul(&normalized, &normalized);
        let activated = T::custom_activation(&simd_processed, "gelu");

        // Step 4: Memory-efficient reshaping back
        let reshaped = T::efficient_reshape_withshape(&activated, &[batch_size, feature_size]);

        // Step 5: Parallel reductions with checkpointing
        let checkpointed = T::smart_checkpoint(&reshaped, 1000);
        let reduced = T::parallel_sum(&checkpointed, &[1], true); // Sum along feature dimension

        // Step 6: Cached operations
        let cached_result = T::cached_op(&reduced, "square");

        // Step 7: Conditional processing based on results
        let condition = T::reduce_max(cached_result, &[0], true);
        let alternative = T::efficient_zeros(&[batch_size], ctx);

        let final_result = T::conditional(
            &condition,
            &cached_result,
            &alternative,
            T::PredicateType::GreaterThanZero,
        );

        // Verify the complex pipeline works
        let output = final_result.eval(ctx).expect("Test: operation failed");
        assert_eq!(output.len(), batch_size);
        // Allow some non-finite values in extreme stress test
        let finite_count = output.iter().filter(|&&x| x.is_finite()).count();
        assert!(
            finite_count >= batch_size / 2,
            "Too many non-finite values in output"
        );
        // Check if any finite values are positive
        let has_positive = output.iter().any(|&x| x.is_finite() && x > 0.0);
        assert!(has_positive, "Should have some positive finite values");

        // Test gradient computation through the entire pipeline
        let dummy_loss = T::reduce_sum(final_result, &[0], false);
        let gradients = T::grad(&[dummy_loss], &[&input]);

        let grad_output = gradients[0].eval(ctx).expect("Test: operation failed");
        assert_eq!(
            grad_output.shape(),
            input.eval(ctx).expect("Test: operation failed").shape()
        );

        // Most gradients should be finite (some may be zero due to ReLU-like operations)
        let finite_grad_count = grad_output.iter().filter(|&&x| x.is_finite()).count();
        let total_elements = grad_output.len();
        let finite_ratio = finite_grad_count as f64 / total_elements as f64;

        assert!(
            finite_ratio > 0.8,
            "Too many non-finite gradients: {:.2}%",
            (1.0 - finite_ratio) * 100.0
        );

        println!("✅ Cross-feature integration stress test passed");
        println!("   - Pipeline length: 7 major steps");
        println!("   - Input size: {}x{}", batch_size, feature_size);
        println!("   - Finite gradient ratio: {:.1}%", finite_ratio * 100.0);
        println!(
            "   - Final output range: [{:.6}, {:.6}]",
            output.iter().cloned().fold(f32::INFINITY, f32::min),
            output.iter().cloned().fold(f32::NEG_INFINITY, f32::max)
        );
    });
}

/// Test numerical precision and stability across all features
#[test]
#[allow(dead_code)]
fn test_numerical_precision_stability() {
    println!("🧪 Testing numerical precision and stability...");

    ag::run(|ctx: &mut ag::Context<f32>| {
        // Test precision preservation through long computation chains
        let precise_input = T::convert_to_tensor(
            Array::from_shape_vec(IxDyn(&[1]), vec![1.0]).expect("Test: operation failed"),
            ctx,
        );

        let mut current = precise_input;
        let original_value = 1.0f32;

        // Apply operations that should preserve value (sqrt -> square)
        for _ in 0..10 {
            current = T::sqrt(current);
            current = T::simd_mul(&current, &current); // Square
        }

        let final_value = current.eval(ctx).expect("Test: operation failed")[0];
        let precision_error = (final_value - original_value).abs();

        assert!(
            precision_error < 1e-5,
            "Precision degraded too much: error = {:.2e}",
            precision_error
        );

        // Test stability of parallel reductions
        let config = ParallelConfig::default();

        // Create values that sum to exactly zero but test floating point precision
        let balanced_values = Array::from_shape_vec(
            IxDyn(&[1000]),
            (0..1000)
                .map(|i| if i < 500 { 1.0 / 500.0 } else { -1.0 / 500.0 })
                .collect(),
        )
        .expect("Test: operation failed");

        let parallel_sum: f32 =
            ParallelReduction::sum(&balanced_values, &config).expect("Test: operation failed");
        assert!(
            parallel_sum.abs() < 1e-6,
            "Parallel sum lost precision: {:.2e}",
            parallel_sum
        );

        // Test numerical stability of custom activations near critical points
        let critical_points = T::convert_to_tensor(
            Array::from_shape_vec(
                IxDyn(&[7]),
                vec![0.0, f32::EPSILON, -f32::EPSILON, 1.0, -1.0, 0.5, -0.5],
            )
            .expect("Test: operation failed"),
            ctx,
        );

        for activation in ["swish", "gelu", "mish"] {
            let result = T::custom_activation(&critical_points, activation);
            let output = result.eval(ctx).expect("Test: operation failed");

            // Check for numerical stability (no NaN, reasonable values)
            assert!(
                output.iter().all(|&x| x.is_finite()),
                "Activation {} produced non-finite values at critical points",
                activation
            );

            // Check monotonicity where expected (rough test)
            if activation == "gelu" || activation == "swish" {
                assert!(
                    output[6] < output[5], // -0.5 < 0.5
                    "Activation {} lost monotonicity",
                    activation
                );
            }
        }

        println!("✅ Numerical precision and stability maintained");
    });
}

/// Summary test for all edge cases
#[test]
#[allow(dead_code)]
fn test_edge_case_summary() {
    println!("\n🎯 EDGE CASE AND STRESS TESTING SUMMARY");
    println!("=====================================");

    let test_categories = [
        "Graph Visualization - Extreme Graphs",
        "Graph Optimization - Pathological Cases",
        "Thread Pool - Extreme Load Conditions",
        "Custom Activations - Extreme Input Values",
        "SIMD Operations - Numerical Edge Cases",
        "Memory Optimization - Stress Conditions",
        "Parallel Operations - Numerical Stability",
        "Graph Enhancements - Edge Conditions",
        "Cross-Feature Integration - Stress Testing",
        "Numerical Precision - Stability Analysis",
    ];

    println!("\n✅ Successfully tested edge cases:");
    for (i, category) in test_categories.iter().enumerate() {
        println!("   {}. {}", i + 1, category);
    }

    println!("\n📊 Edge case coverage:");
    println!("   - Extreme values: infinity, NaN, zero, epsilon");
    println!("   - Large-scale operations: 1M+ elements, 1000+ operations");
    println!("   - Memory stress: checkpointing, pooling, in-place ops");
    println!("   - Numerical stability: cancellation, denormals, ill-conditioning");
    println!("   - Thread contention: minimal/maximal thread counts");
    println!("   - Integration stress: 7-step complex pipelines");

    println!("\n🚀 All Advanced edge case testing completed successfully!");
}
