Random forest and machine learning

Question

I am quite new to using python for machine learning. I come from a background of programming in Fortran, so as you may imagine, python is quite a leap. I work in chemistry and have become involved in chemiformatics (applying data science techniques to chemistry). As such, the application of pythons extensive machine learning libraries is important. I also need my codes to be efficent. I have written a code which runs and seems to work OK. What I would like to know is:

1. How best to improve it/make it more efficient.

2. Any suggestions on alternative formulations to those I have used and if possible a reason why another route maybe superior?

I tend to work with continuous data and regression models.

Edit:

Thank you for all of your comments so far. Apologies for the indentation error this was a copy mistake.

To give a few more details I am aiming to use the code to make prediction of chemical properties such as toxicity, melting points, solubility etc. These sorts of properties are the focus of research efforts in academia and industry to provide pre-screening of target molecules for certain properties.

The data I am providing as input is a csv file. The first column, is a label (molecule name). The last column, is the target value from experiment or quantum chemical calculation. The in between columns are descriptors calculated based on some molecular structure format (2D SMILES, 3D crystal structure etc). An example of a subset of the data is below:

Typically there would be 100 - 150 descriptors which provide information and between 100 and several thousand examples. These examples need to be split into training and test sets.

End Edit

import scipy
import math
import numpy as np
import pandas as pd
import plotly.plotly as py
import os.path
import sys

from time import time
from sklearn import preprocessing, metrics, cross_validation
from sklearn.cross_validation import train_test_split
from sklearn.ensemble import RandomForestRegressor
from sklearn.grid_search import GridSearchCV
from sklearn.cross_validation import KFold

fname = str(raw_input('Please enter the input file name containing total dataset and descriptors (assumes csv file, column headings and first column are labels\n'))
if os.path.isfile(fname) :
    SubFeAll = pd.read_csv(fname, sep=",")
else:
    sys.exit("ERROR: input file does not exist")

#SubFeAll = pd.read_csv(fname, sep=",")
SubFeAll = SubFeAll.fillna(SubFeAll.mean()) # replace the NA values with the mean of the descriptor
header = SubFeAll.columns.values # Use the column headers as the descriptor labels
SubFeAll.head()

# Set the numpy global random number seed (similar effect to random_state) 
np.random.seed(1)  

# Random Forest results initialised
RFr2 = []
RFmse = []
RFrmse = []

# Predictions results initialised 
RFpredictions = []

metcount = 0

# Give the array from pandas to numpy
npArray = np.array(SubFeAll)
print header.shape
npheader = np.array(header[1:-1])
print("Array shape X = %d, Y = %d " % (npArray.shape))
datax, datay =  npArray.shape

# Print specific nparray values to check the data
print("The first element of the input data set, as a minial check please ensure this is as expected = %s" % npArray[0,0])

# Split the data into: names labels of the molecules ; y the True results ; X the descriptors for each data point
names = npArray[:,0]
X = npArray[:,1:-1].astype(float)
y = npArray[:,-1] .astype(float)
X = preprocessing.scale(X)
print X.shape

# Open output files
train_name = "Training.csv"
test_name = "Predictions.csv"
fi_name = "Feature_importance.csv"

with open(train_name,'w') as ftrain, open(test_name,'w') as fpred, open(fi_name,'w') as ffeatimp:
        ftrain.write("This file contains the training information for the Random Forest models\n")
        ftrain.write("The code use a ten fold cross validation 90% training 10% test at each fold so ten training sets are used here,\n")
        ftrain.write("Interation %d ,\n" %(metcount+1))

        fpred.write("This file contains the prediction information for the Random Forest models\n")
        fpred.write("Predictions are made over a ten fold cross validation hence training on 90% test on 10%. The final prediction are return iteratively over this ten fold cros validation once,\n")
        fpred.write("optimised parameters are located via a grid search at each fold,\n")
        fpred.write("Interation %d ,\n" %(metcount+1))

        ffeatimp.write("This file contains the feature importance information for the Random Forest model,\n")
        ffeatimp.write("Interation %d ,\n" %(metcount+1))

        # Begin the K-fold cross validation over ten folds
        kf = KFold(datax, n_folds=10, shuffle=True, random_state=0)
        print "------------------- Begining Ten Fold Cross Validation -------------------"
        for train, test in kf:
            XTrain, XTest, yTrain, yTest = X[train], X[test], y[train], y[test]
            ytestdim = yTest.shape[0]
                print("The test set values are : ")
                i = 0
                if ytestdim%5 == 0:
                        while i < ytestdim:
                                print round(yTest[i],2),'\t', round(yTest[i+1],2),'\t', round(yTest[i+2],2),'\t', round(yTest[i+3],2),'\t', round(yTest[i+4],2)
                                ftrain.write(str(round(yTest[i],2))+','+ str(round(yTest[i+1],2))+','+str(round(yTest[i+2],2))+','+str(round(yTest[i+3],2))+','+str(round(yTest[i+4],2))+',\n')
                                i += 5
                elif ytestdim%4 == 0:
                        while i < ytestdim:
                                print round(yTest[i],2),'\t', round(yTest[i+1],2),'\t', round(yTest[i+2],2),'\t', round(yTest[i+3],2)
                                ftrain.write(str(round(yTest[i],2))+','+str(round(yTest[i+1],2))+','+str(round(yTest[i+2],2))+','+str(round(yTest[i+3],2))+',\n')
                                i += 4
                elif ytestdim%3 == 0 :
                        while i < ytestdim :
                                print round(yTest[i],2),'\t', round(yTest[i+1],2),'\t', round(yTest[i+2],2)
                                ftrain.write(str(round(yTest[i],2))+','+str(round(yTest[i+1],2))+','+str(round(yTest[i+2],2))+',\n')
                                i += 3
                elif ytestdim%2 == 0 :
                        while i < ytestdim :
                                print round(yTest[i],2), '\t', round(yTest[i+1],2)
                                ftrain.write(str(round(yTest[i],2))+','+str(round(yTest[i+1],2))+',\n')
                                i += 2
                        else :
                                while i< ytestdim :
                                        print round(yTest[i],2)
                                        ftrain.write(str(round(yTest[i],2))+',\n')
                                        i += 1        

                print "\n"
                # random forest grid search parameters
            print "------------------- Begining Random Forest Grid Search -------------------"
                rfparamgrid = {"n_estimators": [10], "max_features": ["auto", "sqrt", "log2"], "max_depth": [5,7]}
                rf = RandomForestRegressor(random_state=0,n_jobs=2)
                RfGridSearch = GridSearchCV(rf,param_grid=rfparamgrid,scoring='mean_squared_error',cv=10)
                start = time()
                RfGridSearch.fit(XTrain,yTrain)

                # Get best random forest parameters
                print("GridSearchCV took %.2f seconds for %d candidate parameter settings" %(time() - start,len(RfGridSearch.grid_scores_)))
                RFtime = time() - start,len(RfGridSearch.grid_scores_)
                #print(RfGridSearch.grid_scores_)  # Diagnos
                print("n_estimators = %d " % RfGridSearch.best_params_['n_estimators'])
                ne = RfGridSearch.best_params_['n_estimators']
                print("max_features = %s " % RfGridSearch.best_params_['max_features'])
                mf = RfGridSearch.best_params_['max_features']
                print("max_depth = %d " % RfGridSearch.best_params_['max_depth'])
                md = RfGridSearch.best_params_['max_depth']

                ftrain.write("Random Forest")
                ftrain.write("RF search time, %s ,\n" % (str(RFtime)))
                ftrain.write("Number of Trees, %s ,\n" % str(ne))
                ftrain.write("Number of feature at split, %s ,\n" % str(mf))
                ftrain.write("Max depth of tree, %s ,\n" % str(md))

                # Train random forest and predict with optimised parameters
                print("\n\n------------------- Starting opitimised RF training -------------------")
                optRF = RandomForestRegressor(n_estimators = ne, max_features = mf, max_depth = md, random_state=0)
                optRF.fit(XTrain, yTrain)       # Train the model
                RFfeatimp = optRF.feature_importances_
                indices = np.argsort(RFfeatimp)[::-1]
                print("Training R2 = %5.2f" % optRF.score(XTrain,yTrain))
                print("Starting optimised RF prediction")
                RFpreds = optRF.predict(XTest)
                print("The predicted values now follow :")
                RFpredsdim = RFpreds.shape[0]
                i = 0
                if RFpredsdim%5 == 0:
                        while i < RFpredsdim:
                                print round(RFpreds[i],2),'\t', round(RFpreds[i+1],2),'\t', round(RFpreds[i+2],2),'\t', round(RFpreds[i+3],2),'\t', round(RFpreds[i+4],2)
                                i += 5
                elif RFpredsdim%4 == 0:
                        while i < RFpredsdim:
                                print round(RFpreds[i],2),'\t', round(RFpreds[i+1],2),'\t', round(RFpreds[i+2],2),'\t', round(RFpreds[i+3],2)
                                i += 4
                elif RFpredsdim%3 == 0 :
                        while i < RFpredsdim :
                                print round(RFpreds[i],2),'\t', round(RFpreds[i+1],2),'\t', round(RFpreds[i+2],2)
                                i += 3
                elif RFpredsdim%2 == 0 :
                        while i < RFpredsdim :
                                print round(RFpreds[i],2), '\t', round(RFpreds[i+1],2)
                                i += 2
                else :
                        while i< RFpredsdim :
                                print round(RFpreds[i],2)
                i += 1
                print "\n"
                RFr2.append(optRF.score(XTest, yTest))
                RFmse.append( metrics.mean_squared_error(yTest,RFpreds))
                RFrmse.append(math.sqrt(RFmse[metcount]))
                print ("Random Forest prediction statistics for fold %d are; MSE = %5.2f RMSE = %5.2f R2 = %5.2f\n\n" % (metcount+1, RFmse[metcount], RFrmse[metcount],RFr2[metcount]))

                ftrain.write("Random Forest prediction statistics for fold %d are, MSE =, %5.2f, RMSE =, %5.2f, R2 =, %5.2f,\n\n" % (metcount+1, RFmse[metcount], RFrmse[metcount],RFr2[metcount]))



                ffeatimp.write("Feature importance rankings from random forest,\n")
                for i in range(RFfeatimp.shape[0]) :
                        ffeatimp.write("%d. , feature %d , %s,  (%f),\n" % (i + 1, indices[i], npheader[indices[i]], RFfeatimp[indices[i]]))


                # Store prediction in original order of data (itest) whilst following through the current test set order (j)
            metcount += 1

                ftrain.write("Fold %d, \n" %(metcount))

            print "------------------- Next Fold %d -------------------" %(metcount+1)
            j = 0
            for itest in test :
                RFpredictions.append(RFpreds[j])
                j += 1


        lennames = names.shape[0]
        lenpredictions = len(RFpredictions)
        lentrue = y.shape[0]
        if lennames == lenpredictions == lentrue :
                fpred.write("Names/Label,, Prediction Random Forest,, True Value,\n") 
                for i in range(0,lennames) :
                        fpred.write(str(names[i])+",,"+str(RFpredictions[i])+",,"+str(y[i])+",\n")
        else :
                fpred.write("ERROR - names, prediction and true value array size mismatch. Dumping arrays for manual inspection in predictions.csv\n")
                fpred.write("Array printed in the order names/Labels, predictions RF and true values\n")
                fpred.write(names+"\n")
                fpred.write(RFpredictions+"\n")
                fpred.write(y+"\n")
                sys.exit("ERROR - names, prediction and true value array size mismatch. Dumping arrays for manual inspection in predictions.csv")

        print "Final averaged Random Forest metrics : "
        RFamse  = sum(RFmse)/10
        RFmse_sd = np.std(RFmse)
        RFarmse = sum(RFrmse)/10
        RFrmse_sd = np.std(RFrmse)
        RFslope, RFintercept, RFr_value, RFp_value, RFstd_err = scipy.stats.linregress(RFpredictions, y)
        RFR2 = RFr_value**2
        print "Average Mean Squared Error = ", RFamse, " +/- ", RFmse_sd 
        print "Average Root Mean Squared Error = ", RFarmse, " +/- ", RFrmse_sd
        print "R2 Final prediction against True values = ", RFR2

        fpred.write("\n")
        fpred.write("FINAL PREDICTION STATISTICS,\n")
        fpred.write("Random Forest average MSE, %s, +/-, %s,\n" %(str(RFamse), str(RFmse_sd)))
        fpred.write("Random Forest average RMSE, %s, +/-, %s,\n" %(str(RFarmse), str(RFrmse_sd)))
        fpred.write("Random Forest slope, %s,   Random Forest intercept, %s,\n" %(str(RFslope), str(RFintercept)))
        fpred.write("Random Forest standard error, %s,\n" %(str(RFstd_err)))
        fpred.write("Random Forest R, %s,\n" %(str(RFr_value)))
        fpred.write("Random Forest R2, %s,\n" %(str(RFR2)))

ftrain.close()
fpred.close()
ffeatimp.close()

Answer 1

Disclaimer: I'm not familiar with the work you're doing, so I'll limit my comments to the structure of your code. Generally, I'd make the following changes:

• Split data and code
• Write more functions
• Don't repeat yourself

Split data and code

You have a lot of lines that contain both hardcoded information to display to the user and calculations that are being performed. At the very least, I could imagine splitting this into two classes: a calculation class and a reporter. Put all the logic for extracting information in to the calculation class, pass it to the reporter to display to the user. The calculation class can init all of your empty arrays and contain methods for operating on data.

Perhaps later you won't want to print everything, you'll want to save it to a logfile, put on a website, etc. Perhaps you'll want to report everything the same way, but you want to compare two different algorithms. Separating the logic and how you report the results will make this much simpler.

Much of the information on what to do with what data is tied up in lines like this:

print round(RFpreds[i],2),'\t', round(RFpreds[i+1],2),'\t', round(RFpreds[i+2],2),'\t', round(RFpreds[i+3],2),'\t', round(RFpreds[i+4],2)

Write more functions

An excellent piece of advice I've heard is that if you write a block of code and put a comment at the top, you've probably just written something that should be a function. Something like the collection of if checks related to ytestdim seems like it could be abstracted out into its own function.

This allows you to review what your code is doing at each level of abstraction. Ultimately this makes debugging a lot easier, and it makes it especially easy to compare a mathematical technique you're familiar with in the real world with your code implementation of it.

Don't repeat yourself

I also see plenty of opportunities to combine collections of statements. For example:

print("n_estimators = %d " % RfGridSearch.best_params_['n_estimators'])
ne = RfGridSearch.best_params_['n_estimators']
print("max_features = %s " % RfGridSearch.best_params_['max_features'])
mf = RfGridSearch.best_params_['max_features']
print("max_depth = %d " % RfGridSearch.best_params_['max_depth'])
md = RfGridSearch.best_params_['max_depth']

This is really a single operation performed 3 times. In situations like this, make a list ["n_estimators", "max_features", "max_depth"], and iterate over the list.

Making changes like this prevent you from making simple spelling errors from one line to the next, and make it easier to pull behavior into classes and methods.

Miscellaneous additional items:

If you're using the with...open() as notation, you don't need to worry about close().

You can print multiple things on the same line with this notation:

print a,
print b,
print c

It's easier to read dictionaries that are defined the following way:

rfparamgrid = {
    "n_estimators": [10],
    "max_features": ["auto", "sqrt", "log2"],
    "max_depth": [5,7],
}