bqth29 / simulated-bifurcation-algorithm Goto Github PK

Python CPU/GPU implementation of the Simulated Bifurcation (SB) algorithm to solve quadratic optimization problems (QUBO, Ising, TSP, optimal asset allocations for a portfolio, etc.).

License: MIT License

Python 100.00%

optimal-asset-allocations optimization portfolio python pytorch quantum-algorithms qubo simulated-bifurcation-algorithm

simulated-bifurcation-algorithm's People

Contributors

Stargazers

Watchers

simulated-bifurcation-algorithm's Issues

Error when running on CUDA

Hi Guys
When I tried to run the SB with cuda I get the erro below:
/usr/local/lib/python3.10/dist-packages/simulated_bifurcation/optimizer/stop_window.py in __compare_energies(self, sampled_spins)
111 sampled_spins.t(), sampled_spins.t(), torch.unsqueeze(self.ising_tensor, 0)
112 ).reshape(self.n_agents)
--> 113 torch.eq(
114 energies,
115 self.energies,

RuntimeError: Expected all tensors to be on the same device, but found at least two devices, cuda:0 and cpu!

ValueError: Matrix must be square.

I used this code to perform an opt task and I read on the docs that the An Ising problem, given a null-diagonal square symmetrical matrix.

sb.set_env(time_step=.1, pressure_slope=.01, heat_coefficient=.06)
best_vector, best_value = sb.maximize(matrix, agents=100, device='cuda', 
                                      max_steps=10000, sampling_period=30, ballistic= True, 
                                      convergence_threshold=50, use_window=True, heated=True, `best_only=True)`

The following error appears:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
[<ipython-input-65-6fffbfba281a>](https://localhost:8080/#) in <cell line: 2>()
      1 sb.set_env(time_step=.1, pressure_slope=.01, heat_coefficient=.06)
----> 2 best_vector, best_value = sb.maximize(matrix, agents=100, device='cuda', 
      3                                       max_steps=10000, sampling_period=30, ballistic= True,
      4                                       convergence_threshold=50, use_window=True, heated=True, best_only=True)
      5 

6 frames
[/usr/local/lib/python3.10/dist-packages/simulated_bifurcation/simulated_bifurcation.py](https://localhost:8080/#) in maximize(matrix, vector, constant, domain, dtype, device, agents, max_steps, best_only, ballistic, heated, verbose, use_window, sampling_period, convergence_threshold, timeout, input_type)
    864         domain = input_type
    865 
--> 866     return optimize(
    867         matrix,
    868         vector,

[/usr/local/lib/python3.10/dist-packages/simulated_bifurcation/simulated_bifurcation.py](https://localhost:8080/#) in optimize(matrix, vector, constant, domain, dtype, device, agents, max_steps, best_only, ballistic, heated, minimize, verbose, use_window, sampling_period, convergence_threshold, timeout, input_type)
    315         domain = input_type
    316 
--> 317     model = build_model(
    318         matrix=matrix,
    319         vector=vector,

[/usr/local/lib/python3.10/dist-packages/simulated_bifurcation/simulated_bifurcation.py](https://localhost:8080/#) in build_model(matrix, vector, constant, domain, dtype, device, input_type)
   1043 
   1044     if domain == "spin":
-> 1045         return SpinQuadraticPolynomial(
   1046             matrix=matrix,
   1047             vector=vector,

[/usr/local/lib/python3.10/dist-packages/simulated_bifurcation/polynomial/spin_polynomial.py](https://localhost:8080/#) in __init__(self, matrix, vector, constant, dtype, device, silence_deprecation_warning)
    155             )
    156 
--> 157         super().__init__(
    158             matrix,
    159             vector,

[/usr/local/lib/python3.10/dist-packages/simulated_bifurcation/polynomial/base_multivariate_polynomial.py](https://localhost:8080/#) in __init__(self, matrix, vector, constant, accepted_values, dtype, device, silence_deprecation_warning)
     81             )
     82         self.__check_device(device)
---> 83         self.__init_matrix(matrix, dtype, device)
     84         self.__init_vector(vector, dtype, device)
     85         self.__init_constant(constant, dtype, device)

[/usr/local/lib/python3.10/dist-packages/simulated_bifurcation/polynomial/base_multivariate_polynomial.py](https://localhost:8080/#) in __init_matrix(self, matrix, dtype, device)
    189     ) -> None:
    190         tensor_matrix = self._cast_matrix_to_tensor(matrix, dtype, device)
--> 191         self.__check_square_matrix(tensor_matrix)
    192         self.__matrix = tensor_matrix
    193         self.__dimension = tensor_matrix.shape[0]

[/usr/local/lib/python3.10/dist-packages/simulated_bifurcation/polynomial/base_multivariate_polynomial.py](https://localhost:8080/#) in __check_square_matrix(matrix)
    252             raise ValueError(f"Matrix requires two dimension, got {matrix.ndim}.")
    253         if matrix.shape[0] != matrix.shape[1]:
--> 254             raise ValueError("Matrix must be square.")
    255 
    256     def __check_vector_shape(self, vector: torch.Tensor) -> None:

ValueError: Matrix must be square.

TypeError: 'type' object is not subscriptable when importing simulated_bifurcation

Hi,

When importing simulated_bifurcation, I got en error "TypeError: 'type' object is not subscriptable "

Error message:
class IsingPolynomialInterface(ABC):
11
12 """

simulated_bifurcation/polynomial/ising_polynomial_interface.py in IsingPolynomialInterface()
28 vector: Union[torch.Tensor, np.ndarray, None] = None,
29 constant: Union[int, float, None] = None,
---> 30 accepted_values: Union[torch.Tensor, np.ndarray, list[int], None] = None,
31 dtype: torch.dtype = torch.float32,
32 device: str = "cpu",

Could you help solve this problem? Thanks

[ENH] Replace `logging.warning` with `warnings.warn`

As discussed in #30

It might be more useful to use warnings.warn instead of logging.warning as the former can be caught like an exception and the user can specify some other settings to try that might improve convergence and not require human intervention. Currently, with logging.warning, there's no easy way to programmatically know that the simulation did not converge (there's no way to introspect the agents or the state of the simulator that I am aware of).

Thank you for your consideration!

RuntimeError: expected m1 and m2 to have the same dtype, but got: double != float

Hi, obtaining this error when trying to change the data type to torch.float64 in the minimize() function. Please help! Thanks!

SB Optimizer computation dtype v. Model dtype

Currently, the oscillators in the SB optimizer have the same dtype as the IsingCore model which itself inherits its dtype from the polynomial model defined by the user. Although it makes sense to create a polynomial model with an integer dtype (float32, float64, ...) and to cast the SB results to this integer dtype to allow a full-integer computation, it is counter-productive to use this very dtype for the SB optimization because the oscillators' range of values is [-1, 1] which would not work with integer values.

Thus, it would be nice to allow the user to chose a dtype for the model and a dtype for the optimization.

Several options are availables to remedy this problem:

Option 1: int to float mapping

The dtype provided in the sb.optimize, sb.minimize and sb.maximize functions, is used for the model and the SB computation is derived from it:

if the dtype is a float (float8, float16, float32, float64) it is also used for SB
if the dtype is an integer (int8, int16, int32, int64), SB uses the float dtype encoded on the same number of bits (int8 -> float8, int16 -> float16, etc.)

Option 2: dtype is only for SB computation

The dtype passed is only used for the SB computation (a float dtype is required). If the model to optimize is created first, it can have any dtype, but the equivalent Ising model will have its own dtype. If the polynomial is directly provided in the sb.maximize or sb.minimize function, its dtype will be the SB computation one as well.

Option 3: use two parameters in functions

The optimization functions use 2 parameters: model_dtype and computation_dtype which are respectively used to create the model and run SB

[BUG] MemoryError when creating an instance of IntegerPolynomial for integers with a large number of bits

Description

Creating an instance of IntegerPolynomial over integers with a large number of bits raises a MemoryError.

The IntegerPolynomial.__init__ creates a list of all accepted values which has size 2**number_of_bits.

Code example

import torch
import simulated_bifurcation as sb

sb.build_model(matrix, input_type="int42")

Traceback

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Users\agero\PycharmProjects\simulated-bifurcation-algorithm\.test_env\.pypi_env_no_gpu\lib\site-packages\simulated_bifurcation\simulated_bifurcation.py", line 823, in build_model
    return IntegerPolynomial(
  File "C:\Users\agero\PycharmProjects\simulated-bifurcation-algorithm\.test_env\.pypi_env_no_gpu\lib\site-packages\simulated_bifurcation\polynomial\integer_polynomial.py", line 60, in __init__
    matrix, vector, constant, [*range(2**number_of_bits)], dtype, device
MemoryError

[BUG] Algorithm performance dwindles when the number of agents increases

Considering 6 binary variables $x_1$, $x_2$, $x_3$, $x_4$, $x_5$ and $x_6$, we seek to solve

$$\mbox{Maximize } 2x_1 + 2x_2 + 2x_3 + 2x_4 + 4.5x_5 + 3x_6$$

subjected to the following constraints:

$x_1 + x_2 \leq 1$
$x_1 + x_3 \leq 1$
$x_2 + x_3 + x_4 \leq 1$
$x_3 + x_4 \leq 1$
$x_4 + x_5 \leq 1$
$x_5 + x_6 \leq 1$

Introducing a penalty coefficient $P$, we can write this LP problem as a QUBO problem that we want to maximize:

$$ x_1 + x_2 + x_3 + x_4 - Px_1x_2 - Px_1x_3 - Px_2x_3 - Px_2x_4 - 2Px_3x_4 - Px_4x_5 - Px_5x_6$$

We can chose any value of $P$ as long as one violated constraint results in a negative objective function. Thus, we set $P = 2 + 2 + 2 + 2 + 4.5 + 3 = 15.5$

The optimal situation is met when only $x_1$, $x_4$ and $x_6$ are set to $1$ with an objective function of $7$.

This can be written in terms of PyTorch tensors:

import torch

P = 15.5
Q = torch.tensor(
    [
        [2, -P, -P, 0, 0, 0],
        [0, 2, -P, -P, 0, 0],
        [0, 0, 2, -2 * P, 0, 0],
        [0, 0, 0, 2, -P, 0],
        [0, 0, 0, 0, 4.5, -P],
        [0, 0, 0, 0, 0, 3],
    ]
)

I tried to use the SB algorithm to solve this problem but got an unexpected behavior. The code snippet below allows to reproduce it. In a few words, I used different numbers of agents with 100 runs of the SB algorithm each time and looked at the obtained objective values:

from collections import defaultdict
from simulated_bifurcation.models import QUBO

torch.manual_seed(42) #  Reproducibility
study = {}

for agents in [1, 2, 5, 20, 50, 100]:
    study[agents] = defaultdict(int)
    for _ in range(100):
        _, obj = QUBO(Q).maximize(agents=agents, verbose=False)
        study[agents][obj] += 1

Plotting the study's data as pie charts for a more understandable visualization, I got:

We can clearly see that the optimal value is getting more and more rare as the number of agents increases which is counter-intuitive since using more agents should result in an objective function at least as good as what we get with fewer agents. In the same time, as the number of agents increases, the bad values get more and more frequent.

RuntimeError: expected scalar type Double but found Float

I used the following code to optimize the proposed dataset:

import yfinance as yf
import torch
import random # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
import yfinance as yf
import matplotlib.pyplot as plt


import json
import numpy as np
#from deap import base, creator, tools, algorithms

from datetime import datetime as dt, timedelta as td
#from datetime import datetime
import simulated_bifurcation as sb


asset_name = 'AAPL'

import numpy as np

def generate_weights(data):
    # Get the number of columns in the data
    num_cols = data.shape[1]

    # Generate random weights between 0 and 1
    weights = np.random.rand(num_cols, num_cols)

    # Normalize the weights
    normalized_weights = weights / np.sum(weights)

    return normalized_weights


def generate_normalized_weights(data):
    # Get the number of rows in the data
    num_rows = data.shape[0]

    # Generate random weights between 0 and 1
    weights = np.random.rand(num_rows, 1)

    # Normalize the weights
    normalized_weights = weights / np.sum(weights)
    print(normalized_weights.shape)

    return normalized_weights.flatten()


data = yf.download(tickers=asset_name, period='1y', interval='1d')
data
m_sb = (torch.DoubleTensor(generate_weights(data)))
m_sb = m.double()
m_sb 
sb.set_env(time_step=.1, pressure_slope=.01, heat_coefficient=.06)
best_vector, best_value = sb.maximize(m_sb, #domain='int10',
                                      agents=100, device='cuda',
                                      max_steps=10000, sampling_period=30, ballistic= True,
                                      convergence_threshold=50, use_window=True, heated=True, best_only=True)

I faced the following issue, any suggestions to solve this issue:

🔁 Iterations       :   0%|          | 0/10000 [00:00<?, ? steps/s]
🏁 Bifurcated agents:   0%|          | 0/100 [00:00<?, ? agents/s]
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
[<ipython-input-65-480d0b2ec372>](https://localhost:8080/#) in <cell line: 2>()
      1 sb.set_env(time_step=.1, pressure_slope=.01, heat_coefficient=.06)
----> 2 best_vector, best_value = sb.maximize(matrix, #domain='int10',
      3                                       agents=100, device='cuda',
      4                                       max_steps=10000, sampling_period=30, ballistic= True,
      5                                       convergence_threshold=50, use_window=True, heated=True, best_only=True)

7 frames
[/usr/local/lib/python3.10/dist-packages/simulated_bifurcation/optimizer/stop_window.py](https://localhost:8080/#) in __compare_energies(self, sampled_spins)
    109 
    110     def __compare_energies(self, sampled_spins: torch.Tensor) -> None:
--> 111         energies = torch.nn.functional.bilinear(
    112             sampled_spins.t(), sampled_spins.t(), torch.unsqueeze(self.ising_tensor, 0)
    113         ).reshape(self.n_agents)

RuntimeError: expected scalar type Double but found Float

[ENH] Speed up computation when agents have converged

When the stop window is used, it could be interesting to remove the converged agents to reduce the number of computations and thus speed up the convergence of the other agents.

Currently, the converged agents' spins are stored so their final value is frozen but they are not removed from the oscillators tensors which means that some spins are updated for no reason.

Add GPU tests

Currently, unit tests are only run on CPU but we should add GPU tests to ensure the behavior is the same on CUDA devices.

[ENH] Add Integer Portfolio Optimization

In this paper, they demonstrate the use of SB for solving the portfolio integer optimization problem. I'm sure you are super busy but it would be useful to reproduce this work on different number of assets in order to demonstrate both the quality as well as scalability of this package.

bqth29 / simulated-bifurcation-algorithm Goto Github PK

simulated-bifurcation-algorithm's People

Contributors

Stargazers

Watchers

Forkers

simulated-bifurcation-algorithm's Issues

Description

Code example

Traceback

Recommend Projects

Recommend Topics

Recommend Org