kernridgeCV1.py

#!/usr/bin/python
# Import libraries needed
import sys
import numpy as np
#from sklearn.svm import SVR
#from sklearn.neighbors import KNeighborsRegressor
from sklearn.kernel_ridge import KernelRidge
from sklearn.linear_model import Ridge
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import classification_report
from sklearn.svm import SVC
import argparse
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline
from sklearn import linear_model
from __future__ import print_function
from sklearn import datasets
from sklearn.model_selection import train_test_split

#Create the argument parser
parser = argparse.ArgumentParser (description = "ploting a polynomial fit of x and y data values")
parser.add_argument("input")
parser.add_argument("degree", type = int)
#parser.add_argument("C", type = int)
arguments = parser.parse_args()


# Load the traing and validation sets
validate_data = arguments.input + ".validate.xy"
training_data = arguments.input + ".training.xy"
plot_name = arguments.input + ".simple_polynomial_fit.png"

#Create an empty list of x and y values and load in the text file and interpret it as an array (2 coloumns, x and y, with 100 rows)
tx_y_data = np.loadtxt(training_data, delimiter = '\t')
# This x_y_data array is two columns, just as the text file was. Split the input array into an array for x and a separate one for y
xtdata = tx_y_data[:,0] #all rows in the 0th column
ytdata = tx_y_data[:,1] #all

# X and y data needs to be numpy array [n_samples,n_features]. For simple x-y data, this is [100,1] (if we have 100 data points) [rows,columns]
#print the shape of the x data. As of now it is a list with 100 elements. As read from disk these 1-D arrays are size [100, ]. *Convert a simple list to a [1,100] array, and then transpose to make it [100,1]*
#This is done because the previous code is just a list with no dimensiality. Turns into an array with a long row of 100 elements. The transform makes it 100 rows with 1 column, for both x and y.
xtdata = np.array([xtdata])
xtdata = xtdata.T
#Print the shape of the x data

#repeat for y data
#ytdata = np.array([ytdata])
#ytdata = ytdata.T
#print the shape of y data

#Repeat the process above for validate file
vx_y_data = np.loadtxt(validate_data, delimiter = '\t')
# This x_y_data array is two columns, just as the text file was. Split the input array into an array for x and a separate one for y
x_vdata = vx_y_data[:,0] #all rows in the 0th column
y_vdata = vx_y_data[:,1] #all

# X and y data needs to be numpy array [n_samples,n_features]. For simple x-y data, this is [100,1] (if we have 100 data points) [rows,columns]
#print the shape of the x data. As of now it is a list with 100 elements. As read from disk these 1-D arrays are size [100, ]. *Convert a simple list to a [1,100] array, and then transpose to make it [100,1]*
x_vdata = np.array([x_vdata])
x_vdata = x_vdata.T
#Print the shape of the x data

#repeat for y data
#y_vdata = np.array([y_vdata])
#y_vdata = y_vdata.T


fit_stdev = list()
# Define the model that we will be using to do the fit
#regr = linear_model.LinearRegression(fit_intercept=True)  # just a simple LinearRegression fit
#regr = make_pipeline( PolynomialFeatures(degree=arguments.degree), linear_model.LinearRegression(fit_intercept=True) )
#regr = SVR(degree=arguments.degree, gamma='auto', coef0=0.0, tol=10, C=arguments.C, epsilon=35, shrinking=True, cache_size=200, verbose=False, max_iter=-1, kernel='poly')
#regr = NearestNeighbors(n_neighbors=2, algorithm='auto').fit(x_vdata)
#distances, indices = nbrs.kneighbors(X)
#regr = KNeighborsRegressor(n_neighbors=10, weights='distance', algorithm='auto', leaf_size=30, p=1, metric='minkowski', metric_params=None, n_jobs=1)
kr = KernelRidge(alpha=1.0, coef0=1, kernel='polynomial', gamma=None, degree = arguments.degree, kernel_params=None)
kr.fit (xtdata, ytdata)
fit_stdev.append( np.sqrt( np.mean((kr.predict(x_vdata) - y_vdata) ** 2) ) )
print "The standard deviation of the fit is:", fit_stdev

print(__doc__)

# Loading the Digits dataset
digits = datasets.load_digits()

# To apply an classifier on this data, we need to flatten the image, to
# turn the data in a (samples, feature) matrix:
n_samples = len(digits.images)
X = digits.images.reshape((n_samples, -1))
y = digits.target

# Split the dataset in two equal parts
#X_train, X_test, y_train, y_test = train_test_split(
 #   X, y, test_size=0.5, random_state=0)

# Set the parameters by cross-validation
tuned_parameters = [{'kernel': ['rbf'], 'gamma': [1e-3, 1e-4],
                     'C': [1, 10, 100, 1000]},
                    {'kernel': ['polynomial'], 'C': [1, 10, 100, 1000]}]
scores = ['precision', 'recall']
for score in scores:
    print("# Tuning hyper-parameters for %s" % score)
    print()

    clf = GridSearchCV(SVC(C=1), tuned_parameters, cv=5,
                       scoring='%s_macro' % score)
    clf.fit(xtdata, ytdata)

    print("Best parameters set found on development set:")
    print()
    print(clf.best_params_)
    print()
    print("Grid scores on development set:")
    print()
    means = clf.cv_results_['mean_test_score']
    stds = clf.cv_results_['std_test_score']
    for mean, std, params in zip(means, stds, clf.cv_results_['params']):
        print("%0.3f (+/-%0.03f) for %r"
              % (mean, std * 2, params))
    print()

    print("Detailed classification report:")
    print()
    print("The model is trained on the full development set.")
    print("The scores are computed on the full evaluation set.")
    print()
    y_true, y_pred = y_test, clf.predict(X_test)
    print(classification_report(y_true, y_pred))
    print()

R = Ridge(alpha=1.0,fit_intercept=False)
kr.fit(xtdata, ytdata)
R.fit(xtdata, ytdata)
print np.dot(xtdata.transpose(),kr.dual_coef_)
print R.coef_

#fit_coefficients = regr.coef_ to pull out the linear regression
#fit_coefficients = regr.named_steps['linearregression'].coef_
#print 'The coefficients of the regression are' +str(fit_coefficients)
#fit_intercept = regr.intercept_
#fit_intercept = regr.named_steps['linearregression'].intercept_
#print "The intercept of the line is %d" % (fit_intercept)

# Evaluate how well the model did, usually using the validation data
# Now evaluate how the model performs on the validation data. Calculate the mean squared error, predict and fit will take in data, calculate a y and x and compare. In prediction, use those parameters with given x data to spit out some y data.  Take a difference between that and the true y value to get the errors
mean_sq_error = np.mean((kr.predict(x_vdata) - y_vdata) ** 2)
sqrt_mean_sq_error = (mean_sq_error) ** .5
# The variance score, determined by the ML fitting: 1 is perfect prediction. Score uses the fit results we have to determine a normalized mean squared erro\]\AASASDF"    wwq\`"AS
variance_score = kr.score(x_vdata, y_vdata)
print "The varience score is %d" % (variance_score)
print "The mean sqaured error of the fit is %d" % (mean_sq_error)
print "The square root of the mean squared error is %d" % (sqrt_mean_sq_error)
#print 'regr is:\n', regr

# Plot training data
plt.plot (xtdata, ytdata, 'bo', label="Training Data")
# Plot validation data
plt.plot (x_vdata, y_vdata, 'ro', label="Validation Data")
# Plot fit line
plt.plot (x_vdata, kr.predict(x_vdata), label="Fit Function with degree %d." % (arguments.degree))
#Create a legend
plt.legend(shadow=True, loc="upper left")

#Lable the plot
plt.xlabel('Time (sec)')
plt.ylabel('Signal Strength')
plt.title('Signal Strength vs. Time for Relative Activation Energies')
plt.savefig("kernel_ridge_poly_fit.png")