From f7782b2a689176dcaec9bd66b7044bef011490b1 Mon Sep 17 00:00:00 2001
From: Lukas Novak <95358264+NovakLBUT@users.noreply.github.com>
Date: Fri, 8 Apr 2022 03:35:15 +0200
Subject: [PATCH] Efficient hyperbolic set added

---
 .../polynomials/baseclass/PolynomialBasis.py  | 42 ++++++++++++++++++-
 1 file changed, 41 insertions(+), 1 deletion(-)

diff --git a/src/UQpy/surrogates/polynomial_chaos/polynomials/baseclass/PolynomialBasis.py b/src/UQpy/surrogates/polynomial_chaos/polynomials/baseclass/PolynomialBasis.py
index 8f44cf722..649e0f28c 100644
--- a/src/UQpy/surrogates/polynomial_chaos/polynomials/baseclass/PolynomialBasis.py
+++ b/src/UQpy/surrogates/polynomial_chaos/polynomials/baseclass/PolynomialBasis.py
@@ -104,7 +104,47 @@ def calculate_total_degree_recursive(N, w, rows):
                     row_start = row_end + 1
 
         return subset
-
+    
+    @staticmethod
+    def calculate_hyperbolic_set(inputs_number, degree,q):
+
+        xmono=np.zeros(inputs_number)
+        X=[]
+        X.append(xmono)
+        
+        while np.sum(xmono)<=degree:
+            # generate multi-indices one by one
+            x=np.array(xmono)
+            i = 0
+            for j in range ( inputs_number, 0, -1 ):
+                if ( 0 < x[j-1] ):
+                    i = j
+                    break
+
+            if ( i == 0 ):
+                x[inputs_number-1] = 1
+                xmono=x
+            else:
+                if ( i == 1 ):
+                    t = x[0] + 1
+                    im1 = inputs_number
+                if ( 1 < i ):
+                    t = x[i-1]
+                    im1 = i - 1
+
+                x[i-1] = 0
+                x[im1-1] = x[im1-1] + 1
+                x[inputs_number-1] = x[inputs_number-1] + t - 1
+
+                xmono=x
+                
+            # check the hyperbolic criterion          
+            if (np.round(np.sum(xmono**q)**(1/q), 4) <= degree):
+                X.append(xmono)
+
+
+        return(np.array(X).astype(int))
+    
     @staticmethod
     def calculate_tensor_product_set(inputs_number, degree):
         orders = np.arange(0, degree + 1, 1).tolist()