diff --git a/doc/src/Regression/eos.py b/doc/src/Regression/eos.py
index 60b30551a..26d694670 100644
--- a/doc/src/Regression/eos.py
+++ b/doc/src/Regression/eos.py
@@ -8,7 +8,7 @@ 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.datasets import load_breast_cancer
+from sklearn.svm import SVR
 
 # Where to save the figures and data files
 PROJECT_ROOT_DIR = "Results"
@@ -33,42 +33,52 @@ def data_path(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("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
-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
-X = np.zeros((len(Density),4))
-X[:,3] = Density**(4.0/3.0)
-X[:,2] = Density
-X[:,1] = Density**(2.0/3.0)
+# 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,random_state=1)
+X_train, X_test, y_train, y_test = train_test_split(X,Energies,test_size=0.2)
 
 
-# We use now Scikit-Learn's linear regressor and ridge regressor
-# OLS part
-clf = skl.LinearRegression().fit(X_train, Energies)
-ytilde = clf.predict(X_test)
-EoS['Eols']  = ytilde
-
+svm = SVR(gamma='auto',C=10.0)
+svm.fit(X_train, y_train)
 
 # The mean squared error                               
-print("Mean squared error: %.2f" % mean_squared_error(Energies, ytilde))
-# Explained variance score: 1 is perfect prediction                                 
-print('Variance score: %.2f' % r2_score(Energies, ytilde))
-# Mean absolute error                                                           
-print('Mean absolute error: %.2f' % mean_absolute_error(Energies, ytilde))
-print(clf.coef_, clf.intercept_)
+print("Test set accuracy: {:.2f}".format(svm.score(X_test,y_test)))
+
 
 
-print("Test set accuracy: {:.2f}".format(clf.score(X_test,y_test)))
 
 from sklearn.preprocessing import MinMaxScaler, StandardScaler
 
@@ -80,9 +90,16 @@ 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 before scaling:\n {}".format(X_train_scaled.min(axis=0)))
-print("Feature max values before scaling:\n {}".format(X_train_scaled.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)))
+
 
 
-#svm.fit(X_train_scaled, y_train)
-#print("Test set accuracy scaled data: {:.2f}".format(svm.score(X_test_scaled,y_test)))