Numpy 2 compat, scipy sparse arrays, and python version floor bump

README.md

@@ -12,7 +12,7 @@
 Adelie is a fast and flexible Python package for solving 
 lasso, elastic net, group lasso, and group elastic net problems. 
 
-- **Installation**: [https://jamesyang007.github.io/adelie/notebooks/installation.html](https://jamesyang007.github.io/adelie/notebooks/installation.html)
+- **Installation**: [https://jamesyang007.github.io/adelie/user_guide/notebooks/installation.html](https://jamesyang007.github.io/adelie/user_guide/notebooks/installation.html)
 - **Documentation**: [https://jamesyang007.github.io/adelie](https://jamesyang007.github.io/adelie/)
 - **Source code**: [https://github.com/JamesYang007/adelie](https://github.com/JamesYang007/adelie)
 - **Issue Tracker**: [https://github.com/JamesYang007/adelie/issues](https://github.com/JamesYang007/adelie/issues)

adelie/__init__.py

@@ -1,4 +1,4 @@
-__version__ = "1.1.53.dev"
+__version__ = "2.0.0.dev"
 
 # Set environment flags before loading adelie_core
 import os

adelie/constraint.py

@@ -10,7 +10,7 @@
     MatrixConstraintBase32,
     MatrixConstraintBase64,
 )
-from scipy.sparse import csr_matrix
+from scipy.sparse import csr_array, csr_matrix
 from typing import Union
 import numpy as np
 
@@ -135,7 +135,7 @@ def __init__(self):
 
 
 def linear(
-    A: Union[np.ndarray, csr_matrix, MatrixConstraintBase32, MatrixConstraintBase64],
+    A: Union[np.ndarray, csr_array, csr_matrix, MatrixConstraintBase32, MatrixConstraintBase64],
     lower: np.ndarray,
     upper: np.ndarray,
     *,
@@ -155,15 +155,15 @@ def linear(
 
     Parameters
     ----------
-    A : (m, d) Union[ndarray, csr_matrix, MatrixConstraintBase32, MatrixConstraintBase64]
+    A : (m, d) Union[ndarray, csr_array, csr_matrix, MatrixConstraintBase32, MatrixConstraintBase64]
         Constraint matrix :math:`A`.
     lower : (m,) ndarray
         Lower bound :math:`\\ell`.
     upper : (m,) ndarray
         Upper bound :math:`u`.
     vars : ndarray, optional
         Equivalent to :math:`\\mathrm{diag}(AA^\\top)`.
-        If ``None`` and ``A`` is ``ndarray`` or ``csr_matrix``, it is computed internally.
+        If ``None`` and ``A`` is ``ndarray`` or sparse, it is computed internally.
         Otherwise, it must be explicitly provided by the user.
         Default is ``None``.
     copy : bool, optional
@@ -241,7 +241,7 @@ def linear(
 
         A, _ = _coerce_dtype(A, dtype)
         A = matrix.dense(A, method="constraint", copy=copy)
-    elif isinstance(A, csr_matrix):
+    elif isinstance(A, (csr_array, csr_matrix)):
         if vars is None:
             vars = (A ** 2).sum(axis=1)
 
@@ -293,7 +293,7 @@ def __init__(self):
             self._A = A
             self._lower = np.array(lower, dtype=dtype)
             self._upper = np.array(upper, dtype=dtype)
-            self._vars = np.array(vars, copy=copy, dtype=dtype)
+            self._vars = np.array(vars, copy=copy or None, dtype=dtype)
             core_base.__init__(
                 self,
                 A=self._A,

adelie/diagnostic.py

@@ -20,7 +20,7 @@
 )
 from IPython.display import HTML
 from itertools import cycle
-from scipy.sparse import csr_matrix
+from scipy.sparse import csr_array, csr_matrix
 from typing import Union
 import numpy as np
 import matplotlib.pyplot as plt
@@ -29,7 +29,7 @@
 
 def predict(
     X: Union[np.ndarray, MatrixNaiveBase32, MatrixNaiveBase64],
-    betas: Union[np.ndarray, csr_matrix],
+    betas: Union[np.ndarray, csr_array, csr_matrix],
     intercepts: np.ndarray,
     *,
     offsets: np.ndarray =None,
@@ -63,7 +63,7 @@ def predict(
     X : (n, p) Union[ndarray, MatrixNaiveBase32, MatrixNaiveBase64]
         Feature matrix.
         It is typically one of the matrices defined in :mod:`adelie.matrix` submodule or :class:`numpy.ndarray`.
-    betas : (L, p) or (L, p*K) Union[ndarray, csr_matrix]
+    betas : (L, p) or (L, p*K) Union[ndarray, csr_array, csr_matrix]
         Coefficient vectors :math:`\\beta`.
     intercepts : (L,) or (L, K) ndarray
         Intercepts :math:`\\beta_0`.
@@ -112,10 +112,10 @@ def predict(
     if isinstance(betas, np.ndarray):
         for i in range(etas.shape[0]):
             X.btmul(0, X.cols(), betas[i], etas[i].ravel())
-    elif isinstance(betas, csr_matrix):
+    elif isinstance(betas, (csr_array, csr_matrix)):
         X.sp_tmul(betas, etas.reshape((L, -1))) 
     else:
-        raise RuntimeError("beta is not one of np.ndarray or scipy.sparse.csr_matrix.")
+        raise RuntimeError("beta is not one of np.ndarray or scipy.sparse.csr_array/csr_matrix.")
     etas += intercepts[:, None] + offsets
 
     return etas
@@ -124,7 +124,7 @@ def predict(
 def objective(
     X: Union[np.ndarray, MatrixNaiveBase32, MatrixNaiveBase64], 
     glm: Union[GlmBase32, GlmBase64, GlmMultiBase32, GlmMultiBase64],
-    betas: Union[np.ndarray, csr_matrix], 
+    betas: Union[np.ndarray, csr_array, csr_matrix], 
     intercepts: np.ndarray,
     lmdas: np.ndarray, 
     *,
@@ -148,7 +148,7 @@ def objective(
     glm : Union[GlmBase32, GlmBase64, GlmMultiBase32, GlmMultiBase64] 
         GLM object.
         It is typically one of the GLM classes defined in :mod:`adelie.glm` submodule.
-    betas : (L, p) or (L, p*K) Union[ndarray, csr_matrix]
+    betas : (L, p) or (L, p*K) Union[ndarray, csr_array, csr_matrix]
         Coefficient vectors :math:`\\beta`.
     intercepts : (L,) or (L, K) ndarray
         Intercepts :math:`\\beta_0`.
@@ -259,7 +259,7 @@ def objective(
                 np.float32: core.solver.compute_penalty_dense_32,
                 np.float64: core.solver.compute_penalty_dense_64,
             }[dtype]
-        elif isinstance(betas, csr_matrix): 
+        elif isinstance(betas, (csr_array, csr_matrix)):
             penalty_f = {
                 np.float32: core.solver.compute_penalty_sparse_32,
                 np.float64: core.solver.compute_penalty_sparse_64,
@@ -388,8 +388,8 @@ def gradients(
 
 def gradient_norms(
     grads: np.ndarray,
-    betas: csr_matrix,
-    duals: csr_matrix,
+    betas: Union[csr_array, csr_matrix],
+    duals: Union[csr_array, csr_matrix],
     lmdas: np.ndarray,
     *, 
     constraints: list[Union[ConstraintBase32, ConstraintBase64]] =None,
@@ -423,9 +423,9 @@ def gradient_norms(
     ----------
     grads : (L, p) or (L, p, K) ndarray
         Gradients.
-    betas : (L, p) or (L, p*K) csr_matrix
+    betas : (L, p) or (L, p*K) Union[csr_array, csr_matrix]
         Coefficient vectors :math:`\\beta`.
-    duals : (L, d) csr_matrix
+    duals : (L, d) Union[csr_array, csr_matrix]
         Dual vectors :math:`\\mu`.
     lmdas : (L,) ndarray
         Regularization parameters :math:`\\lambda`.
@@ -497,8 +497,8 @@ def gradient_norms(
             dual_groups = render_dual_groups(constraints)
             mu_grads = np.zeros(grads.shape, dtype=dtype)
             for k in range(L):
-                beta_curr = betas[k].toarray()[0].astype(dtype)
-                mu_curr = duals[k].toarray()[0].astype(dtype)
+                beta_curr = betas[k].toarray().ravel().astype(dtype)
+                mu_curr = duals[k].toarray().ravel().astype(dtype)
                 mu_grads_curr = mu_grads[k].astype(dtype)
                 for constraint, g, gs, dg in zip(
                     constraints,
@@ -576,7 +576,7 @@ def gradient_scores(
 
 def coefficient(
     lmda: float,
-    betas: csr_matrix,
+    betas: Union[csr_array, csr_matrix],
     intercepts: np.ndarray,
     lmdas: np.ndarray,
 ):
@@ -604,7 +604,7 @@ def coefficient(
     ----------
     lmda : float
         New regularization parameter at which to find the solution.
-    betas : (L, p) csr_matrix
+    betas : (L, p) Union[csr_array, csr_matrix]
         Coefficient vectors :math:`\\beta`.
     intercepts : (L,) ndarray
         Intercepts.
@@ -613,15 +613,15 @@ def coefficient(
 
     Returns
     -------
-    beta : (1, p) csr_matrix
+    beta : (1, p) csr_array
         Linearly interpolated coefficient vector at :math:`\\lambda`.
     intercept : float
         Linearly interpolated intercept at :math:`\\lambda`.
     """
     if len(lmdas) == 0:
         raise RuntimeError("lmdas must be non-empty!")
     if len(lmdas) == 1:
-        return betas, lmdas
+        return betas, intercepts[0]
     order = np.argsort(lmdas)
     idx = np.searchsorted(
         lmdas,
@@ -635,9 +635,9 @@ def coefficient(
             "Returning boundary solution."
         )
         idx = np.clip(idx, 0, lmdas.shape[0]-1)
-        return betas[idx], intercepts[idx]
+        return betas[idx:idx+1], intercepts[idx]
 
-    left, right = betas[idx-1], betas[idx]
+    left, right = betas[idx-1:idx], betas[idx:idx+1]
     weight = (lmda - lmdas[idx]) / (lmdas[idx-1] - lmdas[idx])
     beta = left.multiply(weight) + right.multiply(1-weight)
     left, right = intercepts[idx-1], intercepts[idx]
@@ -647,7 +647,7 @@ def coefficient(
 
 
 def plot_coefficients(
-    betas: csr_matrix,
+    betas: Union[csr_array, csr_matrix],
     lmdas: np.ndarray,
     groups: np.ndarray,
     group_sizes: np.ndarray,
@@ -658,7 +658,7 @@ def plot_coefficients(
 
     Parameters
     ----------
-    betas : (L, p) csr_matrix
+    betas : (L, p) Union[csr_array, csr_matrix]
         Coefficient vectors :math:`\\beta`.
     lmdas : (L,) ndarray
         Regularization parameters :math:`\\lambda`.

adelie/glm.py

@@ -15,7 +15,7 @@ def _coerce_dtype(y, dtype):
         np.dtype("float64"): np.float64,
     }
     valid_dtypes = list(dtype_map.keys())
-    y = np.array(y, copy=False, order="C")  # important!
+    y = np.asarray(y, order="C")  # important!
     if dtype is None:
         if not (y.dtype in valid_dtypes):
             raise RuntimeError(