Diffusion solar_potential (master) Browse

RF_uncertainties_test.py
f18b667cc9e6master

RF_uncertainties_test.py
View Options

			import numpy as np
			import pandas as pd
			import xarray as xr
			import os
			import sys
			import hpelm
			import util
			import csv
			from tqdm import tqdm
			from ds import Dataset
			from sklearn.ensemble import RandomForestRegressor
			from sklearn.metrics import mean_squared_error as mse
			from tables import open_file, Atom, Filters
			from meteo_data import Meteo_Reader
			from sklearn.externals import joblib

			def mse_norm(y,t,norm):
			return mse(norm.normalize(y),norm.normalize(t))

			def find_name(filename):
			path, name = os.path.split(filename)
			root, ext = os.path.splitext(name)
			i = 0
			while os.path.exists(filename):
			i += 1
			filename = os.path.join(path, str(i) + '_' + root + ext)
			return filename

			# # Random Forest with uncertainties
			# In the script, a random forest is trained from a training data set, the MSE is computed from the training and test datasets, and the prediction is made for each tree seperately, in order to later compute the mean and the variance of the prediction. This script can only be run if all training data fits into memory. After prediction, the residuals based on the out-of-bag error are computed, and a new forest 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]

			# ## Setup
			# ### Inputs

			data_path = sys.argv[1]
			# data_path='/Users/alinawalch/Documents/EPFL/data/meteo'
			dsname = sys.argv[2]
			queryname = sys.argv[3]

			n_est = int(sys.argv[4])
			m = int(sys.argv[5])
			tree_depth = 100
			save_model = False


			# ### Load dataset and training, validation, testing and query matrices

			my_ds = Dataset(data_path, dsname, queryname)

			# find number of features and number of targets
			nf = my_ds.feature_norm.m
			nt = my_ds.target_norm.m

			# Merge training and validation data for crossvalidation training data and rescale (matrices give normalised data)
			x = util.get_matrix(my_ds.train_X0)
			t = util.get_matrix(my_ds.train_T0)
			test_x = util.get_matrix(my_ds.test_X0)
			test_t = util.get_matrix(my_ds.test_T0)
			if nt == 1:
			t=t.reshape((-1,))
			test_t=test_t.reshape((-1,))

			# ## Random Forest

			# Add model to dataset folder structure
			modelname = 'RF' + str(n_est)
			my_ds.add_model(modelname, queryname)

			forest = RandomForestRegressor(n_estimators = n_est, max_depth = tree_depth, oob_score = True, n_jobs = -1, max_features = m)

			# ### Forest fitting

			timer_0 = util.Timer()
			forest.fit(x, t)
			timer_0.stop()

			logfile = find_name(os.path.join(my_ds.model_path,('log_time.csv')))
			with open(logfile, 'w') as csvfile:
			w = csv.writer(csvfile, delimiter=',')
			w.writerow(['', 'walltime', 'cputime'])
			w.writerow(['model_train', timer_0.walltime, timer_0.cputime])

			if save_model:
			joblib.dump(tree, os.path.join(my_ds.model_path, 'model.pkl'))

			# ### Model uncertainty

			# Declare matrices for sums and sum-of-squares (do not save all n_est outputs)
			N = t.shape[0]
			N_te = test_t.shape[0]
			if nt == 1:
			sigma = np.zeros(N)
			sigma_sq = np.zeros(N)
			else:
			sigma = np.zeros((N,nt))
			sigma_sq = np.zeros((N,nt))

			if nt == 1:
			sigma_te = np.zeros(N_te)
			else:
			sigma_te = np.zeros((N_te,nt))
			k = 1.0

			# Predict for each tree seperately and update the sums and sum-of-squares

			logfile_pred = find_name(os.path.join(my_ds.model_path,('log_pred_model.csv')))
			with open(logfile_pred, 'w') as csvfile:
			w = csv.writer(csvfile, delimiter=',')
			w.writerow(['model', 'mse_tr', 'mse_te'])

			t1 = util.Timer()
			i = 0
			for tree in tqdm(forest.estimators_):
			y = tree.predict(x)
			sigma += y/k
			sigma_sq += (y/k)**2

			y_te = tree.predict(test_x)
			sigma_te += y_te/k

			mse_tr = mse_norm(sigma * k/(i+1),t,my_ds.target_norm)
			mse_te = mse_norm(sigma_te * k/(i+1),test_t,my_ds.target_norm)

			with open(logfile_pred, 'a') as csvfile:
			w = csv.writer(csvfile, delimiter=',')
			w.writerow([i, mse_tr, mse_te])
			i+=1

			#y_q = tree.predict(query_x)
			#sigma_q += y_q/k
			#sigma_sq_q += (y_q/k)**2
			#i += 1

			t1.stop(print_any = False)

			with open(logfile, 'a') as csvfile:
			w = csv.writer(csvfile, delimiter=',')
			w.writerow(['model_pred', t1.walltime, t1.cputime])

			# #### Compute output (mean) and model uncertainty (var) for training data and for query data

			mean = sigma * k/n_est
			var = k*2/n_est(sigma_sq-sigma**2/n_est)
			del sigma
			del sigma_sq

			# #### Obtain out-of-bag and testing outputs and compare mean squared errors (normalized, for comparability with ELM)

			oob = forest.oob_prediction_
			y_te = forest.predict(test_x)

			print('Normalised MSE (testing): %f' %mse_norm(y_te,test_t,my_ds.target_norm))
			print('Normalised MSE (out-of-bag): %f' %mse_norm(oob,t,my_ds.target_norm))
			print('Normalised MSE (training): %f' %mse_norm(mean,t,my_ds.target_norm))

			logpred = find_name(os.path.join(my_ds.model_path,('log_pred_mse.csv')))
			with open(logpred, 'w') as csvfile:
			w = csv.writer(csvfile, delimiter=',')
			w.writerow(['', 'mse'])
			w.writerow(['testing', mse_norm(y_te,test_t,my_ds.target_norm)])
			w.writerow(['out-of-bag', mse_norm(oob,t,my_ds.target_norm)])
			w.writerow(['training', mse_norm(mean,t,my_ds.target_norm)])

solarPOT/Solar/RF_uncertainties_test.pyf18b667cc9e6master

RF_uncertainties_test.pyView Options

solarPOT/Solar/RF_uncertainties_test.py
f18b667cc9e6master

RF_uncertainties_test.py
View Options