diff --git a/doc/src/Regression/eos.py b/doc/src/Regression/eos.py
new file mode 100644
index 000000000..60b30551a
--- /dev/null
+++ b/doc/src/Regression/eos.py
@@ -0,0 +1,88 @@
+
+# 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.datasets import load_breast_cancer
+
+# 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
+X = np.zeros((len(Density),4))
+X[:,3] = Density**(4.0/3.0)
+X[:,2] = Density
+X[:,1] = Density**(2.0/3.0)
+X[:,0] = 1
+
+
+X_train, X_test, y_train, y_test = train_test_split(X,Energies,random_state=1)
+
+
+# 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
+
+
+# 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(clf.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 before scaling:\n {}".format(X_train_scaled.min(axis=0)))
+print("Feature max values before scaling:\n {}".format(X_train_scaled.max(axis=0)))
+
+
+#svm.fit(X_train_scaled, y_train)
+#print("Test set accuracy scaled data: {:.2f}".format(svm.score(X_test_scaled,y_test)))