from__future__importdivisionimportnumpyasnpfrommatplotlibimportpyplotaspltfromsklearn.gaussian_processimportGaussianProcessRegressorfromsklearn.gaussian_process.kernelsimport(RBF,Matern,RationalQuadratic,ExpSineSquared,DotProduct,ConstantKernel)# ----------------------------------------------------------------------number_of_training_samples=1500number_of_testing_samples=500# read coordinates STANDARDIZEDcoords_training_stand=np.loadtxt('coordinates_training_standardized.txt')coords_testing_stand=np.loadtxt('coordinates_testing_standardized.txt')# read time series TRAIN/TESTtimeseries_training=np.loadtxt('timeseries_training.txt')timeseries_testing=np.loadtxt('timeseries_testing.txt')number_of_time_components=np.shape(timeseries_training)[1]# 20# Instantiate a Gaussian Process modelkernel=1.0*Matern(nu=1.5,length_scale=np.ones(coords_training_stand.shape[1]))gp=GaussianProcessRegressor(kernel=kernel)# placeholder for predictionspred_timeseries_training=np.zeros((np.shape(timeseries_training)))pred_timeseries_testing=np.zeros((np.shape(timeseries_testing)))foriinrange(number_of_time_components):print("time component",i)gp.fit(coords_training_stand,timeseries_training[:,i])y_pred,sigma=gp.predict(coords_training_stand,return_std=True)y_pred_test,sigma_test=gp.predict(coords_testing_stand,return_std=True)pred_timeseries_training[:,i]=y_pred
pred_timeseries_testing[:,i]=y_pred_test# plot trainingfig,ax=plt.subplots(5,figsize=(10,20))foriinrange(5):ax[i].plot(timeseries_training[100*i,:20],color='blue',label='Original train')ax[i].plot(pred_timeseries_training[100*i],color='black',label='GP pred train')ax[i].set_xlabel('Time components',fontsize='x-large')ax[i].set_ylabel('Amplitude',fontsize='x-large')ax[i].set_title('Time series n. {:}'.format(100*i+1),fontsize='x-large')ax[i].legend(fontsize='x-large')plt.subplots_adjust(hspace=1)plt.show()plt.close()# plot testingfig,ax=plt.subplots(5,figsize=(10,20))foriinrange(5):ax[i].plot(timeseries_testing[100*i,:20],color='blue',label='Original test')ax[i].plot(pred_timeseries_testing[100*i],color='black',label='GP pred test')ax[i].set_xlabel('Time components',fontsize='x-large')ax[i].set_ylabel('Amplitude',fontsize='x-large')ax[i].set_title('Time series n. {:}'.format(1500+100*i+1),fontsize='x-large')ax[i].legend(fontsize='x-large')plt.subplots_adjust(hspace=1)plt.show()plt.close()