quantrs2_anneal/

chain_break.rs

//! Chain break resolution algorithms for quantum annealing
//!
//! When logical variables are embedded onto physical qubits using chains,
//! the physical qubits in a chain may disagree in the solution. This module
//! provides algorithms to resolve these chain breaks.

use crate::embedding::Embedding;
use crate::ising::{IsingError, IsingResult};
use scirs2_core::random::{thread_rng, ChaCha8Rng, Rng, RngExt, SeedableRng};
use std::collections::{HashMap, HashSet};

/// Represents a solution from quantum annealing hardware
#[derive(Debug, Clone)]
pub struct HardwareSolution {
    /// Values of physical qubits (spin values: +1 or -1)
    pub spins: Vec<i8>,
    /// Energy of this solution
    pub energy: f64,
    /// Number of occurrences (for multiple reads)
    pub occurrences: usize,
}

/// Resolved solution after chain break resolution
#[derive(Debug, Clone)]
pub struct ResolvedSolution {
    /// Values of logical variables
    pub logical_spins: Vec<i8>,
    /// Number of broken chains
    pub chain_breaks: usize,
    /// Energy of the resolved *logical* configuration.
    ///
    /// When a [`LogicalProblem`] is supplied to the resolver this is the exact
    /// energy of `logical_spins` under that problem. When no logical problem is
    /// available the resolver cannot recompute the logical energy, so it reports
    /// the measured hardware energy (a real quantity) as the best available
    /// estimate.
    pub energy: f64,
    /// Original hardware solution
    pub hardware_solution: HardwareSolution,
}

/// Chain break resolution method
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum ResolutionMethod {
    /// Take majority vote within each chain
    MajorityVote,
    /// Minimize energy of the logical problem
    EnergyMinimization,
    /// Use weighted majority based on coupling strengths
    WeightedMajority,
    /// Discard solutions with broken chains
    Discard,
}

/// Chain break resolver
pub struct ChainBreakResolver {
    /// Resolution method to use
    pub method: ResolutionMethod,
    /// Tie-breaking strategy for majority vote
    pub tie_break_random: bool,
    /// Random seed for tie-breaking
    pub seed: Option<u64>,
}

impl Default for ChainBreakResolver {
    fn default() -> Self {
        Self {
            method: ResolutionMethod::MajorityVote,
            tie_break_random: true,
            seed: None,
        }
    }
}

impl ChainBreakResolver {
    /// Resolve chain breaks in a single hardware solution
    pub fn resolve_solution(
        &self,
        hardware_solution: &HardwareSolution,
        embedding: &Embedding,
        logical_problem: Option<&LogicalProblem>,
    ) -> IsingResult<ResolvedSolution> {
        match self.method {
            ResolutionMethod::MajorityVote => {
                self.resolve_majority_vote(hardware_solution, embedding, logical_problem)
            }
            ResolutionMethod::WeightedMajority => {
                self.resolve_weighted_majority(hardware_solution, embedding, logical_problem)
            }
            ResolutionMethod::EnergyMinimization => {
                let problem = logical_problem.ok_or_else(|| {
                    IsingError::InvalidValue(
                        "Energy minimization requires logical problem".to_string(),
                    )
                })?;
                self.resolve_energy_minimization(hardware_solution, embedding, problem)
            }
            ResolutionMethod::Discard => self.resolve_discard(hardware_solution, embedding),
        }
    }

    /// Resolve multiple hardware solutions
    pub fn resolve_solutions(
        &self,
        hardware_solutions: &[HardwareSolution],
        embedding: &Embedding,
        logical_problem: Option<&LogicalProblem>,
    ) -> IsingResult<Vec<ResolvedSolution>> {
        let mut resolved = Vec::new();

        for hw_solution in hardware_solutions {
            match self.resolve_solution(hw_solution, embedding, logical_problem) {
                Ok(solution) => resolved.push(solution),
                Err(_) if self.method == ResolutionMethod::Discard => {
                    // Skip broken solutions when using discard method
                    continue;
                }
                Err(e) => return Err(e),
            }
        }

        // Sort by energy
        resolved.sort_by(|a, b| a.energy.total_cmp(&b.energy));

        Ok(resolved)
    }

    /// Resolve using majority vote
    fn resolve_majority_vote(
        &self,
        hardware_solution: &HardwareSolution,
        embedding: &Embedding,
        logical_problem: Option<&LogicalProblem>,
    ) -> IsingResult<ResolvedSolution> {
        let mut logical_spins = Vec::new();
        let mut chain_breaks = 0;
        let num_vars = embedding.chains.len();
        let mut rng = self.tie_break_rng();

        for var in 0..num_vars {
            let chain = embedding
                .chains
                .get(&var)
                .ok_or_else(|| IsingError::InvalidQubit(var))?;

            // Count votes
            let mut plus_votes = 0;
            let mut minus_votes = 0;

            for &qubit in chain {
                if qubit >= hardware_solution.spins.len() {
                    return Err(IsingError::InvalidQubit(qubit));
                }

                match hardware_solution.spins[qubit] {
                    1 => plus_votes += 1,
                    -1 => minus_votes += 1,
                    _ => return Err(IsingError::InvalidValue("Invalid spin value".to_string())),
                }
            }

            // Determine logical value
            let logical_value = if plus_votes > minus_votes {
                1
            } else if minus_votes > plus_votes {
                -1
            } else {
                // Genuine tie: break it with the (optionally seeded) RNG when
                // random tie-breaking is enabled, otherwise default to +1.
                self.break_tie(&mut rng)
            };

            // Check for chain breaks
            let unanimous = plus_votes == 0 || minus_votes == 0;
            if !unanimous {
                chain_breaks += 1;
            }

            logical_spins.push(logical_value);
        }

        let energy = Self::resolved_energy(&logical_spins, hardware_solution, logical_problem);

        Ok(ResolvedSolution {
            logical_spins,
            chain_breaks,
            energy,
            hardware_solution: hardware_solution.clone(),
        })
    }

    /// Resolve using weighted majority based on coupling strengths
    fn resolve_weighted_majority(
        &self,
        hardware_solution: &HardwareSolution,
        embedding: &Embedding,
        logical_problem: Option<&LogicalProblem>,
    ) -> IsingResult<ResolvedSolution> {
        // Weighted majority voting: weight each qubit's vote by the number of
        // other qubits in the chain that agree with it. This gives more influence
        // to qubits that are part of a larger consensus.

        let num_vars = embedding.chains.len();
        let mut logical_spins = vec![0i8; num_vars];
        let mut chain_breaks = 0;
        let mut rng = self.tie_break_rng();

        for var in 0..num_vars {
            if let Some(chain) = embedding.chains.get(&var) {
                if chain.is_empty() {
                    return Err(IsingError::InvalidValue(format!(
                        "Empty chain for variable {var}"
                    )));
                }

                if chain.len() == 1 {
                    // Single qubit chain - no possibility of chain break
                    logical_spins[var] = hardware_solution.spins[chain[0]];
                    continue;
                }

                // Calculate weighted votes for +1 and -1
                let mut weight_plus = 0.0;
                let mut weight_minus = 0.0;
                let mut has_disagreement = false;

                for &qubit_i in chain {
                    let spin_i = hardware_solution.spins[qubit_i];

                    // Calculate weight: count how many qubits in the chain agree with this one
                    let mut agreement_count = 0.0;
                    for &qubit_j in chain {
                        if qubit_i != qubit_j && hardware_solution.spins[qubit_j] == spin_i {
                            agreement_count += 1.0;
                        }
                    }

                    // Weight is: 1.0 (base) + agreement_count (bonus for consensus)
                    let weight = 1.0 + agreement_count;

                    if spin_i == 1 {
                        weight_plus += weight;
                    } else if spin_i == -1 {
                        weight_minus += weight;
                    }

                    // Check for disagreement
                    if hardware_solution.spins[chain[0]] != spin_i {
                        has_disagreement = true;
                    }
                }

                // Choose the spin value with higher weighted vote
                if weight_plus > weight_minus {
                    logical_spins[var] = 1;
                } else if weight_minus > weight_plus {
                    logical_spins[var] = -1;
                } else {
                    // Genuine tie: break it with the (optionally seeded) RNG
                    // when random tie-breaking is enabled, otherwise fall back
                    // to the first qubit's measured spin.
                    logical_spins[var] = if self.tie_break_random {
                        self.break_tie(&mut rng)
                    } else {
                        hardware_solution.spins[chain[0]]
                    };
                }

                if has_disagreement {
                    chain_breaks += 1;
                }
            }
        }

        let energy = Self::resolved_energy(&logical_spins, hardware_solution, logical_problem);

        Ok(ResolvedSolution {
            logical_spins,
            chain_breaks,
            energy,
            hardware_solution: hardware_solution.clone(),
        })
    }

    /// Resolve by minimizing energy of logical problem
    fn resolve_energy_minimization(
        &self,
        hardware_solution: &HardwareSolution,
        embedding: &Embedding,
        logical_problem: &LogicalProblem,
    ) -> IsingResult<ResolvedSolution> {
        let mut resolved =
            self.resolve_majority_vote(hardware_solution, embedding, Some(logical_problem))?;

        // For each broken chain, try flipping the logical variable
        for var in 0..resolved.logical_spins.len() {
            if self.is_chain_broken(var, hardware_solution, embedding)? {
                // Calculate energy with current value
                let current_energy = logical_problem.calculate_energy(&resolved.logical_spins);

                // Flip and calculate energy
                resolved.logical_spins[var] *= -1;
                let flipped_energy = logical_problem.calculate_energy(&resolved.logical_spins);

                // Keep the flip if it lowers energy
                if flipped_energy >= current_energy {
                    resolved.logical_spins[var] *= -1; // Flip back
                }
            }
        }

        // Recalculate final energy
        resolved.energy = logical_problem.calculate_energy(&resolved.logical_spins);

        Ok(resolved)
    }

    /// Discard solutions with broken chains
    fn resolve_discard(
        &self,
        hardware_solution: &HardwareSolution,
        embedding: &Embedding,
    ) -> IsingResult<ResolvedSolution> {
        let resolved = self.resolve_majority_vote(hardware_solution, embedding, None)?;

        if resolved.chain_breaks > 0 {
            Err(IsingError::HardwareConstraint(format!(
                "Solution has {} broken chains",
                resolved.chain_breaks
            )))
        } else {
            Ok(resolved)
        }
    }

    /// Check if a chain is broken
    fn is_chain_broken(
        &self,
        var: usize,
        hardware_solution: &HardwareSolution,
        embedding: &Embedding,
    ) -> IsingResult<bool> {
        let chain = embedding
            .chains
            .get(&var)
            .ok_or_else(|| IsingError::InvalidQubit(var))?;

        if chain.is_empty() {
            return Ok(false);
        }

        let first_spin = hardware_solution.spins[chain[0]];

        for &qubit in &chain[1..] {
            if hardware_solution.spins[qubit] != first_spin {
                return Ok(true);
            }
        }

        Ok(false)
    }

    /// Construct the RNG used for tie-breaking.
    ///
    /// When a `seed` is configured the generator is deterministic (reproducible
    /// runs); otherwise it is seeded from the thread RNG so ties are broken
    /// genuinely at random.
    fn tie_break_rng(&self) -> ChaCha8Rng {
        match self.seed {
            Some(seed) => ChaCha8Rng::seed_from_u64(seed),
            None => ChaCha8Rng::seed_from_u64(thread_rng().random()),
        }
    }

    /// Break a vote tie. With `tie_break_random` enabled this draws a genuine
    /// random spin from `rng`; otherwise it deterministically defaults to `+1`.
    fn break_tie(&self, rng: &mut ChaCha8Rng) -> i8 {
        if self.tie_break_random {
            if rng.random_bool(0.5) {
                1
            } else {
                -1
            }
        } else {
            1
        }
    }

    /// Energy of the resolved logical configuration.
    ///
    /// If a [`LogicalProblem`] is available the exact logical energy is
    /// computed from the resolved spins. Otherwise the measured hardware energy
    /// is returned, since the logical energy is not derivable without the
    /// problem definition.
    fn resolved_energy(
        logical_spins: &[i8],
        hardware_solution: &HardwareSolution,
        logical_problem: Option<&LogicalProblem>,
    ) -> f64 {
        match logical_problem {
            Some(problem) => problem.calculate_energy(logical_spins),
            None => hardware_solution.energy,
        }
    }
}

/// Represents a logical problem (QUBO or Ising)
#[derive(Debug, Clone)]
pub struct LogicalProblem {
    /// Linear coefficients (`h_i` in Ising, diagonal in QUBO)
    pub linear: Vec<f64>,
    /// Quadratic coefficients as adjacency list
    pub quadratic: HashMap<(usize, usize), f64>,
    /// Constant offset
    pub offset: f64,
}

impl LogicalProblem {
    /// Create a new logical problem
    #[must_use]
    pub fn new(num_vars: usize) -> Self {
        Self {
            linear: vec![0.0; num_vars],
            quadratic: HashMap::new(),
            offset: 0.0,
        }
    }

    /// Calculate energy for a given spin configuration
    #[must_use]
    pub fn calculate_energy(&self, spins: &[i8]) -> f64 {
        let mut energy = self.offset;

        // Linear terms
        for (i, &h) in self.linear.iter().enumerate() {
            if i < spins.len() {
                energy += h * f64::from(spins[i]);
            }
        }

        // Quadratic terms
        for (&(i, j), &J) in &self.quadratic {
            if i < spins.len() && j < spins.len() {
                energy += J * f64::from(spins[i]) * f64::from(spins[j]);
            }
        }

        energy
    }

    /// Convert from QUBO to Ising representation
    pub fn from_qubo(qubo_matrix: &[Vec<f64>], offset: f64) -> IsingResult<Self> {
        let n = qubo_matrix.len();
        let mut problem = Self::new(n);
        problem.offset = offset;

        // Convert QUBO Q_ij to Ising h_i and J_ij
        // x_i = (s_i + 1) / 2
        // Minimize x^T Q x becomes minimize sum_i h_i s_i + sum_{i<j} J_ij s_i s_j

        for i in 0..n {
            for j in i..n {
                let q_ij = qubo_matrix[i][j];
                if q_ij.abs() > 1e-10 {
                    problem.offset += q_ij / 4.0;
                    if i == j {
                        // Diagonal term contributes to linear coefficient
                        problem.linear[i] += q_ij / 2.0;
                    } else {
                        // Off-diagonal term
                        problem.quadratic.insert((i, j), q_ij / 4.0);
                        problem.linear[i] += q_ij / 4.0;
                        problem.linear[j] += q_ij / 4.0;
                    }
                }
            }
        }

        Ok(problem)
    }
}

/// Chain strength optimizer
pub struct ChainStrengthOptimizer {
    /// Minimum chain strength
    pub min_strength: f64,
    /// Maximum chain strength
    pub max_strength: f64,
    /// Number of strength values to try
    pub num_tries: usize,
}

impl Default for ChainStrengthOptimizer {
    fn default() -> Self {
        Self {
            min_strength: 0.1,
            max_strength: 10.0,
            num_tries: 10,
        }
    }
}

impl ChainStrengthOptimizer {
    /// Find optimal chain strength by analyzing the problem
    #[must_use]
    pub fn find_optimal_strength(&self, logical_problem: &LogicalProblem) -> f64 {
        // Calculate statistics of the logical problem coefficients
        let mut all_coeffs = Vec::new();

        // Add linear coefficients
        for &h in &logical_problem.linear {
            if h.abs() > 1e-10 {
                all_coeffs.push(h.abs());
            }
        }

        // Add quadratic coefficients
        for &J in logical_problem.quadratic.values() {
            if J.abs() > 1e-10 {
                all_coeffs.push(J.abs());
            }
        }

        if all_coeffs.is_empty() {
            return 1.0; // Default strength
        }

        // Sort coefficients
        all_coeffs.sort_by(|a, b| a.total_cmp(b));

        // Use median as base strength
        let median = if all_coeffs.len() % 2 == 0 {
            f64::midpoint(
                all_coeffs[all_coeffs.len() / 2 - 1],
                all_coeffs[all_coeffs.len() / 2],
            )
        } else {
            all_coeffs[all_coeffs.len() / 2]
        };

        // Chain strength should be strong enough to keep chains together
        // but not so strong as to dominate the problem
        (median * 1.5).max(self.min_strength).min(self.max_strength)
    }

    /// Optimize chain strength through multiple runs
    #[must_use]
    pub fn optimize_strength(
        &self,
        logical_problem: &LogicalProblem,
        test_solutions: &[Vec<i8>],
    ) -> f64 {
        let mut best_strength = self.find_optimal_strength(logical_problem);
        let mut best_score = f64::INFINITY;

        // Try different strengths
        let step = (self.max_strength - self.min_strength) / (self.num_tries as f64);

        for i in 0..self.num_tries {
            let strength = (i as f64).mul_add(step, self.min_strength);

            // Evaluate this strength
            let score = self.evaluate_strength(strength, logical_problem, test_solutions);

            if score < best_score {
                best_score = score;
                best_strength = strength;
            }
        }

        best_strength
    }

    /// Evaluate a chain strength against a set of representative solutions.
    ///
    /// A chain holds together only if its inter-qubit coupling (the chain
    /// `strength`) is large enough to resist the local field that tries to flip
    /// part of the chain. For each test solution and each logical variable we
    /// compute that local field magnitude
    /// `|h_i + Σ_j J_ij·s_j|` (the energy gradient acting on variable `i`); the
    /// strength must dominate the largest such field to keep every chain intact.
    /// The returned score (lower is better, so it composes with
    /// [`Self::optimize_strength`]'s minimization) penalizes:
    ///
    /// * **under-strength** — `strength` below the required field, scaled by a
    ///   large factor because a broken chain corrupts the embedded solution; and
    /// * **over-strength** — `strength` far above what is needed, which flattens
    ///   the logical problem and degrades solution quality.
    ///
    /// When no test solutions are supplied the worst-case field is estimated
    /// from the problem coefficients alone (every neighbor aligned adversarially).
    fn evaluate_strength(
        &self,
        strength: f64,
        logical_problem: &LogicalProblem,
        test_solutions: &[Vec<i8>],
    ) -> f64 {
        let num_vars = logical_problem.linear.len();

        // Largest local field observed across variables and test solutions.
        let mut required_strength = 0.0_f64;

        if test_solutions.is_empty() {
            // Worst case: every coupling pulls in the breaking direction.
            for i in 0..num_vars {
                let mut field = logical_problem.linear[i].abs();
                for (&(a, b), &j) in &logical_problem.quadratic {
                    if a == i || b == i {
                        field += j.abs();
                    }
                }
                required_strength = required_strength.max(field);
            }
        } else {
            for solution in test_solutions {
                for i in 0..num_vars.min(solution.len()) {
                    let mut field = logical_problem.linear[i];
                    for (&(a, b), &j) in &logical_problem.quadratic {
                        if a == i {
                            if let Some(&s) = solution.get(b) {
                                field += j * f64::from(s);
                            }
                        } else if b == i {
                            if let Some(&s) = solution.get(a) {
                                field += j * f64::from(s);
                            }
                        }
                    }
                    required_strength = required_strength.max(field.abs());
                }
            }
        }

        if required_strength <= 0.0 {
            // No field to resist: any positive strength is fine; prefer the
            // smallest to avoid distorting the (trivial) problem.
            return strength;
        }

        // Penalty for being too weak (chains break) — heavily weighted.
        const BREAK_PENALTY: f64 = 10.0;
        let under = (required_strength - strength).max(0.0) / required_strength;

        // Penalty for being unnecessarily strong (problem distortion).
        let over = (strength - required_strength).max(0.0) / required_strength;

        BREAK_PENALTY.mul_add(under, over)
    }
}

/// Statistics about chain breaks
#[derive(Debug, Clone, Default)]
pub struct ChainBreakStats {
    /// Total number of chains
    pub total_chains: usize,
    /// Number of broken chains per solution
    pub broken_chains: Vec<usize>,
    /// Chain break rate
    pub break_rate: f64,
    /// Most frequently broken variables
    pub frequent_breaks: Vec<(usize, usize)>,
}

impl ChainBreakStats {
    /// Analyze chain breaks across multiple solutions
    pub fn analyze(
        hardware_solutions: &[HardwareSolution],
        embedding: &Embedding,
    ) -> IsingResult<Self> {
        let total_chains = embedding.chains.len();
        let mut broken_chains = Vec::new();
        let mut break_counts: HashMap<usize, usize> = HashMap::new();

        for hw_solution in hardware_solutions {
            let mut breaks_in_solution = 0;

            for (&var, chain) in &embedding.chains {
                if chain.len() > 1 {
                    let first_spin = hw_solution.spins[chain[0]];
                    let is_broken = chain[1..]
                        .iter()
                        .any(|&q| hw_solution.spins[q] != first_spin);

                    if is_broken {
                        breaks_in_solution += 1;
                        *break_counts.entry(var).or_insert(0) += 1;
                    }
                }
            }

            broken_chains.push(breaks_in_solution);
        }

        // Calculate statistics
        let total_breaks: usize = broken_chains.iter().sum();
        let break_rate = if hardware_solutions.is_empty() || total_chains == 0 {
            0.0
        } else {
            total_breaks as f64 / (hardware_solutions.len() * total_chains) as f64
        };

        // Find most frequently broken variables
        let mut frequent_breaks: Vec<(usize, usize)> = break_counts.into_iter().collect();
        frequent_breaks.sort_by_key(|&(_, count)| std::cmp::Reverse(count));
        frequent_breaks.truncate(10); // Keep top 10

        Ok(Self {
            total_chains,
            broken_chains,
            break_rate,
            frequent_breaks,
        })
    }

    /// Get recommendations based on statistics
    #[must_use]
    pub fn get_recommendations(&self) -> Vec<String> {
        let mut recommendations = Vec::new();

        if self.break_rate > 0.5 {
            recommendations.push(
                "High chain break rate detected. Consider increasing chain strength.".to_string(),
            );
        }

        if self.break_rate > 0.2 {
            recommendations.push(
                "Moderate chain breaks. Try optimizing embedding or chain strength.".to_string(),
            );
        }

        if !self.frequent_breaks.is_empty() {
            let vars: Vec<String> = self
                .frequent_breaks
                .iter()
                .take(3)
                .map(|(var, _)| var.to_string())
                .collect();
            recommendations.push(format!(
                "Variables {} frequently have broken chains. Check embedding quality.",
                vars.join(", ")
            ));
        }

        recommendations
    }
}

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

    #[test]
    fn test_majority_vote_resolution() {
        let mut embedding = Embedding::new();
        embedding
            .add_chain(0, vec![0, 1, 2])
            .expect("failed to add chain in test");
        embedding
            .add_chain(1, vec![3, 4, 5])
            .expect("failed to add chain in test");

        let hw_solution = HardwareSolution {
            spins: vec![1, 1, -1, -1, -1, -1], // First chain: 2 vs 1, second: unanimous
            energy: -1.0,
            occurrences: 1,
        };

        let resolver = ChainBreakResolver::default();
        let resolved = resolver
            .resolve_solution(&hw_solution, &embedding, None)
            .expect("failed to resolve solution in test");

        assert_eq!(resolved.logical_spins, vec![1, -1]);
        assert_eq!(resolved.chain_breaks, 1); // First chain is broken
    }

    #[test]
    fn test_chain_strength_optimizer() {
        let mut problem = LogicalProblem::new(3);
        problem.linear = vec![1.0, -0.5, 0.0];
        problem.quadratic.insert((0, 1), -2.0);
        problem.quadratic.insert((1, 2), 1.5);

        let optimizer = ChainStrengthOptimizer::default();
        let strength = optimizer.find_optimal_strength(&problem);

        // Should be around the median of coefficients
        assert!(strength > 0.5 && strength < 5.0);
    }

    #[test]
    fn test_evaluate_strength_uses_required_field() {
        // Single coupling J=2 between var 0 and 1; with the test solution both
        // aligned (+1,+1) the local field on each variable is |J·s| = 2, so the
        // required chain strength is 2.0.
        let mut problem = LogicalProblem::new(2);
        problem.linear = vec![0.0, 0.0];
        problem.quadratic.insert((0, 1), 2.0);

        let optimizer = ChainStrengthOptimizer::default();
        let test_solutions = vec![vec![1_i8, 1]];

        // A strength below the required field is penalized far more heavily than
        // a strength exactly at the required field.
        let score_weak = optimizer.evaluate_strength(0.5, &problem, &test_solutions);
        let score_matched = optimizer.evaluate_strength(2.0, &problem, &test_solutions);
        let score_strong = optimizer.evaluate_strength(6.0, &problem, &test_solutions);

        // The matched strength is the best (lowest score).
        assert!(score_matched < score_weak);
        assert!(score_matched < score_strong);
        // At the required strength both penalties vanish.
        assert!(score_matched.abs() < 1e-12);
        // Under-strength is penalized by the heavy BREAK_PENALTY factor relative
        // to an equal over-strength deviation.
        let score_under = optimizer.evaluate_strength(1.0, &problem, &test_solutions); // 1 below
        let score_over = optimizer.evaluate_strength(3.0, &problem, &test_solutions); // 1 above
        assert!(score_under > score_over);

        // optimize_strength must pick a strength resisting the field (>= ~2).
        let best = optimizer.optimize_strength(&problem, &test_solutions);
        assert!(
            best >= 1.5,
            "expected strength resisting the field, got {best}"
        );
    }

    #[test]
    fn test_chain_break_stats() {
        let mut embedding = Embedding::new();
        embedding
            .add_chain(0, vec![0, 1])
            .expect("failed to add chain in test");
        embedding
            .add_chain(1, vec![2, 3])
            .expect("failed to add chain in test");

        let solutions = vec![
            HardwareSolution {
                spins: vec![1, 1, -1, -1], // No breaks
                energy: -1.0,
                occurrences: 1,
            },
            HardwareSolution {
                spins: vec![1, -1, -1, -1], // First chain broken
                energy: -0.5,
                occurrences: 1,
            },
        ];

        let stats = ChainBreakStats::analyze(&solutions, &embedding)
            .expect("failed to analyze chain break stats in test");

        assert_eq!(stats.total_chains, 2);
        assert_eq!(stats.broken_chains, vec![0, 1]);
        assert_eq!(stats.break_rate, 0.25); // 1 break out of 4 chain instances
    }
}