Add support for Cirq: library for writing, manipulating, and optimizing quantum circuits, and then running them on quantum computers and quantum simulators.

scientific-skills/cirq/SKILL.md

@@ -0,0 +1,340 @@
+---
+name: cirq
+description: Quantum computing framework for building, simulating, optimizing, and executing quantum circuits. Use this skill when working with quantum algorithms, quantum circuit design, quantum simulation (noiseless or noisy), running on quantum hardware (Google, IonQ, AQT, Pasqal), circuit optimization and compilation, noise modeling and characterization, or quantum experiments and benchmarking (VQE, QAOA, QPE, randomized benchmarking).
+---
+
+# Cirq - Quantum Computing with Python
+
+Cirq is Google Quantum AI's open-source framework for designing, simulating, and running quantum circuits on quantum computers and simulators.
+
+## Installation
+
+```bash
+uv pip install cirq
+```
+
+For hardware integration:
+```bash
+# Google Quantum Engine
+uv pip install cirq-google
+
+# IonQ
+uv pip install cirq-ionq
+
+# AQT (Alpine Quantum Technologies)
+uv pip install cirq-aqt
+
+# Pasqal
+uv pip install cirq-pasqal
+
+# Azure Quantum
+uv pip install azure-quantum cirq
+```
+
+## Quick Start
+
+### Basic Circuit
+
+```python
+import cirq
+import numpy as np
+
+# Create qubits
+q0, q1 = cirq.LineQubit.range(2)
+
+# Build circuit
+circuit = cirq.Circuit(
+    cirq.H(q0),              # Hadamard on q0
+    cirq.CNOT(q0, q1),       # CNOT with q0 control, q1 target
+    cirq.measure(q0, q1, key='result')
+)
+
+print(circuit)
+
+# Simulate
+simulator = cirq.Simulator()
+result = simulator.run(circuit, repetitions=1000)
+
+# Display results
+print(result.histogram(key='result'))
+```
+
+### Parameterized Circuit
+
+```python
+import sympy
+
+# Define symbolic parameter
+theta = sympy.Symbol('theta')
+
+# Create parameterized circuit
+circuit = cirq.Circuit(
+    cirq.ry(theta)(q0),
+    cirq.measure(q0, key='m')
+)
+
+# Sweep over parameter values
+sweep = cirq.Linspace('theta', start=0, stop=2*np.pi, length=20)
+results = simulator.run_sweep(circuit, params=sweep, repetitions=1000)
+
+# Process results
+for params, result in zip(sweep, results):
+    theta_val = params['theta']
+    counts = result.histogram(key='m')
+    print(f"θ={theta_val:.2f}: {counts}")
+```
+
+## Core Capabilities
+
+### Circuit Building
+For comprehensive information about building quantum circuits, including qubits, gates, operations, custom gates, and circuit patterns, see:
+- **[references/building.md](references/building.md)** - Complete guide to circuit construction
+
+Common topics:
+- Qubit types (GridQubit, LineQubit, NamedQubit)
+- Single and two-qubit gates
+- Parameterized gates and operations
+- Custom gate decomposition
+- Circuit organization with moments
+- Standard circuit patterns (Bell states, GHZ, QFT)
+- Import/export (OpenQASM, JSON)
+- Working with qudits and observables
+
+### Simulation
+For detailed information about simulating quantum circuits, including exact simulation, noisy simulation, parameter sweeps, and the Quantum Virtual Machine, see:
+- **[references/simulation.md](references/simulation.md)** - Complete guide to quantum simulation
+
+Common topics:
+- Exact simulation (state vector, density matrix)
+- Sampling and measurements
+- Parameter sweeps (single and multiple parameters)
+- Noisy simulation
+- State histograms and visualization
+- Quantum Virtual Machine (QVM)
+- Expectation values and observables
+- Performance optimization
+
+### Circuit Transformation
+For information about optimizing, compiling, and manipulating quantum circuits, see:
+- **[references/transformation.md](references/transformation.md)** - Complete guide to circuit transformations
+
+Common topics:
+- Transformer framework
+- Gate decomposition
+- Circuit optimization (merge gates, eject Z gates, drop negligible operations)
+- Circuit compilation for hardware
+- Qubit routing and SWAP insertion
+- Custom transformers
+- Transformation pipelines
+
+### Hardware Integration
+For information about running circuits on real quantum hardware from various providers, see:
+- **[references/hardware.md](references/hardware.md)** - Complete guide to hardware integration
+
+Supported providers:
+- **Google Quantum AI** (cirq-google) - Sycamore, Weber processors
+- **IonQ** (cirq-ionq) - Trapped ion quantum computers
+- **Azure Quantum** (azure-quantum) - IonQ and Honeywell backends
+- **AQT** (cirq-aqt) - Alpine Quantum Technologies
+- **Pasqal** (cirq-pasqal) - Neutral atom quantum computers
+
+Topics include device representation, qubit selection, authentication, job management, and circuit optimization for hardware.
+
+### Noise Modeling
+For information about modeling noise, noisy simulation, characterization, and error mitigation, see:
+- **[references/noise.md](references/noise.md)** - Complete guide to noise modeling
+
+Common topics:
+- Noise channels (depolarizing, amplitude damping, phase damping)
+- Noise models (constant, gate-specific, qubit-specific, thermal)
+- Adding noise to circuits
+- Readout noise
+- Noise characterization (randomized benchmarking, XEB)
+- Noise visualization (heatmaps)
+- Error mitigation techniques
+
+### Quantum Experiments
+For information about designing experiments, parameter sweeps, data collection, and using the ReCirq framework, see:
+- **[references/experiments.md](references/experiments.md)** - Complete guide to quantum experiments
+
+Common topics:
+- Experiment design patterns
+- Parameter sweeps and data collection
+- ReCirq framework structure
+- Common algorithms (VQE, QAOA, QPE)
+- Data analysis and visualization
+- Statistical analysis and fidelity estimation
+- Parallel data collection
+
+## Common Patterns
+
+### Variational Algorithm Template
+
+```python
+import scipy.optimize
+
+def variational_algorithm(ansatz, cost_function, initial_params):
+    """Template for variational quantum algorithms."""
+
+    def objective(params):
+        circuit = ansatz(params)
+        simulator = cirq.Simulator()
+        result = simulator.simulate(circuit)
+        return cost_function(result)
+
+    # Optimize
+    result = scipy.optimize.minimize(
+        objective,
+        initial_params,
+        method='COBYLA'
+    )
+
+    return result
+
+# Define ansatz
+def my_ansatz(params):
+    q = cirq.LineQubit(0)
+    return cirq.Circuit(
+        cirq.ry(params[0])(q),
+        cirq.rz(params[1])(q)
+    )
+
+# Define cost function
+def my_cost(result):
+    state = result.final_state_vector
+    # Calculate cost based on state
+    return np.real(state[0])
+
+# Run optimization
+result = variational_algorithm(my_ansatz, my_cost, [0.0, 0.0])
+```
+
+### Hardware Execution Template
+
+```python
+def run_on_hardware(circuit, provider='google', device_name='weber', repetitions=1000):
+    """Template for running on quantum hardware."""
+
+    if provider == 'google':
+        import cirq_google
+        engine = cirq_google.get_engine()
+        processor = engine.get_processor(device_name)
+        job = processor.run(circuit, repetitions=repetitions)
+        return job.results()[0]
+
+    elif provider == 'ionq':
+        import cirq_ionq
+        service = cirq_ionq.Service()
+        result = service.run(circuit, repetitions=repetitions, target='qpu')
+        return result
+
+    elif provider == 'azure':
+        from azure.quantum.cirq import AzureQuantumService
+        # Setup workspace...
+        service = AzureQuantumService(workspace)
+        result = service.run(circuit, repetitions=repetitions, target='ionq.qpu')
+        return result
+
+    else:
+        raise ValueError(f"Unknown provider: {provider}")
+```
+
+### Noise Study Template
+
+```python
+def noise_comparison_study(circuit, noise_levels):
+    """Compare circuit performance at different noise levels."""
+
+    results = {}
+
+    for noise_level in noise_levels:
+        # Create noisy circuit
+        noisy_circuit = circuit.with_noise(cirq.depolarize(p=noise_level))
+
+        # Simulate
+        simulator = cirq.DensityMatrixSimulator()
+        result = simulator.run(noisy_circuit, repetitions=1000)
+
+        # Analyze
+        results[noise_level] = {
+            'histogram': result.histogram(key='result'),
+            'dominant_state': max(
+                result.histogram(key='result').items(),
+                key=lambda x: x[1]
+            )
+        }
+
+    return results
+
+# Run study
+noise_levels = [0.0, 0.001, 0.01, 0.05, 0.1]
+results = noise_comparison_study(circuit, noise_levels)
+```
+
+## Best Practices
+
+1. **Circuit Design**
+   - Use appropriate qubit types for your topology
+   - Keep circuits modular and reusable
+   - Label measurements with descriptive keys
+   - Validate circuits against device constraints before execution
+
+2. **Simulation**
+   - Use state vector simulation for pure states (more efficient)
+   - Use density matrix simulation only when needed (mixed states, noise)
+   - Leverage parameter sweeps instead of individual runs
+   - Monitor memory usage for large systems (2^n grows quickly)
+
+3. **Hardware Execution**
+   - Always test on simulators first
+   - Select best qubits using calibration data
+   - Optimize circuits for target hardware gateset
+   - Implement error mitigation for production runs
+   - Store expensive hardware results immediately
+
+4. **Circuit Optimization**
+   - Start with high-level built-in transformers
+   - Chain multiple optimizations in sequence
+   - Track depth and gate count reduction
+   - Validate correctness after transformation
+
+5. **Noise Modeling**
+   - Use realistic noise models from calibration data
+   - Include all error sources (gate, decoherence, readout)
+   - Characterize before mitigating
+   - Keep circuits shallow to minimize noise accumulation
+
+6. **Experiments**
+   - Structure experiments with clear separation (data generation, collection, analysis)
+   - Use ReCirq patterns for reproducibility
+   - Save intermediate results frequently
+   - Parallelize independent tasks
+   - Document thoroughly with metadata
+
+## Additional Resources
+
+- **Official Documentation**: https://quantumai.google/cirq
+- **API Reference**: https://quantumai.google/reference/python/cirq
+- **Tutorials**: https://quantumai.google/cirq/tutorials
+- **Examples**: https://github.com/quantumlib/Cirq/tree/master/examples
+- **ReCirq**: https://github.com/quantumlib/ReCirq
+
+## Common Issues
+
+**Circuit too deep for hardware:**
+- Use circuit optimization transformers to reduce depth
+- See `transformation.md` for optimization techniques
+
+**Memory issues with simulation:**
+- Switch from density matrix to state vector simulator
+- Reduce number of qubits or use stabilizer simulator for Clifford circuits
+
+**Device validation errors:**
+- Check qubit connectivity with device.metadata.nx_graph
+- Decompose gates to device-native gateset
+- See `hardware.md` for device-specific compilation
+
+**Noisy simulation too slow:**
+- Density matrix simulation is O(2^2n) - consider reducing qubits
+- Use noise models selectively on critical operations only
+- See `simulation.md` for performance optimization

scientific-skills/cirq/references/building.md

@@ -0,0 +1,307 @@
+# Building Quantum Circuits
+
+This guide covers circuit construction in Cirq, including qubits, gates, operations, and circuit patterns.
+
+## Basic Circuit Construction
+
+### Creating Circuits
+
+```python
+import cirq
+
+# Create a circuit
+circuit = cirq.Circuit()
+
+# Create qubits
+q0 = cirq.GridQubit(0, 0)
+q1 = cirq.GridQubit(0, 1)
+q2 = cirq.LineQubit(0)
+
+# Add gates to circuit
+circuit.append([
+    cirq.H(q0),
+    cirq.CNOT(q0, q1),
+    cirq.measure(q0, q1, key='result')
+])
+```
+
+### Qubit Types
+
+**GridQubit**: 2D grid topology for hardware-like layouts
+```python
+qubits = cirq.GridQubit.square(2)  # 2x2 grid
+qubit = cirq.GridQubit(row=0, col=1)
+```
+
+**LineQubit**: 1D linear topology
+```python
+qubits = cirq.LineQubit.range(5)  # 5 qubits in a line
+qubit = cirq.LineQubit(3)
+```
+
+**NamedQubit**: Custom-named qubits
+```python
+qubit = cirq.NamedQubit('my_qubit')
+```
+
+## Common Gates and Operations
+
+### Single-Qubit Gates
+
+```python
+# Pauli gates
+cirq.X(qubit)  # NOT gate
+cirq.Y(qubit)
+cirq.Z(qubit)
+
+# Hadamard
+cirq.H(qubit)
+
+# Rotation gates
+cirq.rx(angle)(qubit)  # Rotation around X-axis
+cirq.ry(angle)(qubit)  # Rotation around Y-axis
+cirq.rz(angle)(qubit)  # Rotation around Z-axis
+
+# Phase gates
+cirq.S(qubit)  # √Z gate
+cirq.T(qubit)  # ⁴√Z gate
+```
+
+### Two-Qubit Gates
+
+```python
+# CNOT (Controlled-NOT)
+cirq.CNOT(control, target)
+cirq.CX(control, target)  # Alias
+
+# CZ (Controlled-Z)
+cirq.CZ(q0, q1)
+
+# SWAP
+cirq.SWAP(q0, q1)
+
+# iSWAP
+cirq.ISWAP(q0, q1)
+
+# Controlled rotations
+cirq.CZPowGate(exponent=0.5)(q0, q1)
+```
+
+### Measurement Operations
+
+```python
+# Measure single qubit
+cirq.measure(qubit, key='m')
+
+# Measure multiple qubits
+cirq.measure(q0, q1, q2, key='result')
+
+# Measure all qubits in circuit
+circuit.append(cirq.measure(*qubits, key='final'))
+```
+
+## Advanced Circuit Construction
+
+### Parameterized Gates
+
+```python
+import sympy
+
+# Create symbolic parameters
+theta = sympy.Symbol('theta')
+phi = sympy.Symbol('phi')
+
+# Use in gates
+circuit = cirq.Circuit(
+    cirq.rx(theta)(q0),
+    cirq.ry(phi)(q1),
+    cirq.CNOT(q0, q1)
+)
+
+# Resolve parameters later
+resolved = cirq.resolve_parameters(circuit, {'theta': 0.5, 'phi': 1.2})
+```
+
+### Custom Gates via Unitaries
+
+```python
+import numpy as np
+
+# Define unitary matrix
+unitary = np.array([
+    [1, 0, 0, 0],
+    [0, 1, 0, 0],
+    [0, 0, 0, 1],
+    [0, 0, 1, 0]
+]) / np.sqrt(2)
+
+# Create gate from unitary
+gate = cirq.MatrixGate(unitary)
+operation = gate(q0, q1)
+```
+
+### Gate Decomposition
+
+```python
+# Define custom gate with decomposition
+class MyGate(cirq.Gate):
+    def _num_qubits_(self):
+        return 1
+
+    def _decompose_(self, qubits):
+        q = qubits[0]
+        return [cirq.H(q), cirq.T(q), cirq.H(q)]
+
+    def _circuit_diagram_info_(self, args):
+        return 'MyGate'
+
+# Use the custom gate
+my_gate = MyGate()
+circuit.append(my_gate(q0))
+```
+
+## Circuit Organization
+
+### Moments
+
+Circuits are organized into moments (parallel operations):
+
+```python
+# Explicit moment construction
+circuit = cirq.Circuit(
+    cirq.Moment([cirq.H(q0), cirq.H(q1)]),
+    cirq.Moment([cirq.CNOT(q0, q1)]),
+    cirq.Moment([cirq.measure(q0, key='m0'), cirq.measure(q1, key='m1')])
+)
+
+# Access moments
+for i, moment in enumerate(circuit):
+    print(f"Moment {i}: {moment}")
+```
+
+### Circuit Operations
+
+```python
+# Concatenate circuits
+circuit3 = circuit1 + circuit2
+
+# Insert operations
+circuit.insert(index, operation)
+
+# Append with strategy
+circuit.append(operations, strategy=cirq.InsertStrategy.NEW_THEN_INLINE)
+```
+
+## Circuit Patterns
+
+### Bell State Preparation
+
+```python
+def bell_state_circuit():
+    q0, q1 = cirq.LineQubit.range(2)
+    return cirq.Circuit(
+        cirq.H(q0),
+        cirq.CNOT(q0, q1)
+    )
+```
+
+### GHZ State
+
+```python
+def ghz_circuit(qubits):
+    circuit = cirq.Circuit()
+    circuit.append(cirq.H(qubits[0]))
+    for i in range(len(qubits) - 1):
+        circuit.append(cirq.CNOT(qubits[i], qubits[i+1]))
+    return circuit
+```
+
+### Quantum Fourier Transform
+
+```python
+def qft_circuit(qubits):
+    circuit = cirq.Circuit()
+    for i, q in enumerate(qubits):
+        circuit.append(cirq.H(q))
+        for j in range(i + 1, len(qubits)):
+            circuit.append(cirq.CZPowGate(exponent=1/2**(j-i))(qubits[j], q))
+
+    # Reverse qubit order
+    for i in range(len(qubits) // 2):
+        circuit.append(cirq.SWAP(qubits[i], qubits[len(qubits) - i - 1]))
+
+    return circuit
+```
+
+## Circuit Import/Export
+
+### OpenQASM
+
+```python
+# Export to QASM
+qasm_str = circuit.to_qasm()
+
+# Import from QASM
+from cirq.contrib.qasm_import import circuit_from_qasm
+circuit = circuit_from_qasm(qasm_str)
+```
+
+### Circuit JSON
+
+```python
+import json
+
+# Serialize
+json_str = cirq.to_json(circuit)
+
+# Deserialize
+circuit = cirq.read_json(json_text=json_str)
+```
+
+## Working with Qudits
+
+Qudits are higher-dimensional quantum systems (qutrits, ququarts, etc.):
+
+```python
+# Create qutrit (3-level system)
+qutrit = cirq.LineQid(0, dimension=3)
+
+# Custom qutrit gate
+class QutritXGate(cirq.Gate):
+    def _qid_shape_(self):
+        return (3,)
+
+    def _unitary_(self):
+        return np.array([
+            [0, 0, 1],
+            [1, 0, 0],
+            [0, 1, 0]
+        ])
+
+gate = QutritXGate()
+circuit = cirq.Circuit(gate(qutrit))
+```
+
+## Observables
+
+Create observables from Pauli operators:
+
+```python
+# Single Pauli observable
+obs = cirq.Z(q0)
+
+# Pauli string
+obs = cirq.X(q0) * cirq.Y(q1) * cirq.Z(q2)
+
+# Linear combination
+from cirq import PauliSum
+obs = 0.5 * cirq.X(q0) + 0.3 * cirq.Z(q1)
+```
+
+## Best Practices
+
+1. **Use appropriate qubit types**: GridQubit for hardware-like topologies, LineQubit for 1D problems
+2. **Keep circuits modular**: Build reusable circuit functions
+3. **Use symbolic parameters**: For parameter sweeps and optimization
+4. **Label measurements clearly**: Use descriptive keys for measurement results
+5. **Document custom gates**: Include circuit diagram information for visualization

scientific-skills/cirq/references/experiments.md

@@ -0,0 +1,572 @@
+# Running Quantum Experiments
+
+This guide covers designing and executing quantum experiments, including parameter sweeps, data collection, and using the ReCirq framework.
+
+## Experiment Design
+
+### Basic Experiment Structure
+
+```python
+import cirq
+import numpy as np
+import pandas as pd
+
+class QuantumExperiment:
+    """Base class for quantum experiments."""
+
+    def __init__(self, qubits, simulator=None):
+        self.qubits = qubits
+        self.simulator = simulator or cirq.Simulator()
+        self.results = []
+
+    def build_circuit(self, **params):
+        """Build circuit with given parameters."""
+        raise NotImplementedError
+
+    def run(self, params_list, repetitions=1000):
+        """Run experiment with parameter sweep."""
+        for params in params_list:
+            circuit = self.build_circuit(**params)
+            result = self.simulator.run(circuit, repetitions=repetitions)
+            self.results.append({
+                'params': params,
+                'result': result
+            })
+        return self.results
+
+    def analyze(self):
+        """Analyze experimental results."""
+        raise NotImplementedError
+```
+
+### Parameter Sweeps
+
+```python
+import sympy
+
+# Define parameters
+theta = sympy.Symbol('theta')
+phi = sympy.Symbol('phi')
+
+# Create parameterized circuit
+def parameterized_circuit(qubits, theta, phi):
+    return cirq.Circuit(
+        cirq.ry(theta)(qubits[0]),
+        cirq.rz(phi)(qubits[1]),
+        cirq.CNOT(qubits[0], qubits[1]),
+        cirq.measure(*qubits, key='result')
+    )
+
+# Define sweep
+sweep = cirq.Product(
+    cirq.Linspace('theta', 0, np.pi, 20),
+    cirq.Linspace('phi', 0, 2*np.pi, 20)
+)
+
+# Run sweep
+circuit = parameterized_circuit(cirq.LineQubit.range(2), theta, phi)
+results = cirq.Simulator().run_sweep(circuit, params=sweep, repetitions=1000)
+```
+
+### Data Collection
+
+```python
+def collect_experiment_data(circuit, sweep, simulator, repetitions=1000):
+    """Collect and organize experimental data."""
+
+    data = []
+    results = simulator.run_sweep(circuit, params=sweep, repetitions=repetitions)
+
+    for params, result in zip(sweep, results):
+        # Extract parameters
+        param_dict = {k: v for k, v in params.param_dict.items()}
+
+        # Extract measurements
+        counts = result.histogram(key='result')
+
+        # Store in structured format
+        data.append({
+            **param_dict,
+            'counts': counts,
+            'total': repetitions
+        })
+
+    return pd.DataFrame(data)
+
+# Collect data
+df = collect_experiment_data(circuit, sweep, cirq.Simulator())
+
+# Save to file
+df.to_csv('experiment_results.csv', index=False)
+```
+
+## ReCirq Framework
+
+ReCirq provides a structured framework for reproducible quantum experiments.
+
+### ReCirq Experiment Structure
+
+```python
+"""
+Standard ReCirq experiment structure:
+
+experiment_name/
+├── __init__.py
+├── experiment.py        # Main experiment code
+├── tasks.py            # Data generation tasks
+├── data_collection.py  # Parallel data collection
+├── analysis.py         # Data analysis
+└── plots.py           # Visualization
+"""
+```
+
+### Task-Based Data Collection
+
+```python
+from dataclasses import dataclass
+from typing import List
+import cirq
+
+@dataclass
+class ExperimentTask:
+    """Single task in parameter sweep."""
+    theta: float
+    phi: float
+    repetitions: int = 1000
+
+    def build_circuit(self, qubits):
+        """Build circuit for this task."""
+        return cirq.Circuit(
+            cirq.ry(self.theta)(qubits[0]),
+            cirq.rz(self.phi)(qubits[1]),
+            cirq.CNOT(qubits[0], qubits[1]),
+            cirq.measure(*qubits, key='result')
+        )
+
+    def run(self, qubits, simulator):
+        """Execute task."""
+        circuit = self.build_circuit(qubits)
+        result = simulator.run(circuit, repetitions=self.repetitions)
+        return {
+            'theta': self.theta,
+            'phi': self.phi,
+            'result': result
+        }
+
+# Create tasks
+tasks = [
+    ExperimentTask(theta=t, phi=p)
+    for t in np.linspace(0, np.pi, 10)
+    for p in np.linspace(0, 2*np.pi, 10)
+]
+
+# Execute tasks
+qubits = cirq.LineQubit.range(2)
+simulator = cirq.Simulator()
+results = [task.run(qubits, simulator) for task in tasks]
+```
+
+### Parallel Data Collection
+
+```python
+from multiprocessing import Pool
+import functools
+
+def run_task_parallel(task, qubits, simulator):
+    """Run single task (for parallel execution)."""
+    return task.run(qubits, simulator)
+
+def collect_data_parallel(tasks, qubits, simulator, n_workers=4):
+    """Collect data using parallel processing."""
+
+    # Create partial function with fixed arguments
+    run_func = functools.partial(
+        run_task_parallel,
+        qubits=qubits,
+        simulator=simulator
+    )
+
+    # Run in parallel
+    with Pool(n_workers) as pool:
+        results = pool.map(run_func, tasks)
+
+    return results
+
+# Use parallel collection
+results = collect_data_parallel(tasks, qubits, cirq.Simulator(), n_workers=8)
+```
+
+## Common Quantum Algorithms
+
+### Variational Quantum Eigensolver (VQE)
+
+```python
+import scipy.optimize
+
+def vqe_experiment(hamiltonian, ansatz_func, initial_params):
+    """Run VQE to find ground state energy."""
+
+    def cost_function(params):
+        """Energy expectation value."""
+        circuit = ansatz_func(params)
+
+        # Measure expectation value of Hamiltonian
+        simulator = cirq.Simulator()
+        result = simulator.simulate(circuit)
+        energy = hamiltonian.expectation_from_state_vector(
+            result.final_state_vector,
+            qubit_map={q: i for i, q in enumerate(circuit.all_qubits())}
+        )
+        return energy.real
+
+    # Optimize parameters
+    result = scipy.optimize.minimize(
+        cost_function,
+        initial_params,
+        method='COBYLA'
+    )
+
+    return result
+
+# Example: H2 molecule
+def h2_ansatz(params, qubits):
+    """UCC ansatz for H2."""
+    theta = params[0]
+    return cirq.Circuit(
+        cirq.X(qubits[1]),
+        cirq.ry(theta)(qubits[0]),
+        cirq.CNOT(qubits[0], qubits[1])
+    )
+
+# Define Hamiltonian (simplified)
+qubits = cirq.LineQubit.range(2)
+hamiltonian = cirq.PauliSum.from_pauli_strings([
+    cirq.PauliString({qubits[0]: cirq.Z}),
+    cirq.PauliString({qubits[1]: cirq.Z}),
+    cirq.PauliString({qubits[0]: cirq.Z, qubits[1]: cirq.Z})
+])
+
+# Run VQE
+result = vqe_experiment(
+    hamiltonian,
+    lambda p: h2_ansatz(p, qubits),
+    initial_params=[0.0]
+)
+
+print(f"Ground state energy: {result.fun}")
+print(f"Optimal parameters: {result.x}")
+```
+
+### Quantum Approximate Optimization Algorithm (QAOA)
+
+```python
+def qaoa_circuit(graph, params, p_layers):
+    """QAOA circuit for MaxCut problem."""
+
+    qubits = cirq.LineQubit.range(graph.number_of_nodes())
+    circuit = cirq.Circuit()
+
+    # Initial superposition
+    circuit.append(cirq.H(q) for q in qubits)
+
+    # QAOA layers
+    for layer in range(p_layers):
+        gamma = params[layer]
+        beta = params[p_layers + layer]
+
+        # Problem Hamiltonian (cost)
+        for edge in graph.edges():
+            i, j = edge
+            circuit.append(cirq.ZZPowGate(exponent=gamma)(qubits[i], qubits[j]))
+
+        # Mixer Hamiltonian
+        circuit.append(cirq.rx(2 * beta)(q) for q in qubits)
+
+    circuit.append(cirq.measure(*qubits, key='result'))
+    return circuit
+
+# Run QAOA
+import networkx as nx
+
+graph = nx.cycle_graph(4)
+p_layers = 2
+
+def qaoa_cost(params):
+    """Evaluate QAOA cost function."""
+    circuit = qaoa_circuit(graph, params, p_layers)
+    simulator = cirq.Simulator()
+    result = simulator.run(circuit, repetitions=1000)
+
+    # Calculate MaxCut objective
+    total_cost = 0
+    counts = result.histogram(key='result')
+
+    for bitstring, count in counts.items():
+        cost = 0
+        bits = [(bitstring >> i) & 1 for i in range(graph.number_of_nodes())]
+        for edge in graph.edges():
+            i, j = edge
+            if bits[i] != bits[j]:
+                cost += 1
+        total_cost += cost * count
+
+    return -total_cost / 1000  # Maximize cut
+
+# Optimize
+initial_params = np.random.random(2 * p_layers) * np.pi
+result = scipy.optimize.minimize(qaoa_cost, initial_params, method='COBYLA')
+
+print(f"Optimal cost: {-result.fun}")
+print(f"Optimal parameters: {result.x}")
+```
+
+### Quantum Phase Estimation
+
+```python
+def qpe_circuit(unitary, eigenstate_prep, n_counting_qubits):
+    """Quantum Phase Estimation circuit."""
+
+    counting_qubits = cirq.LineQubit.range(n_counting_qubits)
+    target_qubit = cirq.LineQubit(n_counting_qubits)
+
+    circuit = cirq.Circuit()
+
+    # Prepare eigenstate
+    circuit.append(eigenstate_prep(target_qubit))
+
+    # Apply Hadamard to counting qubits
+    circuit.append(cirq.H(q) for q in counting_qubits)
+
+    # Controlled unitaries
+    for i, q in enumerate(counting_qubits):
+        power = 2 ** (n_counting_qubits - 1 - i)
+        # Apply controlled-U^power
+        for _ in range(power):
+            circuit.append(cirq.ControlledGate(unitary)(q, target_qubit))
+
+    # Inverse QFT on counting qubits
+    circuit.append(inverse_qft(counting_qubits))
+
+    # Measure counting qubits
+    circuit.append(cirq.measure(*counting_qubits, key='phase'))
+
+    return circuit
+
+def inverse_qft(qubits):
+    """Inverse Quantum Fourier Transform."""
+    n = len(qubits)
+    ops = []
+
+    for i in range(n // 2):
+        ops.append(cirq.SWAP(qubits[i], qubits[n - i - 1]))
+
+    for i in range(n):
+        for j in range(i):
+            ops.append(cirq.CZPowGate(exponent=-1/2**(i-j))(qubits[j], qubits[i]))
+        ops.append(cirq.H(qubits[i]))
+
+    return ops
+```
+
+## Data Analysis
+
+### Statistical Analysis
+
+```python
+def analyze_measurement_statistics(results):
+    """Analyze measurement statistics."""
+
+    counts = results.histogram(key='result')
+    total = sum(counts.values())
+
+    # Calculate probabilities
+    probabilities = {state: count/total for state, count in counts.items()}
+
+    # Shannon entropy
+    entropy = -sum(p * np.log2(p) for p in probabilities.values() if p > 0)
+
+    # Most likely outcome
+    most_likely = max(counts.items(), key=lambda x: x[1])
+
+    return {
+        'probabilities': probabilities,
+        'entropy': entropy,
+        'most_likely_state': most_likely[0],
+        'most_likely_probability': most_likely[1] / total
+    }
+```
+
+### Expectation Value Calculation
+
+```python
+def calculate_expectation_value(circuit, observable, simulator):
+    """Calculate expectation value of observable."""
+
+    # Remove measurements
+    circuit_no_measure = cirq.Circuit(
+        m for m in circuit if not isinstance(m, cirq.MeasurementGate)
+    )
+
+    result = simulator.simulate(circuit_no_measure)
+    state_vector = result.final_state_vector
+
+    # Calculate ⟨ψ|O|ψ⟩
+    expectation = observable.expectation_from_state_vector(
+        state_vector,
+        qubit_map={q: i for i, q in enumerate(circuit.all_qubits())}
+    )
+
+    return expectation.real
+```
+
+### Fidelity Estimation
+
+```python
+def state_fidelity(state1, state2):
+    """Calculate fidelity between two states."""
+    return np.abs(np.vdot(state1, state2)) ** 2
+
+def process_fidelity(result1, result2):
+    """Calculate process fidelity from measurement results."""
+
+    counts1 = result1.histogram(key='result')
+    counts2 = result2.histogram(key='result')
+
+    # Normalize to probabilities
+    total1 = sum(counts1.values())
+    total2 = sum(counts2.values())
+
+    probs1 = {k: v/total1 for k, v in counts1.items()}
+    probs2 = {k: v/total2 for k, v in counts2.items()}
+
+    # Classical fidelity (Bhattacharyya coefficient)
+    all_states = set(probs1.keys()) | set(probs2.keys())
+    fidelity = sum(np.sqrt(probs1.get(s, 0) * probs2.get(s, 0))
+                   for s in all_states) ** 2
+
+    return fidelity
+```
+
+## Visualization
+
+### Plot Parameter Landscapes
+
+```python
+import matplotlib.pyplot as plt
+
+def plot_parameter_landscape(theta_vals, phi_vals, energies):
+    """Plot 2D parameter landscape."""
+
+    plt.figure(figsize=(10, 8))
+    plt.contourf(theta_vals, phi_vals, energies, levels=50, cmap='viridis')
+    plt.colorbar(label='Energy')
+    plt.xlabel('θ')
+    plt.ylabel('φ')
+    plt.title('Energy Landscape')
+    plt.show()
+```
+
+### Plot Convergence
+
+```python
+def plot_optimization_convergence(optimization_history):
+    """Plot optimization convergence."""
+
+    iterations = range(len(optimization_history))
+    energies = [result['energy'] for result in optimization_history]
+
+    plt.figure(figsize=(10, 6))
+    plt.plot(iterations, energies, 'b-', linewidth=2)
+    plt.xlabel('Iteration')
+    plt.ylabel('Energy')
+    plt.title('Optimization Convergence')
+    plt.grid(True)
+    plt.show()
+```
+
+### Plot Measurement Distributions
+
+```python
+def plot_measurement_distribution(results):
+    """Plot measurement outcome distribution."""
+
+    counts = results.histogram(key='result')
+
+    plt.figure(figsize=(12, 6))
+    plt.bar(counts.keys(), counts.values())
+    plt.xlabel('Measurement Outcome')
+    plt.ylabel('Counts')
+    plt.title('Measurement Distribution')
+    plt.xticks(rotation=45)
+    plt.tight_layout()
+    plt.show()
+```
+
+## Best Practices
+
+1. **Structure experiments clearly**: Use ReCirq patterns for reproducibility
+2. **Separate tasks**: Divide data generation, collection, and analysis
+3. **Use parameter sweeps**: Explore parameter space systematically
+4. **Save intermediate results**: Don't lose expensive computation
+5. **Parallelize when possible**: Use multiprocessing for independent tasks
+6. **Track metadata**: Record experiment conditions, timestamps, versions
+7. **Validate on simulators**: Test experimental code before hardware
+8. **Implement error handling**: Robust code for long-running experiments
+9. **Version control data**: Track experimental data alongside code
+10. **Document thoroughly**: Clear documentation for reproducibility
+
+## Example: Complete Experiment
+
+```python
+# Full experimental workflow
+class VQEExperiment(QuantumExperiment):
+    """Complete VQE experiment."""
+
+    def __init__(self, hamiltonian, ansatz, qubits):
+        super().__init__(qubits)
+        self.hamiltonian = hamiltonian
+        self.ansatz = ansatz
+        self.history = []
+
+    def build_circuit(self, params):
+        return self.ansatz(params, self.qubits)
+
+    def cost_function(self, params):
+        circuit = self.build_circuit(params)
+        result = self.simulator.simulate(circuit)
+        energy = self.hamiltonian.expectation_from_state_vector(
+            result.final_state_vector,
+            qubit_map={q: i for i, q in enumerate(self.qubits)}
+        )
+        self.history.append({'params': params, 'energy': energy.real})
+        return energy.real
+
+    def run(self, initial_params):
+        result = scipy.optimize.minimize(
+            self.cost_function,
+            initial_params,
+            method='COBYLA',
+            options={'maxiter': 100}
+        )
+        return result
+
+    def analyze(self):
+        # Plot convergence
+        energies = [h['energy'] for h in self.history]
+        plt.plot(energies)
+        plt.xlabel('Iteration')
+        plt.ylabel('Energy')
+        plt.title('VQE Convergence')
+        plt.show()
+
+        return {
+            'final_energy': self.history[-1]['energy'],
+            'optimal_params': self.history[-1]['params'],
+            'num_iterations': len(self.history)
+        }
+
+# Run experiment
+experiment = VQEExperiment(hamiltonian, h2_ansatz, qubits)
+result = experiment.run(initial_params=[0.0])
+analysis = experiment.analyze()
+```

scientific-skills/cirq/references/hardware.md

@@ -0,0 +1,515 @@
+# Hardware Integration
+
+This guide covers running quantum circuits on real quantum hardware through Cirq's device interfaces and service providers.
+
+## Device Representation
+
+### Device Classes
+
+```python
+import cirq
+
+# Define device with connectivity
+class MyDevice(cirq.Device):
+    def __init__(self, qubits, connectivity):
+        self.qubits = qubits
+        self.connectivity = connectivity
+
+    @property
+    def metadata(self):
+        return cirq.DeviceMetadata(
+            self.qubits,
+            self.connectivity
+        )
+
+    def validate_operation(self, operation):
+        # Check if operation is valid on this device
+        if len(operation.qubits) == 2:
+            q0, q1 = operation.qubits
+            if (q0, q1) not in self.connectivity:
+                raise ValueError(f"Qubits {q0} and {q1} not connected")
+```
+
+### Device Constraints
+
+```python
+# Check device metadata
+device = cirq_google.Sycamore
+
+# Get qubit topology
+qubits = device.metadata.qubit_set
+print(f"Available qubits: {len(qubits)}")
+
+# Check connectivity
+for q0 in qubits:
+    neighbors = device.metadata.nx_graph.neighbors(q0)
+    print(f"{q0} connected to: {list(neighbors)}")
+
+# Validate circuit against device
+try:
+    device.validate_circuit(circuit)
+    print("Circuit is valid for device")
+except ValueError as e:
+    print(f"Invalid circuit: {e}")
+```
+
+## Qubit Selection
+
+### Best Qubit Selection
+
+```python
+import cirq_google
+
+# Get calibration metrics
+processor = cirq_google.get_engine().get_processor('weber')
+calibration = processor.get_current_calibration()
+
+# Find qubits with lowest error rates
+def select_best_qubits(calibration, n_qubits):
+    """Select n qubits with best single-qubit gate fidelity."""
+    qubit_fidelities = {}
+
+    for qubit in calibration.keys():
+        if 'single_qubit_rb_average_error_per_gate' in calibration[qubit]:
+            error = calibration[qubit]['single_qubit_rb_average_error_per_gate']
+            qubit_fidelities[qubit] = 1 - error
+
+    # Sort by fidelity
+    best_qubits = sorted(
+        qubit_fidelities.items(),
+        key=lambda x: x[1],
+        reverse=True
+    )[:n_qubits]
+
+    return [q for q, _ in best_qubits]
+
+best_qubits = select_best_qubits(calibration, n_qubits=10)
+```
+
+### Topology-Aware Selection
+
+```python
+def select_connected_qubits(device, n_qubits):
+    """Select connected qubits forming a path or grid."""
+    graph = device.metadata.nx_graph
+
+    # Find connected subgraph
+    import networkx as nx
+    for node in graph.nodes():
+        subgraph = nx.ego_graph(graph, node, radius=n_qubits)
+        if len(subgraph) >= n_qubits:
+            return list(subgraph.nodes())[:n_qubits]
+
+    raise ValueError(f"Could not find {n_qubits} connected qubits")
+```
+
+## Service Providers
+
+### Google Quantum AI (Cirq-Google)
+
+#### Setup
+
+```python
+import cirq_google
+
+# Authenticate (requires Google Cloud project)
+# Set environment variable: GOOGLE_CLOUD_PROJECT=your-project-id
+
+# Get quantum engine
+engine = cirq_google.get_engine()
+
+# List available processors
+processors = engine.list_processors()
+for processor in processors:
+    print(f"Processor: {processor.processor_id}")
+```
+
+#### Running on Google Hardware
+
+```python
+# Create circuit for Google device
+import cirq_google
+
+# Get processor
+processor = engine.get_processor('weber')
+device = processor.get_device()
+
+# Create circuit on device qubits
+qubits = sorted(device.metadata.qubit_set)[:5]
+circuit = cirq.Circuit(
+    cirq.H(qubits[0]),
+    cirq.CZ(qubits[0], qubits[1]),
+    cirq.measure(*qubits, key='result')
+)
+
+# Validate and run
+device.validate_circuit(circuit)
+job = processor.run(circuit, repetitions=1000)
+
+# Get results
+results = job.results()[0]
+print(results.histogram(key='result'))
+```
+
+### IonQ
+
+#### Setup
+
+```python
+import cirq_ionq
+
+# Set API key
+# Option 1: Environment variable
+# export IONQ_API_KEY=your_api_key
+
+# Option 2: In code
+service = cirq_ionq.Service(api_key='your_api_key')
+```
+
+#### Running on IonQ
+
+```python
+import cirq_ionq
+
+# Create service
+service = cirq_ionq.Service(api_key='your_api_key')
+
+# Create circuit (IonQ uses generic qubits)
+qubits = cirq.LineQubit.range(3)
+circuit = cirq.Circuit(
+    cirq.H(qubits[0]),
+    cirq.CNOT(qubits[0], qubits[1]),
+    cirq.CNOT(qubits[1], qubits[2]),
+    cirq.measure(*qubits, key='result')
+)
+
+# Run on simulator
+result = service.run(
+    circuit=circuit,
+    repetitions=1000,
+    target='simulator'
+)
+print(result.histogram(key='result'))
+
+# Run on hardware
+result = service.run(
+    circuit=circuit,
+    repetitions=1000,
+    target='qpu'
+)
+```
+
+#### IonQ Job Management
+
+```python
+# Create job
+job = service.create_job(circuit, repetitions=1000, target='qpu')
+
+# Check job status
+status = job.status()
+print(f"Job status: {status}")
+
+# Wait for completion
+job.wait_until_complete()
+
+# Get results
+results = job.results()
+```
+
+#### IonQ Calibration Data
+
+```python
+# Get current calibration
+calibration = service.get_current_calibration()
+
+# Access metrics
+print(f"Fidelity: {calibration['fidelity']}")
+print(f"Timing: {calibration['timing']}")
+```
+
+### Azure Quantum
+
+#### Setup
+
+```python
+from azure.quantum import Workspace
+from azure.quantum.cirq import AzureQuantumService
+
+# Create workspace connection
+workspace = Workspace(
+    resource_id="/subscriptions/.../resourceGroups/.../providers/Microsoft.Quantum/Workspaces/...",
+    location="eastus"
+)
+
+# Create Cirq service
+service = AzureQuantumService(workspace)
+```
+
+#### Running on Azure Quantum (IonQ Backend)
+
+```python
+# List available targets
+targets = service.targets()
+for target in targets:
+    print(f"Target: {target.name}")
+
+# Run on IonQ simulator
+result = service.run(
+    circuit=circuit,
+    repetitions=1000,
+    target='ionq.simulator'
+)
+
+# Run on IonQ QPU
+result = service.run(
+    circuit=circuit,
+    repetitions=1000,
+    target='ionq.qpu'
+)
+```
+
+#### Running on Azure Quantum (Honeywell Backend)
+
+```python
+# Run on Honeywell System Model H1
+result = service.run(
+    circuit=circuit,
+    repetitions=1000,
+    target='honeywell.hqs-lt-s1'
+)
+
+# Check Honeywell-specific options
+target_info = service.get_target('honeywell.hqs-lt-s1')
+print(f"Target info: {target_info}")
+```
+
+### AQT (Alpine Quantum Technologies)
+
+#### Setup
+
+```python
+import cirq_aqt
+
+# Set API token
+# export AQT_TOKEN=your_token
+
+# Create service
+service = cirq_aqt.AQTSampler(
+    remote_host='https://gateway.aqt.eu',
+    access_token='your_token'
+)
+```
+
+#### Running on AQT
+
+```python
+# Create circuit
+qubits = cirq.LineQubit.range(3)
+circuit = cirq.Circuit(
+    cirq.H(qubits[0]),
+    cirq.CNOT(qubits[0], qubits[1]),
+    cirq.measure(*qubits, key='result')
+)
+
+# Run on simulator
+result = service.run(
+    circuit,
+    repetitions=1000,
+    target='simulator'
+)
+
+# Run on device
+result = service.run(
+    circuit,
+    repetitions=1000,
+    target='device'
+)
+```
+
+### Pasqal
+
+#### Setup
+
+```python
+import cirq_pasqal
+
+# Create Pasqal device
+device = cirq_pasqal.PasqalDevice(qubits=cirq.LineQubit.range(10))
+```
+
+#### Running on Pasqal
+
+```python
+# Create sampler
+sampler = cirq_pasqal.PasqalSampler(
+    remote_host='https://api.pasqal.cloud',
+    access_token='your_token',
+    device=device
+)
+
+# Run circuit
+result = sampler.run(circuit, repetitions=1000)
+```
+
+## Hardware Best Practices
+
+### Circuit Optimization for Hardware
+
+```python
+def optimize_for_hardware(circuit, device):
+    """Optimize circuit for specific hardware."""
+    from cirq.transformers import (
+        optimize_for_target_gateset,
+        merge_single_qubit_gates_to_phxz,
+        drop_negligible_operations
+    )
+
+    # Get device gateset
+    if hasattr(device, 'gateset'):
+        gateset = device.gateset
+    else:
+        gateset = cirq.CZTargetGateset()  # Default
+
+    # Optimize
+    circuit = merge_single_qubit_gates_to_phxz(circuit)
+    circuit = drop_negligible_operations(circuit)
+    circuit = optimize_for_target_gateset(circuit, gateset=gateset)
+
+    return circuit
+```
+
+### Error Mitigation
+
+```python
+def run_with_readout_error_mitigation(circuit, sampler, repetitions):
+    """Mitigate readout errors using calibration."""
+
+    # Measure readout error
+    cal_circuits = []
+    for state in range(2**len(circuit.qubits)):
+        cal_circuit = cirq.Circuit()
+        for i, q in enumerate(circuit.qubits):
+            if state & (1 << i):
+                cal_circuit.append(cirq.X(q))
+        cal_circuit.append(cirq.measure(*circuit.qubits, key='m'))
+        cal_circuits.append(cal_circuit)
+
+    # Run calibration
+    cal_results = [sampler.run(c, repetitions=1000) for c in cal_circuits]
+
+    # Build confusion matrix
+    # ... (implementation details)
+
+    # Run actual circuit
+    result = sampler.run(circuit, repetitions=repetitions)
+
+    # Apply correction
+    # ... (apply inverse of confusion matrix)
+
+    return result
+```
+
+### Job Management
+
+```python
+def submit_jobs_in_batches(circuits, sampler, batch_size=10):
+    """Submit multiple circuits in batches."""
+    jobs = []
+
+    for i in range(0, len(circuits), batch_size):
+        batch = circuits[i:i+batch_size]
+        job_ids = []
+
+        for circuit in batch:
+            job = sampler.run_async(circuit, repetitions=1000)
+            job_ids.append(job)
+
+        jobs.extend(job_ids)
+
+    # Wait for all jobs
+    results = [job.result() for job in jobs]
+    return results
+```
+
+## Device Specifications
+
+### Checking Device Capabilities
+
+```python
+def print_device_info(device):
+    """Print device capabilities and constraints."""
+
+    print(f"Device: {device}")
+    print(f"Number of qubits: {len(device.metadata.qubit_set)}")
+
+    # Gate support
+    print("\nSupported gates:")
+    if hasattr(device, 'gateset'):
+        for gate in device.gateset.gates:
+            print(f"  - {gate}")
+
+    # Connectivity
+    print("\nConnectivity:")
+    graph = device.metadata.nx_graph
+    print(f"  Edges: {graph.number_of_edges()}")
+    print(f"  Average degree: {sum(dict(graph.degree()).values()) / graph.number_of_nodes():.2f}")
+
+    # Duration constraints
+    if hasattr(device, 'gate_durations'):
+        print("\nGate durations:")
+        for gate, duration in device.gate_durations.items():
+            print(f"  {gate}: {duration}")
+```
+
+## Authentication and Access
+
+### Setting Up Credentials
+
+**Google Cloud:**
+```bash
+# Install gcloud CLI
+# Visit: https://cloud.google.com/sdk/docs/install
+
+# Authenticate
+gcloud auth application-default login
+
+# Set project
+export GOOGLE_CLOUD_PROJECT=your-project-id
+```
+
+**IonQ:**
+```bash
+# Set API key
+export IONQ_API_KEY=your_api_key
+```
+
+**Azure Quantum:**
+```python
+# Use Azure CLI or workspace connection string
+# See: https://docs.microsoft.com/azure/quantum/
+```
+
+**AQT:**
+```bash
+# Request access token from AQT
+export AQT_TOKEN=your_token
+```
+
+**Pasqal:**
+```bash
+# Request API access from Pasqal
+export PASQAL_TOKEN=your_token
+```
+
+## Best Practices
+
+1. **Validate circuits before submission**: Use device.validate_circuit()
+2. **Optimize for target hardware**: Decompose to native gates
+3. **Select best qubits**: Use calibration data for qubit selection
+4. **Monitor job status**: Check job completion before retrieving results
+5. **Implement error mitigation**: Use readout error correction
+6. **Batch jobs efficiently**: Submit multiple circuits together
+7. **Respect rate limits**: Follow provider-specific API limits
+8. **Store results**: Save expensive hardware results immediately
+9. **Test on simulators first**: Validate on simulators before hardware
+10. **Keep circuits shallow**: Hardware has limited coherence times

scientific-skills/cirq/references/noise.md

@@ -0,0 +1,515 @@
+# Noise Modeling and Mitigation
+
+This guide covers noise models, noisy simulation, characterization, and error mitigation in Cirq.
+
+## Noise Channels
+
+### Depolarizing Noise
+
+```python
+import cirq
+import numpy as np
+
+# Single-qubit depolarizing channel
+depol_channel = cirq.depolarize(p=0.01)
+
+# Apply to qubit
+q = cirq.LineQubit(0)
+noisy_op = depol_channel(q)
+
+# Add to circuit
+circuit = cirq.Circuit(
+    cirq.H(q),
+    depol_channel(q),
+    cirq.measure(q, key='m')
+)
+```
+
+### Amplitude Damping
+
+```python
+# Amplitude damping (T1 decay)
+gamma = 0.1
+amp_damp = cirq.amplitude_damp(gamma)
+
+# Apply after gate
+circuit = cirq.Circuit(
+    cirq.X(q),
+    amp_damp(q)
+)
+```
+
+### Phase Damping
+
+```python
+# Phase damping (T2 dephasing)
+gamma = 0.1
+phase_damp = cirq.phase_damp(gamma)
+
+circuit = cirq.Circuit(
+    cirq.H(q),
+    phase_damp(q)
+)
+```
+
+### Bit Flip Noise
+
+```python
+# Bit flip channel
+bit_flip_prob = 0.01
+bit_flip = cirq.bit_flip(bit_flip_prob)
+
+circuit = cirq.Circuit(
+    cirq.H(q),
+    bit_flip(q)
+)
+```
+
+### Phase Flip Noise
+
+```python
+# Phase flip channel
+phase_flip_prob = 0.01
+phase_flip = cirq.phase_flip(phase_flip_prob)
+
+circuit = cirq.Circuit(
+    cirq.H(q),
+    phase_flip(q)
+)
+```
+
+### Generalized Amplitude Damping
+
+```python
+# Generalized amplitude damping
+p = 0.1  # Damping probability
+gamma = 0.2  # Excitation probability
+gen_amp_damp = cirq.generalized_amplitude_damp(p=p, gamma=gamma)
+```
+
+### Reset Channel
+
+```python
+# Reset to |0⟩ or |1⟩
+reset_to_zero = cirq.reset(q)
+
+# Reset appears as measurement followed by conditional flip
+circuit = cirq.Circuit(
+    cirq.H(q),
+    reset_to_zero
+)
+```
+
+## Noise Models
+
+### Constant Noise Model
+
+```python
+# Apply same noise to all qubits
+noise = cirq.ConstantQubitNoiseModel(
+    qubit_noise_gate=cirq.depolarize(0.01)
+)
+
+# Simulate with noise
+simulator = cirq.DensityMatrixSimulator(noise=noise)
+result = simulator.run(circuit, repetitions=1000)
+```
+
+### Gate-Specific Noise
+
+```python
+class CustomNoiseModel(cirq.NoiseModel):
+    """Apply different noise to different gate types."""
+
+    def noisy_operation(self, op):
+        # Single-qubit gates: depolarizing noise
+        if len(op.qubits) == 1:
+            return [op, cirq.depolarize(0.001)(op.qubits[0])]
+
+        # Two-qubit gates: higher depolarizing noise
+        elif len(op.qubits) == 2:
+            return [
+                op,
+                cirq.depolarize(0.01)(op.qubits[0]),
+                cirq.depolarize(0.01)(op.qubits[1])
+            ]
+
+        return op
+
+# Use custom noise model
+noise_model = CustomNoiseModel()
+simulator = cirq.DensityMatrixSimulator(noise=noise_model)
+```
+
+### Qubit-Specific Noise
+
+```python
+class QubitSpecificNoise(cirq.NoiseModel):
+    """Different noise for different qubits."""
+
+    def __init__(self, qubit_noise_map):
+        self.qubit_noise_map = qubit_noise_map
+
+    def noisy_operation(self, op):
+        noise_ops = [op]
+        for qubit in op.qubits:
+            if qubit in self.qubit_noise_map:
+                noise = self.qubit_noise_map[qubit]
+                noise_ops.append(noise(qubit))
+        return noise_ops
+
+# Define per-qubit noise
+q0, q1, q2 = cirq.LineQubit.range(3)
+noise_map = {
+    q0: cirq.depolarize(0.001),
+    q1: cirq.depolarize(0.005),
+    q2: cirq.depolarize(0.002)
+}
+
+noise_model = QubitSpecificNoise(noise_map)
+```
+
+### Thermal Noise
+
+```python
+class ThermalNoise(cirq.NoiseModel):
+    """Thermal relaxation noise."""
+
+    def __init__(self, T1, T2, gate_time):
+        self.T1 = T1  # Amplitude damping time
+        self.T2 = T2  # Dephasing time
+        self.gate_time = gate_time
+
+    def noisy_operation(self, op):
+        # Calculate probabilities
+        p_amp = 1 - np.exp(-self.gate_time / self.T1)
+        p_phase = 1 - np.exp(-self.gate_time / self.T2)
+
+        noise_ops = [op]
+        for qubit in op.qubits:
+            noise_ops.append(cirq.amplitude_damp(p_amp)(qubit))
+            noise_ops.append(cirq.phase_damp(p_phase)(qubit))
+
+        return noise_ops
+
+# Typical superconducting qubit parameters
+T1 = 50e-6  # 50 μs
+T2 = 30e-6  # 30 μs
+gate_time = 25e-9  # 25 ns
+
+noise_model = ThermalNoise(T1, T2, gate_time)
+```
+
+## Adding Noise to Circuits
+
+### with_noise Method
+
+```python
+# Add noise to all operations
+noisy_circuit = circuit.with_noise(cirq.depolarize(p=0.01))
+
+# Simulate noisy circuit
+simulator = cirq.DensityMatrixSimulator()
+result = simulator.run(noisy_circuit, repetitions=1000)
+```
+
+### insert_into_circuit Method
+
+```python
+# Manual noise insertion
+def add_noise_to_circuit(circuit, noise_model):
+    noisy_moments = []
+    for moment in circuit:
+        ops = []
+        for op in moment:
+            ops.extend(noise_model.noisy_operation(op))
+        noisy_moments.append(cirq.Moment(ops))
+    return cirq.Circuit(noisy_moments)
+```
+
+## Readout Noise
+
+### Measurement Error Model
+
+```python
+class ReadoutNoiseModel(cirq.NoiseModel):
+    """Model readout/measurement errors."""
+
+    def __init__(self, p0_given_1, p1_given_0):
+        # p0_given_1: Probability of measuring 0 when state is 1
+        # p1_given_0: Probability of measuring 1 when state is 0
+        self.p0_given_1 = p0_given_1
+        self.p1_given_0 = p1_given_0
+
+    def noisy_operation(self, op):
+        if isinstance(op.gate, cirq.MeasurementGate):
+            # Apply bit flip before measurement
+            noise_ops = []
+            for qubit in op.qubits:
+                # Average readout error
+                p_error = (self.p0_given_1 + self.p1_given_0) / 2
+                noise_ops.append(cirq.bit_flip(p_error)(qubit))
+            noise_ops.append(op)
+            return noise_ops
+        return op
+
+# Typical readout errors
+readout_noise = ReadoutNoiseModel(p0_given_1=0.02, p1_given_0=0.01)
+```
+
+## Noise Characterization
+
+### Randomized Benchmarking
+
+```python
+import cirq
+
+def generate_rb_circuit(qubits, depth):
+    """Generate randomized benchmarking circuit."""
+    # Random Clifford gates
+    clifford_gates = [cirq.X, cirq.Y, cirq.Z, cirq.H, cirq.S]
+
+    circuit = cirq.Circuit()
+    for _ in range(depth):
+        for qubit in qubits:
+            gate = np.random.choice(clifford_gates)
+            circuit.append(gate(qubit))
+
+    # Add inverse to return to initial state (ideally)
+    # (simplified - proper RB requires tracking full sequence)
+
+    circuit.append(cirq.measure(*qubits, key='result'))
+    return circuit
+
+# Run RB experiment
+def run_rb_experiment(qubits, depths, repetitions=1000):
+    """Run randomized benchmarking at various depths."""
+    simulator = cirq.DensityMatrixSimulator(
+        noise=cirq.ConstantQubitNoiseModel(cirq.depolarize(0.01))
+    )
+
+    survival_probs = []
+    for depth in depths:
+        circuits = [generate_rb_circuit(qubits, depth) for _ in range(20)]
+
+        total_survival = 0
+        for circuit in circuits:
+            result = simulator.run(circuit, repetitions=repetitions)
+            # Calculate survival probability (returned to |0⟩)
+            counts = result.histogram(key='result')
+            survival = counts.get(0, 0) / repetitions
+            total_survival += survival
+
+        avg_survival = total_survival / len(circuits)
+        survival_probs.append(avg_survival)
+
+    return survival_probs
+
+# Fit to extract error rate
+# p_survival = A * p^depth + B
+# Error per gate ≈ (1 - p) / 2
+```
+
+### Cross-Entropy Benchmarking (XEB)
+
+```python
+def xeb_fidelity(circuit, simulator, ideal_probs, repetitions=10000):
+    """Calculate XEB fidelity."""
+
+    # Run noisy simulation
+    result = simulator.run(circuit, repetitions=repetitions)
+    measured_probs = result.histogram(key='result')
+
+    # Normalize
+    for key in measured_probs:
+        measured_probs[key] /= repetitions
+
+    # Calculate cross-entropy
+    cross_entropy = 0
+    for bitstring, prob in measured_probs.items():
+        if bitstring in ideal_probs:
+            cross_entropy += prob * np.log2(ideal_probs[bitstring])
+
+    # Convert to fidelity
+    n_qubits = len(circuit.all_qubits())
+    fidelity = (2**n_qubits * cross_entropy + 1) / (2**n_qubits - 1)
+
+    return fidelity
+```
+
+## Noise Visualization
+
+### Heatmap Visualization
+
+```python
+import matplotlib.pyplot as plt
+
+def plot_noise_heatmap(device, noise_metric):
+    """Plot noise characteristics across 2D grid device."""
+
+    # Get device qubits (assuming GridQubit)
+    qubits = sorted(device.metadata.qubit_set)
+    rows = max(q.row for q in qubits) + 1
+    cols = max(q.col for q in qubits) + 1
+
+    # Create heatmap data
+    heatmap = np.full((rows, cols), np.nan)
+
+    for qubit in qubits:
+        if isinstance(qubit, cirq.GridQubit):
+            value = noise_metric.get(qubit, 0)
+            heatmap[qubit.row, qubit.col] = value
+
+    # Plot
+    plt.figure(figsize=(10, 8))
+    plt.imshow(heatmap, cmap='RdYlGn_r', interpolation='nearest')
+    plt.colorbar(label='Error Rate')
+    plt.title('Qubit Error Rates')
+    plt.xlabel('Column')
+    plt.ylabel('Row')
+    plt.show()
+
+# Example usage
+noise_metric = {q: np.random.random() * 0.01 for q in device.metadata.qubit_set}
+plot_noise_heatmap(device, noise_metric)
+```
+
+### Gate Fidelity Visualization
+
+```python
+def plot_gate_fidelities(calibration_data):
+    """Plot single- and two-qubit gate fidelities."""
+
+    sq_fidelities = []
+    tq_fidelities = []
+
+    for qubit, metrics in calibration_data.items():
+        if 'single_qubit_rb_fidelity' in metrics:
+            sq_fidelities.append(metrics['single_qubit_rb_fidelity'])
+        if 'two_qubit_rb_fidelity' in metrics:
+            tq_fidelities.append(metrics['two_qubit_rb_fidelity'])
+
+    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
+
+    ax1.hist(sq_fidelities, bins=20)
+    ax1.set_xlabel('Single-Qubit Gate Fidelity')
+    ax1.set_ylabel('Count')
+    ax1.set_title('Single-Qubit Gate Fidelities')
+
+    ax2.hist(tq_fidelities, bins=20)
+    ax2.set_xlabel('Two-Qubit Gate Fidelity')
+    ax2.set_ylabel('Count')
+    ax2.set_title('Two-Qubit Gate Fidelities')
+
+    plt.tight_layout()
+    plt.show()
+```
+
+## Error Mitigation Techniques
+
+### Zero-Noise Extrapolation
+
+```python
+def zero_noise_extrapolation(circuit, noise_levels, simulator):
+    """Extrapolate to zero noise limit."""
+
+    expectation_values = []
+
+    for noise_level in noise_levels:
+        # Scale noise
+        noisy_circuit = circuit.with_noise(
+            cirq.depolarize(p=noise_level)
+        )
+
+        # Measure expectation
+        result = simulator.simulate(noisy_circuit)
+        # ... calculate expectation value
+        exp_val = calculate_expectation(result)
+        expectation_values.append(exp_val)
+
+    # Extrapolate to zero noise
+    from scipy.optimize import curve_fit
+
+    def exponential_fit(x, a, b, c):
+        return a * np.exp(-b * x) + c
+
+    popt, _ = curve_fit(exponential_fit, noise_levels, expectation_values)
+    zero_noise_value = popt[2]
+
+    return zero_noise_value
+```
+
+### Probabilistic Error Cancellation
+
+```python
+def quasi_probability_decomposition(noisy_gate, ideal_gate, noise_model):
+    """Decompose noisy gate into quasi-probability distribution."""
+
+    # Decompose noisy gate as: N = ideal + error
+    # Invert: ideal = (N - error) / (1 - error_rate)
+
+    # This creates a quasi-probability distribution
+    # (some probabilities may be negative)
+
+    # Implementation depends on specific noise model
+    pass
+```
+
+### Readout Error Mitigation
+
+```python
+def mitigate_readout_errors(results, confusion_matrix):
+    """Apply readout error mitigation using confusion matrix."""
+
+    # Invert confusion matrix
+    inv_confusion = np.linalg.inv(confusion_matrix)
+
+    # Get measured counts
+    counts = results.histogram(key='result')
+
+    # Convert to probability vector
+    total_counts = sum(counts.values())
+    measured_probs = np.array([counts.get(i, 0) / total_counts
+                               for i in range(len(confusion_matrix))])
+
+    # Apply inverse
+    corrected_probs = inv_confusion @ measured_probs
+
+    # Convert back to counts
+    corrected_counts = {i: int(p * total_counts)
+                       for i, p in enumerate(corrected_probs) if p > 0}
+
+    return corrected_counts
+```
+
+## Hardware-Based Noise Models
+
+### From Google Calibration
+
+```python
+import cirq_google
+
+# Get calibration data
+processor = cirq_google.get_engine().get_processor('weber')
+noise_props = processor.get_device_specification()
+
+# Create noise model from calibration
+noise_model = cirq_google.NoiseModelFromGoogleNoiseProperties(noise_props)
+
+# Simulate with realistic noise
+simulator = cirq.DensityMatrixSimulator(noise=noise_model)
+result = simulator.run(circuit, repetitions=1000)
+```
+
+## Best Practices
+
+1. **Use density matrix simulator for noisy simulations**: State vector simulators cannot model mixed states
+2. **Match noise model to hardware**: Use calibration data when available
+3. **Include all error sources**: Gate errors, decoherence, readout errors
+4. **Characterize before mitigating**: Understand noise before applying mitigation
+5. **Consider error propagation**: Noise compounds through circuit depth
+6. **Use appropriate benchmarking**: RB for gate errors, XEB for full circuit fidelity
+7. **Visualize noise patterns**: Identify problematic qubits and gates
+8. **Apply targeted mitigation**: Focus on dominant error sources
+9. **Validate mitigation**: Verify that mitigation improves results
+10. **Keep circuits shallow**: Minimize noise accumulation

scientific-skills/cirq/references/simulation.md

@@ -0,0 +1,350 @@
+# Simulation in Cirq
+
+This guide covers quantum circuit simulation, including exact and noisy simulations, parameter sweeps, and the Quantum Virtual Machine (QVM).
+
+## Exact Simulation
+
+### Basic Simulation
+
+```python
+import cirq
+import numpy as np
+
+# Create circuit
+q0, q1 = cirq.LineQubit.range(2)
+circuit = cirq.Circuit(
+    cirq.H(q0),
+    cirq.CNOT(q0, q1),
+    cirq.measure(q0, q1, key='result')
+)
+
+# Simulate
+simulator = cirq.Simulator()
+result = simulator.run(circuit, repetitions=1000)
+
+# Get measurement results
+print(result.histogram(key='result'))
+```
+
+### State Vector Simulation
+
+```python
+# Simulate without measurement to get final state
+simulator = cirq.Simulator()
+result = simulator.simulate(circuit_without_measurement)
+
+# Access state vector
+state_vector = result.final_state_vector
+print(f"State vector: {state_vector}")
+
+# Get amplitudes
+print(f"Amplitude of |00⟩: {state_vector[0]}")
+print(f"Amplitude of |11⟩: {state_vector[3]}")
+```
+
+### Density Matrix Simulation
+
+```python
+# Use density matrix simulator for mixed states
+simulator = cirq.DensityMatrixSimulator()
+result = simulator.simulate(circuit)
+
+# Access density matrix
+density_matrix = result.final_density_matrix
+print(f"Density matrix shape: {density_matrix.shape}")
+```
+
+### Step-by-Step Simulation
+
+```python
+# Simulate moment-by-moment
+simulator = cirq.Simulator()
+for step in simulator.simulate_moment_steps(circuit):
+    print(f"State after moment {step.moment}: {step.state_vector()}")
+```
+
+## Sampling and Measurements
+
+### Run Multiple Shots
+
+```python
+# Run circuit multiple times
+result = simulator.run(circuit, repetitions=10000)
+
+# Access measurement counts
+counts = result.histogram(key='result')
+print(f"Measurement counts: {counts}")
+
+# Get raw measurements
+measurements = result.measurements['result']
+print(f"Shape: {measurements.shape}")  # (repetitions, num_qubits)
+```
+
+### Expectation Values
+
+```python
+# Measure observable expectation value
+from cirq import PauliString
+
+observable = PauliString({q0: cirq.Z, q1: cirq.Z})
+result = simulator.simulate_expectation_values(
+    circuit,
+    observables=[observable]
+)
+print(f"⟨ZZ⟩ = {result[0]}")
+```
+
+## Parameter Sweeps
+
+### Sweep Over Parameters
+
+```python
+import sympy
+
+# Create parameterized circuit
+theta = sympy.Symbol('theta')
+q = cirq.LineQubit(0)
+circuit = cirq.Circuit(
+    cirq.ry(theta)(q),
+    cirq.measure(q, key='m')
+)
+
+# Define parameter sweep
+sweep = cirq.Linspace(key='theta', start=0, stop=2*np.pi, length=50)
+
+# Run sweep
+simulator = cirq.Simulator()
+results = simulator.run_sweep(circuit, params=sweep, repetitions=1000)
+
+# Process results
+for params, result in zip(sweep, results):
+    theta_val = params['theta']
+    counts = result.histogram(key='m')
+    print(f"θ={theta_val:.2f}: {counts}")
+```
+
+### Multiple Parameters
+
+```python
+# Sweep over multiple parameters
+theta = sympy.Symbol('theta')
+phi = sympy.Symbol('phi')
+
+circuit = cirq.Circuit(
+    cirq.ry(theta)(q0),
+    cirq.rz(phi)(q1)
+)
+
+# Product sweep (all combinations)
+sweep = cirq.Product(
+    cirq.Linspace('theta', 0, np.pi, 10),
+    cirq.Linspace('phi', 0, 2*np.pi, 10)
+)
+
+results = simulator.run_sweep(circuit, params=sweep, repetitions=100)
+```
+
+### Zip Sweep (Paired Parameters)
+
+```python
+# Sweep parameters together
+sweep = cirq.Zip(
+    cirq.Linspace('theta', 0, np.pi, 20),
+    cirq.Linspace('phi', 0, 2*np.pi, 20)
+)
+
+results = simulator.run_sweep(circuit, params=sweep, repetitions=100)
+```
+
+## Noisy Simulation
+
+### Adding Noise Channels
+
+```python
+# Create noisy circuit
+noisy_circuit = circuit.with_noise(cirq.depolarize(p=0.01))
+
+# Simulate noisy circuit
+simulator = cirq.DensityMatrixSimulator()
+result = simulator.run(noisy_circuit, repetitions=1000)
+```
+
+### Custom Noise Models
+
+```python
+# Apply different noise to different gates
+noise_model = cirq.NoiseModel.from_noise_model_like(
+    cirq.ConstantQubitNoiseModel(cirq.depolarize(0.01))
+)
+
+# Simulate with noise model
+result = cirq.DensityMatrixSimulator(noise=noise_model).run(
+    circuit, repetitions=1000
+)
+```
+
+See `noise.md` for comprehensive noise modeling details.
+
+## State Histograms
+
+### Visualize Results
+
+```python
+import matplotlib.pyplot as plt
+
+# Get histogram
+result = simulator.run(circuit, repetitions=1000)
+counts = result.histogram(key='result')
+
+# Plot
+plt.bar(counts.keys(), counts.values())
+plt.xlabel('State')
+plt.ylabel('Counts')
+plt.title('Measurement Results')
+plt.show()
+```
+
+### State Probability Distribution
+
+```python
+# Get state vector
+result = simulator.simulate(circuit_without_measurement)
+state_vector = result.final_state_vector
+
+# Compute probabilities
+probabilities = np.abs(state_vector) ** 2
+
+# Plot
+plt.bar(range(len(probabilities)), probabilities)
+plt.xlabel('Basis State Index')
+plt.ylabel('Probability')
+plt.show()
+```
+
+## Quantum Virtual Machine (QVM)
+
+QVM simulates realistic quantum hardware with device-specific constraints and noise.
+
+### Using Virtual Devices
+
+```python
+# Use a virtual Google device
+import cirq_google
+
+# Get virtual device
+device = cirq_google.Sycamore
+
+# Create circuit on device
+qubits = device.metadata.qubit_set
+circuit = cirq.Circuit(device=device)
+
+# Add operations respecting device constraints
+circuit.append(cirq.CZ(qubits[0], qubits[1]))
+
+# Validate circuit against device
+device.validate_circuit(circuit)
+```
+
+### Noisy Virtual Hardware
+
+```python
+# Simulate with device noise
+processor = cirq_google.get_engine().get_processor('weber')
+noise_props = processor.get_device_specification()
+
+# Create realistic noisy simulator
+noisy_sim = cirq.DensityMatrixSimulator(
+    noise=cirq_google.NoiseModelFromGoogleNoiseProperties(noise_props)
+)
+
+result = noisy_sim.run(circuit, repetitions=1000)
+```
+
+## Advanced Simulation Techniques
+
+### Custom Initial State
+
+```python
+# Start from custom state
+initial_state = np.array([1, 0, 0, 1]) / np.sqrt(2)  # |00⟩ + |11⟩
+
+simulator = cirq.Simulator()
+result = simulator.simulate(circuit, initial_state=initial_state)
+```
+
+### Partial Trace
+
+```python
+# Trace out subsystems
+result = simulator.simulate(circuit)
+full_state = result.final_state_vector
+
+# Compute reduced density matrix for first qubit
+from cirq import partial_trace
+reduced_dm = partial_trace(result.final_density_matrix, keep_indices=[0])
+```
+
+### Intermediate State Access
+
+```python
+# Get state at specific moment
+simulator = cirq.Simulator()
+for i, step in enumerate(simulator.simulate_moment_steps(circuit)):
+    if i == 5:  # After 5th moment
+        state = step.state_vector()
+        print(f"State after moment 5: {state}")
+        break
+```
+
+## Simulation Performance
+
+### Optimizing Large Simulations
+
+1. **Use state vector for pure states**: Faster than density matrix
+2. **Avoid density matrix when possible**: Exponentially more expensive
+3. **Batch parameter sweeps**: More efficient than individual runs
+4. **Use appropriate repetitions**: Balance accuracy vs computation time
+
+```python
+# Efficient: Single sweep
+results = simulator.run_sweep(circuit, params=sweep, repetitions=100)
+
+# Inefficient: Multiple individual runs
+results = [simulator.run(circuit, param_resolver=p, repetitions=100)
+           for p in sweep]
+```
+
+### Memory Considerations
+
+```python
+# For large systems, monitor state vector size
+n_qubits = 20
+state_size = 2**n_qubits * 16  # bytes (complex128)
+print(f"State vector size: {state_size / 1e9:.2f} GB")
+```
+
+## Stabilizer Simulation
+
+For circuits with only Clifford gates, use efficient stabilizer simulation:
+
+```python
+# Clifford circuit (H, S, CNOT)
+circuit = cirq.Circuit(
+    cirq.H(q0),
+    cirq.S(q1),
+    cirq.CNOT(q0, q1)
+)
+
+# Use stabilizer simulator (exponentially faster)
+simulator = cirq.CliffordSimulator()
+result = simulator.run(circuit, repetitions=1000)
+```
+
+## Best Practices
+
+1. **Choose appropriate simulator**: Use Simulator for pure states, DensityMatrixSimulator for mixed states
+2. **Use parameter sweeps**: More efficient than running individual circuits
+3. **Validate circuits**: Check circuit validity before long simulations
+4. **Monitor resource usage**: Track memory for large-scale simulations
+5. **Use stabilizer simulation**: When circuits contain only Clifford gates
+6. **Save intermediate results**: For long parameter sweeps or optimization runs

scientific-skills/cirq/references/transformation.md

@@ -0,0 +1,416 @@
+# Circuit Transformations
+
+This guide covers circuit optimization, compilation, and manipulation using Cirq's transformation framework.
+
+## Transformer Framework
+
+### Basic Transformers
+
+```python
+import cirq
+
+# Example circuit
+qubits = cirq.LineQubit.range(3)
+circuit = cirq.Circuit(
+    cirq.H(qubits[0]),
+    cirq.CNOT(qubits[0], qubits[1]),
+    cirq.CNOT(qubits[1], qubits[2])
+)
+
+# Apply built-in transformer
+from cirq.transformers import optimize_for_target_gateset
+
+# Optimize to specific gate set
+optimized = optimize_for_target_gateset(
+    circuit,
+    gateset=cirq.SqrtIswapTargetGateset()
+)
+```
+
+### Merge Single-Qubit Gates
+
+```python
+from cirq.transformers import merge_single_qubit_gates_to_phxz
+
+# Circuit with multiple single-qubit gates
+circuit = cirq.Circuit(
+    cirq.H(q),
+    cirq.T(q),
+    cirq.S(q),
+    cirq.H(q)
+)
+
+# Merge into single operation
+merged = merge_single_qubit_gates_to_phxz(circuit)
+```
+
+### Drop Negligible Operations
+
+```python
+from cirq.transformers import drop_negligible_operations
+
+# Remove gates below threshold
+circuit_with_small_rotations = cirq.Circuit(
+    cirq.rz(1e-10)(q),  # Very small rotation
+    cirq.H(q)
+)
+
+cleaned = drop_negligible_operations(circuit_with_small_rotations, atol=1e-8)
+```
+
+## Custom Transformers
+
+### Transformer Decorator
+
+```python
+from cirq.transformers import transformer_api
+
+@transformer_api.transformer
+def remove_z_gates(circuit: cirq.Circuit) -> cirq.Circuit:
+    """Remove all Z gates from circuit."""
+    new_moments = []
+    for moment in circuit:
+        new_ops = [op for op in moment if not isinstance(op.gate, cirq.ZPowGate)]
+        if new_ops:
+            new_moments.append(cirq.Moment(new_ops))
+    return cirq.Circuit(new_moments)
+
+# Use custom transformer
+transformed = remove_z_gates(circuit)
+```
+
+### Transformer Class
+
+```python
+from cirq.transformers import transformer_primitives
+
+class HToRyTransformer(transformer_primitives.Transformer):
+    """Replace H gates with Ry(π/2)."""
+
+    def __call__(self, circuit: cirq.Circuit, *, context=None) -> cirq.Circuit:
+        def map_op(op: cirq.Operation, _) -> cirq.OP_TREE:
+            if isinstance(op.gate, cirq.HPowGate):
+                return cirq.ry(np.pi/2)(op.qubits[0])
+            return op
+
+        return transformer_primitives.map_operations(
+            circuit,
+            map_op,
+            deep=True
+        ).unfreeze(copy=False)
+
+# Apply transformer
+transformer = HToRyTransformer()
+result = transformer(circuit)
+```
+
+## Gate Decomposition
+
+### Decompose to Target Gateset
+
+```python
+from cirq.transformers import optimize_for_target_gateset
+
+# Decompose to CZ + single-qubit rotations
+target_gateset = cirq.CZTargetGateset()
+decomposed = optimize_for_target_gateset(circuit, gateset=target_gateset)
+
+# Decompose to √iSWAP gates
+sqrt_iswap_gateset = cirq.SqrtIswapTargetGateset()
+decomposed = optimize_for_target_gateset(circuit, gateset=sqrt_iswap_gateset)
+```
+
+### Custom Gate Decomposition
+
+```python
+class Toffoli(cirq.Gate):
+    def _num_qubits_(self):
+        return 3
+
+    def _decompose_(self, qubits):
+        """Decompose Toffoli into basic gates."""
+        c1, c2, t = qubits
+        return [
+            cirq.H(t),
+            cirq.CNOT(c2, t),
+            cirq.T(t)**-1,
+            cirq.CNOT(c1, t),
+            cirq.T(t),
+            cirq.CNOT(c2, t),
+            cirq.T(t)**-1,
+            cirq.CNOT(c1, t),
+            cirq.T(c2),
+            cirq.T(t),
+            cirq.H(t),
+            cirq.CNOT(c1, c2),
+            cirq.T(c1),
+            cirq.T(c2)**-1,
+            cirq.CNOT(c1, c2)
+        ]
+
+# Use decomposition
+circuit = cirq.Circuit(Toffoli()(q0, q1, q2))
+decomposed = cirq.decompose(circuit)
+```
+
+## Circuit Optimization
+
+### Eject Z Gates
+
+```python
+from cirq.transformers import eject_z
+
+# Move Z gates to end of circuit
+circuit = cirq.Circuit(
+    cirq.H(q0),
+    cirq.Z(q0),
+    cirq.CNOT(q0, q1)
+)
+
+ejected = eject_z(circuit)
+```
+
+### Eject Phase Gates
+
+```python
+from cirq.transformers import eject_phased_paulis
+
+# Consolidate phase gates
+optimized = eject_phased_paulis(circuit, atol=1e-8)
+```
+
+### Drop Empty Moments
+
+```python
+from cirq.transformers import drop_empty_moments
+
+# Remove moments with no operations
+cleaned = drop_empty_moments(circuit)
+```
+
+### Align Measurements
+
+```python
+from cirq.transformers import dephase_measurements
+
+# Move measurements to end and remove operations after
+aligned = dephase_measurements(circuit)
+```
+
+## Circuit Compilation
+
+### Compile for Hardware
+
+```python
+import cirq_google
+
+# Get device specification
+device = cirq_google.Sycamore
+
+# Compile circuit to device
+from cirq.transformers import optimize_for_target_gateset
+
+compiled = optimize_for_target_gateset(
+    circuit,
+    gateset=cirq_google.SycamoreTargetGateset()
+)
+
+# Validate compiled circuit
+device.validate_circuit(compiled)
+```
+
+### Two-Qubit Gate Compilation
+
+```python
+# Compile to specific two-qubit gate
+from cirq import two_qubit_to_cz
+
+# Convert all two-qubit gates to CZ
+cz_circuit = cirq.Circuit()
+for moment in circuit:
+    for op in moment:
+        if len(op.qubits) == 2:
+            cz_circuit.append(two_qubit_to_cz(op))
+        else:
+            cz_circuit.append(op)
+```
+
+## Qubit Routing
+
+### Route Circuit to Device Topology
+
+```python
+from cirq.transformers import route_circuit
+
+# Define device connectivity
+device_graph = cirq.NamedTopology(
+    {
+        (0, 0): [(0, 1), (1, 0)],
+        (0, 1): [(0, 0), (1, 1)],
+        (1, 0): [(0, 0), (1, 1)],
+        (1, 1): [(0, 1), (1, 0)]
+    }
+)
+
+# Route logical qubits to physical qubits
+routed_circuit = route_circuit(
+    circuit,
+    device_graph=device_graph,
+    routing_algo=cirq.RouteCQC(device_graph)
+)
+```
+
+### SWAP Network Insertion
+
+```python
+# Manually insert SWAPs for routing
+def insert_swaps(circuit, swap_locations):
+    """Insert SWAP gates at specified locations."""
+    new_circuit = cirq.Circuit()
+    moment_idx = 0
+
+    for i, moment in enumerate(circuit):
+        if i in swap_locations:
+            q0, q1 = swap_locations[i]
+            new_circuit.append(cirq.SWAP(q0, q1))
+        new_circuit.append(moment)
+
+    return new_circuit
+```
+
+## Advanced Transformations
+
+### Unitary Compilation
+
+```python
+import scipy.linalg
+
+# Compile arbitrary unitary to gate sequence
+def compile_unitary(unitary, qubits):
+    """Compile 2x2 unitary using KAK decomposition."""
+    from cirq.linalg import kak_decomposition
+
+    decomp = kak_decomposition(unitary)
+    operations = []
+
+    # Add single-qubit gates before
+    operations.append(cirq.MatrixGate(decomp.single_qubit_operations_before[0])(qubits[0]))
+    operations.append(cirq.MatrixGate(decomp.single_qubit_operations_before[1])(qubits[1]))
+
+    # Add interaction (two-qubit) part
+    x, y, z = decomp.interaction_coefficients
+    operations.append(cirq.XXPowGate(exponent=x/np.pi)(qubits[0], qubits[1]))
+    operations.append(cirq.YYPowGate(exponent=y/np.pi)(qubits[0], qubits[1]))
+    operations.append(cirq.ZZPowGate(exponent=z/np.pi)(qubits[0], qubits[1]))
+
+    # Add single-qubit gates after
+    operations.append(cirq.MatrixGate(decomp.single_qubit_operations_after[0])(qubits[0]))
+    operations.append(cirq.MatrixGate(decomp.single_qubit_operations_after[1])(qubits[1]))
+
+    return operations
+```
+
+### Circuit Simplification
+
+```python
+from cirq.transformers import (
+    merge_k_qubit_unitaries,
+    merge_single_qubit_gates_to_phxz
+)
+
+# Merge adjacent single-qubit gates
+simplified = merge_single_qubit_gates_to_phxz(circuit)
+
+# Merge adjacent k-qubit unitaries
+simplified = merge_k_qubit_unitaries(circuit, k=2)
+```
+
+### Commutation-Based Optimization
+
+```python
+# Commute Z gates through CNOT
+def commute_z_through_cnot(circuit):
+    """Move Z gates through CNOT gates."""
+    new_moments = []
+
+    for moment in circuit:
+        ops = list(moment)
+        # Find Z gates before CNOT
+        z_ops = [op for op in ops if isinstance(op.gate, cirq.ZPowGate)]
+        cnot_ops = [op for op in ops if isinstance(op.gate, cirq.CXPowGate)]
+
+        # Apply commutation rules
+        # Z on control commutes, Z on target anticommutes
+        # (simplified logic here)
+
+        new_moments.append(cirq.Moment(ops))
+
+    return cirq.Circuit(new_moments)
+```
+
+## Transformation Pipelines
+
+### Compose Multiple Transformers
+
+```python
+from cirq.transformers import transformer_api
+
+# Build transformation pipeline
+@transformer_api.transformer
+def optimization_pipeline(circuit: cirq.Circuit) -> cirq.Circuit:
+    # Step 1: Merge single-qubit gates
+    circuit = merge_single_qubit_gates_to_phxz(circuit)
+
+    # Step 2: Drop negligible operations
+    circuit = drop_negligible_operations(circuit)
+
+    # Step 3: Eject Z gates
+    circuit = eject_z(circuit)
+
+    # Step 4: Drop empty moments
+    circuit = drop_empty_moments(circuit)
+
+    return circuit
+
+# Apply pipeline
+optimized = optimization_pipeline(circuit)
+```
+
+## Validation and Analysis
+
+### Circuit Depth Reduction
+
+```python
+# Measure circuit depth before and after
+print(f"Original depth: {len(circuit)}")
+optimized = optimization_pipeline(circuit)
+print(f"Optimized depth: {len(optimized)}")
+```
+
+### Gate Count Analysis
+
+```python
+def count_gates(circuit):
+    """Count gates by type."""
+    counts = {}
+    for moment in circuit:
+        for op in moment:
+            gate_type = type(op.gate).__name__
+            counts[gate_type] = counts.get(gate_type, 0) + 1
+    return counts
+
+original_counts = count_gates(circuit)
+optimized_counts = count_gates(optimized)
+print(f"Original: {original_counts}")
+print(f"Optimized: {optimized_counts}")
+```
+
+## Best Practices
+
+1. **Start with high-level transformers**: Use built-in transformers before writing custom ones
+2. **Chain transformers**: Apply multiple optimizations in sequence
+3. **Validate after transformation**: Ensure circuit correctness and device compatibility
+4. **Measure improvement**: Track depth and gate count reduction
+5. **Use appropriate gatesets**: Match target hardware capabilities
+6. **Consider commutativity**: Exploit gate commutation for optimization
+7. **Test on small circuits first**: Verify transformers work correctly before scaling