quantrs2_tytan/

coherent_ising_machine.rs

//! Coherent Ising Machine (CIM) simulation for quantum-inspired optimization.
//!
//! This module provides a simulation of Coherent Ising Machines, which use
//! optical parametric oscillators to solve optimization problems.

#![allow(dead_code)]

use crate::sampler::{SampleResult, Sampler, SamplerError, SamplerResult};
use scirs2_core::ndarray::{Array, Array1, Array2, IxDyn};
use scirs2_core::random::prelude::*;
use scirs2_core::random::{Distribution, RandNormal, Rng, SeedableRng};
use scirs2_core::Complex64;
use std::collections::HashMap;

type Normal<T> = RandNormal<T>;
use std::f64::consts::PI;

/// Coherent Ising Machine simulator
#[derive(Clone)]
pub struct CIMSimulator {
    /// Number of spins
    pub n_spins: usize,
    /// Pump parameter
    pump_parameter: f64,
    /// Detuning parameter
    detuning: f64,
    /// Time step for evolution
    dt: f64,
    /// Total evolution time
    evolution_time: f64,
    /// Noise strength
    noise_strength: f64,
    /// Coupling strength scaling
    coupling_scale: f64,
    /// Random seed
    seed: Option<u64>,
    /// Use measurement feedback
    use_feedback: bool,
    /// Feedback delay
    feedback_delay: f64,
}

impl CIMSimulator {
    /// Create new CIM simulator
    pub const fn new(n_spins: usize) -> Self {
        Self {
            n_spins,
            pump_parameter: 1.0,
            detuning: 0.0,
            dt: 0.01,
            evolution_time: 10.0,
            noise_strength: 0.1,
            coupling_scale: 1.0,
            seed: None,
            use_feedback: true,
            feedback_delay: 0.1,
        }
    }

    /// Set pump parameter (controls oscillation amplitude)
    pub const fn with_pump_parameter(mut self, pump: f64) -> Self {
        self.pump_parameter = pump;
        self
    }

    /// Set detuning (frequency offset)
    pub const fn with_detuning(mut self, detuning: f64) -> Self {
        self.detuning = detuning;
        self
    }

    /// Set time step
    pub const fn with_time_step(mut self, dt: f64) -> Self {
        self.dt = dt;
        self
    }

    /// Set evolution time
    pub const fn with_evolution_time(mut self, time: f64) -> Self {
        self.evolution_time = time;
        self
    }

    /// Set noise strength
    pub const fn with_noise_strength(mut self, noise: f64) -> Self {
        self.noise_strength = noise;
        self
    }

    /// Set coupling scale
    pub const fn with_coupling_scale(mut self, scale: f64) -> Self {
        self.coupling_scale = scale;
        self
    }

    /// Set random seed
    pub const fn with_seed(mut self, seed: u64) -> Self {
        self.seed = Some(seed);
        self
    }

    /// Enable/disable measurement feedback
    pub const fn with_feedback(mut self, use_feedback: bool) -> Self {
        self.use_feedback = use_feedback;
        self
    }

    /// Simulate CIM evolution
    fn simulate_cim(
        &self,
        coupling_matrix: &Array2<f64>,
        local_fields: &Array1<f64>,
        rng: &mut StdRng,
    ) -> Result<Vec<f64>, String> {
        let n = self.n_spins;
        let steps = (self.evolution_time / self.dt) as usize;

        // Initialize oscillator amplitudes (complex)
        let mut amplitudes: Vec<Complex64> = (0..n)
            .map(|_| {
                let r = rng.random_range(0.0..0.1);
                let theta = rng.random_range(0.0..2.0 * PI);
                Complex64::new(r * theta.cos(), r * theta.sin())
            })
            .collect();

        // Normal distribution for standard Gaussian samples (scaled by noise_strength at use site)
        let standard_normal = Normal::<f64>::new(0.0_f64, 1.0_f64)
            .expect("Normal distribution with mean=0, std=1 is always valid");

        // Evolution loop
        for step in 0..steps {
            let mut new_amplitudes = amplitudes.clone();

            // Construct per-step noise distribution outside inner spin loop for efficiency
            let noise_scale = self.noise_strength;

            for i in 0..n {
                // Compute coupling term
                let mut coupling_term = Complex64::new(0.0, 0.0);
                for j in 0..n {
                    if i != j {
                        let coupling = coupling_matrix[[i, j]] * self.coupling_scale;

                        if self.use_feedback {
                            // Measurement feedback with delay
                            let delayed_step =
                                (step as f64 - self.feedback_delay / self.dt).max(0.0) as usize;
                            let delayed_amp = if delayed_step < step {
                                amplitudes[j]
                            } else {
                                amplitudes[j]
                            };
                            coupling_term += coupling * delayed_amp.re;
                        } else {
                            // Direct coupling
                            coupling_term += coupling * amplitudes[j];
                        }
                    }
                }

                // Add local field
                coupling_term += local_fields[i];

                // Nonlinear evolution equation
                let nonlinear_term = amplitudes[i] * amplitudes[i].norm_sqr();
                let pump_term = self.pump_parameter;
                let detuning_term = Complex64::new(0.0, -self.detuning) * amplitudes[i];

                // Stochastic differential equation
                let deterministic = (pump_term - 1.0) * amplitudes[i] - nonlinear_term
                    + detuning_term
                    + coupling_term;

                // Gaussian noise injection for the SDE Wiener process increments
                // Each spin receives independent real and imaginary noise components
                let noise_re = standard_normal.sample(rng) * noise_scale;
                let noise_im = standard_normal.sample(rng) * noise_scale;
                let noise = Complex64::new(noise_re, noise_im);

                // Update amplitude using Euler-Maruyama scheme: dA = f(A)dt + g*dW
                // where dW ~ sqrt(dt) * N(0,1)
                new_amplitudes[i] =
                    amplitudes[i] + self.dt * deterministic + (self.dt.sqrt()) * noise;
            }

            amplitudes = new_amplitudes;

            // Optional: apply normalization or constraints
            if step % 100 == 0 {
                self.apply_constraints(&mut amplitudes);
            }
        }

        // Extract final spin configuration
        let spins: Vec<f64> = amplitudes.iter().map(|amp| amp.re.signum()).collect();

        Ok(spins)
    }

    /// Apply constraints to maintain physical behavior
    fn apply_constraints(&self, amplitudes: &mut Vec<Complex64>) {
        // Saturation constraint
        let max_amplitude = 2.0;
        for amp in amplitudes.iter_mut() {
            if amp.norm() > max_amplitude {
                *amp = *amp / amp.norm() * max_amplitude;
            }
        }
    }

    /// Convert QUBO to Ising model
    fn qubo_to_ising(&self, qubo_matrix: &Array2<f64>) -> (Array2<f64>, Array1<f64>, f64) {
        let n = qubo_matrix.shape()[0];
        let mut j_matrix = Array2::zeros((n, n));
        let mut h_vector = Array1::zeros(n);
        let mut offset = 0.0;

        // Convert QUBO to Ising: s_i = 2*x_i - 1
        for i in 0..n {
            for j in 0..n {
                if i == j {
                    h_vector[i] += qubo_matrix[[i, i]];
                    offset += qubo_matrix[[i, i]] / 2.0;
                } else if i < j {
                    j_matrix[[i, j]] = qubo_matrix[[i, j]] / 4.0;
                    j_matrix[[j, i]] = qubo_matrix[[i, j]] / 4.0;
                    h_vector[i] += qubo_matrix[[i, j]] / 2.0;
                    h_vector[j] += qubo_matrix[[i, j]] / 2.0;
                    offset += qubo_matrix[[i, j]] / 4.0;
                }
            }
        }

        (j_matrix, h_vector, offset)
    }

    /// Convert Ising spins to binary variables
    fn spins_to_binary(&self, spins: &[f64]) -> Vec<bool> {
        spins.iter().map(|&s| s > 0.0).collect()
    }

    /// Calculate Ising energy
    fn calculate_ising_energy(
        &self,
        spins: &[f64],
        j_matrix: &Array2<f64>,
        h_vector: &Array1<f64>,
    ) -> f64 {
        let n = spins.len();
        let mut energy = 0.0;

        // Quadratic terms
        for i in 0..n {
            for j in i + 1..n {
                energy += j_matrix[[i, j]] * spins[i] * spins[j];
            }
        }

        // Linear terms
        for i in 0..n {
            energy += h_vector[i] * spins[i];
        }

        energy
    }
}

impl Sampler for CIMSimulator {
    fn run_qubo(
        &self,
        qubo: &(Array2<f64>, HashMap<String, usize>),
        shots: usize,
    ) -> SamplerResult<Vec<SampleResult>> {
        let (qubo_matrix, var_map) = qubo;
        let n = qubo_matrix.shape()[0];

        if n != self.n_spins {
            return Err(SamplerError::InvalidParameter(format!(
                "CIM configured for {} spins but QUBO has {} variables",
                self.n_spins, n
            )));
        }

        // Convert QUBO to Ising
        let (j_matrix, h_vector, offset) = self.qubo_to_ising(qubo_matrix);

        // Initialize RNG
        let mut rng = match self.seed {
            Some(seed) => StdRng::seed_from_u64(seed),
            None => StdRng::seed_from_u64(42), // Simple fallback for thread RNG
        };

        let mut results = Vec::new();
        let mut solution_counts: HashMap<Vec<bool>, (f64, usize)> = HashMap::new();

        // Run multiple shots
        for _ in 0..shots {
            // Simulate CIM
            let spins = self.simulate_cim(&j_matrix, &h_vector, &mut rng)?;

            // Convert to binary
            let binary = self.spins_to_binary(&spins);

            // Calculate energy
            let ising_energy = self.calculate_ising_energy(&spins, &j_matrix, &h_vector);
            let qubo_energy = ising_energy + offset;

            // Count occurrences
            let entry = solution_counts
                .entry(binary.clone())
                .or_insert((qubo_energy, 0));
            entry.1 += 1;
        }

        // Convert to sample results
        for (binary, (energy, count)) in solution_counts {
            let assignments: HashMap<String, bool> = var_map
                .iter()
                .map(|(var, &idx)| (var.clone(), binary[idx]))
                .collect();

            results.push(SampleResult {
                assignments,
                energy,
                occurrences: count,
            });
        }

        // Sort by energy (NaN values are treated as equal)
        results.sort_by(|a, b| {
            a.energy
                .partial_cmp(&b.energy)
                .unwrap_or(std::cmp::Ordering::Equal)
        });

        Ok(results)
    }

    fn run_hobo(
        &self,
        hobo: &(Array<f64, IxDyn>, HashMap<String, usize>),
        shots: usize,
    ) -> SamplerResult<Vec<SampleResult>> {
        let (tensor, var_map) = hobo;

        // Quadratize HOBO → QUBO via Rosenberg polynomial reduction.
        let (qubo, ext_var_map) = crate::sampler::energy::hobo_to_qubo(tensor, var_map)
            .map_err(SamplerError::InvalidModel)?;

        // The quadratized QUBO may have more variables than the CIM was
        // originally configured for.  We create a temporary CIM with the
        // correct spin count rather than rejecting the problem.
        let n_qubo = qubo.shape()[0];
        let mut tmp_cim = self.clone();
        tmp_cim.n_spins = n_qubo;

        // Delegate to the QUBO solver.
        let mut results = tmp_cim.run_qubo(&(qubo, ext_var_map), shots)?;

        // Strip auxiliary variables (introduced by quadratization) from each result.
        for result in &mut results {
            result.assignments.retain(|k, _| !k.starts_with("_aux_"));
        }

        Ok(results)
    }
}

/// Advanced CIM with pulse shaping and error correction
pub struct AdvancedCIM {
    /// Base CIM simulator
    pub base_cim: CIMSimulator,
    /// Pulse shaping parameters
    pulse_shape: PulseShape,
    /// Error correction scheme
    error_correction: ErrorCorrectionScheme,
    /// Bifurcation control
    pub bifurcation_control: BifurcationControl,
    /// Multi-round iterations
    pub num_rounds: usize,
}

#[derive(Debug, Clone)]
pub enum PulseShape {
    /// Gaussian pulse
    Gaussian { width: f64, amplitude: f64 },
    /// Hyperbolic secant pulse
    Sech { width: f64, amplitude: f64 },
    /// Custom pulse function
    Custom { name: String, parameters: Vec<f64> },
}

#[derive(Debug, Clone)]
pub enum ErrorCorrectionScheme {
    /// No error correction
    None,
    /// Majority voting
    MajorityVoting { window_size: usize },
    /// Parity check
    ParityCheck { check_matrix: Array2<bool> },
    /// Stabilizer codes
    Stabilizer { generators: Vec<Vec<bool>> },
}

#[derive(Debug, Clone)]
pub struct BifurcationControl {
    /// Initial bifurcation parameter
    pub initial_param: f64,
    /// Final bifurcation parameter
    pub final_param: f64,
    /// Ramp time
    ramp_time: f64,
    /// Ramp function type
    ramp_type: RampType,
}

#[derive(Debug, Clone)]
pub enum RampType {
    Linear,
    Exponential,
    Sigmoid,
    Adaptive,
}

impl AdvancedCIM {
    /// Create new advanced CIM
    pub const fn new(n_spins: usize) -> Self {
        Self {
            base_cim: CIMSimulator::new(n_spins),
            pulse_shape: PulseShape::Gaussian {
                width: 1.0,
                amplitude: 1.0,
            },
            error_correction: ErrorCorrectionScheme::None,
            bifurcation_control: BifurcationControl {
                initial_param: 0.0,
                final_param: 2.0,
                ramp_time: 5.0,
                ramp_type: RampType::Linear,
            },
            num_rounds: 1,
        }
    }

    /// Set pulse shape
    pub fn with_pulse_shape(mut self, shape: PulseShape) -> Self {
        self.pulse_shape = shape;
        self
    }

    /// Set error correction
    pub fn with_error_correction(mut self, scheme: ErrorCorrectionScheme) -> Self {
        self.error_correction = scheme;
        self
    }

    /// Set bifurcation control
    pub const fn with_bifurcation_control(mut self, control: BifurcationControl) -> Self {
        self.bifurcation_control = control;
        self
    }

    /// Set number of rounds
    pub const fn with_num_rounds(mut self, rounds: usize) -> Self {
        self.num_rounds = rounds;
        self
    }

    /// Apply pulse shaping to pump
    fn apply_pulse_shaping(&self, t: f64) -> f64 {
        match &self.pulse_shape {
            PulseShape::Gaussian { width, amplitude } => {
                let sigma = width;
                amplitude * (-t * t / (2.0 * sigma * sigma)).exp()
            }
            PulseShape::Sech { width, amplitude } => amplitude / (t / width).cosh(),
            PulseShape::Custom { .. } => {
                // Custom implementation
                1.0
            }
        }
    }

    /// Apply error correction
    fn apply_error_correction(&self, spins: &mut Vec<f64>, history: &[Vec<f64>]) {
        match &self.error_correction {
            ErrorCorrectionScheme::None => {}
            ErrorCorrectionScheme::MajorityVoting { window_size } => {
                if history.len() >= *window_size {
                    for i in 0..spins.len() {
                        let mut sum = 0.0;
                        for h in history.iter().rev().take(*window_size) {
                            sum += h[i];
                        }
                        spins[i] = if sum > 0.0 { 1.0 } else { -1.0 };
                    }
                }
            }
            ErrorCorrectionScheme::ParityCheck { check_matrix } => {
                // Implement parity check correction
                let n = spins.len();
                let m = check_matrix.shape()[0];

                for i in 0..m {
                    let mut parity = 0;
                    for j in 0..n {
                        if check_matrix[[i, j]] && spins[j] > 0.0 {
                            parity ^= 1;
                        }
                    }
                    // Correct if parity check fails
                    if parity != 0 {
                        // Find minimum weight correction
                        // Simplified: flip first spin in syndrome
                        for j in 0..n {
                            if check_matrix[[i, j]] {
                                spins[j] *= -1.0;
                                break;
                            }
                        }
                    }
                }
            }
            ErrorCorrectionScheme::Stabilizer { .. } => {
                // Stabilizer code implementation
            }
        }
    }

    /// Compute bifurcation parameter
    fn compute_bifurcation_param(&self, t: f64) -> f64 {
        let progress = (t / self.bifurcation_control.ramp_time).min(1.0);
        let initial = self.bifurcation_control.initial_param;
        let final_param = self.bifurcation_control.final_param;

        match self.bifurcation_control.ramp_type {
            RampType::Linear => (final_param - initial).mul_add(progress, initial),
            RampType::Exponential => {
                (final_param - initial).mul_add(1.0 - (-5.0 * progress).exp(), initial)
            }
            RampType::Sigmoid => {
                let x = 10.0 * (progress - 0.5);
                let sigmoid = 1.0 / (1.0 + (-x).exp());
                (final_param - initial).mul_add(sigmoid, initial)
            }
            RampType::Adaptive => {
                // Adaptive based on convergence
                (final_param - initial).mul_add(progress.powi(2), initial)
            }
        }
    }
}

impl Sampler for AdvancedCIM {
    fn run_qubo(
        &self,
        qubo: &(Array2<f64>, HashMap<String, usize>),
        shots: usize,
    ) -> SamplerResult<Vec<SampleResult>> {
        let mut all_results = Vec::new();
        let shots_per_round = shots / self.num_rounds.max(1);

        for round in 0..self.num_rounds {
            // Update pump parameter based on bifurcation control
            let t = round as f64 * self.base_cim.evolution_time / self.num_rounds as f64;
            let pump = self.compute_bifurcation_param(t);

            let mut round_cim = self.base_cim.clone();
            round_cim.pump_parameter = pump * self.apply_pulse_shaping(t);

            // Run CIM for this round
            let round_results = round_cim.run_qubo(qubo, shots_per_round)?;
            all_results.extend(round_results);
        }

        // Aggregate and sort results
        let mut aggregated: HashMap<Vec<bool>, (f64, usize)> = HashMap::new();

        for result in all_results {
            let state: Vec<bool> = qubo.1.keys().map(|var| result.assignments[var]).collect();

            let entry = aggregated.entry(state).or_insert((result.energy, 0));
            entry.1 += result.occurrences;
        }

        let mut final_results: Vec<SampleResult> = aggregated
            .into_iter()
            .map(|(state, (energy, count))| {
                let assignments: HashMap<String, bool> = qubo
                    .1
                    .iter()
                    .zip(state.iter())
                    .map(|((var, _), &val)| (var.clone(), val))
                    .collect();

                SampleResult {
                    assignments,
                    energy,
                    occurrences: count,
                }
            })
            .collect();

        final_results.sort_by(|a, b| {
            a.energy
                .partial_cmp(&b.energy)
                .unwrap_or(std::cmp::Ordering::Equal)
        });

        Ok(final_results)
    }

    fn run_hobo(
        &self,
        hobo: &(Array<f64, IxDyn>, HashMap<String, usize>),
        shots: usize,
    ) -> SamplerResult<Vec<SampleResult>> {
        self.base_cim.run_hobo(hobo, shots)
    }
}

/// Network of coupled CIM modules for large-scale problems
pub struct NetworkedCIM {
    /// Individual CIM modules
    pub modules: Vec<CIMSimulator>,
    /// Inter-module coupling topology
    topology: NetworkTopology,
    /// Synchronization scheme
    sync_scheme: SynchronizationScheme,
    /// Communication delay
    comm_delay: f64,
}

#[derive(Debug, Clone)]
pub enum NetworkTopology {
    /// All-to-all coupling
    FullyConnected,
    /// Ring topology
    Ring,
    /// 2D grid
    Grid2D { rows: usize, cols: usize },
    /// Hierarchical
    Hierarchical { levels: usize },
    /// Custom adjacency
    Custom { adjacency: Array2<bool> },
}

#[derive(Debug, Clone)]
pub enum SynchronizationScheme {
    /// Synchronous updates
    Synchronous,
    /// Asynchronous with random order
    Asynchronous,
    /// Block synchronous
    BlockSynchronous { block_size: usize },
    /// Event-driven
    EventDriven { threshold: f64 },
}

impl NetworkedCIM {
    /// Create new networked CIM
    pub fn new(num_modules: usize, spins_per_module: usize, topology: NetworkTopology) -> Self {
        let modules = (0..num_modules)
            .map(|_| CIMSimulator::new(spins_per_module))
            .collect();

        Self {
            modules,
            topology,
            sync_scheme: SynchronizationScheme::Synchronous,
            comm_delay: 0.0,
        }
    }

    /// Set synchronization scheme
    pub const fn with_sync_scheme(mut self, scheme: SynchronizationScheme) -> Self {
        self.sync_scheme = scheme;
        self
    }

    /// Set communication delay
    pub const fn with_comm_delay(mut self, delay: f64) -> Self {
        self.comm_delay = delay;
        self
    }

    /// Get module neighbors based on topology
    pub fn get_neighbors(&self, module_idx: usize) -> Vec<usize> {
        match &self.topology {
            NetworkTopology::FullyConnected => (0..self.modules.len())
                .filter(|&i| i != module_idx)
                .collect(),
            NetworkTopology::Ring => {
                let n = self.modules.len();
                vec![(module_idx + n - 1) % n, (module_idx + 1) % n]
            }
            NetworkTopology::Grid2D { rows, cols } => {
                let row = module_idx / cols;
                let col = module_idx % cols;
                let mut neighbors = Vec::new();

                if row > 0 {
                    neighbors.push((row - 1) * cols + col);
                }
                if row < rows - 1 {
                    neighbors.push((row + 1) * cols + col);
                }
                if col > 0 {
                    neighbors.push(row * cols + (col - 1));
                }
                if col < cols - 1 {
                    neighbors.push(row * cols + (col + 1));
                }

                neighbors
            }
            _ => Vec::new(),
        }
    }
}

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

    #[test]
    fn test_cim_simulator() {
        let cim = CIMSimulator::new(4)
            .with_pump_parameter(1.5)
            .with_evolution_time(5.0)
            .with_seed(42);

        // Create simple QUBO
        let mut qubo_matrix = Array2::zeros((4, 4));
        qubo_matrix[[0, 1]] = -1.0;
        qubo_matrix[[1, 0]] = -1.0;

        let mut var_map = HashMap::new();
        for i in 0..4 {
            var_map.insert(format!("x{i}"), i);
        }

        let results = cim
            .run_qubo(&(qubo_matrix, var_map), 10)
            .expect("CIM run_qubo should succeed for valid QUBO input");
        assert!(!results.is_empty());
    }

    #[test]
    fn test_advanced_cim() {
        let cim = AdvancedCIM::new(3)
            .with_pulse_shape(PulseShape::Gaussian {
                width: 1.0,
                amplitude: 1.5,
            })
            .with_num_rounds(2);

        assert_eq!(cim.num_rounds, 2);
    }

    #[test]
    fn test_networked_cim() {
        let net_cim = NetworkedCIM::new(4, 2, NetworkTopology::Ring)
            .with_sync_scheme(SynchronizationScheme::Synchronous);

        assert_eq!(net_cim.modules.len(), 4);
        assert_eq!(net_cim.get_neighbors(0), vec![3, 1]);
    }
}