Switched to sparray and np.floating types

robustocs/loaders.py

@@ -48,7 +48,7 @@ def count_sparse_nnz(filename: str) -> int:
 
 
 def load_symmetric_matrix(filename: str, dimension: int
-                          ) -> npt.NDArray[np.float64]:
+                          ) -> npt.NDArray[np.floating]:
     """
     Since NumPy doesn't have a stock way to load symmetric matrices stored in
     symmetric coordinate format, this adds one.
@@ -66,7 +66,7 @@ def load_symmetric_matrix(filename: str, dimension: int
         The matrix represented by the file.
     """
 
-    matrix = np.zeros([dimension, dimension], dtype=float)
+    matrix = np.zeros([dimension, dimension], dtype=np.floating)
 
     with open(filename, 'r') as file:
         for line in file:
@@ -79,7 +79,7 @@ def load_symmetric_matrix(filename: str, dimension: int
 
 
 def load_symmetric_matrix_coo(filename: str, dimension: int, nnz: int
-                              ) -> sparse.spmatrix:
+                              ) -> sparse.sparray:
     """
     Since neither NumPy or SciPy have a stock way to load symmetric matrices
     into sparse coordinate format, this adds one.
@@ -101,9 +101,9 @@ def load_symmetric_matrix_coo(filename: str, dimension: int, nnz: int
     """
 
     # preallocate storage arrays
-    rows: npt.NDArray[np.int8] = np.zeros(nnz)
-    cols: npt.NDArray[np.int8] = np.zeros(nnz)
-    vals: npt.NDArray[np.float64] = np.zeros(nnz)
+    rows: npt.NDArray[np.integer] = np.zeros(nnz, dtype=np.integer)
+    cols: npt.NDArray[np.integer] = np.zeros(nnz, dtype=np.integer)
+    vals: npt.NDArray[np.floating] = np.zeros(nnz, dtype=np.floating)
 
     with open(filename, 'r') as file:
         index: int = 0
@@ -126,15 +126,15 @@ def load_symmetric_matrix_coo(filename: str, dimension: int, nnz: int
 
             index += 1
 
-    return sparse.coo_matrix(
+    return sparse.coo_array(
         (vals, (rows, cols)),
         shape=(dimension, dimension),
-        dtype=np.float64
+        dtype=np.floating
     )
 
 
 def load_symmetric_matrix_csr(filename: str, dimension: int, nnz: int
-                              ) -> sparse.spmatrix:
+                              ) -> sparse.sparray:
     """
     Loads a symmetric matrix into compressed sparse row format. It does this
     by first loading into sparse coordinate format and then converting with
@@ -157,7 +157,7 @@ def load_symmetric_matrix_csr(filename: str, dimension: int, nnz: int
     """
 
     matrix = load_symmetric_matrix_coo(filename, dimension, nnz)
-    return sparse.csr_matrix(matrix)
+    return sparse.csr_array(matrix)
 
 
 def load_ped(filename: str) -> dict[int, list[int]]:
@@ -191,7 +191,7 @@ def load_ped(filename: str) -> dict[int, list[int]]:
 # MATRIX GENERATORS
 # Utility functions for generating matrices from pedigree data.
 
-def makeA(pedigree: dict[int, list[int]]) -> npt.NDArray[np.float64]:
+def makeA(pedigree: dict[int, list[int]]) -> npt.NDArray[np.floating]:
     """
     Constructs Wright's Numerator Relationship Matrix (WNRM) from a given
     pedigree structure.
@@ -210,7 +210,7 @@ def makeA(pedigree: dict[int, list[int]]) -> npt.NDArray[np.float64]:
 
     m = len(pedigree)
     # preallocate memory for A
-    A = np.zeros((m, m), dtype=float)
+    A = np.zeros((m, m), dtype=np.floating)
 
     # iterate over rows
     for i in range(0, m):
@@ -238,9 +238,9 @@ def load_problem(A_filename: str, E_filename: str, S_filename: str,
                  nnzA: int | None = None, nnzS: int | None = None,
                  dimension: int | None = None, pedigree: bool = False,
                  issparse: bool = False
-                 ) -> tuple[npt.NDArray[np.float64] | sparse.spmatrix,
-                            npt.NDArray[np.float64],
-                            npt.NDArray[np.float64] | sparse.spmatrix,
+                 ) -> tuple[npt.NDArray[np.floating] | sparse.sparray,
+                            npt.NDArray[np.floating],
+                            npt.NDArray[np.floating] | sparse.sparray,
                             int]:
     """
     Load a robust genetic selection problem into Python.
@@ -274,12 +274,12 @@ def load_problem(A_filename: str, E_filename: str, S_filename: str,
 
     Returns
     -------
-    ndarray or spmatrix
+    ndarray or sparray
         Covariance matrix of candidates in the cohort.
     ndarray
         Vector of expected values of the expected breeding values of
         candidates in the cohort.
-    ndarray or spmatrix
+    ndarray or sparray
         Covariance matrix of expected breeding values of candidates in the
         cohort.
     int
@@ -307,7 +307,7 @@ def load_problem(A_filename: str, E_filename: str, S_filename: str,
         A = makeA(load_ped(A_filename))
         # HACK this loads the full matrix, then converts it down to sparse
         if issparse:
-            A = sparse.coo_matrix(A)
+            A = sparse.coo_array(A)
     else:
         if issparse:
             if not nnzA:

robustocs/solvers.py

@@ -23,19 +23,19 @@
 
 
 def gurobi_standard_genetics(
-    sigma: npt.NDArray[np.float64] | sparse.spmatrix,
-    mu: npt.NDArray[np.float64],
+    sigma: npt.NDArray[np.floating] | sparse.sparray,
+    mu: npt.NDArray[np.floating],
     sires,  # type could be np.ndarray, sets[ints], lists[int], range, etc
     dams,   # type could be np.ndarray, sets[ints], lists[int], range, etc
     lam: float,  # cannot be called `lambda`, that's reserved in Python
     dimension: int,
-    upper_bound: npt.NDArray[np.float64] | float = 1.0,
-    lower_bound: npt.NDArray[np.float64] | float = 0.0,
+    upper_bound: npt.NDArray[np.floating] | float = 1.0,
+    lower_bound: npt.NDArray[np.floating] | float = 0.0,
     time_limit: float | None = None,
     max_duality_gap: float | None = None,
     model_output: str = '',
     debug: bool = False
-) -> tuple[npt.NDArray[np.float64], float]:
+) -> tuple[npt.NDArray[np.floating], float]:
     """
     Solve the standard genetic selection problem using Gurobi.
 
@@ -118,7 +118,7 @@ def gurobi_standard_genetics(
     )
 
     # set up the two sum-to-half constraints
-    M = np.zeros((2, dimension), dtype=int)
+    M = np.zeros((2, dimension), dtype=np.bool)
     # define the M so that column i is [1;0] if i is a sire and [0;1] otherwise
     M[0, sires] = 1
     M[1, dams] = 1
@@ -141,21 +141,21 @@ def gurobi_standard_genetics(
 
 
 def gurobi_robust_genetics(
-    sigma: npt.NDArray[np.float64] | sparse.spmatrix,
-    mubar: npt.NDArray[np.float64],
-    omega: npt.NDArray[np.float64] | sparse.spmatrix,
+    sigma: npt.NDArray[np.floating] | sparse.sparray,
+    mubar: npt.NDArray[np.floating],
+    omega: npt.NDArray[np.floating] | sparse.sparray,
     sires,  # type could be np.ndarray, sets[ints], lists[int], range, etc
     dams,   # type could be np.ndarray, sets[ints], lists[int], range, etc
     lam: float,  # cannot be called `lambda`, that's reserved in Python
     kappa: float,
     dimension: int,
-    upper_bound: npt.NDArray[np.float64] | float = 1.0,
-    lower_bound: npt.NDArray[np.float64] | float = 0.0,
+    upper_bound: npt.NDArray[np.floating] | float = 1.0,
+    lower_bound: npt.NDArray[np.floating] | float = 0.0,
     time_limit: float | None = None,
     max_duality_gap: float | None = None,
     model_output: str = '',
     debug: bool = False
-) -> tuple[npt.NDArray[np.float64], float, float]:
+) -> tuple[npt.NDArray[np.floating], float, float]:
     """
     Solve the robust genetic selection problem using Gurobi.
 
@@ -227,7 +227,7 @@ def gurobi_robust_genetics(
     ndarray
         Portfolio vector which Gurobi has determined is a solution.
     float
-        Auxillary variable corresponding to uncertainty associated with the
+        Auxiliary variable corresponding to uncertainty associated with the
         portfolio vector which Gurobi has determined is a solution.
     float
         Value of the objective function for returned solution vector.
@@ -254,7 +254,7 @@ def gurobi_robust_genetics(
     )
 
     # set up the two sum-to-half constraints
-    M = np.zeros((2, dimension), dtype=int)
+    M = np.zeros((2, dimension), dtype=np.bool)
     # define the M so that column i is [1;0] if i is a sire and [0;1] otherwise
     M[0, sires] = 1
     M[1, dams] = 1
@@ -280,23 +280,23 @@ def gurobi_robust_genetics(
 
 
 def gurobi_robust_genetics_sqp(
-    sigma: npt.NDArray[np.float64] | sparse.spmatrix,
-    mubar: npt.NDArray[np.float64],
-    omega: npt.NDArray[np.float64] | sparse.spmatrix,
+    sigma: npt.NDArray[np.floating] | sparse.sparray,
+    mubar: npt.NDArray[np.floating],
+    omega: npt.NDArray[np.floating] | sparse.sparray,
     sires,  # type could be np.ndarray, sets[ints], lists[int], range, etc
     dams,   # type could be np.ndarray, sets[ints], lists[int], range, etc
     lam: float,  # cannot be called `lambda`, that's reserved in Python
     kappa: float,
     dimension: int,
-    upper_bound: npt.NDArray[np.float64] | float = 1.0,
-    lower_bound: npt.NDArray[np.float64] | float = 0.0,
+    upper_bound: npt.NDArray[np.floating] | float = 1.0,
+    lower_bound: npt.NDArray[np.floating] | float = 0.0,
     time_limit: float | None = None,
     max_duality_gap: float | None = None,
     max_iterations: int = 1000,
     robust_gap_tol: float = 1e-7,
     model_output: str = '',
     debug: bool = False
-) -> tuple[npt.NDArray[np.float64], float, float]:
+) -> tuple[npt.NDArray[np.floating], float, float]:
     """
     Solve the robust genetic selection problem using SQP in Gurobi.
 
@@ -376,7 +376,7 @@ def gurobi_robust_genetics_sqp(
     ndarray
         Portfolio vector which Gurobi has determined is a solution.
     float
-        Auxillary variable corresponding to uncertainty associated with the
+        Auxiliary variable corresponding to uncertainty associated with the
         portfolio vector which Gurobi has determined is a solution.
     float
         Value of the objective function for returned solution vector.
@@ -403,12 +403,12 @@ def gurobi_robust_genetics_sqp(
     )
 
     # set up the two sum-to-half constraints
-    M = np.zeros((2, dimension), dtype=int)
+    M = np.zeros((2, dimension), dtype=np.bool)
     # define the M so that column i is [1;0] if i is a sire and [0;1] otherwise
     M[0, sires] = 1
     M[1, dams] = 1
     # define the right hand side of the constraint Mx = m
-    m = np.full(2, 0.5, dtype=float)
+    m = np.full(2, 0.5, dtype=np.floating)
     model.addConstr(M@w == m, name="sum-to-half")
 
     # optional controls to stop Gurobi taking too long
@@ -438,7 +438,7 @@ def gurobi_robust_genetics_sqp(
             print("No active constraints!")
 
         # z coefficient for the new constraint
-        w_star: npt.NDArray[np.float64] = np.array(w.X)
+        w_star: npt.NDArray[np.floating] = np.array(w.X)
         alpha: float = sqrt(w_star.transpose()@omega@w_star)
 
         # if gap between z and w'Omega w has converged, done
@@ -452,8 +452,8 @@ def gurobi_robust_genetics_sqp(
 
 
 def highs_bound_like(dimension: int,
-                     value: float | list[float] | npt.NDArray[np.float64]
-                     ):  # -> npt.NDArray[np.float64] | list[float] # BUG broke
+                     value: float | list[float] | npt.NDArray[np.floating]
+                     ):  # BUG -> npt.NDArray[np.floating] | list[float]
     """
     Helper function which allows HiGHS to interpret variable bounds specified
     either as a vector or a single floating point value. If `value` is an array
@@ -466,18 +466,18 @@ def highs_bound_like(dimension: int,
 
 def highs_standard_genetics(
     sigma: sparse.spmatrix,
-    mu: npt.NDArray[np.float64],
+    mu: npt.NDArray[np.floating],
     sires,  # type could be np.ndarray, sets[ints], lists[int], range, etc
     dams,   # type could be np.ndarray, sets[ints], lists[int], range, etc
     lam: float,  # cannot be called `lambda`, that's reserved in Python
     dimension: int,
-    upper_bound: npt.NDArray[np.float64] | list[float] | float = 1.0,
-    lower_bound: npt.NDArray[np.float64] | list[float] | float = 0.0,
+    upper_bound: npt.NDArray[np.floating] | list[float] | float = 1.0,
+    lower_bound: npt.NDArray[np.floating] | list[float] | float = 0.0,
     time_limit: float | None = None,
     max_duality_gap: float | None = None,
     model_output: str = '',
     debug: bool = False
-) -> tuple[npt.NDArray[np.float64], float]:
+) -> tuple[npt.NDArray[np.floating], float]:
     """
     Solve the standard genetic selection problem using HiGHS.
 
@@ -609,7 +609,7 @@ def highs_standard_genetics(
         raise RuntimeError
 
     # by default, col_value is a stock-Python list
-    solution: npt.NDArray[np.float64] = np.array(h.getSolution().col_value)
+    solution: npt.NDArray[np.floating] = np.array(h.getSolution().col_value)
     # we negated the objective function, so negate it back
     objective_value: float = -h.getInfo().objective_function_value
 
@@ -618,22 +618,22 @@ def highs_standard_genetics(
 
 def highs_robust_genetics_sqp(
     sigma: sparse.spmatrix,
-    mubar: npt.NDArray[np.float64],
-    omega: npt.NDArray[np.float64] | sparse.spmatrix,
+    mubar: npt.NDArray[np.floating],
+    omega: npt.NDArray[np.floating] | sparse.spmatrix,
     sires,  # type could be np.ndarray, sets[ints], lists[int], range, etc
     dams,   # type could be np.ndarray, sets[ints], lists[int], range, etc
     lam: float,  # cannot be called `lambda`, that's reserved in Python
     kappa: float,
     dimension: int,
-    upper_bound: npt.NDArray[np.float64] | list[float] | float = 1.0,
-    lower_bound: npt.NDArray[np.float64] | list[float] | float = 0.0,
+    upper_bound: npt.NDArray[np.floating] | list[float] | float = 1.0,
+    lower_bound: npt.NDArray[np.floating] | list[float] | float = 0.0,
     time_limit: float | None = None,
     max_duality_gap: float | None = None,
     max_iterations: int = 1000,
     robust_gap_tol: float = 1e-7,
     model_output: str = '',
     debug: bool = False
-) -> tuple[npt.NDArray[np.float64], float, float]:
+) -> tuple[npt.NDArray[np.floating], float, float]:
     """
     Solve the robust genetic selection problem using SQP in HiGHS.
 
@@ -795,7 +795,7 @@ def highs_robust_genetics_sqp(
 
         # by default, col_value is a stock-Python list
         solution: list[float] = h.getSolution().col_value
-        w_star: npt.NDArray[np.float64] = np.array(solution[:-1])
+        w_star: npt.NDArray[np.floating] = np.array(solution[:-1])
         z_star: float = solution[-1]
 
         # we negated the objective function, so negate it back
@@ -825,7 +825,7 @@ def highs_robust_genetics_sqp(
         # add a new plane to the approximation of the uncertainty cone
         num_nz: int = dimension + 1  # HACK assuming entirely dense
         index: range = range(dimension + 1)
-        value: npt.NDArray[np.float64] = np.append(-omega@w_star, alpha)
+        value: npt.NDArray[np.floating] = np.append(-omega@w_star, alpha)
         h.addRow(0, inf, num_nz, index, value)
 
     # final value of solution is the z value, return separately

robustocs/utils.py

@@ -15,8 +15,8 @@
 
 
 def print_compare_solutions(
-    portfolio1: npt.NDArray[np.float64],
-    portfolio2: npt.NDArray[np.float64],
+    portfolio1: npt.NDArray[np.floating],
+    portfolio2: npt.NDArray[np.floating],
     objective1: float,
     objective2: float,
     precision: int = 5,
@@ -106,8 +106,8 @@ def obj_string(name: str, value: float, precision: int,
 
 def check_uncertainty_constraint(
     z: float,
-    w: npt.NDArray[np.float64],
-    omega: npt.NDArray[np.float64] | sparse.spmatrix,
+    w: npt.NDArray[np.floating],
+    omega: npt.NDArray[np.floating] | sparse.sparray,
     tol: float = 1e-7,
     debug: bool = False
 ) -> bool:
@@ -123,7 +123,7 @@ def check_uncertainty_constraint(
     Parameters
     ----------
     z : float
-        Auxillary variable from a solution to the robust selection problem.
+        Auxiliary variable from a solution to the robust selection problem.
     w : ndarray
         Portfolio vector from a solution to the robust selection problem.
     omega : ndarray