Upload New File

IGARSS/QSVR.py

+# -*- coding: utf-8 -*-
+"""
+@author: Edoardo Pasetto
+
+"""
+
+import os
+import numpy as np
+from sklearn.svm import SVR
+from sklearn.model_selection import train_test_split
+from sklearn.metrics import mean_squared_error
+
+
+
+
+
+import utility
+
+from datetime import datetime
+
+
+
+#to use a different dataset just change the path
+
+path='Dataset\\seabam919.txt' 
+dataset=np.loadtxt(path)
+
+
+
+epsilon=0.1
+
+
+#here the parameters of the models are set
+
+n_samples=25
+K=3         
+B=float(2.5)  
+beta=0
+xi=0
+k0=2   
+MAXRESULTS=40
+
+
+
+#same validation is set to True when the dataset used for training is used for the solutions combination methods
+same_validation=True
+
+
+
+#change the random seed to change how the dataseta are randomly chosen 
+
+random_seed=0
+
+
+exponents_array=np.array(list(range(K)))
+exponents_array=exponents_array-k0
+
+#C_max is the minimum value that C can take
+
+C_max=sum(np.power(B,exponents_array))
+
+#the arrays gamma_values and C_values are used for the hyperparameters validation
+
+gamma_values=[0.1,0.5,1,1.5,2,3,4,5,7,10,20,50]
+C_values=[C_max,2*C_max, 5*C_max,10*C_max,20*C_max,50*C_max]
+
+
+
+
+
+#here the features of the training dataset are selected
+#the chosen features are the one corresponding to the selected wavelengths
+ 
+features=[6,7,8,9,12]
+X=dataset[:,features]
+#%%
+data_indexes=[i for i in range(len(X)) if not -1 in X[i,:]]
+X=X[data_indexes,:]
+
+
+
+
+target=dataset[:,0]
+selected_target=target[data_indexes]
+
+#the train and test dataset are selected with the built-in function
+#it is posible to vary the dataset generation process by modyfing the random seed generator
+X_train, X_test, Y_train, Y_test= train_test_split(X, selected_target, train_size=n_samples, random_state=random_seed) 
+
+
+##################################
+#Data values are changed to the logarithmic domain
+
+X_train=np.log(X_train)
+X_test=np.log(X_test)
+Y_train=np.log(Y_train)
+Y_test=np.log(Y_test)
+
+##############################
+
+
+#from the test set 100 samples are selected that are then used for finding the best hyperparamters(gamma and epsilon) for the SVR
+X_validation, X_test, Y_validation, Y_test=train_test_split(X_test, Y_test,train_size=100, random_state=42)
+
+
+#the 100 samples for validation are then divided into two dataset of dimension 50 that are then used for the hyperparameter optimization
+#for each combination of values for gamma and C the SVR is trained using those hyperparameters and are then tested on the X_val_test dataset
+#the hyperparameters combination that achieves the best MSE is chosen
+
+X_val, X_val_test, Y_val, Y_val_test=train_test_split(X_validation, Y_validation, test_size=0.5, random_state=42)
+
+
+
+
+hyperparameter_matrix, best_gamma, best_C=utility.hyperparameters_validation(X_val, Y_val, X_val_test, Y_val_test, gamma_values, C_values, epsilon)
+
+#the hyperparameter matrix stores the MSE values for all the combination of gamma and C parameters
+
+#IMPORTANT: here the dataset used for the solution combination techniques is selected
+#In order to use the same training samples as the classical model the same dataset used for training is also used for the solutions combination techniques 
+#this is done by setting the variable same_validation to True
+#To use a different dataset for the solutions combinations techniques same_validation should be set to False and the dataset should be initialized in the desired way
+#For the scope of this work the same dataset used for training has been used for the solutions combination techniques in order to keep the same number of samples for both classical and quantum implementations
+
+
+if same_validation==True:
+    X_comb=X_train
+    Y_comb=Y_train
+    
+    
+
+
+#%%
+
+#Here the SVR is initialized using the hyperparameters determined in the previous validation
+
+
+svr=SVR(kernel='rbf',C=best_C,epsilon=epsilon,gamma=best_gamma)
+
+svr.fit(X_train, Y_train)
+
+Y_predicted=svr.predict(X_test)
+
+#predicted test samples in the logarithmic domain
+
+Y_predicted_original=np.exp(Y_predicted) 
+
+
+Y_test_original=np.exp(Y_test)
+#to obtain the final values of the prediction values an exponential is applied to the prediction in the logarithmic domain
+
+
+#%%% Quantum part
+
+
+X_train_reshaped=utility.modify_dataset(X_train)
+X_test_reshaped=utility.modify_dataset(X_test)
+X_comb_reshaped=utility.modify_dataset(X_comb)
+
+
+#%%
+#the hyperparameters of the SVR are also used for its quantum implementation
+
+gamma=best_gamma
+C=best_C
+
+N=n_samples
+
+Q=utility.gen_svm_qubos(X_train,Y_train,B, K, xi, best_gamma,epsilon,beta,k0)
+
+Q_max=np.max(Q)
+
+#Q=Q/Q_max
+
+#%%
+alphas,results=utility.dwave_run(Q,B,K,MAXRESULTS,k0)
+
+#the first N alphas refer to the alpha_n varibales, whereas the last N alphas refer to the alpha_hat_n variables
+#plese refer to the description of the implementation of the QSVR for the details of the formulation
+
+alphas_1=alphas[:,0:N]
+alphas_2=alphas[:,N:]
+
+
+alphas_1_mean=np.mean(alphas_1,axis=0)
+alphas_2_mean=np.mean(alphas_2,axis=0)
+
+
+
+#%% solutions combination
+
+#X_comb and Y_comb are defined to be the training and test datsets for the solutions combination techniques
+
+
+#q_vals_scores_m is a matrix that is used to create the solutions combination
+#in q_val_scores_m at the i-th row and the j-th column it is stored the prediction on the j-th
+#element of the X_comb dataset using the i-th set od alpha coefficients returned by the annealer as solution
+
+
+X_comb_original=np.exp(X_comb)
+Y_comb_original=np.exp(Y_comb)
+
+
+#X_val original and Y_val original store the values of the dataset in the orginal domain
+
+
+q_val_scores_m=np.zeros((alphas.shape[0],X_comb.shape[0]))
+
+
+#q_scores_matrix=np.zeros((MAXRESULTS,Y_test.shape[0])) 
+
+
+
+
+alphas_val_norm_1=np.zeros((alphas_1.shape[1]))
+alphas_val_norm_2=np.zeros((alphas_1.shape[1]))
+alphas_val_norm_lc_1=np.zeros((alphas_1.shape[1]))
+alphas_val_softmax_lc_1=np.zeros((alphas_1.shape[1]))
+
+alphas_val_softmax_1=np.zeros((alphas_1.shape[1]))
+alphas_val_softmax_2=np.zeros((alphas_1.shape[1]))
+alphas_val_norm_lc_2=np.zeros((alphas_1.shape[1]))
+alphas_val_softmax_lc_2=np.zeros((alphas_1.shape[1]))
+
+
+
+alphas_best_1=np.zeros((alphas_1.shape[1]))
+alphas_best_2=np.zeros((alphas_2.shape[1]))
+
+
+
+for i in range(alphas.shape[0]):
+    alp_1=alphas_1[i,:]
+    alp_2=alphas_2[i,:]
+    res=utility.predict(X_comb_reshaped,X_train_reshaped,Y_train,alp_1,alp_2,B,K,epsilon,gamma,C)
+    res=np.reshape(res,(X_comb.shape[0],))
+    q_val_scores_m[i,:]=res
+
+
+q_val_scores_original=np.exp(q_val_scores_m)
+
+
+
+
+#the array validation_coefficients store the mean square error calculated for each set of alphas on the validation set
+#the same procedure applies to log_cosh_coefficients but the log-cosh error function is used 
+
+validation_coefficients=np.zeros(alphas.shape[0])
+for i in range(q_val_scores_m.shape[0]):
+    validation_coefficients[i]=mean_squared_error(Y_comb_original,q_val_scores_original[i,:])
+    
+    
+log_cosh_coefficients=np.zeros(alphas.shape[0])
+for i in range(q_val_scores_m.shape[0]):
+    log_cosh_coefficients[i]=utility.log_cosh(Y_comb_original, q_val_scores_original[i,:])
+    
+    
+#annealer_solutions_mse=validation_coefficients
+
+#the score assigned to each solution is then calculated by considering the multiplicative inverse of the validation coefficients
+
+validation_coefficients=1/validation_coefficients
+log_cosh_coefficients=1/log_cosh_coefficients
+
+
+###################################################
+
+#for both mse and log-cosh the coefficients of the weighted average are calculated
+#they both use a normalization where each score is divided by the sum of all the scores or a softmax
+#by the choice of the loss functions and the combination methods they both yield weights
+#that can be used for a weighted average 
+#the variable best_alpha stores the index of the best solution based on the prediction on the validation set
+
+validation_coefficients_softmax=np.exp(validation_coefficients)/sum(np.exp(validation_coefficients))
+validation_coefficients_norm=validation_coefficients/sum(validation_coefficients)
+
+log_cosh_coefficients_n=log_cosh_coefficients/sum(log_cosh_coefficients)
+log_cosh_coefficients_s=np.exp(log_cosh_coefficients)/sum(np.exp(log_cosh_coefficients))
+
+
+
+
+best_alpha=np.argmax(validation_coefficients)
+validation_coefficients_ba=np.zeros((MAXRESULTS,))
+validation_coefficients_ba[best_alpha]=1
+
+#for the method that selects the best alpha all the weights are set to 0 except for the one that stores the best alpha
+
+
+for i in range(alphas_1.shape[1]):
+    alphas_val_norm_1[i]=np.average(alphas_1[:,i], weights=validation_coefficients_norm)
+    alphas_val_norm_2[i]=np.average(alphas_2[:,i], weights=validation_coefficients_norm)
+    
+    
+    alphas_val_softmax_1[i]=np.average(alphas_1[:,i], weights=validation_coefficients_softmax)
+    alphas_val_softmax_2[i]=np.average(alphas_2[:,i], weights=validation_coefficients_softmax)
+    
+    alphas_val_norm_lc_1[i]=np.average(alphas_1[:,i], weights=log_cosh_coefficients_n)
+    alphas_val_norm_lc_2[i]=np.average(alphas_2[:,i], weights=log_cosh_coefficients_n)
+    
+    alphas_val_softmax_lc_1[i]=np.average(alphas_1[:,i], weights=log_cosh_coefficients_s)
+    alphas_val_softmax_lc_2[i]=np.average(alphas_2[:,i], weights=log_cosh_coefficients_s)
+    
+
+    alphas_best_1[i]=np.average(alphas_1[:,i],weights=validation_coefficients_ba)
+    alphas_best_2[i]=np.average(alphas_2[:,i],weights=validation_coefficients_ba)
+
+
+Y_scores_avg=utility.predict(X_test, X_train, Y_train, alphas_1_mean, alphas_2_mean, B, K, epsilon, gamma, C)
+Y_scores_norm=utility.predict(X_test_reshaped,X_train_reshaped,Y_train, alphas_val_norm_1, alphas_val_norm_2, B,K,epsilon,gamma,C)
+Y_scores_softmax=utility.predict(X_test_reshaped,X_train_reshaped,Y_train, alphas_val_softmax_1, alphas_val_softmax_2, B,K,epsilon,gamma,C)
+Y_scores_lc_norm=utility.predict(X_test_reshaped,X_train_reshaped,Y_train, alphas_val_norm_lc_1, alphas_val_norm_lc_2, B,K,epsilon,gamma,C)
+Y_scores_lc_softmax=utility.predict(X_test_reshaped,X_train_reshaped,Y_train, alphas_val_softmax_lc_1, alphas_val_softmax_lc_2, B,K,epsilon,gamma,C)
+Y_scores_best_alpha=utility.predict(X_test_reshaped,X_train_reshaped,Y_train, alphas_best_1, alphas_best_2, B,K,epsilon,gamma,C)
+
+
+#the prediction are converted in the original domain by applying the function np.exp()
+
+
+Y_scores_avg=np.exp(Y_scores_avg)
+Y_scores_norm=np.exp(Y_scores_norm)
+Y_scores_softmax=np.exp(Y_scores_softmax)
+Y_scores_lc_norm=np.exp(Y_scores_lc_norm)
+Y_scores_lc_softmax=np.exp(Y_scores_lc_softmax)
+Y_scores_best_alpha=np.exp(Y_scores_best_alpha)
+
+
+
+
+Y_scores_norm=np.reshape(Y_scores_norm,(Y_scores_norm.shape[0],))
+Y_scores_softmax=np.reshape(Y_scores_softmax,(Y_scores_softmax.shape[0],))
+Y_scores_lc_norm=np.reshape(Y_scores_lc_norm,(Y_scores_lc_norm.shape[0],))
+Y_scores_lc_softmax=np.reshape(Y_scores_lc_softmax,(Y_scores_lc_softmax.shape[0],))
+Y_scores_avg=np.reshape(Y_scores_avg, (Y_scores_avg.shape[0],))
+Y_scores_best_alpha=np.reshape(Y_scores_best_alpha, (Y_scores_best_alpha.shape[0],))
+
+
+#Y_test are converted back to the original domain
+
+
+
+mse_final_results=np.zeros((6,))
+mse_final_results[0]=mean_squared_error(Y_test_original, Y_scores_norm)
+mse_final_results[1]=mean_squared_error(Y_test_original, Y_scores_softmax)
+mse_final_results[2]=mean_squared_error(Y_test_original, Y_scores_lc_norm)
+mse_final_results[3]=mean_squared_error(Y_test_original, Y_scores_lc_softmax)
+mse_final_results[4]=mean_squared_error(Y_test_original, Y_scores_best_alpha)
+mse_final_results[5]=mean_squared_error(Y_test_original, Y_scores_avg)
+
+
+
+best_method=np.argmin(mse_final_results)
+
+methods={0:(Y_scores_norm,'scores norm'),
+         1:(Y_scores_softmax,'scores softmax'),
+         2:(Y_scores_lc_norm,'scores lc norm'),
+         3:(Y_scores_lc_softmax, 'scores lc softmax'),
+         4:(Y_scores_best_alpha,'best set of alphas'),
+         5:(Y_scores_avg,'simple mean')}
+
+
+mse_best_method=mse_final_results[best_method]
+
+
+SVR_mse=mean_squared_error(Y_test_original, Y_predicted_original)
+
+
+methods_array=np.asarray([methods[i][0] for i in range(len(methods))])
+
+compare=np.vstack((methods_array,Y_predicted_original,Y_test_original))
+
+
+print('SVR mse: \n')
+print(SVR_mse)
+print('QSVR mse: \n')
+print(mse_final_results)
+
+
+
+
+now=datetime.now()
+date_time=str(now)
+
+date_time=date_time[0:19]
+
+date_array=[]
+
+
+for i in range(len(date_time)):
+    date_array.append(date_time[i])
+    if date_array[i]==' ' or date_array[i]=='-' or date_array[i]==':':
+        date_array[i]='_'
+
+date_str=''
+for j in date_array:
+    date_str=date_str+j
+
+
+dir_path='Results\\'+date_str
+os.mkdir(dir_path)
+
+
+
+#IMPORTANT
+#the following lines of code saves the datasets that are already in the logarithmic domain
+
+
+np.save(dir_path+'\\X_train.npy',X_train)
+np.save(dir_path+'\\Y_train.npy',Y_train)
+np.save(dir_path+'\\X_test.npy',X_test)
+np.save(dir_path+'\\Y_test.npy',Y_test)
+np.save(dir_path+'\\Q.npy',Q)
+np.save(dir_path+'\\SVR_predictions.npy',Y_predicted_original)
+np.save(dir_path+'\\QSVR_predictions.npy',methods_array)
+
+
+
+
+
+
+
+
+    
+
+
+
+    
+
+
+
+
+
+
+
+
+
+
+
+