SHARE
TWEET

Untitled

a guest Apr 19th, 2019 104 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import time
  2. import numpy as np
  3. import random
  4. import os
  5. from scipy.spatial import distance
  6. from sklearn.model_selection import train_test_split
  7. from skopt import gp_minimize
  8. from scipy.optimize import minimize
  9. from sklearn.metrics import mean_squared_error
  10. from math import sqrt
  11. from scipy.stats.distributions import norm
  12. from sklearn.metrics.pairwise import euclidean_distances
  13. from scipy.spatial import distance
  14. from numpy import inf
  15. from sklearn.preprocessing import StandardScaler
  16. from sklearn.gaussian_process import GaussianProcessRegressor
  17. from sklearn.gaussian_process.kernels import (RBF, Matern, RationalQuadratic,
  18.                                                       ExpSineSquared, DotProduct,
  19.                                                       ConstantKernel,WhiteKernel)
  20.  
  21. from sklearn import metrics
  22. from sklearn.cluster import KMeans                                           
  23. from multiprocessing import Process                                
  24. from multiprocessing import Pool
  25.        
  26. if __name__ == '__main__':
  27.     with Pool(processes=4) as pool:    
  28.                                                      
  29.         np.random.seed(1)
  30.         np.set_printoptions(suppress=True)
  31.  
  32.  
  33.         ################################ Here we load dataset ###################################################################################################
  34.  
  35.         rmtest=np.loadtxt("test.dat",usecols=range(7))                                        # Load the test dataset
  36.         X=rmtest
  37.         kmeans = KMeans(n_clusters=4)
  38.         kmeans.fit(X)
  39.         y_kmeans = kmeans.predict(X)
  40.         labels = kmeans.labels_
  41.  
  42.         a = X[np.where(labels==0)]
  43.         b = X[np.where(labels==1)]
  44.         c = X[np.where(labels==2)]
  45.         d = X[np.where(labels==3)]
  46.  
  47.         print(a.shape)
  48.         print(b.shape)
  49.         print(c.shape)
  50.         print(d.shape)
  51.  
  52.  
  53.         armpredx=np.delete(a,6,1)                                                         # xvalues, in our case inverse distances
  54.         armpredy=np.delete(a,np.s_[0:6],1)                                                # energy values corresponding to the distances
  55.  
  56.         brmpredx=np.delete(b,6,1)                                                         # xvalues, in our case inverse distances
  57.         brmpredy=np.delete(b,np.s_[0:6],1)                                                # energy values corresponding to the distances
  58.  
  59.         crmpredx=np.delete(c,6,1)                                                         # xvalues, in our case inverse distances
  60.         crmpredy=np.delete(c,np.s_[0:6],1)                                                # energy values corresponding to the distances
  61.  
  62.         drmpredx=np.delete(d,6,1)                                                         # xvalues, in our case inverse distances
  63.         drmpredy=np.delete(d,np.s_[0:6],1)                                                # energy values corresponding to the distances
  64.  
  65.  
  66.  
  67.         sample=input("Enter number of training set points\n")                                 # Number of points in our training set. You will be prompted to choose this.
  68.         sample=int(sample)
  69.  
  70.         np.random.seed(24)
  71.  
  72.         apoints=a[np.random.choice(a.shape[0], sample, replace=False), :]            # Training set based on number of points generated here.
  73.         bpoints=b[np.random.choice(b.shape[0], sample, replace=False), :]  
  74.         cpoints=c[np.random.choice(c.shape[0], sample, replace=False), :]  
  75.         dpoints=d[np.random.choice(d.shape[0], sample, replace=False), :]  
  76.            
  77.         print("Total number of sample points for training is ",apoints.shape)
  78.  
  79.  
  80.  
  81.         if os.path.exists('train.dat'):                                                       # we save the training set in train.dat
  82.           os.remove('train.dat')
  83.  
  84.         with open('train.dat', 'a') as l:
  85.             np.savetxt(l, apoints,fmt='%1.6f',delimiter=' ')  
  86.  
  87.          
  88.         ainpx=np.delete(apoints,6,1)                                                              #inpx is the x-part of the training dataset
  89.         ainpy=np.delete(apoints,np.s_[0:6],1)                                                     #inpy is the y-part of the training dataset
  90.  
  91.         binpx=np.delete(bpoints,6,1)                                                              #inpx is the x-part of the training dataset
  92.         binpy=np.delete(bpoints,np.s_[0:6],1)                                                     #inpy is the y-part of the training dataset
  93.  
  94.         cinpx=np.delete(cpoints,6,1)                                                              #inpx is the x-part of the training dataset
  95.         cinpy=np.delete(cpoints,np.s_[0:6],1)                                                     #inpy is the y-part of the training dataset
  96.  
  97.         dinpx=np.delete(dpoints,6,1)                                                              #inpx is the x-part of the training dataset
  98.         dinpy=np.delete(dpoints,np.s_[0:6],1)                                                     #inpy is the y-part of the training dataset
  99.  
  100.         ########################### GPR  done here#############################################################################################
  101.         ## Instantiate a Gaussian Process model
  102.         ##6 types of kernels represented.
  103.  
  104.         k1 = 1**2* RBF(length_scale=0.01)
  105.         k2=11**2*RationalQuadratic(length_scale=0.1, alpha=1.0, length_scale_bounds=(1e-05, 100000.0), alpha_bounds=(1e-05, 100000.0))
  106.         k4=1**2*ExpSineSquared(length_scale=1.0, periodicity=1.0, length_scale_bounds=(1e-05, 100000.0), periodicity_bounds=(1e-05, 100000.0))
  107.         k3 =Matern(length_scale=[0.1, 0.1, 0.1,0.1, 0.1, 0.1],  length_scale_bounds=(1e-5, 1.0),nu=2.5)
  108.         k5=ConstantKernel(constant_value=1.0, constant_value_bounds=(1e-05, 1000.0))
  109.         k6=DotProduct(sigma_0=10, sigma_0_bounds=(1e-05, 100))
  110.         k7=WhiteKernel(noise_level=0.1)
  111.  
  112.         kernel= k5*k3*k6+k7                                                                 # we choose our kernel here.
  113.         kernel1= k5*k3*k6+k7
  114.         kernel2= k5*k3*k6+k7
  115.         kernel3= k5*k3*k6+k7
  116.  
  117.         seconds = time.time()
  118.        
  119.         #optimisation step.
  120.         gp = (GaussianProcessRegressor(kernel=kernel,n_restarts_optimizer=25))
  121.         gp2 = (GaussianProcessRegressor(kernel=kernel1,n_restarts_optimizer=25))
  122.         gp3 = (GaussianProcessRegressor(kernel=kernel2,n_restarts_optimizer=25))
  123.         gp4 = (GaussianProcessRegressor(kernel=kernel3,n_restarts_optimizer=25))
  124.         input_list=[[ainpx,ainpy.reshape(-1,1)],[binpx,binpy.reshape(-1,1)]]
  125.         f_list=[gp.fit,gp2.fit,gp3.fit,gp4.fit]
  126.         pool.starmap(gp.fit, input_list)
  127.         pool.close()
  128.         # agp=pool.apply(gp.fit,args=(ainpx,ainpy.reshape(-1,1)))
  129.          print(os.getpid())
  130.         # bgp=pool.apply(gp2.fit,args=(binpx,binpy.reshape(-1,1)))
  131.         # print(os.getpid())
  132.         # #cgp=pool.apply(gp.fit,args=(cinpx,cinpy.reshape(-1,1)))
  133.         #dgp=pool.apply(gp.fit,args=(dinpx,dinpy.reshape(-1,1)))
  134.        
  135.         # #fitting step.
  136.         # agp=pool.apply(gp.fit(ainpx, ainpy.reshape(-1,1)))                                                    # Here the GPR finds an optimal value of the kernel based on input x,y
  137.         # bgp=pool.apply(gp2.fit(binpx, binpy.reshape(-1,1)))                                                   # Here the GPR finds an optimal value of the kernel based on input x,y
  138.         # cgp=pool.apply(gp3.fit(cinpx, cinpy.reshape(-1,1)))                                                   # Here the GPR finds an optimal value of the kernel based on input x,y
  139.         # dgp=pool.apply(gp4.fit(dinpx, dinpy.reshape(-1,1)))                                                   # Here the GPR finds an optimal value of the kernel based on input x,y
  140.         secs2=time.time()
  141.         print("Fit done")
  142.         print(secs2-seconds)
  143.         print("GPML kernel: %s" % agp.kernel_)                                              # We print the kernel and the likelihood
  144.         print("Log-marginal-likelihood: %.3f"
  145.                 % agp.log_marginal_likelihood(agp.kernel_.theta))
  146.  
  147.         ######################################################### Generate RMSE values ###############################################################################################3
  148.  
  149.         # name=str(len(ainpx))                                                                  #Generating a filename to save our final observations
  150.         # filename = name+'_rmse'
  151.  
  152.         # if os.path.exists(filename):
  153.           # os.remove(filename)
  154.         # f = open(filename , 'a')
  155.  
  156.         # arm_pred, sigma = pool.apply(gp.predict,args=(armpredx, return_std=True))                            #Generating prediction values for our test database
  157.         # brm_pred, sigma1 = pool.apply(gp2.predict,args=(brmpredx, return_std=True))                              #Generating prediction values for our test database
  158.         # crm_pred, sigma2 = pool.apply(gp3.predict,args=(crmpredx, return_std=True))                              #Generating prediction values for our test database
  159.         # drm_pred, sigma3 = pool.apply(gp4.predict,args=(drmpredx, return_std=True))                              #Generating prediction values for our test database
  160.         # brm_pred, sigma1 = gp.predict(brmpredx, return_std=True)                             #Generating prediction values for our test database
  161.         # crm_pred, sigma2 = gp.predict(crmpredx, return_std=True)                             #Generating prediction values for our test database
  162.         # drm_pred, sigma3 = gp.predict(drmpredx, return_std=True)                             #Generating prediction values for our test database
  163.  
  164.  
  165.         # # print("Observed value --- Predicted value---Absolute difference\n")
  166.  
  167.         # armde=np.c_[armpredy, arm_pred,abs(armpredy-arm_pred)]                                  #Generating an array to compare the original and predictedvalues
  168.         # brmde=np.c_[brmpredy, brm_pred,abs(brmpredy-brm_pred)]
  169.         # crmde=np.c_[crmpredy, crm_pred,abs(crmpredy-crm_pred)]
  170.         # drmde=np.c_[drmpredy, drm_pred,abs(drmpredy-drm_pred)]
  171.  
  172.         # # print (rmde)
  173.  
  174.         # aerr=sqrt(mean_squared_error(armpredy, arm_pred))                                     #Generating an RMSE value for our testset
  175.         # print("RMSE-A ",aerr)
  176.  
  177.         # berr=sqrt(mean_squared_error(brmpredy, brm_pred))                                     #Generating an RMSE value for our testset
  178.         # print("RMSE-B ",berr)
  179.  
  180.         # cerr=sqrt(mean_squared_error(crmpredy, crm_pred))                                     #Generating an RMSE value for our testset
  181.         # print("RMSE-C ",cerr)
  182.  
  183.         # derr=sqrt(mean_squared_error(drmpredy, drm_pred))                                     #Generating an RMSE value for our testset
  184.         # print("RMSE-D ",derr)
  185.  
  186.         # # f.write("RMSE value is ")
  187.         # # f.write(str(err))
  188.         # # f.write("\n")                                                                      # Saving the values to the file. We save the RMSE followed by the comparision array
  189.         # # f.write("Observed value --- Predicted value---Absolute difference\n")
  190.         # # np.savetxt(f, rmde,fmt='%1.6f',delimiter=' ')
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top