kernrrg.py

#!/usr/bin/python

#import modules used for this code. Important to have all of them to prevent future error messages in the terminal

import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
from sklearn.neural_network import MLPRegressor as mlp
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import GridSearchCV
from sklearn.kernel_ridge import KernelRidge

#Tell the program which data file the code will be reading from. These are big data files so they will read slower in the terminal.
raw_data_path = "UW_diffusion_data_for_Brent_with_AN_sorted_by_num.csv" #big csv file
#Tell sage math which row in the data file is the header and to ignore that row. *Not useful for analysis*
raw_data = pd.read_csv(raw_data_path, header=3)
#What are we looking for?
target_value = "E(Impurity) - E(Host) [eV]"

# ************Organize data for fitting************
# Make a copy of the full training frame that has all of the X data columns.
x_data_matrix = raw_data.loc[:,"AREmp (pm) I":"SpaceGroup I"]

# grab just the one column to use as y data. Once as a dataframe and once as a simple list. Create an empty list of y values and load in the text file
target_data = raw_data.loc[:,target_value] #This grabs the column that will be used as the y value for training. This trains the model by pairing the input with expected output.
target_data_list = raw_data[target_value].tolist() #this line of code creates a simplers list version of the pervious line of y values for training data.

# grab just the one column that has the host element labels.
host = raw_data.loc[:,"Host"]  #mini dataframe that has just the host element
hostlist = raw_data['Host'].tolist()  #simple list of the host elements for each row in the x data matrix
impurity = raw_data.loc[:,"Impurity"]  #mini dataframe that as just the impurity element
impuritylist = raw_data['Impurity'].tolist()  #simple list of the impurities for each row in the x data matrix
host_element_list = ["Al","Au","Ca","Cu","Ir","Ni","Pb","Pd","Pt"]  #list of host elements in case it is useful


# ************Scaling and Fitting************
#Kernel Ridge and scaling implemenation.
#Scales the procceessed data frames
scaler = StandardScaler()
scaler.fit(x_data_matrix) #Train the scaler with x training data frame
x_data_matrix = scaler.transform(x_data_matrix)  #Now scale the x training frame


#Kernel Ridge Regression
# This sets up the GridSearch process.
# Here, alpha and gamma will be optimized using CV.
#Note that more than one thing can be optimized simultaneously, just put them all within the curly braces.
regr = GridSearchCV(KernelRidge(kernel='polynomial', coef0=1), cv=20,
                    param_grid={"alpha": [.1,1e-2, 1e-3],
                                "gamma": [10, 25, 50, 75, 100]})
#Without GridSearchCV, the best values for each coefficient can be determined by the outputted R^2 value. Better R^2 values resulted in higher gamma (100) and smaller alpha (.0001).
#regr = KernelRidge(alpha=.000001, coef0=1, kernel='polynomial', gamma=100)
#regr = KernelRidge(coef0=1, kernel='polynomial', kernel_params = {'alpha':[1e10, 0.01, 1e-2, 1e-3, 1e-4], 'gamma': [-10, 1, 10, 100]})


regr.fit(x_data_matrix, target_data)
#predicted_data=regr.predict(validation_frame)
predicted_data=regr.predict(x_data_matrix)
#predicted_data = grid.predict(x_data_matrix)  #This output is just a simple list of numbers
fit_error = predicted_data - target_data

#Print the R^2 value. R-squared is a statistical measure of how close the data are to the fitted regression line. It is also known as the coefficient of determination for multiple regression. 
#The closer R^2 is to 1 indicates that the model explains the variability of the response data around its mean very well.
print "R^2: " +str(regr.score(x_data_matrix, target_data))
#print "R^2: " +str(grid.score(x_data_matrix, target_data))


# ************Plotting the Data************
# We will make one big figure with lots of little subplots -- one for each host element.

# Setup the full figure with nine subplots
fullfig, (axAl, axAu, axCa, axCu, axIr, axNi, axPb, axPd, axPt) = plt.subplots(9)
fullfig.subplots_adjust(hspace=0.5)  #sets the spacing between plots -- default was too tight.
fullfig.suptitle('Activation Energies (eV)')  #overall title for the whole thing
fullfig.set_figheight(20)  #default is to squash all subplots together. This gives more space.

#Host Al, rows 0-41 (42 entries)
axAl.set_ylabel('Host Al')
x = range(0,42)
# Set number of ticks for x-axis
axAl.set_xticks(x)
# Set ticks labels for x-axis
axAl.set_xticklabels(impuritylist[0:42], rotation=270, fontsize=12)
axAl.plot(x, target_data_list[0:42], marker='o', ls='', color='black', label='Test Data')
axAl.plot(x, predicted_data[0:42], color='red', label='Fit')

#Host Au, rows 42-52 (11 entries)
axAu.set_ylabel('Host Au')
#print impuritylist[42:53]   #To double-check that we are plotting the right data
#print target_data_list[42:53]
x = range(0,11)
# Set number of ticks for x-axis
axAu.set_xticks(x)
# Set ticks labels for x-axis
axAu.set_xticklabels(impuritylist[42:53], rotation=270, fontsize=12)
axAu.plot(x, target_data_list[42:53], marker='o', ls='', color='black', label='Test Data')
axAu.plot(x, predicted_data[42:53], color='red', label='Fit')


#Host Ca, rows 53-58 (6 entries)
axCa.set_ylabel('Host Ca')
#print impuritylist[53:59]   #To double-check that we are plotting the right data
#print target_data_list[53:59]
x = range(0,6)
# Set number of ticks for x-axis
axCa.set_xticks(x)
# Set ticks labels for x-axis
axCa.set_xticklabels(impuritylist[53:59], rotation=270, fontsize=12)
axCa.plot(x, target_data_list[53:59], marker='o', ls='', color='black', label='Test Data')
axCa.plot(x, predicted_data[53:59], color='red', label='Fit')

#Host Cu, rows 59-98 (39 entries)
axCu.set_ylabel('Host Cu')
#print impuritylist[59-98]   #To double-check that we are plotting the right data
#print target_data_list[59-98]
x = range(0,39)
# Set number of ticks for x-axis
axCu.set_xticks(x)
# Set ticks labels for x-axis
axCu.set_xticklabels(impuritylist[59:98], rotation=270, fontsize=12)
axCu.plot(x, target_data_list[59:98], marker='o', ls='', color='black', label='Test Data')
axCu.plot(x, predicted_data[59:98], color='red', label='Fit')

#Host Ir, rows 99-109
axIr.set_ylabel('Host Ir')
#print impuritylist[99-110]   #To double-check that we are plotting the right data
#print target_data_list[99-110]
x = range(0,10)
# Set number of ticks for x-axis
axIr.set_xticks(x)
# Set ticks labels for x-axis
axIr.set_xticklabels(impuritylist[99:109], rotation=270, fontsize=12)
axIr.plot(x, target_data_list[99:109], marker='o', ls='', color='black', label='Test Data')
axIr.plot(x, predicted_data[99:109], color='red', label='Fit')

#Host Ni, rows 110-149
axNi.set_ylabel('Host Ni')
#print impuritylist[110-149]   #To double-check that we are plotting the right data
#print target_data_list[110-149]
x = range(0,39)
# Set number of ticks for x-axis
axNi.set_xticks(x)
# Set ticks labels for x-axis
axNi.set_xticklabels(impuritylist[110:149], rotation=270, fontsize=12)
axNi.plot(x, target_data_list[110:149], marker='o', ls='', color='black', label='Test Data')
axNi.plot(x, predicted_data[110:149], color='red', label='Fit')


#Host Pb, rows 150-157
axPb.set_ylabel('Host Pb')
#print impuritylist[150-157]   #To double-check that we are plotting the right data
#print target_data_list[150-157]
x = range(0,7)
# Set number of ticks for x-axis
axPb.set_xticks(x)
# Set ticks labels for x-axis
axPb.set_xticklabels(impuritylist[150:157], rotation=270, fontsize=12)
axPb.plot(x, target_data_list[150:157], marker='o', ls='', color='black', label='Test Data')
axPb.plot(x, predicted_data[150:157], color='red', label='Fit')

#Host Pd, rows 158-187
axPd.set_ylabel('Host Pd')
#print impuritylist[158-187]   #To double-check that we are plotting the right data
#print target_data_list[158-187]
x = range(0,29)
# Set number of ticks for x-axis
axPd.set_xticks(x)
# Set ticks labels for x-axis
axPd.set_xticklabels(impuritylist[158:187], rotation=270, fontsize=12)
axPd.plot(x, target_data_list[158:187], marker='o', ls='', color='black', label='Test Data')
axPd.plot(x, predicted_data[158:187], color='red', label='Fit')

#Host Pt, rows 188-217
axPt.set_ylabel('Host Pt')
#print impuritylist[188-217]   #To double-check that we are plotting the right data
#print target_data_list[188-217]
x = range(0,29)
# Set number of ticks for x-axis
axPt.set_xticks(x)
# Set ticks labels for x-axis
axPt.set_xticklabels(impuritylist[188:217], rotation=270, fontsize=12)
axPt.plot(x, target_data_list[188:217], marker='o', ls='', color='black', label='Test Data')
axPt.plot(x, predicted_data[188:217], color='red', label='Fit')


fullfig.savefig("activation_energies.png")