Homec4science
Diffusion solar_potential (master) Browse

solarPOT/Solar/ELM_compute_residuals.py
f18b667cc9e6master

Last Change
Raw File

ELM_compute_residuals.py
View Options

# coding: utf-8
import numpy as np
import pandas as pd
import xarray as xr
import os
import hpelm
import util
import h5py
import sys
from tqdm import tqdm
from ds import Dataset
from ELM_ensemble import HPELM_Ensemble
from sklearn.metrics import mean_squared_error as mse
from tables import open_file, Atom, Filters
# # Extreme Learning Machine Ensemble with uncertainties
# The script makes use of the class **ELM_ensemble**, which is designed in analogy to the scikit-learn structure.
# An ensemble of extreme learning machines is fitted with the training data (and the validation data if 'V' mode is selected). While training, the out-of-bag samples are stored for later reconstruction of the out-of-bag prediction. After fitting, the model is used to predict on the training data and on the query data (storing the outputs of each run). After prediction, the residuals based on the out-of-bag data are computed, and a new ELM ensemble is trained on the residuals. The variance of the prediction represents the model uncertainty, while the predicted residuals estimate the data uncertainty. [See FG poster]
# _**Notes**_:
# - **Should the validation set be different (random) for each ELM?**
# - We should find a better method to implement the bootstrap - this method does not work for **BIG DATA**
# - Computation of residuals is based on above assumption and requires all data to be loaded into memory
# - Parallelising across the different ELMs (embedded in *ELM_ensemble.py*) should bring significant improvements to expected performance
# - Find out how HPELM implements the parallelisation of the *hdf5_tools* - so that we can use the same approach for the functions in *norms.py*
# - We need a solution for handling very large query datasets!
# ## Setup
############################ Inputs ############################
# data_path = os.path.abspath("/mnt/sda1/hyenergy/data/meteo")
data_path = sys.argv[1]
dsname = sys.argv[2]
modelname = sys.argv[3]
n_est = int(sys.argv[4])
n_nodes = int(sys.argv[5])
accel = str(sys.argv[6])
######### CONSTANTS ############
COMPUTE_RESIDUALS = True
#################################################################
# ### Load dataset and training matrix (for bootstrapping), validation and testing (for direct error computation)
print("\n\nELM residual training with uncertainties\nEnsemble size: %d\nNumber of nodes: %d\n\n" %(n_est, n_nodes))
my_ds = Dataset(data_path, dsname)
nf = my_ds.feature_norm.m
nt = my_ds.target_norm.m
my_ds.get_matrices(['train', 'test', 'val'])
print('\nLoaded dataset %s' %dsname)
if accel is None: accel = 'CPU'
my_ds.add_model(modelname)
# define logile names
logfile_mse = util.find_name(os.path.join(my_ds.model_path,( 'res%d_%d_log_pred_mse.csv' %(n_est, n_nodes) ) ))
logfile_time = util.find_name(os.path.join(my_ds.model_path,( 'res%d_%d_log_time.csv' %(n_est, n_nodes) ) ))
# define file paths for residual files (targets (T) and features (X))
res_T_tr = os.path.join( my_ds.train_path_out, 'train_residuals_t.hdf5' )
res_T_val = os.path.join( my_ds.train_path_out, 'val_residuals_t.hdf5' )
res_T_te = os.path.join( my_ds.test_path_out , 'test_residuals_t.hdf5' )
res_X_tr = os.path.join( my_ds.train_path_out, 'train_residuals_f.hdf5' )
res_X_val = os.path.join( my_ds.train_path_out, 'val_residuals_f.hdf5' )
res_X_te = os.path.join( my_ds.test_path_out , 'test_residuals_f.hdf5' )
# obtain the rescaled out-of-bag prediction, model variance and total model mean based on the predicted training data
if COMPUTE_RESIDUALS:
oob = util.get_matrix( os.path.join(my_ds.train_path_out, 'oob_prediction.hdf5' ) )
mean = util.get_matrix( os.path.join(my_ds.train_path_out, 'train_prediction.hdf5' ) )
var = util.get_matrix( os.path.join(my_ds.train_path_out, 'train_prediction_var.hdf5') )
t0 = my_ds.target_norm.rescale(my_ds.train_t)
res_oob = (oob-t0)**2 - var
res = np.maximum(res_oob, 0)
del res_oob, oob, var, mean, t0
# ### Fit residual model
# #### Prepare residual files
# Split data into training, validation and testing before writing to files in order to obtain test dataset to check model quality:
N = res.shape[0]
ind = np.random.permutation(N)
# keep 10% for testing, and 10% of the remaining data for validation
ind_tr = int( 0.81 * N )
ind_val = int( 0.9 * N )
# In order to be able to train the residual model, we need to write the residuals to a file and compute their normalisation - this is the core functionality of the **Table_Writer** class. The class writes the file and computes mean, norm, max and min, and returns a pandas DataFrame with it, that is then stored in memory:
res_table = util.Table_Writer(my_ds.train_path_out, my_ds.target_norm.names, 'residuals')
res_table.open_hdf5(k = 10.0)
res_table.write_hdf5(res[:ind_tr], table = False)
res_norm = res_table.close_hdf5()
res_norm.to_csv(os.path.join(my_ds.model_path, "norm_residuals.csv"))
# The norm computed above can be loaded into a **Normalizer** object, which contains the functions *normalize* and *rescale* that are useful when working with the data. Below, the residual norm is loaded, and the normalized residuals for training and validation as well as the targets for training, validation and testing are saved in the file given by the path *res_T_XX*. The targets for testing are never needed from memory and do hence not need to be stored.
res_norm = util.Normalizer(os.path.join(my_ds.model_path, "norm_residuals.csv"))
res_norm.set_status('mean')
res = res_norm.normalize(res)
res_tr = res[ : ind_tr ]
res_val = res[ ind_tr : ind_val ]
res_te = res[ ind_val : ]
util.write_hdf5( res_tr , res_T_tr )
util.write_hdf5( res_val, res_T_val )
util.write_hdf5( res_te , res_T_te )
# Now the training dataset has changed, so the features must also be rewritten:
util.write_hdf5( my_ds.train_x[ : ind_tr ], res_X_tr )
util.write_hdf5( my_ds.train_x[ ind_tr : ind_val ], res_X_val )
util.write_hdf5( my_ds.train_x[ ind_val : ], res_X_te )
del res
c4science · Help