quantrs2_anneal/applications/

protein_folding.rs

//! Protein Folding Optimization with Quantum Error Correction
//!
//! This module implements advanced protein folding optimization using quantum annealing
//! with integrated quantum error correction and cutting-edge algorithms. It addresses
//! the fundamental protein folding problem using the HP (Hydrophobic-Polar) model and
//! advanced lattice models with quantum-enhanced optimization.
//!
//! Key Features:
//! - HP model and extended lattice formulations
//! - Quantum error correction for noise-resilient folding predictions
//! - Advanced algorithms (∞-QAOA, Zeno annealing, adiabatic shortcuts)
//! - Neural network guided optimization schedules
//! - Bayesian hyperparameter optimization
//! - Multi-objective optimization (energy, compactness, stability)

use std::collections::{HashMap, VecDeque};
use std::fmt;

use crate::advanced_quantum_algorithms::{
    AdvancedAlgorithmConfig, AdvancedQuantumAlgorithms, AlgorithmSelectionStrategy,
    InfiniteDepthQAOA, InfiniteQAOAConfig, QuantumZenoAnnealer, ZenoConfig,
};
use crate::applications::{
    ApplicationError, ApplicationResult, IndustrySolution, OptimizationProblem,
};
use crate::bayesian_hyperopt::{optimize_annealing_parameters, BayesianHyperoptimizer};
use crate::ising::{IsingModel, QuboModel};
use crate::neural_annealing_schedules::{NeuralAnnealingScheduler, NeuralSchedulerConfig};
use crate::quantum_error_correction::{
    ErrorCorrectionCode, ErrorMitigationConfig, ErrorMitigationManager, LogicalAnnealingEncoder,
    NoiseResilientAnnealingProtocol, SyndromeDetector,
};
use crate::simulator::{AnnealingParams, QuantumAnnealingSimulator};
use std::fmt::Write;

/// Amino acid types in the HP model
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub enum AminoAcidType {
    /// Hydrophobic amino acid
    Hydrophobic,
    /// Polar amino acid
    Polar,
}

impl fmt::Display for AminoAcidType {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            Self::Hydrophobic => write!(f, "H"),
            Self::Polar => write!(f, "P"),
        }
    }
}

/// Protein sequence representation
#[derive(Debug, Clone)]
pub struct ProteinSequence {
    /// Sequence of amino acids
    pub sequence: Vec<AminoAcidType>,
    /// Sequence identifier
    pub id: String,
    /// Additional metadata
    pub metadata: HashMap<String, String>,
}

impl ProteinSequence {
    /// Create new protein sequence
    #[must_use]
    pub fn new(sequence: Vec<AminoAcidType>, id: String) -> Self {
        Self {
            sequence,
            id,
            metadata: HashMap::new(),
        }
    }

    /// Create sequence from string representation
    pub fn from_string(sequence: &str, id: String) -> ApplicationResult<Self> {
        let mut amino_acids = Vec::new();

        for ch in sequence.chars() {
            match ch.to_ascii_uppercase() {
                'H' => amino_acids.push(AminoAcidType::Hydrophobic),
                'P' => amino_acids.push(AminoAcidType::Polar),
                _ => {
                    return Err(ApplicationError::DataValidationError(format!(
                        "Invalid amino acid character: {ch}"
                    )))
                }
            }
        }

        Ok(Self::new(amino_acids, id))
    }

    /// Get sequence length
    #[must_use]
    pub fn length(&self) -> usize {
        self.sequence.len()
    }

    /// Count hydrophobic residues
    #[must_use]
    pub fn hydrophobic_count(&self) -> usize {
        self.sequence
            .iter()
            .filter(|&&aa| aa == AminoAcidType::Hydrophobic)
            .count()
    }

    /// Count polar residues
    #[must_use]
    pub fn polar_count(&self) -> usize {
        self.sequence
            .iter()
            .filter(|&&aa| aa == AminoAcidType::Polar)
            .count()
    }
}

/// Lattice types for protein folding
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum LatticeType {
    /// 2D square lattice
    Square2D,
    /// 2D triangular lattice
    Triangular2D,
    /// 3D cubic lattice
    Cubic3D,
    /// 3D face-centered cubic lattice
    FCC3D,
}

/// Position on a lattice
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub struct LatticePosition {
    pub x: i32,
    pub y: i32,
    pub z: i32,
}

impl LatticePosition {
    #[must_use]
    pub const fn new(x: i32, y: i32, z: i32) -> Self {
        Self { x, y, z }
    }

    #[must_use]
    pub fn distance(&self, other: &Self) -> f64 {
        let dx = f64::from(self.x - other.x);
        let dy = f64::from(self.y - other.y);
        let dz = f64::from(self.z - other.z);
        dz.mul_add(dz, dx.mul_add(dx, dy * dy)).sqrt()
    }

    #[must_use]
    pub const fn manhattan_distance(&self, other: &Self) -> i32 {
        (self.x - other.x).abs() + (self.y - other.y).abs() + (self.z - other.z).abs()
    }
}

/// Protein folding configuration on a lattice
#[derive(Debug, Clone)]
pub struct ProteinFolding {
    /// The protein sequence
    pub sequence: ProteinSequence,
    /// Lattice type
    pub lattice_type: LatticeType,
    /// Positions of amino acids on the lattice
    pub positions: Vec<LatticePosition>,
}

impl ProteinFolding {
    /// Create new protein folding
    #[must_use]
    pub fn new(sequence: ProteinSequence, lattice_type: LatticeType) -> Self {
        let positions = vec![LatticePosition::new(0, 0, 0); sequence.length()];
        Self {
            sequence,
            lattice_type,
            positions,
        }
    }

    /// Set position for amino acid at index
    pub fn set_position(
        &mut self,
        index: usize,
        position: LatticePosition,
    ) -> ApplicationResult<()> {
        if index >= self.sequence.length() {
            return Err(ApplicationError::InvalidConfiguration(format!(
                "Index {} out of bounds for sequence length {}",
                index,
                self.sequence.length()
            )));
        }
        self.positions[index] = position;
        Ok(())
    }

    /// Check if folding is valid (no overlaps, connected chain)
    #[must_use]
    pub fn is_valid(&self) -> bool {
        // Check for position overlaps
        let mut position_set = std::collections::HashSet::new();
        for &pos in &self.positions {
            if !position_set.insert(pos) {
                return false; // Overlap found
            }
        }

        // Check chain connectivity
        for i in 1..self.positions.len() {
            let dist = self.positions[i - 1].manhattan_distance(&self.positions[i]);
            if dist != 1 {
                return false; // Not connected
            }
        }

        true
    }

    /// Calculate hydrophobic-hydrophobic contacts
    #[must_use]
    pub fn hydrophobic_contacts(&self) -> i32 {
        let mut contacts = 0;

        for i in 0..self.sequence.length() {
            if self.sequence.sequence[i] != AminoAcidType::Hydrophobic {
                continue;
            }

            for j in (i + 2)..self.sequence.length() {
                // Skip adjacent in sequence
                if self.sequence.sequence[j] != AminoAcidType::Hydrophobic {
                    continue;
                }

                if self.positions[i].manhattan_distance(&self.positions[j]) == 1 {
                    contacts += 1;
                }
            }
        }

        contacts
    }

    /// Calculate compactness (radius of gyration)
    #[must_use]
    pub fn radius_of_gyration(&self) -> f64 {
        let n = self.positions.len() as f64;

        // Calculate center of mass
        let cx = self.positions.iter().map(|p| f64::from(p.x)).sum::<f64>() / n;
        let cy = self.positions.iter().map(|p| f64::from(p.y)).sum::<f64>() / n;
        let cz = self.positions.iter().map(|p| f64::from(p.z)).sum::<f64>() / n;

        // Calculate radius of gyration
        let rg_sq = self
            .positions
            .iter()
            .map(|p| {
                let dx = f64::from(p.x) - cx;
                let dy = f64::from(p.y) - cy;
                let dz = f64::from(p.z) - cz;
                dz.mul_add(dz, dx.mul_add(dx, dy * dy))
            })
            .sum::<f64>()
            / n;

        rg_sq.sqrt()
    }

    /// Calculate total energy (negative hydrophobic contacts + compactness penalty)
    #[must_use]
    pub fn total_energy(&self) -> f64 {
        let hh_contacts = f64::from(self.hydrophobic_contacts());
        let compactness_penalty = self.radius_of_gyration();

        // Energy = negative contacts (favorable) + compactness penalty
        0.1f64.mul_add(compactness_penalty, -hh_contacts)
    }
}

/// Protein folding optimization problem
#[derive(Debug, Clone)]
pub struct ProteinFoldingProblem {
    /// The protein sequence to fold
    pub sequence: ProteinSequence,
    /// Lattice type for folding
    pub lattice_type: LatticeType,
    /// Optimization objectives
    pub objectives: Vec<FoldingObjective>,
    /// Quantum error correction framework
    pub qec_framework: Option<String>,
    /// Advanced algorithm configuration
    pub advanced_config: AdvancedAlgorithmConfig,
    /// Neural scheduling configuration
    pub neural_config: Option<NeuralSchedulerConfig>,
}

/// Folding optimization objectives
#[derive(Debug, Clone)]
pub enum FoldingObjective {
    /// Maximize hydrophobic-hydrophobic contacts
    MaximizeHHContacts,
    /// Minimize radius of gyration (compactness)
    MinimizeRadiusOfGyration,
    /// Minimize total energy
    MinimizeTotalEnergy,
    /// Maximize structural stability
    MaximizeStability,
}

impl ProteinFoldingProblem {
    /// Create new protein folding problem
    #[must_use]
    pub fn new(sequence: ProteinSequence, lattice_type: LatticeType) -> Self {
        Self {
            sequence,
            lattice_type,
            objectives: vec![FoldingObjective::MaximizeHHContacts],
            qec_framework: None,
            advanced_config: AdvancedAlgorithmConfig {
                enable_infinite_qaoa: true,
                enable_zeno_annealing: true,
                enable_adiabatic_shortcuts: true,
                enable_counterdiabatic: true,
                selection_strategy: AlgorithmSelectionStrategy::ProblemSpecific,
                track_performance: true,
            },
            neural_config: None,
        }
    }

    /// Enable quantum error correction
    #[must_use]
    pub fn with_quantum_error_correction(mut self, config: String) -> Self {
        self.qec_framework = Some(config);
        self
    }

    /// Enable neural annealing schedules
    #[must_use]
    pub fn with_neural_annealing(mut self, config: NeuralSchedulerConfig) -> Self {
        self.neural_config = Some(config);
        self
    }

    /// Add optimization objective
    #[must_use]
    pub fn add_objective(mut self, objective: FoldingObjective) -> Self {
        self.objectives.push(objective);
        self
    }

    /// Solve using advanced quantum algorithms
    pub fn solve_with_advanced_algorithms(&self) -> ApplicationResult<ProteinFolding> {
        println!("Starting protein folding optimization with advanced quantum algorithms");

        // Convert to QUBO formulation
        let (qubo, variable_map) = self.to_qubo()?;

        // Create advanced quantum algorithms coordinator
        let algorithms = AdvancedQuantumAlgorithms::with_config(self.advanced_config.clone());

        // Solve using best algorithm for this problem size
        let result = algorithms.solve(&qubo, None).map_err(|e| {
            ApplicationError::OptimizationError(format!("Advanced algorithm failed: {e:?}"))
        })?;

        let solution = result
            .map_err(|e| ApplicationError::OptimizationError(format!("Solver error: {e}")))?;

        // Convert binary solution back to protein folding
        self.solution_from_binary(&solution, &variable_map)
    }

    /// Solve using quantum error correction
    pub fn solve_with_qec(&self) -> ApplicationResult<ProteinFolding> {
        if let Some(ref qec_framework) = self.qec_framework {
            println!("Starting noise-resilient protein folding optimization");

            // Convert to QUBO and then to Ising
            let (qubo, variable_map) = self.to_qubo()?;
            let ising_model = qubo.to_ising();

            // Use error mitigation for protein folding optimization
            let error_config = ErrorMitigationConfig::default();
            let mut error_manager = ErrorMitigationManager::new(error_config).map_err(|e| {
                ApplicationError::OptimizationError(format!(
                    "Failed to create error manager: {e:?}"
                ))
            })?;

            // First perform standard annealing
            let params = AnnealingParams::default();
            let annealer = QuantumAnnealingSimulator::new(params.clone()).map_err(|e| {
                ApplicationError::OptimizationError(format!("Failed to create annealer: {e:?}"))
            })?;
            let annealing_result = annealer.solve(&ising_model.0).map_err(|e| {
                ApplicationError::OptimizationError(format!("Annealing failed: {e:?}"))
            })?;

            // Convert simulator result to error mitigation format
            let error_mitigation_result =
                crate::quantum_error_correction::error_mitigation::AnnealingResult {
                    solution: annealing_result
                        .best_spins
                        .iter()
                        .map(|&x| i32::from(x))
                        .collect(),
                    energy: annealing_result.best_energy,
                    num_occurrences: 1,
                    chain_break_fraction: 0.0,
                    timing: std::collections::HashMap::new(),
                    info: std::collections::HashMap::new(),
                };

            // Apply error mitigation to improve the result
            let mitigation_result = error_manager
                .apply_mitigation(&ising_model.0, error_mitigation_result, &params)
                .map_err(|e| {
                    ApplicationError::OptimizationError(format!("Error mitigation failed: {e:?}"))
                })?;

            let solution = &mitigation_result.mitigated_result.solution;

            // Convert back to folding
            self.solution_from_binary(&solution, &variable_map)
        } else {
            Err(ApplicationError::InvalidConfiguration(
                "Quantum error correction not enabled".to_string(),
            ))
        }
    }

    /// Solve using neural annealing schedules
    pub fn solve_with_neural_annealing(&self) -> ApplicationResult<ProteinFolding> {
        if let Some(ref neural_config) = self.neural_config {
            println!("Starting protein folding with neural-guided annealing");

            let (qubo, variable_map) = self.to_qubo()?;

            // Create neural annealing scheduler
            let mut neural_scheduler = NeuralAnnealingScheduler::new(neural_config.clone())
                .map_err(|e| {
                    ApplicationError::OptimizationError(format!(
                        "Failed to create neural scheduler: {e:?}"
                    ))
                })?;

            // Optimize using neural-guided schedules
            let result = neural_scheduler.optimize(&qubo).map_err(|e| {
                ApplicationError::OptimizationError(format!("Neural annealing failed: {e:?}"))
            })?;

            let solution = result.map_err(|e| {
                ApplicationError::OptimizationError(format!("Neural solver error: {e}"))
            })?;

            self.solution_from_binary(&solution, &variable_map)
        } else {
            Err(ApplicationError::InvalidConfiguration(
                "Neural annealing not enabled".to_string(),
            ))
        }
    }

    /// Optimize hyperparameters using Bayesian optimization
    pub fn optimize_hyperparameters(&self) -> ApplicationResult<HashMap<String, f64>> {
        println!("Optimizing folding hyperparameters with Bayesian optimization");

        // Define objective function for hyperparameter optimization
        let objective = |params: &[f64]| -> f64 {
            // params[0] = temperature, params[1] = annealing steps, params[2] = algorithm type
            let temperature = params[0];
            let steps = params[1] as usize;
            // let _algorithm_type = params[2] as usize;

            // Simple folding energy approximation for optimization
            let length = self.sequence.length() as f64;
            let hydrophobic_ratio = self.sequence.hydrophobic_count() as f64 / length;

            // Simulate annealing performance
            let energy_scale = -hydrophobic_ratio * length / 4.0; // Approximate optimal HH contacts
            let thermal_noise = temperature * 0.1;
            let convergence_factor = (steps as f64 / 1000.0).min(1.0);

            // Return negative energy (minimize)
            -(energy_scale * convergence_factor - thermal_noise)
        };

        // Run Bayesian optimization
        let best_params = optimize_annealing_parameters(objective, Some(30)).map_err(|e| {
            ApplicationError::OptimizationError(format!("Bayesian optimization failed: {e:?}"))
        })?;

        let mut result = HashMap::new();
        result.insert("optimal_temperature".to_string(), best_params[0]);
        result.insert("optimal_steps".to_string(), best_params[1]);
        result.insert("optimal_algorithm".to_string(), best_params[2]);

        Ok(result)
    }

    /// Convert binary solution vector to protein folding
    fn solution_from_binary(
        &self,
        solution: &[i32],
        variable_map: &HashMap<String, usize>,
    ) -> ApplicationResult<ProteinFolding> {
        let mut folding = ProteinFolding::new(self.sequence.clone(), self.lattice_type);

        // Decode positions from binary solution
        // This is a simplified decoding - in practice would use more sophisticated encoding
        let seq_len = self.sequence.length();

        // Place first amino acid at origin
        folding.set_position(0, LatticePosition::new(0, 0, 0))?;

        // Decode subsequent positions using relative moves
        for i in 1..seq_len {
            let move_var_base = (i - 1) * 2; // 2 bits per move for 2D

            if move_var_base + 1 < solution.len() {
                let move_x = i32::from(solution[move_var_base] > 0);
                let move_y = i32::from(solution[move_var_base + 1] > 0);

                let prev_pos = folding.positions[i - 1];
                let new_pos = match (move_x, move_y) {
                    (0, 0) => LatticePosition::new(prev_pos.x + 1, prev_pos.y, prev_pos.z),
                    (0, 1) => LatticePosition::new(prev_pos.x, prev_pos.y + 1, prev_pos.z),
                    (1, 0) => LatticePosition::new(prev_pos.x - 1, prev_pos.y, prev_pos.z),
                    (1, 1) => LatticePosition::new(prev_pos.x, prev_pos.y - 1, prev_pos.z),
                    _ => LatticePosition::new(prev_pos.x + 1, prev_pos.y, prev_pos.z), // Default case
                };

                folding.set_position(i, new_pos)?;
            }
        }

        // Verify and fix if necessary
        if !folding.is_valid() {
            // Apply simple validation fix
            self.fix_invalid_folding(&mut folding)?;
        }

        Ok(folding)
    }

    /// Fix invalid folding configuration
    fn fix_invalid_folding(&self, folding: &mut ProteinFolding) -> ApplicationResult<()> {
        // Simple SAW (Self-Avoiding Walk) generation
        let mut positions = vec![LatticePosition::new(0, 0, 0)];
        let mut used_positions = std::collections::HashSet::new();
        used_positions.insert(positions[0]);

        // Generate valid chain
        for i in 1..self.sequence.length() {
            let prev_pos = positions[i - 1];

            // Try each direction
            let directions = [
                LatticePosition::new(1, 0, 0),
                LatticePosition::new(-1, 0, 0),
                LatticePosition::new(0, 1, 0),
                LatticePosition::new(0, -1, 0),
            ];

            let mut placed = false;
            for &dir in &directions {
                let new_pos = LatticePosition::new(
                    prev_pos.x + dir.x,
                    prev_pos.y + dir.y,
                    prev_pos.z + dir.z,
                );

                if used_positions.insert(new_pos) {
                    positions.push(new_pos);
                    placed = true;
                    break;
                }
            }

            if !placed {
                return Err(ApplicationError::OptimizationError(
                    "Cannot generate valid folding configuration".to_string(),
                ));
            }
        }

        folding.positions = positions;
        Ok(())
    }
}

impl OptimizationProblem for ProteinFoldingProblem {
    type Solution = ProteinFolding;
    type ObjectiveValue = f64;

    fn description(&self) -> String {
        format!(
            "Protein folding optimization for sequence {} (length: {}) on {:?} lattice",
            self.sequence.id,
            self.sequence.length(),
            self.lattice_type
        )
    }

    fn size_metrics(&self) -> HashMap<String, usize> {
        let mut metrics = HashMap::new();
        metrics.insert("sequence_length".to_string(), self.sequence.length());
        metrics.insert(
            "hydrophobic_count".to_string(),
            self.sequence.hydrophobic_count(),
        );
        metrics.insert("polar_count".to_string(), self.sequence.polar_count());
        metrics.insert("variables".to_string(), (self.sequence.length() - 1) * 2); // Simplified
        metrics
    }

    fn validate(&self) -> ApplicationResult<()> {
        if self.sequence.length() < 2 {
            return Err(ApplicationError::DataValidationError(
                "Sequence must have at least 2 amino acids".to_string(),
            ));
        }

        if self.sequence.length() > 200 {
            return Err(ApplicationError::ResourceLimitExceeded(
                "Sequence too long for current implementation".to_string(),
            ));
        }

        Ok(())
    }

    fn to_qubo(&self) -> ApplicationResult<(QuboModel, HashMap<String, usize>)> {
        self.validate()?;

        let seq_len = self.sequence.length();
        let num_vars = (seq_len - 1) * 2; // 2 bits per move for 2D lattice

        let mut qubo = QuboModel::new(num_vars);
        let mut variable_map = HashMap::new();

        // Map variables
        for i in 0..(seq_len - 1) {
            variable_map.insert(format!("move_{i}_x"), i * 2);
            variable_map.insert(format!("move_{i}_y"), i * 2 + 1);
        }

        // Add objective terms (simplified QUBO formulation)
        for obj in &self.objectives {
            match obj {
                FoldingObjective::MaximizeHHContacts => {
                    self.add_contact_terms(&mut qubo)?;
                }
                FoldingObjective::MinimizeRadiusOfGyration => {
                    self.add_compactness_terms(&mut qubo)?;
                }
                _ => {
                    // Combined terms for other objectives
                    self.add_contact_terms(&mut qubo)?;
                    self.add_compactness_terms(&mut qubo)?;
                }
            }
        }

        Ok((qubo, variable_map))
    }

    fn evaluate_solution(
        &self,
        solution: &Self::Solution,
    ) -> ApplicationResult<Self::ObjectiveValue> {
        if !solution.is_valid() {
            return Err(ApplicationError::DataValidationError(
                "Invalid folding configuration".to_string(),
            ));
        }

        Ok(solution.total_energy())
    }

    fn is_feasible(&self, solution: &Self::Solution) -> bool {
        solution.is_valid()
    }
}

impl ProteinFoldingProblem {
    /// Add hydrophobic contact terms to QUBO
    fn add_contact_terms(&self, qubo: &mut QuboModel) -> ApplicationResult<()> {
        // Simplified contact terms - in practice would be more complex
        let weight = -1.0; // Negative to encourage contacts

        for i in 0..qubo.num_variables {
            for j in (i + 1)..qubo.num_variables {
                // Add interaction terms that promote favorable configurations
                qubo.set_quadratic(i, j, weight * 0.1)?;
            }
        }

        Ok(())
    }

    /// Add compactness terms to QUBO
    fn add_compactness_terms(&self, qubo: &mut QuboModel) -> ApplicationResult<()> {
        // Penalty for extended configurations
        let penalty = 0.5;

        for i in 0..qubo.num_variables {
            qubo.set_linear(i, penalty)?;
        }

        Ok(())
    }
}

impl IndustrySolution for ProteinFolding {
    type Problem = ProteinFoldingProblem;

    fn from_binary(problem: &Self::Problem, binary_solution: &[i8]) -> ApplicationResult<Self> {
        let solution_i32: Vec<i32> = binary_solution.iter().map(|&x| i32::from(x)).collect();
        let variable_map = HashMap::new(); // Simplified
        problem.solution_from_binary(&solution_i32, &variable_map)
    }

    fn summary(&self) -> HashMap<String, String> {
        let mut summary = HashMap::new();
        summary.insert("sequence_id".to_string(), self.sequence.id.clone());
        summary.insert(
            "sequence_length".to_string(),
            self.sequence.length().to_string(),
        );
        summary.insert(
            "lattice_type".to_string(),
            format!("{:?}", self.lattice_type),
        );
        summary.insert(
            "hydrophobic_contacts".to_string(),
            self.hydrophobic_contacts().to_string(),
        );
        summary.insert(
            "radius_of_gyration".to_string(),
            format!("{:.3}", self.radius_of_gyration()),
        );
        summary.insert(
            "total_energy".to_string(),
            format!("{:.3}", self.total_energy()),
        );
        summary.insert("is_valid".to_string(), self.is_valid().to_string());
        summary
    }

    fn metrics(&self) -> HashMap<String, f64> {
        let mut metrics = HashMap::new();
        metrics.insert(
            "hydrophobic_contacts".to_string(),
            f64::from(self.hydrophobic_contacts()),
        );
        metrics.insert("radius_of_gyration".to_string(), self.radius_of_gyration());
        metrics.insert("total_energy".to_string(), self.total_energy());
        metrics.insert(
            "compactness_score".to_string(),
            1.0 / (1.0 + self.radius_of_gyration()),
        );
        metrics
    }

    fn export_format(&self) -> ApplicationResult<String> {
        let mut output = String::new();

        let _ = writeln!(output, "# Protein Folding Result");
        let _ = writeln!(output, "Sequence ID: {}", self.sequence.id);
        let _ = writeln!(
            output,
            "Sequence: {}",
            self.sequence
                .sequence
                .iter()
                .map(|aa| format!("{aa}"))
                .collect::<String>()
        );
        let _ = writeln!(output, "Lattice Type: {:?}", self.lattice_type);
        let _ = write!(output, "Length: {}\n", self.sequence.length());
        let _ = write!(
            output,
            "Hydrophobic Contacts: {}\n",
            self.hydrophobic_contacts()
        );
        let _ = write!(
            output,
            "Radius of Gyration: {:.3}\n",
            self.radius_of_gyration()
        );
        let _ = write!(output, "Total Energy: {:.3}\n", self.total_energy());
        let _ = write!(output, "Valid Configuration: {}\n", self.is_valid());

        output.push_str("\n# Positions\n");
        for (i, pos) in self.positions.iter().enumerate() {
            let _ = write!(
                output,
                "{}: {} ({}, {}, {})\n",
                i, self.sequence.sequence[i], pos.x, pos.y, pos.z
            );
        }

        Ok(output)
    }
}

/// Create benchmark protein folding problems
pub fn create_benchmark_problems(
    size: usize,
) -> ApplicationResult<
    Vec<Box<dyn OptimizationProblem<Solution = ProteinFolding, ObjectiveValue = f64>>>,
> {
    let mut problems = Vec::new();

    // Generate test sequences of different complexity
    let sequences = match size {
        s if s <= 10 => {
            vec!["HPHPPHHPHH", "HPHHHPPHPH", "HPPHHPPHPP"]
        }
        s if s <= 25 => {
            vec![
                "HPHPPHHPHHHPPHPPHHHPHPH",
                "HHHPPHPPHHPPHPPHHHPHPH",
                "HPHPPHPPHHHPPHPPHHHPHP",
                "HPPHHPPHPPHHHPPHPPHHPH",
            ]
        }
        _ => {
            vec![
                "HPHPPHHPHHHPPHPPHHHPHPHHPHPPHHPHHHPPHPPHHHPHPH",
                "HHHPPHPPHHPPHPPHHHPHPHHHHPPHPPHHPPHPPHHHPHPH",
                "HPHPPHPPHHHPPHPPHHHPHPHPHPPHPPHHHPPHPPHHHPHP",
            ]
        }
    };

    for (i, seq_str) in sequences.iter().enumerate() {
        let sequence = ProteinSequence::from_string(seq_str, format!("benchmark_{i}"))?;
        let problem = ProteinFoldingProblem::new(sequence, LatticeType::Square2D);

        // Note: This is a simplified implementation for the trait object
        // In practice, would need proper trait object conversion
        problems.push(Box::new(problem)
            as Box<
                dyn OptimizationProblem<Solution = ProteinFolding, ObjectiveValue = f64>,
            >);
    }

    Ok(problems)
}

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

    #[test]
    fn test_protein_sequence_creation() {
        let sequence = ProteinSequence::from_string("HPHPPHHPHH", "test".to_string())
            .expect("Failed to create protein sequence from valid string");
        assert_eq!(sequence.length(), 10);
        assert_eq!(sequence.hydrophobic_count(), 6);
        assert_eq!(sequence.polar_count(), 4);
    }

    #[test]
    fn test_lattice_position() {
        let pos1 = LatticePosition::new(0, 0, 0);
        let pos2 = LatticePosition::new(1, 0, 0);

        assert_eq!(pos1.manhattan_distance(&pos2), 1);
        assert!((pos1.distance(&pos2) - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_protein_folding_validation() {
        let sequence = ProteinSequence::from_string("HPHH", "test".to_string())
            .expect("Failed to create protein sequence");
        let mut folding = ProteinFolding::new(sequence, LatticeType::Square2D);

        // Set valid positions
        folding
            .set_position(0, LatticePosition::new(0, 0, 0))
            .expect("Failed to set position 0");
        folding
            .set_position(1, LatticePosition::new(1, 0, 0))
            .expect("Failed to set position 1");
        folding
            .set_position(2, LatticePosition::new(2, 0, 0))
            .expect("Failed to set position 2");
        folding
            .set_position(3, LatticePosition::new(2, 1, 0))
            .expect("Failed to set position 3");

        assert!(folding.is_valid());
    }

    #[test]
    fn test_hydrophobic_contacts() {
        let sequence = ProteinSequence::from_string("HPHH", "test".to_string())
            .expect("Failed to create protein sequence");
        let mut folding = ProteinFolding::new(sequence, LatticeType::Square2D);

        // Create L-shape with HH contact
        folding
            .set_position(0, LatticePosition::new(0, 0, 0))
            .expect("Failed to set position 0"); // H
        folding
            .set_position(1, LatticePosition::new(1, 0, 0))
            .expect("Failed to set position 1"); // P
        folding
            .set_position(2, LatticePosition::new(2, 0, 0))
            .expect("Failed to set position 2"); // H
        folding
            .set_position(3, LatticePosition::new(2, 1, 0))
            .expect("Failed to set position 3"); // H

        // No HH contacts in this configuration
        assert_eq!(folding.hydrophobic_contacts(), 0);

        // Modify to create HH contact
        folding
            .set_position(3, LatticePosition::new(1, 1, 0))
            .expect("Failed to update position 3"); // H adjacent to P
                                                    // Still no HH contact between non-adjacent in sequence

        // Create folded structure with HH contact
        folding
            .set_position(0, LatticePosition::new(0, 0, 0))
            .expect("Failed to set folded position 0"); // H
        folding
            .set_position(1, LatticePosition::new(1, 0, 0))
            .expect("Failed to set folded position 1"); // P
        folding
            .set_position(2, LatticePosition::new(1, 1, 0))
            .expect("Failed to set folded position 2"); // H
        folding
            .set_position(3, LatticePosition::new(0, 1, 0))
            .expect("Failed to set folded position 3"); // H

        // Now positions 0 and 3 (both H) are adjacent
        assert_eq!(folding.hydrophobic_contacts(), 1);
    }

    #[test]
    fn test_problem_creation() {
        let sequence = ProteinSequence::from_string("HPHPPHHPHH", "test".to_string())
            .expect("Failed to create protein sequence");
        let problem = ProteinFoldingProblem::new(sequence, LatticeType::Square2D);

        assert!(problem.validate().is_ok());

        let metrics = problem.size_metrics();
        assert_eq!(metrics["sequence_length"], 10);
    }

    #[test]
    fn test_qubo_conversion() {
        let sequence = ProteinSequence::from_string("HPHH", "test".to_string())
            .expect("Failed to create protein sequence");
        let problem = ProteinFoldingProblem::new(sequence, LatticeType::Square2D);

        let (qubo, variable_map) = problem
            .to_qubo()
            .expect("Failed to convert problem to QUBO");
        assert_eq!(qubo.num_variables, 6); // (4-1) * 2 = 6 variables
        assert!(!variable_map.is_empty());
    }
}