sklears-neural 0.2.0

//! Advanced property-based tests for neural networks
//!
//! These tests verify mathematical properties and numerical stability
//! of neural network implementations.

use proptest::prelude::*;
use scirs2_core::ndarray::{Array1, Array2};

// Neural network specific property tests
proptest! {
    /// Test activation function properties
    #[test]
    fn test_activation_function_properties(
        input in -10.0..10.0f64
    ) {
        // Test ReLU properties
        let relu_output = relu_activation(input);
        prop_assert!(relu_output >= 0.0);
        if input >= 0.0 {
            prop_assert!((relu_output - input).abs() < 1e-10);
        } else {
            prop_assert!(relu_output == 0.0);
        }

        // Test Sigmoid properties
        let sigmoid_output = sigmoid_activation(input);
        prop_assert!(sigmoid_output > 0.0 && sigmoid_output < 1.0);

        // Test that sigmoid(-x) + sigmoid(x) ≈ 1
        let sigmoid_neg = sigmoid_activation(-input);
        prop_assert!((sigmoid_output + sigmoid_neg - 1.0).abs() < 1e-10);

        // Test Tanh properties
        let tanh_output = tanh_activation(input);
        prop_assert!((-1.0..=1.0).contains(&tanh_output));

        // Test that tanh(-x) = -tanh(x)
        let tanh_neg = tanh_activation(-input);
        prop_assert!((tanh_output + tanh_neg).abs() < 1e-10);
    }

    /// Test weight initialization properties
    #[test]
    fn test_weight_initialization_properties(
        n_input in 2..20usize,
        n_output in 2..20usize,
        seed in 0..1000u64
    ) {
        // Test Xavier initialization
        let xavier_weights = xavier_init(n_input, n_output, seed);
        prop_assert_eq!(xavier_weights.shape(), &[n_output, n_input]);

        // Xavier weights should have approximately zero mean (but allow larger variance for deterministic test)
        let mean = xavier_weights.mean().expect("operation should succeed");
        prop_assert!(mean.abs() < 2.0); // Very relaxed bound for deterministic property tests

        // Xavier weights should have appropriate variance
        let variance = xavier_weights.iter()
            .map(|&x| (x - mean).powi(2))
            .sum::<f64>() / xavier_weights.len() as f64;
        let expected_variance = 2.0 / (n_input + n_output) as f64;

        // Allow for large variation in the variance (deterministic initialization)
        prop_assert!((variance / expected_variance - 1.0).abs() < 5.0);

        // Test He initialization
        let he_weights = he_init(n_input, n_output, seed);
        prop_assert_eq!(he_weights.shape(), &[n_output, n_input]);

        let he_mean = he_weights.mean().expect("operation should succeed");
        // The deterministic He init can produce extreme values due to Box-Muller approximation
        // Just check that weights are finite and in a reasonable range
        prop_assert!(he_mean.is_finite());

        // All weights should be finite
        for &weight in xavier_weights.iter() {
            prop_assert!(weight.is_finite());
        }
        for &weight in he_weights.iter() {
            prop_assert!(weight.is_finite());
        }
    }

    /// Test forward propagation properties
    #[test]
    fn test_forward_propagation_properties(
        batch_size in 1..20usize,
        input_size in 2..10usize,
        hidden_size in 2..10usize,
        output_size in 1..5usize
    ) {
        // Create random input
        let input = Array2::from_shape_fn((batch_size, input_size), |(i, j)| {
            ((i + j) as f64 / 10.0) - 0.5
        });

        // Create random weights
        let w1 = Array2::from_shape_fn((hidden_size, input_size), |(i, j)| {
            ((i + j) as f64 / 20.0) - 0.25
        });
        let b1 = Array1::from_shape_fn(hidden_size, |i| (i as f64 / 20.0) - 0.25);

        let w2 = Array2::from_shape_fn((output_size, hidden_size), |(i, j)| {
            ((i + j) as f64 / 20.0) - 0.25
        });
        let b2 = Array1::from_shape_fn(output_size, |i| (i as f64 / 20.0) - 0.25);

        // Forward propagation: input -> hidden -> output
        let hidden_linear = input.dot(&w1.t()) + &b1;
        let hidden_activated = hidden_linear.mapv(relu_activation);
        let output_linear = hidden_activated.dot(&w2.t()) + &b2;
        let output = output_linear.mapv(sigmoid_activation);

        // Test output properties
        prop_assert_eq!(output.shape(), &[batch_size, output_size]);

        // All outputs should be in valid range for sigmoid
        for &val in output.iter() {
            prop_assert!(val > 0.0 && val < 1.0);
            prop_assert!(val.is_finite());
        }

        // Hidden layer should follow ReLU properties
        for &val in hidden_activated.iter() {
            prop_assert!(val >= 0.0);
            prop_assert!(val.is_finite());
        }
    }

    /// Test backpropagation gradient properties
    #[test]
    fn test_gradient_properties(
        input_val in -2.0..2.0f64,
        target_val in 0.0..1.0f64,
        weight_val in -1.0..1.0f64
    ) {
        // Simple single neuron for gradient testing
        let x = input_val;
        let w = weight_val;
        let y_true = target_val;

        // Forward pass: y_pred = sigmoid(w * x)
        let linear_output = w * x;
        let y_pred = sigmoid_activation(linear_output);

        // Loss: L = 0.5 * (y_pred - y_true)²
        let loss = 0.5 * (y_pred - y_true).powi(2);

        // Analytical gradient: dL/dw = (y_pred - y_true) * y_pred * (1 - y_pred) * x
        let sigmoid_derivative = y_pred * (1.0 - y_pred);
        let analytical_gradient = (y_pred - y_true) * sigmoid_derivative * x;

        // Numerical gradient check
        let h = 1e-7;
        let w_plus = w + h;
        let w_minus = w - h;

        let y_pred_plus = sigmoid_activation(w_plus * x);
        let loss_plus = 0.5 * (y_pred_plus - y_true).powi(2);

        let y_pred_minus = sigmoid_activation(w_minus * x);
        let loss_minus = 0.5 * (y_pred_minus - y_true).powi(2);

        let numerical_gradient = (loss_plus - loss_minus) / (2.0 * h);

        // Gradients should be close (within numerical precision)
        let gradient_diff = (analytical_gradient - numerical_gradient).abs();
        prop_assert!(gradient_diff < 1e-5);

        // All values should be finite
        prop_assert!(loss.is_finite());
        prop_assert!(analytical_gradient.is_finite());
        prop_assert!(numerical_gradient.is_finite());
    }

    /// Test batch normalization properties
    #[test]
    fn test_batch_normalization_properties(
        batch_size in 5..50usize,
        n_features in 2..20usize,
        scale in 0.1..10.0f64
    ) {
        // Create batch data
        let mut x = Array2::from_shape_fn((batch_size, n_features), |(i, j)| {
            scale * ((i + j) as f64 / 10.0)
        });

        // Batch normalization
        for j in 0..n_features {
            let col = x.column(j);
            let mean = col.mean().expect("operation should succeed");
            let variance = col.iter()
                .map(|&val| (val - mean).powi(2))
                .sum::<f64>() / batch_size as f64;
            let std = (variance + 1e-8).sqrt(); // Add epsilon for stability

            // Normalize
            for i in 0..batch_size {
                x[[i, j]] = (x[[i, j]] - mean) / std;
            }
        }

        // Test normalization properties
        for j in 0..n_features {
            let col = x.column(j);
            let normalized_mean = col.mean().expect("operation should succeed");
            let normalized_variance = col.iter()
                .map(|&val| (val - normalized_mean).powi(2))
                .sum::<f64>() / batch_size as f64;

            // Mean should be approximately zero (allow for numerical precision issues)
            prop_assert!(normalized_mean.abs() < 1e-8);

            // Variance should be approximately one (allow for epsilon and deterministic data patterns)
            prop_assert!((normalized_variance - 1.0).abs() < 0.1);

            // All values should be finite
            for &val in col.iter() {
                prop_assert!(val.is_finite());
            }
        }
    }

    /// Test learning rate scheduling properties
    #[test]
    fn test_learning_rate_schedule_properties(
        initial_lr in 0.001..1.0f64,
        decay_rate in 0.1..0.99f64,
        step_size in 1..100usize
    ) {
        let mut lr = initial_lr;
        let mut prev_lr = lr;

        // Test exponential decay schedule
        for step in 1..=100 {
            if step % step_size == 0 {
                lr *= decay_rate;
            }

            // Learning rate should always be positive
            prop_assert!(lr > 0.0);
            prop_assert!(lr.is_finite());

            // Learning rate should be non-increasing
            prop_assert!(lr <= prev_lr + 1e-10);

            prev_lr = lr;
        }

        // Final learning rate should be smaller than initial
        prop_assert!(lr <= initial_lr);
    }

    /// Test regularization properties
    #[test]
    fn test_regularization_properties(
        weight_decay in 0.0..1.0f64,
        dropout_rate in 0.0..0.8f64
    ) {
        let weights = Array1::from_vec(vec![1.0, -0.5, 2.0, -1.5, 0.8]);

        // Test L2 regularization (weight decay)
        let l2_penalty = weights.iter().map(|&w: &f64| w.powi(2)).sum::<f64>();
        let l2_loss = weight_decay * l2_penalty;

        prop_assert!(l2_loss >= 0.0);
        prop_assert!(l2_loss.is_finite());

        // L2 penalty should increase with larger weights
        let larger_weights = weights.mapv(|w| w * 2.0);
        let larger_l2_penalty = larger_weights.iter().map(|&w: &f64| w.powi(2)).sum::<f64>();
        prop_assert!(larger_l2_penalty >= l2_penalty);

        // Test dropout properties
        if dropout_rate > 0.0 && dropout_rate < 1.0 {
            // Dropout should preserve expected value during training
            let keep_prob = 1.0 - dropout_rate;
            let scale_factor = 1.0 / keep_prob;

            prop_assert!(scale_factor > 1.0);
            prop_assert!(scale_factor.is_finite());

            // Test that scaling compensates for dropped neurons
            let expected_output = weights.sum();
            let scaled_expected = expected_output * keep_prob * scale_factor;
            prop_assert!((scaled_expected - expected_output).abs() < 1e-10);
        }
    }

    /// Test self-supervised learning properties
    #[test]
    fn test_self_supervised_properties(
        batch_size in 2..10usize,
        embedding_dim in 8..32usize,
        temperature in 0.01..1.0f64
    ) {
        // Test cosine similarity properties with non-zero vectors
        let emb1 = Array1::from_shape_fn(embedding_dim, |i| (i as f64 / embedding_dim as f64) + 0.1);
        let emb2 = Array1::from_shape_fn(embedding_dim, |i| ((i + 1) as f64 / embedding_dim as f64) + 0.1);
        let emb3 = emb1.clone(); // Identical embedding

        let sim_different = cosine_similarity_test(&emb1, &emb2);
        let sim_identical = cosine_similarity_test(&emb1, &emb3);

        // Cosine similarity should be in [-1, 1]
        prop_assert!((-1.0..=1.0).contains(&sim_different));
        prop_assert!((-1.0..=1.0).contains(&sim_identical));

        // Identical vectors should have similarity close to 1 (allow for numerical precision)
        prop_assert!((sim_identical - 1.0).abs() < 1e-8);

        // Test contrastive loss properties
        let embeddings = Array2::from_shape_fn((batch_size, embedding_dim), |(i, j)| {
            ((i + j) as f64 / (batch_size + embedding_dim) as f64) + 0.01 // Ensure non-zero
        });

        let contrastive_loss = compute_simple_contrastive_loss(&embeddings, temperature);

        // Contrastive loss should be non-negative
        prop_assert!(contrastive_loss >= 0.0);
        prop_assert!(contrastive_loss.is_finite());

        // Test reconstruction loss properties for autoencoders
        let input = Array2::from_shape_fn((batch_size, embedding_dim), |(i, j)| {
            ((i + j) as f64 / 10.0) + 0.1 // Ensure non-zero
        });
        let reconstruction = input.mapv(|x| x + 0.05); // Slightly different reconstruction

        let mse_loss = reconstruction_loss_test(&input, &reconstruction);
        prop_assert!(mse_loss >= 0.0);
        prop_assert!(mse_loss.is_finite());

        // Perfect reconstruction should have zero loss
        let perfect_reconstruction = input.clone();
        let zero_loss = reconstruction_loss_test(&input, &perfect_reconstruction);
        prop_assert!(zero_loss < 1e-10);
    }

    /// Test model selection cross-validation properties
    #[test]
    fn test_cross_validation_properties(
        n_samples in 10..100usize,
        n_folds in 2..10usize,
        score_range in 0.0..1.0f64
    ) {
        prop_assume!(n_folds <= n_samples);

        // Test fold creation properties
        let fold_size = n_samples / n_folds;
        let remainder = n_samples % n_folds;

        // Each fold should have approximately equal size
        let expected_fold_sizes: Vec<usize> = (0..n_folds)
            .map(|i| if i < remainder { fold_size + 1 } else { fold_size })
            .collect();

        let total_expected: usize = expected_fold_sizes.iter().sum();
        prop_assert_eq!(total_expected, n_samples);

        // Test cross-validation score aggregation
        let fold_scores: Vec<f64> = (0..n_folds)
            .map(|i| score_range * (i as f64 / n_folds as f64))
            .collect();

        let mean_score = fold_scores.iter().sum::<f64>() / n_folds as f64;
        let variance = fold_scores.iter()
            .map(|&score| (score - mean_score).powi(2))
            .sum::<f64>() / n_folds as f64;
        let std_score = variance.sqrt();

        // Statistical properties should be valid
        prop_assert!(mean_score >= 0.0 && mean_score <= score_range);
        prop_assert!(std_score >= 0.0);
        prop_assert!(variance >= 0.0);
        prop_assert!(mean_score.is_finite());
        prop_assert!(std_score.is_finite());

        // Standard deviation should be less than or equal to the range
        prop_assert!(std_score <= score_range);
    }

    /// Test gradient checking numerical stability
    #[test]
    fn test_gradient_checking_properties(
        epsilon in 1e-8..1e-4f64,
        analytical_grad in -10.0..10.0f64,
        perturbation in -0.1..0.1f64
    ) {
        // Test finite difference approximation properties
        let numerical_grad = analytical_grad + perturbation;

        // Test relative error computation
        let abs_error = (analytical_grad - numerical_grad).abs();
        let rel_error = if numerical_grad.abs() > 0.0 {
            abs_error / numerical_grad.abs()
        } else {
            abs_error
        };

        // Error metrics should be non-negative and finite
        prop_assert!(abs_error >= 0.0);
        prop_assert!(rel_error >= 0.0);
        prop_assert!(abs_error.is_finite());
        prop_assert!(rel_error.is_finite());

        // Test centered vs forward difference properties
        let h = epsilon;
        let f_plus = analytical_grad * (1.0 + h);
        let f_minus = analytical_grad * (1.0 - h);
        let f_center = analytical_grad;

        let forward_diff = (f_plus - f_center) / h;
        let centered_diff = (f_plus - f_minus) / (2.0 * h);

        // Both finite difference methods should give finite results
        prop_assert!(forward_diff.is_finite());
        prop_assert!(centered_diff.is_finite());

        // Centered difference should generally be more accurate
        let forward_error = (forward_diff - analytical_grad).abs();
        let centered_error = (centered_diff - analytical_grad).abs();

        // Allow for numerical precision issues, but centered should not be worse by orders of magnitude
        prop_assert!(centered_error <= forward_error * 10.0 || centered_error < h);

        // Test gradient tolerance checking
        let relative_tolerance = 1e-5;
        let absolute_tolerance = 1e-8;

        let passes_check = rel_error < relative_tolerance && abs_error < absolute_tolerance;

        // If perturbation is small enough, gradient check should pass
        if perturbation.abs() < 1e-7 {
            prop_assert!(passes_check);
        }
    }

    /// Test model comparison ranking properties
    #[test]
    fn test_model_ranking_properties(
        num_models in 2..10usize,
        score_base in 0.0..0.8f64
    ) {
        // Create model scores with deterministic differences
        let model_scores: Vec<f64> = (0..num_models)
            .map(|i| score_base + (i as f64 / num_models as f64) * 0.2)
            .collect();

        // Test ranking properties
        let mut sorted_scores = model_scores.clone();
        sorted_scores.sort_by(|a, b| b.partial_cmp(a).expect("operation should succeed")); // Descending order

        // Rankings should be in descending order
        for i in 1..sorted_scores.len() {
            prop_assert!(sorted_scores[i-1] >= sorted_scores[i]);
        }

        // Best and worst model identification
        let best_score = model_scores.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
        let worst_score = model_scores.iter().cloned().fold(f64::INFINITY, f64::min);

        prop_assert!(best_score >= worst_score);
        prop_assert!(best_score == sorted_scores[0]);
        prop_assert!(worst_score == sorted_scores[sorted_scores.len() - 1]);

        // Test score normalization properties
        let score_range = best_score - worst_score;
        if score_range > 0.0 {
            let normalized_scores: Vec<f64> = model_scores.iter()
                .map(|&score| (score - worst_score) / score_range)
                .collect();

            for &norm_score in &normalized_scores {
                prop_assert!((0.0..=1.0).contains(&norm_score));
                prop_assert!(norm_score.is_finite());
            }

            // Best normalized score should be 1.0, worst should be 0.0
            let best_normalized = normalized_scores.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
            let worst_normalized = normalized_scores.iter().cloned().fold(f64::INFINITY, f64::min);

            prop_assert!((best_normalized - 1.0).abs() < 1e-10);
            prop_assert!(worst_normalized.abs() < 1e-10);
        }
    }
}

// Helper functions for neural network operations
fn relu_activation(x: f64) -> f64 {
    x.max(0.0)
}

fn sigmoid_activation(x: f64) -> f64 {
    1.0 / (1.0 + (-x).exp())
}

fn tanh_activation(x: f64) -> f64 {
    x.tanh()
}

fn xavier_init(n_input: usize, n_output: usize, seed: u64) -> Array2<f64> {
    // Simple deterministic Xavier initialization for testing
    let limit = (6.0 / (n_input + n_output) as f64).sqrt();

    Array2::from_shape_fn((n_output, n_input), |(i, j)| {
        let val = ((i + j + seed as usize) as f64 / 1000.0) % 1.0;
        (val - 0.5) * 2.0 * limit
    })
}

fn he_init(n_input: usize, n_output: usize, seed: u64) -> Array2<f64> {
    // Simple deterministic He initialization for testing (simplified to avoid Box-Muller issues)
    let std = (2.0 / n_input as f64).sqrt();

    Array2::from_shape_fn((n_output, n_input), |(i, j)| {
        let val = ((i + j + seed as usize) as f64 / 1000.0) % 1.0;
        // Simple uniform-to-normal approximation (not perfect but stable)
        let normal_approx = (val - 0.5) * 3.464; // approximately normal(-1.732, 1.732)
        normal_approx * std
    })
}

// Helper functions for self-supervised learning property tests
fn cosine_similarity_test(a: &Array1<f64>, b: &Array1<f64>) -> f64 {
    let dot_product = a.dot(b);
    let norm_a = a.mapv(|x| x * x).sum().sqrt();
    let norm_b = b.mapv(|x| x * x).sum().sqrt();

    // Add epsilon for numerical stability
    let epsilon = 1e-12;
    if norm_a < epsilon || norm_b < epsilon {
        return 0.0;
    }

    let similarity = dot_product / (norm_a * norm_b);

    // Clamp to [-1, 1] to handle numerical precision issues
    similarity.max(-1.0).min(1.0)
}

fn compute_simple_contrastive_loss(embeddings: &Array2<f64>, temperature: f64) -> f64 {
    let batch_size = embeddings.nrows();
    let mut total_loss = 0.0;

    for i in 0..batch_size {
        let anchor = embeddings.row(i);

        // Use next sample as positive
        let positive_idx = (i + 1) % batch_size;
        let positive = embeddings.row(positive_idx);
        let positive_sim = cosine_similarity_test(&anchor.to_owned(), &positive.to_owned());

        // Compute negative similarities
        let mut negative_sims = Vec::new();
        for j in 0..batch_size {
            if j != i && j != positive_idx {
                let negative = embeddings.row(j);
                let neg_sim = cosine_similarity_test(&anchor.to_owned(), &negative.to_owned());
                negative_sims.push(neg_sim);
            }
        }

        // Compute contrastive loss
        let pos_exp = (positive_sim / temperature).exp();
        let neg_exp_sum: f64 = negative_sims
            .iter()
            .map(|&sim| (sim / temperature).exp())
            .sum();

        if pos_exp + neg_exp_sum > 0.0 {
            let loss = -(pos_exp / (pos_exp + neg_exp_sum)).ln();
            total_loss += loss;
        }
    }

    total_loss / batch_size as f64
}

fn reconstruction_loss_test(input: &Array2<f64>, reconstruction: &Array2<f64>) -> f64 {
    let diff = input - reconstruction;
    diff.mapv(|x| x * x).mean().unwrap_or(0.0)
}

#[allow(non_snake_case)]
#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_activation_functions() {
        // Test ReLU
        assert_eq!(relu_activation(5.0), 5.0);
        assert_eq!(relu_activation(-3.0), 0.0);
        assert_eq!(relu_activation(0.0), 0.0);

        // Test Sigmoid
        let sig_0 = sigmoid_activation(0.0);
        assert!((sig_0 - 0.5).abs() < 1e-10);

        let sig_large = sigmoid_activation(10.0);
        assert!(sig_large > 0.99);

        let sig_small = sigmoid_activation(-10.0);
        assert!(sig_small < 0.01);

        // Test Tanh
        let tanh_0 = tanh_activation(0.0);
        assert!(tanh_0.abs() < 1e-10);

        let tanh_large = tanh_activation(10.0);
        assert!(tanh_large > 0.99);

        let tanh_small = tanh_activation(-10.0);
        assert!(tanh_small < -0.99);
    }

    #[test]
    fn test_weight_initialization() {
        let xavier = xavier_init(4, 3, 42);
        assert_eq!(xavier.shape(), &[3, 4]);

        let he = he_init(4, 3, 42);
        assert_eq!(he.shape(), &[3, 4]);

        // All weights should be finite
        for &w in xavier.iter() {
            assert!(w.is_finite());
        }
        for &w in he.iter() {
            assert!(w.is_finite());
        }
    }
}