From ed158a8987a8716be2af7094ab6ddc5b41c08f38 Mon Sep 17 00:00:00 2001
From: Dario Azzimonti <dario.azzimonti@gmail.com>
Date: Mon, 27 Nov 2023 09:20:52 +0100
Subject: [PATCH] Added check on matrix A, if it is hierarchical the integer
 linear programmin is not run, fixes issue #7.

---
 R/reconc.R | 35 +++++++++++++++++++++++++----------
 R/utils.R  | 27 +++++++++++++++++++++++++++
 2 files changed, 52 insertions(+), 10 deletions(-)

diff --git a/R/reconc.R b/R/reconc.R
index bcd40aa..16ee4bd 100644
--- 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
index 67db281..cb0916c 100644
--- 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), ")"))