Refactor weekday calc into utils

_delphi_utils_python/.pylintrc

@@ -12,9 +12,9 @@ disable=logging-format-interpolation,
 [BASIC]
 
 # Allow arbitrarily short-named variables.
-variable-rgx=[a-z_][a-z0-9_]*
-argument-rgx=[a-z_][a-z0-9_]*
-attr-rgx=[a-z_][a-z0-9_]*
+variable-rgx=([a-z_][a-z0-9_]*|[a-zA-Z])
+argument-rgx=([a-z_][a-z0-9_]*|[a-zA-Z])
+attr-rgx=([a-z_][a-z0-9_]*|[a-zA-Z])
 
 [DESIGN]
 

_delphi_utils_python/delphi_utils/__init__.py

@@ -13,5 +13,6 @@
 from .smooth import Smoother
 from .signal import add_prefix
 from .nancodes import Nans
+from .weekday import Weekday
 
 __version__ = "0.2.1"

doctor_visits/delphi_doctor_visits/weekday.py → _delphi_utils_python/delphi_utils/weekday.py

@@ -3,23 +3,48 @@
 
 Created: 2020-05-06
 """
-
-
-
-# third party
 import cvxpy as cp
-from cvxpy.error import SolverError
 import numpy as np
-
-# first party
-from .config import Config
+from cvxpy.error import SolverError
 
 
 class Weekday:
     """Class to handle weekday effects."""
 
     @staticmethod
-    def get_params(data, logger):
+    def get_params(data, denominator_col, numerator_cols, date_col, scales, logger):
+        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.
+        """
+        tmp = data.reset_index()
+        denoms = tmp.groupby(date_col).sum()[denominator_col]
+        nums = tmp.groupby(date_col).sum()[numerator_cols]
+
+        # Construct design matrix to have weekday indicator columns and then day
+        # indicators.
+        X = np.zeros((nums.shape[0], 6 + nums.shape[0]))
+        not_sunday = np.where(nums.index.dayofweek != 6)[0]
+        X[not_sunday, np.array(nums.index.dayofweek)[not_sunday]] = 1
+        X[np.where(nums.index.dayofweek == 6)[0], :6] = -1
+        X[:, 6:] = np.eye(X.shape[0])
+
+        npnums, npdenoms = np.array(nums), np.array(denoms)
+        params = np.zeros((nums.shape[1], X.shape[1]))
+
+        # 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)
+            if result is None:
+                logger.error("Unable to calculate weekday correction")
+            else:
+                params[i,:] = result
+
+        return params
+
+    @staticmethod
+    def _fit(X, scales, npnums, npdenoms):
         r"""Correct a signal estimated as numerator/denominator for weekday effects.
 
         The ordinary estimate would be numerator_t/denominator_t for each time point
@@ -53,57 +78,31 @@ def get_params(data, logger):
 
         ll = (numerator * (X*b + log(denominator)) - sum(exp(X*b) + log(denominator)))
                 / num_days
-
-        Return a matrix of parameters: the entire vector of betas, for each time
-        series column in the data.
         """
-        denoms = data.groupby(Config.DATE_COL).sum()["Denominator"]
-        nums = data.groupby(Config.DATE_COL).sum()[Config.CLI_COLS + Config.FLU1_COL]
-
-        # Construct design matrix to have weekday indicator columns and then day
-        # indicators.
-        X = np.zeros((nums.shape[0], 6 + nums.shape[0]))
-        not_sunday = np.where(nums.index.dayofweek != 6)[0]
-        X[not_sunday, np.array(nums.index.dayofweek)[not_sunday]] = 1
-        X[np.where(nums.index.dayofweek == 6)[0], :6] = -1
-        X[:, 6:] = np.eye(X.shape[0])
-
-        npnums, npdenoms = np.array(nums), np.array(denoms)
-        params = np.zeros((nums.shape[1], X.shape[1]))
-
-        # Loop over the available numerator columns and smooth each separately.
-        for i in range(nums.shape[1]):
-            b = cp.Variable((X.shape[1]))
-
-            lmbda = cp.Parameter(nonneg=True)
-            lmbda.value = 10  # Hard-coded for now, seems robust to changes
-
-            ll = (
-                cp.matmul(npnums[:, i], cp.matmul(X, b) + np.log(npdenoms))
-                - cp.sum(cp.exp(cp.matmul(X, b) + np.log(npdenoms)))
-            ) / X.shape[0]
-            penalty = (
-                lmbda * cp.norm(cp.diff(b[6:], 3), 1) / (X.shape[0] - 2)
-            )  # L-1 Norm of third differences, rewards smoothness
-            scales = [1, 1e5, 1e10, 1e15]
-            for scale in scales:
-                try:
-                    prob = cp.Problem(cp.Minimize((-ll + lmbda * penalty) / scale))
-                    _ = prob.solve()
-                    params[i,:] = b.value
-                    break
-                except SolverError:
-                    # If the magnitude of the objective function is too large, an error is
-                    # thrown; Rescale the objective function by going through loop
-                    pass
-            else:
-                # Leaving params[i,:] = 0 is equivalent to not performing weekday correction
-                logger.error("Unable to calculate weekday correction")
-
-        return params
+        b = cp.Variable((X.shape[1]))
+
+        lmbda = cp.Parameter(nonneg=True)
+        lmbda.value = 10  # Hard-coded for now, seems robust to changes
+
+        ll = (
+            cp.matmul(npnums, cp.matmul(X, b) + np.log(npdenoms)) -
+            cp.sum(cp.exp(cp.matmul(X, b) + np.log(npdenoms)))
+        ) / X.shape[0]
+        # L-1 Norm of third differences, rewards smoothness
+        penalty = lmbda * cp.norm(cp.diff(b[6:], 3), 1) / (X.shape[0] - 2)
+        for scale in scales:
+            try:
+                prob = cp.Problem(cp.Minimize((-ll + lmbda * penalty) / scale))
+                _ = prob.solve()
+                return b.value
+            except SolverError:
+                # If the magnitude of the objective function is too large, an error is
+                # thrown; Rescale the objective function by going through loop
+                continue
+        return None
 
     @staticmethod
-    def calc_adjustment(params, sub_data):
+    def calc_adjustment(params, sub_data, cols, date_col):
         """Apply the weekday adjustment to a specific time series.
 
         Extracts the weekday fixed effects from the parameters and uses these to
@@ -122,14 +121,15 @@ def calc_adjustment(params, sub_data):
         -- this has the same effect.
 
         """
-        for i, c in enumerate(Config.CLI_COLS + Config.FLU1_COL):
-            wd_correction = np.zeros((len(sub_data[c])))
+        tmp = sub_data.copy()
+
+        for i, c in enumerate(cols):
+            wd_correction = np.zeros((len(tmp[c])))
 
             for wd in range(7):
-                mask = sub_data[Config.DATE_COL].dt.dayofweek == wd
-                wd_correction[mask] = sub_data.loc[mask, c] / (
+                mask = tmp[date_col].dt.dayofweek == wd
+                wd_correction[mask] = tmp.loc[mask, c] / (
                     np.exp(params[i, wd]) if wd < 6 else np.exp(-np.sum(params[i, :6]))
                 )
-            sub_data.loc[:, c] = wd_correction
-
-        return sub_data
+            tmp.loc[:, c] = wd_correction
+        return tmp

_delphi_utils_python/setup.py

@@ -7,6 +7,7 @@
 required = [
     "boto3",
     "covidcast",
+    "cvxpy",
     "epiweeks",
     "freezegun",
     "gitpython",

_delphi_utils_python/tests/test_weekday.py

@@ -0,0 +1,74 @@
+import logging
+
+import numpy as np
+import pandas as pd
+from delphi_utils.weekday import Weekday
+
+
+class TestWeekday:
+
+    TEST_DATA = pd.DataFrame({
+        "num": np.arange(1, 11, 1),
+        "den": np.arange(11, 21, 1),
+        "date": pd.date_range("2020-01-01", "2020-01-10")
+    })
+
+    def test_get_params(self):
+        TEST_LOGGER = logging.getLogger()
+
+        result = Weekday.get_params(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]])
+
+        result = Weekday.calc_adjustment(params, self.TEST_DATA, ["num"], "date")
+
+        # Data should be unchanged when params are 0's
+        assert np.allclose(result["num"].values, self.TEST_DATA["num"].values)
+        assert np.allclose(result["den"].values, self.TEST_DATA["den"].values)
+        assert np.array_equal(result["date"].values, self.TEST_DATA["date"].values)
+
+    def test_calc_adjustment(self):
+        params = np.array([[1, -1, 1, -1, 1, -1, 1]])
+
+        result = Weekday.calc_adjustment(params, self.TEST_DATA, ["num"], "date")
+
+        print(result["num"].values)
+        print(result["den"].values)
+        expected_nums = [
+            0.36787944,
+            5.43656366,
+            1.10363832,
+            10.87312731,
+            5,
+            2.20727665,
+            19.0279728,
+            2.94303553,
+            24.46453646,
+            3.67879441,
+        ]
+
+        # The date and "den" column are unchanged by this function
+        assert np.allclose(result["num"].values, expected_nums)
+        assert np.allclose(result["den"].values, self.TEST_DATA["den"].values)
+        assert np.array_equal(result["date"].values, self.TEST_DATA["date"].values)

changehc/delphi_changehc/__init__.py

@@ -13,6 +13,5 @@
 from . import run
 from . import sensor
 from . import update_sensor
-from . import weekday
 
 __version__ = "0.0.0"

changehc/delphi_changehc/run.py

@@ -184,7 +184,7 @@ def run_module(params: Dict[str, Dict[str, Any]]):
                              file_dict["flu_like"],file_dict["covid_like"],dropdate_dt,"fips")
                 more_stats = su_inst.update_sensor(
                     data,
-                    params["common"]["export_dir"]
+                    params["common"]["export_dir"],
                 )
                 stats.extend(more_stats)
 

changehc/delphi_changehc/update_sensor.py

@@ -11,13 +11,12 @@
 # third party
 import numpy as np
 import pandas as pd
-from delphi_utils import GeoMapper, add_prefix, create_export_csv
+from delphi_utils import GeoMapper, add_prefix, create_export_csv, Weekday
 
 # first party
 from .config import Config
 from .constants import SMOOTHED, SMOOTHED_ADJ, SMOOTHED_CLI, SMOOTHED_ADJ_CLI, NA
 from .sensor import CHCSensor
-from .weekday import Weekday
 
 
 def write_to_csv(df, geo_level, write_se, day_shift, out_name, logger, output_path=".", start_date=None, end_date=None):
@@ -176,8 +175,8 @@ def geo_reindex(self, data):
 
 
     def update_sensor(self,
-            data,
-            output_path):
+        data,
+        output_path):
         """Generate sensor values, and write to csv format.
 
         Args:
@@ -192,14 +191,22 @@ def update_sensor(self,
         data.reset_index(inplace=True)
         data_frame = self.geo_reindex(data)
         # handle if we need to adjust by weekday
-        wd_params = Weekday.get_params(data_frame) if self.weekday else None
+        wd_params = Weekday.get_params(
+            data_frame,
+            "den",
+            ["num"],
+            Config.DATE_COL,
+            [1, 1e5],
+            self.logger,
+        ) if self.weekday else None
         # run sensor fitting code (maybe in parallel)
         if not self.parallel:
             dfs = []
             for geo_id, sub_data in data_frame.groupby(level=0):
-                sub_data.reset_index(level=0,inplace=True)
+                sub_data.reset_index(inplace=True)
                 if self.weekday:
-                    sub_data = Weekday.calc_adjustment(wd_params, sub_data)
+                    sub_data = Weekday.calc_adjustment(wd_params, sub_data, ["num"], Config.DATE_COL)
+                sub_data.set_index(Config.DATE_COL, inplace=True)
                 res = CHCSensor.fit(sub_data, self.burnindate, geo_id, self.logger)
                 res = pd.DataFrame(res).loc[final_sensor_idxs]
                 dfs.append(res)
@@ -209,9 +216,10 @@ def update_sensor(self,
             with Pool(n_cpu) as pool:
                 pool_results = []
                 for geo_id, sub_data in data_frame.groupby(level=0,as_index=False):
-                    sub_data.reset_index(level=0, inplace=True)
+                    sub_data.reset_index(inplace=True)
                     if self.weekday:
-                        sub_data = Weekday.calc_adjustment(wd_params, sub_data)
+                        sub_data = Weekday.calc_adjustment(wd_params, sub_data, ["num"], Config.DATE_COL)
+                    sub_data.set_index(Config.DATE_COL, inplace=True)
                     pool_results.append(
                         pool.apply_async(
                             CHCSensor.fit, args=(sub_data, self.burnindate, geo_id, self.logger),