fix: handle zero division in top down methods

.github/workflows/python-publish.yml

@@ -33,7 +33,7 @@ jobs:
       - name: Build package
         run: python -m build
       - name: Publish package
-        uses: pypa/gh-action-pypi-publish@67339c736fd9354cd4f8cb0b744f2b82a74b5c70 # v1.12.3
+        uses: pypa/gh-action-pypi-publish@76f52bc884231f62b9a034ebfe128415bbaabdfc # v1.12.4
         with:
           user: __token__
           password: ${{ secrets.PYPI_API_TOKEN }}

README.md

@@ -81,7 +81,6 @@ from hierarchicalforecast.evaluation import evaluate
 from hierarchicalforecast.methods import BottomUp, TopDown, MiddleOut
 from utilsforecast.losses import mse
 
-
 # Load TourismSmall dataset
 Y_df, S, tags = HierarchicalData.load('./data', 'TourismSmall')
 Y_df['ds'] = pd.to_datetime(Y_df['ds'])
@@ -92,10 +91,9 @@ Y_test_df  = Y_df.groupby('unique_id').tail(4)
 Y_train_df = Y_df.drop(Y_test_df.index)
 
 # Compute base auto-ARIMA predictions
-fcst = StatsForecast(df=Y_train_df,
-                     models=[AutoARIMA(season_length=4), Naive()],
-                     freq='Q', n_jobs=-1)
-Y_hat_df = fcst.forecast(h=4)
+fcst = StatsForecast(models=[AutoARIMA(season_length=4), Naive()],
+                     freq='QE', n_jobs=-1)
+Y_hat_df = fcst.forecast(df=Y_train_df, h=4)
 
 # Reconcile the base predictions
 reconcilers = [
@@ -116,7 +114,6 @@ Assumes you have a test dataframe.
 df = Y_rec_df.merge(Y_test_df, on=['unique_id', 'ds'])
 evaluation = evaluate(df = df,
                       tags = tags,
-                      train_df = Y_train_df,
                       metrics = [mse],
                       benchmark = "Naive")
 ```

hierarchicalforecast/core.py

@@ -19,7 +19,7 @@
 import narwhals as nw
 import numpy as np
 
-# %% ../nbs/src/core.ipynb 7
+# %% ../nbs/src/core.ipynb 6
 def _build_fn_name(fn) -> str:
     fn_name = type(fn).__name__
     func_params = fn.__dict__
@@ -42,7 +42,7 @@ def _build_fn_name(fn) -> str:
         fn_name += "_" + "_".join(func_params)
     return fn_name
 
-# %% ../nbs/src/core.ipynb 11
+# %% ../nbs/src/core.ipynb 10
 def _reverse_engineer_sigmah(
     Y_hat_df: Frame,
     y_hat: np.ndarray,
@@ -87,7 +87,7 @@ def _reverse_engineer_sigmah(
 
     return sigmah
 
-# %% ../nbs/src/core.ipynb 12
+# %% ../nbs/src/core.ipynb 11
 class HierarchicalReconciliation:
     """Hierarchical Reconciliation Class.
 

hierarchicalforecast/evaluation.py

@@ -1,7 +1,7 @@
 # AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/src/evaluation.ipynb.
 
 # %% auto 0
-__all__ = ['rel_mse', 'msse', 'scaled_crps', 'energy_score', 'log_score', 'HierarchicalEvaluation']
+__all__ = ['evaluate']
 
 # %% ../nbs/src/evaluation.ipynb 3
 import narwhals as nw
@@ -14,14 +14,14 @@
 from scipy.stats import multivariate_normal
 from typing import Callable, Optional, Union
 
-# %% ../nbs/src/evaluation.ipynb 7
+# %% ../nbs/src/evaluation.ipynb 5
 def _loss_deprecation_notice(loss):
     warnings.warn(
         f"This loss function ({loss}) will be deprecated in future releases. Please use the `utilsforecast.losses` function instead.",
         FutureWarning,
     )
 
-# %% ../nbs/src/evaluation.ipynb 8
+# %% ../nbs/src/evaluation.ipynb 6
 def _metric_protections(
     y: np.ndarray, y_hat: np.ndarray, weights: Optional[np.ndarray]
 ) -> None:
@@ -130,7 +130,7 @@ def mqloss(
 
     return mqloss
 
-# %% ../nbs/src/evaluation.ipynb 10
+# %% ../nbs/src/evaluation.ipynb 7
 def rel_mse(y, y_hat, y_train, mask=None):
     """Relative Mean Squared Error
 
@@ -170,7 +170,7 @@ def rel_mse(y, y_hat, y_train, mask=None):
     loss = loss / (norm + eps)
     return loss
 
-# %% ../nbs/src/evaluation.ipynb 12
+# %% ../nbs/src/evaluation.ipynb 8
 def msse(y, y_hat, y_train, mask=None):
     """Mean Squared Scaled Error
 
@@ -211,7 +211,7 @@ def msse(y, y_hat, y_train, mask=None):
     loss = loss / (norm + eps)
     return loss
 
-# %% ../nbs/src/evaluation.ipynb 14
+# %% ../nbs/src/evaluation.ipynb 9
 def scaled_crps(y, y_hat, quantiles):
     """Scaled Continues Ranked Probability Score
 
@@ -255,7 +255,7 @@ def scaled_crps(y, y_hat, quantiles):
     loss = 2 * loss * np.sum(np.ones(y.shape)) / (norm + eps)
     return loss
 
-# %% ../nbs/src/evaluation.ipynb 16
+# %% ../nbs/src/evaluation.ipynb 10
 def energy_score(y, y_sample1, y_sample2, beta=2):
     """Energy Score
 
@@ -302,7 +302,7 @@ def energy_score(y, y_sample1, y_sample2, beta=2):
     score = np.mean(term2 - 0.5 * term1)
     return score
 
-# %% ../nbs/src/evaluation.ipynb 17
+# %% ../nbs/src/evaluation.ipynb 11
 def log_score(y, y_hat, cov, allow_singular=True):
     """Log Score.
 
@@ -348,7 +348,7 @@ def log_score(y, y_hat, cov, allow_singular=True):
     score = np.mean(scores)
     return score
 
-# %% ../nbs/src/evaluation.ipynb 19
+# %% ../nbs/src/evaluation.ipynb 12
 class HierarchicalEvaluation:
     """Hierarchical Evaluation Class.
 
@@ -493,7 +493,7 @@ def evaluate(
 
         return evaluation
 
-# %% ../nbs/src/evaluation.ipynb 20
+# %% ../nbs/src/evaluation.ipynb 13
 def evaluate(
     df: FrameT,
     metrics: list[Callable],

hierarchicalforecast/methods.py

@@ -337,7 +337,11 @@ def _reconcile_fcst_proportions(
             fcst_parent = reconciled[idx_parent]
             childs_sum = y_hat[idx_childs].sum()
             for idx_child in idx_childs:
-                reconciled[idx_child] = y_hat[idx_child] * fcst_parent / childs_sum
+                if np.abs(childs_sum) < 1e-8:
+                    n_children = len(idx_childs)
+                    reconciled[idx_child] = fcst_parent / n_children
+                else:
+                    reconciled[idx_child] = y_hat[idx_child] * fcst_parent / childs_sum
     return reconciled
 
 # %% ../nbs/src/methods.ipynb 29

hierarchicalforecast/utils.py

@@ -17,14 +17,14 @@
 from sklearn.preprocessing import OneHotEncoder
 from typing import Optional, Union, Sequence
 
-# %% ../nbs/src/utils.ipynb 6
+# %% ../nbs/src/utils.ipynb 5
 # Global variables
 NUMBA_NOGIL = True
 NUMBA_CACHE = True
 NUMBA_PARALLEL = True
 NUMBA_FASTMATH = True
 
-# %% ../nbs/src/utils.ipynb 7
+# %% ../nbs/src/utils.ipynb 6
 class CodeTimer:
     def __init__(self, name=None, verbose=True):
         self.name = " '" + name + "'" if name else ""
@@ -43,7 +43,7 @@ def __exit__(self, exc_type, exc_value, traceback):
                 + " seconds"
             )
 
-# %% ../nbs/src/utils.ipynb 8
+# %% ../nbs/src/utils.ipynb 7
 def is_strictly_hierarchical(S: np.ndarray, tags: dict[str, np.ndarray]) -> bool:
     # main idea:
     # if S represents a strictly hierarchical structure
@@ -60,7 +60,7 @@ def is_strictly_hierarchical(S: np.ndarray, tags: dict[str, np.ndarray]) -> bool
     nodes = levels_.popitem()[1].size
     return paths == nodes
 
-# %% ../nbs/src/utils.ipynb 10
+# %% ../nbs/src/utils.ipynb 9
 def _to_upper_hierarchy(
     bottom_split: list[str], bottom_values: str, upper_key: str
 ) -> list[str]:
@@ -73,7 +73,7 @@ def join_upper(bottom_value):
 
     return [join_upper(val) for val in bottom_values]
 
-# %% ../nbs/src/utils.ipynb 11
+# %% ../nbs/src/utils.ipynb 10
 def aggregate(
     df: Frame,
     spec: list[list[str]],
@@ -238,7 +238,7 @@ def aggregate(
 
     return Y_df, S_df, tags
 
-# %% ../nbs/src/utils.ipynb 25
+# %% ../nbs/src/utils.ipynb 24
 class HierarchicalPlot:
     """Hierarchical Plot
 
@@ -523,7 +523,7 @@ def plot_hierarchical_predictions_gap(
         plt.grid()
         plt.show()
 
-# %% ../nbs/src/utils.ipynb 46
+# %% ../nbs/src/utils.ipynb 45
 # convert levels to output quantile names
 def level_to_outputs(level: list[int]) -> tuple[list[float], list[str]]:
     """Converts list of levels into output names matching StatsForecast and NeuralForecast methods.
@@ -568,7 +568,7 @@ def quantiles_to_outputs(quantiles: list[float]) -> tuple[list[float], list[str]
             output_names.append("-median")
     return quantiles, output_names
 
-# %% ../nbs/src/utils.ipynb 47
+# %% ../nbs/src/utils.ipynb 46
 # given input array of sample forecasts and inptut quantiles/levels,
 # output a Pandas Dataframe with columns of quantile predictions
 def samples_to_quantiles_df(
@@ -653,7 +653,7 @@ def samples_to_quantiles_df(
 
     return _quantiles, df_nw.to_native()
 
-# %% ../nbs/src/utils.ipynb 55
+# %% ../nbs/src/utils.ipynb 53
 # Masked empirical covariance matrix
 @njit(
     "Array(float64, 2, 'F')(Array(float64, 2, 'C'), Array(bool_, 2, 'C'))",
@@ -691,7 +691,7 @@ def _ma_cov(residuals: np.ndarray, not_nan_mask: np.ndarray):
 
     return W
 
-# %% ../nbs/src/utils.ipynb 56
+# %% ../nbs/src/utils.ipynb 54
 # Shrunk covariance matrix using the Schafer-Strimmer method
 
 
@@ -842,7 +842,7 @@ def _shrunk_covariance_schaferstrimmer_with_nans(
 
     return W
 
-# %% ../nbs/src/utils.ipynb 58
+# %% ../nbs/src/utils.ipynb 56
 # Lasso cyclic coordinate descent
 @njit(
     "Array(float64, 1, 'C')(Array(float64, 2, 'C'), Array(float64, 1, 'C'), float64, int64, float64)",

nbs/examples/AustralianDomesticTourism-Intervals.ipynb

@@ -37,9 +37,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import os\n",
-    "os.environ['NIXTLA_ID_AS_COL'] = '1'\n",
-    "\n",
     "import pandas as pd\n",
     "\n",
     "# compute base forecast no coherent\n",

nbs/examples/AustralianDomesticTourism-Permbu-Intervals.ipynb

@@ -37,9 +37,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import os\n",
-    "os.environ['NIXTLA_ID_AS_COL'] = '1'\n",
-    "\n",
     "import pandas as pd\n",
     "\n",
     "# compute base forecast no coherent\n",

nbs/examples/AustralianDomesticTourism.ipynb

@@ -66,9 +66,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import os\n",
-    "os.environ['NIXTLA_ID_AS_COL'] = '1'\n",
-    "\n",
     "import numpy as np\n",
     "import pandas as pd"
    ]

nbs/examples/AustralianPrisonPopulation.ipynb

@@ -63,9 +63,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import os\n",
-    "os.environ['NIXTLA_ID_AS_COL'] = '1'\n",
-    "\n",
     "import numpy as np\n",
     "import pandas as pd"
    ]

nbs/examples/Introduction.ipynb

@@ -161,9 +161,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import os\n",
-    "os.environ['NIXTLA_ID_AS_COL'] = '1'\n",
-    "\n",
     "import numpy as np\n",
     "import pandas as pd"
    ]