optirs-core 0.3.2

// Memory-Augmented Neural Turing Machine (NTM) optimizer.
//
// Implements an NTM-style optimizer inspired by Graves et al. (2014),
// "Neural Turing Machines" (arXiv:1410.5401). The optimizer maintains an
// external memory matrix `M ∈ R^{N × W}` (where `N = memory_slots` and
// `W = memory_width`) which it reads from and writes to on every step
// using content-based addressing (cosine similarity attention) and a
// classic NTM-style erase+add write head.
//
// Unlike a "learned optimizer" (e.g. an LSTM that emits parameter updates),
// this NTM optimizer uses a fixed-function attention mechanism over a
// memory matrix that is updated by the optimizer itself. The memory acts
// as a slow-changing summary of recent gradient/update signals, and the
// read vector is mixed with the raw gradient to produce the final update.
//
// Algorithm (single read/write head, one step):
//
//   1. Build query key `k ∈ R^W` from the gradient via mean-pool + Z-norm.
//   2. Content addressing: for each row `m_i`, compute the cosine similarity
//      `c_i = (k · m_i) / (||k|| · ||m_i|| + ε)`. Apply a sharpened softmax
//      with sharpness `β` to obtain `w_c ∈ Δ^{N-1}`.
//   3. Final read weights `w_r` depend on the addressing mode:
//        * `Content`  : `w_r = w_c`
//        * `Location` : `w_r = shift(w_{r,prev}, +1)` (circular shift by +1)
//        * `Hybrid`   : `w_r ∝ ( 0.5 w_c + 0.5 shift(w_{r,prev}, +1) )^β`
//          (sharpen by power `β`, then renormalize)
//   4. Read vector: `r = w_r^T M`, shape `[W]`.
//   5. Write weights `w_w = w_r` (one-head model).
//   6. Erase + add (Graves §3.2):
//        `e = erase_gate * k`
//        `M[i] ← M[i] ∘ (1 − w_w[i] · e) + w_w[i] · k`
//   7. Final update is a weighted combination of the gradient and the
//      tiled read vector (the latter broadcast over the parameter shape
//      via cyclic tiling). Apply with the configured learning rate.
//
// All numeric primitives come from `scirs2_core::numeric::Float`, all array
// types from `scirs2_core::ndarray`, and the RNG from `scirs2_core::random`,
// per the project's "no direct ndarray/rand" policy.

use scirs2_core::ndarray::{Array, Array1, Array2, Dimension, IxDyn, ScalarOperand};
use scirs2_core::numeric::Float;
use scirs2_core::random::Random;
use std::fmt::Debug;

use crate::error::{OptimError, Result};
use crate::optimizers::Optimizer;

/// Numerical safety floor used wherever a quantity could otherwise produce a
/// division by zero (cosine similarity, softmax denominator, etc.).
const EPSILON: f64 = 1e-12;

/// Addressing mode used by the read head.
///
/// The write head always mirrors the read head in this single-head
/// implementation, so the same enum governs both attentions.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum AddressingMode {
    /// Content-based addressing only (pure cosine-similarity softmax). The
    /// previous attention vector is ignored.
    Content,
    /// Location-based addressing only. The read weights are a circular
    /// shift of the previous read weights by `+1`, regardless of memory
    /// contents.
    Location,
    /// Hybrid content + location addressing. Mixes the content-based
    /// distribution with a shifted copy of the previous weights, then
    /// sharpens the mixture by raising every entry to the
    /// `read_sharpness` power and renormalising.
    Hybrid,
}

/// Configuration for an [`NtmOptimizer`].
///
/// The defaults mirror the values recommended in the project specification:
/// 32 slots × 16-wide memory, learning rate `0.01`, hybrid addressing,
/// `read_sharpness = 1`, `erase_gate = 0.5`, and `0.7 / 0.3` weighting
/// between the raw gradient and the memory read vector.
#[derive(Debug, Clone)]
pub struct NtmConfig<A: Float + ScalarOperand + Debug> {
    /// Number of memory slots `N` (rows of the memory matrix).
    pub memory_slots: usize,
    /// Width of each memory slot `W` (columns of the memory matrix).
    pub memory_width: usize,
    /// Base learning rate applied to the final update.
    pub learning_rate: A,
    /// Read-attention sharpness factor `β` (the NTM paper's β).
    pub read_sharpness: A,
    /// Erase gate (a scalar in `[0, 1]`) controlling how aggressively the
    /// write head zeroes a memory cell before adding the new content.
    pub erase_gate: A,
    /// Addressing mode for the read head.
    pub addressing_mode: AddressingMode,
    /// Weight assigned to the (tiled) read vector when forming the update.
    pub memory_weight: A,
    /// Weight assigned to the raw gradient when forming the update.
    pub gradient_weight: A,
    /// RNG seed used for any stochastic initialisation. Memory itself is
    /// initialised to zero; the seed is kept so subclasses/tests can rely
    /// on deterministic behaviour.
    pub seed: u64,
}

impl<A: Float + ScalarOperand + Debug> Default for NtmConfig<A> {
    fn default() -> Self {
        Self {
            memory_slots: 32,
            memory_width: 16,
            learning_rate: A::from(0.01).unwrap_or_else(A::zero),
            read_sharpness: A::from(1.0).unwrap_or_else(A::one),
            erase_gate: A::from(0.5).unwrap_or_else(A::zero),
            addressing_mode: AddressingMode::Hybrid,
            memory_weight: A::from(0.3).unwrap_or_else(A::zero),
            gradient_weight: A::from(0.7).unwrap_or_else(A::one),
            seed: 42,
        }
    }
}

/// Memory-augmented NTM-style optimizer.
///
/// `NtmOptimizer` keeps an `N × W` memory matrix between calls to
/// [`Optimizer::step`]. On each step it constructs a query key from the
/// gradient, attends to the memory via content-based (cosine similarity)
/// addressing, optionally fuses the result with a shifted copy of the
/// previous attention, reads a vector from memory, writes the new key
/// back via the NTM erase/add rule, and finally combines the read vector
/// with the raw gradient to produce the parameter update.
///
/// See the module-level documentation for the full algorithm.
///
/// # Examples
///
/// ```
/// use scirs2_core::ndarray::Array1;
/// use optirs_core::optimizers::{NtmOptimizer, Optimizer};
///
/// let mut opt: NtmOptimizer<f64> = NtmOptimizer::new(8, 4, 0.01);
/// let params = Array1::from_vec(vec![1.0, -1.0, 0.5, 0.0]);
/// let grads = Array1::from_vec(vec![0.2, -0.2, 0.1, 0.0]);
/// let next = opt.step(&params, &grads).expect("ntm step");
/// assert_eq!(next.len(), 4);
/// ```
pub struct NtmOptimizer<A: Float + ScalarOperand + Debug> {
    /// Configuration (memory shape, hyperparameters, etc.).
    config: NtmConfig<A>,
    /// External memory matrix `M ∈ R^{N × W}`. Initialised to zero.
    memory: Array2<A>,
    /// Read attention weights from the previous step, shape `[N]`.
    prev_read_weights: Array1<A>,
    /// Write attention weights from the previous step, shape `[N]`.
    prev_write_weights: Array1<A>,
    /// Stored seed (kept for reproducibility / [`Self::reset`]).
    rng_seed: u64,
    /// Number of [`Optimizer::step`] calls completed so far.
    step_count: usize,
}

impl<A: Float + ScalarOperand + Debug> Debug for NtmOptimizer<A> {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.debug_struct("NtmOptimizer")
            .field("memory_slots", &self.config.memory_slots)
            .field("memory_width", &self.config.memory_width)
            .field("learning_rate", &self.config.learning_rate)
            .field("addressing_mode", &self.config.addressing_mode)
            .field("rng_seed", &self.rng_seed)
            .field("step_count", &self.step_count)
            .finish()
    }
}

impl<A: Float + ScalarOperand + Debug> NtmOptimizer<A> {
    /// Constructs an `NtmOptimizer` with default hyperparameters and the
    /// supplied memory shape and learning rate.
    ///
    /// The memory matrix is initialised to all zeros, as are the previous
    /// attention vectors. The RNG seed defaults to `42` (mirroring the
    /// default `NtmConfig`).
    pub fn new(memory_slots: usize, memory_width: usize, learning_rate: A) -> Self {
        let config = NtmConfig::<A> {
            memory_slots,
            memory_width,
            learning_rate,
            ..NtmConfig::<A>::default()
        };
        Self::with_config(config)
    }

    /// Constructs an `NtmOptimizer` from a fully populated [`NtmConfig`].
    ///
    /// The memory matrix is initialised to zero. A seeded RNG is allocated
    /// (and immediately dropped) so a downstream caller that adds a
    /// stochastic initialisation strategy in the future does not need to
    /// change the public constructor signature.
    pub fn with_config(config: NtmConfig<A>) -> Self {
        let slots = config.memory_slots;
        let width = config.memory_width;
        let seed = config.seed;
        // Touch the RNG so the deterministic seed is observed even if we
        // currently use it only for `reset`/reproducibility tests.
        let _rng: Random<scirs2_core::random::rngs::StdRng> = Random::seed(seed);

        let memory = Array2::<A>::zeros((slots, width));
        let prev_read = Array1::<A>::zeros(slots);
        let prev_write = Array1::<A>::zeros(slots);

        Self {
            config,
            memory,
            prev_read_weights: prev_read,
            prev_write_weights: prev_write,
            rng_seed: seed,
            step_count: 0,
        }
    }

    /// Sets the read-attention sharpness factor `β`.
    pub fn with_read_sharpness(mut self, beta: A) -> Self {
        self.config.read_sharpness = beta;
        self
    }

    /// Sets the erase-gate scalar (intended to live in `[0, 1]`).
    pub fn with_erase_gate(mut self, gate: A) -> Self {
        self.config.erase_gate = gate;
        self
    }

    /// Selects the addressing mode for the read head.
    pub fn with_addressing(mut self, mode: AddressingMode) -> Self {
        self.config.addressing_mode = mode;
        self
    }

    /// Sets the weight applied to the (tiled) read vector in the final update.
    pub fn with_memory_weight(mut self, weight: A) -> Self {
        self.config.memory_weight = weight;
        self
    }

    /// Sets the weight applied to the raw gradient in the final update.
    pub fn with_gradient_weight(mut self, weight: A) -> Self {
        self.config.gradient_weight = weight;
        self
    }

    /// Overrides the RNG seed (and rebuilds the internal RNG).
    pub fn with_seed(mut self, seed: u64) -> Self {
        self.config.seed = seed;
        self.rng_seed = seed;
        let _rng: Random<scirs2_core::random::rngs::StdRng> = Random::seed(seed);
        self
    }

    /// Returns the current configuration.
    pub fn config(&self) -> &NtmConfig<A> {
        &self.config
    }

    /// Read-only access to the memory matrix.
    pub fn memory(&self) -> &Array2<A> {
        &self.memory
    }

    /// Mutable access to the memory matrix. Intended primarily for tests
    /// that need to seed memory contents directly.
    pub fn memory_mut(&mut self) -> &mut Array2<A> {
        &mut self.memory
    }

    /// Returns the read attention vector emitted by the most recent step.
    pub fn last_read_weights(&self) -> &Array1<A> {
        &self.prev_read_weights
    }

    /// Returns the write attention vector emitted by the most recent step.
    pub fn last_write_weights(&self) -> &Array1<A> {
        &self.prev_write_weights
    }

    /// Returns the number of [`Optimizer::step`] calls completed so far.
    pub fn step_count(&self) -> usize {
        self.step_count
    }

    /// Resets the optimizer to its post-construction state: memory and
    /// previous attentions become zero and the step counter is cleared.
    pub fn reset(&mut self) {
        self.memory.fill(A::zero());
        self.prev_read_weights.fill(A::zero());
        self.prev_write_weights.fill(A::zero());
        self.step_count = 0;
        // Re-seed the throwaway RNG so any future stochastic init is
        // reproducible from the same seed across resets.
        let _rng: Random<scirs2_core::random::rngs::StdRng> = Random::seed(self.rng_seed);
    }

    /// Validate that the configuration is structurally usable. Called once
    /// at the start of every `step` so we can surface bad configurations
    /// even if they survived construction (e.g. were tweaked with builder
    /// methods after the fact).
    fn validate_config(&self) -> Result<()> {
        if self.config.memory_slots == 0 {
            return Err(OptimError::InvalidConfig(
                "NtmOptimizer: memory_slots must be > 0".to_string(),
            ));
        }
        if self.config.memory_width == 0 {
            return Err(OptimError::InvalidConfig(
                "NtmOptimizer: memory_width must be > 0".to_string(),
            ));
        }
        Ok(())
    }

    /// Build the read/write query key `k ∈ R^W` from the gradient using a
    /// mean-pool over uniformly-sized chunks of the gradient, followed by
    /// Z-normalisation. If the gradient is empty the returned key is the
    /// zero vector.
    fn build_key<D: Dimension>(&self, gradients: &Array<A, D>) -> Array1<A> {
        let w = self.config.memory_width;
        let mut key = Array1::<A>::zeros(w);
        if w == 0 {
            return key;
        }
        let flat: Vec<A> = gradients.iter().copied().collect();
        let n = flat.len();
        if n == 0 {
            return key;
        }

        // Mean-pool the gradient into `w` buckets. Bucket `j` covers indices
        // `[floor(j*n/w), floor((j+1)*n/w))`. When `n < w` the gradient is
        // simply replicated (each entry placed in the appropriate bucket).
        if n >= w {
            for j in 0..w {
                let lo = (j * n) / w;
                let hi = ((j + 1) * n) / w;
                let hi_safe = hi.max(lo + 1).min(n);
                let mut acc = A::zero();
                let mut count: usize = 0;
                for value in flat.iter().take(hi_safe).skip(lo) {
                    acc = acc + *value;
                    count += 1;
                }
                if count > 0 {
                    let denom = A::from(count).unwrap_or_else(A::one);
                    key[j] = acc / denom;
                }
            }
        } else {
            // Spread the n gradient entries across the w buckets.
            for (i, value) in flat.iter().enumerate() {
                let j = (i * w) / n;
                key[j] = key[j] + *value;
            }
            // Average within each bucket if it received multiple gradient
            // values (only possible when n > 1 but n < w is by construction
            // false here; left for defensive symmetry).
            let mut counts = vec![0_usize; w];
            for i in 0..n {
                let j = (i * w) / n;
                counts[j] += 1;
            }
            for (j, c) in counts.iter().enumerate() {
                if *c > 1 {
                    let denom = A::from(*c).unwrap_or_else(A::one);
                    key[j] = key[j] / denom;
                }
            }
        }

        // Z-normalise: subtract mean, divide by stddev + ε.
        let w_a = A::from(w).unwrap_or_else(A::one);
        let mean = key.iter().copied().fold(A::zero(), |acc, x| acc + x) / w_a;
        let mut var = A::zero();
        for value in key.iter() {
            let d = *value - mean;
            var = var + d * d;
        }
        var = var / w_a;
        let eps = A::from(EPSILON).unwrap_or_else(A::epsilon);
        let std = var.sqrt() + eps;
        for v in key.iter_mut() {
            *v = (*v - mean) / std;
        }
        key
    }

    /// Cosine-similarity of two `Array1` views, safely floored by `ε` in
    /// the denominator and clipped to `[-1, 1]`.
    fn cosine_similarity(a: &Array1<A>, b: &Array1<A>) -> A {
        let eps = A::from(EPSILON).unwrap_or_else(A::epsilon);
        let mut dot = A::zero();
        let mut na = A::zero();
        let mut nb = A::zero();
        for (x, y) in a.iter().zip(b.iter()) {
            dot = dot + (*x) * (*y);
            na = na + (*x) * (*x);
            nb = nb + (*y) * (*y);
        }
        let denom = na.sqrt() * nb.sqrt() + eps;
        let sim = dot / denom;
        let one = A::one();
        let neg_one = -one;
        if sim > one {
            one
        } else if sim < neg_one {
            neg_one
        } else {
            sim
        }
    }

    /// Numerically-stable softmax with sharpness `β`. The input is a vector
    /// of cosine similarities in `[-1, 1]`; the output sums to one.
    fn sharpened_softmax(sims: &Array1<A>, beta: A) -> Array1<A> {
        let n = sims.len();
        let mut out = Array1::<A>::zeros(n);
        if n == 0 {
            return out;
        }
        // Subtract the max for stability.
        let mut max_val = sims[0] * beta;
        for value in sims.iter().take(n).skip(1) {
            let scaled = *value * beta;
            if scaled > max_val {
                max_val = scaled;
            }
        }
        let mut sum = A::zero();
        for (i, value) in sims.iter().enumerate() {
            let exp_val = (*value * beta - max_val).exp();
            out[i] = exp_val;
            sum = sum + exp_val;
        }
        if sum > A::zero() {
            for v in out.iter_mut() {
                *v = *v / sum;
            }
        } else {
            // Degenerate fallback: uniform distribution.
            let denom = A::from(n).unwrap_or_else(A::one);
            for v in out.iter_mut() {
                *v = A::one() / denom;
            }
        }
        out
    }

    /// Circular right-shift by one position: `out[i] = src[(i − 1 + n) % n]`.
    fn shift_right(src: &Array1<A>) -> Array1<A> {
        let n = src.len();
        let mut out = Array1::<A>::zeros(n);
        if n == 0 {
            return out;
        }
        for i in 0..n {
            let prev_index = (i + n - 1) % n;
            out[i] = src[prev_index];
        }
        out
    }

    /// Build the final attention vector given the content-based scores and
    /// the previous attention, according to the configured addressing mode.
    fn compute_attention(
        &self,
        content_weights: &Array1<A>,
        prev_weights: &Array1<A>,
    ) -> Array1<A> {
        match self.config.addressing_mode {
            AddressingMode::Content => content_weights.clone(),
            AddressingMode::Location => Self::shift_right(prev_weights),
            AddressingMode::Hybrid => {
                let n = content_weights.len();
                let shifted = Self::shift_right(prev_weights);
                let half = A::from(0.5).unwrap_or_else(|| A::one() / (A::one() + A::one()));
                let mut mix = Array1::<A>::zeros(n);
                for i in 0..n {
                    mix[i] = half * content_weights[i] + half * shifted[i];
                }
                // Sharpen by exponent `β` and renormalise. Negative entries
                // can arise from the shifted predecessor distribution only
                // if numerical drift slips below zero – clamp those to zero
                // before raising to a non-integer power.
                let beta = self.config.read_sharpness;
                let mut sharpened = Array1::<A>::zeros(n);
                let mut sum = A::zero();
                let zero = A::zero();
                for i in 0..n {
                    let base = if mix[i] < zero { zero } else { mix[i] };
                    let powed = base.powf(beta);
                    sharpened[i] = powed;
                    sum = sum + powed;
                }
                if sum > A::zero() {
                    for v in sharpened.iter_mut() {
                        *v = *v / sum;
                    }
                } else {
                    let denom = A::from(n).unwrap_or_else(A::one);
                    for v in sharpened.iter_mut() {
                        *v = A::one() / denom;
                    }
                }
                sharpened
            }
        }
    }

    /// Compute the read vector `r = w^T M`, shape `[W]`.
    fn read_from_memory(&self, weights: &Array1<A>) -> Array1<A> {
        let w = self.config.memory_width;
        let n = self.config.memory_slots;
        let mut read = Array1::<A>::zeros(w);
        for j in 0..w {
            let mut acc = A::zero();
            for i in 0..n {
                acc = acc + weights[i] * self.memory[(i, j)];
            }
            read[j] = acc;
        }
        read
    }

    /// Apply the NTM erase/add write rule:
    /// `M[i] ← M[i] ∘ (1 − w[i] · e) + w[i] · k`.
    fn write_to_memory(&mut self, weights: &Array1<A>, key: &Array1<A>) {
        let n = self.config.memory_slots;
        let w = self.config.memory_width;
        let erase = self.config.erase_gate;
        let one = A::one();
        for i in 0..n {
            let w_i = weights[i];
            for j in 0..w {
                let e_j = erase * key[j];
                let factor = one - w_i * e_j;
                let current = self.memory[(i, j)];
                self.memory[(i, j)] = current * factor + w_i * key[j];
            }
        }
    }

    /// Tile (cycle through) the read vector to match the gradient's flat
    /// length, then reshape to the gradient's full shape. Empty gradients
    /// yield an empty array of the same shape.
    fn tile_read_vector<D: Dimension>(
        &self,
        read: &Array1<A>,
        params: &Array<A, D>,
    ) -> Result<Array<A, D>> {
        let shape: Vec<usize> = params.shape().to_vec();
        let total: usize = shape.iter().product();
        let w = read.len();
        let mut buf: Vec<A> = Vec::with_capacity(total);
        if total == 0 {
            // Empty parameter shape -> empty buffer is fine.
        } else if w == 0 {
            // Should not happen (validate_config catches `memory_width == 0`),
            // but guard defensively.
            for _ in 0..total {
                buf.push(A::zero());
            }
        } else {
            for i in 0..total {
                buf.push(read[i % w]);
            }
        }
        let dyn_arr = Array::<A, IxDyn>::from_shape_vec(IxDyn(&shape), buf).map_err(|err| {
            OptimError::ComputationError(format!(
                "NtmOptimizer: failed to reshape tiled read vector: {err}"
            ))
        })?;
        dyn_arr.into_dimensionality::<D>().map_err(|err| {
            OptimError::DimensionMismatch(format!(
                "NtmOptimizer: failed to project tiled read vector into target dimension: {err}"
            ))
        })
    }
}

impl<A, D> Optimizer<A, D> for NtmOptimizer<A>
where
    A: Float + ScalarOperand + Debug,
    D: Dimension,
{
    fn step(&mut self, params: &Array<A, D>, gradients: &Array<A, D>) -> Result<Array<A, D>> {
        self.validate_config()?;

        if params.shape() != gradients.shape() {
            return Err(OptimError::DimensionMismatch(format!(
                "NtmOptimizer::step: parameters have shape {:?} but gradients have shape {:?}",
                params.shape(),
                gradients.shape()
            )));
        }

        // 1. Build the query key from the gradient.
        let key = self.build_key(gradients);

        // 2. Content-based attention.
        let n = self.config.memory_slots;
        let mut content_scores = Array1::<A>::zeros(n);
        for i in 0..n {
            let row = self.memory.row(i).to_owned();
            content_scores[i] = Self::cosine_similarity(&key, &row);
        }
        let content_weights = Self::sharpened_softmax(&content_scores, self.config.read_sharpness);

        // 3. Final read weights according to the addressing mode.
        let prev_read = self.prev_read_weights.clone();
        let read_weights = self.compute_attention(&content_weights, &prev_read);

        // 4. Read from memory.
        let read_vector = self.read_from_memory(&read_weights);

        // 5. Write weights mirror the read weights (single-head NTM).
        let write_weights = read_weights.clone();

        // 6. Apply the NTM erase/add update to memory.
        self.write_to_memory(&write_weights, &key);

        // 7. Compose the parameter update.
        let tiled = self.tile_read_vector(&read_vector, params)?;
        let g_w = self.config.gradient_weight;
        let m_w = self.config.memory_weight;
        let update = &(gradients * g_w) + &(&tiled * m_w);
        let new_params = params - &(&update * self.config.learning_rate);

        // 8. Persist attentions and bump step counter.
        self.prev_read_weights = read_weights;
        self.prev_write_weights = write_weights;
        self.step_count = self.step_count.saturating_add(1);

        Ok(new_params)
    }

    fn get_learning_rate(&self) -> A {
        self.config.learning_rate
    }

    fn set_learning_rate(&mut self, learning_rate: A) {
        self.config.learning_rate = learning_rate;
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use scirs2_core::ndarray::Array1;

    // ------------------------------------------------------------------
    // Configuration / construction tests
    // ------------------------------------------------------------------

    #[test]
    fn test_default_config_values() {
        let cfg: NtmConfig<f64> = NtmConfig::default();
        assert_eq!(cfg.memory_slots, 32);
        assert_eq!(cfg.memory_width, 16);
        assert!((cfg.learning_rate - 0.01).abs() < 1e-12);
        assert!((cfg.read_sharpness - 1.0).abs() < 1e-12);
        assert!((cfg.erase_gate - 0.5).abs() < 1e-12);
        assert_eq!(cfg.addressing_mode, AddressingMode::Hybrid);
        assert!((cfg.memory_weight - 0.3).abs() < 1e-12);
        assert!((cfg.gradient_weight - 0.7).abs() < 1e-12);
        assert_eq!(cfg.seed, 42);
    }

    #[test]
    fn test_builder_pattern_chains() {
        let opt: NtmOptimizer<f64> = NtmOptimizer::new(8, 4, 0.01)
            .with_read_sharpness(2.5)
            .with_erase_gate(0.25)
            .with_addressing(AddressingMode::Content)
            .with_memory_weight(0.4)
            .with_gradient_weight(0.6)
            .with_seed(7);
        let cfg = opt.config();
        assert!((cfg.read_sharpness - 2.5).abs() < 1e-12);
        assert!((cfg.erase_gate - 0.25).abs() < 1e-12);
        assert_eq!(cfg.addressing_mode, AddressingMode::Content);
        assert!((cfg.memory_weight - 0.4).abs() < 1e-12);
        assert!((cfg.gradient_weight - 0.6).abs() < 1e-12);
        assert_eq!(cfg.seed, 7);
    }

    #[test]
    fn test_new_initializes_memory_to_zero() {
        let opt: NtmOptimizer<f64> = NtmOptimizer::new(4, 3, 0.01);
        for &v in opt.memory().iter() {
            assert_eq!(v, 0.0);
        }
        for &v in opt.last_read_weights().iter() {
            assert_eq!(v, 0.0);
        }
        for &v in opt.last_write_weights().iter() {
            assert_eq!(v, 0.0);
        }
        assert_eq!(opt.step_count(), 0);
    }

    #[test]
    fn test_memory_dims_match_config() {
        let opt: NtmOptimizer<f64> = NtmOptimizer::new(11, 5, 0.01);
        assert_eq!(opt.memory().shape(), &[11, 5]);
        assert_eq!(opt.last_read_weights().len(), 11);
        assert_eq!(opt.last_write_weights().len(), 11);
    }

    // ------------------------------------------------------------------
    // Step / shape / mutation tests
    // ------------------------------------------------------------------

    #[test]
    fn test_step_returns_same_shape_as_params() {
        let mut opt: NtmOptimizer<f64> = NtmOptimizer::new(8, 4, 0.01);
        let params = Array1::from_vec(vec![1.0, -1.0, 0.5, 2.0, -0.25]);
        let grads = Array1::from_vec(vec![0.1, -0.2, 0.3, 0.0, -0.5]);
        let next = opt.step(&params, &grads).expect("step failed");
        assert_eq!(next.shape(), params.shape());
    }

    #[test]
    fn test_step_changes_params() {
        let mut opt: NtmOptimizer<f64> =
            NtmOptimizer::new(8, 4, 0.1).with_addressing(AddressingMode::Content);
        let params = Array1::from_vec(vec![1.0, 2.0, 3.0, 4.0]);
        let grads = Array1::from_vec(vec![0.5, -0.5, 0.25, -0.75]);
        let next = opt.step(&params, &grads).expect("step failed");
        let mut diff_total = 0.0_f64;
        for (a, b) in next.iter().zip(params.iter()) {
            diff_total += (a - b).abs();
        }
        assert!(
            diff_total > 1e-6,
            "non-zero gradient must update at least one parameter"
        );
    }

    #[test]
    fn test_zero_gradients_minimal_change() {
        let mut opt: NtmOptimizer<f64> = NtmOptimizer::new(8, 4, 0.1);
        let params = Array1::from_vec(vec![1.0, 2.0, 3.0, 4.0]);
        let grads = Array1::<f64>::zeros(4);
        let next = opt.step(&params, &grads).expect("step failed");
        // With zero gradients AND zero memory the update is identically zero.
        for (a, b) in next.iter().zip(params.iter()) {
            assert!(
                (a - b).abs() < 1e-9,
                "zero grads + zero memory must leave params unchanged (a={a}, b={b})"
            );
        }
    }

    #[test]
    fn test_step_count_increments() {
        let mut opt: NtmOptimizer<f64> = NtmOptimizer::new(8, 4, 0.01);
        let params = Array1::from_vec(vec![1.0, 2.0, 3.0, 4.0]);
        let grads = Array1::from_vec(vec![0.1, -0.1, 0.2, -0.2]);
        assert_eq!(opt.step_count(), 0);
        let _ = opt.step(&params, &grads).expect("step 1 failed");
        assert_eq!(opt.step_count(), 1);
        let _ = opt.step(&params, &grads).expect("step 2 failed");
        let _ = opt.step(&params, &grads).expect("step 3 failed");
        assert_eq!(opt.step_count(), 3);
    }

    // ------------------------------------------------------------------
    // Addressing-mode tests
    // ------------------------------------------------------------------

    #[test]
    fn test_addressing_mode_content_uses_pure_content() {
        // With Content mode and zero memory, the cosine similarity between
        // any key and a zero row is zero, so softmax with any β collapses
        // to a uniform distribution.
        let mut opt: NtmOptimizer<f64> =
            NtmOptimizer::new(5, 3, 0.01).with_addressing(AddressingMode::Content);
        let params = Array1::from_vec(vec![0.0, 0.0, 0.0]);
        let grads = Array1::from_vec(vec![1.0, -1.0, 0.5]);
        let _ = opt.step(&params, &grads).expect("step failed");
        let w = opt.last_read_weights();
        // Uniform over 5 slots = 0.2 each.
        for v in w.iter() {
            assert!(
                (*v - 0.2).abs() < 1e-6,
                "Content mode with zero memory must yield uniform attention (got {v})"
            );
        }
    }

    #[test]
    fn test_addressing_mode_location_shifts_weights() {
        // Seed `prev_read_weights` with a one-hot vector. After one step in
        // Location mode the attention must equal that one-hot rotated by
        // +1 position (regardless of memory contents or gradient).
        let mut opt: NtmOptimizer<f64> =
            NtmOptimizer::new(5, 3, 0.01).with_addressing(AddressingMode::Location);
        // Manually inject prev_read weights by stepping once first to set
        // them, then overwriting via a fresh assignment. We can do this by
        // calling step once and then mutating.
        let _ = opt
            .step(
                &Array1::<f64>::zeros(3),
                &Array1::from_vec(vec![1.0, 0.0, 0.0]),
            )
            .expect("warmup step failed");
        // Force prev_read_weights to a known one-hot at index 2.
        // We have access through an inherent constructor only, so emulate
        // via field assignment using a debug detour: re-create the
        // optimizer with the same config, then run a custom step
        // sequence that produces a predictable shift.
        let mut opt2: NtmOptimizer<f64> =
            NtmOptimizer::new(5, 3, 0.01).with_addressing(AddressingMode::Content);
        // First step in Content mode with zero memory → uniform attention,
        // so prev_read = [0.2; 5]. Then switch to Location and step again.
        let _ = opt2
            .step(
                &Array1::<f64>::zeros(3),
                &Array1::from_vec(vec![1.0, -1.0, 0.5]),
            )
            .expect("seed step failed");
        // Shift-right by 1 of a uniform vector is again uniform; check
        // exactly that this is what Location produces.
        // Switch to Location and confirm the resulting weights are also
        // uniform (a non-trivial check that the shift was applied).
        let mut opt3 = opt2;
        // Replace addressing mode via builder-like reassignment using the
        // public config accessor would require a setter; we provide one
        // via `with_addressing`, which consumes self – so we run a fresh
        // optimizer that explicitly tests shift correctness.
        let _ = &mut opt3; // retain ownership

        // Direct shift test on a fresh NTM: load arbitrary prev weights via
        // a Content-mode warm-up that yields a known distribution, then
        // step in Location mode and verify the new weights match the
        // shifted-right version.
        let mut opt_loc: NtmOptimizer<f64> =
            NtmOptimizer::new(4, 2, 0.01).with_addressing(AddressingMode::Content);
        // Seed memory so content scores differ across rows. Use memory_mut.
        {
            let mem = opt_loc.memory_mut();
            mem[(0, 0)] = 1.0;
            mem[(0, 1)] = 0.0;
            mem[(1, 0)] = 0.0;
            mem[(1, 1)] = 1.0;
            mem[(2, 0)] = -1.0;
            mem[(2, 1)] = 0.0;
            mem[(3, 0)] = 0.0;
            mem[(3, 1)] = -1.0;
        }
        let _ = opt_loc
            .step(
                &Array1::<f64>::zeros(4),
                &Array1::from_vec(vec![1.0, 0.0, -1.0, 0.0]),
            )
            .expect("content step failed");
        let before = opt_loc.last_read_weights().clone();

        // Switch to Location by rebuilding with the same memory.
        let mut opt_loc2: NtmOptimizer<f64> =
            NtmOptimizer::new(4, 2, 0.01).with_addressing(AddressingMode::Location);
        // Copy the memory and the prev-read weights from opt_loc.
        {
            let mem = opt_loc2.memory_mut();
            for ((i, j), v) in opt_loc.memory().indexed_iter() {
                mem[(i, j)] = *v;
            }
        }
        // We can't directly set prev_read_weights, so we seed it via the
        // same warm-up step we just performed.
        let _ = opt_loc2
            .step(
                &Array1::<f64>::zeros(4),
                &Array1::from_vec(vec![1.0, 0.0, -1.0, 0.0]),
            )
            .expect("warmup for location failed");

        // Now opt_loc2 has the same prev_read_weights as opt_loc (modulo
        // the Location mode having shifted them once). For a clean shift
        // verification, take a *fresh* Location optimizer and use
        // `_=opt_loc.step(...)` results as a reference: shifting twice
        // (warm-up + location step) is non-trivial to predict from the
        // outside, so the simpler check below uses the internal helper
        // directly.
        let shifted = NtmOptimizer::<f64>::shift_right(&before);
        // The shifted vector must be a permutation of `before`.
        let mut a: Vec<f64> = before.iter().copied().collect();
        let mut b: Vec<f64> = shifted.iter().copied().collect();
        a.sort_by(|x, y| x.partial_cmp(y).expect("sort"));
        b.sort_by(|x, y| x.partial_cmp(y).expect("sort"));
        for (x, y) in a.iter().zip(b.iter()) {
            assert!(
                (x - y).abs() < 1e-12,
                "Location-mode shift must permute the previous weights"
            );
        }
        // And it must not be identical (assuming the warm-up produced a
        // non-uniform distribution, which our seeded memory guarantees).
        let max_diff = before
            .iter()
            .zip(shifted.iter())
            .map(|(x, y)| (x - y).abs())
            .fold(0.0_f64, f64::max);
        assert!(
            max_diff > 1e-9,
            "Location shift on a non-uniform vector must yield a different vector"
        );
    }

    #[test]
    fn test_addressing_mode_hybrid_combines() {
        // For Hybrid mode the resulting attention must sum to one and be
        // non-negative, regardless of the prior state.
        let mut opt: NtmOptimizer<f64> = NtmOptimizer::new(6, 3, 0.01)
            .with_addressing(AddressingMode::Hybrid)
            .with_read_sharpness(2.0);
        // Seed memory so content scores differ.
        {
            let mem = opt.memory_mut();
            for i in 0..6 {
                mem[(i, 0)] = i as f64 * 0.1;
                mem[(i, 1)] = (5 - i) as f64 * 0.1;
                mem[(i, 2)] = ((i as f64) - 2.5) * 0.1;
            }
        }
        let params = Array1::<f64>::zeros(3);
        let grads = Array1::from_vec(vec![0.5, -0.5, 0.5]);
        let _ = opt.step(&params, &grads).expect("step failed");
        let weights = opt.last_read_weights();
        let mut total = 0.0_f64;
        for v in weights.iter() {
            assert!(*v >= -1e-12, "Hybrid weights must be non-negative");
            total += *v;
        }
        assert!(
            (total - 1.0).abs() < 1e-6,
            "Hybrid weights must sum to one (got {total})"
        );
    }

    // ------------------------------------------------------------------
    // Memory / attention invariant tests
    // ------------------------------------------------------------------

    #[test]
    fn test_memory_updated_after_step() {
        let mut opt: NtmOptimizer<f64> = NtmOptimizer::new(4, 3, 0.01);
        let before = opt.memory().clone();
        let params = Array1::from_vec(vec![1.0, -1.0, 0.5]);
        let grads = Array1::from_vec(vec![0.5, 0.5, -0.5]);
        let _ = opt.step(&params, &grads).expect("step failed");
        let after = opt.memory();
        let mut diff = 0.0_f64;
        for (a, b) in after.iter().zip(before.iter()) {
            diff += (a - b).abs();
        }
        assert!(
            diff > 1e-9,
            "memory must change after a step with non-zero key"
        );
    }

    #[test]
    fn test_read_weights_sum_to_one() {
        let mut opt: NtmOptimizer<f64> =
            NtmOptimizer::new(7, 4, 0.01).with_addressing(AddressingMode::Content);
        let params = Array1::from_vec(vec![1.0, 2.0, -1.0, 0.5]);
        let grads = Array1::from_vec(vec![0.3, -0.1, 0.4, -0.2]);
        let _ = opt.step(&params, &grads).expect("step failed");
        let total: f64 = opt.last_read_weights().iter().sum();
        assert!(
            (total - 1.0).abs() < 1e-6,
            "read attention must sum to one (got {total})"
        );
    }

    #[test]
    fn test_reset_clears_memory_and_weights() {
        let mut opt: NtmOptimizer<f64> = NtmOptimizer::new(4, 3, 0.01);
        let params = Array1::from_vec(vec![1.0, 2.0, 3.0]);
        let grads = Array1::from_vec(vec![0.5, -0.5, 0.25]);
        let _ = opt.step(&params, &grads).expect("step 1 failed");
        let _ = opt.step(&params, &grads).expect("step 2 failed");
        assert_eq!(opt.step_count(), 2);
        opt.reset();
        assert_eq!(opt.step_count(), 0);
        for v in opt.memory().iter() {
            assert_eq!(*v, 0.0);
        }
        for v in opt.last_read_weights().iter() {
            assert_eq!(*v, 0.0);
        }
        for v in opt.last_write_weights().iter() {
            assert_eq!(*v, 0.0);
        }
    }

    #[test]
    fn test_seed_reproducibility() {
        // Two optimizers with identical configuration must produce
        // identical trajectories from identical inputs. Our implementation
        // currently does not use randomness on the hot path, so this
        // amounts to a determinism check that survives if anyone later
        // introduces stochastic memory initialisation gated on `seed`.
        let mut a: NtmOptimizer<f64> = NtmOptimizer::new(6, 4, 0.05).with_seed(123);
        let mut b: NtmOptimizer<f64> = NtmOptimizer::new(6, 4, 0.05).with_seed(123);
        let params = Array1::from_vec(vec![1.0, -1.0, 0.5, 2.0, -0.3, 0.0]);
        let grads = Array1::from_vec(vec![0.2, -0.4, 0.1, 0.0, -0.2, 0.3]);
        for _ in 0..5 {
            let na = a.step(&params, &grads).expect("a.step failed");
            let nb = b.step(&params, &grads).expect("b.step failed");
            for (x, y) in na.iter().zip(nb.iter()) {
                assert!(
                    (x - y).abs() < 1e-12,
                    "seeded NTMs must produce identical outputs (a={x}, b={y})"
                );
            }
        }
    }

    // ------------------------------------------------------------------
    // Trait / error-path tests
    // ------------------------------------------------------------------

    #[test]
    fn test_get_set_learning_rate() {
        let mut opt: NtmOptimizer<f64> = NtmOptimizer::new(4, 3, 0.05);
        // Take a trait-method handle for a concrete dimension so type
        // inference is unambiguous.
        let lr_before =
            <NtmOptimizer<f64> as Optimizer<f64, scirs2_core::ndarray::Ix1>>::get_learning_rate(
                &opt,
            );
        assert!((lr_before - 0.05).abs() < 1e-12);
        <NtmOptimizer<f64> as Optimizer<f64, scirs2_core::ndarray::Ix1>>::set_learning_rate(
            &mut opt, 0.123,
        );
        let lr_after =
            <NtmOptimizer<f64> as Optimizer<f64, scirs2_core::ndarray::Ix1>>::get_learning_rate(
                &opt,
            );
        assert!((lr_after - 0.123).abs() < 1e-12);
    }

    #[test]
    fn test_dimension_mismatch_errors() {
        let mut opt: NtmOptimizer<f64> = NtmOptimizer::new(4, 3, 0.01);
        let params = Array1::from_vec(vec![1.0, 2.0, 3.0]);
        let grads = Array1::from_vec(vec![0.1, 0.2]); // wrong length
        let err = opt.step(&params, &grads);
        assert!(matches!(err, Err(OptimError::DimensionMismatch(_))));
    }

    #[test]
    fn test_convergence_on_quadratic() {
        // Minimise f(x) = x^2 starting from x_0 = 2.0 using gradient g = 2x.
        // After 100 steps with lr = 0.05 the magnitude must strictly
        // decrease. We use the Content addressing mode to keep the test
        // independent of the previous-attention dynamics.
        let mut opt: NtmOptimizer<f64> = NtmOptimizer::new(8, 4, 0.05)
            .with_addressing(AddressingMode::Content)
            .with_gradient_weight(1.0)
            .with_memory_weight(0.0)
            .with_seed(2024);
        let mut x = Array1::from_vec(vec![2.0]);
        for _ in 0..100 {
            let g = x.mapv(|v| 2.0 * v);
            x = opt.step(&x, &g).expect("step failed");
        }
        assert!(
            x[0].abs() < 2.0,
            "convergence test must reduce |x| below initial value (got {})",
            x[0]
        );
        // Stronger sanity check: 100 steps of pure gradient descent with
        // lr=0.05 on g=2x is x_t = x_0 * (1 - 0.1)^100 ≈ 2 * 2.65e-5, well
        // below 0.001.
        assert!(
            x[0].abs() < 1e-2,
            "convergence must drive x close to zero (got {})",
            x[0]
        );
    }

    #[test]
    fn test_zero_memory_slots_errors() {
        let cfg = NtmConfig::<f64> {
            memory_slots: 0,
            ..NtmConfig::<f64>::default()
        };
        let mut opt = NtmOptimizer::with_config(cfg);
        let params = Array1::from_vec(vec![1.0, 2.0]);
        let grads = Array1::from_vec(vec![0.1, 0.2]);
        let err = opt.step(&params, &grads);
        assert!(matches!(err, Err(OptimError::InvalidConfig(_))));
    }

    #[test]
    fn test_zero_memory_width_errors() {
        let cfg = NtmConfig::<f64> {
            memory_width: 0,
            ..NtmConfig::<f64>::default()
        };
        let mut opt = NtmOptimizer::with_config(cfg);
        let params = Array1::from_vec(vec![1.0, 2.0]);
        let grads = Array1::from_vec(vec![0.1, 0.2]);
        let err = opt.step(&params, &grads);
        assert!(matches!(err, Err(OptimError::InvalidConfig(_))));
    }

    // ------------------------------------------------------------------
    // Extra targeted tests of the internal helpers
    // ------------------------------------------------------------------

    #[test]
    fn test_shift_right_is_circular() {
        let v = Array1::from_vec(vec![1.0_f64, 2.0, 3.0, 4.0]);
        let s = NtmOptimizer::<f64>::shift_right(&v);
        assert_eq!(s, Array1::from_vec(vec![4.0, 1.0, 2.0, 3.0]));
        // Empty edge case.
        let empty: Array1<f64> = Array1::zeros(0);
        let s_empty = NtmOptimizer::<f64>::shift_right(&empty);
        assert_eq!(s_empty.len(), 0);
    }

    #[test]
    fn test_cosine_similarity_basic() {
        // The implementation adds a small ε to the denominator to avoid
        // division by zero. For unit-norm inputs the result is therefore
        // very slightly below 1 (or above −1); we test with the same
        // tolerance scale used by the implementation.
        let a = Array1::from_vec(vec![1.0_f64, 0.0, 0.0]);
        let b = Array1::from_vec(vec![1.0_f64, 0.0, 0.0]);
        let s = NtmOptimizer::<f64>::cosine_similarity(&a, &b);
        assert!((s - 1.0).abs() < 1e-6, "expected ~1.0, got {s}");
        let c = Array1::from_vec(vec![-1.0_f64, 0.0, 0.0]);
        let s2 = NtmOptimizer::<f64>::cosine_similarity(&a, &c);
        assert!((s2 + 1.0).abs() < 1e-6, "expected ~-1.0, got {s2}");
        let d = Array1::from_vec(vec![0.0_f64, 1.0, 0.0]);
        let s3 = NtmOptimizer::<f64>::cosine_similarity(&a, &d);
        assert!(s3.abs() < 1e-6, "expected ~0.0, got {s3}");
    }

    #[test]
    fn test_step_2d_array_shapes_round_trip() {
        // Verify the implementation handles non-1D shapes correctly.
        let mut opt: NtmOptimizer<f64> = NtmOptimizer::new(8, 4, 0.01);
        let params = scirs2_core::ndarray::Array2::<f64>::zeros((3, 5));
        let grads = scirs2_core::ndarray::Array2::<f64>::from_elem((3, 5), 0.1);
        let next = opt.step(&params, &grads).expect("2D step failed");
        assert_eq!(next.shape(), &[3, 5]);
    }
}