Homec4science
Diffusion solar_potential (rc_energy_potential) Browse

solarPOT/Solar/ELM_uncertainties_query.py
dd1c0248d00drc_energy_potential

Last Change
Raw File

ELM_uncertainties_query.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]
queryname = sys.argv[3]
n_est = int(sys.argv[4])
n_nodes = int(sys.argv[5])
accel = sys.argv[6]
#################################################################
# ### Load dataset and training matrix (for bootstrapping), validation and testing (for direct error computation)
my_ds = Dataset(data_path, dsname, queryname)
nf = my_ds.feature_norm.m
nt = my_ds.target_norm.m
if accel is None: accel = 'CPU'
modelname = 'ELM' + str(n_est) + '_' + str(n_nodes) + '_' + accel
my_ds.add_model(modelname, queryname)
modelpath = os.path.join(my_ds.model_path,'model')
ens = HPELM_Ensemble(path = modelpath, n_estimators = n_est, n_nodes = n_nodes, n_features = nf, n_targets = nt, oob = True, accelerator = accel)
# ### Model fitting
# If model has been fitted before and to save computational time, uncomment *ens.load()* and comment out the fitting.
# Loading model
ens.load()
# #### Prediction on query dataset
# predict on query set:
ens.predict(my_ds.query_X, my_ds.query_Y, label = 'query' , eval = True, norm = my_ds.target_norm)
query_time = ens.prediction_times_
# ## Residuals
# #### Declare and load model
residualpath = os.path.join(my_ds.model_path,'residuals')
res_ens = HPELM_Ensemble(path = residualpath, n_estimators = n_est, n_nodes = n_nodes, n_features = nf, n_targets = nt, oob = False, accelerator = accel)
res_ens.load()
# ### Predict on residual model
res_norm = util.Normalizer(os.path.join(my_ds.model_path, "norm_residuals.csv"))
# #### Predict on query grid
res_query_Y = os.path.join(my_ds.query_path_out, 'residual_prediction.hdf5')
res_ens.predict(my_ds.query_X, res_query_Y, eval = True, norm = res_norm, label = 'query') # eval computes directly the mean and the variance of the model
res_query_time = res_ens.prediction_times_
query_y_var_f = os.path.splitext(my_ds.query_Y)[0] + '_var.hdf5'
try:
X = util.get_matrix(my_ds.query_X0)
Y = util.get_matrix(my_ds.query_Y)
var_model = util.get_matrix(query_y_var_f)
var_data = util.get_matrix(os.path.join(my_ds.query_path_out, 'residual_prediction.hdf5'))
var_tot = var_model + var_data
# ## Evaluation
my_ds.make_xarray(X, Y, rescale = False).to_netcdf(os.path.join(my_ds.query_path_out, 'prediction.nc'))
my_ds.make_xarray(X, np.sqrt(np.maximum(var_model,0)), rescale = False).to_netcdf(os.path.join(my_ds.query_path_out, 'model_unc.nc'))
my_ds.make_xarray(X, np.sqrt(np.maximum(var_data,0)), rescale = False).to_netcdf(os.path.join(my_ds.query_path_out, 'data_unc.nc'))
my_ds.make_xarray(X, np.sqrt(np.maximum(var_tot,0)), rescale = False).to_netcdf(os.path.join(my_ds.query_path_out, 'total_unc.nc'))
except:
print('failed to load query data and to create xarray datasets - may be due to out-of-memory issue')
c4science · Help