This notebook is an example of how to generate a working submission for the Kaggle Plankton competition based on the offical competition tutorial.
This Notebook was adapted from Ehud Ben-Reuven's file on GitHub. The file's original name was 141215-tutorial-submission.ipynb.
#Import libraries for doing image analysis
from skimage import measure
from skimage import morphology
from skimage.io import imread
from skimage.transform import resize
from sklearn.ensemble import RandomForestClassifier as RF
from sklearn.cross_validation import StratifiedKFold as KFold
import glob
import numpy as np
import os
import pandas
import warnings
warnings.filterwarnings("ignore")
import datetime
start_time = datetime.datetime.now()
print start_time
Set the random number seed so results are reproducible.
# 3.716 for 17 for 25-by-25
# 3.51 for 71 for 64-by-64
# 3.405 for 19937 for 128-by-128 [used all 32 GB RAM]
np.random.seed(71)
cd C:\Kaggle\2015\Plankton
The training data is organized in a series of subdirectories that contain examples for the each class of interest. We will store the list of directory names to aid in labelling the data classes for training and testing purposes.
# The directory structure gives the plankton classes
directory_names = os.listdir("train")
# find the largest nonzero region
def getLargestRegion(props, labelmap, imagethres):
regionmaxprop = None
for regionprop in props:
# check to see if the region is at least 50% nonzero
if sum(imagethres[labelmap == regionprop.label])*1.0/regionprop.area < 0.50:
continue
if regionmaxprop is None:
regionmaxprop = regionprop
if regionmaxprop.filled_area < regionprop.filled_area:
regionmaxprop = regionprop
return regionmaxprop
def getMinorMajorRatio(image):
image = image.copy()
# Create the thresholded image to eliminate some of the background.
# Use median as more robust measure of central tendency than mean.
imagethr = np.where(image > np.mean(image),0.,1.0)
#Dilate the image
imdilated = morphology.dilation(imagethr, np.ones((4,4)))
# Create the label list
label_list = measure.label(imdilated)
label_list = imagethr*label_list
label_list = label_list.astype(int)
region_list = measure.regionprops(label_list)
maxregion = getLargestRegion(region_list, label_list, imagethr)
# guard against cases where the segmentation fails by providing zeros
ratio = 0.0
if ((not maxregion is None) and (maxregion.major_axis_length != 0.0)):
ratio = 0.0 if maxregion is None else maxregion.minor_axis_length*1.0 / maxregion.major_axis_length
return ratio
With our code for the ratio of minor to major axis, let's add the raw pixel values to the list of features for our dataset. In order to use the pixel values in a model for our classifier, we need a fixed length feature vector, so we will rescale the images to be constant size and add the fixed number of pixels to the feature vector.
To create the feature vectors, we will loop through each of the directories in our training data set and then loop over each image within that class. For each image, we will rescale it to 25 x 25 pixels and then add the rescaled pixel values to a feature vector, X. The last feature we include will be our width-to-length ratio. We will also create the class label in the vector y, which will have the true class label for each row of the feature vector, X.
folder = "acantharia_protist"
x = [fileNameDir for fileNameDir in os.walk(os.path.join("train", folder))]
x[0][0]
fileName = x[0][2][0] # first image filename
print fileName
nameFileImage = os.path.join(fileNameDir[0], fileName)
print nameFileImage
image = imread(nameFileImage, as_grey=True)
print image.shape
print np.mean(image)
print np.median(image)
axisratio = getMinorMajorRatio(image)
axisratio
#from matplotlib import pyplot as plt
#from pylab import cm
#plt.imshow(image, cmap=cm.gray)
#plt.show()
# Rescale the images and create the combined metrics and training labels
# Count .jpgs in train folder
numberofImages = 0
for folder in directory_names:
# os.walk gives (dirpath, dirnames, filenames)
for fileNameDir in os.walk(os.path.join("train", folder)):
for fileName in fileNameDir[2]: # filenames
numberofImages += (fileName[-4:] == ".jpg")
# We'll rescale the images to be 1: 25x25, 2:64x64, 3:128x128
maxPixel = 25
imageSize = maxPixel * maxPixel
num_rows = numberofImages # one row for each image in the training dataset
num_features = imageSize + 1 # for ratio feature
# X is the feature vector with one row of features per image
# consisting of the pixel values and our metric
X = np.zeros((num_rows, num_features), dtype=float)
print X.shape
# y is the numeric class label
y = np.zeros((num_rows))
print y.shape
files = []
# Generate training data
i = 0
label = 0
# List of string of class names
namesClasses = list()
print "Reading images"
# Navigate through the list of directories
for folder in directory_names:
# Append the string class name for each class
currentClass = folder.split(os.pathsep)[-1]
namesClasses.append(currentClass)
for fileNameDir in os.walk(os.path.join("train",folder)):
for fileName in fileNameDir[2]: # filenames
# Only read in the images
if fileName[-4:] == ".jpg":
# Read in the images and create the features
nameFileImage = os.path.join(fileNameDir[0], fileName)
image = imread(nameFileImage, as_grey=True)
files.append(nameFileImage)
axisratio = getMinorMajorRatio(image)
image = resize(image, (maxPixel, maxPixel))
# Store the rescaled image pixels and the axis ratio
X[i, 0:imageSize] = np.reshape(image, (1, imageSize))
X[i, imageSize] = axisratio
# Store the classlabel
y[i] = label
i += 1
# report progress for each 5% done
report = [int((j+1)*num_rows/20.) for j in range(20)]
if i in report: print np.ceil(i *100.0 / num_rows), "% done"
label += 1
numberofImages
We choose a random forest model to classify the images. Random forests perform well in many classification tasks and have robust default settings. We will give a brief description of a random forest model so that you can understand its two main free parameters: n_estimators and max_features.
A random forest model is an ensemble model of n_estimators number of decision trees. During the training process, each decision tree is grown automatically by making a series of conditional splits on the data. At each split in the decision tree, a random sample of max_features number of features is chosen and used to make a conditional decision on which of the two nodes that the data will be grouped in. The best condition for the split is determined by the split that maximizes the class purity of the nodes directly below. The tree continues to grow by making additional splits until the leaves are pure or the leaves have less than the minimum number of samples for a split (in sklearn default for min_samples_split is two data points). The final majority class purity of the terminal nodes of the decision tree are used for making predictions on what class a new data point will belong. Then, the aggregate vote across the forest determines the class prediction for new samples.
The competition scoring uses a multiclass log-loss metric to compute your overall score. In the next steps, we define the multiclass log-loss function and compute your estimated score on the training dataset.
def multiclass_log_loss(y_true, y_pred, eps=1e-15):
"""Multi class version of Logarithmic Loss metric.
https://www.kaggle.com/wiki/MultiClassLogLoss
Parameters
----------
y_true : array, shape = [n_samples]
true class, intergers in [0, n_classes - 1)
y_pred : array, shape = [n_samples, n_classes]
Returns
-------
loss : float
"""
predictions = np.clip(y_pred, eps, 1 - eps)
# normalize row sums to 1
predictions /= predictions.sum(axis=1)[:, np.newaxis]
actual = np.zeros(y_pred.shape)
n_samples = actual.shape[0]
actual[np.arange(n_samples), y_true.astype(int)] = 1
vectsum = np.sum(actual * np.log(predictions))
loss = -1.0 / n_samples * vectsum
return loss
# Get the probability predictions for computing the log-loss function
kf = KFold(y, n_folds=5)
# prediction probabilities number of samples, by number of classes
y_pred = np.zeros((len(y),len(set(y))))
for train, test in kf:
X_train, X_test, y_train, y_test = X[train,:], X[test,:], y[train], y[test]
clf = RF(n_estimators=100, n_jobs=3)
clf.fit(X_train, y_train)
y_pred[test] = clf.predict_proba(X_test)
multiclass_log_loss(y, y_pred)
The multiclass log loss function is an classification error metric that heavily penalizes you for being both confident (either predicting very high or very low class probability) and wrong. Throughout the competition you will want to check that your model improvements are driving this loss metric lower.
Now that you've made a simple metric, created a model, and examined the model's performance on the training data, the next step is to make improvements to your model to make it more competitive. The random forest model we created does not perform evenly across all classes and in some cases fails completely. By creating new features and looking at some of your distributions for the problem classes directly, you can identify features that specifically help separate those classes from the others. You can add new metrics by considering other image properties, stratified sampling, transformations, or other models for the classification.
header = "acantharia_protist_big_center,acantharia_protist_halo,acantharia_protist,amphipods,appendicularian_fritillaridae,appendicularian_s_shape,appendicularian_slight_curve,appendicularian_straight,artifacts_edge,artifacts,chaetognath_non_sagitta,chaetognath_other,chaetognath_sagitta,chordate_type1,copepod_calanoid_eggs,copepod_calanoid_eucalanus,copepod_calanoid_flatheads,copepod_calanoid_frillyAntennae,copepod_calanoid_large_side_antennatucked,copepod_calanoid_large,copepod_calanoid_octomoms,copepod_calanoid_small_longantennae,copepod_calanoid,copepod_cyclopoid_copilia,copepod_cyclopoid_oithona_eggs,copepod_cyclopoid_oithona,copepod_other,crustacean_other,ctenophore_cestid,ctenophore_cydippid_no_tentacles,ctenophore_cydippid_tentacles,ctenophore_lobate,decapods,detritus_blob,detritus_filamentous,detritus_other,diatom_chain_string,diatom_chain_tube,echinoderm_larva_pluteus_brittlestar,echinoderm_larva_pluteus_early,echinoderm_larva_pluteus_typeC,echinoderm_larva_pluteus_urchin,echinoderm_larva_seastar_bipinnaria,echinoderm_larva_seastar_brachiolaria,echinoderm_seacucumber_auricularia_larva,echinopluteus,ephyra,euphausiids_young,euphausiids,fecal_pellet,fish_larvae_deep_body,fish_larvae_leptocephali,fish_larvae_medium_body,fish_larvae_myctophids,fish_larvae_thin_body,fish_larvae_very_thin_body,heteropod,hydromedusae_aglaura,hydromedusae_bell_and_tentacles,hydromedusae_h15,hydromedusae_haliscera_small_sideview,hydromedusae_haliscera,hydromedusae_liriope,hydromedusae_narco_dark,hydromedusae_narco_young,hydromedusae_narcomedusae,hydromedusae_other,hydromedusae_partial_dark,hydromedusae_shapeA_sideview_small,hydromedusae_shapeA,hydromedusae_shapeB,hydromedusae_sideview_big,hydromedusae_solmaris,hydromedusae_solmundella,hydromedusae_typeD_bell_and_tentacles,hydromedusae_typeD,hydromedusae_typeE,hydromedusae_typeF,invertebrate_larvae_other_A,invertebrate_larvae_other_B,jellies_tentacles,polychaete,protist_dark_center,protist_fuzzy_olive,protist_noctiluca,protist_other,protist_star,pteropod_butterfly,pteropod_theco_dev_seq,pteropod_triangle,radiolarian_chain,radiolarian_colony,shrimp_caridean,shrimp_sergestidae,shrimp_zoea,shrimp-like_other,siphonophore_calycophoran_abylidae,siphonophore_calycophoran_rocketship_adult,siphonophore_calycophoran_rocketship_young,siphonophore_calycophoran_sphaeronectes_stem,siphonophore_calycophoran_sphaeronectes_young,siphonophore_calycophoran_sphaeronectes,siphonophore_other_parts,siphonophore_partial,siphonophore_physonect_young,siphonophore_physonect,stomatopod,tornaria_acorn_worm_larvae,trichodesmium_bowtie,trichodesmium_multiple,trichodesmium_puff,trichodesmium_tuft,trochophore_larvae,tunicate_doliolid_nurse,tunicate_doliolid,tunicate_partial,tunicate_salp_chains,tunicate_salp,unknown_blobs_and_smudges,unknown_sticks,unknown_unclassified".split(',')
labels = map(lambda s: s.split('/')[-1], namesClasses)
#get the total test images
fnames = glob.glob(os.path.join("test", "*.jpg"))
numberofTestImages = len(fnames)
numberofTestImages
X_test = np.zeros((numberofTestImages, num_features), dtype=float)
images = map(lambda fileName: fileName.split('/')[-1], fnames)
i = 0
# report progress for each 5% done
report = [int((j+1)*numberofTestImages/20.) for j in range(20)]
for fileName in fnames:
# Read in the images and create the features
image = imread(fileName, as_grey=True)
axisratio = getMinorMajorRatio(image)
image = resize(image, (maxPixel, maxPixel))
# Store the rescaled image pixels and the axis ratio
X_test[i, 0:imageSize] = np.reshape(image, (1, imageSize))
X_test[i, imageSize] = axisratio
i += 1
if i in report: print np.ceil(i *100.0 / numberofTestImages), "% done"
y_pred = clf.predict_proba(X_test)
y_pred.shape
df = pandas.DataFrame(y_pred, columns=labels, index=images)
df.index.name = 'image'
df = df[header]
df.to_csv('submissions/submission.csv')
!gzip submissions/submission.csv
!ls -l submissions/submission.csv.gz
stop_time = datetime.datetime.now()
print stop_time
print (stop_time - start_time), "elapsed time"