mhjensen
106 lines
3.0 KiB
Python
106 lines
3.0 KiB
Python
|
|
# Common imports
|
|
import os
|
|
import numpy as np
|
|
import pandas as pd
|
|
import matplotlib.pyplot as plt
|
|
import matplotlib.pyplot as plt
|
|
import sklearn.linear_model as skl
|
|
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
|
|
from sklearn.model_selection import train_test_split
|
|
from sklearn.svm import SVR
|
|
|
|
# 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')
|
|
infile = open(data_path("MassEval2016.dat"),'r')
|
|
# Read the EoS data as csv file and organize the data into two arrays with density and energies
|
|
# Read the experimental data with Pandas
|
|
Masses = pd.read_fwf(infile, usecols=(2,3,4,6,11),
|
|
names=('N', 'Z', 'A', 'Element', 'Ebinding'),
|
|
widths=(1,3,5,5,5,1,3,4,1,13,11,11,9,1,2,11,9,1,3,1,12,11,1),
|
|
header=39,
|
|
index_col=False)
|
|
|
|
# Extrapolated values are indicated by '#' in place of the decimal place, so
|
|
# the Ebinding column won't be numeric. Coerce to float and drop these entries.
|
|
Masses['Ebinding'] = pd.to_numeric(Masses['Ebinding'], errors='coerce')
|
|
Masses = Masses.dropna()
|
|
# Convert from keV to MeV.
|
|
Masses['Ebinding'] /= 1000
|
|
|
|
# Group the DataFrame by nucleon number, A.
|
|
Masses = Masses.groupby('A')
|
|
# Find the rows of the grouped DataFrame with the maximum binding energy.
|
|
Masses = Masses.apply(lambda t: t[t.Ebinding==t.Ebinding.max()])
|
|
A = Masses['A']
|
|
Z = Masses['Z']
|
|
N = Masses['N']
|
|
Element = Masses['Element']
|
|
Energies = Masses['Ebinding']
|
|
# Now we set up the design matrix X
|
|
X = np.zeros((len(A),5))
|
|
X[:,0] = 1
|
|
X[:,1] = A
|
|
X[:,2] = A**(2.0/3.0)
|
|
X[:,3] = A**(-1.0/3.0)
|
|
X[:,4] = A**(-1.0)
|
|
|
|
|
|
X_train, X_test, y_train, y_test = train_test_split(X,Energies,test_size=0.2)
|
|
|
|
|
|
svm = SVR(gamma='auto',C=10.0)
|
|
svm.fit(X_train, y_train)
|
|
|
|
# The mean squared error
|
|
print("Test set accuracy: {:.2f}".format(svm.score(X_test,y_test)))
|
|
|
|
|
|
|
|
|
|
from sklearn.preprocessing import MinMaxScaler, StandardScaler
|
|
|
|
scaler = MinMaxScaler()
|
|
scaler.fit(X_train)
|
|
X_train_scaled = scaler.transform(X_train)
|
|
X_test_scaled = scaler.transform(X_test)
|
|
|
|
print("Feature min values before scaling:\n {}".format(X_train.min(axis=0)))
|
|
print("Feature max values before scaling:\n {}".format(X_train.max(axis=0)))
|
|
|
|
print("Feature min values after scaling:\n {}".format(X_train_scaled.min(axis=0)))
|
|
print("Feature max values after scaling:\n {}".format(X_train_scaled.max(axis=0)))
|
|
|
|
|
|
|
|
svm = SVR(gamma='auto',C=10.0)
|
|
svm.fit(X_train_scaled, y_train)
|
|
|
|
|
|
print("Test set accuracy scaled data: {:.2f}".format(svm.score(X_test_scaled,y_test)))
|
|
|
|
|
|
|