From 517a2fe1d52af9bef45a736caa63b5b39d237a2d Mon Sep 17 00:00:00 2001
From: farid <farid.mohammadi@iws.uni-stuttgart.de>
Date: Tue, 3 Aug 2021 14:59:37 +0200
Subject: [PATCH] [surrogate] new normalization function.

---
 .../surrogate_models/RegressionFastARD.py     | 94 +++++++++++--------
 1 file changed, 56 insertions(+), 38 deletions(-)

diff --git a/BayesValidRox/surrogate_models/RegressionFastARD.py b/BayesValidRox/surrogate_models/RegressionFastARD.py
index f1a26bcfa..460381b1b 100755
--- a/BayesValidRox/surrogate_models/RegressionFastARD.py
+++ b/BayesValidRox/surrogate_models/RegressionFastARD.py
@@ -13,7 +13,7 @@ from sklearn.linear_model._base import LinearModel
 import warnings
 from sklearn.utils import check_X_y,check_array,as_float_array
 from scipy.linalg import pinvh
-
+from sklearn.preprocessing import normalize as f_normalize
 
 def update_precisions(Q,S,q,s,A,active,tol,n_samples,clf_bias):
     '''
@@ -54,7 +54,8 @@ def update_precisions(Q,S,q,s,A,active,tol,n_samples,clf_bias):
     
     # changes in precision for features already in model is below threshold
     no_delta       = np.sum( abs( Anew - Arec ) > tol ) == 0
-    
+    # if same_features: print(abs( Anew - Arec ))
+    # print("same_features = {} no_delta = {}".format(same_features,no_delta))
     # check convergence: if no features to add or delete and small change in 
     #                    precision for current features then terminate
     converged = False
@@ -153,41 +154,52 @@ class RegressionFastARD(LinearModel,RegressorMixin):
     '''
     
     def __init__( self, n_iter=300, start=None, tol=1e-3, fit_intercept=True, 
-                  copy_X=True, compute_score=False, verbose=False):
+                  normalize=False, copy_X=True, compute_score=False, verbose=False):
         self.n_iter          = n_iter
         self.start           = start
         self.tol             = tol
         self.scores_         = list()
         self.fit_intercept   = fit_intercept
+        self.normalize       = normalize
         self.copy_X          = copy_X
         self.compute_score   = compute_score
         self.verbose         = verbose
         
         
-    def _center_data(self,X,y):
-        ''' Centers data'''
-        X     = as_float_array(X,self.copy_X)
-        # normalisation should be done in preprocessing!
-        X_std = np.ones(X.shape[1], dtype = X.dtype)
+    def _preprocess_data(self,X, y):
+        """Center and scale data.
+        Centers data to have mean zero along axis 0. If fit_intercept=False or if
+        the X is a sparse matrix, no centering is done, but normalization can still
+        be applied. The function returns the statistics necessary to reconstruct
+        the input data, which are X_offset, y_offset, X_scale, such that the output
+            X = (X - X_offset) / X_scale
+        X_scale is the L2 norm of X - X_offset. 
+        """
+        
+        if self.copy_X:
+            X = X.copy(order='K')
+    
+        y = np.asarray(y, dtype=X.dtype)
+    
         if self.fit_intercept:
-            X_mean = np.average(X,axis = 0)
-            y_mean = np.average(y,axis = 0)
-            X     -= X_mean
-            y     -= y_mean
+            X_offset = np.average(X, axis=0)
+            X -= X_offset
+            if self.normalize:
+                X, X_scale = f_normalize(X, axis=0, copy=False,
+                                         return_norm=True)
+            else:
+                X_scale = np.ones(X.shape[1], dtype=X.dtype)
+            y_offset = np.average(y, axis=0)
+            y = y - y_offset
         else:
-            X_mean = np.zeros(X.shape[1],dtype = X.dtype)
-            y_mean = 0. if y.ndim == 1 else np.zeros(y.shape[1], dtype=X.dtype)
-            Y_std = np.ones(y.shape, dtype = y.dtype)
-            # TODO
-            # X_mean = np.average(X,axis = 0)
-            # X_std = np.std(X,axis = 0)
-            # X_std[0]=np.finfo(np.float32).eps
-            # X     = (X - X_mean) / X_std
-            # y_mean = np.average(y,axis = 0)
-            # Y_std = np.std(y,axis = 0)
-            # y     = (y - y_mean) / Y_std
-        return X,y, X_mean, y_mean, X_std, Y_std
-        
+            X_offset = np.zeros(X.shape[1], dtype=X.dtype)
+            X_scale = np.ones(X.shape[1], dtype=X.dtype)
+            if y.ndim == 1:
+                y_offset = X.dtype.type(0)
+            else:
+                y_offset = np.zeros(y.shape[1], dtype=X.dtype)
+    
+        return X, y, X_offset, y_offset, X_scale  
         
     def fit(self,X,y):
         '''
@@ -209,11 +221,10 @@ class RegressionFastARD(LinearModel,RegressorMixin):
         X, y = check_X_y(X, y, dtype=np.float64, y_numeric=True)
         n_samples, n_features = X.shape
         
-        X, y, X_mean, y_mean, X_std, Y_std = self._center_data(X, y)
+        X, y, X_mean, y_mean, X_std = self._preprocess_data(X, y)
         self._x_mean_ = X_mean
         self._y_mean  = y_mean
         self._x_std   = X_std
-        self._y_std   = Y_std
 
         #  precompute X'*Y , X'*X for faster iterations & allocate memory for
         #  sparsity & quality vectors
@@ -285,13 +296,13 @@ class RegressionFastARD(LinearModel,RegressorMixin):
             rss     = np.sum( ( y - np.dot(X[:,active] , Mn) )**2 )
 
             # if near perfect fit , then terminate
-            if (rss / n_samples/var_y) < self.tol:
-                warnings.warn('Early termination due to near perfect fit')
-                converged = True
-                break
+            # if (rss / n_samples/var_y) < self.tol:
+            #     warnings.warn('Early termination due to near perfect fit')
+            #     converged = True
+            #     break
             beta    = n_samples - np.sum(active) + np.sum(Aa * Sdiag )
             beta   /= ( rss + np.finfo(np.float32).eps )
-            # print("Active terms: {0} Beta: {1:.3e} error: {2:.3e}".format(np.sum(active), beta, rss / n_samples/var_y))
+
             # update precision parameters of coefficients
             A,converged  = update_precisions(Q,S,q,s,A,active,self.tol,
                                              n_samples,False)
@@ -323,7 +334,14 @@ class RegressionFastARD(LinearModel,RegressorMixin):
         self.converged     = converged
         if self.compute_score:
             self.scores_       = np.array(scores_)
-        self._set_intercept(X_mean,y_mean,X_std)
+        
+        # set intercept_
+        if self.fit_intercept:
+            self.coef_ = self.coef_ / X_std
+            self.intercept_ = y_mean - np.dot(X_mean, self.coef_.T)
+        else:
+            self.intercept_ = 0.
+        
         return self
     
     def log_marginal_likelihood(self, XXa,XYa, Aa, beta):
@@ -368,15 +386,15 @@ class RegressionFastARD(LinearModel,RegressorMixin):
         ----------
         [1] Bishop, C. M. (2006). Pattern recognition and machine learning. springer.
         '''
-        x         = (X - self._x_mean_) / self._x_std
-        y_hat     = np.dot(x,self.coef_) * self._y_std + self._y_mean
-        # y_hat     = np.dot(x,self.coef_) * self._y_mean
+        
+        y_hat     = np.dot(X,self.coef_) + self.intercept_
         
         if return_std:
+            if self.normalize:
+                x   = (X - self._x_mean_) / self._x_std
             var_hat   = 1./self.alpha_
             var_hat  += np.sum( np.dot(x[:,self.active_],self.sigma_) * x[:,self.active_], axis = 1)
-            var_hat  *= self._y_std**2
-            std_hat = np.sqrt(var_hat)
+            std_hat   = np.sqrt(var_hat)
             return y_hat, std_hat
         else:
             return y_hat
-- 
GitLab