Merge pull request #18 from tamnva/development

add feature important

Author: tamnva

hydroecolstm/__init__.py

@@ -1,6 +1,6 @@
-__version__ = "0.3.4"
+__version__ = "0.3.5"
 __author__ = 'Tam V. Nguyen'
 __credits__ = 'Helmholtz Centre for Environmental Research (UFZ)'
 
 from .import data, interface, model, utility, train
-__all__ = ["data", "interface", "model", "utility", "train"]
+__all__ = ["data", "interface", "model", "utility", "train", "feat_importance"]

hydroecolstm/data/read_data.py

@@ -225,4 +225,4 @@ def read_scale_data(config):
     data['y_train_scale'] = data['y_scaler'].transform(x=data["y_train"])
     data['y_valid_scale'] = data['y_scaler'].transform(x=data["y_valid"])  
 
-    return data
+    return data

hydroecolstm/data/scaler.py

@@ -23,37 +23,55 @@ def fit(self, x=None, method=None):
 
         for i, method_name in zip(range(len(method)), method):
             if method_name == "MinMaxScaler":
+                
+                # If cannot scale by MinMaxScaler then do not scale
+                if torch.any(self.maxs[i] - self.mins[i] == 0.0):
+                    print("max-min is zero, cannot use MinMaxScaler, no transform")
+                    self.mins[self.maxs[i] - self.mins[i] == 0.0] = 0.0
+                    self.maxs[self.maxs[i] - self.mins[i] == 0.0] = 1.0
+                
                 scaler_a.append(self.mins[i])
                 scaler_b.append(self.maxs[i] - self.mins[i])
+                
             elif method_name=="Z-score":
+                
+                if torch.any(self.stds[i] == 0.0):
+                    print("standard deviation is zero, cannot use Z-score, no transform")
+                    self.means[self.stds[i] == 0.] = 0.0
+                    self.stds[self.stds[i] == 0.] = 1.0
+                
                 scaler_a.append(self.means[i])
                 scaler_b.append(self.stds[i])
+                
             elif method_name=="None":
                 scaler_a.append(0.0)
                 scaler_b.append(1.0)
+                
             else:
                 print("Error: unknown scaler")
                 SystemExit("Program stop, please change scaler")
 
         scaler_ab = torch.cat((torch.tensor(scaler_a, dtype=torch.float32),
                                torch.tensor(scaler_b, dtype=torch.float32)), 0)
 
-        self.scaler_parameter = torch.reshape(scaler_ab, 
-                                              (2,len(scaler_a)))
+        self.scaler_parameter = torch.reshape(
+            scaler_ab,(2,len(scaler_a)))
 
-    def transform(self, x:dict[str:torch.tensor]=None) -> list: 
+    def transform(self, x:dict[str:torch.tensor]=None) -> list:
         x_scale = {}
         for object_id in x:
-            x_scale[object_id] =  torch.div(torch.sub(x[object_id], 
-                                                      self.scaler_parameter[0,:]), 
-                                            self.scaler_parameter[1,:])               
+            x_scale[object_id] =  torch.div(
+                torch.sub(x[object_id],self.scaler_parameter[0,:]),
+                self.scaler_parameter[1,:])               
         return x_scale
 
     def inverse(self, x:list=None) -> list:        
         x_inverse = {}
         for object_id in x:
-            x_inverse[object_id] =  torch.add(self.scaler_parameter[0,:],
-                                              x[object_id]*self.scaler_parameter[1,:])
+            x_inverse[object_id] =  torch.add(
+                self.scaler_parameter[0,:],
+                x[object_id]*self.scaler_parameter[1,:])
+
         return x_inverse
 
 def _column_mins(input_tensor: torch.tensor=None):
@@ -105,7 +123,8 @@ def get_scaler_name(config):
                 scaler_name_input.append(name)
 
     # scaler name target
-    scaler_name_target = config["scaler_target_features"]*len(config["target_features"])
+    scaler_name_target = config["scaler_target_features"]*len(
+        config["target_features"])
 
     return scaler_name_input, scaler_name_target
 

hydroecolstm/feat_importance/__init__.py

@@ -0,0 +1,3 @@
+from .perm_feat_importance import pfib
+
+__all__= ["pfib"]

hydroecolstm/feat_importance/perm_feat_importance.py

@@ -0,0 +1,68 @@
+
+import copy
+import torch
+import pandas as pd
+from hydroecolstm.utility.evaluation_function import EvaluationFunction
+
+#features=features
+#x_test_scale=data["x_test_scale"]
+#y_test=data["y_test"]
+#y_scaler=data["y_scaler"]
+#trained_model=model
+#objective_function_name="NSE"
+#y_column_name=data["y_column_name"]
+#nskip=config["warmup_length"]
+#seed=100
+            
+# Permutation feature important basin wise
+def pfib(features: str, x_test_scale:dict[str, torch.Tensor], 
+         y_test:dict[str, torch.Tensor], y_scaler, 
+         trained_model, objective_function_name:str, 
+         nskip:int, y_column_name:str, seed:int=None):
+    
+    # Evaluation function
+    objective = EvaluationFunction(function_name=objective_function_name, 
+                                   nskip=nskip, y_column_name=y_column_name)
+    
+    #obj = objective(y_test, 
+    #                y_scaler.inverse(trained_model.evaluate(x_test_scale)))
+    
+    # Loop over features
+    for i in range(len(features)): 
+        x_perm = {}
+        
+        for key, x in zip(x_test_scale.keys(), x_test_scale.values()):
+
+            # Shuffle index of feature i
+            if seed is not None: 
+                torch.manual_seed(0)
+                
+            idx = torch.randperm(x.shape[0])
+            
+            # Shuffle data 
+            x_copy = copy.deepcopy(x)
+            x_copy[:,i] = x_copy[idx, i]
+            
+            # Save permutated data for each key
+            x_perm[key] = copy.deepcopy(x_copy)
+            
+        prediction = y_scaler.inverse(trained_model.evaluate(x_perm))
+        
+        if i == 0: 
+            output = objective(y_test, prediction)
+            output.columns = features[i] + "_" + output.columns 
+            #output.rename(columns={output.columns[0]: features[i]}, 
+            #          inplace=True)
+        else: 
+            temp = objective(y_test, prediction)
+            temp.columns = features[i] + "_" + temp.columns 
+            
+            output = pd.concat([output, temp], axis=1)
+            #output[features[i]] =  objective(
+            #    y_test, prediction)["objective_function_value"]
+        
+        #output.columns = ["s"]
+        
+    return output  #.subtract(obj['objective_function_value'], axis=0)
+
+

hydroecolstm/interface/utility.py

@@ -341,9 +341,12 @@ def config_to_text(config):
                     elif key == "valid_period":
                         out_text.append("  - " + str(config["valid_period"][0])[:16] + "\n")
                         out_text.append("  - " + str(config["valid_period"][1])[:16] + "\n")
-                    else:
+                    elif key == "test_period":
                         out_text.append("  - " + str(config["test_period"][0])[:16] + "\n")
-                        out_text.append("  - " + str(config["test_period"][1])[:16] + "\n")                           
+                        out_text.append("  - " + str(config["test_period"][1])[:16] + "\n") 
+                    else:
+                        out_text.append("  - " + str(config["forecast_period"][0])[:16] + "\n")
+                        out_text.append("  - " + str(config["forecast_period"][1])[:16] + "\n")                           
             except:
                 # Non list object writte in 1 line
                 out_text.append(key +": " + str(config[key]) + "\n")

hydroecolstm/train/custom_loss.py

@@ -8,11 +8,11 @@ def __init__(self, loss_function:str):
         # Dict of all available loss functions
         loss_functions = {"MSE": self.MSE, 
                           "RMSE": self.RMSE,
+                          "RMSE_normalize": self.RMSE_normalize,
                           "MAE": self.MAE,
                           "NSE_complement": self.NSE_complement}
 
         # Use this loss function
-
         self.loss_function = loss_functions[loss_function]
 
     def forward(self, y_true:torch.Tensor, y_predict:torch.Tensor) -> torch.Tensor:
@@ -23,7 +23,10 @@ def forward(self, y_true:torch.Tensor, y_predict:torch.Tensor) -> torch.Tensor:
 
         mask = ~torch.isnan(y_true)
         loss = self.loss_function(y_true, y_predict, mask)
-            
+        
+        if torch.isnan(loss).any():
+            raise ValueError("loss is nan, cannot train the model, check training data")
+
         return loss
 
     # Mean square error
@@ -52,8 +55,18 @@ def RMSE(self, y_true:torch.Tensor, y_predict:torch.Tensor,
 
         # Root Mean Square Error
         rmse = self.MSE(y_true, y_predict, mask)**0.5
+        
         return rmse
 
+    def RMSE_normalize(self, y_true:torch.Tensor, y_predict:torch.Tensor,
+            mask:torch.Tensor)-> torch.Tensor: 
+        
+        # Root Mean Square Error
+        rmse_normalize = self.MSE(y_true, y_predict, mask)**0.5/ torch.mean(
+            y_true[mask])
+        
+        return rmse_normalize
+    
     # Complement to 1 of the Nash-Sutcliffe (or 1- Nash sutcliffe)
     def NSE_complement(self, y_true:torch.Tensor, y_predict:torch.Tensor,
             mask:torch.Tensor)-> torch.Tensor: 
@@ -66,14 +79,3 @@ def NSE_complement(self, y_true:torch.Tensor, y_predict:torch.Tensor,
 
         # Minimize loss, so output should be sse/ssd, which is 1 - NSE
         return sse/ssd
-
-
-#x = CustomLoss(config["loss_function"])
-
-
-
-
-
-
-
-

hydroecolstm/train/trainer.py

@@ -1,12 +1,12 @@
 import numpy as np
 import pandas as pd
 from ray import train
+import ray
 import tempfile
 import os
 import copy
 import torch
 from torch.utils.data import DataLoader
-from pathlib import Path
 from hydroecolstm.train.custom_loss import CustomLoss
 from hydroecolstm.data.custom_dataset import CustomDataset
 
@@ -124,7 +124,7 @@ def train(self,
 
             print(f"Epoch [{epoch+1}/{self.n_epochs}], ", 
                   f"average_train_loss = {train_loss_epoch[-1]:.8f}, ",
-                  f"avearge_valid_loss = {valid_loss_epoch[-1]:.8f}")
+                  f"average_valid_loss = {valid_loss_epoch[-1]:.8f}")
 
             # Early stopping based on validation loss and make checkpoint
             flag = early_stopping(valid_loss_epoch[-1], self.model)
@@ -157,11 +157,6 @@ def train(self,
                                                      valid_loss_epoch, 
                                                      check_point)
 
-        # Save loss_epoch incase of automatic hyperparam optim with tune
-        #self.loss_epoch.to_csv(
-        #    Path(self.out_dir, str(np.random.randint(1, 1e9)) + ".txt"),
-        #    sep='\t')
-
         return self.model
 
     # Save intermediate result at check points
@@ -183,9 +178,10 @@ def _save_check_point(self, train_loss_epoch, valid_loss_epoch,
                                                     valid_loss_epoch, 
                                                     check_point)
 
-            train.report({'loss': train_loss_epoch[-1], 
-                          'loss_epoch': loss_epoch},
-                         checkpoint=checkpoint)
+            if ray.train._internal.session.get_session(): 
+                train.report({'loss': train_loss_epoch[-1],
+                              'loss_epoch': loss_epoch},
+                             checkpoint=checkpoint)
 
     # Create data frame of epoch number, train loss, valid loss
     def _create_train_loss_df(self, train_loss_epoch, valid_loss_epoch, 

hydroecolstm/utility/__init__.py

@@ -1,4 +1,5 @@
 from .evaluation_function import EvaluationFunction
 from .plot import plot
+from .format_conversion import tensor_to_pandas_df
 
 __all__= ["EvaluationFunction", "plot"]

hydroecolstm/utility/evaluation_function.py

@@ -1,7 +1,8 @@
 import torch
+import pandas as pd
 
 class EvaluationFunction():
-    def __init__(self, function_name:str, nskip:int):
+    def __init__(self, function_name:str, nskip:int, y_column_name:str):
 
         # Dict of all available evaluation functions
         evaluation_functions = {"MSE": self.MSE, "RMSE": self.RMSE,
@@ -10,6 +11,8 @@ def __init__(self, function_name:str, nskip:int):
         # Selected evaluation function
         self.eval_function = evaluation_functions[function_name]
         self.nskip = nskip
+        self.function_name = function_name
+        self.y_column_name = y_column_name
 
     def __call__(self, y_true:torch.Tensor, y_predict:torch.Tensor) -> torch.Tensor:
 
@@ -19,16 +22,21 @@ def __call__(self, y_true:torch.Tensor, y_predict:torch.Tensor) -> torch.Tensor:
         for key in y_true.keys():
             eval_values[key] = self.eval_function(y_true[key][self.nskip:,],
                                                   y_predict[key][self.nskip:,])
-
-        avg_eval_values = sum(sum(eval_values.values()))/((len(eval_values))*eval_values[next(iter(eval_values))].shape[0])
 
-        return eval_values, avg_eval_values
+        
+        df = pd.DataFrame(torch.stack(list(eval_values.values())).numpy())
+        df.index = eval_values.keys()
+        df.columns = [self.function_name + "_" + name 
+                      for name in self.y_column_name]
+        
+        return df
 
     def MSE(self, ytrue:torch.Tensor, ypredict:torch.Tensor):
         mask = ~torch.isnan(ytrue)
         mse = []
         for i in range(ytrue.shape[1]):
-            mse.append(torch.mean((ytrue[:,i][mask[:,i]] - ypredict[:,i][mask[:,i]])**2))
+            mse.append(torch.mean((ytrue[:,i][mask[:,i]] - 
+                                   ypredict[:,i][mask[:,i]])**2))
         mse = torch.stack(mse)
         return mse
 
@@ -38,23 +46,20 @@ def RMSE(self, ytrue:torch.Tensor, ypredict:torch.Tensor):
         rmse = mse**0.5
         return rmse
 
-    # 1 - Nash–Sutcliffe efficiency (NSE)
+    # Nash–Sutcliffe efficiency (NSE)
     def NSE(self, ytrue:torch.Tensor, ypredict:torch.Tensor):
         mask = ~torch.isnan(ytrue)
 
         # Sum of Square Error (sse) = sum((true-predict)^2)
         # Sum of Square Difference around mean (ssd) = sum((true-mean_true)^2)
         sse = []        
         ssd = []
+        
         for i in range(ytrue.shape[1]):
             sse.append(torch.sum((ytrue[:,i][mask[:,i]] - ypredict[:,i][mask[:,i]])**2))
             ssd.append(torch.sum((ytrue[:,i][mask[:,i]] - torch.nanmean(ytrue[:,i]))**2))
 
-        # get 1 - nse, here I call it as nse
         nse = 1.0 - torch.stack(sse)/torch.stack(ssd)
-        
-        if torch.isnan(nse).any():
-            raise ValueError("nan values found when calculating NSE - zero division")
 
         return nse