aiming-lab/autoresearchclaw/quantum-qiskit

import numpy as np
from qiskit import QuantumCircuit
from qiskit.circuit import ParameterVector
from qiskit.circuit.library import (
    ZFeatureMap,
    ZZFeatureMap,
    StatePreparation,
    EfficientSU2,
)
from qiskit.primitives import StatevectorSampler, StatevectorEstimator  # V2 primitives
from qiskit.quantum_info import Statevector, SparsePauliOp
from qiskit_aer import AerSimulator
from qiskit_algorithms.optimizers import SPSA, COBYLA, L_BFGS_B, ADAM
from qiskit_algorithms.utils import algorithm_globals
from qiskit_machine_learning.algorithms.classifiers import VQC

from qiskit_nature.units import DistanceUnit
from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.second_q.mappers import ParityMapper, JordanWignerMapper

def build_angle_encoding(num_features: int) -> QuantumCircuit:
    """Hadamard plus single-qubit Z-rotation per feature.
    Mathematically equivalent to ZFeatureMap(reps=1).
    """
    return ZFeatureMap(feature_dimension=num_features, reps=1)
def build_amplitude_encoding(num_features: int):
    """Load an L2-normalized, zero-padded input as the amplitudes of a
    quantum state. The encoding uses ceil(log2(num_features)) qubits.
    Returns (circuit, parameter_vector, num_qubits). The caller binds
    parameters per-sample via the helper below.
    """
    num_qubits = int(np.ceil(np.log2(max(num_features, 2))))
    full_dim = 2 ** num_qubits
    params = ParameterVector("x_amp", full_dim)
    qc = QuantumCircuit(num_qubits)
    qc.append(StatePreparation(list(params)), range(num_qubits))
    return qc, params, num_qubits
def amplitude_binding(x: np.ndarray, params, num_qubits: int) -> dict:
    """Build the parameter-value dict for a single input sample."""
    x_norm = x / max(float(np.linalg.norm(x)), 1e-12)
    padded = np.zeros(2 ** num_qubits, dtype=np.float64)
    padded[: len(x_norm)] = x_norm
    padded = padded / max(float(np.linalg.norm(padded)), 1e-12)
    return {params[i]: float(padded[i]) for i in range(len(padded))}
def build_zz_feature_map(num_features: int) -> QuantumCircuit:
    """Two repetitions of Hadamard plus pairwise ZZ entangling rotations."""
    return ZZFeatureMap(
        feature_dimension=num_features, reps=2, entanglement="linear"
    )

def assert_different_output_states(qc_a, qc_b, x, tol: float = 1e-6):
    sv_a = Statevector(qc_a.assign_parameters(x))
    sv_b = Statevector(qc_b.assign_parameters(x))
    diff = float(np.linalg.norm(sv_a.data - sv_b.data))
    assert diff > tol, f"encoders produced identical states (diff={diff})"

def build_ansatz(num_qubits: int, reps: int = 2) -> QuantumCircuit:
    """Hardware-efficient ansatz with alternating Pauli rotations
    and a linear chain of CNOT entanglers. Trainable parameter count
    is (reps + 1) * num_qubits for the default su2_gates = ['ry']."""
    return EfficientSU2(
        num_qubits=num_qubits, reps=reps, entanglement="linear"
    )

def train_vqc(
    feature_map: QuantumCircuit,
    ansatz: QuantumCircuit,
    X_train: np.ndarray,
    y_train: np.ndarray,
    seed: int,
    maxiter: int = 200,
) -> VQC:
    algorithm_globals.random_seed = seed
    vqc = VQC(
        feature_map=feature_map,
        ansatz=ansatz,
        loss="cross_entropy",
        optimizer=COBYLA(maxiter=maxiter),
        sampler=StatevectorSampler(seed=seed),
    )
    vqc.fit(X_train, y_train)
    return vqc

def build_h2_hamiltonian(bond_length_angstrom: float):
    driver = PySCFDriver(
        atom=f"H 0 0 0; H 0 0 {bond_length_angstrom}",
        basis="sto3g",
        charge=0,
        spin=0,
        unit=DistanceUnit.ANGSTROM,
    )
    problem = driver.run()
    num_particles = tuple(problem.num_particles)         # (1, 1) for H2
    mapper = ParityMapper(num_particles=num_particles)   # 2-qubit reduction
    qubit_op = mapper.map(problem.hamiltonian.second_q_op())
    e_nuclear = float(problem.nuclear_repulsion_energy)
    # H2 in STO-3G with parity mapping plus 2-qubit reduction produces a
    # 2-qubit Hamiltonian (not 4-qubit).
    return qubit_op, e_nuclear
def run_vqe(qubit_op, e_nuclear, optimizer_name: str, seed: int):
    algorithm_globals.random_seed = seed
    rng = np.random.RandomState(seed)
    ansatz = build_ansatz(num_qubits=qubit_op.num_qubits, reps=2)
    initial_point = rng.normal(0.0, 0.1, ansatz.num_parameters)
    estimator = StatevectorEstimator(seed=seed)
    energy_history: list[tuple[int, float]] = []
    def energy(theta: np.ndarray) -> float:
        bound = ansatz.assign_parameters(theta)
        result = estimator.run([(bound, qubit_op)]).result()
        e = float(result[0].data.evs) + e_nuclear
        energy_history.append((len(energy_history) + 1, e))
        return e
    optimizers = {
        "spsa": SPSA(maxiter=200),
        "cobyla": COBYLA(maxiter=200, rhobeg=0.1, tol=1e-4),
        "lbfgsb": L_BFGS_B(maxiter=100, ftol=1e-6),
        "adam": ADAM(maxiter=200, lr=0.05, beta_1=0.9, beta_2=0.999),
    }
    if optimizer_name not in optimizers:
        raise ValueError(f"unknown optimizer: {optimizer_name}")
    result = optimizers[optimizer_name].minimize(energy, initial_point)
    return result, energy_history

from collections import deque
def cumulative_shots_to_threshold(
    energy_history: list[tuple[int, float]],
    e_target: float,
    threshold: float = 0.0016,    # 1.6 mHa
    shots_per_eval: int = 1024,
    window: int = 5,
) -> int | None:
    """Return cumulative shots at the first point where the running mean
    over `window` evaluations stays within `threshold` of `e_target` for
    `window` consecutive windows. Return None if never reached."""
    buf = deque(maxlen=window)
    streak = 0
    for eval_count, energy in energy_history:
        buf.append(energy)
        if len(buf) < window:
            continue
        if abs(sum(buf) / window - e_target) <= threshold:
            streak += 1
            if streak >= window:
                return eval_count * shots_per_eval
        else:
            streak = 0
    return None

def encode_features_to_circuit(x: np.ndarray, n_qubits: int) -> QuantumCircuit:
    """Per-feature embedding equivalent to phi(x) = [cos(pi*x/2), sin(pi*x/2)].
    Apply RY(pi * x_i) on qubit i so |0> maps to cos(pi*x_i/2)|0> + sin(pi*x_i/2)|1>."""
    qc = QuantumCircuit(n_qubits)
    for i in range(n_qubits):
        qc.ry(float(x[i]) * np.pi, i)
    return qc
def build_mps_ansatz(n_qubits: int, reps_for_chi: int) -> QuantumCircuit:
    """Linear-chain entangling ansatz; effective bond dimension <= 2 ** reps_for_chi.
    reps_for_chi=4 covers chi up to 16."""
    return EfficientSU2(
        num_qubits=n_qubits, reps=reps_for_chi, entanglement="linear"
    )
def mps_class_logits(
    x: np.ndarray,
    theta: np.ndarray,
    ansatz: QuantumCircuit,
    n_classes: int,
    max_bond: int = 16,
) -> np.ndarray:
    """Return one logit per class, computed via AerSimulator MPS method.
    Each class c corresponds to a Pauli observable acting on the first
    ceil(log2(n_classes)) qubits with sign pattern fixed by the bits of c."""
    import math
    n_qubits = ansatz.num_qubits
    sim = AerSimulator(
        method="matrix_product_state",
        matrix_product_state_max_bond_dimension=int(max_bond),
    )
    bound_ansatz = ansatz.assign_parameters(theta)
    qc = encode_features_to_circuit(x, n_qubits)
    qc.compose(bound_ansatz, inplace=True)
    n_label_qubits = max(1, math.ceil(math.log2(n_classes)))
    logits = []
    for c in range(n_classes):
        pauli = list("I" * n_qubits)
        for bit_idx in range(n_label_qubits):
            if (c >> bit_idx) & 1:
                pauli[bit_idx] = "Z"
        obs = SparsePauliOp.from_list([("".join(pauli[::-1]), 1.0)])
        qc_with_save = qc.copy()
        qc_with_save.save_expectation_value(obs, list(range(n_qubits)))
        result = sim.run(qc_with_save).result()
        logits.append(float(result.data(0)["expectation_value"]))
    return np.array(logits)