Task2_6/Task2_6.py - master - timshk/Shalaev_Stepik

Task2_6/Task2_6.py
407 строк18 KB
timshk
Stepik задания
01 июн 2026, 18:58
01 июн 2026, 18:582e75834
"""
БЕНЧМАРК ZNE: СРАВНЕНИЕ МИТИГАЦИИ НА ШУМОВОМ СИМУЛЯТОРЕ
Тестовая схема: состояние Белла |Φ⁺⟩ = (|00⟩ + |11⟩)/√2
Наблюдаемая: ⟨Z₀Z₁⟩ (идеальное значение = 1.0)

Метрики:
- Точность: абсолютная ошибка |1.0 - ⟨ZZ⟩|
- Стоимость: общее число запусков и время (очередь + выполнение)
- Устойчивость: R² линейной регрессии зашумленных значений по scale_factors
"""

import numpy as np
import time
import matplotlib.pyplot as plt
from qiskit import QuantumCircuit, transpile
from qiskit.quantum_info import SparsePauliOp
from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel, depolarizing_error, ReadoutError
import warnings
warnings.filterwarnings('ignore')

# ============================================================================
# КОНФИГУРАЦИЯ
# ============================================================================
SHOTS = 4096                     # Количество выстрелов для каждого запуска
SCALE_FACTORS = [1, 2, 3]        # Коэффициенты масштабирования шума
OBSERVABLE = SparsePauliOp.from_list([("ZZ", 1.0)])  # Наблюдаемая ⟨Z₀Z₁⟩


# ============================================================================
# 1. ВСПОМОГАТЕЛЬНАЯ ФУНКЦИЯ ДЛЯ РАСЧЁТА R² (БЕЗ SKLEARN)
# ============================================================================
def r2_score_manual(y_true, y_pred):
    """
    Вычисляет коэффициент детерминации R² вручную.
    R² = 1 - (SS_res / SS_tot)
    """
    y_true = np.array(y_true)
    y_pred = np.array(y_pred)
    
    ss_res = np.sum((y_true - y_pred) ** 2)
    ss_tot = np.sum((y_true - np.mean(y_true)) ** 2)
    
    if ss_tot == 0:
        return 1.0
    
    return 1 - (ss_res / ss_tot)


# ============================================================================
# 2. СОЗДАНИЕ ТЕСТОВОЙ СХЕМЫ (СОСТОЯНИЕ БЕЛЛА)
# ============================================================================
def create_bell_circuit():
    """Создаёт схему состояния Белла на 2 кубитах (без измерений)."""
    circuit = QuantumCircuit(2, name="Bell_State")
    circuit.h(0)
    circuit.cx(0, 1)
    return circuit


# ============================================================================
# 3. ГЛОБАЛЬНЫЙ ФОЛДИНГ ДЛЯ УСИЛЕНИЯ ШУМА
# ============================================================================
def fold_circuit_global(circuit, scale_factor):
    """
    Глобальный фолдинг схемы.
    При scale_factor=1 возвращает исходную схему.
    При scale_factor=2: схема + обратная схема.
    При scale_factor=3: схема + обратная + схема.
    """
    if scale_factor <= 1.0:
        return circuit.copy()
    
    k = int(round(scale_factor))
    if k < 1:
        k = 1
    
    folded = circuit.copy()
    
    if k == 1:
        return folded
    
    for i in range(k - 1):
        if i % 2 == 0:
            folded = folded.compose(circuit.inverse(), qubits=range(circuit.num_qubits))
        else:
            folded = folded.compose(circuit, qubits=range(circuit.num_qubits))
    
    return folded


# ============================================================================
# 4. ВЫЧИСЛЕНИЕ ⟨ZZ⟩ ИЗ РЕЗУЛЬТАТОВ ИЗМЕРЕНИЙ
# ============================================================================
def compute_expectation_zz(counts):
    """
    Вычисляет ожидаемое значение ⟨Z₀Z₁⟩ из измерений.
    Формула: ⟨ZZ⟩ = (n00 + n11 - n01 - n10) / total_shots
    """
    n00 = counts.get('00', 0)
    n01 = counts.get('01', 0)
    n10 = counts.get('10', 0)
    n11 = counts.get('11', 0)
    
    total = n00 + n01 + n10 + n11
    if total == 0:
        return 0.0
    
    expectation = (n00 + n11 - n01 - n10) / total
    return expectation


# ============================================================================
# 5. ЗАПУСК СХЕМЫ НА БЭКЕНДЕ И СБОР МЕТРИК
# ============================================================================
def run_and_get_expectation(circuit, backend, shots=SHOTS):
    """
    Запускает схему на указанном бэкенде и возвращает:
    - expectation ⟨ZZ⟩
    - counts
    - execution_time
    """
    # Добавляем измерения в вычислительном базисе
    meas_circuit = circuit.copy()
    meas_circuit.measure_all()
    
    # Транспиляция под бэкенд
    transpiled = transpile(meas_circuit, backend)
    
    # Замер времени
    start_time = time.time()
    
    # Запуск
    result = backend.run(transpiled, shots=shots).result()
    
    execution_time = time.time() - start_time
    
    counts = result.get_counts()
    expectation = compute_expectation_zz(counts)
    
    return {
        'expectation': expectation,
        'counts': counts,
        'execution_time': execution_time
    }


# ============================================================================
# 6. СОЗДАНИЕ ШУМОВОГО БЭКЕНДА
# ============================================================================
def create_noisy_backend():
    """
    Создаёт шумовой симулятор с моделью деполяризации.
    Ошибки: 1% на однокубитные гейты, 2% на двухкубитные гейты.
    """
    noise_model = NoiseModel()
    
    # Деполяризующий шум для однокубитных гейтов (1%)
    depol_1q = depolarizing_error(0.01, 1)
    noise_model.add_all_qubit_quantum_error(depol_1q, ['h', 'x', 'rz', 'sx'])
    
    # Деполяризующий шум для двухкубитных гейтов (2%)
    depol_2q = depolarizing_error(0.02, 2)
    noise_model.add_all_qubit_quantum_error(depol_2q, ['cx', 'cz'])
    
    # Ошибки считывания (2%)
    readout_err = ReadoutError([[0.98, 0.02], [0.02, 0.98]])
    for i in range(2):
        noise_model.add_readout_error(readout_err, [i])
    
    backend = AerSimulator(noise_model=noise_model)
    
    print("✅ Создан шумовой симулятор с моделью деполяризации:")
    print("   - Однокубитные гейты: ошибка 1%")
    print("   - Двухкубитные гейты: ошибка 2%")
    print("   - Считывание: ошибка 2%")
    
    return backend


# ============================================================================
# 7. ОСНОВНАЯ ФУНКЦИЯ БЕНЧМАРКА
# ============================================================================
def run_benchmark():
    """Запускает бенчмарк ZNE и выводит таблицу результатов."""
    
    print("\n" + "="*70)
    print("БЕНЧМАРК ZNE: СОСТОЯНИЕ БЕЛЛА |Φ⁺⟩")
    print("="*70)
    
    # Создание шумового бэкенда
    backend = create_noisy_backend()
    print(f"\nShots: {SHOTS}")
    print(f"Scale factors: {SCALE_FACTORS}")
    
    # Идеальное значение
    IDEAL = 1.0
    
    # Создание исходной схемы
    base_circuit = create_bell_circuit()
    print(f"\nИсходная схема (глубина: {base_circuit.depth()} гейтов):")
    print(base_circuit.draw(output='text'))
    
    # ========================================================================
    # ЗАПУСК БЕЗ МИТИГАЦИИ (BASELINE)
    # ========================================================================
    print("\n" + "-"*70)
    print("ЗАПУСК 1: БЕЗ МИТИГАЦИИ (scale_factor = 1)")
    print("-"*70)
    
    baseline_result = run_and_get_expectation(base_circuit, backend, shots=SHOTS)
    baseline_error = abs(IDEAL - baseline_result['expectation'])
    
    print(f"   ⟨ZZ⟩_raw   = {baseline_result['expectation']:.6f}")
    print(f"   Ошибка     = {baseline_error:.6f}")
    print(f"   exec_time  = {baseline_result['execution_time']:.2f} с")
    print(f"   counts     = {baseline_result['counts']}")
    
    # ========================================================================
    # ЗАПУСК С ZNE ДЛЯ КАЖДОГО SCALE_FACTOR
    # ========================================================================
    print("\n" + "-"*70)
    print("ЗАПУСК 2: ZNE МИТИГАЦИЯ (global folding)")
    print("-"*70)
    
    zne_results = []
    for sf in SCALE_FACTORS:
        print(f"\n  Scale factor = {sf}...")
        
        # Создание растянутой схемы
        folded = fold_circuit_global(base_circuit, sf)
        print(f"    Глубина схемы: {folded.depth()} гейтов")
        
        # Запуск
        result = run_and_get_expectation(folded, backend, shots=SHOTS)
        zne_results.append(result)
        
        print(f"    ⟨ZZ⟩_noisy = {result['expectation']:.6f}")
        print(f"    exec_time  = {result['execution_time']:.2f} с")
    
    # ========================================================================
    # ЭКСТРАПОЛЯЦИЯ К НУЛЕВОМУ ШУМУ
    # ========================================================================
    noisy_values = [r['expectation'] for r in zne_results]
    
    # Линейная экстраполяция
    coeffs_linear = np.polyfit(SCALE_FACTORS, noisy_values, deg=1)
    mitigated_linear = coeffs_linear[1]  # значение при scale=0
    error_linear = abs(IDEAL - mitigated_linear)
    
    # Полиномиальная экстраполяция (2-й степени)
    coeffs_poly = np.polyfit(SCALE_FACTORS, noisy_values, deg=2)
    mitigated_poly = coeffs_poly[2]  # значение при scale=0
    error_poly = abs(IDEAL - mitigated_poly)
    
    # R² для линейной модели (без sklearn)
    linear_pred = [coeffs_linear[0]*sf + coeffs_linear[1] for sf in SCALE_FACTORS]
    r2 = r2_score_manual(noisy_values, linear_pred)
    
    # ========================================================================
    # РАСЧЁТ СТОИМОСТИ
    # ========================================================================
    total_jobs = 1 + len(SCALE_FACTORS)
    total_exec_time = baseline_result['execution_time'] + sum(r['execution_time'] for r in zne_results)
    
    # ========================================================================
    # ВЫВОД ИТОГОВОЙ ТАБЛИЦЫ
    # ========================================================================
    print("\n" + "="*70)
    print("РЕЗУЛЬТАТЫ ЭКСТРАПОЛЯЦИИ")
    print("="*70)
    print(f"Линейная модель:     ⟨ZZ⟩ = {mitigated_linear:.6f}, ошибка = {error_linear:.6f}")
    print(f"Полиномиальная (deg2): ⟨ZZ⟩ = {mitigated_poly:.6f}, ошибка = {error_poly:.6f}")
    print(f"R² линейной регрессии: {r2:.6f}")
    
    improvement = baseline_error / error_linear if error_linear > 0 else float('inf')
    
    print("\n" + "="*70)
    print("ИТОГОВАЯ ТАБЛИЦА")
    print("="*70)
    
    print(f"\n{'Метод':<25} {'⟨ZZ⟩':<12} {'Ошибка':<12} {'Запусков':<10} {'Exec (с)':<10}")
    print("-" * 70)
    print(f"{'Без митигации':<25} {baseline_result['expectation']:<12.6f} {baseline_error:<12.6f} {1:<10} {baseline_result['execution_time']:<10.2f}")
    print(f"{'ZNE (linear)':<25} {mitigated_linear:<12.6f} {error_linear:<12.6f} {total_jobs:<10} {total_exec_time:<10.2f}")
    
    print("\n" + "-"*70)
    print("ПОДРОБНЫЕ РЕЗУЛЬТАТЫ ZNE")
    print("-"*70)
    print(f"{'scale_factor':<15} {'⟨ZZ⟩':<12} {'exec_time (с)':<12}")
    print("-" * 40)
    for sf, res in zip(SCALE_FACTORS, zne_results):
        print(f"{sf:<15} {res['expectation']:<12.6f} {res['execution_time']:<12.2f}")
    
    # ========================================================================
    # ПОСТРОЕНИЕ ГРАФИКОВ
    # ========================================================================
    # Точки для построения гладкой кривой экстраполяции
    x_smooth = np.linspace(0, max(SCALE_FACTORS) + 0.5, 100)
    y_linear = coeffs_linear[0] * x_smooth + coeffs_linear[1]
    y_poly = coeffs_poly[0] * x_smooth**2 + coeffs_poly[1] * x_smooth + coeffs_poly[2]
    
    # Создание фигуры с двумя подграфиками
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
    
    # График 1: Экстраполяция
    ax1.scatter(SCALE_FACTORS, noisy_values, color='red', s=100, zorder=5, label='Измеренные значения')
    ax1.plot(x_smooth, y_linear, 'b--', linewidth=2, label=f'Линейная экстраполяция (R²={r2:.4f})')
    ax1.plot(x_smooth, y_poly, 'g-', linewidth=2, label='Полиномиальная (deg=2)')
    ax1.scatter(0, mitigated_linear, color='blue', s=150, marker='s', zorder=5, label=f'ZNE результат (линейная) = {mitigated_linear:.4f}')
    ax1.scatter(0, mitigated_poly, color='green', s=150, marker='^', zorder=5, label=f'ZNE результат (полином) = {mitigated_poly:.4f}')
    ax1.axhline(y=IDEAL, color='black', linestyle=':', linewidth=2, label=f'Идеальное значение = {IDEAL}')
    ax1.set_xlabel('Scale factor (λ)', fontsize=12)
    ax1.set_ylabel('⟨Z₀Z₁⟩', fontsize=12)
    ax1.set_title('ZNE экстраполяция к нулевому шуму', fontsize=14)
    ax1.legend(loc='best', fontsize=10)
    ax1.grid(True, alpha=0.3)
    ax1.set_xlim(-0.5, max(SCALE_FACTORS) + 0.5)
    ax1.set_ylim(min(noisy_values) - 0.1, 1.05)
    
    # График 2: Сравнение ошибок
    methods = ['Без митигации', 'ZNE (линейная)']
    errors = [baseline_error, error_linear]
    colors_errors = ['red', 'blue']
    
    bars = ax2.bar(methods, errors, color=colors_errors, edgecolor='black', linewidth=1.5)
    ax2.set_ylabel('Абсолютная ошибка', fontsize=12)
    ax2.set_title(f'Улучшение точности: {improvement:.2f}x', fontsize=14)
    ax2.grid(True, alpha=0.3, axis='y')
    
    # Добавление значений на столбцы
    for bar, err in zip(bars, errors):
        ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.002, 
                f'{err:.5f}', ha='center', va='bottom', fontsize=11, fontweight='bold')
    
    plt.tight_layout()
    plt.savefig('zne_benchmark_results.png', dpi=150, bbox_inches='tight')
    plt.show()
    
    print("\n📊 График сохранён как 'zne_benchmark_results.png'")
    
    # ========================================================================
    # ВЫВОДЫ
    # ========================================================================
    print("\n" + "="*70)
    print("ВЫВОДЫ")
    print("="*70)
    print(f"• Точность: улучшение ошибки в {improvement:.2f}x")
    print(f"• Стоимость: {total_jobs}x больше запусков, время увеличено в {total_exec_time / baseline_result['execution_time']:.1f}x")
    print(f"• Устойчивость: R² = {r2:.6f} ({'хорошая линейная зависимость' if r2 > 0.9 else 'слабая линейная зависимость'})")
    
    if r2 < 0.7:
        print("\n⚠️ Низкий R² указывает на немарковский или когерентный шум.")
        print("   Для такого типа шума стандартный ZNE с глобальным фолдингом малоэффективен.")
    else:
        print("\n✅ Высокий R² указывает на марковский шум, хорошо поддающийся экстраполяции.")
    
    # Возвращаем данные для возможного сохранения
    return {
        'baseline': {
            'expectation': baseline_result['expectation'],
            'error': baseline_error,
            'execution_time': baseline_result['execution_time']
        },
        'zne': {
            'mitigated_linear': mitigated_linear,
            'error_linear': error_linear,
            'improvement': improvement,
            'r2': r2,
            'total_jobs': total_jobs,
            'total_exec_time': total_exec_time,
            'noisy_values': noisy_values,
            'scale_factors': SCALE_FACTORS,
            'ideal': IDEAL
        }
    }


# ============================================================================
# 8. ТОЧКА ВХОДА
# ============================================================================
if __name__ == "__main__":
    results = run_benchmark()
    
    # Текстовая визуализация ASCII
    print("\n" + "="*70)
    print("ТЕКСТОВАЯ ВИЗУАЛИЗАЦИЯ ЭКСТРАПОЛЯЦИИ")
    print("="*70)
    
    scale_factors = results['zne']['scale_factors']
    noisy_vals = results['zne']['noisy_values']
    mitigated = results['zne']['mitigated_linear']
    ideal = results['zne']['ideal']
    
    print("\nЗависимость ⟨ZZ⟩ от scale_factor:")
    print("   1.00 ┤")
    for sf, val in zip(scale_factors, noisy_vals):
        bar_len = int(val * 50)
        bar = "█" * bar_len
        print(f"sf={sf}  ┤ {bar} {val:.4f}")
    
    # Экстраполяция
    mitigated_bar = int(mitigated * 50)
    print(f"extrap┤ {'█' * mitigated_bar} {mitigated:.4f} (при sf=0)")
    print(f"ideal ┤ {'█' * 50} {ideal:.4f}")
    
    print("\n✅ БЕНЧМАРК ЗАВЕРШЁН")