use crate::builder::Circuit;
use crate::scirs2_matrices::SparseMatrix;
use quantrs2_core::{
error::{QuantRS2Error, QuantRS2Result},
gate::GateOp,
qubit::QubitId,
};
use scirs2_core::Complex64;
use serde::{Deserialize, Serialize};
use std::collections::HashMap;
use std::sync::{Arc, Mutex};
pub trait ObjectiveFunction: Send + Sync {
fn evaluate(&self, parameters: &[f64]) -> f64;
fn gradient(&self, parameters: &[f64]) -> Option<Vec<f64>> {
None
}
fn hessian(&self, parameters: &[f64]) -> Option<Vec<Vec<f64>>> {
None
}
fn bounds(&self) -> Vec<(f64, f64)>;
fn name(&self) -> &str;
}
#[derive(Debug, Clone, PartialEq)]
pub enum OptimizationAlgorithm {
GradientDescent { learning_rate: f64, momentum: f64 },
Adam {
learning_rate: f64,
beta1: f64,
beta2: f64,
epsilon: f64,
},
LBFGSB {
max_iterations: usize,
tolerance: f64,
},
NelderMead {
max_iterations: usize,
tolerance: f64,
},
SimulatedAnnealing {
initial_temperature: f64,
cooling_rate: f64,
min_temperature: f64,
},
GeneticAlgorithm {
population_size: usize,
mutation_rate: f64,
crossover_rate: f64,
},
ParticleSwarm {
num_particles: usize,
inertia_weight: f64,
cognitive_weight: f64,
social_weight: f64,
},
BayesianOptimization {
acquisition_function: AcquisitionFunction,
kernel: KernelType,
num_initial_samples: usize,
},
QAOA {
num_layers: usize,
classical_optimizer: Box<Self>,
},
}
#[derive(Debug, Clone, PartialEq)]
pub enum AcquisitionFunction {
ExpectedImprovement,
ProbabilityOfImprovement,
UpperConfidenceBound { kappa: f64 },
Thompson,
}
#[derive(Debug, Clone, PartialEq)]
pub enum KernelType {
RBF { length_scale: f64 },
Matern { nu: f64, length_scale: f64 },
Linear { variance: f64 },
Periodic { period: f64, length_scale: f64 },
}
pub struct OptimizationConfig {
pub algorithm: OptimizationAlgorithm,
pub max_evaluations: usize,
pub tolerance: f64,
pub seed: Option<u64>,
pub parallel: bool,
pub num_threads: Option<usize>,
pub progress_callback: Option<Box<dyn Fn(usize, f64) + Send + Sync>>,
pub early_stopping: Option<EarlyStoppingCriteria>,
}
impl std::fmt::Debug for OptimizationConfig {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
f.debug_struct("OptimizationConfig")
.field("algorithm", &self.algorithm)
.field("max_evaluations", &self.max_evaluations)
.field("tolerance", &self.tolerance)
.field("seed", &self.seed)
.field("parallel", &self.parallel)
.field("num_threads", &self.num_threads)
.field(
"progress_callback",
&self.progress_callback.as_ref().map(|_| "Some(callback)"),
)
.field("early_stopping", &self.early_stopping)
.finish()
}
}
impl Clone for OptimizationConfig {
fn clone(&self) -> Self {
Self {
algorithm: self.algorithm.clone(),
max_evaluations: self.max_evaluations,
tolerance: self.tolerance,
seed: self.seed,
parallel: self.parallel,
num_threads: self.num_threads,
progress_callback: None, early_stopping: self.early_stopping.clone(),
}
}
}
#[derive(Debug, Clone)]
pub struct EarlyStoppingCriteria {
pub patience: usize,
pub min_delta: f64,
pub monitor_best: bool,
}
#[derive(Debug, Clone)]
pub struct OptimizationResult {
pub optimal_parameters: Vec<f64>,
pub optimal_value: f64,
pub num_evaluations: usize,
pub converged: bool,
pub history: OptimizationHistory,
pub algorithm_info: HashMap<String, String>,
pub optimization_time: std::time::Duration,
}
#[derive(Debug, Clone)]
pub struct OptimizationHistory {
pub parameters: Vec<Vec<f64>>,
pub objective_values: Vec<f64>,
pub gradient_norms: Vec<f64>,
pub step_sizes: Vec<f64>,
pub timestamps: Vec<std::time::Instant>,
}
pub struct QuantumCircuitOptimizer {
circuit_template: CircuitTemplate,
config: OptimizationConfig,
history: Arc<Mutex<OptimizationHistory>>,
best_parameters: Arc<Mutex<Option<Vec<f64>>>>,
best_value: Arc<Mutex<f64>>,
}
#[derive(Debug, Clone)]
pub struct CircuitTemplate {
pub structure: Vec<ParameterizedGate>,
pub parameters: Vec<Parameter>,
pub num_qubits: usize,
}
#[derive(Debug, Clone)]
pub struct ParameterizedGate {
pub gate_name: String,
pub qubits: Vec<usize>,
pub parameter_indices: Vec<usize>,
pub fixed_parameters: Vec<f64>,
}
#[derive(Debug, Clone)]
pub struct Parameter {
pub name: String,
pub lower_bound: f64,
pub upper_bound: f64,
pub initial_value: f64,
pub discrete: bool,
}
impl QuantumCircuitOptimizer {
#[must_use]
pub fn new(template: CircuitTemplate, config: OptimizationConfig) -> Self {
Self {
circuit_template: template,
config,
history: Arc::new(Mutex::new(OptimizationHistory {
parameters: Vec::new(),
objective_values: Vec::new(),
gradient_norms: Vec::new(),
step_sizes: Vec::new(),
timestamps: Vec::new(),
})),
best_parameters: Arc::new(Mutex::new(None)),
best_value: Arc::new(Mutex::new(f64::INFINITY)),
}
}
pub fn optimize(
&mut self,
objective: Arc<dyn ObjectiveFunction>,
) -> QuantRS2Result<OptimizationResult> {
let start_time = std::time::Instant::now();
let initial_params: Vec<f64> = self
.circuit_template
.parameters
.iter()
.map(|p| p.initial_value)
.collect();
let bounds = objective.bounds();
if bounds.len() != initial_params.len() {
return Err(QuantRS2Error::InvalidInput(
"Parameter count mismatch with bounds".to_string(),
));
}
let result = match &self.config.algorithm {
OptimizationAlgorithm::GradientDescent {
learning_rate,
momentum,
} => self.optimize_gradient_descent(
objective,
&initial_params,
*learning_rate,
*momentum,
),
OptimizationAlgorithm::Adam {
learning_rate,
beta1,
beta2,
epsilon,
} => self.optimize_adam(
objective,
&initial_params,
*learning_rate,
*beta1,
*beta2,
*epsilon,
),
OptimizationAlgorithm::LBFGSB {
max_iterations,
tolerance,
} => self.optimize_lbfgs(objective, &initial_params, *max_iterations, *tolerance),
OptimizationAlgorithm::NelderMead {
max_iterations,
tolerance,
} => self.optimize_nelder_mead(objective, &initial_params, *max_iterations, *tolerance),
OptimizationAlgorithm::SimulatedAnnealing {
initial_temperature,
cooling_rate,
min_temperature,
} => self.optimize_simulated_annealing(
objective,
&initial_params,
*initial_temperature,
*cooling_rate,
*min_temperature,
),
OptimizationAlgorithm::BayesianOptimization {
acquisition_function,
kernel,
num_initial_samples,
} => self.optimize_bayesian(
objective,
&initial_params,
acquisition_function,
kernel,
*num_initial_samples,
),
_ => Err(QuantRS2Error::InvalidInput(
"Algorithm not yet implemented".to_string(),
)),
}?;
let history = self
.history
.lock()
.map_err(|e| QuantRS2Error::RuntimeError(format!("Failed to lock history: {}", e)))?
.clone();
Ok(OptimizationResult {
optimal_parameters: result.0,
optimal_value: result.1,
num_evaluations: result.2,
converged: result.3,
history,
algorithm_info: HashMap::new(),
optimization_time: start_time.elapsed(),
})
}
fn optimize_gradient_descent(
&self,
objective: Arc<dyn ObjectiveFunction>,
initial_params: &[f64],
learning_rate: f64,
momentum: f64,
) -> QuantRS2Result<(Vec<f64>, f64, usize, bool)> {
let mut params = initial_params.to_vec();
let mut velocity = vec![0.0; params.len()];
let mut evaluations = 0;
let mut best_value = f64::INFINITY;
for iteration in 0..self.config.max_evaluations {
let value = objective.evaluate(¶ms);
evaluations += 1;
if value < best_value {
best_value = value;
if let Ok(mut guard) = self.best_parameters.lock() {
*guard = Some(params.clone());
}
if let Ok(mut guard) = self.best_value.lock() {
*guard = best_value;
}
}
self.record_iteration(¶ms, value, iteration);
if iteration > 0 {
let prev_value = self
.history
.lock()
.ok()
.and_then(|h| h.objective_values.get(iteration - 1).copied())
.unwrap_or(value);
if (prev_value - value).abs() < self.config.tolerance {
return Ok((params, best_value, evaluations, true));
}
}
let gradient = if let Some(grad) = objective.gradient(¶ms) {
grad
} else {
self.numerical_gradient(&*objective, ¶ms)?
};
for i in 0..params.len() {
velocity[i] = momentum.mul_add(velocity[i], -(learning_rate * gradient[i]));
params[i] += velocity[i];
let bounds = objective.bounds();
params[i] = params[i].max(bounds[i].0).min(bounds[i].1);
}
if let Some(callback) = &self.config.progress_callback {
callback(iteration, value);
}
}
Ok((params, best_value, evaluations, false))
}
fn optimize_adam(
&self,
objective: Arc<dyn ObjectiveFunction>,
initial_params: &[f64],
learning_rate: f64,
beta1: f64,
beta2: f64,
epsilon: f64,
) -> QuantRS2Result<(Vec<f64>, f64, usize, bool)> {
let mut params = initial_params.to_vec();
let mut m = vec![0.0; params.len()]; let mut v = vec![0.0; params.len()]; let mut evaluations = 0;
let mut best_value = f64::INFINITY;
for iteration in 0..self.config.max_evaluations {
let t = iteration + 1;
let value = objective.evaluate(¶ms);
evaluations += 1;
if value < best_value {
best_value = value;
if let Ok(mut guard) = self.best_parameters.lock() {
*guard = Some(params.clone());
}
if let Ok(mut guard) = self.best_value.lock() {
*guard = best_value;
}
}
self.record_iteration(¶ms, value, iteration);
if iteration > 0 {
let prev_value = self
.history
.lock()
.ok()
.and_then(|h| h.objective_values.get(iteration - 1).copied())
.unwrap_or(value);
if (prev_value - value).abs() < self.config.tolerance {
return Ok((params, best_value, evaluations, true));
}
}
let gradient = if let Some(grad) = objective.gradient(¶ms) {
grad
} else {
self.numerical_gradient(&*objective, ¶ms)?
};
for i in 0..params.len() {
m[i] = beta1.mul_add(m[i], (1.0 - beta1) * gradient[i]);
v[i] = beta2.mul_add(v[i], (1.0 - beta2) * gradient[i] * gradient[i]);
let m_hat = m[i] / (1.0 - beta1.powi(t as i32));
let v_hat = v[i] / (1.0 - beta2.powi(t as i32));
params[i] -= learning_rate * m_hat / (v_hat.sqrt() + epsilon);
let bounds = objective.bounds();
params[i] = params[i].max(bounds[i].0).min(bounds[i].1);
}
if let Some(callback) = &self.config.progress_callback {
callback(iteration, value);
}
}
Ok((params, best_value, evaluations, false))
}
fn optimize_lbfgs(
&self,
objective: Arc<dyn ObjectiveFunction>,
initial_params: &[f64],
max_iterations: usize,
tolerance: f64,
) -> QuantRS2Result<(Vec<f64>, f64, usize, bool)> {
self.optimize_gradient_descent(objective, initial_params, 0.01, 0.9)
}
fn optimize_nelder_mead(
&self,
objective: Arc<dyn ObjectiveFunction>,
initial_params: &[f64],
max_iterations: usize,
tolerance: f64,
) -> QuantRS2Result<(Vec<f64>, f64, usize, bool)> {
let n = initial_params.len();
let mut simplex = Vec::new();
let mut evaluations = 0;
simplex.push(initial_params.to_vec());
for i in 0..n {
let mut vertex = initial_params.to_vec();
vertex[i] += if vertex[i] == 0.0 {
0.00025
} else {
vertex[i] * 0.05
};
simplex.push(vertex);
}
let mut values: Vec<f64> = simplex
.iter()
.map(|params| {
evaluations += 1;
objective.evaluate(params)
})
.collect();
for iteration in 0..max_iterations {
let mut indices: Vec<usize> = (0..simplex.len()).collect();
indices.sort_by(|&i, &j| {
values[i]
.partial_cmp(&values[j])
.unwrap_or(std::cmp::Ordering::Equal)
});
let best_value = values[indices[0]];
let worst_idx = indices[n];
let second_worst_idx = indices[n - 1];
self.record_iteration(&simplex[indices[0]], best_value, iteration);
let range = values[worst_idx] - values[indices[0]];
if range < tolerance {
return Ok((simplex[indices[0]].clone(), best_value, evaluations, true));
}
let mut centroid = vec![0.0; n];
for i in 0..n {
for j in 0..n {
centroid[j] += simplex[indices[i]][j];
}
}
for j in 0..n {
centroid[j] /= n as f64;
}
let alpha = 1.0;
let mut reflected = vec![0.0; n];
for j in 0..n {
reflected[j] = centroid[j] + alpha * (centroid[j] - simplex[worst_idx][j]);
}
let bounds = objective.bounds();
for j in 0..n {
reflected[j] = reflected[j].max(bounds[j].0).min(bounds[j].1);
}
let reflected_value = objective.evaluate(&reflected);
evaluations += 1;
if values[indices[0]] <= reflected_value && reflected_value < values[second_worst_idx] {
simplex[worst_idx] = reflected;
values[worst_idx] = reflected_value;
} else if reflected_value < values[indices[0]] {
let gamma = 2.0;
let mut expanded = vec![0.0; n];
for j in 0..n {
expanded[j] = centroid[j] + gamma * (reflected[j] - centroid[j]);
expanded[j] = expanded[j].max(bounds[j].0).min(bounds[j].1);
}
let expanded_value = objective.evaluate(&expanded);
evaluations += 1;
if expanded_value < reflected_value {
simplex[worst_idx] = expanded;
values[worst_idx] = expanded_value;
} else {
simplex[worst_idx] = reflected;
values[worst_idx] = reflected_value;
}
} else {
let rho = 0.5;
let mut contracted = vec![0.0; n];
for j in 0..n {
contracted[j] = centroid[j] + rho * (simplex[worst_idx][j] - centroid[j]);
contracted[j] = contracted[j].max(bounds[j].0).min(bounds[j].1);
}
let contracted_value = objective.evaluate(&contracted);
evaluations += 1;
if contracted_value < values[worst_idx] {
simplex[worst_idx] = contracted;
values[worst_idx] = contracted_value;
} else {
let sigma = 0.5;
for i in 1..=n {
for j in 0..n {
simplex[i][j] = simplex[indices[0]][j]
+ sigma * (simplex[i][j] - simplex[indices[0]][j]);
simplex[i][j] = simplex[i][j].max(bounds[j].0).min(bounds[j].1);
}
values[i] = objective.evaluate(&simplex[i]);
evaluations += 1;
}
}
}
if let Some(callback) = &self.config.progress_callback {
callback(iteration, best_value);
}
}
let mut best_idx = 0;
let mut best_value = values[0];
for i in 1..values.len() {
if values[i] < best_value {
best_value = values[i];
best_idx = i;
}
}
Ok((simplex[best_idx].clone(), best_value, evaluations, false))
}
fn optimize_simulated_annealing(
&self,
objective: Arc<dyn ObjectiveFunction>,
initial_params: &[f64],
initial_temperature: f64,
cooling_rate: f64,
min_temperature: f64,
) -> QuantRS2Result<(Vec<f64>, f64, usize, bool)> {
use scirs2_core::random::prelude::*;
let mut rng = thread_rng();
let mut current_params = initial_params.to_vec();
let mut current_value = objective.evaluate(¤t_params);
let mut best_params = current_params.clone();
let mut best_value = current_value;
let mut temperature = initial_temperature;
let mut evaluations = 1;
let bounds = objective.bounds();
for iteration in 0..self.config.max_evaluations {
if temperature < min_temperature {
break;
}
let mut neighbor_params = current_params.clone();
for i in 0..neighbor_params.len() {
let range = bounds[i].1 - bounds[i].0;
let step = rng.random_range(-0.1..0.1) * range * temperature / initial_temperature;
neighbor_params[i] = (neighbor_params[i] + step)
.max(bounds[i].0)
.min(bounds[i].1);
}
let neighbor_value = objective.evaluate(&neighbor_params);
evaluations += 1;
let delta = neighbor_value - current_value;
if delta < 0.0 || rng.random::<f64>() < (-delta / temperature).exp() {
current_params = neighbor_params;
current_value = neighbor_value;
if current_value < best_value {
best_params.clone_from(¤t_params);
best_value = current_value;
}
}
self.record_iteration(¤t_params, current_value, iteration);
temperature *= cooling_rate;
if let Some(callback) = &self.config.progress_callback {
callback(iteration, best_value);
}
}
Ok((
best_params,
best_value,
evaluations,
temperature < min_temperature,
))
}
fn optimize_bayesian(
&self,
objective: Arc<dyn ObjectiveFunction>,
initial_params: &[f64],
acquisition_function: &AcquisitionFunction,
kernel: &KernelType,
num_initial_samples: usize,
) -> QuantRS2Result<(Vec<f64>, f64, usize, bool)> {
self.optimize_nelder_mead(
objective,
initial_params,
self.config.max_evaluations,
self.config.tolerance,
)
}
fn numerical_gradient(
&self,
objective: &dyn ObjectiveFunction,
params: &[f64],
) -> QuantRS2Result<Vec<f64>> {
let epsilon = 1e-8;
let mut gradient = vec![0.0; params.len()];
for i in 0..params.len() {
let mut params_plus = params.to_vec();
let mut params_minus = params.to_vec();
params_plus[i] += epsilon;
params_minus[i] -= epsilon;
let f_plus = objective.evaluate(¶ms_plus);
let f_minus = objective.evaluate(¶ms_minus);
gradient[i] = (f_plus - f_minus) / (2.0 * epsilon);
}
Ok(gradient)
}
fn record_iteration(&self, params: &[f64], value: f64, iteration: usize) {
if let Ok(mut history) = self.history.lock() {
history.parameters.push(params.to_vec());
history.objective_values.push(value);
history.gradient_norms.push(0.0); history.step_sizes.push(0.0); history.timestamps.push(std::time::Instant::now());
}
}
#[must_use]
pub fn get_best_parameters(&self) -> Option<Vec<f64>> {
self.best_parameters.lock().ok().and_then(|g| g.clone())
}
#[must_use]
pub fn get_best_value(&self) -> f64 {
self.best_value.lock().ok().map_or(f64::INFINITY, |g| *g)
}
pub fn build_circuit(&self, parameters: &[f64]) -> QuantRS2Result<Circuit<32>> {
if parameters.len() != self.circuit_template.parameters.len() {
return Err(QuantRS2Error::InvalidInput(
"Parameter count mismatch".to_string(),
));
}
let mut circuit = Circuit::<32>::new();
for gate_template in &self.circuit_template.structure {
}
Ok(circuit)
}
}
pub struct VQEObjective {
hamiltonian: SparseMatrix,
circuit_template: CircuitTemplate,
bounds: Vec<(f64, f64)>,
}
impl VQEObjective {
#[must_use]
pub fn new(hamiltonian: SparseMatrix, circuit_template: CircuitTemplate) -> Self {
let bounds = circuit_template
.parameters
.iter()
.map(|p| (p.lower_bound, p.upper_bound))
.collect();
Self {
hamiltonian,
circuit_template,
bounds,
}
}
}
impl ObjectiveFunction for VQEObjective {
fn evaluate(&self, parameters: &[f64]) -> f64 {
parameters.iter().map(|x| x * x).sum::<f64>()
}
fn bounds(&self) -> Vec<(f64, f64)> {
self.bounds.clone()
}
fn name(&self) -> &'static str {
"VQE"
}
}
pub struct QAOAObjective {
problem_hamiltonian: SparseMatrix,
mixer_hamiltonian: SparseMatrix,
num_layers: usize,
bounds: Vec<(f64, f64)>,
}
impl QAOAObjective {
#[must_use]
pub fn new(
problem_hamiltonian: SparseMatrix,
mixer_hamiltonian: SparseMatrix,
num_layers: usize,
) -> Self {
let bounds = vec![(0.0, 2.0 * std::f64::consts::PI); 2 * num_layers];
Self {
problem_hamiltonian,
mixer_hamiltonian,
num_layers,
bounds,
}
}
}
impl ObjectiveFunction for QAOAObjective {
fn evaluate(&self, parameters: &[f64]) -> f64 {
parameters.iter().map(|x| x.sin().powi(2)).sum::<f64>()
}
fn bounds(&self) -> Vec<(f64, f64)> {
self.bounds.clone()
}
fn name(&self) -> &'static str {
"QAOA"
}
}
impl Default for OptimizationConfig {
fn default() -> Self {
Self {
algorithm: OptimizationAlgorithm::Adam {
learning_rate: 0.01,
beta1: 0.9,
beta2: 0.999,
epsilon: 1e-8,
},
max_evaluations: 1000,
tolerance: 1e-6,
seed: None,
parallel: false,
num_threads: None,
progress_callback: None,
early_stopping: None,
}
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_optimization_config_creation() {
let config = OptimizationConfig::default();
assert_eq!(config.max_evaluations, 1000);
assert_eq!(config.tolerance, 1e-6);
}
#[test]
fn test_vqe_objective() {
let hamiltonian = SparseMatrix::identity(4);
let template = CircuitTemplate {
structure: Vec::new(),
parameters: vec![Parameter {
name: "theta".to_string(),
lower_bound: 0.0,
upper_bound: 2.0 * std::f64::consts::PI,
initial_value: 0.5,
discrete: false,
}],
num_qubits: 2,
};
let objective = VQEObjective::new(hamiltonian, template);
let value = objective.evaluate(&[0.5]);
assert!(value >= 0.0);
}
#[test]
fn test_qaoa_objective() {
let problem_h = SparseMatrix::identity(4);
let mixer_h = SparseMatrix::identity(4);
let objective = QAOAObjective::new(problem_h, mixer_h, 2);
assert_eq!(objective.bounds().len(), 4);
let value = objective.evaluate(&[0.5, 1.0, 1.5, 2.0]);
assert!(value >= 0.0);
}
#[test]
fn test_circuit_template() {
let template = CircuitTemplate {
structure: vec![ParameterizedGate {
gate_name: "RY".to_string(),
qubits: vec![0],
parameter_indices: vec![0],
fixed_parameters: Vec::new(),
}],
parameters: vec![Parameter {
name: "theta".to_string(),
lower_bound: 0.0,
upper_bound: 2.0 * std::f64::consts::PI,
initial_value: 0.0,
discrete: false,
}],
num_qubits: 1,
};
assert_eq!(template.parameters.len(), 1);
assert_eq!(template.structure.len(), 1);
}
struct TestObjective;
impl ObjectiveFunction for TestObjective {
fn evaluate(&self, parameters: &[f64]) -> f64 {
parameters.iter().map(|x| (x - 1.0).powi(2)).sum()
}
fn bounds(&self) -> Vec<(f64, f64)> {
vec![(-5.0, 5.0); 2]
}
fn name(&self) -> &'static str {
"test"
}
}
#[test]
fn test_optimizer_creation() {
let template = CircuitTemplate {
structure: Vec::new(),
parameters: vec![
Parameter {
name: "x1".to_string(),
lower_bound: -5.0,
upper_bound: 5.0,
initial_value: 0.0,
discrete: false,
},
Parameter {
name: "x2".to_string(),
lower_bound: -5.0,
upper_bound: 5.0,
initial_value: 0.0,
discrete: false,
},
],
num_qubits: 1,
};
let config = OptimizationConfig::default();
let optimizer = QuantumCircuitOptimizer::new(template, config);
assert_eq!(optimizer.circuit_template.parameters.len(), 2);
}
}