Lab2.py - master - LadyChat2015/stepik

Lab2.py
555 строк24 KB
LadyChat2015
upload files
01 мар 2026, 16:53
01 мар 2026, 16:530eb0088
"""
Лабораторная работа: Валидация алгоритма Гровера на трёх платформах
- Идеальный симулятор (AerSimulator без шума)
- Симулятор с моделями шума (fake_backend)
- Реальный квантовый процессор (IBM Quantum)
"""

import numpy as np
from qiskit import QuantumCircuit, transpile
from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler
from qiskit.visualization import plot_histogram, plot_distribution
import matplotlib.pyplot as plt
from qiskit_ibm_runtime.fake_provider import FakeManila, FakeManilaV2
from qiskit.circuit.library import GroverOperator, MCMT, ZGate
from qiskit.circuit import QuantumRegister, ClassicalRegister, AncillaRegister
import time
import os
from dotenv import load_dotenv

# Загружаем переменные окружения
load_dotenv()


class GroverAlgorithmTester:
    """
    Класс для тестирования алгоритма Гровера на различных платформах
    """

    def __init__(self, n_qubits=3, marked_states=['101']):
        """
        Инициализация тестера

        Args:
            n_qubits: количество кубитов
            marked_states: список помеченных состояний (в двоичной строке)
        """
        self.n_qubits = n_qubits
        self.marked_states = marked_states
        self.circuit = None

        print(f"\n{'=' * 60}")
        print(f"АЛГОРИТМ ГРОВЕРА - ТЕСТИРОВАНИЕ НА {n_qubits} КУБИТАХ")
        print(f"Помеченные состояния: {marked_states}")
        print(f"{'=' * 60}")

    def create_oracle(self):
        """
        Создаёт оракул для помеченных состояний
        """
        oracle = QuantumCircuit(self.n_qubits)

        for marked_state in self.marked_states:
            # Применяем X гейты к нулевым битам в помеченном состоянии
            for i, bit in enumerate(reversed(marked_state)):
                if bit == '0':
                    oracle.x(i)

            # Мульти-контролируемый Z гейт
            oracle.h(self.n_qubits - 1)
            oracle.mcx(list(range(self.n_qubits - 1)), self.n_qubits - 1)
            oracle.h(self.n_qubits - 1)

            # Возвращаем X гейты обратно
            for i, bit in enumerate(reversed(marked_state)):
                if bit == '0':
                    oracle.x(i)

        return oracle.to_gate(label='Oracle')

    def create_diffusion(self):
        """
        Создаёт оператор диффузии Гровера
        """
        diffusion = QuantumCircuit(self.n_qubits)

        # H на всех кубитах
        diffusion.h(range(self.n_qubits))

        # X на всех кубитах
        diffusion.x(range(self.n_qubits))

        # Мульти-контролируемый Z
        diffusion.h(self.n_qubits - 1)
        diffusion.mcx(list(range(self.n_qubits - 1)), self.n_qubits - 1)
        diffusion.h(self.n_qubits - 1)

        # X на всех кубитах
        diffusion.x(range(self.n_qubits))

        # H на всех кубитах
        diffusion.h(range(self.n_qubits))

        return diffusion.to_gate(label='Diffusion')

    def build_circuit(self, iterations=None):
        """
        Строит полную схему алгоритма Гровера

        Args:
            iterations: количество итераций Гровера (если None, вычисляется оптимально)
        """
        # Вычисляем оптимальное количество итераций
        if iterations is None:
            N = 2 ** self.n_qubits
            M = len(self.marked_states)
            # Оптимальное число итераций ≈ (π/4) * sqrt(N/M)
            iterations = int(np.pi / 4 * np.sqrt(N / M))
            print(f"\nОптимальное число итераций: {iterations}")

        # Создаём схему
        circuit = QuantumCircuit(self.n_qubits, self.n_qubits)

        # Инициализация: H на всех кубитах
        circuit.h(range(self.n_qubits))
        circuit.barrier()

        # Получаем оракул и оператор диффузии
        oracle = self.create_oracle()
        diffusion = self.create_diffusion()

        # Итерации Гровера
        for i in range(iterations):
            circuit.append(oracle, range(self.n_qubits))
            circuit.barrier()
            circuit.append(diffusion, range(self.n_qubits))
            circuit.barrier()

        # Измерение
        circuit.measure(range(self.n_qubits), range(self.n_qubits))

        self.circuit = circuit
        return circuit

    def run_on_ideal_simulator(self, shots=4096):
        """
        Запуск на идеальном симуляторе (без шума)
        """
        print("\n" + "=" * 60)
        print("1. ИДЕАЛЬНЫЙ СИМУЛЯТОР (AerSimulator без шума)")
        print("=" * 60)

        # Создаём симулятор
        simulator = AerSimulator(method='statevector')

        # Транспилируем схему
        compiled_circuit = transpile(self.circuit, simulator)

        # Запускаем
        start_time = time.time()
        job = simulator.run(compiled_circuit, shots=shots)
        result = job.result()
        end_time = time.time()

        # Получаем результаты
        counts = result.get_counts(compiled_circuit)

        print(f"Время выполнения: {end_time - start_time:.3f} сек")
        print(f"Результаты ({shots} выстрелов):")
        self._print_results(counts)

        # Сохраняем визуализацию
        self._save_histogram(counts, 'ideal_simulator', 'Идеальный симулятор')

        return counts

    def run_on_noisy_simulator(self, shots=4096):
        """
        Запуск на симуляторе с моделью шума (на основе реального устройства)
        """
        print("\n" + "=" * 60)
        print("2. СИМУЛЯТОР С МОДЕЛЬЮ ШУМА (Fake Manila)")
        print("=" * 60)

        # Используем fake backend для моделирования шума
        fake_backend = FakeManilaV2()
        noise_model = NoiseModel.from_backend(fake_backend)

        # Создаём симулятор с шумом
        simulator = AerSimulator(noise_model=noise_model)

        # Транспилируем схему под архитектуру fake backend
        compiled_circuit = transpile(self.circuit, fake_backend)

        # Запускаем
        start_time = time.time()
        job = simulator.run(compiled_circuit, shots=shots)
        result = job.result()
        end_time = time.time()

        # Получаем результаты
        counts = result.get_counts(compiled_circuit)

        print(f"Время выполнения: {end_time - start_time:.3f} сек")
        print(f"Модель шума: {fake_backend.name}")
        print(f"Результаты ({shots} выстрелов):")
        self._print_results(counts)

        # Сохраняем визуализацию
        self._save_histogram(counts, 'noisy_simulator', 'Симулятор с шумом (Fake Manila)')

        return counts

    def run_on_real_device(self, shots=4096):
        """
        Запуск на реальном квантовом устройстве IBM Quantum
        """
        print("\n" + "=" * 60)
        print("3. РЕАЛЬНОЕ КВАНТОВОЕ УСТРОЙСТВО")
        print("=" * 60)

        # Получаем токен
        token = os.getenv('IBM_QUANTUM_TOKEN')
        if not token:
            print("Токен IBM Quantum не найден!")
            token = input("Введите ваш токен IBM Quantum: ").strip()
            if not token:
                print("Пропускаем запуск на реальном устройстве")
                return None

        try:
            # Подключаемся к IBM Quantum
            service = QiskitRuntimeService(channel="ibm_quantum", token=token)

            # Выбираем наименее загруженное устройство
            backend = service.least_busy(simulator=False, operational=True, min_num_qubits=self.n_qubits)
            print(f"Выбрано устройство: {backend.name}")
            print(f"Количество кубитов: {backend.num_qubits}")
            print(f"Вероятность ошибки считывания: {backend.measure_lo}")  # Для информации

            # Транспилируем схему
            compiled_circuit = transpile(self.circuit, backend, optimization_level=3)

            # Создаём сэмплер и запускаем
            sampler = Sampler(backend=backend)

            print(f"Отправка задачи... ID: будет получен")
            start_time = time.time()
            job = sampler.run([compiled_circuit], shots=shots)
            print(f"ID задачи: {job.job_id()}")
            print("Ожидание выполнения... (это может занять несколько минут)")

            # Ждём результат
            result = job.result()
            end_time = time.time()

            # Получаем результаты
            counts = result[0].data.meas.get_counts()

            print(f"Время выполнения (включая ожидание): {end_time - start_time:.3f} сек")
            print(f"Результаты ({shots} выстрелов):")
            self._print_results(counts)

            # Сохраняем визуализацию
            self._save_histogram(counts, 'real_device', f'Реальное устройство ({backend.name})')

            return counts, backend.name

        except Exception as e:
            print(f"Ошибка при запуске на реальном устройстве: {e}")
            return None

    def _print_results(self, counts):
        """
        Печатает результаты в отформатированном виде
        """
        total = sum(counts.values())

        # Сортируем по убыванию
        sorted_counts = sorted(counts.items(), key=lambda x: x[1], reverse=True)

        print(f"{'Состояние':<10} {'Количество':<12} {'Вероятность':<12}")
        print("-" * 34)
        for state, count in sorted_counts[:5]:  # Показываем топ-5
            prob = count / total
            mark = "✓" if state in self.marked_states else " "
            print(f"{state}{mark:<3} {count:<12} {prob:<12.4f}")

        if len(counts) > 5:
            print(f"... и ещё {len(counts) - 5} состояний")

    def _save_histogram(self, counts, filename, title):
        """
        Сохраняет гистограмму результатов
        """
        try:
            # Сортируем для лучшей визуализации
            sorted_states = sorted(counts.keys())

            # Выделяем помеченные состояния
            colors = ['red' if state in self.marked_states else 'blue' for state in sorted_states]

            fig, ax = plt.subplots(figsize=(12, 6))
            plot_histogram(counts, ax=ax, color=colors, bar_labels=False)
            ax.set_title(f"{title}\nПомеченные состояния: {self.marked_states}", fontsize=14)
            ax.set_xlabel("Состояния", fontsize=12)
            ax.set_ylabel("Количество", fontsize=12)

            # Добавляем легенду
            from matplotlib.patches import Patch
            legend_elements = [
                Patch(facecolor='red', label='Помеченное состояние'),
                Patch(facecolor='blue', label='Другое состояние')
            ]
            ax.legend(handles=legend_elements, loc='upper right')

            plt.tight_layout()
            plt.savefig(f'{filename}.png', dpi=150)
            print(f"Гистограмма сохранена в {filename}.png")
            plt.show()
        except Exception as e:
            print(f"Не удалось сохранить гистограмму: {e}")

    def analyze_results(self, ideal_counts, noisy_counts, real_result):
        """
        Сравнительный анализ результатов с трёх платформ
        """
        print("\n" + "=" * 60)
        print("СРАВНИТЕЛЬНЫЙ АНАЛИЗ РЕЗУЛЬТАТОВ")
        print("=" * 60)

        if real_result:
            real_counts, real_backend = real_result
        else:
            real_counts = None
            real_backend = "Недоступно"

        # Нормализуем вероятности
        def get_probs(counts):
            total = sum(counts.values())
            return {k: v / total for k, v in counts.items()}

        ideal_probs = get_probs(ideal_counts)
        noisy_probs = get_probs(noisy_counts) if noisy_counts else {}
        real_probs = get_probs(real_counts) if real_counts else {}

        # Собираем все возможные состояния
        all_states = set(ideal_probs.keys())
        if noisy_counts:
            all_states |= set(noisy_counts.keys())
        if real_counts:
            all_states |= set(real_counts.keys())

        all_states = sorted(all_states)

        # Создаём таблицу сравнения
        print(f"\n{'Состояние':<12} {'Идеал':<10} {'Шумный':<10} {'Реальный':<10} {'Маркер':<8}")
        print("-" * 50)

        # Статистика по помеченным состояниям
        marked_probs_ideal = 0
        marked_probs_noisy = 0
        marked_probs_real = 0

        for state in all_states:
            ideal_prob = ideal_probs.get(state, 0)
            noisy_prob = noisy_probs.get(state, 0) if noisy_counts else 0
            real_prob = real_probs.get(state, 0) if real_counts else 0

            is_marked = state in self.marked_states
            marker = "✓" if is_marked else " "

            print(f"{state:<12} {ideal_prob:<10.4f} {noisy_prob:<10.4f} {real_prob:<10.4f} {marker:<8}")

            if is_marked:
                marked_probs_ideal += ideal_prob
                marked_probs_noisy += noisy_prob
                marked_probs_real += real_prob

        print("\n" + "=" * 60)
        print("КЛЮЧЕВЫЕ МЕТРИКЫ")
        print("=" * 60)

        # Суммарная вероятность помеченных состояний
        print(f"\n1. Суммарная вероятность помеченных состояний:")
        print(f"   - Идеальный симулятор: {marked_probs_ideal:.4f} (должно быть близко к 1.0)")
        print(f"   - Шумный симулятор:    {marked_probs_noisy:.4f}")
        if real_counts:
            print(f"   - Реальное устройство:  {marked_probs_real:.4f}")

        # Соотношение сигнал/шум
        if real_counts:
            print(f"\n2. Соотношение сигнал/шум (реальное устройство):")
            total_prob_marked = marked_probs_real
            total_prob_noise = 1 - total_prob_marked
            print(f"   - Сигнал (помеченные): {total_prob_marked:.4f}")
            print(f"   - Шум (остальные):    {total_prob_noise:.4f}")
            print(f"   - SNR:                 {total_prob_marked / total_prob_noise:.4f}")

        # Анализ ошибок
        print("\n3. АНАЛИЗ ИСТОЧНИКОВ ОТКЛОНЕНИЙ:")

        # Расхождение между идеальным и шумным симулятором
        mse_noise = sum((ideal_probs.get(s, 0) - noisy_probs.get(s, 0)) ** 2 for s in all_states) / len(all_states)
        print(f"\n   a) Влияние модели шума (MSE идеал vs шумный): {mse_noise:.6f}")

        if real_counts:
            mse_real = sum((ideal_probs.get(s, 0) - real_probs.get(s, 0)) ** 2 for s in all_states) / len(all_states)
            print(f"   b) Реальное устройство vs идеал (MSE): {mse_real:.6f}")

            # Сравнение шумного симулятора с реальным устройством
            mse_model_vs_real = sum((noisy_probs.get(s, 0) - real_probs.get(s, 0)) ** 2 for s in all_states) / len(
                all_states)
            print(f"   c) Точность модели шума (MSE шумный vs реальный): {mse_model_vs_real:.6f}")

        # Диагностика ошибок
        print("\n4. ДИАГНОСТИКА ИСТОЧНИКОВ ОШИБОК:")

        if real_counts:
            # Анализ нежелательных состояний
            unwanted_states = [s for s in all_states if s not in self.marked_states]

            # Состояния с одной ошибкой (ближайшие по Хэммингову расстоянию)
            one_bit_errors = []
            for marked in self.marked_states:
                for unwanted in unwanted_states:
                    # Считаем расстояние Хэмминга
                    distance = sum(b1 != b2 for b1, b2 in zip(marked, unwanted))
                    if distance == 1:
                        one_bit_errors.append(unwanted)

            one_bit_error_prob = sum(real_probs.get(s, 0) for s in one_bit_errors)
            other_error_prob = sum(real_probs.get(s, 0) for s in unwanted_states if s not in one_bit_errors)

            print(f"\n   - Вероятность одно-битовых ошибок: {one_bit_error_prob:.4f}")
            print(f"   - Вероятность прочих ошибок: {other_error_prob:.4f}")

            # Интерпретация
            if one_bit_error_prob > other_error_prob * 2:
                print("\n   ВЫВОД: Доминируют одно-битовые ошибки → проблемы с декогеренцией")
            elif other_error_prob > one_bit_error_prob * 2:
                print("\n   ВЫВОД: Доминируют много-битовые ошибки → проблемы с гейтами")
            else:
                print("\n   ВЫВОД: Смешанный тип ошибок → комбинация факторов")

        print("\n" + "=" * 60)
        print("ИТОГОВОЕ ЗАКЛЮЧЕНИЕ")
        print("=" * 60)

        if real_counts:
            success_threshold = 0.7
            if marked_probs_real >= success_threshold:
                print(f"\n✅ Алгоритм успешно выполнен на реальном устройстве!")
                print(f"   Вероятность успеха {marked_probs_real:.2%} > {success_threshold:.0%}")
            else:
                print(f"\n⚠️ Алгоритм выполнен с низкой вероятностью успеха {marked_probs_real:.2%}")
                print(f"   Требуется оптимизация или использование устройств с меньшим шумом")

            # Сравнение с симуляцией
            if abs(marked_probs_real - marked_probs_noisy) < 0.1:
                print(f"\n   Модель шума хорошо предсказывает поведение реального устройства")
            else:
                print(f"\n   Модель шума не точно предсказывает поведение реального устройства")
                print(f"   Разница: {abs(marked_probs_real - marked_probs_noisy):.2%}")

        print("\n" + "=" * 60)


def test_with_different_parameters():
    """
    Тестирование с разными параметрами
    """
    print("\n" + "=" * 60)
    print("ТЕСТИРОВАНИЕ С РАЗНЫМИ ПАРАМЕТРАМИ")
    print("=" * 60)

    test_cases = [
        (2, ['11'], "2 кубита, одно помеченное состояние"),
        (3, ['101'], "3 кубита, одно помеченное состояние"),
        (3, ['101', '110'], "3 кубита, два помеченных состояния"),
        (4, ['1010'], "4 кубита, одно помеченное состояние")
    ]

    for n_qubits, marked_states, description in test_cases:
        print(f"\n▶ Тест: {description}")

        tester = GroverAlgorithmTester(n_qubits=n_qubits, marked_states=marked_states)
        circuit = tester.build_circuit()

        # Показываем схему
        print(f"\nРазмер схемы: {circuit.depth()} глубина, {circuit.size()} операций")

        # Запускаем на идеальном симуляторе
        ideal_counts = tester.run_on_ideal_simulator(shots=2048)

        # Пауза для визуального разделения тестов
        print("\n" + "-" * 40)


def main():
    """
    Основная функция
    """
    print("=" * 60)
    print("ЛАБОРАТОРНАЯ РАБОТА: ВАЛИДАЦИЯ АЛГОРИТМА ГРОВЕРА")
    print("=" * 60)

    # Выбор режима работы
    print("\nВыберите режим работы:")
    print("1. Полное тестирование (все три платформы)")
    print("2. Только симуляторы (идеальный + шумный)")
    print("3. Тестирование с разными параметрами")

    choice = input("\nВаш выбор (1/2/3): ").strip()

    if choice == '3':
        test_with_different_parameters()
        return

    # Основная конфигурация
    n_qubits = 3
    marked_states = ['101']  # Можно изменить

    # Создаём тестер
    tester = GroverAlgorithmTester(n_qubits=n_qubits, marked_states=marked_states)

    # Строим схему
    circuit = tester.build_circuit()

    # Показываем схему
    print("\nКВАНТОВАЯ СХЕМА:")
    print(circuit.draw(output='text'))

    try:
        circuit.draw('mpl')
        plt.savefig('grover_circuit.png')
        print("Схема сохранена в 'grover_circuit.png'")
    except:
        pass

    # Запуск на идеальном симуляторе
    ideal_counts = tester.run_on_ideal_simulator(shots=4096)

    # Запуск на шумном симуляторе
    noisy_counts = tester.run_on_noisy_simulator(shots=4096)

    # Запуск на реальном устройстве (если выбрано)
    real_result = None
    if choice == '1':
        print("\n" + "=" * 60)
        response = input("Запустить на реальном устройстве IBM Quantum? (да/нет): ").strip().lower()
        if response in ['да', 'yes', 'y', 'д']:
            real_result = tester.run_on_real_device(shots=4096)

    # Анализ результатов
    tester.analyze_results(ideal_counts, noisy_counts, real_result)

    print("\n" + "=" * 60)
    print("ЛАБОРАТОРНАЯ РАБОТА ЗАВЕРШЕНА")
    print("=" * 60)


if __name__ == "__main__":
    main()