test_error.py

import numpy as np
from matplotlib import pyplot as plt
from qiskit_aer import AerSimulator

from errors import model_depolarizing_error, _construct_grover_solver, model_thermal_error_2q, model_thermal_error_1q, \
    build_combined_noise_model
from subset_sum_oracle import SubsetSumOracle
from verify import solve_subset_sum_with_bitstrings, check_counts_statistically_accurate


def run_2d_error_sweep(numbers: list[int], target: int, shots: int = 1024):
    """
    Run Grover's algorithm across a grid of depolarizing and thermal error values,
    then plot a contour map of accuracy.
    """
    # Define error ranges
    depolarizing_errs = [None, 0.0001, 0.00025, 0.001, 0.002, 0.003, 0.005, 0.0075, 0.01]
    thermal_t1_means = [200e3, 150e3, 100e3, 75e3, 50e3, 25e3, 15e3, 10e3, 5e3]  # Lower = more thermal noise

    # Get correct solutions for verification
    correct_bitstrings = [bs for _, bs in solve_subset_sum_with_bitstrings(numbers, target)]
    n_solutions = len(correct_bitstrings)

    # Create meshgrid for plotting
    # Convert None to 0 for plotting purposes
    depol_values = [0 if e is None else e for e in depolarizing_errs]
    thermal_values = thermal_t1_means

    # Store accuracy results
    accuracy_grid = np.zeros((len(thermal_t1_means), len(depolarizing_errs)))

    total_runs = len(depolarizing_errs) * len(thermal_t1_means)
    run_count = 0

    for i, t1_mean in enumerate(thermal_t1_means):
        for j, depol_err in enumerate(depolarizing_errs):
            run_count += 1
            print(f"Run {run_count}/{total_runs}: depol={depol_err}, T1_mean={t1_mean}")

            oracle = SubsetSumOracle(numbers, target)
            n_qubits = oracle.n_items + oracle.sum_bits + 1
            # Build noise model
            noise_model = build_combined_noise_model(depol_err, t1_mean, n_qubits)

            simulator = AerSimulator(noise_model=noise_model)
            solver = _construct_grover_solver(numbers, target, simulator)

            result = solver.run_known_solutions(n_solutions=n_solutions, shots=shots)

            if result is None:
                accuracy_grid[i, j] = 0  # No solution found
            else:
                accuracy, _ = check_counts_statistically_accurate(result, correct_bitstrings)
                # Cap infinity at a reasonable value for plotting
                if accuracy == float('inf'):
                    accuracy = 500
                if accuracy < 0:
                    accuracy = 0
                accuracy_grid[i, j] = accuracy

    # Plot contour map
    plot_error_contour(depol_values, thermal_values, accuracy_grid)


def plot_error_contour(depol_values, thermal_values, accuracy_grid):
    """Plot 2D contour map of accuracy vs error parameters."""
    fig, ax = plt.subplots(figsize=(12, 9))

    X, Y = np.meshgrid(depol_values, thermal_values)

    # Create filled contour plot
    levels = np.linspace(accuracy_grid.min(), min(accuracy_grid.max(), 500), 20)
    contourf = ax.contourf(X, Y, accuracy_grid, levels=levels, cmap='RdYlGn')

    # Add contour lines
    contour = ax.contour(X, Y, accuracy_grid, levels=levels, colors='black', linewidths=0.5, alpha=0.5)
    ax.clabel(contour, inline=True, fontsize=8, fmt='%.0f%%')

    # Colorbar
    cbar = plt.colorbar(contourf, ax=ax)
    cbar.set_label('Accuracy Gap (%)', fontsize=12)

    # Labels
    ax.set_xlabel('Depolarizing Error Rate', fontsize=12)
    ax.set_ylabel('Thermal T1 Mean (ns)', fontsize=12)
    ax.set_title('Grover Algorithm Accuracy vs Noise Parameters', fontsize=14)

    # Use log scale for better visualization
    ax.set_yscale('log')

    plt.tight_layout()
    plt.show()


def find_error_rate(numbers, target, min_percentage_gap=30.0, shots=1024,
                    err_min=0.0, err_max=0.005, tolerance=0.0001, samples_per_test=3):
    """
    Find the maximum depolarizing error rate that still achieves the minimum accuracy gap.
    Uses binary search with multiple samples to handle randomness in quantum simulation.

    Args:
        numbers: List of integers for subset sum problem
        target: Target sum
        min_percentage_gap: Minimum required accuracy gap (%)
        shots: Number of shots per simulation
        err_min: Lower bound of error rate search range
        err_max: Upper bound of error rate search range
        tolerance: Stop when search range is smaller than this
        samples_per_test: Number of samples to take at each error rate to reduce noise

    Returns:
        The estimated maximum error rate that maintains accuracy above min_percentage_gap
    """
    correct_bitstrings = [bs for _, bs in solve_subset_sum_with_bitstrings(numbers, target)]
    n_solutions = len(correct_bitstrings)

    def test_error_rate(depol_err):
        """Test an error rate multiple times and return the fraction of passes."""
        passes = 0
        total_accuracy = 0.0

        for _ in range(samples_per_test):
            noise_model = model_depolarizing_error(err=depol_err if depol_err > 0 else None)
            simulator = AerSimulator(noise_model=noise_model)
            solver = _construct_grover_solver(numbers, target, simulator)
            result = solver.run_known_solutions(n_solutions=n_solutions, shots=shots)

            accuracy, _ = check_counts_statistically_accurate(result, correct_bitstrings)
            total_accuracy += accuracy if accuracy != float('inf') else 500

            if accuracy > min_percentage_gap:
                passes += 1

        avg_accuracy = total_accuracy / samples_per_test
        pass_rate = passes / samples_per_test
        return pass_rate, avg_accuracy

    low = err_min
    high = err_max
    best_passing_err = None

    # First, verify that low error rate passes (sanity check)
    pass_rate, avg_acc = test_error_rate(low)
    print(f"Initial check at err={low:.6f}: pass_rate={pass_rate:.2f}, avg_accuracy={avg_acc:.1f}%")
    if pass_rate >= 0.5:
        best_passing_err = low

    iteration = 0
    max_iterations = 20  # Safety limit

    while high - low > tolerance and iteration < max_iterations:
        iteration += 1
        mid = (low + high) / 2

        pass_rate, avg_accuracy = test_error_rate(mid)
        print(f"Iteration {iteration}: err={mid:.6f}, pass_rate={pass_rate:.2f}, avg_accuracy={avg_accuracy:.1f}%")

        # Decision logic with hysteresis to handle noise
        if pass_rate >= 0.5:  # Majority of samples pass
            # This error rate works, try higher
            best_passing_err = mid
            low = mid
        else:
            # This error rate fails, try lower
            high = mid

    # Return the best passing error rate found, or do a final refinement
    if best_passing_err is not None:
        print(f"Found error rate: {best_passing_err:.6f} (range: [{low:.6f}, {high:.6f}])")
        return best_passing_err
    else:
        print("No passing error rate found in range")
        return None


if __name__ == '__main__':
    numbers = [1, 2, 3]
    target = 5
    # find_error_rate(numbers, target)
    run_2d_error_sweep(numbers, target)