Added check on matrix A, if it is hierarchical the integer linear pro…

…grammin is not run, fixes issue #7.
IDSIA · Nov 27, 2023 · ed158a8 · ed158a8
1 parent 600549b
commit ed158a8
Show file tree

Hide file tree

Showing 2 changed files with 52 additions and 10 deletions.
diff --git a/R/reconc.R b/R/reconc.R
@@ -218,12 +218,27 @@ reconc_BUIS <- function(S,
   bottom_base_forecasts = split_hierarchy.res$bottom
 
   # H, G
-  get_HG.res = .get_HG(A, upper_base_forecasts, distr[split_hierarchy.res$upper_idxs], in_type[split_hierarchy.res$upper_idxs])
-  H = get_HG.res$H
-  upper_base_forecasts_H = get_HG.res$Hv
-  G = get_HG.res$G
-  upper_base_forecasts_G = get_HG.res$Gv
-
+  if(.check_hierarchical(A)){
+    H = A
+    G = NULL
+    upper_base_forecasts_H = upper_base_forecasts
+    upper_base_forecasts_G = NULL
+    in_typeH = in_type[split_hierarchy.res$upper_idxs]
+    distr_H  = distr[split_hierarchy.res$upper_idxs]
+    in_typeG = NULL
+    distr_G  = NULL
+  }else{
+    get_HG.res = .get_HG(A, upper_base_forecasts, distr[split_hierarchy.res$upper_idxs], in_type[split_hierarchy.res$upper_idxs])
+    H = get_HG.res$H
+    upper_base_forecasts_H = get_HG.res$Hv
+    G = get_HG.res$G
+    upper_base_forecasts_G = get_HG.res$Gv
+    in_typeH = get_HG.res$Hin_type
+    distr_H  = get_HG.res$Hdistr
+    in_typeG = get_HG.res$Gin_type
+    distr_G  = get_HG.res$Gdistr
+  }
+
   # Reconciliation using BUIS
   n_upper = nrow(A)
   n_bottom = ncol(A)
@@ -249,8 +264,8 @@ reconc_BUIS <- function(S,
       b = (B %*% c),
       # (num_samples x 1)
       u = unlist(upper_base_forecasts_H[[hi]]),
-      in_type_ = get_HG.res$Hin_type[[hi]],
-      distr_ = get_HG.res$Hdistr[[hi]]
+      in_type_ = in_typeH[[hi]],
+      distr_ = distr_H[[hi]]
     )
     B[, b_mask] = .resample(B[, b_mask], weights)
   }
@@ -263,8 +278,8 @@ reconc_BUIS <- function(S,
       weights = weights * .compute_weights(
         b = (B %*% c),
         u = unlist(upper_base_forecasts_G[[gi]]),
-        in_type_ = get_HG.res$Gin_type[[gi]],
-        distr_ = get_HG.res$Gdistr[[gi]]
+        in_type_ = in_typeG[[gi]],
+        distr_ = distr_G[[gi]]
       )
     }
     B = .resample(B, weights)

diff --git a/R/utils.R b/R/utils.R
@@ -2,6 +2,7 @@
 .DISTR_SET = c("gaussian", "poisson", "nbinom")
 .DISTR_SET2 = c("continuous", "discrete")
 .check_input <- function(S, base_forecasts, in_type, distr) {
+
   if (!(nrow(S) == length(base_forecasts))) {
     stop("Input error: nrow(S) != length(base_forecasts)")
   }
@@ -84,6 +85,32 @@
   }
 }
 
+# Returns TRUE if A is a hierarchy matrix 
+# (according to Definition 1 in "Find Maximal Hierarchy")
+# If this returns TRUE we avoid solving the integer linear programming problem
+.check_hierarchical <- function(A) {
+
+  k <- nrow(A)
+  m <- ncol(A)
+
+  for (i in 1:k) {
+    for (j in 1:k) {
+      if (i < j) {
+        cond1 = A[i,] %*% A[j,] != 0     # Upper i and j have some common descendants
+        cond2 = sum(A[i,] - A[j,] >= 0) < m  # Upper j is not a descendant of upper i
+        cond3 = sum(A[i,] - A[j,] <= 0) < m  # Upper i is not a descendant of upper j
+        if (cond1 & cond2 & cond3) {
+          return(FALSE)
+        }
+      }
+    }
+  }
+
+  return(TRUE)
+
+}
+
+
 # Misc
 .shape <- function(m) {
   print(paste0("(", nrow(m), ",", ncol(m), ")"))