New models API and working knapsack problem

docs/modules/models/knapsack.md

@@ -1,5 +1,11 @@
 # Knapsack problem
 
+$$\text{Minimize } - \sum_{i = 1}^{N} c_{i}x_{i} + A_{\text{totalWeight}} \left ( \sum_{i = 1}^{N} w_{i}x_{i} - \sum_{k = 1}^{W} ky_{k} \right )^{2} + A_{\text{onlyOneWeight}} \left ( 1 - \sum_{k = 0}^{W} y_{k} \right )^{2}$$
+
+with:
+- $A_{\text{totalWeight}} = sum_{i = 1}^{N} c_{i} + 1$
+- $A_{\text{onlyOneWeight}} = \left ( sum_{i = 1}^{N} c_{i} + 1 \right ) \left ( W + sum_{i = 1}^{N} w_{i} + 1 \right)$
+
 ```{eval-rst}
 .. autoclass:: simulated_bifurcation.models.Knapsack
     :members:

src/simulated_bifurcation/models/__init__.py

@@ -3,14 +3,12 @@
 solved using the Simulated Bifurcation algorithm. This module is a
 pre-defined models library to extend the use of SB.
 
-It also comes with an abstract class called `ABCModel` which acts as a
+It also comes with an abstract class called `ABCQuadraticModel` which acts as a
 basis for defining custom optimization models that are not implemented
 yet.
 
 """
 
-from .abc_model import ABCModel
-from .ising import Ising
-from .markowitz import Markowitz, SequentialMarkowitz
-from .number_partitioning import NumberPartitioning
-from .qubo import QUBO
+from .knapsack import knapsack
+from .markowitz import markowitz, sequential_markowitz
+from .number_partitioning import number_partitioning

src/simulated_bifurcation/models/abc_quadratic_model.py

@@ -0,0 +1,60 @@
+from abc import ABC, abstractmethod
+from typing import Any, Literal, Optional, Union
+
+import torch
+
+from ..core.quadratic_polynomial import QuadraticPolynomial
+
+
+class ABCQuadraticModel(ABC):
+    domain: str
+    sense: str
+
+    @abstractmethod
+    def _as_quadratic_polynomial(
+        self, dtype: Optional[torch.dtype], device: Optional[Union[str, torch.device]]
+    ) -> QuadraticPolynomial:
+        raise NotImplementedError() # pragma: no cover
+
+    @abstractmethod
+    def _from_optimized_tensor(
+        self, optimized_tensor: torch.Tensor, optimized_cost: torch.Tensor
+    ) -> Any:
+        raise NotImplementedError() # pragma: no cover
+
+    def solve(
+        self,
+        dtype: Optional[torch.dtype] = None,
+        device: Optional[Union[str, torch.device]] = None,
+        agents: int = 128,
+        max_steps: int = 10_000,
+        mode: Literal["ballistic", "discrete"] = "ballistic",
+        heated: bool = False,
+        verbose: bool = True,
+        early_stopping: bool = True,
+        sampling_period: int = 50,
+        convergence_threshold: int = 50,
+        timeout: Optional[float] = None,
+    ) -> Any:
+        quadratic_model = self._as_quadratic_polynomial(dtype=dtype, device=device)
+        kwargs = {
+            "domain": self.domain,
+            "agents": agents,
+            "max_steps": max_steps,
+            "mode": mode,
+            "heated": heated,
+            "best_only": True,
+            "verbose": verbose,
+            "early_stopping": early_stopping,
+            "sampling_period": sampling_period,
+            "convergence_threshold": convergence_threshold,
+            "timeout": timeout,
+        }
+        if self.sense == "minimize":
+            return self._from_optimized_tensor(*quadratic_model.minimize(**kwargs))
+        elif self.sense == "maximize":
+            return self._from_optimized_tensor(*quadratic_model.maximize(**kwargs))
+        else: # pragma: no cover
+            raise ValueError(
+                f'Unknown optimization sense {self.sense}. Expected "maximize" or "minimize".'
+            )

src/simulated_bifurcation/models/knapsack.py

@@ -0,0 +1,88 @@
+from typing import List, Optional, Union
+
+import numpy as np
+import torch
+
+from ..core.quadratic_polynomial import QuadraticPolynomial
+from .abc_quadratic_model import ABCQuadraticModel
+
+
+class Knapsack(ABCQuadraticModel):
+    """
+    Note
+    ----
+    It is advised to run `solve` in discrete mode with `torch.float32` dtype for an optimal behavior.
+    """
+
+    domain = "binary"
+    sense = "minimize"
+
+    def __init__(
+        self, costs: List[Union[int, float]], weights: List[int], max_weight: int
+    ):
+        self.__costs = list(map(float, costs))
+        self.__weights = list(map(float, weights))
+        self.__max_weight = max_weight
+
+    def _as_quadratic_polynomial(
+        self, dtype: Optional[torch.dtype], device: Optional[Union[str, torch.device]]
+    ) -> QuadraticPolynomial:
+        # Define penalties for quadratic constraints
+        weight_constraint_penalty = float(np.sum(self.__costs) + 1)
+        only_one_weight_constraint_penalty = weight_constraint_penalty * (
+            self.__max_weight + 1 + float(np.sum(self.__weights))
+        )
+
+        # Define arrays and tensors
+        costs_array = np.array(self.__costs).reshape(-1, 1)
+        weights_array = np.array(self.__weights).reshape(-1, 1)
+        possible_total_weights_array = np.arange(self.__max_weight + 1).reshape(-1, 1)
+
+        quadratic = np.block(
+            [
+                [
+                    weight_constraint_penalty * weights_array @ weights_array.T,
+                    -weight_constraint_penalty
+                    * weights_array
+                    @ possible_total_weights_array.T,
+                ],
+                [
+                    -weight_constraint_penalty
+                    * possible_total_weights_array
+                    @ weights_array.T,
+                    weight_constraint_penalty
+                    * (possible_total_weights_array @ possible_total_weights_array.T)
+                    + only_one_weight_constraint_penalty,
+                ],
+            ]
+        )
+
+        linear = np.concat(
+            (
+                -costs_array,
+                -2
+                * only_one_weight_constraint_penalty
+                * np.ones(self.__max_weight + 1).reshape(-1, 1),
+            )
+        ).reshape(
+            -1,
+        )
+
+        bias = only_one_weight_constraint_penalty
+
+        return QuadraticPolynomial(quadratic, linear, bias, dtype=dtype, device=device)
+
+    def _from_optimized_tensor(
+        self, optimized_tensor: torch.Tensor, optimized_cost: torch.Tensor
+    ) -> List[int]:
+        return [
+            index
+            for index in range(len(self.__costs))
+            if optimized_tensor[index].item() == 1.0
+        ]
+
+
+def knapsack(
+    costs: List[Union[int, float]], weights: List[int], max_weight: int
+) -> Knapsack:
+    return Knapsack(costs, weights, max_weight)