Task2_3/Task2_3.py - master - timshk/Shalaev_Stepik

Task2_3/Task2_3.py
365 строк15 KB
timshk
Stepik задания
01 июн 2026, 18:58
01 июн 2026, 18:582e75834
"""
ВЕБ-ПРИЛОЖЕНИЕ ДЛЯ ГИБРИДНОЙ ОПТИМИЗАЦИИ С QAOA
Задача: Max-Cut на графе из 4 вершин (учебный прототип)
Фреймворк: Streamlit + Qiskit Aer
"""

import streamlit as st
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.quantum_info import SparsePauliOp, Pauli

# ============================================================================
# 1. НАСТРОЙКА СТРАНИЦЫ И СТИЛЯ
# ============================================================================
st.set_page_config(
    page_title="QAOA Max-Cut Optimizer",
    page_icon="✂️",
    layout="wide",
    initial_sidebar_state="collapsed"
)

st.title("✂️ QAOA Max-Cut Optimizer")
st.markdown("Гибридный квантово-классический алгоритм для задачи максимального разреза графа")
st.divider()

# ============================================================================
# 2. БОКОВАЯ ПАНЕЛЬ: ВВОД ДАННЫХ
# ============================================================================
with st.sidebar:
    st.header("⚙️ Параметры задачи")
    
    # Ввод матрицы смежности (4x4)
    st.subheader("Матрица смежности (4x4)")
    st.caption("Введите веса рёбер (0 = нет ребра)")
    
    # Создание матрицы через виджет number_input для каждого элемента
    cols = st.columns(4)
    matrix = np.zeros((4, 4))
    
    for i in range(4):
        for j in range(4):
            if i <= j:
                # Только верхний треугольник (симметричная матрица)
                with cols[j]:
                    val = st.number_input(
                        f"({i+1},{j+1})",
                        min_value=0.0,
                        max_value=10.0,
                        value=1.0 if (i != j and (i,j) in [(0,1), (1,2), (2,3), (0,3)]) else 0.0,
                        step=0.5,
                        key=f"m_{i}_{j}",
                        label_visibility="collapsed"
                    )
                    matrix[i, j] = val
                    matrix[j, i] = val
            else:
                # Пропускаем нижний треугольник (уже заполнен)
                pass
    
    st.divider()
    
    # Слайдер для количества слоёв p
    p = st.slider("Количество слоёв QAOA (p)", min_value=1, max_value=5, value=1, help="Больше слоёв = точнее, но сложнее схема")
    
    # Кнопка запуска
    run_button = st.button("🚀 Запустить оптимизацию", type="primary", use_container_width=True)

# ============================================================================
# 3. ОТОБРАЖЕНИЕ ГРАФА
# ============================================================================
col1, col2 = st.columns([1, 1])

with col1:
    st.subheader("📊 Исходный граф")
    
    # Построение графа для визуализации
    G = nx.Graph()
    for i in range(4):
        G.add_node(i, label=f"Вершина {i+1}")
    
    edges = []
    for i in range(4):
        for j in range(i+1, 4):
            if matrix[i, j] > 0:
                G.add_edge(i, j, weight=matrix[i, j])
                edges.append((i, j, matrix[i, j]))
    
    # Рисуем граф
    fig, ax = plt.subplots(figsize=(4, 4))
    pos = nx.circular_layout(G)
    nx.draw_networkx_nodes(G, pos, node_size=500, node_color='lightblue', ax=ax)
    nx.draw_networkx_labels(G, pos, labels={i: f"{i+1}" for i in range(4)}, ax=ax, font_size=12, font_weight='bold')
    
    # Рёбра с весами
    nx.draw_networkx_edges(G, pos, width=2, edge_color='gray', ax=ax)
    edge_labels = {(i, j): f"{w:.1f}" for (i, j, w) in edges}
    nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, font_size=10, ax=ax)
    
    ax.set_title("Граф задачи Max-Cut")
    ax.axis('off')
    st.pyplot(fig)

# ============================================================================
# 4. ОПРЕДЕЛЕНИЕ BACKEND-ФУНКЦИЙ
# ============================================================================
def maxcut_hamiltonian(matrix):
    """
    Преобразует матрицу смежности в гамильтониан Max-Cut.
    
    Для Max-Cut: H = Σ_{i<j} w_ij * (1 - Z_i Z_j) / 2
    Для QAOA минимизируем ⟨H⟩ → максимизируем вес разреза.
    """
    n = matrix.shape[0]
    pauli_list = []
    coeffs = []
    
    for i in range(n):
        for j in range(i+1, n):
            w = matrix[i, j]
            if w != 0:
                # Член (1 - Z_i Z_j) / 2
                # Константа 1/2 даст сдвиг энергии, не влияющий на оптимизацию
                # Добавляем -w/2 * Z_i Z_j
                pauli_str = ['I'] * n
                pauli_str[i] = 'Z'
                pauli_str[j] = 'Z'
                pauli_list.append(''.join(pauli_str))
                coeffs.append(-w / 2)
    
    return SparsePauliOp(pauli_list, coeffs)


def build_qaoa_circuit(hamiltonian, gamma, beta):
    """
    Строит параметризованную схему QAOA для одного слоя (p=1).
    
    Схема:
    1. H⊗n для создания суперпозиции
    2. Гамильтонова эволюция: exp(-iγ H) = произведение exp(-iγ * Z_i Z_j)
    3. Смешивающий член: exp(-iβ Σ X_i) = произведение Rx(2β)
    """
    n = hamiltonian.num_qubits
    circuit = QuantumCircuit(n, n)
    
    # Суперпозиция
    circuit.h(range(n))
    
    # Гамильтонова эволюция exp(-iγ H)
    # Для Max-Cut: H состоит из членов Z_i Z_j
    for pauli, coeff in zip(hamiltonian.paulis, hamiltonian.coeffs):
        # coeff = -w/2, но для exp(iγ H) нужен γ * coeff
        # Преобразуем в exp(iγ * (w/2) * Z_i Z_j) сделаем через Rzz
        # Стандартная схема: Rzz(θ) = exp(-iθ/2 * Z_i Z_j)
        # Нам нужно exp(iγ * w/2 * Z_i Z_j) → θ = -γ * w
        if coeff.imag != 0:
            continue
        w_effective = -2 * coeff.real  # восстанавливаем вес
        theta = -gamma * w_effective
        
        # Находим индексы кубитов, где стоит Z
        indices = [i for i, p in enumerate(pauli.to_label()) if p == 'Z']
        if len(indices) == 2:
            i, j = indices
            circuit.rzz(theta, i, j)
    
    # Смешивающий член exp(-iβ Σ X_i) → Rx(2β)
    for i in range(n):
        circuit.rx(2 * beta, i)
    
    # Измерение
    circuit.measure(range(n), range(n))
    
    return circuit


def evaluate_energy(hamiltonian, counts, shots):
    """
    Вычисляет ожидаемое значение ⟨H⟩ из результатов измерений.
    """
    energy = 0.0
    n = hamiltonian.num_qubits
    
    for bitstring, count in counts.items():
        # Преобразуем битовую строку в массив ±1 (little-endian как в Qiskit)
        bits = np.array([1 if c == '1' else -1 for c in bitstring[::-1]])
        
        # Вычисляем энергию: Σ coeff * Π (Z)
        for pauli, coeff in zip(hamiltonian.paulis, hamiltonian.coeffs):
            term_val = 1.0
            for i, p in enumerate(pauli.to_label()):
                if p == 'Z':
                    term_val *= bits[i]
            energy += coeff * term_val * count
    
    return energy / shots


def run_qaoa(matrix, p_layers, shots=1024, max_iter=50):
    """
    Запускает гибридный цикл QAOA.
    
    Args:
        matrix: матрица смежности (4x4)
        p_layers: количество слоёв QAOA (p)
        shots: количество измерений
        max_iter: максимальное число итераций оптимизатора
    
    Returns:
        optimal_bitstring: лучшая найденная битовая строка
        optimal_energy: соответствующее значение энергии
        history: список энергий по итерациям
    """
    # Получаем гамильтониан
    hamiltonian = maxcut_hamiltonian(matrix)
    n_qubits = hamiltonian.num_qubits
    
    # Инициализация симулятора
    simulator = AerSimulator()
    
    # Начальные параметры (γ и β для каждого слоя)
    # Для p=1: gamma и beta — скаляры
    initial_params = np.random.uniform(0, np.pi, 2 * p_layers)
    
    # История для визуализации
    history = []
    
    # Целевая функция для оптимизатора
    def objective(params):
        nonlocal history
        
        # Разделяем параметры на гаммы и беты
        gammas = params[:p_layers]
        betas = params[p_layers:]
        
        # Для p=1 просто берём gamma и beta
        gamma = gammas[0] if p_layers == 1 else gammas
        beta = betas[0] if p_layers == 1 else betas
        
        # Строим схему
        circuit = build_qaoa_circuit(hamiltonian, gamma, beta)
        
        # Транспиляция и запуск
        pass_manager = generate_preset_pass_manager(optimization_level=1, backend=simulator)
        compiled = pass_manager.run(circuit)
        result = simulator.run(compiled, shots=shots).result()
        counts = result.get_counts()
        
        # Вычисляем энергию
        energy = evaluate_energy(hamiltonian, counts, shots)
        history.append(energy)
        
        return energy
    
    # Запуск оптимизации
    with st.spinner("🔄 Выполняется гибридный цикл QAOA..."):
        result = minimize(
            objective,
            initial_params,
            method='COBYLA',
            options={'maxiter': max_iter, 'disp': False}
        )
    
    # После оптимизации получаем лучшее решение
    optimal_params = result.x
    gammas = optimal_params[:p_layers]
    betas = optimal_params[p_layers:]
    
    # Финальный запуск с оптимальными параметрами
    gamma_final = gammas[0] if p_layers == 1 else gammas
    beta_final = betas[0] if p_layers == 1 else betas
    final_circuit = build_qaoa_circuit(hamiltonian, gamma_final, beta_final)
    
    pass_manager = generate_preset_pass_manager(optimization_level=1, backend=simulator)
    compiled = pass_manager.run(final_circuit)
    final_result = simulator.run(compiled, shots=1024).result()
    final_counts = final_result.get_counts()
    
    # Находим битовую строку с минимальной энергией
    min_energy = float('inf')
    optimal_bitstring = None
    for bitstring, count in final_counts.items():
        bits = np.array([1 if c == '1' else -1 for c in bitstring[::-1]])
        energy = 0.0
        for pauli, coeff in zip(hamiltonian.paulis, hamiltonian.coeffs):
            term_val = 1.0
            for i, p in enumerate(pauli.to_label()):
                if p == 'Z':
                    term_val *= bits[i]
            energy += coeff * term_val
        if energy < min_energy:
            min_energy = energy
            optimal_bitstring = bitstring
    
    return optimal_bitstring, min_energy, history


# ============================================================================
# 5. ЗАПУСК ОПТИМИЗАЦИИ
# ============================================================================
if run_button:
    with col2:
        st.subheader("🔬 Результаты оптимизации")
        
        # Запуск QAOA
        optimal_bitstring, energy, history = run_qaoa(matrix, p, shots=1024, max_iter=30)
        
        # Декодирование решения
        # В Qiskit битовая строка: первый символ — последний кубит (little-endian)
        # Для наглядности переворачиваем
        bitstring_reversed = optimal_bitstring[::-1]
        partition = [int(b) for b in bitstring_reversed]
        
        # Вычисляем вес разреза
        cut_weight = 0
        for i in range(4):
            for j in range(i+1, 4):
                if matrix[i, j] > 0 and partition[i] != partition[j]:
                    cut_weight += matrix[i, j]
        
        # Отображение результатов
        st.metric("Вес разреза (Max-Cut)", f"{cut_weight:.2f}")
        st.metric("Полученное разбиение", f"Группа 0: {[i+1 for i, p in enumerate(partition) if p == 0]}, Группа 1: {[i+1 for i, p in enumerate(partition) if p == 1]}")
        
        # Визуализация графа с цветами по решению
        st.subheader("🎨 Граф с разбиением")
        fig2, ax2 = plt.subplots(figsize=(4, 4))
        
        colors = ['lightcoral' if partition[i] == 0 else 'lightgreen' for i in range(4)]
        nx.draw_networkx_nodes(G, pos, node_size=500, node_color=colors, ax=ax2)
        nx.draw_networkx_labels(G, pos, labels={i: f"{i+1}" for i in range(4)}, ax=ax2, font_size=12, font_weight='bold')
        nx.draw_networkx_edges(G, pos, width=2, edge_color='gray', ax=ax2)
        nx.draw_networkx_edge_labels(G, pos, edge_labels={(i, j): f"{w:.1f}" for (i, j, w) in edges}, font_size=10, ax=ax2)
        
        ax2.set_title(f"Решение Max-Cut (вес = {cut_weight:.2f})")
        ax2.axis('off')
        st.pyplot(fig2)
        
        # График сходимости
        st.subheader("📈 Сходимость оптимизатора")
        fig3, ax3 = plt.subplots(figsize=(8, 4))
        ax3.plot(range(1, len(history)+1), history, 'b-o', linewidth=2, markersize=4)
        ax3.set_xlabel('Итерация')
        ax3.set_ylabel('⟨H⟩ (энергия)')
        ax3.set_title('Сходимость классического оптимизатора COBYLA')
        ax3.grid(True, alpha=0.3)
        st.pyplot(fig3)
        
        # Дополнительная информация
        with st.expander("ℹ️ Подробная информация"):
            st.write(f"**Количество слоёв QAOA (p):** {p}")
            st.write(f"**Оптимальная битовая строка (little-endian):** `{optimal_bitstring}`")
            st.write(f"**Оптимальная битовая строка (big-endian):** `{bitstring_reversed}`")
            st.write(f"**Минимальная энергия⟨H⟩:** {energy:.4f}")
            st.write(f"**Количество итераций оптимизатора:** {len(history)}")
            st.write(f"**Матрица смежности:**")
            st.dataframe(matrix.astype(float))

else:
    with col2:
        st.info("👈 Настройте параметры в боковой панели и нажмите 'Запустить оптимизацию'")