Predicting Citrus varieties from Reflectance Data

This is yet another exercise from the Caltech/JPL Summer School on Big Data Analytics. This one is for the lectures by David Thompson on Dimensionality Reduction. The problem is to differentiate citrus plants from non-citrus plants (and then between different varieties of citrus plants) from the light spectrum they reflect as observed by an airborne instrument. The exercise requires the student to do some basic data analysis, then perform Principal Component Analysis (PCA) to find the best features to differentiate between Citrus and Non-Citrus plants. In the next part of the exercise, the student is asked to use Linear Discriminant Analysis (LDA) to predict one of six varieties of Citrus plants.

The dataset consists of 4 comma-separated (CSV) text files. The first two are for the first (PCA) exercise, and the second two are for the second (LDA) exercise. Each pair consists of a data file with the class (citrus or non-citrus in the first pair and the citrus species in the second case), the plant number and a set of reflectance numbers, each corresponding to a different wavelength. There are 426 wavelengths in the first case and 388 in the second. Each plant has multiple records, possibly corresponding to observations at different times. The actual wavelength values are supplied in the second file (for each pair), one wavelength per line. The first dataset has 5,944 records and the second dataset has 2,109 records.

In this post I describe my attempt at solving the exercise. I used Python with Numpy and Scikit-Learn, and Matplotlib for the visualizations.

We first attempt to visualize the data and see if we can observe a difference in the reflectance properties between citrus and non-citrus plants. This is

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

from sklearn.cross_validation import KFold from sklearn.svm import SVC from sklearn.metrics import accuracy_score from sklearn.preprocessing import StandardScaler import matplotlib.pyplot as plt import numpy as np import os import pandas as pd import random DATA_DIR = "../../data/citrus_pca" def plot_by_class(X, y, w): """Plot reflectance across wavelengths for citrus vs non-citrus""" viz_df = pd.DataFrame() viz_df["wavelength"] = w colnames = {1: "citrus", 2: "noncit"} for yval in range(1,3): Xc = X[np.where(y == yval)[0], :] xmin = np.min(Xc, axis=0) xavg = np.mean(Xc, axis=0) xmax = np.max(Xc, axis=0) viz_df[colnames[yval] + "_min"] = xmin.T viz_df[colnames[yval] + "_avg"] = xavg.T viz_df[colnames[yval] + "_max"] = xmax.T plt.subplot(211) viz_df.plot(x="wavelength", y=["citrus_min", "citrus_avg", "citrus_max"], title="Citrus Plants", style=["r--", "b", "g--"]) plt.xticks([]) plt.xlabel("") plt.ylabel("reflectance") plt.subplot(212) viz_df.plot(x="wavelength", y=["noncit_min", "noncit_avg", "noncit_max"], title="Non-Citrus Plants", style=["r--", "b", "g--"]) plt.ylabel("reflectance") plt.show() cid_df = pd.read_csv(os.path.join(DATA_DIR, "Citrus_Identification_Data.txt")) y_class = np.asarray(cid_df["class"]) y_plant = np.asarray(cid_df["plant"]) X = np.matrix(cid_df[cid_df.columns[3:]]) # capture the bNNN cols (representing wavelengths) wavelengths = np.loadtxt(os.path.join(DATA_DIR, "Citrus_Identification_Wavelengths.txt")) # How many classes? print "#-classes:", len(np.unique(y_class)) print "#-plants:", len(np.unique(y_plant)) # Plot reflectance vs wavelength for citrus vs non-citrus plants plot_by_class(X, y_class, wavelengths)

The code above tells us that the data describes two classes (Citrus, Non-Citrus) and 160 unique plants. The plot shown below shows the variation in reflectance for each class of plant for different wavelengths. The dashed green and red lines represent the maximum and minimum reflectance for all plants in that class, and the blue line represents the mean. As you can see, the graphs look quite similar, although there are small differences in the reflectance.

Our next attempt is to see if there are differences across individual plants. The following code calculates the mean spectrum for each of 10 randomly chosen plants. Once more, differences appear to be in the amount of reflectance of each plant, the shape of the spectrum looks fairly uniform across plants.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

def plot_random_plants(X, y, w, n): """Plot reflectance vs wavelengths for n random plants""" viz_df = pd.DataFrame() viz_df["wavelength"] = w # generate n random plant ids plant_ids = set() num_unique = len(np.unique(y)) while len(plant_ids) < n: plant_ids.add(int(num_unique * random.random())) # for each plant plot mean value of wavelengths plant_ids_lst = list(plant_ids) pcols = [] for i in range(n): Xp = X[np.where(y == plant_ids_lst[i])[0], :] xavg = np.mean(Xp, axis=0) pcols.append("p" + str(i)) viz_df["p" + str(i)] = xavg.T viz_df.plot(x="wavelength", y=pcols, title="Mean for %d random plants" % (n)) plt.ylabel("reflectance") plt.show() ## Plot Reflectace vs Wavelengths for 10 random plants on same chart plot_random_plants(X, y_plant, wavelengths, 10)

We then attempt to calculate the correlation between features. As this blog post points out, high correlation between features may be an indication that there is one or more underlying factors that drives their similarity, and some less important factors drive their differences. When features are highly correlated, it may be possible to describe them adequately with a small number of principal components.

The code below computes an all-pairs correlations across all pairs of features, and plots the correlation matrix as a heatmap. I scaled the data first (mean of 0 and standard deviation of 1) first because PCA (next step) needs the input to be scaled. This does not have an effect on the correlation coefficients. As you can see from the heatmap, there is some correlation around the lower end of the spectrum and once again some (weaker) correlation at the higher end.

1 2 3 4 5 6 7 8 9 10 11 12 13 14

def plot_correlations_as_heatmap(X): """Plots heatmap from blue (low) to pink (high)""" corr = np.corrcoef(X) plt.imshow(corr, interpolation="nearest", cmap=plt.cm.cool) plt.xticks(()) plt.yticks(()) plt.title("Correlation Matrix") plt.colorbar() plt.show() # plot heat map of correlations between features scaler = StandardScaler() Xscaled = scaler.fit_transform(X) plot_correlations_as_heatmap(Xscaled)

We then perform PCA on the scaled dataset and see how the eigenvalues decay across the number of components. We restrict our plot to the first 50 principal components (to keep the chart readable), which explains 85.6% of the total variance. The first few components account for the bulk of the variance, so the behavior of the system can be approximated by considering only the first few components. Since the correlations are not very strong, the decay is not that dramatic.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

def plot_fraction_of_variance_explained(s, n): """Plot the first n values of eigenvalue s as fraction of total""" fve = s / np.cumsum(s)[-1] plt.plot(np.arange(len(fve))[0:n], fve[0:n]) plt.title("Eigenvalue Decay") plt.xlabel("Principal Eigenvalues") plt.ylabel("Fraction of Variance Explained") plt.show() return fve # Perform a PCA on the (scaled) matrix U, s, V = np.linalg.svd(Xscaled) # Plot eigenvalue decay (as fraction of variance explained) fve = plot_fraction_of_variance_explained(s, 50)

We then project our data to the first 10 components (accounting for 69.4% of the total variance) and plot the correlations between the first few components. We do this for both the original data and the projected data. As expected, the original data exhibit higher correlation than the projected data. However, even in the projected data, there is no clear separation between citrus (blue) and non-citrus (red).

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

def project_data(X, U, k): Ured = U[:, 0:k] Sred = np.diag(s[0:k]) return np.dot(Ured, Sred) def recover_data(Xred, V, k): Vred = V[0:k, :] return np.dot(Xred, Vred) def plot_scatter_plots(X, y, title): pairs = [(0, 1), (1, 2), (0, 2)] Xparts = [X[np.where(y == 1)[0], :], X[np.where(y == 2)[0], :]] for i in range(len(pairs)): plt.subplot(1, 3, i + 1) for j in range(len(Xparts)): ecol = 'b' if j == 0 else 'r' plt.scatter(Xparts[j][:, pairs[i][0]], Xparts[j][:, pairs[i][1]], s=10, edgecolors=ecol, facecolors='none') if i == 1: plt.title(title) plt.xlabel("X" + str(pairs[i][0])) plt.ylabel("X" + str(pairs[i][1])) plt.xticks(()) plt.yticks(()) i = i + 1 plt.show() # Project dataset onto the first k components k = 10 print "Variance explained with %d principal components = %0.3f" % (k, np.sum(fve[0:k])) # Project Xscaled to first k dimensions Xprojected = project_data(Xscaled, U, k) plot_by_class(Xprojected, y_class, np.arange(k)) # Scatter plot between 1/2. 2/3, 1/3 plot_scatter_plots(Xscaled, y_class, "Original Data") plot_scatter_plots(Xprojected, y_class, "Projected Data")

The final step in this part of the exercise was to build a classifier to predict between citrus and non-citrus plants. Because it appeared that the reflectance of citrus plants are spread out higher than non-citrus plants (also borne out by the original chart of citrus vs non-citrus reflectances), I decided to use a Support Vector Machine with a Radial (RBF) kernel. The results were a bit counter-intuitive - with 3-fold cross-validation over the original data, the accuracy was 85.1% but after PCA only 80.1%. I also tried using a RandomForest with 200 trees (code not shown here), and obtained an accuracy of 84.3% with the original data and 85.7% after PCA (data projected to the first 10 components).

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

def evaluate_classification(X, y): """Evaluate classifier using 10 fold cross validation""" kfold = KFold(X.shape[0], 3) scores = [] for train, test in kfold: Xtrain, Xtest, ytrain, ytest = X[train], X[test], y_class[train], y_class[test] clf = SVC() # rbf kernel default clf.fit(Xtrain, ytrain) ypred = clf.predict(Xtest) scores.append(accuracy_score(ypred, ytest)) return 1.0 * sum(scores) / len(scores) # Train classifier on original data with 10-fold cross validation score_orig = evaluate_classification(Xscaled, y_class) score_pca = evaluate_classification(Xprojected, y_class) print "Classification accuracy: original: %.3f, after PCA: %.3f" % (score_orig, score_pca)

The next part of the exercise was to build an LDA classifier to differentiate between different species (varieties) of citrus plants, again based on similar reflectance data. The first step is to do PCA to project the data onto its two principal components. As you can see, there is some separability between the different varieties.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39

from sklearn.cross_validation import train_test_split, KFold from sklearn.lda import LDA from sklearn.metrics import accuracy_score from sklearn.preprocessing import StandardScaler import matplotlib.pyplot as plt import numpy as np import os import pandas as pd DATA_DIR = "../../data/citrus_pca" cvdf = pd.read_csv(os.path.join(DATA_DIR, "Citrus_Varieties.txt")) y_species = np.asarray(cvdf["class"]) y_plant = np.asarray(cvdf["plant"]) X = np.matrix(cvdf[cvdf.columns[3:]]) wavelengths = np.loadtxt(os.path.join(DATA_DIR, "Citrus_Varieties_Wavelengths.txt")) # Reduce to 2 dimensions and visualize scaler = StandardScaler() Xscaled = scaler.fit_transform(X) U, s, V = np.linalg.svd(Xscaled) Ured = U[:, 0:2] Sred = np.diag(s[0:2]) Xproj = np.dot(Ured, Sred) species_ids = np.unique(y_species) print "species_ids:", species_ids colors = ['b', 'g', 'r', 'c', 'm', 'y'] i = 0 for species_id in species_ids: Xpart = Xproj[np.where(y_species == species_id)[0], :] plt.scatter(Xpart[:, 0], Xpart[:, 1], color=colors[i]) i = i + 1 plt.title("Citrus Species (first 2 Principal Components)") plt.xlabel("X0") plt.ylabel("X1") plt.show()

The next step is to train a multiclass LDA classifier and test its accuracy on a holdout test dataset. We split our dataset 75/25 train/test, train a Scikit-Learn LDA classifier on the training data, and obtain a pretty impressive score of 97.3% accuracy on the hold out test set. LDA provides the "mean" vector for each of its classes. In order to check which species are most similar, we build a correlation matrix and represent it as a heatmap. As can be easily seen from the heatmap, the most similar species (in terms of reflectance spectrum) are the Mandarin, Pummelo and the Hybrid.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

# Perform multiclass LDA Xtrain, Xtest, ytrain, ytest = train_test_split(Xscaled, y_species, test_size=0.25, random_state=42) clf = LDA(len(species_ids)) clf.fit(Xtrain, ytrain) ypred = clf.predict(Xtest) print "LDA Accuracy Score: %.3f" % (accuracy_score(ypred, ytest)) # What varieties are most spectrally similar? corr = np.corrcoef(clf.means_) plt.imshow(corr, interpolation="nearest", cmap=plt.cm.cool) plt.xticks(np.arange(len(species_ids)), species_ids, rotation=45) plt.yticks(np.arange(len(species_ids)), species_ids) plt.colorbar() plt.show()

Finally, in an effort to test if we can predict species (and then plant) effectively from the reflectance data, we build small harnesses to do 10-fold cross validation with the LDA classifier. The cross-validation accuracy score for the species predictor is a slightly less impressive 78.6% and for the plant predictor an abysmal 3.2%.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

# Find LDA classifier accuracy using cross validation kfold = KFold(Xscaled.shape[0], 10) scores = [] for train, test in kfold: Xtrain, Xtest, ytrain, ytest = Xscaled[train], Xscaled[test], \ y_species[train], y_species[test] clf = LDA(len(species_ids)) clf.fit(Xtrain, ytrain) ypred = clf.predict(Xtest) scores.append(accuracy_score(ypred, ytest)) print "LDA accuracy predicting species: %.3f" % (1.0 * sum(scores) / len(scores)) # Find LDA classifier accuracy for predicting specific plants kfold = KFold(Xscaled.shape[0], 10) scores = [] for train, test in kfold: Xtrain, Xtest, ytrain, ytest = Xscaled[train], Xscaled[test], \ y_plant[train], y_plant[test] clf = LDA(len(species_ids)) clf.fit(Xtrain, ytrain) ypred = clf.predict(Xtest) scores.append(accuracy_score(ypred, ytest)) print "LDA accuracy predicting plants: %.3f" % (1.0 * sum(scores) / len(scores))

This is all I have for today. As always (at least for most of my recent posts), this code is available on github here. I had a lot of fun doing this exercise even though the results are probably not that great. I did not know about the LDA classifier before this, got to learn about this. I also learned quite a few Matplotlib tricks. I hope you enjoyed it too.

IR Math in Java : Cluster Visualization

I've been trying to learn clustering algorithms lately. I was planning to write about them this week, but some last minute refactoring to remove redundancies and make the code more readable resulted in everything going to hell. So I guess I will have to write about them next week.

Almost all clustering algorithms (at least the ones I have seen) seem to be non-deterministic, mainly because they select documents randomly from the collection to build the initial clusters. As a result, they can come up with wildly different clusters depending on how the initial clusters were formed. In my previous (un-refactored) code, for example, the K-Means algorithm converged to the same set of clusters most of the time, but with the changes, they no longer do.

Working through this for some time, I decided I needed to see for myself what the "correct" clusters were. So if I could visualize the documents as points in the n-dimensional term space, clumps of points would correspond to clusters. The problem was that I had only 2 (or maximum 3) dimensions of visualization to work with.

Luckily for me, smarter people than I have faced and solved the same problem, and they have been kind enough to write about it on the web. The solution is to do Dimensionality Reduction, extracting from the term-document matrix the first 2 or 3 most interesting components (or Principal Components) and use them as the values for a 2-dimensional or 3-dimensional scatter chart.

The mathematical background for Principal Component Analysis (PCA) is explained very nicely in this tutorial, which I quote verbatim below.

The mathematical technique used in PCA is called eigen analysis: we solve for the eigenvalues and eigenvectors of a square symmetric matrix with sums of squares and cross products. The eigenvector associated with the largest eigenvalue has the same direction as the first principal component. The eigenvector associated with the second largest eigenvalue determines the direction of the second principal component. The sum of the eigenvalues equals the trace of the square matrix and the maximum number of eigenvectors equals the number of rows (or columns) of this matrix.

It then goes on to explain the algorithm that should be used for reducing and extracting the most interesting dimensions (See Section 6, Algorithms) for a non-square matrix such as our term-document matrix. Essentially, it decomposes the term-document matrix A using Singular Value Decomposition (SVD) into 3 matrices, U, S and V, where the following equation holds.

  A = U * S * VT

Here S is a square diagonal matrix, where values are in descending order down the diagonal. So for a 2 dimensional reduction, the principal components correspond to the first 2 columns of V, and for a 3 dimensional reduction, the principal components correspond to the first 3 columns of V.

The Java code to generate data for drawing the charts is trivial, mainly because we use the Jama matrix library, which does all the heavy lifting for SVD calculations.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39

// Source: src/main/java/com/mycompany/myapp/clustering/PcaClusterVisualizer.java package com.mycompany.myapp.clustering; import java.io.FileWriter; import java.io.IOException; import java.io.PrintWriter; import Jama.Matrix; import Jama.SingularValueDecomposition; public class PcaClusterVisualizer { private final String PLOT_2D_OUTPUT = "plot2d.dat"; private final String PLOT_3D_OUTPUT = "plot3d.dat"; public void reduce(Matrix tdMatrix, String[] docNames) throws IOException { PrintWriter plot2dWriter = new PrintWriter(new FileWriter(PLOT_2D_OUTPUT)); PrintWriter plot3dWriter = new PrintWriter(new FileWriter(PLOT_3D_OUTPUT)); SingularValueDecomposition svd = new SingularValueDecomposition(tdMatrix); Matrix v = svd.getV(); // we know that the diagonal of S is ordered, so we can take the // first 3 cols from V, for use in plot2d and plot3d Matrix vRed = v.getMatrix(0, v.getRowDimension() - 1, 0, 2); for (int i = 0; i < v.getRowDimension(); i++) { // rows plot2dWriter.printf("%6.4f %6.4f %s%n", Math.abs(vRed.get(i, 0)), Math.abs(vRed.get(i, 1)), docNames[i]); plot3dWriter.printf("%6.4f %6.4f %6.4f %s%n", Math.abs(vRed.get(i, 0)), Math.abs(vRed.get(i, 1)), Math.abs(vRed.get(i, 2)), docNames[i]); } plot2dWriter.flush(); plot3dWriter.flush(); plot2dWriter.close(); plot3dWriter.close(); } }

The term-document matrix is generated from my 7 document title collection that I have been using for my experiments, using the following snippet of JUnit code. See one of my earlier posts titled IR Math in Java : TF, IDF and LSI for the actual data and details on the classes being used.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

@Test public void testPcaClusterVisualization() throws Exception { // for brevity, this block is in a @Before method in the actual // code, it has been globbed together here for readability VectorGenerator vectorGenerator = new VectorGenerator(); vectorGenerator.setDataSource(new DriverManagerDataSource( "com.mysql.jdbc.Driver", "jdbc:mysql://localhost:3306/tmdb", "irstuff", "irstuff")); Map<String,Reader> documents = new LinkedHashMap<String,Reader>(); BufferedReader reader = new BufferedReader( new FileReader("src/test/resources/data/indexing_sample_data.txt")); String line = null; while ((line = reader.readLine()) != null) { String[] docTitleParts = StringUtils.split(line, ";"); documents.put(docTitleParts[0], new StringReader(docTitleParts[1])); } vectorGenerator.generateVector(documents); IdfIndexer indexer = new IdfIndexer(); tdMatrix = indexer.transform(vectorGenerator.getMatrix()); documentNames = vectorGenerator.getDocumentNames(); documentCollection = new DocumentCollection(tdMatrix, documentNames); // this is my actual @Test block PCAClusterVisualizer visualizer = new PCAClusterVisualizer(); visualizer.reduce(tdMatrix, documentNames); }

This generates 2 data files which are used as inputs to gnuplot to generate 2D and 3D scatter charts. The data, chart, and the gnuplot code to generate the chart is shown in the table below:

1 2 3 4 5 6 7 8 # plot2d.dat 0.0000 0.2261 D1 0.0468 0.0000 D2 0.0000 0.7363 D3 0.0000 0.6378 D4 0.0000 0.0000 D5 0.8751 0.0000 D6 0.4817 0.0000 D7	1 2 3 4 5 6 7 8 # plot3d.dat 0.0000 0.2261 0.0000 D1 0.0468 0.0000 0.2997 D2 0.0000 0.7363 0.0000 D3 0.0000 0.6378 0.0000 D4 0.0000 0.0000 0.0000 D5 0.8751 0.0000 0.4723 D6 0.4817 0.0000 0.8289 D7
1 2 3 4 5 6 # plot2d.gp set style data labels unset key plot 'plot2d.dat' using 1:2:3 \ with labels font "arial,11" \ textcolor lt 1	1 2 3 4 5 6 # plot3d.gp set style data labels unset key splot 'plot3d.dat' using 1:2:3:4 \ with labels font "arial,11" \ textcolor lt 1

From the charts above, it appears that the following clusters may be valid for our test document set. Notice that although D3 and D7 appear really close (overlapped) on the 3D chart, they don't seem to be close going by the 2D chart or the data. In any case, the results look believable, although not perfect, but that could be due to dimensionality reduction and/or the small data set.

1 2 3 4 5 6 7 8 9 10 11

C0: [D1, D2, D5] D1 Human machine interface for <b>computer</b> applications D2 A survey of user opinion of <b>computer</b> system response time D5 The generation of random, binary and ordered trees C1: [D3, D4] D3 The <b>EPS</b> user interface management <b>system</b> D4 <b>System</b> and human system engineering testing of <b>EPS</b> C2: [D7] D7 Graph minors: A survey C3: [D6] D6 The intersection graph of paths in trees

I think this post may be helpful to programmers like me who are just getting into IR (most people who are heavily into IR would probably know this stuff already). Text mining algorithms, by their very nature, need to deal with n-dimensional data, and the ability to visualize the data in 2D or 3D can be quite enlightening, so this is a useful tool to have in one's text mining toolbox.

Update 2009-04-26: In recent posts, I have been building on code written and described in previous posts, so there were (and rightly so) quite a few requests for the code. So I've created a project on Sourceforge to host the code. You will find the complete source code built so far in the project's SVN repository.