scirs2_core/array_protocol/

training.rs

// Copyright (c) 2025, `SciRS2` Team
//
// Licensed under either of
//
// * Apache License, Version 2.0
//   (LICENSE-APACHE or http://www.apache.org/licenses/LICENSE-2.0)
// * MIT license
//   (LICENSE-MIT or http://opensource.org/licenses/MIT)
//
// at your option.
//

//! Training utilities for neural networks using the array protocol.
//!
//! This module provides utilities for training neural networks using the
//! array protocol, including datasets, dataloaders, loss functions, and
//! training loops.

use std::fmt;
use std::time::Instant;

use ndarray::{Array, Array0, Dimension};
use rand::seq::SliceRandom;
use rand::Rng;
use rand::SeedableRng;

use crate::array_protocol::grad::{GradientDict, Optimizer};
use crate::array_protocol::ml_ops::ActivationFunc;
use crate::array_protocol::neural::Sequential;
use crate::array_protocol::operations::{multiply, subtract};
use crate::array_protocol::{activation, ArrayProtocol, NdarrayWrapper};
use crate::error::{CoreError, CoreResult, ErrorContext};

/// Type alias for batch data
pub type BatchData = (Vec<Box<dyn ArrayProtocol>>, Vec<Box<dyn ArrayProtocol>>);

/// Dataset trait for providing data samples.
pub trait Dataset {
    /// Get the number of samples in the dataset.
    fn len(&self) -> usize;

    /// Check if the dataset is empty.
    fn is_empty(&self) -> bool {
        self.len() == 0
    }

    /// Get a sample from the dataset by index.
    fn get(&self, index: usize) -> Option<(Box<dyn ArrayProtocol>, Box<dyn ArrayProtocol>)>;

    /// Get the input shape of the dataset.
    fn inputshape(&self) -> Vec<usize>;

    /// Get the output shape of the dataset.
    fn outputshape(&self) -> Vec<usize>;
}

/// In-memory dataset with arrays.
pub struct InMemoryDataset {
    /// Input data samples.
    inputs: Vec<Box<dyn ArrayProtocol>>,

    /// Target output samples.
    targets: Vec<Box<dyn ArrayProtocol>>,

    /// Input shape.
    inputshape: Vec<usize>,

    /// Output shape.
    outputshape: Vec<usize>,
}

impl InMemoryDataset {
    /// Create a new in-memory dataset.
    pub fn new(
        inputs: Vec<Box<dyn ArrayProtocol>>,
        targets: Vec<Box<dyn ArrayProtocol>>,
        inputshape: Vec<usize>,
        outputshape: Vec<usize>,
    ) -> Self {
        assert_eq!(
            inputs.len(),
            targets.len(),
            "Inputs and targets must have the same length"
        );

        Self {
            inputs,
            targets,
            inputshape,
            outputshape,
        }
    }

    /// Create an in-memory dataset from arrays.
    pub fn from_arrays<T, D1, D2>(inputs: Array<T, D1>, targets: Array<T, D2>) -> Self
    where
        T: Clone + Send + Sync + 'static,
        D1: Dimension + Send + Sync,
        D2: Dimension + Send + Sync,
    {
        let inputshape = inputs.shape().to_vec();
        let outputshape = targets.shape().to_vec();

        // Handle batched _inputs
        let num_samples = inputshape[0];
        assert_eq!(
            num_samples, outputshape[0],
            "Inputs and targets must have the same number of samples"
        );

        let mut input_samples = Vec::with_capacity(num_samples);
        let mut target_samples = Vec::with_capacity(num_samples);

        // Create dynamic arrays with the appropriate shape to handle arbitrary dimensions
        let to_dyn_inputs = inputs.into_dyn();
        let to_dyn_targets = targets.into_dyn();

        for i in 0..num_samples {
            // Use index_axis instead of slice for better compatibility with different dimensions
            let input_view = to_dyn_inputs.index_axis(ndarray::Axis(0), i);
            let inputarray = input_view.to_owned();
            input_samples.push(Box::new(NdarrayWrapper::new(inputarray)) as Box<dyn ArrayProtocol>);

            let target_view = to_dyn_targets.index_axis(ndarray::Axis(0), i);
            let target_array = target_view.to_owned();
            target_samples
                .push(Box::new(NdarrayWrapper::new(target_array)) as Box<dyn ArrayProtocol>);
        }

        Self {
            inputs: input_samples,
            targets: target_samples,
            inputshape: inputshape[1..].to_vec(),
            outputshape: outputshape[1..].to_vec(),
        }
    }
}

impl Dataset for InMemoryDataset {
    fn len(&self) -> usize {
        self.inputs.len()
    }

    fn get(&self, index: usize) -> Option<(Box<dyn ArrayProtocol>, Box<dyn ArrayProtocol>)> {
        if index >= self.len() {
            return None;
        }

        Some((self.inputs[index].clone(), self.targets[index].clone()))
    }

    fn inputshape(&self) -> Vec<usize> {
        self.inputshape.clone()
    }

    fn outputshape(&self) -> Vec<usize> {
        self.outputshape.clone()
    }
}

/// Data loader for batching and shuffling datasets.
pub struct DataLoader {
    /// The dataset to load from.
    dataset: Box<dyn Dataset>,

    /// Batch size.
    batch_size: usize,

    /// Whether to shuffle the dataset.
    shuffle: bool,

    /// Random number generator seed.
    seed: Option<u64>,

    /// Indices of the dataset.
    indices: Vec<usize>,

    /// Current position in the dataset.
    position: usize,
}

impl DataLoader {
    /// Create a new data loader.
    pub fn new(
        dataset: Box<dyn Dataset>,
        batch_size: usize,
        shuffle: bool,
        seed: Option<u64>,
    ) -> Self {
        let indices = (0..dataset.len()).collect();

        Self {
            dataset,
            batch_size,
            shuffle,
            seed,
            indices,
            position: 0,
        }
    }

    /// Reset the data loader.
    pub fn reset(&mut self) {
        self.position = 0;

        if self.shuffle {
            let mut rng = match self.seed {
                Some(s) => rand::rngs::StdRng::seed_from_u64(s),
                None => {
                    let mut rng = rand::rng();
                    // Get a random seed from rng and create a new StdRng
                    let random_seed: u64 = rng.random();
                    rand::rngs::StdRng::seed_from_u64(random_seed)
                }
            };

            self.indices.shuffle(&mut rng);
        }
    }

    /// Get the next batch from the dataset.
    pub fn next_batch(&mut self) -> Option<BatchData> {
        if self.position >= self.dataset.len() {
            return None;
        }

        // Determine how many samples to take
        let remaining = self.dataset.len() - self.position;
        let batch_size = std::cmp::min(self.batch_size, remaining);

        // Get the batch
        let mut inputs = Vec::with_capacity(batch_size);
        let mut targets = Vec::with_capacity(batch_size);

        for i in 0..batch_size {
            let index = self.indices[self.position + i];
            if let Some((input, target)) = self.dataset.get(index) {
                inputs.push(input);
                targets.push(target);
            }
        }

        // Update position
        self.position += batch_size;

        Some((inputs, targets))
    }

    /// Get the number of batches in the dataset.
    pub fn numbatches(&self) -> usize {
        self.dataset.len().div_ceil(self.batch_size)
    }

    /// Get a reference to the dataset.
    pub fn dataset(&self) -> &dyn Dataset {
        self.dataset.as_ref()
    }
}

/// Iterator implementation for DataLoader.
impl Iterator for DataLoader {
    type Item = BatchData;

    fn next(&mut self) -> Option<Self::Item> {
        self.next_batch()
    }
}

/// Loss function trait.
pub trait Loss {
    /// Compute the loss between predictions and targets.
    fn forward(
        &self,
        predictions: &dyn ArrayProtocol,
        targets: &dyn ArrayProtocol,
    ) -> CoreResult<Box<dyn ArrayProtocol>>;

    /// Compute the gradient of the loss with respect to predictions.
    fn backward(
        &self,
        predictions: &dyn ArrayProtocol,
        targets: &dyn ArrayProtocol,
    ) -> CoreResult<Box<dyn ArrayProtocol>>;

    /// Get the name of the loss function.
    fn name(&self) -> &str;
}

/// Mean squared error loss.
pub struct MSELoss {
    /// Name of the loss function.
    name: String,

    /// Whether to reduce the loss.
    reduction: String,
}

impl MSELoss {
    /// Create a new MSE loss.
    pub fn new(reduction: Option<&str>) -> Self {
        Self {
            name: "MSELoss".to_string(),
            reduction: reduction.unwrap_or("mean").to_string(),
        }
    }
}

impl Loss for MSELoss {
    fn forward(
        &self,
        predictions: &dyn ArrayProtocol,
        targets: &dyn ArrayProtocol,
    ) -> CoreResult<Box<dyn ArrayProtocol>> {
        // Compute squared difference
        let diff = subtract(predictions, targets)?;
        let squared = multiply(diff.as_ref(), diff.as_ref())?;

        // Apply reduction
        match self.reduction.as_str() {
            "none" => Ok(squared),
            "mean" => {
                // Compute mean of all elements
                if let Some(array) = squared
                    .as_any()
                    .downcast_ref::<NdarrayWrapper<f64, ndarray::IxDyn>>()
                {
                    let mean = array.as_array().mean().unwrap();
                    let result = Array0::<f64>::from_elem((), mean);
                    Ok(Box::new(NdarrayWrapper::new(result)))
                } else {
                    Err(CoreError::NotImplementedError(ErrorContext::new(
                        "Mean reduction not implemented for this array type".to_string(),
                    )))
                }
            }
            "sum" => {
                // Compute sum of all elements
                if let Some(array) = squared
                    .as_any()
                    .downcast_ref::<NdarrayWrapper<f64, ndarray::IxDyn>>()
                {
                    let sum = array.as_array().sum();
                    let result = Array0::<f64>::from_elem((), sum);
                    Ok(Box::new(NdarrayWrapper::new(result)))
                } else {
                    Err(CoreError::NotImplementedError(ErrorContext::new(
                        "Sum reduction not implemented for this array type".to_string(),
                    )))
                }
            }
            _ => Err(CoreError::InvalidArgument(ErrorContext::new(format!(
                "Unknown reduction: {reduction}",
                reduction = self.reduction
            )))),
        }
    }

    fn backward(
        &self,
        predictions: &dyn ArrayProtocol,
        targets: &dyn ArrayProtocol,
    ) -> CoreResult<Box<dyn ArrayProtocol>> {
        // Gradient of MSE loss: 2 * (predictions - targets)
        let diff = subtract(predictions, targets)?;
        let factor = Box::new(NdarrayWrapper::new(ndarray::Array0::<f64>::from_elem(
            (),
            2.0,
        )));
        let grad = multiply(factor.as_ref(), diff.as_ref())?;

        // Apply reduction scaling if needed
        match self.reduction.as_str() {
            "none" => Ok(grad),
            "mean" => {
                // For mean reduction, scale by 1/N
                if let Some(array) = grad
                    .as_any()
                    .downcast_ref::<NdarrayWrapper<f64, ndarray::IxDyn>>()
                {
                    let n = array.as_array().len() as f64;
                    let scale_factor = Box::new(NdarrayWrapper::new(
                        ndarray::Array0::<f64>::from_elem((), 1.0 / n),
                    ));
                    Ok(multiply(scale_factor.as_ref(), grad.as_ref())?)
                } else {
                    Ok(grad)
                }
            }
            "sum" => Ok(grad),
            _ => Err(CoreError::InvalidArgument(ErrorContext::new(format!(
                "Unknown reduction: {reduction}",
                reduction = self.reduction
            )))),
        }
    }

    fn name(&self) -> &str {
        &self.name
    }
}

/// Cross-entropy loss.
pub struct CrossEntropyLoss {
    /// Name of the loss function.
    name: String,

    /// Whether to reduce the loss.
    reduction: String,
}

impl CrossEntropyLoss {
    /// Create a new cross-entropy loss.
    pub fn new(reduction: Option<&str>) -> Self {
        Self {
            name: "CrossEntropyLoss".to_string(),
            reduction: reduction.unwrap_or("mean").to_string(),
        }
    }
}

impl Loss for CrossEntropyLoss {
    fn forward(
        &self,
        predictions: &dyn ArrayProtocol,
        targets: &dyn ArrayProtocol,
    ) -> CoreResult<Box<dyn ArrayProtocol>> {
        // Apply softmax to predictions
        let softmax_preds = activation(predictions, ActivationFunc::Softmax)?;

        // Compute cross-entropy
        if let (Some(preds_array), Some(targets_array)) = (
            softmax_preds
                .as_any()
                .downcast_ref::<NdarrayWrapper<f64, ndarray::IxDyn>>(),
            targets
                .as_any()
                .downcast_ref::<NdarrayWrapper<f64, ndarray::IxDyn>>(),
        ) {
            let preds = preds_array.as_array();
            let targets = targets_array.as_array();

            // Compute -targets * log(preds)
            let log_preds = preds.mapv(|x| x.max(1e-10).ln());

            // Compute element-wise multiplication and then negate
            let mut losses = targets.clone();
            losses.zip_mut_with(&log_preds, |t, l| *t = -(*t * *l));

            // Apply reduction
            match self.reduction.as_str() {
                "none" => Ok(Box::new(NdarrayWrapper::new(losses))),
                "mean" => {
                    let mean = losses.mean().unwrap();
                    let result = Array0::<f64>::from_elem((), mean);
                    Ok(Box::new(NdarrayWrapper::new(result)))
                }
                "sum" => {
                    let sum = losses.sum();
                    let result = Array0::<f64>::from_elem((), sum);
                    Ok(Box::new(NdarrayWrapper::new(result)))
                }
                _ => Err(CoreError::InvalidArgument(ErrorContext::new(format!(
                    "Unknown reduction: {reduction}",
                    reduction = self.reduction
                )))),
            }
        } else {
            Err(CoreError::NotImplementedError(ErrorContext::new(
                "CrossEntropy not implemented for these array types".to_string(),
            )))
        }
    }

    fn backward(
        &self,
        predictions: &dyn ArrayProtocol,
        targets: &dyn ArrayProtocol,
    ) -> CoreResult<Box<dyn ArrayProtocol>> {
        // For cross-entropy with softmax: gradient is softmax(predictions) - targets
        let softmax_preds = activation(predictions, ActivationFunc::Softmax)?;
        let grad = subtract(softmax_preds.as_ref(), targets)?;

        // Apply reduction scaling if needed
        match self.reduction.as_str() {
            "none" => Ok(grad),
            "mean" => {
                // For mean reduction, scale by 1/N
                if let Some(array) = grad
                    .as_any()
                    .downcast_ref::<NdarrayWrapper<f64, ndarray::IxDyn>>()
                {
                    let n = array.as_array().len() as f64;
                    let scale_factor = Box::new(NdarrayWrapper::new(
                        ndarray::Array0::<f64>::from_elem((), 1.0 / n),
                    ));
                    Ok(multiply(scale_factor.as_ref(), grad.as_ref())?)
                } else {
                    Ok(grad)
                }
            }
            "sum" => Ok(grad),
            _ => Err(CoreError::InvalidArgument(ErrorContext::new(format!(
                "Unknown reduction: {reduction}",
                reduction = self.reduction
            )))),
        }
    }

    fn name(&self) -> &str {
        &self.name
    }
}

/// Metrics for evaluating model performance.
pub struct Metrics {
    /// Loss values.
    losses: Vec<f64>,

    /// Accuracy values (if applicable).
    accuracies: Option<Vec<f64>>,

    /// Name of the metrics object.
    name: String,
}

impl Metrics {
    /// Create a new metrics object.
    pub fn new(name: &str) -> Self {
        Self {
            losses: Vec::new(),
            accuracies: None,
            name: name.to_string(),
        }
    }

    /// Add a loss value.
    pub fn add_loss(&mut self, loss: f64) {
        self.losses.push(loss);
    }

    /// Add an accuracy value.
    pub fn add_accuracy(&mut self, accuracy: f64) {
        if self.accuracies.is_none() {
            self.accuracies = Some(Vec::new());
        }

        if let Some(accuracies) = &mut self.accuracies {
            accuracies.push(accuracy);
        }
    }

    /// Get the mean loss.
    pub fn mean_loss(&self) -> Option<f64> {
        if self.losses.is_empty() {
            return None;
        }

        let sum: f64 = self.losses.iter().sum();
        Some(sum / self.losses.len() as f64)
    }

    /// Get the mean accuracy.
    pub fn mean_accuracy(&self) -> Option<f64> {
        if let Some(accuracies) = &self.accuracies {
            if accuracies.is_empty() {
                return None;
            }

            let sum: f64 = accuracies.iter().sum();
            Some(sum / accuracies.len() as f64)
        } else {
            None
        }
    }

    /// Reset the metrics.
    pub fn reset(&mut self) {
        self.losses.clear();
        if let Some(accuracies) = &mut self.accuracies {
            accuracies.clear();
        }
    }

    /// Get the name of the metrics object.
    pub fn name(&self) -> &str {
        &self.name
    }
}

impl fmt::Display for Metrics {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(
            f,
            "{}: loss = {:.4}",
            self.name,
            self.mean_loss().unwrap_or(0.0)
        )?;

        if let Some(acc) = self.mean_accuracy() {
            write!(f, ", accuracy = {acc:.4}")?;
        }

        Ok(())
    }
}

/// Training progress callback trait.
pub trait TrainingCallback {
    /// Called at the start of each epoch.
    fn on_epoch_start(&mut self, epoch: usize, numepochs: usize);

    /// Called at the end of each epoch.
    fn on_epoch_end(&mut self, epoch: usize, numepochs: usize, metrics: &Metrics);

    /// Called at the start of each batch.
    fn on_batch_start(&mut self, batch: usize, numbatches: usize);

    /// Called at the end of each batch.
    fn on_batch_end(&mut self, batch: usize, numbatches: usize, loss: f64);

    /// Called at the start of training.
    fn on_train_start(&mut self, numepochs: usize);

    /// Called at the end of training.
    fn on_train_end(&mut self, metrics: &Metrics);
}

/// Progress bar callback for displaying training progress.
pub struct ProgressCallback {
    /// Whether to display a progress bar.
    verbose: bool,

    /// Start time of the current epoch.
    epoch_start: Option<Instant>,

    /// Start time of training.
    train_start: Option<Instant>,
}

impl ProgressCallback {
    /// Create a new progress callback.
    pub fn new(verbose: bool) -> Self {
        Self {
            verbose,
            epoch_start: None,
            train_start: None,
        }
    }
}

impl TrainingCallback for ProgressCallback {
    fn on_epoch_start(&mut self, epoch: usize, numepochs: usize) {
        if self.verbose {
            println!("Epoch {}/{}", epoch + 1, numepochs);
        }

        self.epoch_start = Some(Instant::now());
    }

    fn on_epoch_end(&mut self, _epoch: usize, numepochs: usize, metrics: &Metrics) {
        if self.verbose {
            if let Some(start) = self.epoch_start {
                let duration = start.elapsed();
                println!("{} - {}ms", metrics, duration.as_millis());
            } else {
                println!("{metrics}");
            }
        }
    }

    fn on_batch_start(&mut self, _batch: usize, _numbatches: usize) {
        // No-op for this callback
    }

    fn on_batch_end(&mut self, batch: usize, numbatches: usize, loss: f64) {
        if self.verbose && (batch + 1).is_multiple_of((numbatches / 10).max(1)) {
            print!("\rBatch {}/{} - loss: {:.4}", batch + 1, numbatches, loss);
            if batch + 1 == numbatches {
                println!();
            }
        }
    }

    fn on_train_start(&mut self, numepochs: usize) {
        if self.verbose {
            println!("Starting training for {numepochs} epochs");
        }

        self.train_start = Some(Instant::now());
    }

    fn on_train_end(&mut self, metrics: &Metrics) {
        if self.verbose {
            if let Some(start) = self.train_start {
                let duration = start.elapsed();
                println!("Training completed in {}s", duration.as_secs());
            } else {
                println!("Training completed");
            }

            if let Some(acc) = metrics.mean_accuracy() {
                println!("Final accuracy: {acc:.4}");
            }
        }
    }
}

/// Model trainer for neural networks.
pub struct Trainer {
    /// The model to train.
    model: Sequential,

    /// The optimizer to use.
    optimizer: Box<dyn Optimizer>,

    /// The loss function to use.
    lossfn: Box<dyn Loss>,

    /// The callbacks to use during training.
    callbacks: Vec<Box<dyn TrainingCallback>>,

    /// Training metrics.
    train_metrics: Metrics,

    /// Validation metrics.
    val_metrics: Option<Metrics>,
}

impl Trainer {
    /// Create a new trainer.
    pub fn new(model: Sequential, optimizer: Box<dyn Optimizer>, lossfn: Box<dyn Loss>) -> Self {
        Self {
            model,
            optimizer,
            lossfn,
            callbacks: Vec::new(),
            train_metrics: Metrics::new("train"),
            val_metrics: None,
        }
    }

    /// Add a callback to the trainer.
    pub fn add_callback(&mut self, callback: Box<dyn TrainingCallback>) {
        self.callbacks.push(callback);
    }

    /// Train the model.
    pub fn train(
        &mut self,
        train_loader: &mut DataLoader,
        numepochs: usize,
        mut val_loader: Option<&mut DataLoader>,
    ) -> CoreResult<()> {
        // Notify callbacks that training is starting
        for callback in &mut self.callbacks {
            callback.on_train_start(numepochs);
        }

        // Initialize validation metrics if needed
        if val_loader.is_some() && self.val_metrics.is_none() {
            self.val_metrics = Some(Metrics::new("val"));
        }

        // Train for the specified number of epochs
        for epoch in 0..numepochs {
            // Reset metrics
            self.train_metrics.reset();
            if let Some(metrics) = &mut self.val_metrics {
                metrics.reset();
            }

            // Notify callbacks that epoch is starting
            for callback in &mut self.callbacks {
                callback.on_epoch_start(epoch, numepochs);
            }

            // Train on the training set
            self.train_epoch(train_loader)?;

            // Validate on the validation set if provided
            if let Some(ref mut val_loader) = val_loader {
                self.validate(val_loader)?;
            }

            // Notify callbacks that epoch is ending
            for callback in &mut self.callbacks {
                callback.on_epoch_end(
                    epoch,
                    numepochs,
                    if let Some(val_metrics) = &self.val_metrics {
                        val_metrics
                    } else {
                        &self.train_metrics
                    },
                );
            }
        }

        // Notify callbacks that training is ending
        for callback in &mut self.callbacks {
            callback.on_train_end(if let Some(val_metrics) = &self.val_metrics {
                val_metrics
            } else {
                &self.train_metrics
            });
        }

        Ok(())
    }

    /// Train for one epoch.
    fn train_epoch(&mut self, dataloader: &mut DataLoader) -> CoreResult<()> {
        // Set model to training mode
        self.model.train();

        // Reset data loader
        dataloader.reset();

        let numbatches = dataloader.numbatches();

        // Train on batches
        for batch_idx in 0..numbatches {
            let (inputs, targets) = dataloader.next_batch().unwrap();
            // Notify callbacks that batch is starting
            for callback in &mut self.callbacks {
                callback.on_batch_start(batch_idx, numbatches);
            }

            // Forward pass
            let batch_loss = self.train_batch(&inputs, &targets)?;

            // Update metrics
            self.train_metrics.add_loss(batch_loss);

            // Notify callbacks that batch is ending
            for callback in &mut self.callbacks {
                callback.on_batch_end(batch_idx, numbatches, batch_loss);
            }
        }

        Ok(())
    }

    /// Train on a single batch.
    fn train_batch(
        &mut self,
        inputs: &[Box<dyn ArrayProtocol>],
        targets: &[Box<dyn ArrayProtocol>],
    ) -> CoreResult<f64> {
        // Zero gradients
        self.optimizer.zero_grad();

        // Forward pass
        let mut batch_loss = 0.0;

        for (input, target) in inputs.iter().zip(targets.iter()) {
            // Forward pass through model
            let output = self.model.forward(input.as_ref())?;

            // Compute loss
            let loss = self.lossfn.forward(output.as_ref(), target.as_ref())?;

            // Get loss value
            if let Some(loss_array) = loss
                .as_any()
                .downcast_ref::<NdarrayWrapper<f64, ndarray::IxDyn>>()
            {
                let loss_value = loss_array.as_array().sum();
                batch_loss += loss_value;
            }

            // Backward pass - compute gradients
            // For now, implement a simple gradient approximation using finite differences
            // In a full implementation, this would be automatic differentiation

            let learningrate = 0.001; // Default learning rate

            // Simple gradient estimation for demonstration
            // This computes numerical gradients for the model parameters

            // Get current output for gradient computation
            let current_output = self.model.forward(input.as_ref())?;
            let current_loss = self
                .lossfn
                .forward(current_output.as_ref(), target.as_ref())?;
            let _current_loss_value = if let Some(loss_array) = current_loss
                .as_any()
                .downcast_ref::<NdarrayWrapper<f64, ndarray::IxDyn>>()
            {
                loss_array.as_array().sum()
            } else {
                0.0
            };

            // Compute gradients via backpropagation
            let gradients = self.compute_gradients(
                input.as_ref(),
                target.as_ref(),
                current_output.as_ref(),
                current_loss.as_ref(),
            )?;

            // Apply gradients to model parameters
            self.apply_gradients(&gradients, learningrate)?;

            // Store gradients in optimizer for momentum-based optimizers
            self.optimizer.accumulate_gradients(&gradients)?;
        }

        // Compute average loss
        let batch_loss = batch_loss / inputs.len() as f64;

        // Update weights
        self.optimizer.step()?;

        Ok(batch_loss)
    }

    /// Compute gradients via backpropagation
    fn compute_gradients(
        &self,
        input: &dyn ArrayProtocol,
        target: &dyn ArrayProtocol,
        output: &dyn ArrayProtocol,
        _loss: &dyn ArrayProtocol,
    ) -> CoreResult<GradientDict> {
        // Start backpropagation from loss
        let mut gradients = GradientDict::new();

        // Compute gradient of loss with respect to output
        let loss_grad = self.lossfn.backward(output, target)?;

        // Backpropagate through the model
        let model_gradients = self.model.backward(input, loss_grad.as_ref())?;

        // Merge gradients
        gradients.merge(model_gradients);

        Ok(gradients)
    }

    /// Apply computed gradients to model parameters
    fn apply_gradients(&mut self, gradients: &GradientDict, learningrate: f64) -> CoreResult<()> {
        // Apply gradients to each parameter in the model
        for (param_name, gradient) in gradients.iter() {
            self.model
                .update_parameter(param_name, gradient.as_ref(), learningrate)?;
        }

        Ok(())
    }

    /// Validate the model.
    fn validate(&mut self, dataloader: &mut DataLoader) -> CoreResult<()> {
        // Set model to evaluation mode
        self.model.eval();

        // Reset validation metrics
        if let Some(metrics) = &mut self.val_metrics {
            metrics.reset();
        } else {
            return Ok(());
        }

        // Reset data loader
        dataloader.reset();

        let numbatches = dataloader.numbatches();

        // Validate on batches
        for _ in 0..numbatches {
            let (inputs, targets) = dataloader.next_batch().unwrap();
            // Forward pass without gradient tracking
            let mut batch_loss = 0.0;
            let mut batch_correct = 0;
            let mut batch_total = 0;

            for (input, target) in inputs.iter().zip(targets.iter()) {
                // Forward pass through model
                let output = self.model.forward(input.as_ref())?;

                // Compute loss
                let loss = self.lossfn.forward(output.as_ref(), target.as_ref())?;

                // Get loss value
                if let Some(loss_array) = loss
                    .as_any()
                    .downcast_ref::<NdarrayWrapper<f64, ndarray::IxDyn>>()
                {
                    let loss_value = loss_array.as_array().sum();
                    batch_loss += loss_value;
                }

                // Compute accuracy for classification problems
                if let (Some(output_array), Some(target_array)) = (
                    output
                        .as_any()
                        .downcast_ref::<NdarrayWrapper<f64, ndarray::Ix2>>(),
                    target
                        .as_any()
                        .downcast_ref::<NdarrayWrapper<f64, ndarray::Ix2>>(),
                ) {
                    // Get predictions (argmax)
                    let output_vec = output_array.as_array();
                    let target_vec = target_array.as_array();

                    // For simplicity, assume 2D arrays [batch_size, num_classes]
                    if output_vec.ndim() == 2 && target_vec.ndim() == 2 {
                        for (out_row, target_row) in
                            output_vec.outer_iter().zip(target_vec.outer_iter())
                        {
                            // Find the index of the maximum value in the output row
                            let mut max_idx = 0;
                            let mut max_val = out_row[0];

                            for (i, &val) in out_row.iter().enumerate().skip(1) {
                                if val > max_val {
                                    max_idx = i;
                                    max_val = val;
                                }
                            }

                            // Find the index of 1 in the target row (one-hot encoding)
                            if let Some(target_idx) = target_row.iter().position(|&x| x == 1.0) {
                                if max_idx == target_idx {
                                    batch_correct += 1;
                                }
                            }

                            batch_total += 1;
                        }
                    }
                }
            }

            // Compute average loss and accuracy
            let batch_loss = batch_loss / inputs.len() as f64;
            let batch_accuracy = if batch_total > 0 {
                batch_correct as f64 / batch_total as f64
            } else {
                0.0
            };

            // Update validation metrics
            if let Some(metrics) = &mut self.val_metrics {
                metrics.add_loss(batch_loss);
                metrics.add_accuracy(batch_accuracy);
            }
        }

        Ok(())
    }

    /// Get training metrics.
    pub const fn train_metrics(&self) -> &Metrics {
        &self.train_metrics
    }

    /// Get validation metrics.
    pub fn val_metrics(&self) -> Option<&Metrics> {
        self.val_metrics.as_ref()
    }
}

// Helper functions

#[cfg(test)]
mod tests {
    use super::*;
    use crate::array_protocol::{self, NdarrayWrapper};
    use ndarray::Array2;

    #[test]
    fn test_in_memory_dataset() {
        // Create input and target arrays
        let inputs = Array2::<f64>::ones((10, 5));
        let targets = Array2::<f64>::zeros((10, 2));

        // Create dataset
        let dataset = InMemoryDataset::from_arrays(inputs, targets);

        // Check properties
        assert_eq!(dataset.len(), 10);
        assert_eq!(dataset.inputshape(), vec![5]);
        assert_eq!(dataset.outputshape(), vec![2]);

        // Get a sample
        let (input, target) = dataset.get(0).unwrap();
        assert!(input
            .as_any()
            .downcast_ref::<NdarrayWrapper<f64, ndarray::IxDyn>>()
            .is_some());
        assert!(target
            .as_any()
            .downcast_ref::<NdarrayWrapper<f64, ndarray::IxDyn>>()
            .is_some());
    }

    #[test]
    fn test_dataloader() {
        // Create input and target arrays
        let inputs = Array2::<f64>::ones((10, 5));
        let targets = Array2::<f64>::zeros((10, 2));

        // Create dataset and data loader
        let dataset = Box::new(InMemoryDataset::from_arrays(inputs, targets));
        let mut loader = DataLoader::new(dataset, 4, true, Some(42));

        // Check properties
        assert_eq!(loader.numbatches(), 3);

        // Get batches
        let (batch1_inputs, batch1_targets) = loader.next_batch().unwrap();
        assert_eq!(batch1_inputs.len(), 4);
        assert_eq!(batch1_targets.len(), 4);

        let (batch2_inputs, batch2_targets) = loader.next_batch().unwrap();
        assert_eq!(batch2_inputs.len(), 4);
        assert_eq!(batch2_targets.len(), 4);

        let (batch3_inputs, batch3_targets) = loader.next_batch().unwrap();
        assert_eq!(batch3_inputs.len(), 2);
        assert_eq!(batch3_targets.len(), 2);

        // Reset and get another batch
        loader.reset();
        let (batch1_inputs, batch1_targets) = loader.next_batch().unwrap();
        assert_eq!(batch1_inputs.len(), 4);
        assert_eq!(batch1_targets.len(), 4);
    }

    #[test]
    fn test_mse_loss() {
        // Initialize the array protocol system
        array_protocol::init();

        // Create prediction and target arrays
        let predictions = Array2::<f64>::ones((2, 3));
        let targets = Array2::<f64>::zeros((2, 3));

        let predictions_wrapped = NdarrayWrapper::new(predictions);
        let targets_wrapped = NdarrayWrapper::new(targets);

        // Create loss function
        let mse = MSELoss::new(Some("mean"));

        // Compute loss with proper error handling
        match mse.forward(&predictions_wrapped, &targets_wrapped) {
            Ok(loss) => {
                if let Some(loss_array) = loss
                    .as_any()
                    .downcast_ref::<NdarrayWrapper<f64, ndarray::Ix0>>()
                {
                    // Expected: mean((1 - 0)^2) = 1.0
                    assert_eq!(loss_array.as_array()[()], 1.0);
                } else {
                    println!("Loss not of expected type NdarrayWrapper<f64, Ix0>");
                }
            }
            Err(e) => {
                println!("MSE Loss forward not fully implemented: {e}");
            }
        }
    }

    #[test]
    fn test_metrics() {
        // Create metrics
        let mut metrics = Metrics::new("test");

        // Add loss values
        metrics.add_loss(1.0);
        metrics.add_loss(2.0);
        metrics.add_loss(3.0);

        // Add accuracy values
        metrics.add_accuracy(0.5);
        metrics.add_accuracy(0.6);
        metrics.add_accuracy(0.7);

        // Check mean values
        assert_eq!(metrics.mean_loss().unwrap(), 2.0);
        assert_eq!(metrics.mean_accuracy().unwrap(), 0.6);

        // Reset metrics
        metrics.reset();
        assert!(metrics.mean_loss().is_none());
        assert!(metrics.mean_accuracy().is_none());
    }
}