Algorytm symulowanej bifurkacjiΒΆ
Jest to heurystyczny algorytm zainspirowany fizykΔ , zaproponowany przez Goto et al., w 2019 roku. Jak pokazujΔ ostatnie badania, jest to SOTA algorytm do rozwiΔ zywania problemu Isinga (Pawlowski 2025).
Intuicyjnie, moΕΌemy o nim myΕleΔ jako zakodowaniu naszego problemu uoptymalizacjyjnego jako pewien ukΕad dynamiczny opisany rΓ³waniami rΓ³ΕΌniczkowymi. Te rΓ³wnania wywodzΔ siΔ z rΓ³wniaΕ Hamiltona z mechaniki klasycznej. Intuicuyjnie opisujΔ ewolucjΔ poΕoΕΌenia $x$ i pΔdu $y$ pewnej grupy czΔ steczek.
$$ \dot{y}_i = -[a_0 - a(t)]x_i + c_0f_i $$ $$ \dot{x}_i = a_0 y_i $$
gdzie kropka opisuje pochodnΔ po czasie, $a_0$ jest staΕΔ , $a(t)$ jest parametrem bifurkacji rosnΔ cym w czasie $t$ od 0 do $a_0$, a $c_0$ jest parametrem normalizacyjnym zaleΕΌym od siΕ $f_i$. SiΕy $f_i$ sΔ de facto gradientami energi. Co jest waΕΌne, statnie stacjonarnym, po zbinaryzowaniu te rΓ³wnania opisujΔ lokalne minimum modelu Isinga.
RozwiΔ zujΔ c te rΓ³wnania rΓ³ΕΌniczkowe numerycznie (np. metodΔ symplektycznΔ Eulera) otrzymujemy rozwiΔ zanie oryginalnego problemu. Dla naszych zastosowaΕ moΕΌemy go traktowaΔ jako czarnΔ skrzynkΔ.
ZostanΔ Pokazane dwie wersje tego algorytmu ktΓ³re nieznacznie siΔ od siebie rΓ³ΕΌniΔ , balistyczna symulowana bifurkacja i dyskretna symulowana bifurkacja. Te algorytmy w bardzo naturalny sposΓ³b siΔ zrΓ³wnoleglajΔ .
Balistyczna Symulowana BifurkacjaΒΆ
W kaΕΌdym kroku $t_k$ "pozycja" $x(t_k)$ i "pΔd" $y(t_k)$ sΔ aktualizowane za pomocΔ nastΔpujΔ cych wzorΓ³w:
$$ y_i(t_{k+1}) = y_i(t_k) + \{-[a_0 - a(t_k)]x_i(t_k) + c_0 \sum_{j=1}^N J_{i,j}x_j(t_k) \}\Delta_t $$ $$ x_i(t_{k+1}) = x_i(t_k) + a_0 y_i(t_{k+1})\Delta_t $$
Oznaczenia sΔ identyczne jak wczeΕniej. Przez $\Delta_t$ oznaczamy dΕugoΕΔ kroku czasowego. Ponadto przyjmujemy tzw. "nieelastycznΔ scianΔ" dla $\vert x \vert = 1$. Oznacza to, ΕΌe jeΕΌeli w dowolnym momΔcie algorytmu $\vert x \vert > 1$ to obcinamy $x$ do odpowiedniej wartoΕci (-1 lub 1) i ustawiamy $y = 0$. MoΕΌna o tym myΕleΔ jak o czΔ steczce uderzajΔ cej w Εcianie i siΔ do niej przyklejajΔ cΔ . Matematycznie moΕΌemy tΔ operacje wyraziΔ funkcjΔ $\text{Wall}(x, y)$.
$$ \text{Wall}(x, y) = \left\{ \begin{array}{ll} (\text{sign}(x), 0), & \text{jeΕΌeli } \vert x \vert > 1 \\ (x, y) & \text{w pozostaΕych przypadkach}\\ \end{array} \right. $$
Na koΕcu obliczeΕ naleΕΌy zbinaryzowaΔ "pozycjΔ" by otrzymaΔ stan Isinga w rozwiΔ zaniu. UΕΌyjemy do tego funkcji znaku $\text{sign}(x)$. Dla przypomnienia:
$$ \text{sign}(x) = \left\{ \begin{array}{ll} 1, & \text{jeΕΌeli } x > 0 \\ 0, & \text{jeΕΌeli } x = 0\\ -1, & \text{jeΕΌeli } x < 0\\ \end{array} \right. $$
PseudokodΒΆ
from setup import setup_paths
setup_paths()
import numpy as np
from math import sqrt
from typing import Optional
from tqdm import tqdm
from funkcje_pomocnicze import calculate_energy
def wall(x: np.ndarray, y: np.ndarray):
mask = np.abs(x) > 1
x[mask] = np.sign(x[mask])
y[mask] = 0
return x, y
def balistic_simulated_bifurcation_single(J, h, num_steps, time_step,
a_0: Optional[float] = None, c_0_scaling: Optional[float] = None):
if a_0 is None:
a_0 = 1
N, _ = J.shape
mean_J = np.sqrt(np.sum(np.square(J)) / (N * (N - 1)) )
c_0 = 0.5 / (mean_J * sqrt(N))
if c_0_scaling is not None:
c_0 *= c_0_scaling
a = np.linspace(0, a_0, num=num_steps)
x = np.zeros(N)
y = np.random.uniform(-0.1, 0.1, N)
for t in tqdm(range(num_steps), desc="Symulowana Bifurkacja"):
y += (-1 * (a_0 - a[t]) * x + c_0 * (J @ x + h )) * time_step # y(t+1); x(t), x(t)
x += a_0 * y * time_step # x(t + 1); y(t+1)
x, y = wall(x, y)
x = np.sign(x)
return x, calculate_energy(J, h, x)
# test dla maΕej instancji
import numpy as np
import matplotlib.pyplot as plt
from funkcje_pomocnicze import read_instance, test_pegasus
J, h = read_instance(test_pegasus.path, convention="minus_half")
state, energy = balistic_simulated_bifurcation_single(J, h, num_steps=200, time_step=0.25, a_0=1, c_0_scaling=0.7)
print(f"Otrzymana energia: {energy}")
print(f"Stan podstawowy: {test_pegasus.best_energy}")
Symulowana Bifurkacja: 100%|ββββββββββ| 200/200 [00:00<00:00, 51536.57it/s]
Otrzymana energia: -467.0 Stan podstawowy: -469.0
ZrΓ³wnoleglenieΒΆ
Przez to ΕΌe algorytm sam w sobie jest stochastyczny nie mamy gwarancji ΕΌe w trakcie jednego wywoΕania otrzymamy najlepszy moΕΌliwy wynik. Jedym z rozwiΔ zaΔ tego problemu jest zrΓ³wnoleglenie algorytmu. Oznacza to, ΕΌe jednoczeΕnie Εledzimy wiele "trajektorii" i na koΕcu obliczeΕ moΕΌemy wybraΔ najlepszΔ z nich. MoΕΌna tego dokonaΔ wywoΕujΔ c algorytm wielokrotnie, ale jest to nieefektywne. Lepszym rozwiΔ zaniem jest trzymanie wielu stanΓ³w uΕoΕΌonych w macierz. Niech $\bm{x}$ bΔdzie wektorem zmiennych $x_i$, wtedy majΔ c $M$ trajektorii:
$$ \sigma = \begin{bmatrix} \vert & \vert & \cdots & \vert \\ \bm{x}_1 & \bm{x}_2 & \cdots & \bm{x}_M \\ \vert & \vert & \cdots & \vert \end{bmatrix} $$
WybΓ³r czy trzymamy wektory w kolumnach czy w wierszach jest arbitralny. To podejΕcie pozwala nam rΓ³wnolegle liczyΔ wiele rozwiΔ zaΕ i jest peΕnΔ formΔ tego algorytmu. Wtedy
W podobny sposΓ³b trzymamy $y$.. BΔdzie to wymagaΕo pewnych technicznych zmian w kodzie, ale gΕΓ³wna idea pozostaje niezmieniona.
# Implementacja Balistycznej Symulowanej Bifurkacji dla wielu trajektorii
from setup import setup_paths
setup_paths()
import numpy as np
from math import sqrt
from typing import Optional
from tqdm import tqdm
from funkcje_pomocnicze import calculate_energy_matrix
def wall(x: np.ndarray, y: np.ndarray):
mask = np.abs(x) > 1
x[mask] = np.sign(x[mask])
y[mask] = 0
return x, y
def balistic_simulated_bifurcation(J, h, num_steps, time_step, num_trajectories: int,
a_0: Optional[float] = None, c_0_scaling: Optional[float] = None):
if a_0 is None:
a_0 = 1
N, _ = J.shape
mean_J = np.sqrt(np.sum(np.square(J)) / (N * (N - 1)) )
c_0 = 0.5 / (mean_J * sqrt(N))
if c_0_scaling is not None:
c_0 *= c_0_scaling
a = np.linspace(0, a_0, num=num_steps)
x = np.zeros((N, num_trajectories))
y = np.random.uniform(-0.1, 0.1, (N, num_trajectories))
for t in tqdm(range(num_steps), desc="Symulowana Bifurkacja"):
y += (-1 * (a_0 - a[t]) * x + c_0 * (J @ x + h.reshape((N, 1)))) * time_step # y(t+1); x(t), x(t)
x += a_0 * y * time_step # x(t + 1); y(t+1)
x, y = wall(x, y)
x = np.sign(x)
return x, calculate_energy_matrix(J, h, x)
# Test na maΕym przykΕadzie z 216 spinami
# Balistyczna symulowana Bifurkacja
import numpy as np
import matplotlib.pyplot as plt
from funkcje_pomocnicze import read_instance, test_pegasus, Grid100
instance = test_pegasus
J, h = read_instance(instance.path, convention="minus_half")
states, energies = balistic_simulated_bifurcation(J, h, num_steps=200, time_step=0.25, num_trajectories=100, a_0=1, c_0_scaling=0.7)
print(f"Otrzymana energia: {min(energies)}")
print(f"Stan podstawowy: {instance.best_energy}")
print(f"Luka energetyczna: {((instance.best_energy - min(energies))/(2 * instance.best_energy) * 100):.2f}%")
unique_values, counts = np.unique(energies, return_counts=True)
plt.figure(figsize=(9, 6))
plt.bar(unique_values, counts, color='blue', edgecolor='black', alpha=0.7)
plt.xticks(energies)
plt.axvline(x=instance.best_energy, color='red', linestyle='--', linewidth=2)
plt.text(instance.best_energy, plt.ylim()[1]*0.8, 'Stan podstawowy', color='black', rotation=90,
verticalalignment='center', horizontalalignment='right')
plt.xlabel("Energia")
plt.ylabel("CzΔstoΕΔ wystΔpowania rozwiΔ
zania")
plt.title("RozkΕad znalezionych rozwiΔ
zaΕ dla instancji testowej o 216 spinach")
plt.grid(axis="y")
plt.show()
Symulowana Bifurkacja: 0%| | 0/200 [00:00<?, ?it/s]
Symulowana Bifurkacja: 100%|ββββββββββ| 200/200 [00:00<00:00, 3496.08it/s]
Otrzymana energia: -467.0 Stan podstawowy: -469.0 Luka energetyczna: 0.21%
Dyskretna Symulowana BifurkacjaΒΆ
Jest podobna do balistycznej, jedyna rΓ³ΕΌnica jest taka, ΕΌe podczas update parametru $y$ czΔΕΔ $J_{i,j}x_j(t_k)$ zastΔpujemy na $J_{i,j}\text{sign}(x_j(t_k))$, gdzie $\text{sign}$ jest funkcjΔ znaku.
Wtedy rΓ³wnania ktΓ³re opisujΔ ewolucjΔ parametrΓ³w $x$ i $y$ wyglΔ dajΔ nastΔpujΔ co:
$$ y_i(t_{k+1}) = y_i(t_k) + \{-[a_0 - a(t_k)]x_i(t_k) + c_0 \sum_{j=1}^N J_{i,j}\text{sign}(x_j(t_k)) \}\Delta_t $$ $$ x_i(t_{k+1}) = x_i(t_k) + a_0 y_i(t_{k+1})\Delta_t $$
Reszta algorytmu pozostaje bez zmian.
# Implementacja dyskretnej Symulowanej Bifurkacji
from setup import setup_paths
setup_paths()
import numpy as np
from math import sqrt
from typing import Optional
from tqdm import tqdm
from funkcje_pomocnicze import calculate_energy_matrix
def wall(x: np.ndarray, y: np.ndarray):
mask = np.abs(x) > 1
x[mask] = np.sign(x[mask])
y[mask] = 0
return x, y
def discrete_simulated_bifurcation(J, h, num_steps, time_step, num_trajectories: int,
a_0: Optional[float] = None, c_0_scaling: Optional[float] = None):
if a_0 is None:
a_0 = 1
N, _ = J.shape
mean_J = np.sqrt(np.sum(np.square(J)) / (N * (N - 1)) )
c_0 = 0.5 / (mean_J * sqrt(N))
if c_0_scaling is not None:
c_0 *= c_0_scaling
a = np.linspace(0, a_0, num=num_steps)
x = np.zeros((N, num_trajectories))
y = np.random.uniform(-0.1, 0.1, (N, num_trajectories))
for t in tqdm(range(num_steps), desc="Symulowana Bifurkacja"):
y += (-1 * (a_0 - a[t]) * x + c_0 * (J @ np.sign(x) + h.reshape((N, 1)))) * time_step # y(t+1); x(t), x(t)
x += a_0 * y * time_step # x(t + 1); y(t+1)
x, y = wall(x, y)
x = np.sign(x)
return x, calculate_energy_matrix(J, h, x)
# Test na maΕym przykΕadzie z 216 spinami
# Balistyczna symulowana Bifurkacja
import numpy as np
import matplotlib.pyplot as plt
from funkcje_pomocnicze import read_instance, test_pegasus
J, h = read_instance(test_pegasus.path, convention="minus_half")
states, energies = discrete_simulated_bifurcation(J, h, num_steps=200, time_step=0.5, num_trajectories=100, a_0=1, c_0_scaling=0.7)
print(f"Otrzymana energia: {min(energies)}")
print(f"Stan podstawowy: {test_pegasus.best_energy}")
unique_values, counts = np.unique(energies, return_counts=True)
plt.figure(figsize=(9, 6))
plt.bar(unique_values, counts, color='blue', edgecolor='black', alpha=0.7)
#plt.xticks(energies)
plt.axvline(x=test_pegasus.best_energy, color='red', linestyle='--', linewidth=2)
plt.text(test_pegasus.best_energy, plt.ylim()[1]*0.8, 'Stan podstawowy', color='black', rotation=90,
verticalalignment='center', horizontalalignment='right')
plt.xlabel("Energia")
plt.ylabel("CzΔstoΕΔ wystΔpowania rozwiΔ
zania")
plt.title("RozkΕad znalezionych rozwiΔ
zaΕ dla instancji testowej o 216 spinach")
plt.grid(axis="y")
plt.show()
Symulowana Bifurkacja: 0%| | 0/200 [00:00<?, ?it/s]
Symulowana Bifurkacja: 100%|ββββββββββ| 200/200 [00:00<00:00, 1999.01it/s]
Otrzymana energia: -469.0 Stan podstawowy: -469.0
# PorΓ³wnanie jakoΕci rozwiΔ
zaΕ i czasu obliczeΕ dla obydwu wersji
import time
import matplotlib.pylab as plt
from funkcje_pomocnicze import read_instance, P2, P4, P6, P8, P12
from IPython.utils.io import capture_output
num_traj = 2**7
num_steps = 200
time_step = 0.75
instances = [P2, P4, P6, P8, P12]
gaps = {}
times = {}
for instance in tqdm(instances, desc="Zbieranie danych"):
J, h = read_instance(instance.path, convention="minus_half")
with capture_output() as captured:
begin = time.time()
_, energies1 = balistic_simulated_bifurcation(J, h, num_steps=num_steps,
time_step=time_step, num_trajectories=num_traj)
end = time.time()
elapsed_balistic = end - begin
begin = time.time()
_, energies2 = discrete_simulated_bifurcation(J, h, num_steps=num_steps,
time_step=time_step, num_trajectories=num_traj)
end = time.time()
elapsed_discrete = end - begin
gap_balistic = (instance.best_energy - min(energies1)) / (2*instance.best_energy) * 100
gap_discrete = (instance.best_energy - min(energies2)) /(2*instance.best_energy) * 100
gaps[(instance.name, "balistyczna")] = gap_balistic
gaps[(instance.name, "dyskretna")] = gap_discrete
times[(instance.name, "balistyczna")] = elapsed_balistic
times[(instance.name, "dyskretna")] = elapsed_discrete
instancje = [instance.name for instance in instances]
x = np.arange(len(instancje))
width = 0.25
offset = 0.15
fig, ax = plt.subplots(figsize=(12, 6))
barset = []
param = [-1, 1]
for idx, s in enumerate(["balistyczna", "dyskretna"]):
a = [gaps[(i, s)] for i in instancje]
bar = ax.bar(x + offset * param[idx], a, width, label=rf"{s} symulowana bifurkacja")
barset.append(bar)
(y_min, y_max) = ax.get_ylim()
if y_max <1.1:
y_max = 1.1
ax.set_ylim((-0.2, y_max))
ax.set_xlabel("Rozmiar instancji")
ax.set_xticks(x)
ax.set_xticklabels(instancje, rotation=45)
ax.set_ylabel("Luka energetyczna [%]")
ax.set_title("Luka energetyczna dla rΓ³ΕΌnych rozmiarΓ³w instancji uΕΌywajΔ
c tych samych wartoΕci parametrΓ³w")
ax.axhline(0, color='black', linestyle='--')
ax.legend()
for bars in barset:
for bar in bars:
if bar.get_height() == 0:
center = bar.get_x() + bar.get_width()/2
ax.plot(center, 0, marker='*', markersize=10, color=bar.get_facecolor())
plt.grid(axis="y")
plt.show()
fig2, ax2 = plt.subplots(figsize=(12, 6))
for idx, s in enumerate(["balistyczna", "dyskretna"]):
a = [times[(i, s)] for i in instancje]
bar = ax2.bar(x + offset * param[idx], a, width, label=rf"{s} symulowana bifurkacja")
ax2.set_xlabel("Rozmiar instancji")
ax2.set_xticks(x)
ax2.set_xticklabels(instancje, rotation=45)
ax2.set_ylabel("Czas [s]")
ax2.set_title("Czas do znalezienia rozwiΔ
zania dla rΓ³ΕΌnych rozmiarΓ³w instancji uΕΌywajΔ
c tych samych wartoΕci parametrΓ³w")
ax2.legend()
plt.grid(axis="y")
plt.show()
Zbieranie danych: 100%|ββββββββββ| 5/5 [00:09<00:00, 1.90s/it]
WpΕyw parametrΓ³w na zachowanie algorytmuΒΆ
Skupimy siΔ na dyskretnej symulowanej bifurkacji.
# krok czasowy delta t
import time
import matplotlib.pylab as plt
from funkcje_pomocnicze import read_instance, P2, P4, P6, P8
from IPython.utils.io import capture_output
num_traj = 2**7
num_steps = 200
time_steps = [0.1, 0.25, 0.5, 0.75, 1.0]
instances = [P2, P4, P6, P8]
gaps = {}
times = {}
with tqdm(total=len(instances)*len(time_steps), desc="zbieranie danych") as pbar:
for t in time_steps:
for instance in instances:
J, h = read_instance(instance.path, convention="minus_half")
with capture_output() as captured:
begin = time.time()
_, energies = discrete_simulated_bifurcation(J, h, num_steps=num_steps,
time_step=t, num_trajectories=num_traj)
end = time.time()
elapsed_discrete = end - begin
gap = (instance.best_energy - min(energies)) / (2*instance.best_energy) * 100
if gap < 0:
print(instance.name, min(energies))
gaps[(instance.name, t)] = gap
times[(instance.name, t)] = elapsed_discrete
pbar.update(1)
instancje = [instance.name for instance in instances]
x = np.arange(len(instancje))
width = 0.1
offset = 0.15
fig, ax = plt.subplots(figsize=(12, 6))
barset = []
param = [-2, -1, 0, 1, 2]
for idx, s in enumerate(time_steps):
a = [gaps[(i, s)] for i in instancje]
bar = ax.bar(x + offset * param[idx], a, width, label=rf"$\Delta t =${s}")
barset.append(bar)
(y_min, y_max) = ax.get_ylim()
if y_max <1.1:
y_max = 1.1
ax.set_ylim((-0.2, y_max))
ax.set_xlabel("Rozmiar instancji")
ax.set_xticks(x)
ax.set_xticklabels(instancje, rotation=45)
ax.set_ylabel("Luka energetyczna [%]")
ax.set_title("Luka energetyczna dla rΓ³ΕΌnych rozmiarΓ³w instancji i rΓ³ΕΌnych wartoΕci kroku czasowego")
ax.axhline(0, color='black', linestyle='--')
ax.legend()
for bars in barset:
for bar in bars:
if bar.get_height() == 0:
center = bar.get_x() + bar.get_width()/2
ax.plot(center, 0, marker='*', markersize=10, color=bar.get_facecolor())
plt.grid(axis="y")
plt.show()
zbieranie danych: 100%|ββββββββββ| 20/20 [00:10<00:00, 1.92it/s]
# IloΕΔ trajektorii
import matplotlib.pylab as plt
from funkcje_pomocnicze import read_instance, P2, P4, P6, P8
from IPython.utils.io import capture_output
from tqdm import tqdm
time_step = 0.75
num_steps = 200
instances = [P2, P4, P6, P8]
trajectories = [2**6, 2**7, 2**8, 2**9, 2**10]
gaps = {}
times = {}
with tqdm(total=len(instances)*len(trajectories), desc="zbieranie danych") as pbar:
for t in trajectories:
for instance in instances:
J, h = read_instance(instance.path, convention="minus_half")
with capture_output() as captured:
begin = time.time()
_, energies = discrete_simulated_bifurcation(J, h, num_steps=num_steps,
time_step=time_step, num_trajectories=t)
end = time.time()
elapsed_discrete = end - begin
gap = (instance.best_energy - min(energies)) / (2*instance.best_energy) * 100
if gap < 0:
print(instance.name, min(energies))
gaps[(instance.name, t)] = gap
times[(instance.name, t)] = elapsed_discrete
pbar.update(1)
instancje = [instance.name for instance in instances]
x = np.arange(len(instancje))
width = 0.1
offset = 0.15
fig, ax = plt.subplots(figsize=(12, 6))
barset = []
param = [-2, -1, 0, 1, 2]
for idx, s in enumerate(trajectories):
a = [gaps[(i, s)] for i in instancje]
bar = ax.bar(x + offset * param[idx], a, width, label=rf"{s} trajektorii")
barset.append(bar)
(y_min, y_max) = ax.get_ylim()
if y_max <1.1:
y_max = 1.1
ax.set_ylim((-0.2, y_max))
ax.set_xlabel("Rozmiar instancji")
ax.set_xticks(x)
ax.set_xticklabels(instancje, rotation=45)
ax.set_ylabel("Luka energetyczna [%]")
ax.set_title("Luka energetyczna dla rΓ³ΕΌnych rozmiarΓ³w instancji i rΓ³ΕΌnej iloΕci Εledzonych trajektorii")
ax.axhline(0, color='black', linestyle='--')
ax.legend()
for bars in barset:
for bar in bars:
if bar.get_height() == 0:
center = bar.get_x() + bar.get_width()/2
ax.plot(center, 0, marker='*', markersize=10, color=bar.get_facecolor())
plt.grid(axis="y")
plt.show()
fig2, ax2 = plt.subplots(figsize=(12, 6))
for idx, s in enumerate(trajectories):
a = [times[(i, s)] for i in instancje]
bar = ax2.bar(x + offset * param[idx], a, width, label=rf"{s} trajektorii")
ax2.set_xlabel("Rozmiar instancji")
ax2.set_xticks(x)
ax2.set_xticklabels(instancje, rotation=45)
ax2.set_ylabel("Czas [s]")
ax2.set_title("Czas do znalezienia rozwiΔ
zania dla rΓ³ΕΌnych rozmiarΓ³w instancji i rΓ³ΕΌnej iloΕci Εledzonych trajektorii")
ax2.legend()
plt.grid(axis="y")
plt.show()
zbieranie danych: 100%|ββββββββββ| 20/20 [00:29<00:00, 1.48s/it]
# IloΕΔ krokΓ³w
import matplotlib.pylab as plt
from funkcje_pomocnicze import read_instance, P2, P4, P6, P8
from IPython.utils.io import capture_output
from tqdm import tqdm
time_step = 0.75
num_trajectories = 2**7
instances = [P2, P4, P6, P8]
steps = [50, 100, 200, 500, 1000]
gaps = {}
times = {}
with tqdm(total=len(instances)*len(steps), desc="zbieranie danych") as pbar:
for t in steps:
for instance in instances:
J, h = read_instance(instance.path, convention="minus_half")
with capture_output() as captured:
begin = time.time()
_, energies = discrete_simulated_bifurcation(J, h, num_steps=t,
time_step=time_step, num_trajectories=num_trajectories)
end = time.time()
elapsed_discrete = end - begin
gap = (instance.best_energy - min(energies)) / (2*instance.best_energy) * 100
if gap < 0:
print(instance.name, min(energies))
gaps[(instance.name, t)] = gap
times[(instance.name, t)] = elapsed_discrete
pbar.update(1)
instancje = [instance.name for instance in instances]
x = np.arange(len(instancje))
width = 0.1
offset = 0.15
fig, ax = plt.subplots(figsize=(12, 6))
barset = []
param = [-2, -1, 0, 1, 2]
for idx, s in enumerate(steps):
a = [gaps[(i, s)] for i in instancje]
bar = ax.bar(x + offset * param[idx], a, width, label=rf"{s} krokΓ³w")
barset.append(bar)
(y_min, y_max) = ax.get_ylim()
if y_max <1.1:
y_max = 1.1
ax.set_ylim((-0.2, y_max))
ax.set_xlabel("Rozmiar instancji")
ax.set_xticks(x)
ax.set_xticklabels(instancje, rotation=45)
ax.set_ylabel("Luka energetyczna [%]")
ax.set_title("Luka energetyczna dla rΓ³ΕΌnych rozmiarΓ³w instancji i rΓ³ΕΌnej iloΕci krokΓ³w")
ax.axhline(0, color='black', linestyle='--')
ax.legend()
for bars in barset:
for bar in bars:
if bar.get_height() == 0:
center = bar.get_x() + bar.get_width()/2
ax.plot(center, 0, marker='*', markersize=10, color=bar.get_facecolor())
plt.grid(axis="y")
plt.show()
fig2, ax2 = plt.subplots(figsize=(12, 6))
for idx, s in enumerate(steps):
a = [times[(i, s)] for i in instancje]
bar = ax2.bar(x + offset * param[idx], a, width, label=rf"{s} krokΓ³w")
ax2.set_xlabel("Rozmiar instancji")
ax2.set_xticks(x)
ax2.set_xticklabels(instancje, rotation=45)
ax2.set_ylabel("Czas [s]")
ax2.set_title("Czas do znalezienia rozwiΔ
zania dla rΓ³ΕΌnych rozmiarΓ³w instancji i rΓ³ΕΌnej iloΕci krokΓ³w")
ax2.legend()
plt.grid(axis="y")
plt.show()
zbieranie danych: 100%|ββββββββββ| 20/20 [00:17<00:00, 1.17it/s]
BibliografiaΒΆ
Goto, H., et al., Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems. Sci. Adv. 5, eaav2372 (2019). DOI:10.1126/sciadv.aav2372
Goto, H., et al., High-performance combinatorial optimization based on classical mechanics. Sci. Adv. 7, eabe7953 (2021).
DOI:10.1126/sciadv.abe7953Pawlowski, J., et al. "Closing the Quantum-Classical Scaling Gap in Approximate Optimization." arXiv preprint arXiv:2505.22514 (2025). DOI:10.48550/arXiv.2505.22514