mhjensen
69 lines
2.0 KiB
Python
69 lines
2.0 KiB
Python
# Common imports
|
|
import os
|
|
import numpy as np
|
|
import pandas as pd
|
|
import matplotlib.pyplot as plt
|
|
from sklearn.linear_model import LinearRegression, Ridge, Lasso
|
|
from sklearn.metrics import mean_squared_error
|
|
from sklearn.model_selection import KFold
|
|
from sklearn.model_selection import cross_val_score
|
|
|
|
|
|
# Where to save the figures and data files
|
|
PROJECT_ROOT_DIR = "Results"
|
|
FIGURE_ID = "Results/FigureFiles"
|
|
DATA_ID = "DataFiles/"
|
|
|
|
if not os.path.exists(PROJECT_ROOT_DIR):
|
|
os.mkdir(PROJECT_ROOT_DIR)
|
|
|
|
if not os.path.exists(FIGURE_ID):
|
|
os.makedirs(FIGURE_ID)
|
|
|
|
if not os.path.exists(DATA_ID):
|
|
os.makedirs(DATA_ID)
|
|
|
|
def image_path(fig_id):
|
|
return os.path.join(FIGURE_ID, fig_id)
|
|
|
|
def data_path(dat_id):
|
|
return os.path.join(DATA_ID, dat_id)
|
|
|
|
def save_fig(fig_id):
|
|
plt.savefig(image_path(fig_id) + ".png", format='png')
|
|
|
|
infile = open(data_path("EoS.csv"),'r')
|
|
|
|
# Read the EoS data as csv file and organize the data into two arrays with density and energies
|
|
EoS = pd.read_csv(infile, names=('Density', 'Energy'))
|
|
EoS['Energy'] = pd.to_numeric(EoS['Energy'], errors='coerce')
|
|
EoS = EoS.dropna()
|
|
Energies = EoS['Energy']
|
|
Density = EoS['Density']
|
|
# The design matrix now as function of various polytrops
|
|
|
|
Maxpolydegree = 30
|
|
X = np.zeros((len(Density),Maxpolydegree))
|
|
X[:,0] = 1.0
|
|
estimated_mse_sklearn = np.zeros(Maxpolydegree)
|
|
polynomial = np.zeros(Maxpolydegree)
|
|
k =5
|
|
kfold = KFold(n_splits = k)
|
|
|
|
for polydegree in range(1, Maxpolydegree):
|
|
polynomial[polydegree] = polydegree
|
|
for degree in range(polydegree):
|
|
X[:,degree] = Density**(degree/3.0)
|
|
OLS = LinearRegression()
|
|
# loop over trials in order to estimate the expectation value of the MSE
|
|
estimated_mse_folds = cross_val_score(OLS, X, Energies, scoring='neg_mean_squared_error', cv=kfold)
|
|
#[:, np.newaxis]
|
|
estimated_mse_sklearn[polydegree] = np.mean(-estimated_mse_folds)
|
|
|
|
plt.plot(polynomial, np.log10(estimated_mse_sklearn), label='Test Error')
|
|
plt.xlabel('Polynomial degree')
|
|
plt.ylabel('log10[MSE]')
|
|
plt.legend()
|
|
plt.show()
|
|
|