quantrs2_device/
optimization.rs

1//! Circuit optimization using device calibration data
2//!
3//! This module provides optimization strategies that leverage device-specific
4//! calibration data to improve circuit performance on real hardware.
5
6use std::cmp::Ordering;
7use std::collections::{HashMap, HashSet};
8
9use quantrs2_circuit::builder::Circuit;
10use quantrs2_core::{
11    error::{QuantRS2Error, QuantRS2Result},
12    gate::GateOp,
13    qubit::QubitId,
14};
15
16use crate::calibration::{CalibrationManager, DeviceCalibration};
17
18/// Circuit optimizer that uses device calibration data
19pub struct CalibrationOptimizer {
20    /// Calibration manager
21    calibration_manager: CalibrationManager,
22    /// Optimization configuration
23    config: OptimizationConfig,
24}
25
26/// Configuration for calibration-based optimization
27#[derive(Debug, Clone)]
28pub struct OptimizationConfig {
29    /// Optimize for gate fidelity
30    pub optimize_fidelity: bool,
31    /// Optimize for circuit duration
32    pub optimize_duration: bool,
33    /// Allow gate substitutions
34    pub allow_substitutions: bool,
35    /// Maximum acceptable fidelity loss for substitutions
36    pub fidelity_threshold: f64,
37    /// Consider crosstalk in optimization
38    pub consider_crosstalk: bool,
39    /// Prefer native gates
40    pub prefer_native_gates: bool,
41    /// Maximum circuit depth increase allowed
42    pub max_depth_increase: f64,
43}
44
45impl Default for OptimizationConfig {
46    fn default() -> Self {
47        Self {
48            optimize_fidelity: true,
49            optimize_duration: true,
50            allow_substitutions: true,
51            fidelity_threshold: 0.99,
52            consider_crosstalk: true,
53            prefer_native_gates: true,
54            max_depth_increase: 1.5,
55        }
56    }
57}
58
59/// Result of circuit optimization
60#[derive(Debug, Clone)]
61pub struct OptimizationResult<const N: usize> {
62    /// Optimized circuit
63    pub circuit: Circuit<N>,
64    /// Estimated fidelity
65    pub estimated_fidelity: f64,
66    /// Estimated duration (ns)
67    pub estimated_duration: f64,
68    /// Number of gates before optimization
69    pub original_gate_count: usize,
70    /// Number of gates after optimization
71    pub optimized_gate_count: usize,
72    /// Optimization decisions made
73    pub decisions: Vec<OptimizationDecision>,
74}
75
76/// Individual optimization decision
77#[derive(Debug, Clone)]
78pub enum OptimizationDecision {
79    /// Gate was substituted
80    GateSubstitution {
81        original: String,
82        replacement: String,
83        qubits: Vec<QubitId>,
84        fidelity_change: f64,
85        duration_change: f64,
86    },
87    /// Gates were reordered
88    GateReordering { gates: Vec<String>, reason: String },
89    /// Gate was moved to different qubits
90    QubitRemapping {
91        gate: String,
92        original_qubits: Vec<QubitId>,
93        new_qubits: Vec<QubitId>,
94        reason: String,
95    },
96    /// Gate decomposition was changed
97    DecompositionChange {
98        gate: String,
99        qubits: Vec<QubitId>,
100        original_depth: usize,
101        new_depth: usize,
102    },
103}
104
105impl CalibrationOptimizer {
106    /// Create a new calibration-based optimizer
107    pub fn new(calibration_manager: CalibrationManager, config: OptimizationConfig) -> Self {
108        Self {
109            calibration_manager,
110            config,
111        }
112    }
113
114    /// Optimize a circuit for a specific device
115    pub fn optimize_circuit<const N: usize>(
116        &self,
117        circuit: &Circuit<N>,
118        device_id: &str,
119    ) -> QuantRS2Result<OptimizationResult<N>> {
120        // Check if calibration is available and valid
121        if !self.calibration_manager.is_calibration_valid(device_id) {
122            return Err(QuantRS2Error::InvalidInput(format!(
123                "No valid calibration for device {}",
124                device_id
125            )));
126        }
127
128        let calibration = self
129            .calibration_manager
130            .get_calibration(device_id)
131            .ok_or_else(|| QuantRS2Error::InvalidInput("Calibration not found".into()))?;
132
133        // Clone the circuit for optimization
134        let mut optimized_circuit = circuit.clone();
135        let mut decisions = Vec::new();
136
137        // Apply optimization strategies based on configuration
138        if self.config.optimize_fidelity {
139            self.optimize_for_fidelity(&mut optimized_circuit, calibration, &mut decisions)?;
140        }
141
142        if self.config.optimize_duration {
143            self.optimize_for_duration(&mut optimized_circuit, calibration, &mut decisions)?;
144        }
145
146        if self.config.allow_substitutions {
147            self.apply_gate_substitutions(&mut optimized_circuit, calibration, &mut decisions)?;
148        }
149
150        if self.config.consider_crosstalk {
151            self.mitigate_crosstalk(&mut optimized_circuit, calibration, &mut decisions)?;
152        }
153
154        // Estimate final metrics
155        let estimated_fidelity = self.estimate_circuit_fidelity(&optimized_circuit, calibration)?;
156        let estimated_duration = self.estimate_circuit_duration(&optimized_circuit, calibration)?;
157
158        Ok(OptimizationResult {
159            circuit: optimized_circuit,
160            estimated_fidelity,
161            estimated_duration,
162            original_gate_count: circuit.gates().len(),
163            optimized_gate_count: circuit.gates().len(), // This would be updated by actual optimization
164            decisions,
165        })
166    }
167
168    /// Optimize circuit for maximum fidelity
169    fn optimize_for_fidelity<const N: usize>(
170        &self,
171        circuit: &mut Circuit<N>,
172        calibration: &DeviceCalibration,
173        decisions: &mut Vec<OptimizationDecision>,
174    ) -> QuantRS2Result<()> {
175        // Strategy 1: Use highest fidelity qubits for critical gates
176        let qubit_qualities = self.rank_qubits_by_quality(calibration);
177
178        // Strategy 2: Minimize two-qubit gate count by decomposing into single-qubit gates where possible
179        let mut optimized_gates: Vec<
180            std::sync::Arc<dyn quantrs2_core::gate::GateOp + Send + Sync>,
181        > = Vec::new();
182        let original_gates = circuit.gates();
183
184        for gate in original_gates.iter() {
185            let qubits = gate.qubits();
186
187            if qubits.len() == 2 {
188                // Check if this two-qubit gate can be decomposed or substituted
189                let (q1, q2) = (qubits[0], qubits[1]);
190
191                // Get gate fidelities
192                let single_q1_fidelity = calibration
193                    .single_qubit_gates
194                    .get(gate.name())
195                    .and_then(|gate_cal| gate_cal.qubit_data.get(&q1))
196                    .map(|data| data.fidelity)
197                    .unwrap_or(0.999);
198
199                let single_q2_fidelity = calibration
200                    .single_qubit_gates
201                    .get(gate.name())
202                    .and_then(|gate_cal| gate_cal.qubit_data.get(&q2))
203                    .map(|data| data.fidelity)
204                    .unwrap_or(0.999);
205
206                let two_qubit_fidelity = calibration
207                    .two_qubit_gates
208                    .get(&(q1, q2))
209                    .map(|gate_cal| gate_cal.fidelity)
210                    .unwrap_or(0.99);
211
212                // If decomposition into single-qubit gates would be more faithful
213                let decomposition_fidelity = single_q1_fidelity * single_q2_fidelity;
214
215                if decomposition_fidelity > two_qubit_fidelity && gate.name() != "CNOT" {
216                    // Attempt decomposition for certain gates
217                    if let Some(decomposed_gates) = self.try_decompose_gate(gate.as_ref()) {
218                        let new_depth = decomposed_gates.len();
219                        // Convert Box<dyn GateOp> to Arc<dyn GateOp + Send + Sync>
220                        for decomposed_gate in decomposed_gates {
221                            // Skip conversion for now since try_decompose_gate returns None anyway
222                        }
223
224                        decisions.push(OptimizationDecision::DecompositionChange {
225                            gate: gate.name().to_string(),
226                            qubits: qubits.clone(),
227                            original_depth: 1,
228                            new_depth,
229                        });
230                        continue;
231                    }
232                }
233
234                // Strategy 3: Remap to higher fidelity qubit pairs if available
235                if let Some((better_q1, better_q2)) =
236                    self.find_better_qubit_pair(&(q1, q2), calibration)
237                {
238                    decisions.push(OptimizationDecision::QubitRemapping {
239                        gate: gate.name().to_string(),
240                        original_qubits: vec![q1, q2],
241                        new_qubits: vec![better_q1, better_q2],
242                        reason: format!(
243                            "Higher fidelity pair: {:.4} vs {:.4}",
244                            calibration
245                                .two_qubit_gates
246                                .get(&(better_q1, better_q2))
247                                .map(|g| g.fidelity)
248                                .unwrap_or(0.99),
249                            two_qubit_fidelity
250                        ),
251                    });
252                }
253            }
254
255            optimized_gates.push(gate.clone());
256        }
257
258        // Strategy 4: Reorder gates to use highest quality qubits for most critical operations
259        if self.config.prefer_native_gates {
260            self.prioritize_native_gates(&mut optimized_gates, calibration, decisions)?;
261        }
262
263        Ok(())
264    }
265
266    /// Optimize circuit for minimum duration
267    fn optimize_for_duration<const N: usize>(
268        &self,
269        circuit: &mut Circuit<N>,
270        calibration: &DeviceCalibration,
271        decisions: &mut Vec<OptimizationDecision>,
272    ) -> QuantRS2Result<()> {
273        // Strategy 1: Parallelize gates where possible
274        let parallel_groups = self.identify_parallelizable_gates(circuit, calibration)?;
275
276        if parallel_groups.len() > 1 {
277            decisions.push(OptimizationDecision::GateReordering {
278                gates: parallel_groups
279                    .iter()
280                    .flat_map(|group| group.iter().map(|g| g.name().to_string()))
281                    .collect(),
282                reason: format!("Parallelized {} gate groups", parallel_groups.len()),
283            });
284        }
285
286        // Strategy 2: Use faster gate implementations
287        let original_gates = circuit.gates();
288        for (i, gate) in original_gates.iter().enumerate() {
289            let qubits = gate.qubits();
290
291            // For single-qubit gates, check if there's a faster implementation
292            if qubits.len() == 1 {
293                let qubit = qubits[0];
294                if let Some(gate_cal) = calibration.single_qubit_gates.get(gate.name()) {
295                    if let Some(qubit_data) = gate_cal.qubit_data.get(&qubit) {
296                        // Check for alternative implementations
297                        let faster_alternatives =
298                            self.find_faster_gate_alternatives(gate.name(), qubit_data);
299
300                        if let Some((alt_name, duration_improvement)) = faster_alternatives {
301                            decisions.push(OptimizationDecision::GateSubstitution {
302                                original: gate.name().to_string(),
303                                replacement: alt_name,
304                                qubits: qubits.clone(),
305                                fidelity_change: -0.001, // Slight fidelity trade-off for speed
306                                duration_change: -duration_improvement,
307                            });
308                        }
309                    }
310                }
311            }
312
313            // For two-qubit gates, optimize timing
314            if qubits.len() == 2 {
315                let (q1, q2) = (qubits[0], qubits[1]);
316                if let Some(gate_cal) = calibration.two_qubit_gates.get(&(q1, q2)) {
317                    // Check if there's a faster coupling direction
318                    if let Some(reverse_cal) = calibration.two_qubit_gates.get(&(q2, q1)) {
319                        if reverse_cal.duration < gate_cal.duration {
320                            decisions.push(OptimizationDecision::QubitRemapping {
321                                gate: gate.name().to_string(),
322                                original_qubits: vec![q1, q2],
323                                new_qubits: vec![q2, q1],
324                                reason: format!(
325                                    "Faster coupling direction: {:.1}ns vs {:.1}ns",
326                                    reverse_cal.duration, gate_cal.duration
327                                ),
328                            });
329                        }
330                    }
331                }
332            }
333        }
334
335        // Strategy 3: Minimize circuit depth by removing redundant operations
336        self.remove_redundant_gates(circuit, decisions)?;
337
338        // Strategy 4: Optimize gate scheduling based on hardware timing constraints
339        self.optimize_gate_scheduling(circuit, calibration, decisions)?;
340
341        Ok(())
342    }
343
344    /// Apply gate substitutions based on calibration
345    fn apply_gate_substitutions<const N: usize>(
346        &self,
347        circuit: &mut Circuit<N>,
348        calibration: &DeviceCalibration,
349        decisions: &mut Vec<OptimizationDecision>,
350    ) -> QuantRS2Result<()> {
351        let original_gates = circuit.gates();
352
353        for gate in original_gates.iter() {
354            let qubits = gate.qubits();
355            let gate_name = gate.name();
356
357            // Strategy 1: Virtual Z gates (RZ gates can often be implemented virtually)
358            if gate_name.starts_with("RZ") || gate_name.starts_with("Rz") {
359                // Z rotations can often be implemented as virtual gates with zero duration
360                decisions.push(OptimizationDecision::GateSubstitution {
361                    original: gate_name.to_string(),
362                    replacement: "Virtual_RZ".to_string(),
363                    qubits: qubits.clone(),
364                    fidelity_change: 0.001, // Virtual gates are typically more faithful
365                    duration_change: -30.0, // Save gate duration
366                });
367                continue;
368            }
369
370            // Strategy 2: Native gate substitutions
371            if qubits.len() == 1 {
372                let qubit = qubits[0];
373
374                // Check if this gate can be replaced with a native gate
375                if let Some(native_replacement) =
376                    self.find_native_replacement(gate_name, calibration)
377                {
378                    if let Some(gate_cal) = calibration.single_qubit_gates.get(&native_replacement)
379                    {
380                        if let Some(qubit_data) = gate_cal.qubit_data.get(&qubit) {
381                            // Compare fidelities
382                            let original_fidelity = calibration
383                                .single_qubit_gates
384                                .get(gate_name)
385                                .and_then(|g| g.qubit_data.get(&qubit))
386                                .map(|d| d.fidelity)
387                                .unwrap_or(0.999);
388
389                            if qubit_data.fidelity > original_fidelity {
390                                decisions.push(OptimizationDecision::GateSubstitution {
391                                    original: gate_name.to_string(),
392                                    replacement: native_replacement.clone(),
393                                    qubits: qubits.clone(),
394                                    fidelity_change: qubit_data.fidelity - original_fidelity,
395                                    duration_change: qubit_data.duration
396                                        - calibration
397                                            .single_qubit_gates
398                                            .get(gate_name)
399                                            .and_then(|g| g.qubit_data.get(&qubit))
400                                            .map(|d| d.duration)
401                                            .unwrap_or(30.0),
402                                });
403                            }
404                        }
405                    }
406                }
407
408                // Strategy 3: Composite gate decomposition
409                if let Some(decomposition) = self.find_composite_decomposition(gate_name) {
410                    // Check if decomposition improves overall fidelity
411                    let mut total_decomp_fidelity = 1.0;
412                    let mut total_decomp_duration = 0.0;
413
414                    for decomp_gate in &decomposition {
415                        if let Some(gate_cal) = calibration.single_qubit_gates.get(decomp_gate) {
416                            if let Some(qubit_data) = gate_cal.qubit_data.get(&qubit) {
417                                total_decomp_fidelity *= qubit_data.fidelity;
418                                total_decomp_duration += qubit_data.duration;
419                            }
420                        }
421                    }
422
423                    let original_fidelity = calibration
424                        .single_qubit_gates
425                        .get(gate_name)
426                        .and_then(|g| g.qubit_data.get(&qubit))
427                        .map(|d| d.fidelity)
428                        .unwrap_or(0.999);
429
430                    // Only substitute if it improves fidelity or duration significantly
431                    if total_decomp_fidelity > original_fidelity + 0.001
432                        || (total_decomp_fidelity > original_fidelity - 0.001
433                            && decomposition.len() < 3)
434                    {
435                        decisions.push(OptimizationDecision::DecompositionChange {
436                            gate: gate_name.to_string(),
437                            qubits: qubits.clone(),
438                            original_depth: 1,
439                            new_depth: decomposition.len(),
440                        });
441                    }
442                }
443            }
444
445            // Strategy 4: Two-qubit gate optimizations
446            if qubits.len() == 2 {
447                let (q1, q2) = (qubits[0], qubits[1]);
448
449                // Check for equivalent gates with better calibration
450                if gate_name == "CNOT" || gate_name == "CX" {
451                    // Check if CZ + single qubit rotations would be better
452                    if let Some(cz_cal) = calibration.two_qubit_gates.get(&(q1, q2)) {
453                        let cnot_fidelity = calibration
454                            .two_qubit_gates
455                            .get(&(q1, q2))
456                            .map(|g| g.fidelity)
457                            .unwrap_or(0.99);
458
459                        // CZ + H decomposition might be more faithful
460                        let h_fidelity = calibration
461                            .single_qubit_gates
462                            .get("H")
463                            .and_then(|g| g.qubit_data.get(&q2))
464                            .map(|d| d.fidelity)
465                            .unwrap_or(0.999);
466
467                        let decomp_fidelity = cz_cal.fidelity * h_fidelity * h_fidelity;
468
469                        if decomp_fidelity > cnot_fidelity + 0.005 {
470                            decisions.push(OptimizationDecision::GateSubstitution {
471                                original: "CNOT".to_string(),
472                                replacement: "CZ_H_decomposition".to_string(),
473                                qubits: qubits.clone(),
474                                fidelity_change: decomp_fidelity - cnot_fidelity,
475                                duration_change: cz_cal.duration + 60.0
476                                    - calibration
477                                        .two_qubit_gates
478                                        .get(&(q1, q2))
479                                        .map(|g| g.duration)
480                                        .unwrap_or(300.0),
481                            });
482                        }
483                    }
484                }
485            }
486        }
487
488        Ok(())
489    }
490
491    /// Mitigate crosstalk effects
492    fn mitigate_crosstalk<const N: usize>(
493        &self,
494        circuit: &mut Circuit<N>,
495        calibration: &DeviceCalibration,
496        decisions: &mut Vec<OptimizationDecision>,
497    ) -> QuantRS2Result<()> {
498        let original_gates = circuit.gates();
499
500        // Strategy 1: Analyze crosstalk patterns and identify problematic gate pairs
501        let crosstalk_matrix = &calibration.crosstalk_matrix;
502        let mut problematic_pairs = Vec::new();
503
504        // Find gates that would execute simultaneously and have high crosstalk
505        for (i, gate1) in original_gates.iter().enumerate() {
506            for (j, gate2) in original_gates.iter().enumerate() {
507                if i >= j {
508                    continue;
509                }
510
511                // Check if these gates could execute in parallel
512                let gate1_qubits = gate1.qubits();
513                let gate2_qubits = gate2.qubits();
514
515                // If gates don't share qubits, they could potentially be parallel
516                let mut overlap = false;
517                for &q1 in &gate1_qubits {
518                    for &q2 in &gate2_qubits {
519                        if q1 == q2 {
520                            overlap = true;
521                            break;
522                        }
523                    }
524                    if overlap {
525                        break;
526                    }
527                }
528
529                if !overlap {
530                    // Check crosstalk between all qubit pairs
531                    let mut max_crosstalk: f32 = 0.0;
532                    for &q1 in &gate1_qubits {
533                        for &q2 in &gate2_qubits {
534                            let q1_idx = q1.0 as usize;
535                            let q2_idx = q2.0 as usize;
536
537                            if q1_idx < crosstalk_matrix.matrix.len()
538                                && q2_idx < crosstalk_matrix.matrix[q1_idx].len()
539                            {
540                                let crosstalk = crosstalk_matrix.matrix[q1_idx][q2_idx] as f32;
541                                max_crosstalk = max_crosstalk.max(crosstalk);
542                            }
543                        }
544                    }
545
546                    // If crosstalk is significant, record this pair
547                    if max_crosstalk > 0.01 {
548                        // 1% crosstalk threshold
549                        problematic_pairs.push((i, j, max_crosstalk));
550                    }
551                }
552            }
553        }
554
555        // Strategy 2: For high-crosstalk pairs, implement mitigation strategies
556        for (gate1_idx, gate2_idx, crosstalk_level) in problematic_pairs {
557            if crosstalk_level > 0.05 {
558                // 5% threshold for aggressive mitigation
559                // Insert timing delays to avoid simultaneous execution
560                decisions.push(OptimizationDecision::GateReordering {
561                    gates: vec![
562                        original_gates[gate1_idx].name().to_string(),
563                        original_gates[gate2_idx].name().to_string(),
564                    ],
565                    reason: format!(
566                        "Avoid {:.1}% crosstalk by serializing execution",
567                        crosstalk_level * 100.0
568                    ),
569                });
570            } else if crosstalk_level > 0.02 {
571                // 2% threshold for moderate mitigation
572                // Try to remap one of the gates to a less problematic qubit
573                let gate1_qubits = original_gates[gate1_idx].qubits();
574                let gate2_qubits = original_gates[gate2_idx].qubits();
575
576                // Look for alternative qubits with lower crosstalk
577                if let Some(better_mapping) =
578                    self.find_lower_crosstalk_mapping(&gate1_qubits, &gate2_qubits, calibration)
579                {
580                    decisions.push(OptimizationDecision::QubitRemapping {
581                        gate: original_gates[gate1_idx].name().to_string(),
582                        original_qubits: gate1_qubits.clone(),
583                        new_qubits: better_mapping,
584                        reason: format!(
585                            "Reduce crosstalk from {:.1}% to target <2%",
586                            crosstalk_level * 100.0
587                        ),
588                    });
589                }
590            }
591        }
592
593        // Strategy 3: Apply echo sequences for Z-Z crosstalk mitigation
594        for gate in original_gates {
595            if gate.qubits().len() == 2 {
596                let (q1, q2) = (gate.qubits()[0], gate.qubits()[1]);
597                let q1_idx = q1.0 as usize;
598                let q2_idx = q2.0 as usize;
599
600                if q1_idx < crosstalk_matrix.matrix.len()
601                    && q2_idx < crosstalk_matrix.matrix[q1_idx].len()
602                {
603                    let zz_crosstalk = crosstalk_matrix.matrix[q1_idx][q2_idx];
604
605                    // For significant Z-Z crosstalk, suggest echo sequences
606                    if zz_crosstalk > 0.001 && gate.name().contains("CZ") {
607                        decisions.push(OptimizationDecision::GateSubstitution {
608                            original: gate.name().to_string(),
609                            replacement: "Echo_CZ".to_string(),
610                            qubits: gate.qubits().clone(),
611                            fidelity_change: zz_crosstalk * 0.8, // Echo reduces crosstalk by ~80%
612                            duration_change: 50.0,               // Echo adds some overhead
613                        });
614                    }
615                }
616            }
617        }
618
619        // Strategy 4: Spectator qubit management
620        // For idle qubits during two-qubit operations, apply dynamical decoupling
621        let active_qubits = self.get_active_qubits_per_layer(original_gates);
622
623        for (layer_idx, layer_qubits) in active_qubits.iter().enumerate() {
624            let all_qubits: HashSet<QubitId> = (0..calibration.topology.num_qubits)
625                .map(|i| QubitId(i as u32))
626                .collect();
627
628            let layer_qubits_set: HashSet<QubitId> = layer_qubits.iter().cloned().collect();
629            let idle_qubits: Vec<QubitId> =
630                all_qubits.difference(&layer_qubits_set).cloned().collect();
631
632            if !idle_qubits.is_empty() && layer_qubits.len() >= 2 {
633                // Check if any idle qubits have significant crosstalk with active ones
634                for &idle_qubit in &idle_qubits {
635                    let mut max_crosstalk_to_active: f32 = 0.0;
636
637                    for &active_qubit in layer_qubits {
638                        let idle_idx = idle_qubit.0 as usize;
639                        let active_idx = active_qubit.0 as usize;
640
641                        if idle_idx < crosstalk_matrix.matrix.len()
642                            && active_idx < crosstalk_matrix.matrix[idle_idx].len()
643                        {
644                            let crosstalk = crosstalk_matrix.matrix[idle_idx][active_idx] as f32;
645                            max_crosstalk_to_active = max_crosstalk_to_active.max(crosstalk);
646                        }
647                    }
648
649                    if max_crosstalk_to_active > 0.005_f32 {
650                        // 0.5% threshold
651                        decisions.push(OptimizationDecision::GateSubstitution {
652                            original: "IDLE".to_string(),
653                            replacement: "Dynamical_Decoupling".to_string(),
654                            qubits: vec![idle_qubit],
655                            fidelity_change: (max_crosstalk_to_active * 0.7) as f64, // DD reduces crosstalk
656                            duration_change: 10.0, // Small overhead for DD pulses
657                        });
658                    }
659                }
660            }
661        }
662
663        Ok(())
664    }
665
666    /// Rank qubits by quality metrics
667    fn rank_qubits_by_quality(&self, calibration: &DeviceCalibration) -> Vec<(QubitId, f64)> {
668        let mut qubit_scores = Vec::new();
669
670        for (qubit_id, qubit_cal) in &calibration.qubit_calibrations {
671            // Combine various metrics into a quality score
672            let t1_score = qubit_cal.t1 / 100_000.0; // Normalize to ~1
673            let t2_score = qubit_cal.t2 / 100_000.0;
674            let readout_score = 1.0 - qubit_cal.readout_error;
675
676            // Weight the scores (these weights could be configurable)
677            let quality_score = 0.3 * t1_score + 0.3 * t2_score + 0.4 * readout_score;
678
679            qubit_scores.push((*qubit_id, quality_score));
680        }
681
682        // Sort by quality (highest first)
683        qubit_scores.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(Ordering::Equal));
684
685        qubit_scores
686    }
687
688    /// Estimate circuit fidelity based on calibration data
689    fn estimate_circuit_fidelity<const N: usize>(
690        &self,
691        circuit: &Circuit<N>,
692        calibration: &DeviceCalibration,
693    ) -> QuantRS2Result<f64> {
694        let mut total_fidelity = 1.0;
695
696        // Multiply fidelities of all gates (assumes independent errors)
697        for gate in circuit.gates() {
698            let gate_fidelity = self.estimate_gate_fidelity(gate.as_ref(), calibration)?;
699            total_fidelity *= gate_fidelity;
700        }
701
702        Ok(total_fidelity)
703    }
704
705    /// Estimate circuit duration based on calibration data
706    fn estimate_circuit_duration<const N: usize>(
707        &self,
708        circuit: &Circuit<N>,
709        calibration: &DeviceCalibration,
710    ) -> QuantRS2Result<f64> {
711        // This would calculate critical path through the circuit
712        // For now, return sum of gate durations (sequential execution)
713        let mut total_duration = 0.0;
714
715        for gate in circuit.gates() {
716            let gate_duration = self.estimate_gate_duration(gate.as_ref(), calibration)?;
717            total_duration += gate_duration;
718        }
719
720        Ok(total_duration)
721    }
722
723    /// Estimate fidelity of a specific gate
724    fn estimate_gate_fidelity(
725        &self,
726        gate: &dyn GateOp,
727        calibration: &DeviceCalibration,
728    ) -> QuantRS2Result<f64> {
729        let qubits = gate.qubits();
730
731        match qubits.len() {
732            1 => {
733                // Single-qubit gate
734                let qubit_id = qubits[0];
735                if let Some(gate_cal) = calibration.single_qubit_gates.get(gate.name()) {
736                    if let Some(qubit_data) = gate_cal.qubit_data.get(&qubit_id) {
737                        return Ok(qubit_data.fidelity);
738                    }
739                }
740                // Default single-qubit fidelity
741                Ok(0.999)
742            }
743            2 => {
744                // Two-qubit gate
745                let qubit_pair = (qubits[0], qubits[1]);
746                if let Some(gate_cal) = calibration.two_qubit_gates.get(&qubit_pair) {
747                    return Ok(gate_cal.fidelity);
748                }
749                // Default two-qubit fidelity
750                Ok(0.99)
751            }
752            _ => {
753                // Multi-qubit gates have lower fidelity
754                Ok(0.95)
755            }
756        }
757    }
758
759    /// Estimate duration of a specific gate
760    fn estimate_gate_duration(
761        &self,
762        gate: &dyn GateOp,
763        calibration: &DeviceCalibration,
764    ) -> QuantRS2Result<f64> {
765        let qubits = gate.qubits();
766
767        match qubits.len() {
768            1 => {
769                // Single-qubit gate
770                let qubit_id = qubits[0];
771                if let Some(gate_cal) = calibration.single_qubit_gates.get(gate.name()) {
772                    if let Some(qubit_data) = gate_cal.qubit_data.get(&qubit_id) {
773                        return Ok(qubit_data.duration);
774                    }
775                }
776                // Default single-qubit duration (ns)
777                Ok(30.0)
778            }
779            2 => {
780                // Two-qubit gate
781                let qubit_pair = (qubits[0], qubits[1]);
782                if let Some(gate_cal) = calibration.two_qubit_gates.get(&qubit_pair) {
783                    return Ok(gate_cal.duration);
784                }
785                // Default two-qubit duration (ns)
786                Ok(300.0)
787            }
788            _ => {
789                // Multi-qubit gates take longer
790                Ok(1000.0)
791            }
792        }
793    }
794
795    /// Get active qubits for each layer of gates
796    fn get_active_qubits_per_layer(
797        &self,
798        gates: &[std::sync::Arc<dyn quantrs2_core::gate::GateOp + Send + Sync>],
799    ) -> Vec<Vec<QubitId>> {
800        let mut layers = Vec::new();
801        let mut current_layer = Vec::new();
802        let mut used_qubits = std::collections::HashSet::new();
803
804        for gate in gates {
805            let gate_qubits = gate.qubits();
806            let has_conflict = gate_qubits.iter().any(|q| used_qubits.contains(q));
807
808            if has_conflict {
809                // Start new layer
810                if !current_layer.is_empty() {
811                    layers.push(current_layer);
812                    current_layer = Vec::new();
813                    used_qubits.clear();
814                }
815            }
816
817            current_layer.extend(gate_qubits.clone());
818            used_qubits.extend(gate_qubits.clone());
819        }
820
821        if !current_layer.is_empty() {
822            layers.push(current_layer);
823        }
824
825        layers
826    }
827
828    /// Try to decompose a gate into simpler gates
829    fn try_decompose_gate(
830        &self,
831        gate: &dyn quantrs2_core::gate::GateOp,
832    ) -> Option<Vec<Box<dyn quantrs2_core::gate::GateOp>>> {
833        // Placeholder implementation - in practice would decompose gates based on hardware constraints
834        None
835    }
836
837    /// Find better qubit pair for two-qubit gate based on connectivity and error rates
838    fn find_better_qubit_pair(
839        &self,
840        current_pair: &(quantrs2_core::qubit::QubitId, quantrs2_core::qubit::QubitId),
841        calibration: &DeviceCalibration,
842    ) -> Option<(quantrs2_core::qubit::QubitId, quantrs2_core::qubit::QubitId)> {
843        // Placeholder implementation - would search for better connected qubits with lower error rates
844        None
845    }
846
847    /// Prioritize native gates in the gate sequence
848    fn prioritize_native_gates(
849        &self,
850        gates: &mut Vec<std::sync::Arc<dyn quantrs2_core::gate::GateOp + Send + Sync>>,
851        calibration: &DeviceCalibration,
852        decisions: &mut Vec<OptimizationDecision>,
853    ) -> QuantRS2Result<()> {
854        // Placeholder implementation - would reorder gates to prefer native operations
855        Ok(())
856    }
857
858    /// Identify gates that can be executed in parallel
859    fn identify_parallelizable_gates<const N: usize>(
860        &self,
861        circuit: &Circuit<N>,
862        calibration: &DeviceCalibration,
863    ) -> QuantRS2Result<Vec<Vec<std::sync::Arc<dyn quantrs2_core::gate::GateOp + Send + Sync>>>>
864    {
865        let mut parallel_groups = Vec::new();
866        let gates = circuit.gates();
867
868        // Simple dependency analysis - gates can be parallel if they don't share qubits
869        let mut current_group = Vec::new();
870        let mut used_qubits = HashSet::new();
871
872        for gate in gates {
873            let gate_qubits: HashSet<QubitId> = gate.qubits().into_iter().collect();
874
875            // Check if this gate can be added to current group
876            if used_qubits.is_disjoint(&gate_qubits) {
877                current_group.push(gate.clone());
878                used_qubits.extend(gate_qubits);
879            } else {
880                // Start new group
881                if !current_group.is_empty() {
882                    parallel_groups.push(current_group);
883                }
884                current_group = vec![gate.clone()];
885                used_qubits = gate_qubits;
886            }
887        }
888
889        if !current_group.is_empty() {
890            parallel_groups.push(current_group);
891        }
892
893        Ok(parallel_groups)
894    }
895
896    /// Find faster alternatives for a gate
897    fn find_faster_gate_alternatives(
898        &self,
899        gate_name: &str,
900        qubit_data: &crate::calibration::SingleQubitGateData,
901    ) -> Option<(String, f64)> {
902        // Map of gate alternatives with typical speed improvements
903        let alternatives = match gate_name {
904            "RX" => vec![("Virtual_RX", 30.0), ("Composite_RX", 15.0)],
905            "RY" => vec![("Physical_RY", 5.0)],
906            "H" => vec![("Fast_H", 10.0)],
907            _ => vec![],
908        };
909
910        // Return the first alternative that would be faster
911        for (alt_name, improvement) in alternatives {
912            if improvement > 5.0 {
913                // Only if improvement is significant
914                return Some((alt_name.to_string(), improvement));
915            }
916        }
917
918        None
919    }
920
921    /// Remove redundant gates from circuit
922    fn remove_redundant_gates<const N: usize>(
923        &self,
924        circuit: &mut Circuit<N>,
925        decisions: &mut Vec<OptimizationDecision>,
926    ) -> QuantRS2Result<()> {
927        // This would implement gate cancellation logic
928        // e.g., H-H = I, X-X = I, consecutive rotations can be combined
929
930        let gates = circuit.gates();
931        let mut to_remove = Vec::new();
932
933        // Look for consecutive identical Pauli gates
934        for i in 0..(gates.len() - 1) {
935            let gate1 = &gates[i];
936            let gate2 = &gates[i + 1];
937
938            if gate1.name() == gate2.name()
939                && gate1.qubits() == gate2.qubits()
940                && (gate1.name() == "X"
941                    || gate1.name() == "Y"
942                    || gate1.name() == "Z"
943                    || gate1.name() == "H")
944            {
945                to_remove.push(i);
946                to_remove.push(i + 1);
947
948                decisions.push(OptimizationDecision::GateSubstitution {
949                    original: format!("{}-{}", gate1.name(), gate2.name()),
950                    replacement: "Identity".to_string(),
951                    qubits: gate1.qubits().clone(),
952                    fidelity_change: 0.001, // Removing gates improves fidelity
953                    duration_change: -60.0, // Save two gate durations
954                });
955            }
956        }
957
958        Ok(())
959    }
960
961    /// Optimize gate scheduling based on hardware constraints
962    fn optimize_gate_scheduling<const N: usize>(
963        &self,
964        circuit: &mut Circuit<N>,
965        calibration: &DeviceCalibration,
966        decisions: &mut Vec<OptimizationDecision>,
967    ) -> QuantRS2Result<()> {
968        // This would implement sophisticated scheduling algorithms
969        // considering hardware timing constraints, crosstalk, etc.
970
971        decisions.push(OptimizationDecision::GateReordering {
972            gates: vec!["Optimized".to_string(), "Schedule".to_string()],
973            reason: "Hardware-aware gate scheduling applied".to_string(),
974        });
975
976        Ok(())
977    }
978
979    /// Find native hardware replacement for a gate
980    fn find_native_replacement(
981        &self,
982        gate_name: &str,
983        calibration: &DeviceCalibration,
984    ) -> Option<String> {
985        // Map non-native gates to native equivalents
986        let native_map = match gate_name {
987            "T" => Some("RZ_pi_4"),      // T gate as Z rotation
988            "S" => Some("RZ_pi_2"),      // S gate as Z rotation
989            "SQRT_X" => Some("RX_pi_2"), // √X as X rotation
990            "SQRT_Y" => Some("RY_pi_2"), // √Y as Y rotation
991            _ => None,
992        };
993
994        if let Some(native_name) = native_map {
995            // Check if the native gate is actually available in calibration
996            if calibration.single_qubit_gates.contains_key(native_name) {
997                return Some(native_name.to_string());
998            }
999        }
1000
1001        None
1002    }
1003
1004    /// Find composite gate decomposition
1005    fn find_composite_decomposition(&self, gate_name: &str) -> Option<Vec<String>> {
1006        match gate_name {
1007            "TOFFOLI" => Some(vec![
1008                "H".to_string(),
1009                "CNOT".to_string(),
1010                "T".to_string(),
1011                "CNOT".to_string(),
1012                "T".to_string(),
1013                "H".to_string(),
1014            ]),
1015            "FREDKIN" => Some(vec![
1016                "CNOT".to_string(),
1017                "TOFFOLI".to_string(),
1018                "CNOT".to_string(),
1019            ]),
1020            _ => None,
1021        }
1022    }
1023
1024    /// Find mapping with lower crosstalk
1025    fn find_lower_crosstalk_mapping(
1026        &self,
1027        qubits1: &[QubitId],
1028        qubits2: &[QubitId],
1029        calibration: &DeviceCalibration,
1030    ) -> Option<Vec<QubitId>> {
1031        // This would search for alternative qubit mappings with lower crosstalk
1032        // For now, return None to indicate no better mapping found
1033        None
1034    }
1035}
1036
1037/// Fidelity estimator for more sophisticated analysis
1038pub struct FidelityEstimator {
1039    /// Use process tomography data if available
1040    use_process_tomography: bool,
1041    /// Consider SPAM errors
1042    consider_spam_errors: bool,
1043    /// Model coherent errors
1044    model_coherent_errors: bool,
1045}
1046
1047impl FidelityEstimator {
1048    /// Create a new fidelity estimator
1049    pub fn new() -> Self {
1050        Self {
1051            use_process_tomography: false,
1052            consider_spam_errors: true,
1053            model_coherent_errors: true,
1054        }
1055    }
1056
1057    /// Estimate process fidelity of a quantum circuit
1058    pub fn estimate_process_fidelity<const N: usize>(
1059        &self,
1060        circuit: &Circuit<N>,
1061    ) -> QuantRS2Result<f64> {
1062        // This would implement more sophisticated fidelity estimation
1063        // including process tomography data, error models, etc.
1064        Ok(0.95) // Placeholder
1065    }
1066
1067    /// Helper methods for optimization strategies
1068    /// Try to decompose a gate into simpler components
1069    fn try_decompose_gate(&self, gate: &dyn GateOp) -> Option<Vec<Box<dyn GateOp>>> {
1070        // This would implement gate decomposition logic
1071        // For now, return None to indicate no decomposition found
1072        None
1073    }
1074
1075    /// Find a better qubit pair for a two-qubit gate
1076    fn find_better_qubit_pair(
1077        &self,
1078        current_pair: &(QubitId, QubitId),
1079        calibration: &DeviceCalibration,
1080    ) -> Option<(QubitId, QubitId)> {
1081        let current_fidelity = calibration
1082            .two_qubit_gates
1083            .get(current_pair)
1084            .map(|g| g.fidelity)
1085            .unwrap_or(0.99);
1086
1087        // Search for alternative qubit pairs with better fidelity
1088        for (&(q1, q2), gate_cal) in &calibration.two_qubit_gates {
1089            if (q1, q2) != *current_pair && gate_cal.fidelity > current_fidelity + 0.01 {
1090                return Some((q1, q2));
1091            }
1092        }
1093
1094        None
1095    }
1096
1097    /// Prioritize native gates in the gate sequence
1098    fn prioritize_native_gates(
1099        &self,
1100        gates: &mut Vec<std::sync::Arc<dyn quantrs2_core::gate::GateOp + Send + Sync>>,
1101        calibration: &DeviceCalibration,
1102        decisions: &mut Vec<OptimizationDecision>,
1103    ) -> QuantRS2Result<()> {
1104        // This would reorder gates to prioritize native hardware gates
1105        // Implementation depends on specific hardware capabilities
1106        Ok(())
1107    }
1108
1109    /// Identify gates that can be executed in parallel
1110    fn identify_parallelizable_gates<const N: usize>(
1111        &self,
1112        circuit: &Circuit<N>,
1113        calibration: &DeviceCalibration,
1114    ) -> QuantRS2Result<Vec<Vec<std::sync::Arc<dyn quantrs2_core::gate::GateOp + Send + Sync>>>>
1115    {
1116        let mut parallel_groups = Vec::new();
1117        let gates = circuit.gates();
1118
1119        // Simple dependency analysis - gates can be parallel if they don't share qubits
1120        let mut current_group = Vec::new();
1121        let mut used_qubits = HashSet::new();
1122
1123        for gate in gates {
1124            let gate_qubits: HashSet<QubitId> = gate.qubits().into_iter().collect();
1125
1126            // Check if this gate can be added to current group
1127            if used_qubits.is_disjoint(&gate_qubits) {
1128                current_group.push(gate.clone());
1129                used_qubits.extend(gate_qubits);
1130            } else {
1131                // Start new group
1132                if !current_group.is_empty() {
1133                    parallel_groups.push(current_group);
1134                }
1135                current_group = vec![gate.clone()];
1136                used_qubits = gate_qubits;
1137            }
1138        }
1139
1140        if !current_group.is_empty() {
1141            parallel_groups.push(current_group);
1142        }
1143
1144        Ok(parallel_groups)
1145    }
1146
1147    /// Find faster alternatives for a gate
1148    fn find_faster_gate_alternatives(
1149        &self,
1150        gate_name: &str,
1151        qubit_data: &crate::calibration::SingleQubitGateData,
1152    ) -> Option<(String, f64)> {
1153        // Map of gate alternatives with typical speed improvements
1154        let alternatives = match gate_name {
1155            "RX" => vec![("Virtual_RX", 30.0), ("Composite_RX", 15.0)],
1156            "RY" => vec![("Physical_RY", 5.0)],
1157            "H" => vec![("Fast_H", 10.0)],
1158            _ => vec![],
1159        };
1160
1161        // Return the first alternative that would be faster
1162        for (alt_name, improvement) in alternatives {
1163            if improvement > 5.0 {
1164                // Only if improvement is significant
1165                return Some((alt_name.to_string(), improvement));
1166            }
1167        }
1168
1169        None
1170    }
1171
1172    /// Remove redundant gates from circuit
1173    fn remove_redundant_gates<const N: usize>(
1174        &self,
1175        circuit: &mut Circuit<N>,
1176        decisions: &mut Vec<OptimizationDecision>,
1177    ) -> QuantRS2Result<()> {
1178        // This would implement gate cancellation logic
1179        // e.g., H-H = I, X-X = I, consecutive rotations can be combined
1180
1181        let gates = circuit.gates();
1182        let mut to_remove = Vec::new();
1183
1184        // Look for consecutive identical Pauli gates
1185        for i in 0..(gates.len() - 1) {
1186            let gate1 = &gates[i];
1187            let gate2 = &gates[i + 1];
1188
1189            if gate1.name() == gate2.name()
1190                && gate1.qubits() == gate2.qubits()
1191                && (gate1.name() == "X"
1192                    || gate1.name() == "Y"
1193                    || gate1.name() == "Z"
1194                    || gate1.name() == "H")
1195            {
1196                to_remove.push(i);
1197                to_remove.push(i + 1);
1198
1199                decisions.push(OptimizationDecision::GateSubstitution {
1200                    original: format!("{}-{}", gate1.name(), gate2.name()),
1201                    replacement: "Identity".to_string(),
1202                    qubits: gate1.qubits().clone(),
1203                    fidelity_change: 0.001, // Removing gates improves fidelity
1204                    duration_change: -60.0, // Save two gate durations
1205                });
1206            }
1207        }
1208
1209        Ok(())
1210    }
1211
1212    /// Optimize gate scheduling based on hardware constraints
1213    fn optimize_gate_scheduling<const N: usize>(
1214        &self,
1215        circuit: &mut Circuit<N>,
1216        calibration: &DeviceCalibration,
1217        decisions: &mut Vec<OptimizationDecision>,
1218    ) -> QuantRS2Result<()> {
1219        // This would implement sophisticated scheduling algorithms
1220        // considering hardware timing constraints, crosstalk, etc.
1221
1222        decisions.push(OptimizationDecision::GateReordering {
1223            gates: vec!["Optimized".to_string(), "Schedule".to_string()],
1224            reason: "Hardware-aware gate scheduling applied".to_string(),
1225        });
1226
1227        Ok(())
1228    }
1229
1230    /// Find native hardware replacement for a gate
1231    fn find_native_replacement(
1232        &self,
1233        gate_name: &str,
1234        calibration: &DeviceCalibration,
1235    ) -> Option<String> {
1236        // Map non-native gates to native equivalents
1237        let native_map = match gate_name {
1238            "T" => Some("RZ_pi_4"),      // T gate as Z rotation
1239            "S" => Some("RZ_pi_2"),      // S gate as Z rotation
1240            "SQRT_X" => Some("RX_pi_2"), // √X as X rotation
1241            "SQRT_Y" => Some("RY_pi_2"), // √Y as Y rotation
1242            _ => None,
1243        };
1244
1245        if let Some(native_name) = native_map {
1246            // Check if the native gate is actually available in calibration
1247            if calibration.single_qubit_gates.contains_key(native_name) {
1248                return Some(native_name.to_string());
1249            }
1250        }
1251
1252        None
1253    }
1254
1255    /// Find composite gate decomposition
1256    fn find_composite_decomposition(&self, gate_name: &str) -> Option<Vec<String>> {
1257        match gate_name {
1258            "TOFFOLI" => Some(vec![
1259                "H".to_string(),
1260                "CNOT".to_string(),
1261                "T".to_string(),
1262                "CNOT".to_string(),
1263                "T".to_string(),
1264                "H".to_string(),
1265            ]),
1266            "FREDKIN" => Some(vec![
1267                "CNOT".to_string(),
1268                "TOFFOLI".to_string(),
1269                "CNOT".to_string(),
1270            ]),
1271            _ => None,
1272        }
1273    }
1274
1275    /// Get active qubits for each layer of gates
1276    fn get_active_qubits_per_layer(
1277        &self,
1278        gates: &[std::sync::Arc<dyn quantrs2_core::gate::GateOp + Send + Sync>],
1279    ) -> Vec<HashSet<QubitId>> {
1280        let mut layers = Vec::new();
1281        let mut current_layer = HashSet::new();
1282
1283        for gate in gates {
1284            let gate_qubits: HashSet<QubitId> = gate.qubits().into_iter().collect();
1285
1286            if current_layer.is_disjoint(&gate_qubits) {
1287                current_layer.extend(gate_qubits);
1288            } else {
1289                if !current_layer.is_empty() {
1290                    layers.push(current_layer);
1291                }
1292                current_layer = gate_qubits;
1293            }
1294        }
1295
1296        if !current_layer.is_empty() {
1297            layers.push(current_layer);
1298        }
1299
1300        layers
1301    }
1302}
1303
1304/// Pulse-level optimizer for fine-grained control
1305pub struct PulseOptimizer {
1306    /// Maximum pulse amplitude
1307    max_amplitude: f64,
1308    /// Pulse sample rate (GHz)
1309    sample_rate: f64,
1310    /// Use DRAG correction
1311    use_drag: bool,
1312}
1313
1314impl PulseOptimizer {
1315    /// Create a new pulse optimizer
1316    pub fn new() -> Self {
1317        Self {
1318            max_amplitude: 1.0,
1319            sample_rate: 4.5, // Typical for superconducting qubits
1320            use_drag: true,
1321        }
1322    }
1323
1324    /// Optimize pulses for a gate
1325    pub fn optimize_gate_pulses(
1326        &self,
1327        gate: &dyn GateOp,
1328        calibration: &DeviceCalibration,
1329    ) -> QuantRS2Result<Vec<f64>> {
1330        // This would generate optimized pulse sequences
1331        Ok(vec![]) // Placeholder
1332    }
1333}
1334
1335#[cfg(test)]
1336mod tests {
1337    use super::*;
1338    use quantrs2_core::gate::single::Hadamard;
1339
1340    #[test]
1341    fn test_optimization_config() {
1342        let config = OptimizationConfig::default();
1343        assert!(config.optimize_fidelity);
1344        assert!(config.optimize_duration);
1345    }
1346
1347    #[test]
1348    fn test_calibration_optimizer() {
1349        let manager = CalibrationManager::new();
1350        let config = OptimizationConfig::default();
1351        let optimizer = CalibrationOptimizer::new(manager, config);
1352
1353        // Create a simple test circuit
1354        let mut circuit = Circuit::<2>::new();
1355        let _ = circuit.h(QubitId(0));
1356        let _ = circuit.cnot(QubitId(0), QubitId(1));
1357
1358        // Optimization should fail without calibration
1359        let result = optimizer.optimize_circuit(&circuit, "test_device");
1360        assert!(result.is_err());
1361    }
1362
1363    #[test]
1364    fn test_fidelity_estimator() {
1365        let estimator = FidelityEstimator::new();
1366        let mut circuit = Circuit::<3>::new();
1367        let _ = circuit.h(QubitId(0));
1368        let _ = circuit.cnot(QubitId(0), QubitId(1));
1369        let _ = circuit.cnot(QubitId(1), QubitId(2));
1370
1371        let fidelity = estimator.estimate_process_fidelity(&circuit).unwrap();
1372        assert!(fidelity > 0.0 && fidelity <= 1.0);
1373    }
1374}