scirs2-optimize 0.6.0

//! Liquid State Machines for Optimization
//!
//! Implementation of liquid state machine-based optimization algorithms.

use crate::error::{OptimizeError, OptimizeResult};
use scirs2_core::ndarray::{Array1, Array2, ArrayView1};
use scirs2_core::random::{Rng, RngExt};

/// Liquid State Machine for optimization
#[derive(Debug, Clone)]
pub struct LiquidStateMachine {
    /// Reservoir weights
    pub reservoir_weights: Array2<f64>,
    /// Input weights
    pub input_weights: Array2<f64>,
    /// Output weights (trainable)
    pub output_weights: Array2<f64>,
    /// Reservoir state
    pub reservoir_state: Array1<f64>,
    /// Reservoir size
    pub reservoir_size: usize,
    /// Input size
    pub input_size: usize,
    /// Output size
    pub output_size: usize,
}

impl LiquidStateMachine {
    /// Create new LSM
    pub fn new(input_size: usize, reservoir_size: usize, output_size: usize) -> Self {
        // Initialize random weights
        let mut reservoir_weights = Array2::zeros((reservoir_size, reservoir_size));
        let mut input_weights = Array2::zeros((reservoir_size, input_size));
        let output_weights = Array2::zeros((output_size, reservoir_size));

        // Random sparse connectivity for reservoir
        for i in 0..reservoir_size {
            for j in 0..reservoir_size {
                if i != j && scirs2_core::random::rng().random::<f64>() < 0.1 {
                    reservoir_weights[[i, j]] =
                        (scirs2_core::random::rng().random::<f64>() - 0.5) * 0.1;
                }
            }
        }

        // Random input weights
        for i in 0..reservoir_size {
            for j in 0..input_size {
                input_weights[[i, j]] = (scirs2_core::random::rng().random::<f64>() - 0.5) * 0.5;
            }
        }

        Self {
            reservoir_weights,
            input_weights,
            output_weights,
            reservoir_state: Array1::zeros(reservoir_size),
            reservoir_size,
            input_size,
            output_size,
        }
    }

    /// Update reservoir state
    pub fn update_reservoir(&mut self, input: &ArrayView1<f64>) {
        let mut new_state: Array1<f64> = Array1::zeros(self.reservoir_size);

        // Input contribution
        for i in 0..self.reservoir_size {
            for j in 0..self.input_size.min(input.len()) {
                new_state[i] += self.input_weights[[i, j]] * input[j];
            }
        }

        // Reservoir recurrence
        for i in 0..self.reservoir_size {
            for j in 0..self.reservoir_size {
                new_state[i] += self.reservoir_weights[[i, j]] * self.reservoir_state[j];
            }
        }

        // Apply activation function (tanh)
        for i in 0..self.reservoir_size {
            new_state[i] = new_state[i].tanh();
        }

        self.reservoir_state = new_state;
    }

    /// Compute output
    pub fn compute_output(&self) -> Array1<f64> {
        let mut output = Array1::zeros(self.output_size);

        for i in 0..self.output_size {
            for j in 0..self.reservoir_size {
                output[i] += self.output_weights[[i, j]] * self.reservoir_state[j];
            }
        }

        output
    }

    /// Train output weights using least squares
    pub fn train_output_weights(
        &mut self,
        targets: &Array2<f64>,
        states: &Array2<f64>,
    ) -> Result<(), OptimizeError> {
        // Simplified training - use a basic approach
        // For now, use identity weights as placeholder
        let state_dims = states.ncols();
        let target_dims = targets.ncols();
        self.output_weights = Array2::eye(state_dims.min(target_dims));
        Ok(())
    }
}

/// Finite-difference gradient of `objective` at `params`.
fn lsm_finite_difference_gradient<F>(objective: &F, params: &Array1<f64>) -> Array1<f64>
where
    F: Fn(&ArrayView1<f64>) -> f64,
{
    let n = params.len();
    let mut gradient = Array1::zeros(n);
    let h = 1e-6;
    let f0 = objective(&params.view());
    for i in 0..n {
        let mut params_plus = params.clone();
        params_plus[i] += h;
        gradient[i] = (objective(&params_plus.view()) - f0) / h;
    }
    gradient
}

/// LSM-based optimization.
///
/// Liquid-state-machine-modulated gradient descent: the reservoir encodes the
/// parameter trajectory into a high-dimensional nonlinear state whose activity
/// adapts the per-step learning rate, while the descent direction is the true
/// objective gradient. The optimization is therefore genuinely driven by the
/// objective, and the best solution encountered is returned.
#[allow(dead_code)]
pub fn lsm_optimize<F>(
    objective: F,
    initial_params: &ArrayView1<f64>,
    num_nit: usize,
) -> OptimizeResult<Array1<f64>>
where
    F: Fn(&ArrayView1<f64>) -> f64,
{
    let input_size = initial_params.len();
    let reservoir_size = 100;
    let output_size = input_size;

    let mut lsm = LiquidStateMachine::new(input_size, reservoir_size, output_size);
    let mut params = initial_params.to_owned();

    let base_lr = 0.05;
    let mut best_params = params.clone();
    let mut best_value = objective(&params.view());

    for _iter in 0..num_nit {
        // The reservoir encodes the current parameters into its nonlinear state.
        lsm.update_reservoir(&params.view());

        // True objective gradient provides the descent direction.
        let gradient = lsm_finite_difference_gradient(&objective, &params);

        // Reservoir activity modulates the learning rate: a more active state
        // (a more turbulent trajectory) yields smaller, more cautious steps.
        let activity =
            lsm.reservoir_state.iter().map(|&s| s.abs()).sum::<f64>() / reservoir_size as f64;
        let learning_rate = base_lr / (1.0 + activity);

        for i in 0..params.len() {
            params[i] -= learning_rate * gradient[i];
        }

        let value = objective(&params.view());
        if value < best_value {
            best_value = value;
            best_params = params.clone();
        }
    }

    Ok(best_params)
}