text: # Regression analysis using scikit-learn functions
  # Common imports
  import numpy as np
  import pandas as pd
  import matplotlib.pyplot as plt
  import os
  from pylab import plt, mpl
  plt.style.use('seaborn')
  mpl.rcParams['font.family'] = 'serif'
  
  def MakePlot(x,y, styles, labels, axlabels):
      plt.figure(figsize=(10,6))
      for i in range(len(x)):
          plt.plot(x[i], y[i], styles[i], label = labels[i])
          plt.xlabel(axlabels[0])
          plt.ylabel(axlabels[1])
      plt.legend(loc=0)
  
  
  
  # 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("MassEval2016.dat"),'r')
  
  
  # 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
  A = Masses['A']
  Z = Masses['Z']
  N = Masses['N']
  Element = Masses['Element']
  Energies = Masses['Ebinding']*A
  
  df = pd.DataFrame({'A':A,'Z':Z, 'N':N,'Element':Element,'Energies':Energies})
  Nucleus = df.loc[lambda df: df.Z==8, :]
  # drop cases with no number
  Nucleus = Nucleus.dropna()
  print(Nucleus)
  Nucleus['NeutronSeparationEnergies'] = Nucleus['Energies'].diff(+1)
  MakePlot([Nucleus.A], [Nucleus.NeutronSeparationEnergies], ['b'], ['Neutron Separation Energy'], ['$A$','$S_n$'])
  save_fig('Nucleus')
  plt.show()
  