added extended version to manocca

python/src/explainer.py

@@ -219,7 +219,7 @@ class Explainer :
         if return_raw_contrib == True :
             return res
         else :
-            df_loadings = res.iloc[:,:-3]
+            df_loadings = res.iloc[:,:-3] #-3 because we added p, beta and chi2 at the end
             df_loadings = df_loadings*df_loadings
             df_prod_contrib = res["chi2"].values.reshape(-1,1)*df_loadings
             return df_prod_contrib.sum().sort_values(ascending = False)#.to_dict()

python/src/manocca.py

@@ -99,7 +99,7 @@ class MANOCCA:
 
     """
     def __init__(self, predictors, outputs, covariates=None, cols_outputs = None,
-                 cols_predictors = None, cols_covariates = None, prodV_red=None, n_comp = None, prod_to_keep = None, use_resid = True, use_pca = True, n_jobs = 1):
+                 cols_predictors = None, cols_covariates = None, prodV_red=None, n_comp = None, prod_to_keep = None, use_resid = True, use_pca = True, use_extended = False, n_jobs = 1):
 
         ### Initializing 
         self.outputs = outputs
@@ -126,6 +126,7 @@ class MANOCCA:
         self.n_jobs = n_jobs
         self.use_resid = use_resid
         self.use_pca = use_pca
+        self.use_extended = use_extended
 
         # Filled later
         self.prodV = None
@@ -148,7 +149,10 @@ class MANOCCA:
                 self.predictors = np.apply_along_axis(lambda x : pt.adjust_covariates(x,covariates), axis = 0, arr = self.predictors)
 
         else : # If not we preprocess the data and compute prodV and prodV_red
-            self.prodV = self.get_prodV_para_wrap(self.outputs)
+            if self.use_extended : 
+                self.prodV = self.get_prodV_extended_wrap(self.outputs)
+            else :
+                self.prodV = self.get_prodV_para_wrap(self.outputs)
             # self.prodV = scale(self.prodV)
             if not isinstance(self.prod_to_keep, type(None)): # Filtering out some columns
                 self.filter_prodV_columns()
@@ -267,6 +271,35 @@ class MANOCCA:
         else :
             return tmp
 
+    ### For extended ###
+    def get_prodV_extended(self, DD0, job_id, nb_compute):
+        L_prodV = []
+        for i in range(job_id, DD0.shape[1], nb_compute):
+    #         tmp = np.transpose(np.transpose(DD0[:,(i+1):DD0.shape[1]])*DD0[:,i])
+            tmp = (DD0[:,(i):]*DD0[:,i].reshape(-1,1))
+            tmp = pt.get_qt(tmp)
+            L_prodV += [tmp]
+        return L_prodV
+
+    def get_prodV_extended_wrap(self, DD0, verbose = 10):
+        n_jobs = self.n_jobs
+        if n_jobs == -1:
+            nb_compute = cpu_count()
+        else :
+            nb_compute = min(cpu_count(),n_jobs)
+        print("Computing prodV with %i cpu" %nb_compute)
+
+        res = Parallel(n_jobs=n_jobs, verbose = verbose)(delayed(self.get_prodV_extended)(DD0,j, nb_compute) for j in range(nb_compute))
+        
+        # reordering
+        res_ordered = []
+        for i in range(len(res[0])):
+            for L in res :
+                if L != [] : 
+                    res_ordered += [L.pop(0)]
+        return np.hstack([DD0]+res_ordered)
+
+    ### Other functions ###
     def get_prod_cols(self, cols, sep = '|'):
         return [cols[i]+sep+cols[j] for i in range(len(cols)-1) for j in range(i+1,len(cols))]
 

python/src/tools/compute_manova.py

@@ -20,6 +20,7 @@ def custom_manova(Y,X,C=None, return_beta = False):
 
     # (sign_num, num) = np.linalg.slogdet(dot_Y-X.shape[0]*np.dot(beta,beta.T))
     (sign_num, num) = np.linalg.slogdet(dot_Y - np.dot( np.dot(beta, X.T) , np.dot(X, beta.T)))
+    # (sign_num, num) = np.linalg.slogdet((Y - X @ beta.T).T@(Y - X @beta.T))
     (sign_denom, denom) = np.linalg.slogdet(dot_Y)
     lamb = np.exp(sign_num*num-(sign_denom*denom))
     # print(lamb)
@@ -81,6 +82,8 @@ def linear_regression(Y, X, C=None):
         beta = np.dot(np.linalg.inv(np.dot(X.T,X)) , np.dot(X.T,Y))
         beta = beta.T
         # beta = np.dot(X.T,Y)/X.shape[0]
+        # print(beta.shape)
+        # beta = beta.T
     return beta