+"""This module implements I/O"""
+from __future__ import annotations
+
+from pathlib import Path
+from typing import TypeVar, cast
+
+import numpy as np
+import pandas as pd
+from numpy.typing import NDArray
+from pandas import DataFrame
+
+DEFAULT_REQUIRED_COLUMN_NAMES = [
+ "asset",
+ "outstanding_now",
+ "min_outstanding_future",
+ "max_outstanding_future",
+ "income_now",
+ "regcap_now",
+]
+
+PortfolioDataT = TypeVar("PortfolioDataT", bound="PortfolioData")
+
+
+[docs]class PortfolioData:
+ """The ``PortfolioData`` stores the data used for portfolio optimization."""
+
+[docs] def __init__(
+ self,
+ portfolio_dataframe: DataFrame,
+ columns_rename: dict[str, str] | None = None,
+ ):
+ """Creates a ``PortfolioData`` object from a pandas ``DataFrame``.
+
+ The portfolio data is expected to contain at least the following columns names:
+
+ - ``"assets"``: The name of the asset.
+ - ``"outstanding_now_now"``: Current outstanding amount per asset.
+ - ``"min_outstanding_future"``: Lower bound outstanding amount in the future per asset.
+ - ``"max_outstanding_future"``: Upper bound outstanding amount in the future per asset.
+ - ``"income_now"``: Current income per asset, corresponds to return multiplied by the current outstanding amount.
+ - ``"regcap_now"``: Current regulatory capital per asset.
+
+ Different column names in the dataset can be used but need to be provided as a
+ renaming dictionary to the ``columns_rename`` argument.
+
+ Args:
+ portfolio_dataframe: Pandas ``DataFrame`` containing the portfolio data.
+ column_rename: to rename columns provided as dict with new column names as keys
+ and to replace column name as value. Example
+ ``{"outstanding_2021": "outstanding_now"}``.
+
+ Raises:
+ ValueError if required columns are not present in dataset.
+ """
+ if columns_rename is not None:
+ portfolio_dataframe.rename(columns=columns_rename, inplace=True)
+
+ # Validate dataset to contain required column names
+ for required_column_name in DEFAULT_REQUIRED_COLUMN_NAMES:
+ if required_column_name not in portfolio_dataframe.columns:
+ raise ValueError(
+ f"Required column name {required_column_name} is not in dataset."
+ )
+
+ self.portfolio_df = portfolio_dataframe
+
+[docs] @classmethod
+ def from_file(
+ cls: type[PortfolioDataT],
+ filename: str | Path,
+ columns_rename: dict[str, str] | None = None,
+ ) -> PortfolioDataT:
+ """Reads portfolio data object into ``PortfolioData``.
+
+ The portfolio data is expected to contain at least the following columns names:
+
+ - ``"assets"``: The name of the asset.
+ - ``"outstanding_now_now"``: Current outstanding amount per asset.
+ - ``"min_outstanding_future"``: Lower bound outstanding amount in the future per asset.
+ - ``"max_outstanding_future"``: Upper bound outstanding amount in the future per asset.
+ - ``"income_now"``: Current income per asset, corresponds to return multiplied by the current outstanding amount.
+ - ``"regcap_now"``: Current regulatory capital per asset.
+
+ Different column names in the dataset can be used but need to be provided as a
+ renaming dictionary to the ``columns_rename`` argument.
+
+ Args:
+ filename: path to portfolio data. If instead ``benchmark_dataset`` is
+ provided, a default benchmark dataset containing 52 assets will be used.
+ column_rename: to rename columns provided as dict with new column names as keys
+ and to replace column name as value. Example
+ ``{"outstanding_2021": "outstanding_now"}``.
+
+ Raises:
+ ValueError if required columns are not present in dataset.
+ """
+ if str(filename) == "benchmark_dataset":
+ filename = Path(__file__).parents[1] / "datasets" / "benchmark_dataset.xlsx"
+
+ filename = str(filename)
+ if filename.endswith(".xlsx"):
+ portfolio_data = pd.read_excel(filename)
+ elif filename.endswith(".csv"):
+ portfolio_data = pd.read_csv(filename)
+ elif filename.endswith(".json"):
+ portfolio_data = pd.read_json(filename)
+ else:
+ raise ValueError("Datatype not supported.")
+ return cls(portfolio_data, columns_rename)
+
+[docs] def get_outstanding_now(self) -> NDArray[np.float_]:
+ """Gets the `outstanding_now` data from the dataset.
+
+ Returns:
+ The `outstanding_now` column from the dataset as a numpy array.
+ """
+ return cast(NDArray[np.float_], self.portfolio_df["outstanding_now"].to_numpy())
+
+[docs] def get_l_bound(self) -> NDArray[np.float_]:
+ """Gets the `l_bound` data from the dataset.
+
+ Returns:
+ The `min_outstanding_future` column from the dataset as a numpy array.
+ """
+ return cast(
+ NDArray[np.float_], self.portfolio_df["min_outstanding_future"].to_numpy()
+ )
+
+[docs] def get_u_bound(self) -> NDArray[np.float_]:
+ """Gets the `u_bound` data from the dataset.
+
+ Returns:
+ The `max_outstanding_future` column from the dataset as a numpy array.
+ """
+ return cast(
+ NDArray[np.float_], self.portfolio_df["max_outstanding_future"].to_numpy()
+ )
+
+[docs] def get_income(self) -> NDArray[np.float_]:
+ """Gets the `income` data from the dataset.
+
+ Returns:
+ The `income_now` column from the dataset as a numpy array.
+ """
+ return cast(NDArray[np.float_], self.portfolio_df["income_now"].to_numpy())
+
+[docs] def get_capital(self) -> NDArray[np.float_]:
+ """Gets the `capital` data from the dataset.
+
+ Returns:
+ The `regcap_now` column from the dataset as a numpy array.
+ """
+ return cast(NDArray[np.float_], self.portfolio_df["regcap_now"].to_numpy())
+
+[docs] def get_returns(self) -> NDArray[np.float_]:
+ """Gets the `returns` data from the dataset.
+
+ Returns:
+ Returns is defined as income / outstanding_now
+ """
+ income = self.get_income()
+ outstanding_now = self.get_outstanding_now()
+ return cast(NDArray[np.float_], income / outstanding_now)
+
+[docs] def get_column(self, column_name: str) -> NDArray[np.float_]:
+ """Gets the specified column from the dataset.
+
+ Args:
+ column_name: Name of the column to get.
+
+ Returns:
+ The `regcap_now` columns from the dataset as a numpy array.
+ """
+ return cast(NDArray[np.float_], self.portfolio_df[column_name].to_numpy())
+
+[docs] def print_portfolio_info(self) -> None:
+ """Prints information about portfolio data to terminal."""
+ outstanding_now = self.get_outstanding_now()
+ l_bound = self.get_l_bound()
+ u_bound = self.get_u_bound()
+ income = self.get_income()
+ capital = self.get_capital()
+
+ # Calculate the total outstanding amount in now
+ total_outstanding_now = np.sum(outstanding_now)
+ print(f"Total outstanding now: {total_outstanding_now:.2f}")
+
+ # Calculate the ROC for now
+ roc_now = np.sum(income) / np.sum(capital)
+ print(f"ROC now: {roc_now:.4f}")
+
+ # Calculate the HHI diversification for now
+ hhi_now = (
+ np.sum(total_outstanding_now**2) / np.sum(total_outstanding_now) ** 2
+ )
+ print(f"HHI now: {hhi_now:.4f}")
+
+ if "emis_intens_now" in self.portfolio_df:
+ # Calculate the total emissions for now
+ total_emission_now = np.sum(
+ self.get_column("emis_intens_now") * total_outstanding_now
+ )
+ print(f"Total Emission now: {total_emission_now:.2f}")
+
+ # Calculate the average emission intensity now
+ relative_total_emission = total_emission_now / total_outstanding_now
+ print(f"Relative emission intensity now: {relative_total_emission:.2f}")
+
+ # Estimate the total outstanding amount and its standard deviation for future. This
+ # follows from the assumption of a symmetric probability distribution on the
+ # interval [l_bound, u_bound] and the central limit theorem.
+ expected_total_outstanding_future = np.sum(u_bound + l_bound) / 2
+ expected_stddev_total_outstanding_future = np.linalg.norm(
+ (u_bound - l_bound) / 2
+ )
+ print(
+ f"Expected total outstanding future: {expected_total_outstanding_future:.2f}",
+ f"Std dev: {expected_stddev_total_outstanding_future:.2f}",
+ )
+
+ # Estimate a average growth factor and its standard deviation for now-future. This
+ # consists of the (averaged) amount per asset in future, which is the outcome of the
+ # optimization, divided by the amount for now.
+ expected_average_growth_fac = np.sum(
+ (u_bound + l_bound) / (2 * total_outstanding_now)
+ )
+ expected_stddev_average_growth_fac = np.linalg.norm(
+ (u_bound - l_bound) / (2 * total_outstanding_now)
+ )
+ print(
+ f"Expected average growth factor: {expected_average_growth_fac:.4f}",
+ f"Std dev: {expected_stddev_average_growth_fac:.4f}",
+ )
+
+
+
+[docs] def __repr__(self) -> str:
+ """Representation for debugging."""
+ txt = f"{self.__class__.__name__} object containing the following data:\n"
+ txt += repr(self.portfolio_df)
+ return txt
+
+[docs] def __str__(self) -> str:
+ """String representation of the ``PortfolioData`` object."""
+ txt = f"{self.__class__.__name__} object containing the following data:\n"
+ txt += str(self.portfolio_df)
+ return txt
+
+[docs] def __contains__(self, other: object) -> bool:
+ """Check if ``other`` is part of the dataset."""
+ return other in self.portfolio_df
+
+"""This module implements required post processing steps."""
+from __future__ import annotations
+
+from collections.abc import Mapping
+
+import numpy as np
+from dimod import SampleSet
+from numpy.typing import ArrayLike, NDArray
+from scipy.spatial import ConvexHull # pylint: disable=no-name-in-module
+
+from tno.quantum.problems.portfolio_optimization.components.io import PortfolioData
+
+
+[docs]class Decoder:
+ """``Decoder`` class for decoding samples and samplesets."""
+
+[docs] def __init__(
+ self,
+ portfolio_data: PortfolioData,
+ k: int,
+ ) -> None:
+ """
+ Init for the ``Decoder`` Class.
+
+ Args:
+ portfolio_data: A ``PortfolioData`` object containing the portfolio to
+ optimize.
+ k: The number of bits that are used to represent the outstanding amount for
+ each asset. A fixed point representation is used to represent `$2^k$`
+ different equidistant values in the range `$[LB_i, UB_i]$` for asset i.
+ """
+ self.number_of_assets = len(portfolio_data)
+ self.k = k
+ self.mantissa = np.power(2, np.arange(self.k))
+ self.l_bound = portfolio_data.get_l_bound()
+ self.u_bound = portfolio_data.get_u_bound()
+ self.multiplier = (self.u_bound - self.l_bound) / (2**self.k - 1)
+
+[docs] def decode_sample(self, sample: Mapping[int, int]) -> NDArray[np.float_]:
+ """Decodes a sample to the `oustanding_future` array.
+
+ Args:
+ sample: Sample as returned by D-Wave.
+
+ Returns:
+ Array containing all `outstanding future` values.
+ """
+ sample_array = np.array(
+ [sample[i] for i in range(self.number_of_assets * self.k)], dtype=np.uint8
+ )
+ sample_reshaped = sample_array.reshape((self.number_of_assets, self.k))
+ ints = np.sum(sample_reshaped * self.mantissa, axis=1)
+ outstanding_future = self.l_bound + self.multiplier * ints
+ if (self.l_bound > outstanding_future).any() or (
+ self.u_bound < outstanding_future
+ ).any():
+ raise ValueError("Bounds not obeyed.")
+
+ return np.asarray(outstanding_future, dtype=np.float_)
+
+[docs] def decode_sampleset(self, sampleset: SampleSet) -> NDArray[np.float_]:
+ """Efficiently decodes a `sampleset` create a matrix of `oustanding_future`
+ values.
+
+ Each row in the matrix corresponds to a different sample in the `sampleset`.
+
+ Args:
+ sampleset: ``SampleSet`` as returned by D-Wave.
+
+ Returns:
+ Matrix containing all `outstanding future` values.
+ """
+ samples_matrix = sampleset.record.sample[:, : self.number_of_assets * self.k]
+ samples_reshaped = samples_matrix.reshape(
+ (len(sampleset), self.number_of_assets, self.k)
+ )
+
+ ints = np.sum(samples_reshaped * self.mantissa, axis=2)
+ outstanding_future = self.l_bound + self.multiplier * ints
+
+ if (self.l_bound > outstanding_future).any() or (
+ self.u_bound < outstanding_future
+ ).any():
+ raise ValueError("Bounds not obeyed.")
+
+ return np.asarray(outstanding_future, dtype=np.float_)
+
+
+[docs]def pareto_front(
+ xvals: ArrayLike,
+ yvals: ArrayLike,
+ min_points: int = 50,
+ upper_right_quadrant: bool = True,
+) -> tuple[NDArray[np.float_], NDArray[np.float_]]:
+ """Calculates the pareto front with at least min_points data points by repeatedly
+ creating a convex hull around data points.
+
+ Args:
+ xvals: x-values of data points
+ yvals: y-values of data points
+ min_points: minimum number of points to be selected
+ upper_right_quadrant: If ``True``, only show the upper right quadrant of the
+ pareto front.
+
+ Returns:
+ x, y values of the points that are on the pareto front
+ """
+ points = np.vstack((xvals, yvals)).T
+ points = np.unique(points, axis=0)
+ if len(points) < 3:
+ return points.T[0], points.T[1]
+
+ hull = ConvexHull(points)
+ pareto_points = points[hull.vertices]
+
+ # Expand the pareto front so that it contains at least min_points.
+ for _ in range(min_points):
+ if len(pareto_points) >= min_points:
+ break
+ # Remove current hull vertices from data and create a new hull
+ points = np.delete(points, hull.vertices, axis=0)
+ if len(points) < 3:
+ break
+ hull = ConvexHull(points)
+ # Add the new hull vertices to the pareto front
+ new_points = points[hull.vertices]
+ pareto_points = np.vstack((pareto_points, new_points))
+ if upper_right_quadrant:
+ pareto_points = _get_upper_quadrant(pareto_points)
+
+ if upper_right_quadrant:
+ pareto_points = _get_upper_quadrant(pareto_points)
+
+ return pareto_points.T[0], pareto_points.T[1]
+
+
+def _get_upper_quadrant(points: NDArray[np.float_]) -> NDArray[np.float_]:
+ """Removes all values that are not in the upper right quadrant of the pareto front."""
+ x_values = points.T[0]
+ y_values = points.T[1]
+
+ x_bound = x_values[np.argmax(y_values)]
+ y_bound = y_values[np.argmax(x_values)]
+ mask = (x_values >= x_bound) & (y_values >= y_bound)
+
+ x_values = x_values[mask]
+ y_values = y_values[mask]
+
+ return np.vstack((x_values, y_values)).T
+
+"""This module contains the ``QuboCompiler`` class.
+
+The ``QuboCompiler`` can create a variety of QUBO formulation by combining different
+objectives and constraints with their corresponding penalty or preference parameters.
+"""
+from __future__ import annotations
+
+from collections.abc import Callable
+from functools import partial
+from typing import TypeVar
+
+import numpy as np
+from numpy.typing import NDArray
+
+from tno.quantum.problems.portfolio_optimization.components.io import PortfolioData
+
+from ._qubo_factory import QuboFactory
+
+QuboCompilerT = TypeVar("QuboCompilerT", bound="QuboCompiler")
+
+
+[docs]class QuboCompiler:
+ """QuboCompiler - A compiler class for creating QUBO instances.
+
+ This class provides a convenient interface for combining different QUBO formulations
+ without needing to worry about the qubo size.
+
+ Methods:
+
+ - `add_minimize_hhi`: Adds the to minimize HHI QUBO to the compile list.
+ - `add_maximize_roc`: Adds a ROC and optionally a stabilizing QUBO to the compile
+ list.
+ - `add_emission_constraint`: Adds an emission constraint QUBO to the compile list.
+ - `add_growth_factor_constraint`: Adds the growth factor constraint QUBO to the
+ compile list.
+
+ """
+
+[docs] def __init__(self, portfolio_data: PortfolioData, k: int) -> None:
+ """Init of the ``QuboCompiler`` class.
+
+ The ``QuboCompiler`` can create a variety of QUBO formulation by combining
+ different objectives and constraints with penalty or preference parameters.
+
+ Args:
+ portfolio_data: A ``PortfolioData`` object containing the portfolio to
+ optimize.
+ k: The number of bits that are used to represent the outstanding amount for
+ each asset. A fixed point representation is used to represent `$2^k$`
+ different equidistant values in the range `$[LB_i, UB_i]$` for asset i.
+ """
+ self._qubo_factory = QuboFactory(portfolio_data, k)
+
+ self._to_compile: list[Callable[[], tuple[NDArray[np.float_], float]]] = []
+ self._compiled_qubos: list[NDArray[np.float_]] = []
+
+[docs] def add_minimize_hhi(
+ self: QuboCompilerT,
+ ) -> QuboCompilerT:
+ r"""Adds the minimize HHI objective to the compile list.
+
+ The HHI objective is given by
+
+ $$HHI(x) = \sum_{i=1}^N\left(\frac{x_i}{\sum_{j=1}^N x_j}\right)^2,$$
+
+ where
+
+ - `$N$` is the total number of assets,
+ - `$x_i$` is the future outstanding amount for asset `$i$`.
+
+ For the QUBO formulation, see the docs of
+ :py:meth:`~portfolio_optimization.components.qubos.QuboFactory.calc_minimize_hhi`.
+
+ Returns:
+ Self.
+ """
+ self._to_compile.append(self._qubo_factory.calc_minimize_hhi)
+ return self
+
+[docs] def add_maximize_roc(
+ self: QuboCompilerT,
+ formulation: int,
+ ancilla_variables: int = 0,
+ ) -> QuboCompilerT:
+ """Adds the maximize ROC objective and based on the chosen formulation a
+ stabilize c constraint.
+
+ Args:
+ formulation: Integer representing which formulation to pick. If formulation
+ is ``1``, then one QUBO term will be added. If formulation is ``2``,
+ then 2 QUBO terms will be added as well, but the argument
+ `ancilla_variables` must be provided.
+ ancilla_variables: Number of ancilla variables to use for formulation ``2``.
+
+ Returns:
+ Self.
+ """
+ if formulation == 1:
+ self._to_compile.append(self._qubo_factory.calc_maximize_roc1)
+ elif formulation == 2:
+ self._qubo_factory.n_vars += ancilla_variables
+ self._to_compile.append(self._qubo_factory.calc_maximize_roc2)
+ self._to_compile.append(self._qubo_factory.calc_stabilize_c)
+
+ return self
+
+[docs] def add_emission_constraint(
+ self: QuboCompilerT,
+ emission_now: str,
+ emission_future: str | None = None,
+ reduction_percentage_target: float = 0.7,
+ ) -> QuboCompilerT:
+ r"""Adds the emission constraint to the compile list.
+
+ The constraint is given by
+
+ .. math::
+
+ \frac{\sum_{i=1}^Nf_i \cdot x_i}{\sum_{i=1}^N x_i}
+ =
+ g_e \frac{\sum_{i=1}^Ne_i \cdot y_i}{\sum_{i=1}^N y_i},
+
+ where:
+
+ - `$x_i$` is the future outstanding amount for asset `$i$`,
+ - `$y_i$` is the current outstanding amount for asset `$i$`,
+ - `$e_i$` is the current emission intensity for asset `$i$`,
+ - `$f_i$` is the expected emission intensity at the future for asset `$i$`,
+ - `$g_e$` is the target value for the relative emission reduction.
+
+ For the QUBO formulation, see the docs of
+ :py:meth:`~portfolio_optimization.components.qubos.QuboFactory.calc_emission_constraint`.
+
+ Args:
+ emission_now: Name of the column in the portfolio dataset corresponding to
+ the variables at current time.
+ emission_future: Name of the column in the portfolio dataset corresponding
+ to the variables at future time. If no value is provided, it is assumed
+ that the value is constant over time, i.e., the variable
+ ``emission_now`` will be used.
+ reduction_percentage_target: target value for reduction percentage amount.
+
+ Returns:
+ Self.
+ """
+ method = partial(
+ self._qubo_factory.calc_emission_constraint,
+ emission_now=emission_now,
+ emission_future=emission_future,
+ reduction_percentage_target=reduction_percentage_target,
+ )
+ self._to_compile.append(method)
+ return self
+
+[docs] def add_growth_factor_constraint(
+ self: QuboCompilerT, growth_target: float
+ ) -> QuboCompilerT:
+ # pylint: disable=line-too-long
+ r"""Adds the capital growth factor constraint to the compile list.
+
+ The constraint is given by
+
+ $$\frac{\sum_{i=1}^N x_i}{\sum_{i=1}^N y_i} = g_c,$$
+
+ where
+
+ - `$N$` is the total number of assets,
+ - `$x_i$` is the future outstanding amount for asset `$i$`,
+ - `$y_i$` is the current outstanding amount for asset `$i$`,
+ - `$g_c$` is the target value for the total growth factor.
+
+ For the QUBO formulation, see the docs of
+ :py:meth:`~portfolio_optimization.components.qubos.QuboFactory.calc_growth_factor_constraint`.
+
+ Args:
+ growth_target: target value for growth factor total outstanding amount.
+
+ Returns:
+ Self.
+ """
+ # pylint: enable=line-too-long
+ method = partial(
+ self._qubo_factory.calc_growth_factor_constraint,
+ growth_target=growth_target,
+ )
+ self._to_compile.append(method)
+ return self
+
+[docs] def compile(self: QuboCompilerT) -> QuboCompilerT:
+ """Compiles all QUBOs in the compile list.
+
+ Returns:
+ Self."""
+ self._compiled_qubos = []
+ for constructor in self._to_compile:
+ qubo, _ = constructor()
+ self._compiled_qubos.append(qubo)
+ return self
+
+[docs] def make_qubo(self, *lambdas: float) -> tuple[NDArray[np.float_], float]:
+ """Makes a QUBO of the entire problem with the given lambdas.
+
+ Args:
+ lambdas: Scaling parameters for each QUBO in the formulation.
+
+ Returns:
+ Tuple containing the QUBO matrix and offset.
+ """
+ if len(lambdas) != len(self._compiled_qubos):
+ raise ValueError(
+ "Number of lambdas does not correspond with the number of Hamiltonians."
+ )
+ qubo = sum(
+ (
+ lambda_i * qubo_i
+ for lambda_i, qubo_i in zip(lambdas, self._compiled_qubos)
+ ),
+ start=np.zeros_like(self._compiled_qubos[0]),
+ )
+
+ return qubo, float("nan")
+
+"""This module contains the ``QuboFactory`` class.
+
+The ``QuboFactory`` class provides a convenient interface for constructing intermediate
+QUBO matrices for different objectives and constraints.
+"""
+from __future__ import annotations
+
+import numpy as np
+from numpy.typing import NDArray
+
+from tno.quantum.problems.portfolio_optimization.components.io import PortfolioData
+
+# pylint: disable=too-many-instance-attributes
+
+
+[docs]class QuboFactory:
+ """QuboFactory - A factory class for creating QUBO instances.
+
+ This class provides a convenient interface for constructing intermediate QUBO
+ matrices for different objectives and constraints.
+
+ Methods:
+
+ - `calc_minimize_hhi`: Calculates the to minimize HHI QUBO
+ - `calc_maximize_roc1`: Calculates the to maximize return on capital QUBO variant 1
+ - `calc_maximize_roc2`: Calculates the to maximize return on capital QUBO variant 2
+ - `calc_emission_constraint`: Calculates the emission constraint QUBO
+ - `calc_growth_factor_constraint`: Calculates the growth factor constraint QUBO
+ - `calc_stabilize_c`: Calculates the constraint QUBO that stabilizes growth factor.
+ """
+
+[docs] def __init__(self, portfolio_data: PortfolioData, k: int) -> None:
+ """Init of the ``QuboFactory``.
+
+ Args:
+ portfolio_data: A ``PortfolioData`` object containing the portfolio to
+ optimize.
+ k: The number of bits that are used to represent the outstanding amount for
+ each asset. A fixed point representation is used to represent `$2^k$`
+ different equidistant values in the range `$[LB_i, UB_i]$` for asset i.
+ """
+ self.portfolio_data = portfolio_data
+ self.number_of_assets = len(portfolio_data)
+ self.n_vars = self.number_of_assets * k
+ self.outstanding_now = portfolio_data.get_outstanding_now()
+ self.l_bound = portfolio_data.get_l_bound()
+ self.u_bound = portfolio_data.get_u_bound()
+ self.income = portfolio_data.get_income()
+ self.returns = portfolio_data.get_returns()
+ self.capital = portfolio_data.get_capital()
+ self.k = k
+
+[docs] def calc_minimize_hhi(self) -> tuple[NDArray[np.float_], float]:
+ r"""Calculates the to minimize HHI QUBO.
+
+ The QUBO formulation is given by
+
+ .. math::
+
+ QUBO(x)
+ =
+ \sum_{i=1}^N\left(LB_i + \frac{UB_i-LB_i}{2^k-1}\sum_{j=0}^{k-1}2^j\cdot x_{i,j}\right)^2,
+
+ where:
+
+ - `$LB_i$` is the lower bound for asset `$i$`,
+ - `$UB_i$` is the upper bound for asset `$i$`,
+ - `$k$` is the number of bits,
+ - and `$x_{i,j}$` are the $k$ binary variables for asset `$i$` with $j<k$.
+
+ Returns:
+ qubo matrix and its offset
+ """
+ qubo = np.zeros((self.n_vars, self.n_vars))
+ offset = np.sum(self.l_bound**2)
+ multiplier = (self.u_bound - self.l_bound) / (2**self.k - 1)
+ for asset_i, (l_bound_i, multiplier_i) in enumerate(
+ zip(self.l_bound, multiplier)
+ ):
+ # For asset i: (LB + multiplier * sum_j 2**j x_j) ** 2
+ for bit_j in range(self.k):
+ idx_1 = asset_i * self.k + bit_j
+
+ # Diagonal elements: 2 * LB * multiplier * sum_j 2**j x_j
+ qubo[idx_1, idx_1] += l_bound_i * multiplier_i * 2 ** (bit_j + 1)
+ qubo[idx_1, idx_1] += multiplier_i**2 * 2 ** (2 * bit_j)
+
+ # Elements: multiplier**2 * (sum_j 2**j x_j) * (sum_j' 2**j' x_j')
+ for bit_j_prime in range(bit_j + 1, self.k):
+ idx_2 = asset_i * self.k + bit_j_prime
+ qubo[idx_1, idx_2] += multiplier_i**2 * 2 ** (
+ bit_j + bit_j_prime + 1
+ )
+
+ qubo = qubo
+ return qubo, offset
+
+[docs] def calc_emission_constraint(
+ self,
+ emission_now: str,
+ emission_future: str | None = None,
+ reduction_percentage_target: float = 0.7,
+ ) -> tuple[NDArray[np.float_], float]:
+ r"""Calculates the emission constraint QUBO for arbitrary target reduction target
+
+ The QUBO formulation is given by
+
+ .. math::
+
+ QUBO(x)
+ &=
+ \left(
+ \sum_{i=1}^N f_i x_i
+ - g_e E
+ \sum_{i=1}^N x_i
+ \right)^2,
+
+ x_i & = LB_i+\frac{UB_i-LB_i}{2^k-1}\sum_{j=0}^{k-1}2^j\cdot x_{i,j},
+
+ E &= \frac{\sum_{i=1}^N e_i \cdot y_i}{\sum_{i=1}^N y_i},
+
+ where:
+
+ - `$LB_i$` is the lower bound for asset `$i$`,
+ - `$UB_i$` is the upper bound for asset `$i$`,
+ - `$k$` is the number of bits,
+ - `$e_i$` is the current emission intensity for asset `$i$`,
+ - `$f_i$` is the expected emission intensity at the future for asset `$i$`,
+ - `$y_i$` is the current outstanding amount for asset `$i$`,
+ - `$g_e$` is the target value for the relative emission reduction,
+ - and `$x_{i,j}$` are the $k$ binary variables for asset `$i$` with $j<k$.
+
+ Args:
+ emission_now: Name of the column in the portfolio dataset corresponding to
+ the variables at current time.
+ emission_future: Name of the column in the portfolio dataset corresponding
+ to the variables at future time. If no value is provided, it is assumed
+ that the value is constant over time, i.e., the variable
+ ``emission_now`` will be used.
+ reduction_percentage_target: Target value for reduction percentage amount.
+
+ Raises:
+ KeyError: if the provided column names are not in the portfolio_data.
+
+ Returns:
+ qubo matrix and its offset
+ """
+ if emission_future is None:
+ emission_future = emission_now
+
+ if emission_now not in self.portfolio_data:
+ raise KeyError(
+ f"Column name {emission_now} not present in portfolio dataset."
+ )
+ if emission_future not in self.portfolio_data:
+ raise KeyError(
+ f"Column name {emission_future} not present in portfolio dataset."
+ )
+
+ emission_intensity_now = self.portfolio_data.get_column(emission_now)
+ emission_intensity_future = self.portfolio_data.get_column(emission_future)
+
+ total_emission_now = np.sum(emission_intensity_now * self.outstanding_now)
+ relelative_total_emission_now = total_emission_now / np.sum(
+ self.outstanding_now
+ )
+
+ alpha = np.sum(
+ (
+ emission_intensity_future
+ - reduction_percentage_target * relelative_total_emission_now
+ )
+ * self.l_bound
+ )
+
+ mantisse = np.power(2, np.arange(self.k))
+ multiplier = (
+ (
+ emission_intensity_future
+ - reduction_percentage_target * relelative_total_emission_now
+ )
+ * (self.u_bound - self.l_bound)
+ / (2**self.k - 1)
+ )
+ beta = np.kron(multiplier, mantisse)
+ qubo = np.zeros((self.n_vars, self.n_vars))
+
+ offset = alpha**2
+ for idx1 in range(self.number_of_assets * self.k):
+ qubo[idx1, idx1] += 2 * alpha * beta[idx1] + beta[idx1] ** 2
+ for idx2 in range(idx1 + 1, self.number_of_assets * self.k):
+ qubo[idx1, idx2] += 2 * beta[idx1] * beta[idx2]
+
+ qubo = qubo
+ return qubo, offset
+
+[docs] def calc_growth_factor_constraint(
+ self, growth_target: float
+ ) -> tuple[NDArray[np.float_], float]:
+ r"""Calculates the growth factor constraint QUBO
+
+ The QUBO formulation is given by
+
+ .. math::
+
+ QUBO(x)
+ =
+ \left(
+ \frac{\sum_{i=1}^N LB_i + \frac{UB_i-LB_i}{2^k-1}\sum_{j=0}^{k-1} 2^j\cdot x_{i,j}}{\sum_{i=1}^N y_i}
+ - g_c
+ \right)^2
+
+ where:
+
+ - `$LB_i$` is the lower bound for asset `$i$`,
+ - `$UB_i$` is the upper bound for asset `$i$`,
+ - `$k$` is the number of bits,
+ - `$g_c$` is the target value for the total growth factor,
+ - `$y_i$` is the current outstanding amount for asset `$i$`,
+ - and `$x_{i,j}$` are the $k$ binary variables for asset `$i$` with $j<k$.
+
+ Args:
+ growth_target: target value for growth factor total outstanding amount.
+
+ Returns:
+ qubo matrix and its offset
+ """
+ total_outstanding_now = np.sum(self.outstanding_now)
+ alpha = np.sum(self.l_bound) / total_outstanding_now - growth_target
+
+ mantisse = np.power(2, np.arange(self.k))
+ multiplier = (self.u_bound - self.l_bound) / (
+ (2**self.k - 1) * total_outstanding_now
+ )
+ beta = np.kron(multiplier, mantisse)
+
+ qubo = np.zeros((self.n_vars, self.n_vars))
+ # We only fill the upper left part of the matrix, since we don't use ancilla
+ # variables
+ qubo_slice = qubo[
+ : self.k * self.number_of_assets, : self.k * self.number_of_assets
+ ]
+ # Add the off diagonal elements
+ qubo_slice += np.triu(2 * np.outer(beta, beta), k=1)
+ # Add the diagonal elements
+ np.fill_diagonal(qubo_slice, beta**2 + 2 * alpha * beta)
+ offset = alpha**2
+ return qubo, float(offset)
+
+[docs] def calc_maximize_roc1(self) -> tuple[NDArray[np.float_], float]:
+ r"""Calculates the to maximize ROC QUBO for variant 1.
+
+ The QUBO formulation is given by
+
+ .. math::
+
+ QUBO(x)
+ =
+ -\sum_{i=1}^N\frac{r_i}{c_i\cdot y_i}
+ \left(LB_i + \frac{UB_i-LB_i}{2^k-1}\sum_{j=0}^{k-1}2^j\cdot x_{i,j}\right),
+
+ where
+
+ - `$LB_i$` is the lower bound for asset `$i$`,
+ - `$UB_i$` is the upper bound for asset `$i$`,
+ - `$k$` is the number of bits,
+ - `$y_i$` is the current outstanding amount for asset `$i$`,
+ - `$r_i$` is the current return for asset `$i$`,
+ - `$c_i$` is the regulatory capital for asset `$i$`,
+ - and `$x_{i,j}$` are the $k$ binary variables for asset `$i$` with $j<k$.
+
+ Returns:
+ qubo matrix and its offset
+ """
+ theta = self.returns / (self.outstanding_now * self.capital)
+ offset = np.sum(theta * self.l_bound)
+ mantisse = np.power(2, np.arange(self.k))
+ multiplier = theta * (self.u_bound - self.l_bound) / (2**self.k - 1)
+ qubo_diag = np.kron(multiplier, mantisse)
+
+ qubo = np.diag(qubo_diag)
+ return -qubo, -offset
+
+[docs] def calc_maximize_roc2(self) -> tuple[NDArray[np.float_], float]:
+ r"""Calculates the to maximize ROC QUBO for variant 2.
+
+ The QUBO formulation is given by
+
+ .. math::
+
+ QUBO(x,g)
+ &=
+ -
+ G_{inv}(g) \cdot
+ \sum_{i=1}^N\frac{r_i}{y_i}
+ \left(LB_i + \frac{UB_i-LB_i}{2^k-1}\sum_{j=0}^{k-1}2^j\cdot x_{i,j}\right),
+
+ G_{inv}(g) &=
+ 1 + \sum_{j=0}^{k-1} 2^{-j-1}(2^{-j-1} - 1)\cdot g_{j},
+
+ where
+
+ - `$LB_i$` is the lower bound for asset `$i$`,
+ - `$UB_i$` is the upper bound for asset `$i$`,
+ - `$k$` is the number of bits,
+ - `$a$` is the number of ancilla variables,
+ - `$y_i$` is the current outstanding amount for asset `$i$`,
+ - `$r_i$` is the return for asset `$i$`,
+ - `$c_i$` is the regulatory capital for asset `$i$`,
+ - `$x_{i,j}$` are the $k$ binary variables for asset `$i$` with $j<k$.
+ - `$g_{j}$` are the $a$ binary ancilla variables with $j<a$.
+
+ Returns:
+ qubo matrix and its offset
+ """
+ ancilla_variables = self.n_vars - self.k * self.number_of_assets
+
+ alpha = np.sum(self.l_bound * self.returns / self.outstanding_now)
+ mantisse = mantisse = np.power(2, np.arange(self.k))
+ multiplier = (
+ self.returns
+ * (self.u_bound - self.l_bound)
+ / (self.outstanding_now * (2**self.k - 1))
+ )
+ beta = np.kron(multiplier, mantisse)
+
+ gamma = np.power(2.0, np.arange(-1, -ancilla_variables - 1, -1))
+ gamma = gamma**2 - gamma
+
+ qubo = np.zeros((self.n_vars, self.n_vars))
+ np.fill_diagonal(
+ qubo[: self.number_of_assets * self.k, : self.number_of_assets * self.k],
+ beta,
+ )
+ qubo[
+ : self.number_of_assets * self.k, self.number_of_assets * self.k :
+ ] += np.outer(beta, gamma)
+ np.fill_diagonal(
+ qubo[self.number_of_assets * self.k :, self.number_of_assets * self.k :],
+ alpha * gamma,
+ )
+ offset = alpha
+ return -qubo, -offset
+
+[docs] def calc_stabilize_c(self) -> tuple[NDArray[np.float_], float]:
+ r"""Calculates the QUBO that stabilizes the growth factor in the second ROC
+ formulation.
+
+ The QUBO formulation is given by
+
+ .. math::
+
+ QUBO(x,g)
+ &=
+ \left(
+ \sum_{i=1}^N\frac{c_i}{y_i}
+ \left(LB_i + \frac{UB_i-LB_i}{2^k-1}\sum_{j=0}^{k-1}2^j\cdot x_{i,j}\right)
+ - G_C(g)\sum_{i=1}^N c_i
+ \right)^2,
+
+ G_C &= 1 + \sum_{j=0}^{k-1} 2^{-j - 1} \cdot g_j,
+
+ where
+
+ - `$LB_i$` is the lower bound for asset `$i$`,
+ - `$UB_i$` is the upper bound for asset `$i$`,
+ - `$k$` is the number of bits,
+ - `$a$` is the number of ancilla variables,
+ - `$y_i$` is the current outstanding amount for asset `$i$`,
+ - `$c_i$` is the regulatory capital for asset `$i$`,
+ - `$x_{i,j}$` are the $k$ binary variables for asset `$i$` with $j<k$,
+ - `$g_j$` are the $a$ ancillary binary variables with $j<a$.
+
+ Returns:
+ qubo matrix and its offset
+ """
+ ancilla_variables = self.n_vars - self.k * self.number_of_assets
+ alpha = np.sum(self.capital * self.l_bound / self.outstanding_now) - np.sum(
+ self.capital
+ )
+
+ mantisse = mantisse = np.power(2, np.arange(self.k))
+ multiplier = (
+ self.capital
+ * (self.u_bound - self.l_bound)
+ / (self.outstanding_now * (2**self.k - 1))
+ )
+ beta = np.kron(multiplier, mantisse)
+
+ gamma = -np.power(2.0, np.arange(-1, -ancilla_variables - 1, -1)) * np.sum(
+ self.capital
+ )
+
+ qubo = np.zeros((self.n_vars, self.n_vars))
+ qubo_upper_left = qubo[
+ : self.number_of_assets * self.k, : self.number_of_assets * self.k
+ ]
+ qubo_upper_left += np.triu(2 * np.outer(beta, beta), k=1)
+ np.fill_diagonal(qubo_upper_left, 2 * alpha * beta + beta**2)
+
+ qubo_upper_right = qubo[
+ : self.number_of_assets * self.k, self.number_of_assets * self.k :
+ ]
+ qubo_upper_right += 2 * np.outer(beta, gamma)
+
+ qubo_lower_right = qubo[
+ self.number_of_assets * self.k :, self.number_of_assets * self.k :
+ ]
+ qubo_lower_right += np.triu(2 * np.outer(gamma, gamma), k=1)
+ np.fill_diagonal(qubo_lower_right, 2 * alpha * gamma + gamma**2)
+ offset = alpha**2
+
+ return qubo, offset
+
+"""This module contains a container for Results object."""
+from __future__ import annotations
+
+import numpy as np
+import pandas as pd
+from numpy.typing import NDArray
+
+from tno.quantum.problems.portfolio_optimization.components.io import PortfolioData
+
+
+[docs]class Results:
+ """Results container"""
+
+[docs] def __init__(
+ self,
+ portfolio_data: PortfolioData,
+ provided_emission_constraints: list[tuple[str, str, float, str]] | None = None,
+ provided_growth_target: float | None = None,
+ ) -> None:
+ """Init of Results container.
+
+ Args:
+ portfolio_data: the portfolio data
+ provided_emission_constraints: list of all the emission constraints that are
+ provided. Each list element contains the ``emission_now``,
+ ``emission_future`` and ``reduction_percentage_target`` input.
+ provided_growth_target: target outstanding amount growth factor if the
+ growth factor constraint is set, otherwise None.
+ """
+ self.portfolio_data = portfolio_data
+ if provided_emission_constraints is None:
+ self.provided_emission_constraints = []
+ else:
+ self.provided_emission_constraints = provided_emission_constraints
+ self.provided_growth_target = provided_growth_target
+
+ self._outstanding_now = portfolio_data.get_outstanding_now()
+ self._total_outstanding_now = np.sum(self._outstanding_now)
+
+ self._returns = portfolio_data.get_income() / self._outstanding_now
+ self._capital = portfolio_data.get_capital()
+ self._roc_now = np.sum(portfolio_data.get_income()) / np.sum(self._capital)
+ self._hhi_now = (
+ np.sum(self._outstanding_now**2) / np.sum(self._outstanding_now) ** 2
+ )
+
+ self.columns = [
+ "outstanding amount",
+ "diff ROC",
+ "diff diversification",
+ "diff outstanding",
+ ] + [
+ "diff " + constraint[3] for constraint in self.provided_emission_constraints
+ ]
+ self.results_df = pd.DataFrame(columns=self.columns)
+
+[docs] def __len__(self) -> int:
+ """Return the number of samples stored in the ``Results`` object."""
+ return len(self.results_df)
+
+[docs] def add_result(self, outstanding_future_samples: NDArray[np.float_]) -> None:
+ """Adds a new outstanding_future data point to results container.
+
+ Args:
+ outstanding_future_samples: outstanding amounts in the future for each
+ sample of the dataset.
+ """
+ for oustanding_future in outstanding_future_samples:
+ total_outstanding_future = np.sum(oustanding_future)
+ # Compute the ROC growth
+ roc = np.sum(oustanding_future * self._returns) / np.sum(
+ oustanding_future * self._capital / self._outstanding_now
+ )
+ roc_growth = 100 * (roc / self._roc_now - 1)
+
+ # Compute the diversification.
+ hhi = np.sum(oustanding_future**2) / total_outstanding_future**2
+ diff_diversification = 100 * (1 - (hhi / self._hhi_now))
+
+ # Compute the growth outstanding in outstanding amount
+ growth_outstanding = total_outstanding_future / self._total_outstanding_now
+
+ new_data = [
+ tuple(oustanding_future),
+ roc_growth,
+ diff_diversification,
+ growth_outstanding,
+ ]
+
+ # Compute the emission constraint growths
+ for (
+ column_name_now,
+ column_name_future,
+ _,
+ _,
+ ) in self.provided_emission_constraints:
+ total_relative_emission_now = (
+ np.sum(
+ self._outstanding_now
+ * self.portfolio_data.get_column(column_name_now)
+ )
+ / self._total_outstanding_now
+ )
+ total_relative_emission_future = np.sum(
+ oustanding_future
+ * self.portfolio_data.get_column(column_name_future)
+ ) / np.sum(oustanding_future)
+
+ new_data.append(
+ 100
+ * (total_relative_emission_future / total_relative_emission_now - 1)
+ )
+
+ # Write results
+ self.results_df.loc[len(self.results_df)] = new_data
+
+[docs] def head(self, n: int = 5) -> pd.DataFrame:
+ """Returns first n rows of self.results_df DataFrame
+
+ Args:
+ selected_columns: By default all columns
+ n: number of results to return
+ """
+ selected_columns = [
+ column for column in self.columns if column != "Outstanding amount"
+ ]
+ return self.results_df[selected_columns].head(n)
+
+[docs] def drop_duplicates(self) -> None:
+ """Drops duplicates in results DataFrame"""
+ self.results_df.drop_duplicates(subset=["outstanding amount"], inplace=True)
+
+[docs] def slice_results(
+ self, tolerance: float = 0.0
+ ) -> tuple[
+ tuple[NDArray[np.float_], NDArray[np.float_]],
+ tuple[NDArray[np.float_], NDArray[np.float_]],
+ ]:
+ """Helper function that slices the results in two groups, those results that
+ satisfy all constraints and those that violate at least one of the growth factor
+ or emission constraints.
+
+ Args:
+ tolerance: tolerance on how strict the constraints need to be satisfied (in
+ percentage point). Example: if the desired target growth rate is 1.2, if
+ the tolerance is set to 0.05 (5%). Solutions that increase outstanding
+ amount by a factor of 1.15 are considered to satisfy the constraints
+ given the tolerance.
+
+ Returns:
+ Relative difference (diversification, roc) coordinates for solutions that
+ satisfy all constraints, and for those that do not satisfy all constraints.
+
+ Raises:
+ ValueError: when there are no emission or growth factor constraints set.
+ """
+ if (
+ self.provided_growth_target is None
+ and len(self.provided_emission_constraints) == 0
+ ):
+ raise ValueError("There are no emission or growth constraints set.")
+
+ mask_growth_target = (
+ True
+ if self.provided_growth_target is None
+ else (
+ self.results_df["diff outstanding"]
+ >= 100 * (self.provided_growth_target - tolerance - 1)
+ )
+ )
+
+ mask_emission_constraint = True
+ for _, _, value, name in self.provided_emission_constraints:
+ mask_emission_constraint &= self.results_df["diff " + name] <= 100 * (
+ (value + tolerance) - 1
+ )
+
+ combined_mask = mask_growth_target & mask_emission_constraint
+
+ # Filter data based on masks
+ filtered_data_met = self.results_df[combined_mask]
+ filtered_data_violated = self.results_df[~combined_mask]
+
+ x_met = filtered_data_met["diff diversification"].to_numpy()
+ y_met = filtered_data_met["diff ROC"].to_numpy()
+ x_violated = filtered_data_violated["diff diversification"].to_numpy()
+ y_violated = filtered_data_violated["diff ROC"].to_numpy()
+ return (x_met, y_met), (x_violated, y_violated)
+
+"""This module contains visualization tools."""
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+import matplotlib.pyplot as plt
+
+from tno.quantum.problems.portfolio_optimization.components.postprocess import (
+ pareto_front,
+)
+
+if TYPE_CHECKING:
+ from collections.abc import Sequence
+ from typing import Any
+
+ from matplotlib.axes import Axes
+ from matplotlib.collections import PatchCollection
+ from matplotlib.colors import Colormap
+ from matplotlib.typing import ColorType
+ from numpy.typing import ArrayLike
+
+
+[docs]def plot_points(
+ diversification_values: ArrayLike,
+ roc_values: ArrayLike,
+ color: str | None = None,
+ label: str | None = None,
+ c: ArrayLike | Sequence[ColorType] | ColorType | None = None,
+ vmin: float | None = None,
+ vmax: float | None = None,
+ alpha: float | None = None,
+ cmap: str | Colormap | None = None,
+ ax: Axes | None = None, # pylint: disable=invalid-name
+) -> PatchCollection | Any:
+ """Plots the given data-points in a Diversification-ROC plot.
+
+ Args:
+ diversification_values: 1-D ``ArrayLike`` containing the x values of the plot.
+ roc_values: 1-D ``ArrayLike`` containing the y values of the plot.
+ color: Optional color to use for the points. For an overview of allowed colors
+ see the `Matplotlib Documentation`_. If ``None`` is given, a default color
+ will be assigned by ``matplotlib``. Default is ``None``.
+ label: Label to use in the legend. If ``None`` is given, no label will be used.
+ Default is ``None``.
+ c: The marker colors as used by ``matplotlib``.
+ vmin: min value of data range that colormap covers as used by ``matplotlib``.
+ vmax: max value of data range that colormap covers as used by ``matplotlib``.
+ alpha: The alpha blending value as used by ``matplotlib``.
+ cmap: The Colormap instance or registered colormap name as used by ``matplotlib``.
+ ax: ``Axes`` to plot on. If ``None``, a new figure with one ``Axes`` will be
+ created.
+
+ Returns:
+ The ``matplotlib`` PathCollection object created by scatter.
+
+ .. _Matplotlib Documentation: https://matplotlib.org/stable/gallery/color/named_colors.html
+ """
+ if ax is None:
+ _, ax = plt.subplots()
+ collection = ax.scatter(
+ diversification_values,
+ roc_values,
+ color=color,
+ c=c,
+ alpha=alpha,
+ vmin=vmin,
+ vmax=vmax,
+ cmap=cmap,
+ label=label,
+ )
+ ax.set_xlabel("Diversification Change (%)")
+ ax.set_ylabel("ROC Change (%)")
+ ax.grid()
+ if label is not None:
+ ax.legend()
+
+ xlim = ax.get_xlim()
+ ylim = ax.get_ylim()
+ ax.axhline(0, color="black", linewidth=1)
+ ax.axvline(0, color="black", linewidth=1)
+ ax.set_xlim(*xlim, auto=True)
+ ax.set_ylim(*ylim, auto=True)
+ return collection
+
+
+[docs]def plot_front(
+ diversification_values: ArrayLike,
+ roc_values: ArrayLike,
+ color: str | None = None,
+ label: str | None = None,
+ c: ArrayLike | Sequence[ColorType] | ColorType | None = None,
+ vmin: float | None = None,
+ vmax: float | None = None,
+ alpha: float | None = None,
+ cmap: str | Colormap | None = None,
+ ax: Axes | None = None, # pylint: disable=invalid-name
+) -> PatchCollection:
+ """Plots a pareto front of the given data-points in a Diversification-ROC plot.
+
+ Args:
+ diversification_values: 1-D ``ArrayLike`` containing the x values of the plot.
+ roc_values: 1-D ``ArrayLike`` containing the y values of the plot.
+ color: Optional color to use for the points. For an overview of allowed colors
+ see the `Matplotlib Documentation`_. If ``None`` is given, a default color
+ will be assigned by ``matplotlib``. Default is ``None``.
+ label: Label to use in the legend. If ``None`` is given, no label will be used.
+ Default is ``None``.
+ c: The marker colors as used by ``matplotlib``.
+ vmin: min value of data range that colormap covers as used by ``matplotlib``.
+ vmax: max value of data range that colormap covers as used by ``matplotlib``.
+ alpha: The alpha blending value as used by ``matplotlib``.
+ cmap: The Colormap instance or registered colormap name as used by ``matplotlib``.
+ ax: ``Axes`` to plot on. If ``None``, a new figure with one ``Axes`` will be
+ created.
+
+ Returns:
+ The ``matplotlib`` PathCollection object created by scatter.
+
+ .. _Matplotlib Documentation: https://matplotlib.org/stable/gallery/color/named_colors.html
+ """
+ x_values, y_values = pareto_front(diversification_values, roc_values)
+ return plot_points(
+ x_values,
+ y_values,
+ color=color,
+ c=c,
+ alpha=alpha,
+ vmin=vmin,
+ vmax=vmax,
+ cmap=cmap,
+ label=label,
+ ax=ax,
+ )
+
+"""This module contains the ``PortfolioOptimizer`` class."""
+from __future__ import annotations
+
+import itertools
+import math
+from datetime import datetime
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+from dimod.core.sampler import Sampler
+from dwave.samplers import SimulatedAnnealingSampler
+from numpy.typing import ArrayLike, NDArray
+from pandas import DataFrame
+from tqdm import tqdm
+
+from tno.quantum.problems.portfolio_optimization.components import (
+ Decoder,
+ PortfolioData,
+ QuboCompiler,
+ Results,
+)
+
+
+[docs]class PortfolioOptimizer:
+ """The ``PortfolioOptimizer`` class is used to convert multi-objective portfolio
+ optimization problems into QUBO problems which can then be solved using QUBO solving
+ techniques such as simulated or quantum annealing.
+
+ The following objectives can be considered
+
+ - `return on capital`, indicated by ROC,
+ - `diversification`, indicated by the `Herfindahl-Hirschman Index`_ HHI.
+
+ The following constraints can be added
+
+ - `capital growth`, demand a minimum increase in outstanding assets.
+ - `emission reduction`, demand a minimum reduction for an arbitrary emission type.
+
+ Usage example:
+
+ .. code-block::
+
+ import numpy as np
+ from dwave.samplers import SimulatedAnnealingSampler
+
+ from tno.quantum.problems.portfolio_optimization import PortfolioOptimizer
+
+ # Choose sampler for solving qubo
+ sampler = SimulatedAnnealingSampler()
+ sampler_kwargs = {"num_reads": 20, "num_sweeps": 200}
+
+ # Set up penalty coefficients for the constraints
+ lambdas1 = np.logspace(-16, 1, 25, endpoint=False, base=10.0)
+ lambdas2 = np.logspace(-16, 1, 25, endpoint=False, base=10.0)
+ lambdas3 = np.array([1])
+
+ # Create portfolio optimization problem
+ portfolio_optimizer = PortfolioOptimizer("benchmark_dataset")
+ portfolio_optimizer.add_minimize_hhi(weights=lambdas1)
+ portfolio_optimizer.add_maximize_roc(formulation=1, weights_roc=lambdas1)
+ portfolio_optimizer.add_emission_constraint(
+ weights=lambdas3,
+ emission_now="emis_intens_now",
+ emission_future="emis_intens_future",
+ name="emission",
+ )
+
+ # Solve the portfolio optimization problem
+ results = portfolio_optimizer.run(sampler, sampler_kwargs)
+ print(results.head())
+
+ .. _Herfindahl-Hirschman Index: https://nl.wikipedia.org/wiki/Herfindahl-index
+ """
+
+[docs] def __init__(
+ self,
+ portfolio_data: PortfolioData | DataFrame | str | Path,
+ k: int = 2,
+ columns_rename: dict[str, str] | None = None,
+ ) -> None:
+ """Init ``PortfolioOptimizer``.
+
+ Args:
+ portfolio_data: Portfolio data represented by a ``PortfolioData`` object,
+ a pandas ``DataFrame`` or a path to where portfolio data is stored. See
+ the docstring of
+ :py:class:`~portfolio_optimization.components.io.PortfolioData` for data
+ input conventions.
+ k: The number of bits that are used to represent the outstanding amount for
+ each asset. A fixed point representation is used to represent `$2^k$`
+ different equidistant values in the range `$[LB_i, UB_i]$` for asset i.
+ column_rename: can be used to rename data columns. See the docstring of
+ :py:class:`~portfolio_optimization.components.io.PortfolioData` for
+ example.
+
+ Raises:
+ TypeError: If the provided ``portfolio_data`` input has the wrong type.
+ """
+ if isinstance(portfolio_data, PortfolioData):
+ self.portfolio_data = portfolio_data
+ elif isinstance(portfolio_data, DataFrame):
+ self.portfolio_data = PortfolioData(portfolio_data, columns_rename)
+ elif isinstance(portfolio_data, (Path, str)):
+ self.portfolio_data = PortfolioData.from_file(
+ portfolio_data, columns_rename
+ )
+ else:
+ raise TypeError(
+ "`portfolio_data` must be of type `PortfolioData`, `DataFrame`, `Path` "
+ f"or `str`, but was of type {type(portfolio_data)}"
+ )
+ self._qubo_compiler = QuboCompiler(self.portfolio_data, k)
+ self.decoder = Decoder(self.portfolio_data, k)
+ self._all_lambdas: list[NDArray[np.float_]] = []
+ self._provided_emission_constraints: list[tuple[str, str, float, str]] = []
+ self._provided_growth_target: float | None = None
+
+[docs] def add_minimize_hhi(self, weights: ArrayLike | None = None) -> None:
+ r"""Adds the minimize HHI objective to the portfolio optimization problem.
+
+ The HHI objective is given by
+
+ $$HHI(x) = \sum_{i=1}^N\left(\frac{x_i}{\sum_{j=1}^N x_j}\right)^2,$$
+
+ where
+
+ - `$N$` is the total number of assets,
+ - `$x_i$` is the future outstanding amount for asset `$i$`.
+
+ As the objective contains non-quadratic terms, a QUBO formulation requires
+ approximations. For the QUBO formulation, see the docs of
+ :py:class:`~portfolio_optimization.components.qubos.QuboFactory`.
+ :py:meth:`~portfolio_optimization.components.qubos.QuboFactory.calc_minimize_hhi`.
+
+ Usage example:
+
+ >>> from tno.quantum.problems.portfolio_optimization import PortfolioOptimizer
+ >>> import numpy as np
+ >>> portfolio_optimizer = PortfolioOptimizer(filename="benchmark_dataset")
+ >>> lambdas = np.logspace(-16, 1, 25, endpoint=False, base=10.0)
+ >>> portfolio_optimizer.add_minimize_hhi(weights=lambdas)
+
+ Args:
+ weights: The coefficients that are considered as penalty parameter.
+ """
+ self._all_lambdas.append(self._parse_weight(weights))
+ self._qubo_compiler.add_minimize_hhi()
+
+[docs] def add_maximize_roc(
+ self,
+ formulation: int,
+ weights_roc: ArrayLike | None = None,
+ ancilla_variables: int = 0,
+ weights_stabilize: ArrayLike | None = None,
+ ) -> None:
+ r"""Adds the maximize ROC objective to the portfolio optimization problem.
+
+ The ROC objective is given by
+
+ .. math::
+
+ ROC(x) =
+ \frac{\sum_{i=1}^N \frac{x_i \cdot r_i}{y_i}}
+ {\sum_{i=1}^N \frac{x_i \cdot c_i}{y_i}},
+
+ where
+
+ - `$N$` is the total number of assets,
+ - `$x_i$` is the future outstanding amount for asset `$i$`,
+ - `$y_i$` is the current outstanding amount for asset `$i$`,
+ - `$r_i$` is the return for asset `$i$`,
+ - `$c_i$` is the regulatory capital for asset `$i$`.
+
+ As the ROC is not a quadratic function, it is approximated using two different
+ formulations:
+
+ formulation 1:
+ $$ROC_1(x)=\sum_{i=1}^N\frac{x_i\cdot r_i}{c_i\cdot y_i}$$
+
+ Adds 1 qubo term, use ``weights_roc`` to scale.
+
+ formulation 2:
+ $$ROC_2(x)=\frac{1}{G_C \cdot C_{21}}\sum_{i=1}^N x_i\frac{r_i}{y_i}$$
+
+ In this formulation, $G_C \cdot C_{21}$ approximates a fixed regulatory
+ capital growth which is equal for all assets, where
+
+ - `$1≤G_C<2$` is a growth factor to be estimated using ancilla variables,
+ - `$C_{21} = \sum_{i=1}^N c_{i}$` is the sum of all assets' regulatory capital.
+
+ This formulation adds 2 qubo terms, one for the ROC term, and one to stabilize the
+ capital growth. The stabilize qubo requires an extra argument ``ancilla_variables``.
+ Use ``weights_roc`` and ``weights_stabilize`` to scale both qubo's accordingly.
+
+ For the different QUBO formulations, see the docs of
+ :py:class:`~portfolio_optimization.components.qubos.QuboFactory`.
+
+ Usage example:
+
+ >>> from tno.quantum.problems.portfolio_optimization import PortfolioOptimizer
+ >>> import numpy as np
+ >>> portfolio_optimizer = PortfolioOptimizer(filename="benchmark_dataset")
+ >>> lambdas = np.logspace(-16, 1, 25, endpoint=False, base=10.0)
+ >>> portfolio_optimizer.add_maximize_roc(...)
+
+ Args:
+ formulation: the ROC QUBO formulation that is being used.
+ Possible options are: [1, 2].
+ weights_roc: The coefficients that are considered as penalty parameter for
+ maximizing the roc objective.
+ ancilla_variables: The number of ancillary variables that are used to
+ represent ``G_C`` using fixed point representation. Only relevant for
+ roc formulation ``2``.
+ weights_stabilize: The coefficients that are considered as penalty parameter
+ for the stabilizing constraint. Only relevant for roc formulation ``2``.
+
+ Raises:
+ ValueError: If invalid formulation is provided.
+ """
+ allowed_formulation_options = [1, 2]
+ if formulation not in allowed_formulation_options:
+ raise ValueError(
+ "Invalid formulation input provided, "
+ f"choose from {allowed_formulation_options}."
+ )
+
+ self._all_lambdas.append(self._parse_weight(weights_roc))
+ if formulation == 2:
+ self._all_lambdas.append(self._parse_weight(weights_stabilize))
+ self._qubo_compiler.add_maximize_roc(formulation, ancilla_variables)
+
+[docs] def add_emission_constraint(
+ self,
+ emission_now: str,
+ emission_future: str | None = None,
+ reduction_percentage_target: float = 0.7,
+ name: str | None = None,
+ weights: ArrayLike | None = None,
+ ) -> None:
+ r"""Adds emission constraint to the portfolio optimization problem.
+
+ The constraint is given by
+
+ .. math::
+
+ \frac{\sum_{i=1}^Nf_i \cdot x_i}{\sum_{i=1}^N x_i}
+ =
+ g_e \frac{\sum_{i=1}^Ne_i \cdot y_i}{\sum_{i=1}^N y_i},
+
+ where:
+
+ - `$N$` is the total number of assets,
+ - `$x_i$` is the future outstanding amount for asset `$i$`,
+ - `$y_i$` is the current outstanding amount for asset `$i$`,
+ - `$e_i$` is the current emission intensity for asset `$i$`,
+ - `$f_i$` is the expected emission intensity at the future for asset `$i$`,
+ - `$g_e$` is the target value for the relative emission reduction.
+
+ Usage example:
+
+ >>> from tno.quantum.problems.portfolio_optimization import PortfolioOptimizer
+ >>> import numpy as np
+ >>> portfolio_optimizer = PortfolioOptimizer(filename="benchmark_dataset")
+ >>> lambdas = np.logspace(-16, 1, 25, endpoint=False, base=10.0)
+ >>> portfolio_optimizer.add_emission_constraint(
+ ... emission_now="emis_intens_now", weights=lambdas
+ ... )
+
+ For the QUBO formulation, see the docs of
+ :py:class:`~portfolio_optimization.components.qubos.QuboFactory`.
+ :py:meth:`~portfolio_optimization.components.qubos.QuboFactory.calc_emission_constraint`.
+
+ Args:
+ emission_now: Name of the column in the portfolio dataset corresponding to
+ the variables emission intensity at current time.
+ emission_future: Name of the column in the portfolio dataset corresponding
+ to the variables emission intensity at future time. If no value is
+ provided, it is assumed that the emission intensity is constant over
+ time, i.e., the variable ``emission_now`` will be used.
+ reduction_percentage_target: target value for reduction percentage amount.
+ name: Name that will be used for emission constraint in the results df.
+ weights: The coefficients that are considered as penalty parameter.
+ """
+ # Store emission constraint information
+ if name is None:
+ name = emission_now
+ if emission_future is None:
+ emission_future = emission_now
+ self._provided_emission_constraints.append(
+ (emission_now, emission_future, reduction_percentage_target, name)
+ )
+ self._all_lambdas.append(self._parse_weight(weights))
+
+ self._qubo_compiler.add_emission_constraint(
+ emission_now=emission_now,
+ emission_future=emission_future,
+ reduction_percentage_target=reduction_percentage_target,
+ )
+
+[docs] def add_growth_factor_constraint(
+ self, growth_target: float, weights: ArrayLike | None = None
+ ) -> None:
+ # pylint: disable=line-too-long
+ r"""Adds an outstanding amount growth factor constraint to the portfolio
+ optimization problem.
+
+ The constraint is given by
+
+ $$\frac{\sum_{i=1}^N x_i}{\sum_{i=1}^N y_i} = g_c,$$
+
+ where
+
+ - `$N$` is the total number of assets,
+ - `$x_i$` is the future outstanding amount for asset `$i$`,
+ - `$y_i$` is the current outstanding amount for asset `$i$`,
+ - `$g_c$` is the target value for the total growth factor.
+
+ This constraint can only be added once.
+
+ Usage example:
+
+ >>> from tno.quantum.problems.portfolio_optimization import PortfolioOptimizer
+ >>> import numpy as np
+ >>> portfolio_optimizer = PortfolioOptimizer(filename="benchmark_dataset")
+ >>> lambdas = np.logspace(-16, 1, 25, endpoint=False, base=10.0)
+ >>> portfolio_optimizer.add_emission_constraint(growth_target=1.2, weights=lambdas)
+
+ For the QUBO formulation, see the docs of
+ :py:class:`~portfolio_optimization.components.qubos.QuboFactory`.
+ :py:meth:`~portfolio_optimization.components.qubos.QuboFactory.calc_growth_factor_constraint`.
+
+ Args:
+ growth_target: target value for growth factor total outstanding amount.
+ weights: The coefficients that are considered as penalty parameter.
+
+ Raises:
+ ValueError: If constraint has been added before.
+ """
+ # pylint: enable=line-too-long
+ if self._provided_growth_target is not None:
+ raise ValueError("Growth factor constraint has been set before.")
+
+ self._all_lambdas.append(self._parse_weight(weights))
+ self._provided_growth_target = growth_target
+ self._qubo_compiler.add_growth_factor_constraint(growth_target)
+
+[docs] def run(
+ self,
+ sampler: Sampler | None = None,
+ sampler_kwargs: dict[str, Any] | None = None,
+ verbose: bool = True,
+ ) -> Results:
+ # pylint: disable=line-too-long
+ """
+ Optimizes a portfolio given the set of provided constraints.
+
+ Usage example:
+
+ >>> from tno.quantum.problems.portfolio_optimization import PortfolioOptimizer
+ >>> from dwave.samplers import SimulatedAnnealingSampler
+ >>> portfolio_optimizer = PortfolioOptimizer(filename="benchmark_dataset")
+ >>> portfolio_optimizer.add_minimize_HHI()
+ >>> portfolio_optimizer.run(sampler=SimulatedAnnealingSampler(), verbose=False)
+
+ Args:
+ sampler: Instance of a D-Wave Sampler that can be used to solve the QUBO.
+ More information can be found in the `D-Wave Ocean Documentation`_.
+ By default the ``SimulatedAnnealingSampler`` is being used.
+ sampler_kwargs: The sampler specific key-word arguments.
+ verbose: If True, print detailed information during execution
+
+ Returns:
+ results
+
+ Raises:
+ ValueError: if constraints are not set
+
+ .. _D-Wave Ocean Documentation: https://docs.ocean.dwavesys.com/projects/system/en/stable/reference/samplers.html
+
+ """
+ # pylint: enable=line-too-long
+ if verbose:
+ self.portfolio_data.print_portfolio_info()
+
+ sampler = SimulatedAnnealingSampler() if sampler is None else sampler
+ sampler_kwargs = {} if sampler_kwargs is None else sampler_kwargs
+
+ if verbose:
+ print("Status: creating model")
+ if self._provided_growth_target is not None:
+ print(f"Growth target: {self._provided_growth_target - 1:.1%}")
+
+ for _, _, target_value, name in self._provided_emission_constraints:
+ print(
+ f"Emission constraint: {name}, "
+ f"target reduction percentage: {target_value - 1:.1%}"
+ )
+
+ self._qubo_compiler.compile()
+
+ results = Results(
+ portfolio_data=self.portfolio_data,
+ provided_emission_constraints=self._provided_emission_constraints,
+ provided_growth_target=self._provided_growth_target,
+ )
+
+ if verbose:
+ print("Status: calculating")
+ starttime = datetime.now()
+
+ total_steps = math.prod(map(len, self._all_lambdas))
+ lambdas_iterator = tqdm(
+ itertools.product(*self._all_lambdas), total=total_steps
+ )
+
+ for lambdas in lambdas_iterator:
+ # Compile the model and generate QUBO
+ qubo, _ = self._qubo_compiler.make_qubo(*lambdas)
+ # Solve the QUBO
+ response = sampler.sample_qubo(qubo, **sampler_kwargs)
+ # Postprocess solution. Iterate over all found solutions. (Compute future portfolios)
+ outstanding_future_samples = self.decoder.decode_sampleset(response)
+ results.add_result(outstanding_future_samples)
+
+ if verbose:
+ print("Drop duplicate samples in results.")
+ results.drop_duplicates()
+
+ if verbose:
+ print("Number of unique samples: ", len(results))
+ print("Time consumed:", datetime.now() - starttime)
+
+ return results
+
+ @staticmethod
+ def _parse_weight(weights: ArrayLike | None = None) -> NDArray[np.float_]:
+ """Converts weights into NumPy array and if needed set default weights to [1.0]
+
+ Args:
+ weights: penalty coefficients.
+
+ Returns:
+ Numpy array of weights
+ """
+ if weights is None:
+ return np.array([1.0])
+ return np.asarray(weights, dtype=np.float_)
+
+"""This module contains visualization tools."""
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+import matplotlib.pyplot as plt
+
+from tno.quantum.problems.portfolio_optimization.components.postprocess import (
+ pareto_front,
+)
+
+if TYPE_CHECKING:
+ from collections.abc import Sequence
+ from typing import Any
+
+ from matplotlib.axes import Axes
+ from matplotlib.collections import PatchCollection
+ from matplotlib.colors import Colormap
+ from matplotlib.typing import ColorType
+ from numpy.typing import ArrayLike
+
+
+[docs]def plot_points(
+ diversification_values: ArrayLike,
+ roc_values: ArrayLike,
+ color: str | None = None,
+ label: str | None = None,
+ c: ArrayLike | Sequence[ColorType] | ColorType | None = None,
+ vmin: float | None = None,
+ vmax: float | None = None,
+ alpha: float | None = None,
+ cmap: str | Colormap | None = None,
+ ax: Axes | None = None, # pylint: disable=invalid-name
+) -> PatchCollection | Any:
+ """Plots the given data-points in a Diversification-ROC plot.
+
+ Args:
+ diversification_values: 1-D ``ArrayLike`` containing the x values of the plot.
+ roc_values: 1-D ``ArrayLike`` containing the y values of the plot.
+ color: Optional color to use for the points. For an overview of allowed colors
+ see the `Matplotlib Documentation`_. If ``None`` is given, a default color
+ will be assigned by ``matplotlib``. Default is ``None``.
+ label: Label to use in the legend. If ``None`` is given, no label will be used.
+ Default is ``None``.
+ c: The marker colors as used by ``matplotlib``.
+ vmin: min value of data range that colormap covers as used by ``matplotlib``.
+ vmax: max value of data range that colormap covers as used by ``matplotlib``.
+ alpha: The alpha blending value as used by ``matplotlib``.
+ cmap: The Colormap instance or registered colormap name as used by ``matplotlib``.
+ ax: ``Axes`` to plot on. If ``None``, a new figure with one ``Axes`` will be
+ created.
+
+ Returns:
+ The ``matplotlib`` PathCollection object created by scatter.
+
+ .. _Matplotlib Documentation: https://matplotlib.org/stable/gallery/color/named_colors.html
+ """
+ if ax is None:
+ _, ax = plt.subplots()
+ collection = ax.scatter(
+ diversification_values,
+ roc_values,
+ color=color,
+ c=c,
+ alpha=alpha,
+ vmin=vmin,
+ vmax=vmax,
+ cmap=cmap,
+ label=label,
+ )
+ ax.set_xlabel("Diversification Change (%)")
+ ax.set_ylabel("ROC Change (%)")
+ ax.grid()
+ if label is not None:
+ ax.legend()
+
+ xlim = ax.get_xlim()
+ ylim = ax.get_ylim()
+ ax.axhline(0, color="black", linewidth=1)
+ ax.axvline(0, color="black", linewidth=1)
+ ax.set_xlim(*xlim, auto=True)
+ ax.set_ylim(*ylim, auto=True)
+ return collection
+
+
+[docs]def plot_front(
+ diversification_values: ArrayLike,
+ roc_values: ArrayLike,
+ color: str | None = None,
+ label: str | None = None,
+ c: ArrayLike | Sequence[ColorType] | ColorType | None = None,
+ vmin: float | None = None,
+ vmax: float | None = None,
+ alpha: float | None = None,
+ cmap: str | Colormap | None = None,
+ ax: Axes | None = None, # pylint: disable=invalid-name
+) -> PatchCollection:
+ """Plots a pareto front of the given data-points in a Diversification-ROC plot.
+
+ Args:
+ diversification_values: 1-D ``ArrayLike`` containing the x values of the plot.
+ roc_values: 1-D ``ArrayLike`` containing the y values of the plot.
+ color: Optional color to use for the points. For an overview of allowed colors
+ see the `Matplotlib Documentation`_. If ``None`` is given, a default color
+ will be assigned by ``matplotlib``. Default is ``None``.
+ label: Label to use in the legend. If ``None`` is given, no label will be used.
+ Default is ``None``.
+ c: The marker colors as used by ``matplotlib``.
+ vmin: min value of data range that colormap covers as used by ``matplotlib``.
+ vmax: max value of data range that colormap covers as used by ``matplotlib``.
+ alpha: The alpha blending value as used by ``matplotlib``.
+ cmap: The Colormap instance or registered colormap name as used by ``matplotlib``.
+ ax: ``Axes`` to plot on. If ``None``, a new figure with one ``Axes`` will be
+ created.
+
+ Returns:
+ The ``matplotlib`` PathCollection object created by scatter.
+
+ .. _Matplotlib Documentation: https://matplotlib.org/stable/gallery/color/named_colors.html
+ """
+ x_values, y_values = pareto_front(diversification_values, roc_values)
+ return plot_points(
+ x_values,
+ y_values,
+ color=color,
+ c=c,
+ alpha=alpha,
+ vmin=vmin,
+ vmax=vmax,
+ cmap=cmap,
+ label=label,
+ ax=ax,
+ )
+