1.5 MiB
1.5 MiB
In [1]:
print ("Notebook 1: Why Is ML Difficult?")
#This is Python Notebook to walk through polynomial regression examples
#We will use this to think about regression
import numpy as np
import sklearn as sk
from sklearn import datasets, linear_model
from sklearn.preprocessing import PolynomialFeatures
import matplotlib as mpl
from matplotlib import pyplot as plt
#import ml_style as sty # optional style file, can be commented
#mpl.rcParams.update(sty.style) # optional style file, can be commented
%matplotlib notebook
# The Training Data
N_train=100
sigma_train=1;
# Train on integers
x=np.linspace(0.05,0.95,N_train)
# Draw random noise
s = sigma_train*np.random.randn(N_train)
#linear
y=2*x+s
#Tenth Order
#y=2*x-10*x**5+15*x**10+s
p1=plt.plot(x,y, "o",ms=15, label='Training')
#Linear Regression
# Create linear regression object
clf = linear_model.LinearRegression()
# Train the model using the training sets
clf.fit(x[:, np.newaxis], y)
# The coefficients
xplot=np.linspace(0.02,0.98,200)
linear_plot=plt.plot(xplot, clf.predict(xplot[:, np.newaxis]),label='Linear')
#Polynomial Regression
poly3 = PolynomialFeatures(degree=3)
X = poly3.fit_transform(x[:,np.newaxis])
clf3 = linear_model.LinearRegression()
clf3.fit(X,y)
Xplot=poly3.fit_transform(xplot[:,np.newaxis])
poly3_plot=plt.plot(xplot, clf3.predict(Xplot), label='Poly 3')
#poly5 = PolynomialFeatures(degree=5)
#X = poly5.fit_transform(x[:,np.newaxis])
#clf5 = linear_model.LinearRegression()
#clf5.fit(X,y)
#Xplot=poly5.fit_transform(xplot[:,np.newaxis])
#plt.plot(xplot, clf5.predict(Xplot), 'r--',linewidth=1)
poly10 = PolynomialFeatures(degree=10)
X = poly10.fit_transform(x[:,np.newaxis])
clf10 = linear_model.LinearRegression()
clf10.fit(X,y)
Xplot=poly10.fit_transform(xplot[:,np.newaxis])
poly10_plot=plt.plot(xplot, clf10.predict(Xplot), label='Poly 10')
axes = plt.gca()
axes.set_ylim([-7,7])
handles, labels=axes.get_legend_handles_labels()
plt.legend(handles,labels, loc='lower center')
plt.xlabel("$x$")
plt.ylabel("$y$")
Title="$N=$"+str(N_train)+", $\sigma=$"+str(sigma_train)
plt.title(Title+" (train)")
plt.tight_layout()
plt.show()
#Linear Filename
filename_train=Title+"train-linear.pdf"
#Tenth Order Filename
#filename_train=Title+"train-o10.pdf"
plt.savefig(filename_train)
# Generate Test Data
#Number of test data
N_test=20
sigma_test=sigma_train
max_x=1.2
x_test=max_x*np.random.random(N_test)
# Draw random noise
s_test = sigma_test*np.random.randn(N_test)
#Linear
y_test=2*x_test+s_test
#Tenth order
#y_test=2*x_test-10*x_test**5+15*x_test**10+s_test
#Make design matrices for prediction
x_plot=np.linspace(0,max_x, 200)
X3 = poly3.fit_transform(x_plot[:,np.newaxis])
X10 = poly10.fit_transform(x_plot[:,np.newaxis])
%matplotlib notebook
fig = plt.figure()
p1=plt.plot(x_test,y_test.transpose(), 'o', ms=12, label='data')
p2=plt.plot(x_plot,clf.predict(x_plot[:,np.newaxis]), label='linear')
p3=plt.plot(x_plot,clf3.predict(X3), label='3rd order')
p10=plt.plot(x_plot,clf10.predict(X10), label='10th order')
plt.legend(loc=2)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.legend(loc='best')
plt.title(Title+" (pred.)")
plt.tight_layout()
plt.show()
#Linear Filename
filename_test=Title+"pred-linear.pdf"
#Tenth Order Filename
#filename_test=Title+"pred-o10.pdf"
plt.savefig(filename_test)
plt.ylim((-6,12))Out [1]:
Notebook 1: Why Is ML Difficult?
<IPython.core.display.Javascript object>
/usr/local/lib/python3.7/site-packages/sklearn/linear_model/base.py:509: RuntimeWarning: internal gelsd driver lwork query error, required iwork dimension not returned. This is likely the result of LAPACK bug 0038, fixed in LAPACK 3.2.2 (released July 21, 2010). Falling back to 'gelss' driver. linalg.lstsq(X, y)
<IPython.core.display.Javascript object>
(-6, 12)
In [5]:
print ("Notebook 2: Gradient Descent")
#This cell sets up basic plotting functions awe
#we will use to visualize the gradient descent routines.
#Make plots inline
%matplotlib inline
#Make 3D plots
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
#from matplotlib import animation
#from IPython.display import HTML
from matplotlib.colors import LogNorm
#from itertools import zip_longest
#Import Numpy
import numpy as np
#Define function for plotting
def plot_surface(x,y,z,azim=-60,elev=40, dist=10, cmap="RdYlBu_r"):
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
plot_args = {'rstride': 1, 'cstride': 1, 'cmap':cmap,
'linewidth': 20, 'antialiased': True,
'vmin': -2, 'vmax': 2}
ax.plot_surface(x, y, z, **plot_args)
ax.view_init(azim=azim, elev=elev)
ax.dist=dist
ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
ax.set_zlim(-2, 2)
plt.xticks([-1, -0.5, 0, 0.5, 1],
[r"$-1$", r"$-1/2$", r"$0$", r"$1/2$", r"$1$"])
plt.yticks([-1, -0.5, 0, 0.5, 1],
[r"$-1$", r"$-1/2$", r"$0$", r"$1/2$", r"$1$"])
ax.set_zticks([-2, -1, 0, 1, 2])
ax.set_zticklabels([r"$-2$", r"$-1$", r"$0$", r"$1$", r"$2$"])
ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.set_xlabel(r"$x$", fontsize=18)
ax.set_ylabel(r"$y$", fontsize=18)
ax.set_zlabel(r"z", fontsize=18)
return fig,ax;
def overlay_trajectory_quiver(ax,obj_func,trajectory, color='k'):
xs=trajectory[:,0]
ys=trajectory[:,1]
zs=obj_func(xs,ys)
ax.quiver(xs[:-1], ys[:-1], zs[:-1], xs[1:]-xs[:-1], ys[1:]-ys[:-1],zs[1:]-zs[:-1],color=color,arrow_length_ratio=0.3)
return ax;
def overlay_trajectory(ax,obj_func,trajectory,label,color='k'):
xs=trajectory[:,0]
ys=trajectory[:,1]
zs=obj_func(xs,ys)
ax.plot(xs,ys,zs, color, label=label)
return ax;
def overlay_trajectory_contour_M(ax,trajectory, label,color='k',lw=2):
xs=trajectory[:,0]
ys=trajectory[:,1]
ax.plot(xs,ys, color, label=label,lw=lw)
ax.plot(xs[-1],ys[-1],color+'>', markersize=14)
return ax;
def overlay_trajectory_contour(ax,trajectory, label,color='k',lw=2):
xs=trajectory[:,0]
ys=trajectory[:,1]
ax.plot(xs,ys, color, label=label,lw=lw)
return ax;
#DEFINE SURFACES WE WILL WORK WITH
#Define monkey saddle and gradient
def monkey_saddle(x,y):
return x**3 - 3*x*y**2
def grad_monkey_saddle(params):
x=params[0]
y=params[1]
grad_x= 3*x**2-3*y**2
grad_y= -6*x*y
return [grad_x,grad_y]
#Define saddle surface
def saddle_surface(x,y,a=1,b=1):
return a*x**2-b*y**2
def grad_saddle_surface(params,a=1,b=1):
x=params[0]
y=params[1]
grad_x= a*x
grad_y= -1*b*y
return [grad_x,grad_y]
# Define minima_surface
def minima_surface(x,y,a=1,b=1):
return a*x**2+b*y**2-1
def grad_minima_surface(params,a=1,b=1):
x=params[0]
y=params[1]
grad_x= 2*a*x
grad_y= 2*b*y
return [grad_x,grad_y]
def beales_function(x,y):
f=np.square(1.5-x+x*y)+np.square(2.25-x+x*y*y)+np.square(2.625-x+x*y**3)
return f
def grad_beales_function(params):
x=params[0]
y=params[1]
grad_x=2*(1.5-x+x*y)*(-1+y)+2*(2.25-x+x*y**2)*(-1+y**2)+2*(2.625-x+x*y**3)*(-1+y**3)
grad_y=2*(1.5-x+x*y)*x+4*(2.25-x+x*y**2)*x*y+6*(2.625-x+x*y**3)*x*y**2
return [grad_x,grad_y]
def contour_beales_function():
#plot beales function
x, y = np.meshgrid(np.arange(-4.5, 4.5, 0.2), np.arange(-4.5, 4.5, 0.2))
fig, ax = plt.subplots(figsize=(10, 6))
z=beales_function(x,y)
ax.contour(x, y, z, levels=np.logspace(0, 5, 35), norm=LogNorm(), cmap="RdYlBu_r")
ax.plot(3,0.5, 'r*', markersize=18)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_xlim((-4.5, 4.5))
ax.set_ylim((-4.5, 4.5))
return fig,ax
#Make plots of surfaces
x, y = np.mgrid[-1:1:31j, -1:1:31j]
fig1,ax1=plot_surface(x,y,monkey_saddle(x,y))
fig2,ax2=plot_surface(x,y,saddle_surface(x,y))
fig3,ax3=plot_surface(x,y,minima_surface(x,y,5),0)
#Contour plot of Beale's Function
fig4,ax4 =contour_beales_function()
plt.show()
#This writes a simple gradient descent, gradient descent+ momentum,
#nesterov.
#Mean-gradient based methods
def gd(grad, init, n_epochs=1000, eta=10**-4, noise_strength=0):
#This is a simple optimizer
params=np.array(init)
param_traj=np.zeros([n_epochs+1,2])
param_traj[0,]=init
v=0;
for j in range(n_epochs):
noise=noise_strength*np.random.randn(params.size)
v=eta*(np.array(grad(params))+noise)
params=params-v
param_traj[j+1,]=params
return param_traj
def gd_with_mom(grad, init, n_epochs=5000, eta=10**-4, gamma=0.9,noise_strength=0):
params=np.array(init)
param_traj=np.zeros([n_epochs+1,2])
param_traj[0,]=init
v=0
for j in range(n_epochs):
noise=noise_strength*np.random.randn(params.size)
v=gamma*v+eta*(np.array(grad(params))+noise)
params=params-v
param_traj[j+1,]=params
return param_traj
def NAG(grad, init, n_epochs=5000, eta=10**-4, gamma=0.9,noise_strength=0):
params=np.array(init)
param_traj=np.zeros([n_epochs+1,2])
param_traj[0,]=init
v=0
for j in range(n_epochs):
noise=noise_strength*np.random.randn(params.size)
params_nesterov=params-gamma*v
v=gamma*v+eta*(np.array(grad(params_nesterov))+noise)
params=params-v
param_traj[j+1,]=params
return param_traj
# Investigate effect of learning rate in GD
x, y = np.meshgrid(np.arange(-4.5, 4.5, 0.2), np.arange(-4.5, 4.5, 0.2))
fig, ax = plt.subplots(figsize=(10, 6))
z=minima_surface(x,y,1,10)
ax.contour(x, y, z, levels=np.logspace(0, 5, 35), norm=LogNorm(), cmap="RdYlBu_r")
ax.plot(0,0, 'r*', markersize=18)
#initial point
init1=[-2,4]
init2=[-1.7,4]
init3=[-1.5,4]
init4=[-3,4.5]
eta1=0.1
eta2=0.5
eta3=1
eta4=1.01
gd_1=gd(grad_minima_surface,init1, n_epochs=100, eta=eta1)
gd_2=gd(grad_minima_surface,init2, n_epochs=100, eta=eta2)
gd_3=gd(grad_minima_surface,init3, n_epochs=100, eta=eta3)
gd_4=gd(grad_minima_surface,init4, n_epochs=10, eta=eta4)
#print(gd_1)
overlay_trajectory_contour(ax,gd_1,'$\eta=$%s'% eta1,'g--*', lw=0.5)
overlay_trajectory_contour(ax,gd_2,'$\eta=$%s'% eta2,'b-<', lw=0.5)
overlay_trajectory_contour(ax,gd_3,'$\eta=$%s'% eta3,'->', lw=0.5)
overlay_trajectory_contour(ax,gd_4,'$\eta=$%s'% eta4,'c-o', lw=0.5)
plt.legend(loc=2)
plt.show()
fig.savefig("GD3regimes.pdf", bbox_inches='tight')
#Methods that exploit first and second moments of gradient: RMS-PROP and ADAMS
def rms_prop(grad, init, n_epochs=5000, eta=10**-3, beta=0.9,epsilon=10**-8,noise_strength=0):
params=np.array(init)
param_traj=np.zeros([n_epochs+1,2])
param_traj[0,]=init#Import relevant packages
grad_sq=0;
for j in range(n_epochs):
noise=noise_strength*np.random.randn(params.size)
g=np.array(grad(params))+noise
grad_sq=beta*grad_sq+(1-beta)*g*g
v=eta*np.divide(g,np.sqrt(grad_sq+epsilon))
params= params-v
param_traj[j+1,]=params
return param_traj
def adams(grad, init, n_epochs=5000, eta=10**-4, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0):
params=np.array(init)
param_traj=np.zeros([n_epochs+1,2])
param_traj[0,]=init
v=0;
grad_sq=0;
for j in range(n_epochs):
noise=noise_strength*np.random.randn(params.size)
g=np.array(grad(params))+noise
v=gamma*v+(1-gamma)*g
grad_sq=beta*grad_sq+(1-beta)*g*g
v_hat=v/(1-gamma)
grad_sq_hat=grad_sq/(1-beta)
params=params-eta*np.divide(v_hat,np.sqrt(grad_sq_hat+epsilon))
param_traj[j+1,]=params
return param_traj
#Make static plot of the results
Nsteps=10**4
lr_l=10**-3
lr_s=10**-6
init1=np.array([4,3])
fig1, ax1=contour_beales_function()
gd_trajectory1=gd(grad_beales_function,init1,Nsteps, eta=lr_s, noise_strength=0)
gdm_trajectory1=gd_with_mom(grad_beales_function,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
NAG_trajectory1=NAG(grad_beales_function,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory1=rms_prop(grad_beales_function,init1,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=0)
adam_trajectory1=adams(grad_beales_function,init1,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)
overlay_trajectory_contour_M(ax1,gd_trajectory1, 'GD','k')
overlay_trajectory_contour_M(ax1,gd_trajectory1, 'GDM','m')
overlay_trajectory_contour_M(ax1,NAG_trajectory1, 'NAG','c--')
overlay_trajectory_contour_M(ax1,rms_prop_trajectory1,'RMS', 'b-.')
overlay_trajectory_contour_M(ax1,adam_trajectory1,'ADAMS', 'r')
plt.legend(loc=2)
#init2=np.array([1.5,1.5])
#gd_trajectory2=gd(grad_beales_function,init2,Nsteps, eta=10**-6, noise_strength=0)
#gdm_trajectory2=gd_with_mom(grad_beales_function,init2,Nsteps,eta=10**-6, gamma=0.9,noise_strength=0)
#NAG_trajectory2=NAG(grad_beales_function,init2,Nsteps,eta=10**-6, gamma=0.9,noise_strength=0)
#rms_prop_trajectory2=rms_prop(grad_beales_function,init2,Nsteps,eta=10**-3, beta=0.9,epsilon=10**-8,noise_strength=0)
#adam_trajectory2=adams(grad_beales_function,init2,Nsteps,eta=10**-3, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)
#overlay_trajectory_contour_M(ax1,gdm_trajectory2, 'GDM','m')
#overlay_trajectory_contour_M(ax1,NAG_trajectory2, 'NAG','c--')
#overlay_trajectory_contour_M(ax1,rms_prop_trajectory2,'RMS', 'b-.')
#overlay_trajectory_contour_M(ax1,adam_trajectory2,'ADAMS', 'r')
init3=np.array([-1,4])
gd_trajectory3=gd(grad_beales_function,init3,10**5, eta=lr_s, noise_strength=0)
gdm_trajectory3=gd_with_mom(grad_beales_function,init3,10**5,eta=lr_s, gamma=0.9,noise_strength=0)
NAG_trajectory3=NAG(grad_beales_function,init3,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory3=rms_prop(grad_beales_function,init3,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=0)
adam_trajectory3=adams(grad_beales_function,init3,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)
overlay_trajectory_contour_M(ax1,gd_trajectory3, 'GD','k')
overlay_trajectory_contour_M(ax1,gdm_trajectory3, 'GDM','m')
overlay_trajectory_contour_M(ax1,NAG_trajectory3, 'NAG','c--')
overlay_trajectory_contour_M(ax1,rms_prop_trajectory3,'RMS', 'b-.')
overlay_trajectory_contour_M(ax1,adam_trajectory3,'ADAMS', 'r')
init4=np.array([-2,-4])
gd_trajectory4=gd(grad_beales_function,init4,Nsteps, eta=lr_s, noise_strength=0)
gdm_trajectory4=gd_with_mom(grad_beales_function,init4,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
NAG_trajectory4=NAG(grad_beales_function,init4,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory4=rms_prop(grad_beales_function,init4,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=0)
adam_trajectory4=adams(grad_beales_function,init4,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)
overlay_trajectory_contour_M(ax1,gd_trajectory4, 'GD','k')
overlay_trajectory_contour_M(ax1,gdm_trajectory4, 'GDM','m')
overlay_trajectory_contour_M(ax1,NAG_trajectory4, 'NAG','c--')
overlay_trajectory_contour_M(ax1,rms_prop_trajectory4,'RMS', 'b-.')
overlay_trajectory_contour_M(ax1,adam_trajectory4,'ADAMS', 'r')
plt.show()Notebook 2: Gradient Descent
/anaconda3/lib/python3.6/site-packages/matplotlib/contour.py:1537: UserWarning: Log scale: values of z <= 0 have been masked
warnings.warn('Log scale: values of z <= 0 have been masked')
In [1]:
print ("Section 3: Linear Regression")
%matplotlib inline
#This code is modified from plot_cv_diabetes.py in the skit-learn documentation
#and plot_ridge_path.py
from __future__ import print_function
print(__doc__)
import numpy as np
import matplotlib.pyplot as plt
import seaborn
from sklearn import datasets, linear_model
#Load Training Data set with 200 examples
number_examples=200
diabetes = datasets.load_diabetes()
X = diabetes.data[:number_examples]
y = diabetes.target[:number_examples]
#Set up Lasso and Ridge Regression models
ridge=linear_model.Ridge()
lasso = linear_model.Lasso()
#Chooose paths
alphas = np.logspace(-2, 2, 10)
# To see how well we learn, we partition the dataset into a training set with 150
# as well as a test set with 50 examples. We record their errors respectively.
n_samples = 150
n_samples_train = 100
X_train, X_test = X[:n_samples_train], X[n_samples_train:]
y_train, y_test = y[:n_samples_train], y[n_samples_train:]
train_errors_ridge = list()
test_errors_ridge = list()
train_errors_lasso = list()
test_errors_lasso = list()
#Initialize coeffficients for ridge regression and Lasso
coefs_ridge = []
coefs_lasso=[]
for a in alphas:
ridge.set_params(alpha=a)
ridge.fit(X_train, y_train)
coefs_ridge.append(ridge.coef_)
# Use the coefficient of determination R^2 as the performance of prediction.
train_errors_ridge.append(ridge.score(X_train, y_train))
test_errors_ridge.append(ridge.score(X_test, y_test))
lasso.set_params(alpha=a)
lasso.fit(X_train, y_train)
coefs_lasso.append(lasso.coef_)
train_errors_lasso.append(lasso.score(X_train, y_train))
test_errors_lasso.append(lasso.score(X_test, y_test))
###############################################################################
# Display results
# First see how the 10 features we learned scale as we change the regularization parameter
plt.subplot(1,2,1)
plt.semilogx(alphas, np.abs(coefs_ridge))
axes = plt.gca()
#ax.set_xscale('log')
#ax.set_xlim(ax.get_xlim()[::-1]) # reverse axis
plt.xlabel(r'$\lambda$',fontsize=18)
plt.ylabel('$|w_i|$',fontsize=18)
plt.title('Ridge')
#plt.savefig("Ridge_sparsity_scale.pdf.pdf")
plt.subplot(1,2,2)
plt.semilogx(alphas, np.abs(coefs_lasso))
axes = plt.gca()
#ax.set_xscale('log')
#ax.set_xlim(ax.get_xlim()[::-1]) # reverse axis
plt.xlabel(r'$\lambda$',fontsize=18)
#plt.ylabel('$|\mathbf{w}|$',fontsize=18)
plt.title('LASSO')
#plt.savefig("LASSO_sparsity_scale.pdf")
plt.show()
# Plot our performance on both the training and test data
plt.semilogx(alphas, train_errors_ridge, 'b',label='Train (Ridge)')
plt.semilogx(alphas, test_errors_ridge, '--b',label='Test (Ridge)')
plt.semilogx(alphas, train_errors_lasso, 'g',label='Train (LASSO)')
plt.semilogx(alphas, test_errors_lasso, '--g',label='Test (LASSO)')
#plt.vlines(alpha_optim, plt.ylim()[0], np.max(test_errors), color='k',
# linewidth=3, label='Optimum on test')
plt.legend(loc='upper right')
plt.ylim([0, 1.0])
plt.xlabel(r'$\lambda$',fontsize=18)
plt.ylabel('Performance')
#plt.savefig("Ridge_LASSO_sparsity_performance.pdf")
plt.show()Section 3: Linear Regression Automatically created module for IPython interactive environment
In [1]:
print ("Section 4: Hamiltonian")
import numpy as np
import scipy.sparse as sp
np.random.seed(12)
import warnings
#Comment this to turn on warnings
warnings.filterwarnings('ignore')
### define Ising model aprams
# system size
L=40
# create 10000 random Ising states
states=np.random.choice([-1, 1], size=(10000,L))
def ising_energies(states,L):
"""
This function calculates the energies of the states in the nn Ising Hamiltonian
"""
J=np.zeros((L,L),)
for i in range(L):
J[i,(i+1)%L]-=1.0
# compute energies
E = np.einsum('...i,ij,...j->...',states,J,states)
return E
# calculate Ising energies
energies=ising_energies(states,L)
# reshape Ising states into RL samples: S_iS_j --> X_p
states=np.einsum('...i,...j->...ij', states, states)
shape=states.shape
states=states.reshape((shape[0],shape[1]*shape[2]))
# build final data set
Data=[states,energies]
# define number of samples
n_samples=400
# define train and test data sets
X_train=Data[0][:n_samples]
Y_train=Data[1][:n_samples] #+ np.random.normal(0,4.0,size=X_train.shape[0])
X_test=Data[0][n_samples:3*n_samples//2]
Y_test=Data[1][n_samples:3*n_samples//2] #+ np.random.normal(0,4.0,size=X_test.shape[0])
from sklearn import linear_model
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import seaborn
%matplotlib inline
# set up Lasso and Ridge Regression models
leastsq=linear_model.LinearRegression()
ridge=linear_model.Ridge()
lasso = linear_model.Lasso()
# define error lists
train_errors_leastsq = []
test_errors_leastsq = []
train_errors_ridge = []
test_errors_ridge = []
train_errors_lasso = []
test_errors_lasso = []
# set refularisations trength values
lmbdas = np.logspace(-4, 5, 10)
#Initialize coeffficients for ridge regression and Lasso
coefs_leastsq = []
coefs_ridge = []
coefs_lasso=[]
for lmbda in lmbdas:
### ordinary least squares
leastsq.fit(X_train, Y_train) # fit model
coefs_leastsq.append(leastsq.coef_) # store weights
# use the coefficient of determination R^2 as the performance of prediction.
train_errors_leastsq.append(leastsq.score(X_train, Y_train))
test_errors_leastsq.append(leastsq.score(X_test,Y_test))
### apply Ridge regression
ridge.set_params(alpha=lmbda) # set regularisation parameter
ridge.fit(X_train, Y_train) # fit model
coefs_ridge.append(ridge.coef_) # store weights
# use the coefficient of determination R^2 as the performance of prediction.
train_errors_ridge.append(ridge.score(X_train, Y_train))
test_errors_ridge.append(ridge.score(X_test,Y_test))
### apply Ridge regression
lasso.set_params(alpha=lmbda) # set regularisation parameter
lasso.fit(X_train, Y_train) # fit model
coefs_lasso.append(lasso.coef_) # store weights
# use the coefficient of determination R^2 as the performance of prediction.
train_errors_lasso.append(lasso.score(X_train, Y_train))
test_errors_lasso.append(lasso.score(X_test,Y_test))
### plot Ising interaction J
J_leastsq=np.array(leastsq.coef_).reshape((L,L))
J_ridge=np.array(ridge.coef_).reshape((L,L))
J_lasso=np.array(lasso.coef_).reshape((L,L))
cmap_args=dict(vmin=-1., vmax=1., cmap='seismic')
fig, axarr = plt.subplots(nrows=1, ncols=3)
axarr[0].imshow(J_leastsq,**cmap_args)
axarr[0].set_title('$\\mathrm{OLS}$',fontsize=16)
axarr[0].tick_params(labelsize=16)
axarr[1].imshow(J_ridge,**cmap_args)
axarr[1].set_title('$\\mathrm{Ridge},\ \\lambda=%.4f$' %(lmbda),fontsize=16)
axarr[1].tick_params(labelsize=16)
im=axarr[2].imshow(J_lasso,**cmap_args)
axarr[2].set_title('$\\mathrm{LASSO},\ \\lambda=%.4f$' %(lmbda),fontsize=16)
axarr[2].tick_params(labelsize=16)
divider = make_axes_locatable(axarr[2])
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar=fig.colorbar(im, cax=cax)
cbar.ax.set_yticklabels(np.arange(-1.0, 1.0+0.25, 0.25),fontsize=14)
cbar.set_label('$J_{i,j}$',labelpad=-40, y=1.12,fontsize=16,rotation=0)
fig.subplots_adjust(right=2.0)
plt.show()
# Plot our performance on both the training and test data
plt.semilogx(lmbdas, train_errors_leastsq, 'b',label='Train (OLS)')
plt.semilogx(lmbdas, test_errors_leastsq,'--b',label='Test (OLS)')
plt.semilogx(lmbdas, train_errors_ridge,'r',label='Train (Ridge)',linewidth=1)
plt.semilogx(lmbdas, test_errors_ridge,'--r',label='Test (Ridge)',linewidth=1)
plt.semilogx(lmbdas, train_errors_lasso, 'g',label='Train (LASSO)')
plt.semilogx(lmbdas, test_errors_lasso, '--g',label='Test (LASSO)')
fig = plt.gcf()
fig.set_size_inches(10.0, 6.0)
#plt.vlines(alpha_optim, plt.ylim()[0], np.max(test_errors), color='k',
# linewidth=3, label='Optimum on test')
plt.legend(loc='lower left',fontsize=16)
plt.ylim([-0.01, 1.01])
plt.xlim([min(lmbdas), max(lmbdas)])
plt.xlabel(r'$\lambda$',fontsize=16)
plt.ylabel('Performance',fontsize=16)
plt.tick_params(labelsize=16)
plt.show()Section 4: Hamiltonian
In [3]:
print ("notebook 6: Logistic Regression, Ising")
import numpy as np
import warnings
#Comment this to turn on warnings
warnings.filterwarnings('ignore')
np.random.seed() # shuffle random seed generator
# Ising model parameters
L=40 # linear system size
J=-1.0 # Ising interaction
T=np.linspace(0.25,4.0,16) # set of temperatures
T_c=2.26 # Onsager critical temperature in the TD limit
##### prepare training and test data sets
import pickle,os
from sklearn.model_selection import train_test_split
###### define ML parameters
num_classes=2
train_to_test_ratio=0.5 # training samples
# path to data directory
path_to_data=os.path.expanduser('~')+'/Dropbox/MachineLearningReview/Datasets/isingMC/'
# load data
file_name = "Ising2DFM_reSample_L40_T=All.pkl" # this file contains 16*10000 samples taken in T=np.arange(0.25,4.0001,0.25)
data = pickle.load(open(path_to_data+file_name,'rb')) # pickle reads the file and returns the Python object (1D array, compressed bits)
data = np.unpackbits(data).reshape(-1, 1600) # Decompress array and reshape for convenience
data=data.astype('int')
data[np.where(data==0)]=-1 # map 0 state to -1 (Ising variable can take values +/-1)
file_name = "Ising2DFM_reSample_L40_T=All_labels.pkl" # this file contains 16*10000 samples taken in T=np.arange(0.25,4.0001,0.25)
labels = pickle.load(open(path_to_data+file_name,'rb')) # pickle reads the file and returns the Python object (here just a 1D array with the binary labels)
# divide data into ordered, critical and disordered
X_ordered=data[:70000,:]
Y_ordered=labels[:70000]
X_critical=data[70000:100000,:]
Y_critical=labels[70000:100000]
X_disordered=data[100000:,:]
Y_disordered=labels[100000:]
del data,labels
# define training and test data sets
X=np.concatenate((X_ordered,X_disordered))
Y=np.concatenate((Y_ordered,Y_disordered))
# pick random data points from ordered and disordered states
# to create the training and test sets
X_train,X_test,Y_train,Y_test=train_test_split(X,Y,train_size=train_to_test_ratio)
# full data set
X=np.concatenate((X_critical,X))
Y=np.concatenate((Y_critical,Y))
print('X_train shape:', X_train.shape)
print('Y_train shape:', Y_train.shape)
print()
print(X_train.shape[0], 'train samples')
print(X_critical.shape[0], 'critical samples')
print(X_test.shape[0], 'test samples')notebook 6: Logistic Regression, Ising
[0;31m---------------------------------------------------------------------------[0m [0;31mFileNotFoundError[0m Traceback (most recent call last) [0;32m<ipython-input-3-fd8d2f3c3b85>[0m in [0;36m<module>[0;34m()[0m [1;32m 28[0m [0;31m# load data[0m[0;34m[0m[0;34m[0m[0m [1;32m 29[0m [0mfile_name[0m [0;34m=[0m [0;34m"Ising2DFM_reSample_L40_T=All.pkl"[0m [0;31m# this file contains 16*10000 samples taken in T=np.arange(0.25,4.0001,0.25)[0m[0;34m[0m[0m [0;32m---> 30[0;31m [0mdata[0m [0;34m=[0m [0mpickle[0m[0;34m.[0m[0mload[0m[0;34m([0m[0mopen[0m[0;34m([0m[0mpath_to_data[0m[0;34m+[0m[0mfile_name[0m[0;34m,[0m[0;34m'rb'[0m[0;34m)[0m[0;34m)[0m [0;31m# pickle reads the file and returns the Python object (1D array, compressed bits)[0m[0;34m[0m[0m [0m[1;32m 31[0m [0mdata[0m [0;34m=[0m [0mnp[0m[0;34m.[0m[0munpackbits[0m[0;34m([0m[0mdata[0m[0;34m)[0m[0;34m.[0m[0mreshape[0m[0;34m([0m[0;34m-[0m[0;36m1[0m[0;34m,[0m [0;36m1600[0m[0;34m)[0m [0;31m# Decompress array and reshape for convenience[0m[0;34m[0m[0m [1;32m 32[0m [0mdata[0m[0;34m=[0m[0mdata[0m[0;34m.[0m[0mastype[0m[0;34m([0m[0;34m'int'[0m[0;34m)[0m[0;34m[0m[0m [0;31mFileNotFoundError[0m: [Errno 2] No such file or directory: '/Users/dylansmith/Dropbox/MachineLearningReview/Datasets/isingMC/Ising2DFM_reSample_L40_T=All.pkl'
In [ ]:
print ("notebook 5: Logistic Regression, SUSY")
# Importing the SUSY Data set
import sys, os
import pandas as pd
import numpy as np
import warnings
#Commnet this to turn on warnings
warnings.filterwarnings('ignore')
seed=12
np.random.seed(seed)
import tensorflow as tf
# suppress tflow compilation warnings
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'
tf.set_random_seed(seed)
columns=["signal", "lepton 1 pT", "lepton 1 eta", "lepton 1 phi", "lepton 2 pT", "lepton 2 eta",
"lepton 2 phi", "missing energy magnitude", "missing energy phi", "MET_rel",
"axial MET", "M_R", "M_TR_2", "R", "MT2", "S_R", "M_Delta_R", "dPhi_r_b", "cos(theta_r1)"]
#Load 1,500,000 rows as train data, 50,000 as test data
filename="/Users/dylansmith/Downloads/SUSY.csv"
df_train=pd.read_csv(filename,names=columns,nrows=1500000,engine='python')
df_test=pd.read_csv(filename,names=columns,nrows=50000, skiprows=1500000,engine='python')
%matplotlib inline
#import ml_style as style
import matplotlib as mpl
import matplotlib.pyplot as plt
#mpl.rcParams.update(style.style)
def getTrainData(nVar):
ExamplesTrain = df_train.iloc[:,1:nVar+1].as_matrix()
#now the signal
ResultsTrain = df_train.iloc[:,0:1].as_matrix()
return (ExamplesTrain,ResultsTrain)
def getTestData(nVar):
ExamplesTest = df_test.iloc[:,1:nVar+1].as_matrix()
#now the signal
ResultsTest = df_test.iloc[:,0:1].as_matrix()
return (ExamplesTest,ResultsTest)
#let's define this as a function so we can call it easily
def runTensorFlowRegression(nVar,alpha):
#make data array placeholder for just first 8 simple features
x = tf.placeholder(tf.float32,[None,nVar])
#make weights and bias
W = tf.Variable(tf.zeros([nVar,2])) #we will make y 'onehot' 0 bit is bkg, 1 bit is signal
b = tf.Variable(tf.zeros([2]))
#make 'answer variable'
y = tf.nn.softmax(tf.matmul(x, W) + b)
#placeholder for correct answer
y_ = tf.placeholder(tf.float32, [None, 2])
#cross entropy
cross_entropy = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits(logits=y,labels=y_)+alpha*tf.nn.l2_loss(W))
#define training step
train_step = tf.train.GradientDescentOptimizer(0.5).minimize(cross_entropy)
#initialize variables
init = tf.global_variables_initializer()
#setup session
sess = tf.Session()
sess.run(init)
#ok now everything is setup for tensorflow, but we need the data in a useful form
#first let's get the variables
Var_train,Sig_train_bit1 = getTrainData(nVar)
#now the signal
Sig_train_bit0 = Sig_train_bit1.copy()
Sig_train_bit0 = 1 - Sig_train_bit0
Sig_train = np.column_stack((Sig_train_bit0,Sig_train_bit1))
#now run with batches of 100 data points
for i in range(0,15000):
start = i*100
end = (i+1)*100-1
batch_x = Var_train[start:end]
batch_y = Sig_train[start:end]
sess.run(train_step, feed_dict={x: batch_x, y_: batch_y})
#now test accuracy
correct_prediction = tf.equal(tf.argmax(y,1), tf.argmax(y_,1))
accuracy = tf.reduce_mean(tf.cast(correct_prediction, tf.float32))
#setup test data
Var_test = df_test.iloc[:,1:nVar+1].as_matrix()
#now the signal
Sig_test_bit1 = df_test.iloc[:,0:1].as_matrix()
Sig_test_bit0 = Sig_test_bit1.copy()
Sig_test_bit0 = 1 - Sig_test_bit0
Sig_test = np.column_stack((Sig_test_bit0,Sig_test_bit1))
print("Accuracy for alpha %.1E : %.3f" %(alpha,sess.run(accuracy, feed_dict={x: Var_test, y_: Sig_test})))
#get the weights
weights = W.eval(session=sess)
#get probabilities assigned (i.e. evaluate y on test data)
probs = y.eval(feed_dict = {x: Var_test}, session = sess)
#print probs
#now let's get the signal efficiency and background rejection on the test data
Acceptance = []
Rejection = []
theshes = np.arange(0,1,0.01)
for thresh in theshes:
it=0
nBkg = 0.0
nBkgCor = 0.0
nSig = 0.0
nSigCor = 0.0
for prob in probs:
if prob[0] > thresh:
sig=False
else:
sig=True
if Sig_test_bit1[it][0]: #actual signal
nSig+=1
if sig:
nSigCor+=1
else: #actual background
nBkg+=1
if not sig:
nBkgCor+=1
it+=1
#now append signal efficiency and bakcground rejection
Acceptance.append(nSigCor/nSig)
Rejection.append(nBkgCor/nBkg)
return (probs,Acceptance,Rejection)
alphas = np.logspace(-10,1,11)
fig = plt.figure()
ax = fig.add_subplot(111)
it=0
for alpha in alphas:
c1 = 1.*( float(it) % 3.)/3.0
c2 = 1.*( float(it) % 9.)/9.0
c3 = 1.*( float(it) % 27.)/27.0
probsSimple,accep,rej = runTensorFlowRegression(8,alpha)
ax.scatter(accep,rej,c=[c1,c2,c3],label='Alpha: %.1E' %alpha)
it+=1
ax.set_xlabel('signal efficiency')
ax.set_ylabel('background rejection')
plt.legend(loc='lower left', fontsize = 'small');
plt.show()
#now let's investigate how mixed the events are
probsSimple,accep,rej = runTensorFlowRegression(8,.00001)
Signal = df_test.iloc[:,0:1]
#print probsSimple[:,1]
df_test_acc = pd.DataFrame({'PROB':probsSimple[:,1]})
df_test_acc['SIG']=Signal
df_test_acc_sig = df_test_acc.query('SIG==1')
df_test_acc_bkg = df_test_acc.query('SIG==0')
df_test_acc_sig.plot(kind='hist',y='PROB',color='blue',alpha=0.5,bins=np.linspace(0,1,10),label='Signal')
df_test_acc_bkg.plot(kind='hist',y='PROB',color='red',label='Background')
alphas = np.logspace(-10,1,11)
fig = plt.figure()
ax = fig.add_subplot(111)
it=0
for alpha in alphas:
c1 = 1.*( float(it) % 3.)/3.0
c2 = 1.*( float(it) % 9.)/9.0
c3 = 1.*( float(it) % 27.)/27.0
probsSimple,accep,rej = runTensorFlowRegression(18,alpha)
ax.scatter(accep,rej,c=[c1,c2,c3],label='Alpha: %.1E' %alpha)
it+=1
ax.set_xlabel('signal efficiency')
ax.set_ylabel('background rejection')
plt.legend(loc='lower left', fontsize = 'small');
plt.show()
from sklearn.neural_network import MLPClassifier
from sklearn.linear_model import SGDClassifier
def runSciKitRegressionL2(nVar,alpha):
X_train, y_train = getTrainData(nVar)
X_test, y_test = getTestData(nVar)
clf = SGDClassifier(loss="log", penalty="l2",alpha=alpha,max_iter=5,tol=None)
clf.fit(X_train,y_train.ravel())
predictions = clf.predict(X_test)
print('Accuracy on test data with alpha %.2E : %.3f' %(alpha,clf.score(X_test,y_test)) )
probs = clf.predict_proba(X_test)
#print probs
#get signal acceptance and background rejection
thresholds = np.arange(0,1,.01)
Acceptance = []
BkgRejection = []
for thresh in thresholds:
it=0
nPredSig=0.0
nPredBkg=0.0
nTotSig=0.0
nTotBkg=0.0
for prob in probs:
if prob[1]>thresh:
predSig=True
else:
predSig=False
if y_test[it][0]:
Sig = True
else:
Sig = False
if Sig:
if predSig:
nPredSig+=1
nTotSig+=1
else:
if not predSig:
nPredBkg+=1
nTotBkg+=1
it+=1
Acceptance.append(nPredSig/nTotSig)
BkgRejection.append(nPredBkg/nTotBkg)
return (probs,Acceptance,BkgRejection)
def runSciKitRegressionL1(nVar,alpha):
X_train, y_train = getTrainData(nVar)
X_test, y_test = getTestData(nVar)
clf = SGDClassifier(loss="log", penalty="l1",alpha=alpha,max_iter=5,tol=None)
clf.fit(X_train,y_train.ravel())
predictions = clf.predict(X_test)
print('Accuracy on test data with alpha %.2E : %.3f' %(alpha,clf.score(X_test,y_test)) )
probs = clf.predict_proba(X_test)
#print probs
#get signal acceptance and background rejection
thresholds = np.arange(0,1,.01)
Acceptance = []
BkgRejection = []
for thresh in thresholds:
it=0
nPredSig=0.0
nPredBkg=0.0
nTotSig=0.0
nTotBkg=0.0
for prob in probs:
if prob[1]>thresh:
predSig=True
else:
predSig=False
if y_test[it][0]:
Sig = True
else:
Sig = False
if Sig:
if predSig:
nPredSig+=1
nTotSig+=1
else:
if not predSig:
nPredBkg+=1
nTotBkg+=1
it+=1
Acceptance.append(nPredSig/nTotSig)
BkgRejection.append(nPredBkg/nTotBkg)
return (probs,Acceptance,BkgRejection)
alphas = np.logspace(-10,1,11)
fig = plt.figure()
ax = fig.add_subplot(111)
it=0
for alpha in alphas:
c1 = 1.*( float(it) % 3.)/3.0
c2 = 1.*( float(it) % 9.)/9.0
c3 = 1.*( float(it) % 27.)/27.0
probs,accept,rej = runSciKitRegressionL1(8,alpha)
ax.scatter(accept,rej,c=[c1,c2,c3],label='Alpha: %.1E' %alpha)
it+=1
ax.set_xlabel('signal efficiency')
ax.set_ylabel('background rejection')
plt.legend(loc='lower left', fontsize = 'small');
plt.show()
#now let's investigate how mixed the events are
probsSimple,accep,rej = runSciKitRegressionL1(8,.5)
Signal = df_test.iloc[:,0:1]
#print probsSimple[:,1]
df_test_acc = pd.DataFrame({'PROB':probsSimple[:,1]})
df_test_acc['SIG']=Signal
df_test_acc_sig = df_test_acc.query('SIG==1')
df_test_acc_bkg = df_test_acc.query('SIG==0')
df_test_acc_sig.plot(kind='hist',y='PROB',color='blue',alpha=0.5,bins=np.linspace(0,1,10),label='Signal')
df_test_acc_bkg.plot(kind='hist',y='PROB',color='red',label='Background')
alphas = np.logspace(-10,1,11)
fig = plt.figure()
ax = fig.add_subplot(111)
it=0
for alpha in alphas:
c1 = 1.*( float(it) % 3.)/3.0
c2 = 1.*( float(it) % 9.)/9.0
c3 = 1.*( float(it) % 27.)/27.0
probs,accept,rej = runSciKitRegressionL1(18,alpha)
ax.scatter(accept,rej,c=[c1,c2,c3],label='Alpha: %.1E' %alpha)
it+=1
ax.set_xlabel('signal efficiency')
ax.set_ylabel('background rejection')
plt.legend(loc='lower left', fontsize = 'small');
plt.show()
alphas = np.logspace(-10,1,11)
fig = plt.figure()
ax = fig.add_subplot(111)
it=0
for alpha in alphas:
c1 = 1.*( float(it) % 3.)/3.0
c2 = 1.*( float(it) % 9.)/9.0
c3 = 1.*( float(it) % 27.)/27.0
probs,accept,rej = runSciKitRegressionL2(8,alpha)
ax.scatter(accept,rej,c=[c1,c2,c3],label='Alpha: %.1E' %alpha)
it+=1
ax.set_xlabel('signal efficiency')
ax.set_ylabel('background rejection')
plt.legend(loc='lower left', fontsize = 'small');
plt.show()
alphas = np.logspace(-10,1,11)
fig = plt.figure()
ax = fig.add_subplot(111)
it=0
for alpha in alphas:
c1 = 1.*( float(it) % 3.)/3.0
c2 = 1.*( float(it) % 9.)/9.0
c3 = 1.*( float(it) % 27.)/27.0
probs,accept,rej = runSciKitRegressionL2(18,alpha)
ax.scatter(accept,rej,c=[c1,c2,c3],label='Alpha: %.1E' %alpha)
it+=1
ax.set_xlabel('signal efficiency')
ax.set_ylabel('background rejection')
plt.legend(loc='lower left', fontsize = 'small');
plt.show()notebook 5: Logistic Regression, SUSY WARNING:tensorflow:From <ipython-input-9-57e47b7861fc>:61: softmax_cross_entropy_with_logits (from tensorflow.python.ops.nn_ops) is deprecated and will be removed in a future version. Instructions for updating: Future major versions of TensorFlow will allow gradients to flow into the labels input on backprop by default. See tf.nn.softmax_cross_entropy_with_logits_v2. Accuracy for alpha 1.0E-10 : 0.770 Accuracy for alpha 1.3E-09 : 0.770 Accuracy for alpha 1.6E-08 : 0.770 Accuracy for alpha 2.0E-07 : 0.770 Accuracy for alpha 2.5E-06 : 0.770 Accuracy for alpha 3.2E-05 : 0.770 Accuracy for alpha 4.0E-04 : 0.769 Accuracy for alpha 5.0E-03 : 0.764 Accuracy for alpha 6.3E-02 : 0.748 Accuracy for alpha 7.9E-01 : 0.659 Accuracy for alpha 1.0E+01 : 0.543
Accuracy for alpha 1.0E-05 : 0.770 Accuracy for alpha 1.0E-10 : 0.790 Accuracy for alpha 1.3E-09 : 0.790 Accuracy for alpha 1.6E-08 : 0.790 Accuracy for alpha 2.0E-07 : 0.790 Accuracy for alpha 2.5E-06 : 0.790
In [ ]: