From f5630543ac7bc78994ec9319cae73c95cda4a5d0 Mon Sep 17 00:00:00 2001 From: Thomas Bouquet Date: Sat, 30 Dec 2023 18:47:43 +0100 Subject: [PATCH 1/5] Operational presolver --- .../optimizer/__init__.py | 1 + .../optimizer/preprocessing.py | 70 ++++++++ tests/optimizer/test_preprocessing.py | 154 ++++++++++++++++++ 3 files changed, 225 insertions(+) create mode 100644 src/simulated_bifurcation/optimizer/preprocessing.py create mode 100644 tests/optimizer/test_preprocessing.py diff --git a/src/simulated_bifurcation/optimizer/__init__.py b/src/simulated_bifurcation/optimizer/__init__.py index bcb5a3ce..d461dbc3 100644 --- a/src/simulated_bifurcation/optimizer/__init__.py +++ b/src/simulated_bifurcation/optimizer/__init__.py @@ -14,6 +14,7 @@ """ from .environment import get_env, reset_env, set_env +from .preprocessing import Preprocessing from .simulated_bifurcation_engine import SimulatedBifurcationEngine from .simulated_bifurcation_optimizer import ( ConvergenceWarning, diff --git a/src/simulated_bifurcation/optimizer/preprocessing.py b/src/simulated_bifurcation/optimizer/preprocessing.py new file mode 100644 index 00000000..0a917278 --- /dev/null +++ b/src/simulated_bifurcation/optimizer/preprocessing.py @@ -0,0 +1,70 @@ +from typing import Optional, Tuple + +import torch + + +class Preprocessing: + def __init__(self, J: torch.Tensor, h: torch.Tensor) -> None: + self.J = J.clone() + self.h = h.clone() + self.optimized_spins = torch.zeros(self.J.shape[0]) + self.shifted_indices = list(range(self.J.shape[0])) + self.dimension = self.J.shape[0] + + def _remove_row(self, index: int): + self.J = torch.cat((self.J[:index], self.J[index + 1 :]), dim=0) + self.h = torch.cat((self.h[:index], self.h[index + 1 :])) + + def _remove_column(self, index: int): + self.J = torch.cat((self.J[:, :index], self.J[:, index + 1 :]), dim=1) + + def _remove_all_coefficients(self, index: int): + self._remove_row(index) + self._remove_column(index) + + def _project_coefficients_in_linear_part(self, row_index: int, sign: int): + self.h += sign * self.J[row_index] + + def _project_coefficients_and_delete_row_and_column( + self, spin_index: int, spin_sign: int + ): + self._project_coefficients_in_linear_part(spin_index, spin_sign) + self._remove_all_coefficients(spin_index) + + def _get_optimizable_spins(self) -> torch.Tensor: + return torch.abs(self.h) >= torch.sum(torch.abs(self.J), dim=0) + + def _get_first_optimizable_spin(self) -> Optional[int]: + optimizable_spins = self._get_optimizable_spins() + if torch.any(optimizable_spins).item(): + return torch.nonzero(optimizable_spins)[0].item() + return None + + def _delete_index(self, index: int): + self.shifted_indices = ( + self.shifted_indices[:index] + self.shifted_indices[index + 1 :] + ) + + def _get_original_index(self, new_index: int) -> int: + return self.shifted_indices[new_index] + + def _set_spin_value(self, spin_index: int, spin_sign: int): + self.optimized_spins[spin_index] = spin_sign + + def _set_first_optimal_spin(self) -> bool: + optimal_spin_index = self._get_first_optimizable_spin() + if optimal_spin_index is None: + return False + sign = -1 if self.h[optimal_spin_index] > 0 else 1 + original_index = self._get_original_index(optimal_spin_index) + self._set_spin_value(original_index, sign) + self._project_coefficients_and_delete_row_and_column(optimal_spin_index, sign) + self._delete_index(optimal_spin_index) + return True + + def presolve(self) -> Tuple[torch.Tensor, torch.Tensor, torch.Tensor]: + _continue = True + while _continue: + optimized_spin = self._set_first_optimal_spin() + _continue = _continue and optimized_spin and self.J.shape[1] > 0 + return self.optimized_spins.clone(), self.J.clone(), self.h.clone() diff --git a/tests/optimizer/test_preprocessing.py b/tests/optimizer/test_preprocessing.py new file mode 100644 index 00000000..bff2983b --- /dev/null +++ b/tests/optimizer/test_preprocessing.py @@ -0,0 +1,154 @@ +import pytest +import torch + +from src.simulated_bifurcation.optimizer import Preprocessing + +J = torch.tensor( + [ + [1, 2, 3], + [2, 1, 4], + [3, 4, 1], + ], + dtype=torch.float32, +) +h = torch.tensor([1, 0, -2], dtype=torch.float32) + +optimizable_J = torch.tensor( + [ + [0.0, 0.5, -1.0], + [0.5, 0.0, 2.0], + [-1.0, 2.0, 0.0], + ], + dtype=torch.float32, +) +optimizable_h = torch.tensor([2.0, 1.0, -4.0], dtype=torch.float32) + + +def test_remove_row(): + preprocesser = Preprocessing(J, h) + preprocesser._remove_row(1) + assert torch.equal( + preprocesser.J, torch.tensor([[1, 2, 3], [3, 4, 1]], dtype=torch.float32) + ) + assert torch.equal(preprocesser.h, torch.tensor([1, -2], dtype=torch.float32)) + preprocesser._remove_row(0) + assert torch.equal(preprocesser.J, torch.tensor([[3, 4, 1]], dtype=torch.float32)) + assert torch.equal(preprocesser.h, torch.tensor([-2], dtype=torch.float32)) + + +def test_remove_column(): + preprocesser = Preprocessing(J, h) + preprocesser._remove_column(1) + assert torch.equal( + preprocesser.J, torch.tensor([[1, 3], [2, 4], [3, 1]], dtype=torch.float32) + ) + assert torch.equal(preprocesser.h, h) + preprocesser._remove_column(0) + assert torch.equal( + preprocesser.J, torch.tensor([[3], [4], [1]], dtype=torch.float32) + ) + assert torch.equal(preprocesser.h, h) + + +def test_remove_all_coefficients(): + preprocesser = Preprocessing(J, h) + preprocesser._remove_all_coefficients(1) + assert torch.equal( + preprocesser.J, torch.tensor([[1, 3], [3, 1]], dtype=torch.float32) + ) + assert torch.equal(preprocesser.h, torch.tensor([1, -2], dtype=torch.float32)) + preprocesser._remove_all_coefficients(0) + assert torch.equal(preprocesser.J, torch.tensor([[1]], dtype=torch.float32)) + assert torch.equal(preprocesser.h, torch.tensor([-2], dtype=torch.float32)) + + +def test_project_coefficients_in_linear_part(): + preprocesser = Preprocessing(J, h) + preprocesser._project_coefficients_in_linear_part(1, 1) + assert torch.equal(preprocesser.J, J) + assert torch.equal(preprocesser.h, torch.tensor([3, 1, 2], dtype=torch.float32)) + preprocesser._project_coefficients_in_linear_part(0, -1) + assert torch.equal(preprocesser.J, J) + assert torch.equal(preprocesser.h, torch.tensor([2, -1, -1], dtype=torch.float32)) + + +def test_project_coefficients_and_delete_row_and_column(): + preprocesser = Preprocessing(J, h) + preprocesser._project_coefficients_and_delete_row_and_column(2, -1) + assert torch.equal( + preprocesser.J, torch.tensor([[1, 2], [2, 1]], dtype=torch.float32) + ) + assert torch.equal(preprocesser.h, torch.tensor([-2, -4], dtype=torch.float32)) + + +def test_get_optimizable_spins(): + preprocesser = Preprocessing(optimizable_J, optimizable_h) + assert torch.equal( + preprocesser._get_optimizable_spins(), torch.Tensor([True, False, True]) + ) + + +def test_get_first_optimizable_spin(): + optimizable_preprocesser = Preprocessing(optimizable_J, optimizable_h) + assert 0 == optimizable_preprocesser._get_first_optimizable_spin() + non_optimizable_preprocesser = Preprocessing(J, h) + assert non_optimizable_preprocesser._get_first_optimizable_spin() is None + + +def test_delete_index(): + preprocesser = Preprocessing(J, h) + assert [0, 1, 2] == preprocesser.shifted_indices + preprocesser._delete_index(1) + assert [0, 2] == preprocesser.shifted_indices + preprocesser._delete_index(0) + assert [2] == preprocesser.shifted_indices + + +def test_get_original_index(): + preprocesser = Preprocessing(J, h) + preprocesser._delete_index(1) + assert 0 == preprocesser._get_original_index(0) + assert 2 == preprocesser._get_original_index(1) + + +def test_set_spin_value(): + preprocesser = Preprocessing(J, h) + assert torch.equal(preprocesser.optimized_spins, torch.tensor([0, 0, 0])) + preprocesser._set_spin_value(2, -1) + assert torch.equal(preprocesser.optimized_spins, torch.tensor([0, 0, -1])) + preprocesser._set_spin_value(0, 1) + assert torch.equal(preprocesser.optimized_spins, torch.tensor([1, 0, -1])) + + +def test_set_first_optimal_spin(): + non_optimizable_preprocesser = Preprocessing(J, h) + assert non_optimizable_preprocesser._set_first_optimal_spin() is False + optimizable_preprocesser = Preprocessing(optimizable_J, optimizable_h) + assert optimizable_preprocesser._set_first_optimal_spin() is True + assert torch.equal( + optimizable_preprocesser.J, torch.tensor([[0, 2], [2, 0]], dtype=torch.float32) + ) + assert torch.equal( + optimizable_preprocesser.h, torch.tensor([0.5, -3], dtype=torch.float32) + ) + assert torch.equal( + optimizable_preprocesser.optimized_spins, torch.tensor([-1, 0, 0]) + ) + assert optimizable_preprocesser.shifted_indices == [1, 2] + + +def test_presolve(): + non_optimizable_preprocesser = Preprocessing(J, h) + ( + non_optimized_spins, + non_optimized_J, + non_optimized_h, + ) = non_optimizable_preprocesser.presolve() + assert torch.equal(non_optimized_spins, torch.tensor([0, 0, 0])) + assert torch.equal(non_optimized_J, J) + assert torch.equal(non_optimized_h, h) + optimizable_preprocesser = Preprocessing(optimizable_J, optimizable_h) + optimized_spins, optimized_J, optimized_h = optimizable_preprocesser.presolve() + assert torch.equal(optimized_spins, torch.tensor([-1, -1, 1])) + assert torch.equal(optimized_J, torch.tensor([]).reshape(0, 0)) + assert torch.equal(optimized_h, torch.tensor([])) From 60a1680a22dd76bd04935054640df203f1600526 Mon Sep 17 00:00:00 2001 From: Thomas Bouquet Date: Sun, 31 Dec 2023 01:23:58 +0100 Subject: [PATCH 2/5] presolve parameter for optimization --- src/simulated_bifurcation/core/ising.py | 32 ++++++++++++++++++- .../core/quadratic_polynomial.py | 18 +++++++++++ .../optimizer/__init__.py | 1 + .../optimizer/postprocessing.py | 21 ++++++++++++ .../optimizer/preprocessing.py | 1 - .../simulated_bifurcation.py | 18 +++++++++++ tests/core/test_ising.py | 23 +++++++++++++ tests/optimizer/test_postprocessing.py | 20 ++++++++++++ tests/optimizer/test_preprocessing.py | 2 +- 9 files changed, 133 insertions(+), 3 deletions(-) create mode 100644 src/simulated_bifurcation/optimizer/postprocessing.py create mode 100644 tests/optimizer/test_postprocessing.py diff --git a/src/simulated_bifurcation/core/ising.py b/src/simulated_bifurcation/core/ising.py index 543a7e84..760bfa23 100644 --- a/src/simulated_bifurcation/core/ising.py +++ b/src/simulated_bifurcation/core/ising.py @@ -24,7 +24,12 @@ import torch from numpy import ndarray -from ..optimizer import SimulatedBifurcationEngine, SimulatedBifurcationOptimizer +from ..optimizer import ( + Postprocessing, + Preprocessing, + SimulatedBifurcationEngine, + SimulatedBifurcationOptimizer, +) # Workaround because `Self` type is only available in Python >= 3.11 SelfIsing = TypeVar("SelfIsing", bound="Ising") @@ -251,6 +256,7 @@ def minimize( sampling_period: int = 50, convergence_threshold: int = 50, timeout: Optional[float] = None, + presolve: bool = False ) -> None: """ Minimize the energy of the Ising model using the Simulated Bifurcation @@ -402,6 +408,30 @@ def minimize( sampling_period, convergence_threshold, ) + if presolve: + presolved_spins, reduced_J, reduced_h = Preprocessing( + self.J, self.h + ).presolve() + if reduced_J.shape[1] == 0: + self.computed_spins = presolved_spins.repeat(agents, 1).t() + return + reduced_model = Ising(reduced_J, reduced_h, self.dtype, self.device) + reduced_model.minimize( + agents=agents, + max_steps=max_steps, + ballistic=ballistic, + heated=heated, + verbose=verbose, + use_window=use_window, + sampling_period=sampling_period, + convergence_threshold=convergence_threshold, + timeout=timeout, + presolve=False, + ) + self.computed_spins = Postprocessing.reconstruct_spins( + reduced_model.computed_spins, presolved_spins + ) + return tensor = self.as_simulated_bifurcation_tensor() spins = optimizer.run_integrator(tensor, use_window) if self.linear_term: diff --git a/src/simulated_bifurcation/core/quadratic_polynomial.py b/src/simulated_bifurcation/core/quadratic_polynomial.py index cbe34865..09fbd0b5 100644 --- a/src/simulated_bifurcation/core/quadratic_polynomial.py +++ b/src/simulated_bifurcation/core/quadratic_polynomial.py @@ -322,6 +322,7 @@ def optimize( sampling_period: int = 50, convergence_threshold: int = 50, timeout: Optional[float] = None, + presolve: bool = False, ) -> Tuple[torch.Tensor, torch.Tensor]: """ Computes a local extremum of the model by optimizing @@ -395,6 +396,10 @@ def optimize( minimize : bool, optional if `True` the optimization direction is minimization, otherwise it is maximization (default is True) + presolve : bool, optional + if `True`, the model will be preprocessed to find spins that can be + optimized before entering the SB algorithm, reducing the + computation size (default is False) Returns ------- @@ -414,6 +419,7 @@ def optimize( sampling_period=sampling_period, convergence_threshold=convergence_threshold, timeout=timeout, + presolve=presolve, ) self.sb_result = self.convert_spins(ising_equivalent, domain) result = self.sb_result.t().to(dtype=self.dtype) @@ -438,6 +444,7 @@ def minimize( sampling_period: int = 50, convergence_threshold: int = 50, timeout: Optional[float] = None, + presolve: bool = False, ) -> Tuple[torch.Tensor, torch.Tensor]: """ Computes a local minimum of the model by optimizing @@ -508,6 +515,10 @@ def minimize( if `True` only the best found solution to the optimization problem is returned, otherwise all the solutions found by the simulated bifurcation algorithm. + presolve : bool, optional + if `True`, the model will be preprocessed to find spins that can be + optimized before entering the SB algorithm, reducing the + computation size (default is False) Returns ------- @@ -526,6 +537,7 @@ def minimize( sampling_period=sampling_period, convergence_threshold=convergence_threshold, timeout=timeout, + presolve=presolve, ) def maximize( @@ -542,6 +554,7 @@ def maximize( sampling_period: int = 50, convergence_threshold: int = 50, timeout: Optional[float] = None, + presolve: bool = False, ) -> Tuple[torch.Tensor, torch.Tensor]: """ Computes a local maximum of the model by optimizing @@ -611,6 +624,10 @@ def maximize( if `True` only the best found solution to the optimization problem is returned, otherwise all the solutions found by the simulated bifurcation algorithm. + presolve : bool, optional + if `True`, the model will be preprocessed to find spins that can be + optimized before entering the SB algorithm, reducing the + computation size (default is False) Returns ------- @@ -629,6 +646,7 @@ def maximize( sampling_period=sampling_period, convergence_threshold=convergence_threshold, timeout=timeout, + presolve=presolve, ) @staticmethod diff --git a/src/simulated_bifurcation/optimizer/__init__.py b/src/simulated_bifurcation/optimizer/__init__.py index d461dbc3..fd29bcfb 100644 --- a/src/simulated_bifurcation/optimizer/__init__.py +++ b/src/simulated_bifurcation/optimizer/__init__.py @@ -14,6 +14,7 @@ """ from .environment import get_env, reset_env, set_env +from .postprocessing import Postprocessing from .preprocessing import Preprocessing from .simulated_bifurcation_engine import SimulatedBifurcationEngine from .simulated_bifurcation_optimizer import ( diff --git a/src/simulated_bifurcation/optimizer/postprocessing.py b/src/simulated_bifurcation/optimizer/postprocessing.py new file mode 100644 index 00000000..c157ae5b --- /dev/null +++ b/src/simulated_bifurcation/optimizer/postprocessing.py @@ -0,0 +1,21 @@ +import torch + + +class Postprocessing: + @staticmethod + def reconstruct_spins( + optimized_spins: torch.Tensor, pre_solved_spins: torch.Tensor + ) -> torch.Tensor: + original_dimension = pre_solved_spins.shape[0] + agents = optimized_spins.shape[1] + reconstructed_spins = torch.zeros(original_dimension, agents, dtype=torch.int32) + reconstructed_spins[pre_solved_spins == 0] = optimized_spins.clone().to( + dtype=torch.int32 + ) + reconstructed_spins[torch.abs(pre_solved_spins) == 1] = ( + pre_solved_spins[torch.abs(pre_solved_spins) == 1] + .repeat(agents, 1) + .t() + .to(dtype=torch.int32) + ) + return reconstructed_spins diff --git a/src/simulated_bifurcation/optimizer/preprocessing.py b/src/simulated_bifurcation/optimizer/preprocessing.py index 0a917278..c6ca4798 100644 --- a/src/simulated_bifurcation/optimizer/preprocessing.py +++ b/src/simulated_bifurcation/optimizer/preprocessing.py @@ -9,7 +9,6 @@ def __init__(self, J: torch.Tensor, h: torch.Tensor) -> None: self.h = h.clone() self.optimized_spins = torch.zeros(self.J.shape[0]) self.shifted_indices = list(range(self.J.shape[0])) - self.dimension = self.J.shape[0] def _remove_row(self, index: int): self.J = torch.cat((self.J[:index], self.J[index + 1 :]), dim=0) diff --git a/src/simulated_bifurcation/simulated_bifurcation.py b/src/simulated_bifurcation/simulated_bifurcation.py index 0e55e8e7..d207a2b4 100644 --- a/src/simulated_bifurcation/simulated_bifurcation.py +++ b/src/simulated_bifurcation/simulated_bifurcation.py @@ -172,6 +172,7 @@ def optimize( sampling_period: int = 50, convergence_threshold: int = 50, timeout: Optional[float] = None, + presolve: bool = False, ) -> Tuple[torch.Tensor, torch.Tensor]: """ Optimize a multivariate quadratic polynomial using the @@ -256,6 +257,10 @@ def optimize( timeout : float | None, default=None, keyword-only Time, in seconds, after which the simulation will be stopped. None means no timeout. + presolve : bool, optional + if `True`, the model will be preprocessed to find spins that can be + optimized before entering the SB algorithm, reducing the + computation size (default is False) Returns ------- @@ -430,6 +435,7 @@ def optimize( sampling_period=sampling_period, convergence_threshold=convergence_threshold, timeout=timeout, + presolve=presolve, ) return result, evaluation @@ -449,6 +455,7 @@ def minimize( sampling_period: int = 50, convergence_threshold: int = 50, timeout: Optional[float] = None, + presolve: bool = False, ) -> Tuple[torch.Tensor, torch.Tensor]: """ Minimize a multivariate quadratic polynomial using the @@ -530,6 +537,10 @@ def minimize( timeout : float | None, default=None, keyword-only Time, in seconds, after which the simulation will be stopped. None means no timeout. + presolve : bool, optional + if `True`, the model will be preprocessed to find spins that can be + optimized before entering the SB algorithm, reducing the + computation size (default is False) Returns ------- @@ -694,6 +705,7 @@ def minimize( sampling_period=sampling_period, convergence_threshold=convergence_threshold, timeout=timeout, + presolve=presolve, ) @@ -712,6 +724,7 @@ def maximize( sampling_period: int = 50, convergence_threshold: int = 50, timeout: Optional[float] = None, + presolve: bool = False, ) -> Tuple[torch.Tensor, torch.Tensor]: """ Maximize a multivariate quadratic polynomial using the @@ -793,6 +806,10 @@ def maximize( timeout : float | None, default=None, keyword-only Time, in seconds, after which the simulation will be stopped. None means no timeout. + presolve : bool, optional + if `True`, the model will be preprocessed to find spins that can be + optimized before entering the SB algorithm, reducing the + computation size (default is False) Returns ------- @@ -957,4 +974,5 @@ def maximize( sampling_period=sampling_period, convergence_threshold=convergence_threshold, timeout=timeout, + presolve=presolve, ) diff --git a/tests/core/test_ising.py b/tests/core/test_ising.py index e88944a3..6ddf1f06 100644 --- a/tests/core/test_ising.py +++ b/tests/core/test_ising.py @@ -95,3 +95,26 @@ def test_negative_ising(): assert negative_ising.linear_term assert len(negative_ising) == 3 assert negative_ising.dimension == 3 + + +def test_presolve_ising(): + ising = Ising(J, h) + ising.minimize(agents=10, presolve=True) + assert tuple(ising.computed_spins.shape) == (3, 10) + + +def test_presolve_solves_all_spins(): + presolvable_J = torch.tensor( + [ + [0.0, 0.5, -1.0], + [0.5, 0.0, 2.0], + [-1.0, 2.0, 0.0], + ], + dtype=torch.float32, + ) + presolvable_h = torch.tensor([2.0, 1.0, -4.0], dtype=torch.float32) + ising = Ising(presolvable_J, presolvable_h) + ising.minimize(agents=3, presolve=True) + assert torch.equal( + ising.computed_spins, torch.tensor([[-1, -1, -1], [-1, -1, -1], [1, 1, 1]]) + ) diff --git a/tests/optimizer/test_postprocessing.py b/tests/optimizer/test_postprocessing.py new file mode 100644 index 00000000..c5cb32eb --- /dev/null +++ b/tests/optimizer/test_postprocessing.py @@ -0,0 +1,20 @@ +import torch + +from src.simulated_bifurcation.optimizer.postprocessing import Postprocessing + + +def test_reconstruct_spins(): + optimized_spins = torch.tensor([[1, 1, -1, 1], [-1, 1, -1, -1], [1, 1, -1, 1]]) + pre_solved_spins = torch.tensor([1, 0, 0, -1, 0]) + assert torch.equal( + torch.tensor( + [ + [1, 1, 1, 1], + [1, 1, -1, 1], + [-1, 1, -1, -1], + [-1, -1, -1, -1], + [1, 1, -1, 1], + ] + ), + Postprocessing.reconstruct_spins(optimized_spins, pre_solved_spins), + ) diff --git a/tests/optimizer/test_preprocessing.py b/tests/optimizer/test_preprocessing.py index bff2983b..edc6f57e 100644 --- a/tests/optimizer/test_preprocessing.py +++ b/tests/optimizer/test_preprocessing.py @@ -1,7 +1,7 @@ import pytest import torch -from src.simulated_bifurcation.optimizer import Preprocessing +from src.simulated_bifurcation.optimizer.preprocessing import Preprocessing J = torch.tensor( [ From aac22f47c232d51b0974a5bee4211e1b2a443f12 Mon Sep 17 00:00:00 2001 From: Thomas Bouquet Date: Sun, 31 Dec 2023 10:35:40 +0100 Subject: [PATCH 3/5] fixed wrong sign in projection --- src/simulated_bifurcation/models/abc_model.py | 12 +++++++++--- src/simulated_bifurcation/optimizer/preprocessing.py | 2 +- tests/core/test_ising.py | 2 +- tests/optimizer/test_preprocessing.py | 10 +++++----- 4 files changed, 16 insertions(+), 10 deletions(-) diff --git a/src/simulated_bifurcation/models/abc_model.py b/src/simulated_bifurcation/models/abc_model.py index 690a7407..0a55848d 100644 --- a/src/simulated_bifurcation/models/abc_model.py +++ b/src/simulated_bifurcation/models/abc_model.py @@ -34,7 +34,8 @@ def optimize( use_window: bool = True, sampling_period: int = 50, convergence_threshold: int = 50, - timeout: Optional[float] = None + timeout: Optional[float] = None, + presolve: bool = False, ) -> Tuple[torch.Tensor, torch.Tensor]: return super().optimize( self.domain, @@ -49,6 +50,7 @@ def optimize( sampling_period=sampling_period, convergence_threshold=convergence_threshold, timeout=timeout, + presolve=presolve, ) def minimize( @@ -63,7 +65,8 @@ def minimize( use_window: bool = True, sampling_period: int = 50, convergence_threshold: int = 50, - timeout: Optional[float] = None + timeout: Optional[float] = None, + presolve: bool = False, ) -> Tuple[torch.Tensor, torch.Tensor]: return self.optimize( agents, @@ -77,6 +80,7 @@ def minimize( sampling_period=sampling_period, convergence_threshold=convergence_threshold, timeout=timeout, + presolve=presolve, ) def maximize( @@ -91,7 +95,8 @@ def maximize( use_window: bool = True, sampling_period: int = 50, convergence_threshold: int = 50, - timeout: Optional[float] = None + timeout: Optional[float] = None, + presolve: bool = False, ) -> Tuple[torch.Tensor, torch.Tensor]: return self.optimize( agents, @@ -105,4 +110,5 @@ def maximize( sampling_period=sampling_period, convergence_threshold=convergence_threshold, timeout=timeout, + presolve=presolve, ) diff --git a/src/simulated_bifurcation/optimizer/preprocessing.py b/src/simulated_bifurcation/optimizer/preprocessing.py index c6ca4798..47f5b495 100644 --- a/src/simulated_bifurcation/optimizer/preprocessing.py +++ b/src/simulated_bifurcation/optimizer/preprocessing.py @@ -22,7 +22,7 @@ def _remove_all_coefficients(self, index: int): self._remove_column(index) def _project_coefficients_in_linear_part(self, row_index: int, sign: int): - self.h += sign * self.J[row_index] + self.h -= sign * self.J[row_index] def _project_coefficients_and_delete_row_and_column( self, spin_index: int, spin_sign: int diff --git a/tests/core/test_ising.py b/tests/core/test_ising.py index 6ddf1f06..7a88540d 100644 --- a/tests/core/test_ising.py +++ b/tests/core/test_ising.py @@ -116,5 +116,5 @@ def test_presolve_solves_all_spins(): ising = Ising(presolvable_J, presolvable_h) ising.minimize(agents=3, presolve=True) assert torch.equal( - ising.computed_spins, torch.tensor([[-1, -1, -1], [-1, -1, -1], [1, 1, 1]]) + ising.computed_spins, torch.tensor([[-1, -1, -1], [1, 1, 1], [1, 1, 1]]) ) diff --git a/tests/optimizer/test_preprocessing.py b/tests/optimizer/test_preprocessing.py index edc6f57e..f4a3655f 100644 --- a/tests/optimizer/test_preprocessing.py +++ b/tests/optimizer/test_preprocessing.py @@ -66,10 +66,10 @@ def test_project_coefficients_in_linear_part(): preprocesser = Preprocessing(J, h) preprocesser._project_coefficients_in_linear_part(1, 1) assert torch.equal(preprocesser.J, J) - assert torch.equal(preprocesser.h, torch.tensor([3, 1, 2], dtype=torch.float32)) + assert torch.equal(preprocesser.h, torch.tensor([-1, -1, -6], dtype=torch.float32)) preprocesser._project_coefficients_in_linear_part(0, -1) assert torch.equal(preprocesser.J, J) - assert torch.equal(preprocesser.h, torch.tensor([2, -1, -1], dtype=torch.float32)) + assert torch.equal(preprocesser.h, torch.tensor([0, 1, -3], dtype=torch.float32)) def test_project_coefficients_and_delete_row_and_column(): @@ -78,7 +78,7 @@ def test_project_coefficients_and_delete_row_and_column(): assert torch.equal( preprocesser.J, torch.tensor([[1, 2], [2, 1]], dtype=torch.float32) ) - assert torch.equal(preprocesser.h, torch.tensor([-2, -4], dtype=torch.float32)) + assert torch.equal(preprocesser.h, torch.tensor([4, 4], dtype=torch.float32)) def test_get_optimizable_spins(): @@ -129,7 +129,7 @@ def test_set_first_optimal_spin(): optimizable_preprocesser.J, torch.tensor([[0, 2], [2, 0]], dtype=torch.float32) ) assert torch.equal( - optimizable_preprocesser.h, torch.tensor([0.5, -3], dtype=torch.float32) + optimizable_preprocesser.h, torch.tensor([1.5, -5], dtype=torch.float32) ) assert torch.equal( optimizable_preprocesser.optimized_spins, torch.tensor([-1, 0, 0]) @@ -149,6 +149,6 @@ def test_presolve(): assert torch.equal(non_optimized_h, h) optimizable_preprocesser = Preprocessing(optimizable_J, optimizable_h) optimized_spins, optimized_J, optimized_h = optimizable_preprocesser.presolve() - assert torch.equal(optimized_spins, torch.tensor([-1, -1, 1])) + assert torch.equal(optimized_spins, torch.tensor([-1, 1, 1])) assert torch.equal(optimized_J, torch.tensor([]).reshape(0, 0)) assert torch.equal(optimized_h, torch.tensor([])) From 834eecccee03a2a69285bc282ddc2c91464f2e25 Mon Sep 17 00:00:00 2001 From: Thomas Bouquet Date: Sun, 31 Dec 2023 11:00:42 +0100 Subject: [PATCH 4/5] docstring --- .../optimizer/postprocessing.py | 25 ++++ .../optimizer/preprocessing.py | 133 +++++++++++++++++- tests/optimizer/test_preprocessing.py | 6 +- 3 files changed, 159 insertions(+), 5 deletions(-) diff --git a/src/simulated_bifurcation/optimizer/postprocessing.py b/src/simulated_bifurcation/optimizer/postprocessing.py index c157ae5b..8235d6e9 100644 --- a/src/simulated_bifurcation/optimizer/postprocessing.py +++ b/src/simulated_bifurcation/optimizer/postprocessing.py @@ -2,10 +2,35 @@ class Postprocessing: + """ + Utilility class to reconstrcut the spin agents from presolved spins + on the original Ising model and spins optimized by the Simulated + Bifurcation algorithm on the reduced Ising model. + """ + @staticmethod def reconstruct_spins( optimized_spins: torch.Tensor, pre_solved_spins: torch.Tensor ) -> torch.Tensor: + """ + Reconstruct the spin vectors by merging the presolved spins and + the spins optimized by the Simulated Bifurcation algorithm on the + reduced Ising model. + + Parameters + ---------- + optimized_spins : torch.Tensor + Spins returned by the Simulated Bifurcation algorithm + that correspond to the presolved reduced Ising model. + pre_solved_spins : torch.Tensor + Presolved spins of the original Ising model. + + Returns + ------- + torch.Tensor + Reconstructed spins taking in account both the SB-optimized + spins and the presolved spins. + """ original_dimension = pre_solved_spins.shape[0] agents = optimized_spins.shape[1] reconstructed_spins = torch.zeros(original_dimension, agents, dtype=torch.int32) diff --git a/src/simulated_bifurcation/optimizer/preprocessing.py b/src/simulated_bifurcation/optimizer/preprocessing.py index 47f5b495..ebdab2b3 100644 --- a/src/simulated_bifurcation/optimizer/preprocessing.py +++ b/src/simulated_bifurcation/optimizer/preprocessing.py @@ -4,6 +4,14 @@ class Preprocessing: + """ + Utilility class to presolve Ising models. It identifies spins + the optimal sign of which can be determined before running the + Simulated Bifurcation algorithm in order to only use a reduced + version of the original Ising model, saving computation resources + and time. + """ + def __init__(self, J: torch.Tensor, h: torch.Tensor) -> None: self.J = J.clone() self.h = h.clone() @@ -11,46 +19,153 @@ def __init__(self, J: torch.Tensor, h: torch.Tensor) -> None: self.shifted_indices = list(range(self.J.shape[0])) def _remove_row(self, index: int): + """ + Remove a given row from both the matrix J and the + vector h. + + Parameters + ---------- + index : int + Index of the row to remove. + """ self.J = torch.cat((self.J[:index], self.J[index + 1 :]), dim=0) self.h = torch.cat((self.h[:index], self.h[index + 1 :])) def _remove_column(self, index: int): + """ + Remove a given column from the matrix J. + + Parameters + ---------- + index : int + Index of the column to remove. + """ self.J = torch.cat((self.J[:, :index], self.J[:, index + 1 :]), dim=1) def _remove_all_coefficients(self, index: int): + """ + Remove the row and the column with the same given index + from both the matrix J and the vector h. + + Parameters + ---------- + index : int + Index of the row and column to remove. + """ self._remove_row(index) self._remove_column(index) def _project_coefficients_in_linear_part(self, row_index: int, sign: int): + """ + Project a row from the matrix J into the vector h. + + Parameters + ---------- + row_index : int + Index of the row to project. + sign : int + Sign of the spin that acts as a projection coefficient. + """ self.h -= sign * self.J[row_index] def _project_coefficients_and_delete_row_and_column( self, spin_index: int, spin_sign: int ): + """ + Project quadratic coefficients in the vector h and remove + the associated row and column for a presolved optimized spin. + + Parameters + ---------- + spin_index : int + Index of the optimized spin. + spin_sign : int + Sign of the optimized spin. + """ self._project_coefficients_in_linear_part(spin_index, spin_sign) self._remove_all_coefficients(spin_index) def _get_optimizable_spins(self) -> torch.Tensor: + """ + Identify all presolvable spins. + + Returns + ------- + torch.Tensor + Boolean tensor that indicates which spins are presolvable. + """ return torch.abs(self.h) >= torch.sum(torch.abs(self.J), dim=0) def _get_first_optimizable_spin(self) -> Optional[int]: + """ + Get the index of the first presolvable spin. + + Returns + ------- + Optional[int] + Index of the first presolvable spin. If no spin was found, + None is returned instead. + """ optimizable_spins = self._get_optimizable_spins() if torch.any(optimizable_spins).item(): return torch.nonzero(optimizable_spins)[0].item() return None - def _delete_index(self, index: int): + def _drop_index(self, index: int): + """ + Remove an index from the shifted indices buffer. + + Parameters + ---------- + index : int + Index to drop. + """ self.shifted_indices = ( self.shifted_indices[:index] + self.shifted_indices[index + 1 :] ) def _get_original_index(self, new_index: int) -> int: + """ + Retrieve the original index of a given index after + indeces have been shifted because of dropped rows + and columns. + + Parameters + ---------- + new_index : int + Index in the reduced model. + + Returns + ------- + int + Original index. + """ return self.shifted_indices[new_index] def _set_spin_value(self, spin_index: int, spin_sign: int): + """ + Set the sign of an optimal spin. + + Parameters + ---------- + spin_index : int + Index of the optimal spin to set. + spin_sign : int + Sign of the optimal spin to set. + """ self.optimized_spins[spin_index] = spin_sign def _set_first_optimal_spin(self) -> bool: + """ + Get the first presolvable spin and set its value. + The coefficients in the tensors are updated and the associated + row and columns are dropped. + + Returns + ------- + bool + Whether a spin has been set or not. + """ optimal_spin_index = self._get_first_optimizable_spin() if optimal_spin_index is None: return False @@ -58,10 +173,24 @@ def _set_first_optimal_spin(self) -> bool: original_index = self._get_original_index(optimal_spin_index) self._set_spin_value(original_index, sign) self._project_coefficients_and_delete_row_and_column(optimal_spin_index, sign) - self._delete_index(optimal_spin_index) + self._drop_index(optimal_spin_index) return True def presolve(self) -> Tuple[torch.Tensor, torch.Tensor, torch.Tensor]: + """ + Presolves an Ising model. + + Returns + ------- + optimized_spins : torch.Tensor + The presolved spins. + reduced_J : torch.Tensor + The reduced version of the original J matrix. Only the coefficients + associated to not presolved spins are remaining. + reduced_h : torch.Tensor + The reduced version of the original h vector. Only the coefficients + associated to not presolved spins are remaining. + """ _continue = True while _continue: optimized_spin = self._set_first_optimal_spin() diff --git a/tests/optimizer/test_preprocessing.py b/tests/optimizer/test_preprocessing.py index f4a3655f..fa9a5f15 100644 --- a/tests/optimizer/test_preprocessing.py +++ b/tests/optimizer/test_preprocessing.py @@ -98,15 +98,15 @@ def test_get_first_optimizable_spin(): def test_delete_index(): preprocesser = Preprocessing(J, h) assert [0, 1, 2] == preprocesser.shifted_indices - preprocesser._delete_index(1) + preprocesser._drop_index(1) assert [0, 2] == preprocesser.shifted_indices - preprocesser._delete_index(0) + preprocesser._drop_index(0) assert [2] == preprocesser.shifted_indices def test_get_original_index(): preprocesser = Preprocessing(J, h) - preprocesser._delete_index(1) + preprocesser._drop_index(1) assert 0 == preprocesser._get_original_index(0) assert 2 == preprocesser._get_original_index(1) From 5f2a57f7d9e83eba44270b362b3c5ab57e87689c Mon Sep 17 00:00:00 2001 From: Thomas Bouquet Date: Thu, 8 Feb 2024 21:02:32 +0100 Subject: [PATCH 5/5] lint black 24.1.1 --- src/simulated_bifurcation/__init__.py | 1 - src/simulated_bifurcation/core/__init__.py | 1 - src/simulated_bifurcation/core/ising.py | 2 -- src/simulated_bifurcation/models/ising.py | 1 - src/simulated_bifurcation/models/markowitz.py | 2 -- .../models/number_partitioning.py | 1 - src/simulated_bifurcation/models/qubo.py | 1 - src/simulated_bifurcation/optimizer/__init__.py | 1 + .../optimizer/simulated_bifurcation_engine.py | 1 - .../optimizer/simulated_bifurcation_optimizer.py | 1 - src/simulated_bifurcation/optimizer/stop_window.py | 1 - .../optimizer/symplectic_integrator.py | 1 - tests/core/test_quadratic_polynomial.py | 11 +---------- tests/polynomial/test_polynomial_map.py | 11 +---------- 14 files changed, 3 insertions(+), 33 deletions(-) diff --git a/src/simulated_bifurcation/__init__.py b/src/simulated_bifurcation/__init__.py index 1c01e114..7721b03b 100644 --- a/src/simulated_bifurcation/__init__.py +++ b/src/simulated_bifurcation/__init__.py @@ -175,7 +175,6 @@ """ - from . import models from .core import Ising, QuadraticPolynomial from .optimizer import ConvergenceWarning, get_env, reset_env, set_env diff --git a/src/simulated_bifurcation/core/__init__.py b/src/simulated_bifurcation/core/__init__.py index 247bf638..9b95ba25 100644 --- a/src/simulated_bifurcation/core/__init__.py +++ b/src/simulated_bifurcation/core/__init__.py @@ -22,7 +22,6 @@ """ - from .ising import Ising from .quadratic_polynomial import ( PolynomialLike, diff --git a/src/simulated_bifurcation/core/ising.py b/src/simulated_bifurcation/core/ising.py index 760bfa23..a06df2e7 100644 --- a/src/simulated_bifurcation/core/ising.py +++ b/src/simulated_bifurcation/core/ising.py @@ -18,7 +18,6 @@ """ - from typing import Optional, TypeVar, Union import torch @@ -36,7 +35,6 @@ class Ising: - """ Internal implementation of the Ising model. diff --git a/src/simulated_bifurcation/models/ising.py b/src/simulated_bifurcation/models/ising.py index 56f3e4ce..62399a34 100644 --- a/src/simulated_bifurcation/models/ising.py +++ b/src/simulated_bifurcation/models/ising.py @@ -7,7 +7,6 @@ class Ising(ABCModel): - """ Implementation of the Ising model. diff --git a/src/simulated_bifurcation/models/markowitz.py b/src/simulated_bifurcation/models/markowitz.py index 2fa3cf32..80425bc3 100644 --- a/src/simulated_bifurcation/models/markowitz.py +++ b/src/simulated_bifurcation/models/markowitz.py @@ -7,7 +7,6 @@ class SequentialMarkowitz(ABCModel): - """ Implementation of the Markowitz model for the integer trading trajectory optimization problem. @@ -154,7 +153,6 @@ def gains(self) -> float: class Markowitz(SequentialMarkowitz): - """ A representation of the Markowitz model for portfolio optimization. Portfolio only takes integer stocks. diff --git a/src/simulated_bifurcation/models/number_partitioning.py b/src/simulated_bifurcation/models/number_partitioning.py index 3ad3d1f2..d15e05cf 100644 --- a/src/simulated_bifurcation/models/number_partitioning.py +++ b/src/simulated_bifurcation/models/number_partitioning.py @@ -7,7 +7,6 @@ class NumberPartitioning(ABCModel): - """ A solver that separates a set of numbers into two subsets, the respective sums of which are as close as possible. diff --git a/src/simulated_bifurcation/models/qubo.py b/src/simulated_bifurcation/models/qubo.py index faa69f22..fe5a3aaa 100644 --- a/src/simulated_bifurcation/models/qubo.py +++ b/src/simulated_bifurcation/models/qubo.py @@ -7,7 +7,6 @@ class QUBO(ABCModel): - """ Quadratic Unconstrained Binary Optimization diff --git a/src/simulated_bifurcation/optimizer/__init__.py b/src/simulated_bifurcation/optimizer/__init__.py index fd29bcfb..3d0534fd 100644 --- a/src/simulated_bifurcation/optimizer/__init__.py +++ b/src/simulated_bifurcation/optimizer/__init__.py @@ -13,6 +13,7 @@ optimization problems. """ + from .environment import get_env, reset_env, set_env from .postprocessing import Postprocessing from .preprocessing import Preprocessing diff --git a/src/simulated_bifurcation/optimizer/simulated_bifurcation_engine.py b/src/simulated_bifurcation/optimizer/simulated_bifurcation_engine.py index 1d2ab5e4..232860f7 100644 --- a/src/simulated_bifurcation/optimizer/simulated_bifurcation_engine.py +++ b/src/simulated_bifurcation/optimizer/simulated_bifurcation_engine.py @@ -4,7 +4,6 @@ class SimulatedBifurcationEngine(Enum): - """ Enum class that gathers the 4 variants of the Simulated Bifurcation algorithm: diff --git a/src/simulated_bifurcation/optimizer/simulated_bifurcation_optimizer.py b/src/simulated_bifurcation/optimizer/simulated_bifurcation_optimizer.py index 4acee27d..3df26b8b 100644 --- a/src/simulated_bifurcation/optimizer/simulated_bifurcation_optimizer.py +++ b/src/simulated_bifurcation/optimizer/simulated_bifurcation_optimizer.py @@ -27,7 +27,6 @@ def __str__(self) -> str: class SimulatedBifurcationOptimizer: - """ The Simulated Bifurcation (SB) algorithm relies on Hamiltonian/quantum mechanics to find local minima of diff --git a/src/simulated_bifurcation/optimizer/stop_window.py b/src/simulated_bifurcation/optimizer/stop_window.py index 1f8ffa61..29edc62a 100644 --- a/src/simulated_bifurcation/optimizer/stop_window.py +++ b/src/simulated_bifurcation/optimizer/stop_window.py @@ -5,7 +5,6 @@ class StopWindow: - """ Optimization tool to monitor spins bifurcation and convergence for the Simulated Bifurcation (SB) algorithm. diff --git a/src/simulated_bifurcation/optimizer/symplectic_integrator.py b/src/simulated_bifurcation/optimizer/symplectic_integrator.py index 90e03849..c3b4e4c6 100644 --- a/src/simulated_bifurcation/optimizer/symplectic_integrator.py +++ b/src/simulated_bifurcation/optimizer/symplectic_integrator.py @@ -4,7 +4,6 @@ class SymplecticIntegrator: - """ Simulates the evolution of spins' momentum and position following the Hamiltonian quantum mechanics equations that drive the diff --git a/tests/core/test_quadratic_polynomial.py b/tests/core/test_quadratic_polynomial.py index 9c718dcf..193810f2 100644 --- a/tests/core/test_quadratic_polynomial.py +++ b/tests/core/test_quadratic_polynomial.py @@ -20,16 +20,7 @@ def test_build_polynomial_from_expression(): x, y, z = symbols("x y z") expression = poly( - x**2 - + 2 * y**2 - + 3 * z**2 - - 2 * x * y - - x * z - - 3 * y * z - - x - - 2 * y - + z - + 2 + x**2 + 2 * y**2 + 3 * z**2 - 2 * x * y - x * z - 3 * y * z - x - 2 * y + z + 2 ) # Valid definitions diff --git a/tests/polynomial/test_polynomial_map.py b/tests/polynomial/test_polynomial_map.py index 40b58190..bd29af51 100644 --- a/tests/polynomial/test_polynomial_map.py +++ b/tests/polynomial/test_polynomial_map.py @@ -167,16 +167,7 @@ def test_init_polynomial_map_with_inconsistent_dimension_raises_error(): def test_init_polynomial_map_from_expression(): x, y, z = symbols("x y z") expression = poly( - x**2 - + 2 * y**2 - + 3 * z**2 - - 2 * x * y - - x * z - - 3 * y * z - - x - - 2 * y - + z - + 2 + x**2 + 2 * y**2 + 3 * z**2 - 2 * x * y - x * z - 3 * y * z - x - 2 * y + z + 2 ) polynomial_map = PolynomialMap.from_expression(expression) assert_expected_polynomial_map(polynomial_map)