Add params to control solver, default Clarabel, pin cvxpy version

_delphi_utils_python/delphi_utils/weekday.py

@@ -12,12 +12,19 @@ class Weekday:
     """Class to handle weekday effects."""
 
     @staticmethod
-    def get_params(data, denominator_col, numerator_cols, date_col, scales, logger):
+    def get_params(data, denominator_col, numerator_cols, date_col, scales, logger, solver_override=None):
         r"""Fit weekday correction for each col in numerator_cols.
 
         Return a matrix of parameters: the entire vector of betas, for each time
         series column in the data.
+
+        solver: Historically used "ECOS" but due to numerical stability issues, "CLARABEL"
+        (introduced in cvxpy 1.3)is now the default solver in cvxpy 1.5.
         """
+        if solver_override is None:
+            solver = cp.CLARABEL
+        else:
+            solver = solver_override
         tmp = data.reset_index()
         denoms = tmp.groupby(date_col).sum()[denominator_col]
         nums = tmp.groupby(date_col).sum()[numerator_cols]
@@ -35,7 +42,7 @@ def get_params(data, denominator_col, numerator_cols, date_col, scales, logger):
 
         # Loop over the available numerator columns and smooth each separately.
         for i in range(nums.shape[1]):
-            result = Weekday._fit(X, scales, npnums[:, i], npdenoms)
+            result = Weekday._fit(X, scales, npnums[:, i], npdenoms, solver)
             if result is None:
                 logger.error("Unable to calculate weekday correction")
             else:
@@ -44,7 +51,18 @@ def get_params(data, denominator_col, numerator_cols, date_col, scales, logger):
         return params
 
     @staticmethod
-    def _fit(X, scales, npnums, npdenoms):
+    def get_params_legacy(data, denominator_col, numerator_cols, date_col, scales, logger):
+        r"""
+        Preserves older default behavior of using the ECOS solver.
+
+        NOTE: "ECOS" solver will not be installed by default as of cvxpy 1.6
+        """
+        return Weekday.get_params(
+            data, denominator_col, numerator_cols, date_col, scales, logger, solver_override=cp.ECOS
+        )
+
+    @staticmethod
+    def _fit(X, scales, npnums, npdenoms, solver):
         r"""Correct a signal estimated as numerator/denominator for weekday effects.
 
         The ordinary estimate would be numerator_t/denominator_t for each time point
@@ -78,6 +96,8 @@ def _fit(X, scales, npnums, npdenoms):
 
         ll = (numerator * (X*b + log(denominator)) - sum(exp(X*b) + log(denominator)))
                 / num_days
+
+        solver: Historically use "ECOS" but due to numerical issues, "CLARABEL" is now default.
         """
         b = cp.Variable((X.shape[1]))
 
@@ -93,7 +113,7 @@ def _fit(X, scales, npnums, npdenoms):
         for scale in scales:
             try:
                 prob = cp.Problem(cp.Minimize((-ll + lmbda * penalty) / scale))
-                _ = prob.solve(solver=cp.CLARABEL)
+                _ = prob.solve(solver=solver)
                 return b.value
             except SolverError:
                 # If the magnitude of the objective function is too large, an error is

_delphi_utils_python/tests/test_weekday.py

@@ -25,6 +25,31 @@ def test_get_params(self):
                                     -0.81464459, -0.76322013, -0.7667211,-0.8251475]])
         assert np.allclose(result, expected_result)
 
+    def test_get_params_legacy(self):
+        TEST_LOGGER = logging.getLogger()
+
+        result = Weekday.get_params_legacy(self.TEST_DATA, "den", ["num"], "date", [1], TEST_LOGGER)
+        print(result)
+        expected_result = [
+            -0.05993665,
+            -0.0727396,
+            -0.05618517,
+            0.0343405,
+            0.12534997,
+            0.04561813,
+            -2.27669028,
+            -1.89564374,
+            -1.5695407,
+            -1.29838116,
+            -1.08216513,
+            -0.92089259,
+            -0.81456355,
+            -0.76317802,
+            -0.76673598,
+            -0.82523745,
+        ]
+        assert np.allclose(result, expected_result)
+
     def test_calc_adjustment_with_zero_parameters(self):
         params = np.array([[0, 0, 0, 0, 0, 0, 0]])
 

claims_hosp/delphi_claims_hosp/update_indicator.py

@@ -13,13 +13,13 @@
 # third party
 import numpy as np
 import pandas as pd
-from delphi_utils import GeoMapper
 
 # first party
-from delphi_utils import Weekday
+from delphi_utils import GeoMapper, Weekday
+
 from .config import Config, GeoConstants
-from .load_data import load_data
 from .indicator import ClaimsHospIndicator
+from .load_data import load_data
 
 
 class ClaimsHospIndicatorUpdater:
@@ -152,15 +152,18 @@ def update_indicator(self, input_filepath, outpath, logger):
         data_frame = self.geo_reindex(data)
 
         # handle if we need to adjust by weekday
-        wd_params = Weekday.get_params(
-            data_frame,
-            "den",
-            ["num"],
-            Config.DATE_COL,
-            [1, 1e5],
-            logger,
-        ) if self.weekday else None
-
+        wd_params = (
+            Weekday.get_params_legacy(
+                data_frame,
+                "den",
+                ["num"],
+                Config.DATE_COL,
+                [1, 1e5],
+                logger,
+            )
+            if self.weekday
+            else None
+        )
         # run fitting code (maybe in parallel)
         rates = {}
         std_errs = {}

claims_hosp/setup.py

@@ -13,6 +13,7 @@
     "pylint==2.8.3",
     "pytest-cov",
     "pytest",
+    "cvxpy<1.6",
 ]
 
 setup(

doctor_visits/delphi_doctor_visits/update_sensor.py

@@ -18,6 +18,7 @@
 
 # first party
 from delphi_utils import Weekday
+
 from .config import Config
 from .geo_maps import GeoMaps
 from .sensor import DoctorVisitsSensor
@@ -125,15 +126,19 @@ def update_sensor(
         (burn_in_dates >= startdate) & (burn_in_dates <= enddate))[0][:len(sensor_dates)]
 
     # handle if we need to adjust by weekday
-    params = Weekday.get_params(
-        data,
-        "Denominator",
-        Config.CLI_COLS + Config.FLU1_COL,
-        Config.DATE_COL,
-        [1, 1e5, 1e10, 1e15],
-        logger,
-    ) if weekday else None
-    if weekday and np.any(np.all(params == 0,axis=1)):
+    params = (
+        Weekday.get_params(
+            data,
+            "Denominator",
+            Config.CLI_COLS + Config.FLU1_COL,
+            Config.DATE_COL,
+            [1, 1e5, 1e10, 1e15],
+            logger,
+        )
+        if weekday
+        else None
+    )
+    if weekday and np.any(np.all(params == 0, axis=1)):
         # Weekday correction failed for at least one count type
         return None
 

doctor_visits/setup.py

@@ -11,6 +11,7 @@
     "pytest-cov",
     "pytest",
     "scikit-learn",
+    "cvxpy>=1.5",
 ]
 
 setup(