Galaxy Image Classification¶
Importing Libraries¶
In [1]:
import numpy as np
from sklearn.model_selection import train_test_split
from keras.preprocessing.image import ImageDataGenerator
from tqdm import tqdm #displays progress bar within for loops
import matplotlib.pyplot as plt
from skimage import io
import pickle
from skimage.transform import resize
import glob
import os
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')
from keras import regularizers, optimizers
from functions.model_functions import plot_confusion_matrix, predict_one_image_cnn
import tensorflow as tf
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3'
Size of all images¶
size = set()
for filename in glob.glob('data/*.jpg'):
im=io.imread(filename).shape
size.add(im)
size for all images are: {(424, 424, 3)}
Visualizing some of the training data¶
In [2]:
#Let's look at four of our galaxies from the 61578 images in our training set.
plt.figure(1, figsize=(12,9))
plt.subplot(231)
plt.imshow(plt.imread('data/train/spiral/100023.jpg'))
plt.title('Galaxy 100023')
plt.subplot(232)
plt.imshow(plt.imread('data/train/elliptical/100078.jpg'))
plt.title('Galaxy 100078')
plt.subplot(233)
plt.imshow(plt.imread('data/train/elliptical/100123.jpg'))
plt.title('Galaxy 100123')
plt.subplot(234)
plt.imshow(plt.imread('data/train/spiral/100143.jpg'))
plt.title('Galaxy 100143')
plt.subplot(235)
plt.imshow(plt.imread('data/train/spiral/233622.jpg'))
plt.title('Galaxy 233622')
plt.subplot(236)
plt.imshow(plt.imread('data/train/spiral/177755.jpg'))
plt.title('Galaxy 177755')
plt.show()
plt.tight_layout()
From the images above, as a scientist I can already tell what morphology these galaxies are of, but now the question is will the computer be able to classify them correctly. Let's first build a baseline model then compare it to a CNN model.
Predictions.csv¶
In [3]:
import pandas as pd
import numpy as np
predictions = pd.read_csv('predictions.csv')
predictions.head()
Out[3]:
There seems to be ambigious headings for our columns and the resources on Kaggle don't seem to help much. Essentially, the classes refer to the morphology on the galaxy.
In [4]:
#All of the classes in our Class1 add up to 1. These classes refer to the probaility of
#the galaxy being either the shape in Class1.1, Class1.2, or in Class1.3.
(predictions['Class1.1'] + predictions['Class1.2'] + predictions['Class1.3']).head()
Out[4]:
In [5]:
#Dropping all the other columns which do not refer to the morphology.
predictions = predictions[['GalaxyID','Class1.1', 'Class1.2', 'Class1.3']]
predictions.head(10)
Out[5]:
Based on the images printed above, we can use our physics brain to determine what the columns stand for and rename the columns.
In [6]:
predictions.columns = ['GalaxyID', 'Elliptical', 'Spiral', 'Irregular']
predictions.head()
Out[6]:
If the value in the irregular was the max value for that row, we called it an irregular.
Note: The following code below no longer works because the irregular folders are deleted.
print('num in train elliptical = {}'.format(len(os.listdir('data/train/elliptical/'))))
print('num in train spiral = {}'.format(len(os.listdir('data/train/spiral/'))))
print('num in train irregular = {}'.format(len(os.listdir('data/train/irregular/'))))
print('num in test elliptical = {}'.format(len(os.listdir('data/test/elliptical/'))))
print('num in test spiral = {}'.format(len(os.listdir('data/test/spiral/'))))
print('num in test irregular = {}'.format(len(os.listdir('data/test/irregular/'))))
num in train elliptical = 13237
num in train spiral = 17526
num in train irregular = 26
num in test elliptical = 13456
num in test spiral = 17300
num in test irregular = 33
For this, we can remove the irregular images since they make such a small portion of our over all dataset.
Now, let's delete the irregular galaxies from the predictions (because we deleted their images). Then we will drop the "Irregular" column.
In [7]:
len(predictions[(predictions.Irregular > predictions.Spiral) & (predictions.Irregular > predictions.Elliptical)])
Out[7]:
In [8]:
predictions[(predictions.Irregular > predictions.Spiral) & (predictions.Irregular > predictions.Elliptical)].head()
Out[8]:
In [9]:
#Dropped all the irregular galaxies.
predictions.drop(predictions[(predictions.Irregular > predictions.Spiral) & (predictions.Irregular > predictions.Elliptical)].index,inplace=True)
In [10]:
predictions.drop(['Irregular'], axis=1, inplace=True)
In [11]:
predictions.head()
Out[11]:
In [2]:
import matplotlib.pyplot as plt
import numpy as np
import os, shutil
from keras import models
from keras import layers
from keras.layers import Dropout
from sklearn.metrics import confusion_matrix, f1_score
np.random.seed(123)
from keras.models import load_model
from keras.preprocessing.image import ImageDataGenerator, array_to_img, img_to_array, load_img
In [3]:
# get all the data in the directory split/test, and reshape them
datagen = ImageDataGenerator(rescale=1./255,
rotation_range=90,
brightness_range=(0.5,2),
)
#Take 30000 images from the training folder and resize images down to 106 x 106 px
data_tr=datagen.flow_from_directory(
'data/train',
target_size=(106, 106), #actual image size
batch_size = 30000,
class_mode='binary',
seed = 123)
In [4]:
#Take 20000 images from the testing folder and resize images down to 106 x 106 px
data_te = ImageDataGenerator(rescale=1./255).flow_from_directory(
'data/test',
target_size=(106, 106),
batch_size = 20000,
class_mode='binary',
seed = 123)
In [5]:
# Defining x_train, y_train, and x_test, y_test by grabbing the first 30,000 batch for train
# and 20,000 image batch for test
x_tr, y_tr = next(data_tr)
x_te, y_te = next(data_te)
In [6]:
#Displaying that our training data is binary
y_tr[:5]
Out[6]:
In [7]:
#Splitting our data using sklearn. 80% of our data is X_train, y_train - what we are training
#our model on. 20% of our data is the X_val, y_val what we are testing our data on.
from sklearn.model_selection import train_test_split
X_train, X_val, y_train, y_val = train_test_split(x_tr, y_tr, test_size=0.20, random_state=123)
In [ ]:
len(x_te)
In [ ]:
X_val.shape
CNN Architecture¶
In [8]:
cnn = models.Sequential()
cnn.add(layers.Conv2D(64, (3, 3), kernel_regularizer=regularizers.l2(0.1), input_shape=(106, 106, 3), padding='valid',strides=(2,2)))
cnn.add(layers.Conv2D(64, (3, 3), input_shape=(106, 106, 3), padding='valid',strides=(2,2)))
cnn.add(layers.BatchNormalization())
cnn.add(layers.Activation('relu'))
cnn.add(layers.MaxPooling2D((2, 2)))
cnn.add(Dropout(0.5))
cnn.add(layers.Conv2D(64, (3, 3), kernel_regularizer=regularizers.l2(0.1)))
cnn.add(layers.BatchNormalization())
cnn.add(layers.Activation('relu'))
cnn.add(layers.MaxPooling2D((2, 2)))
cnn.add(layers.Flatten())
cnn.add(layers.Dense(128))
cnn.add(layers.BatchNormalization())
cnn.add(layers.Activation('relu'))
cnn.add(layers.Dense(32, activation='relu'))
cnn.add(layers.Dense(1, activation='sigmoid'))
# decay_rate = learning_rate / epochs
sgd = optimizers.SGD(lr=0.01, decay=0.0002, momentum=0.9, nesterov=False)
cnn.compile(loss='binary_crossentropy',
optimizer=sgd,
metrics=['acc'])
cnn.load_weights('model_weights/cnn11_only_lr.h5')
Fitting CNN Model¶
In [ ]:
cnn1 = cnn.fit(X_train, y_train,
epochs=50,
validation_data=(X_val, y_val),
batch_size=500)
In [ ]:
print(cnn.summary())
Let's see what our CNN is doing behind the scenes when it is fed this image:¶
In [14]:
# Printing one random image to compare with the hidden layers below.
plt.imshow(X_train[50])
plt.show()
Grabbing the first hidden 8 layers from CNN to visualize the activation from the first layer¶
In [65]:
# layer_outputs = [layer.output for layer in cnn.layers[:8]]
In [9]:
layer_outputs = [layer.output for layer in cnn.layers]
In [20]:
cnn.layers
Out[20]:
In [10]:
activation_model = models.Model(inputs=cnn.input, outputs=layer_outputs)
In [11]:
from keras import models
In [41]:
# from keras import models
from keras.preprocessing import image
activations = activation_model.predict(X_train)
fig, axes = plt.subplots(5, 4, figsize=(12,24))
for i in range(20):
row = i//4
column = i%4
ax = axes[row, column]
first_layer_activation = activations[1]
ax.matshow(first_layer_activation[50, :, :, i], cmap='viridis')
In [17]:
# from keras import models
from keras.preprocessing import image
activations = activation_model.predict(X_train)
fig, axes = plt.subplots(5, 4, figsize=(12,24))
for i in range(20):
row = i//4
column = i%4
ax = axes[row, column]
first_layer_activation = activations[7]
ax.matshow(first_layer_activation[50, :, :, i], cmap='viridis')
In [18]:
# from keras import models
from keras.preprocessing import image
activations = activation_model.predict(X_train)
fig, axes = plt.subplots(5, 4, figsize=(12,24))
for i in range(20):
row = i//4
column = i%4
ax = axes[row, column]
first_layer_activation = activations[0]
ax.matshow(first_layer_activation[50, :, :, i], cmap='viridis')
In [19]:
# from keras import models
from keras.preprocessing import image
activations = activation_model.predict(X_train)
fig, axes = plt.subplots(5, 4, figsize=(12,24))
for i in range(20):
row = i//4
column = i%4
ax = axes[row, column]
first_layer_activation = activations[5]
ax.matshow(first_layer_activation[50, :, :, i], cmap='viridis')
In [20]:
# from keras import models
from keras.preprocessing import image
activations = activation_model.predict(X_train)
fig, axes = plt.subplots(5, 4, figsize=(12,24))
for i in range(20):
row = i//4
column = i%4
ax = axes[row, column]
first_layer_activation = activations[9]
ax.matshow(first_layer_activation[50, :, :, i], cmap='viridis')
In [12]:
# from keras import models
from keras.preprocessing import image
activations = activation_model.predict(X_train)
fig, axes = plt.subplots(5, 4, figsize=(12,24))
for i in range(20):
row = i//4
column = i%4
ax = axes[row, column]
first_layer_activation = activations[3]
ax.matshow(first_layer_activation[50, :, :, i], cmap='viridis')
In [13]:
# from keras import models
from keras.preprocessing import image
activations = activation_model.predict(X_train)
fig, axes = plt.subplots(5, 4, figsize=(12,24))
for i in range(20):
row = i//4
column = i%4
ax = axes[row, column]
first_layer_activation = activations[4]
ax.matshow(first_layer_activation[50, :, :, i], cmap='viridis')
Confusion Matrix for CNN¶
In [69]:
cnn.evaluate(x_te, y_te)
Out[69]:
In [70]:
predictions_cnn = cnn.predict(x_te)
predictions_cnn = np.around(predictions_cnn)
In [71]:
f1_score(y_te, predictions_cnn, average='macro')
Out[71]:
In [72]:
plt.figure()
plot_confusion_matrix(confusion_matrix(y_te, predictions_cnn), classes=['Elliptical', 'Spiral'])
plt.title('Confusion Matrix for CNN Model')
plt.savefig('CNN_ConfusionMatrix.png')
Let's test our model on some images of galaxies it has not seen yet.¶
In [ ]:
def predictoneimage_cnn(cnn, path):
img = load_img(path, target_size=(106, 106))
plt.imshow(img)
img = img_to_array(img)
img = img/255
img = np.expand_dims(img, axis=0)
predict = cnn.predict(img)
return predict
Keras' flow_from_directory method orders folder structure alphabetical, elliptical is first, and 0 is first. So, 0 is elliptical and 1 is spiral.¶
In [73]:
predict_one_image_cnn(cnn, 'data/train/elliptical/100078.jpg')
Out[73]:
In [74]:
predict_one_image_cnn(cnn, 'data/test/elliptical/564639.jpg')
Out[74]:
The CNN model classified the above galaxy as an elliptical. Which is correct since this test image came from the elliptical folder.¶
ROC Curve for CNN Model¶
In [75]:
from sklearn.metrics import roc_curve
In [76]:
fpr_cnn, tpr_cnn, thresholds = roc_curve(y_te, cnn.predict(x_te))
plt.plot(fpr_cnn,tpr_cnn)
plt.plot([0, 1], [0, 1],'r--')
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
plt.title('ROC Curve for CNN Model')
# plt.show()
plt.savefig('CNN_ROC.png')
In [77]:
from sklearn.metrics import auc
roc_auc = auc(fpr_cnn, tpr_cnn)
In [78]:
print(roc_auc)
Decision Tree¶
In [18]:
from sklearn.tree import DecisionTreeClassifier
from sklearn.metrics import accuracy_score
tree_model = DecisionTreeClassifier(max_depth=10)
In [19]:
x_tr.shape
Out[19]:
In [20]:
X_train.shape
Out[20]:
In [21]:
y_train.shape
Out[21]:
In [22]:
x_tr_flat = X_train.reshape((24000,X_train.shape[1]*X_train.shape[2]*X_train.shape[3]))
In [23]:
tree_model.fit(x_tr_flat, y_train)
Out[23]:
In [24]:
x_te.shape
Out[24]:
In [25]:
x_te_flat = x_te.reshape((20000,x_te.shape[1]*x_te.shape[2]*x_te.shape[3]))
In [33]:
y_predict = tree_model.predict(x_te_flat)
In [27]:
y_predict_proba = tree_model.predict_proba(x_te_flat)
In [28]:
y_predict_proba
Out[28]:
In [34]:
accuracy_score(y_te, y_predict)
Out[34]:
In [35]:
from sklearn.metrics import confusion_matrix
confusion_matrix(y_te, y_predict)
Out[35]:
In [36]:
plt.figure()
plot_confusion_matrix(confusion_matrix(y_te, y_predict), classes=['Elliptical', 'Spiral'])
plt.title('Confusion Matrix for Decision Tree')
plt.savefig('DecisionTree_ConfusionMatrix.png')
In [51]:
y_te[:10]
Out[51]:
In [58]:
# np.unique(y_predict_proba)[:5]
y_predict_proba[:10]
Out[58]:
In [59]:
[x[1] for x in y_predict_proba][:10]
Out[59]:
In [105]:
from sklearn.metrics import auc, roc_curve
fpr_decision, tpr_decision, thresholds = roc_curve(y_te, [x[1] for x in y_predict_proba])
plt.plot(fpr_decision,tpr_decision)
plt.plot([0, 1], [0, 1],'r--')
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
plt.title('ROC Curve for Decision Tree')
# plt.show()
plt.savefig('DecisionTree_ROC.png')
In [106]:
roc_auc = auc(fpr_decision, tpr_decision)
In [110]:
# fpr_decision
In [111]:
# tpr_decision
In [109]:
print(roc_auc)
Inception¶
In [79]:
#we can use any of the inception stuff even if they are trained on anything
from keras.applications import inception_v3
from keras.layers import Dense,GlobalAveragePooling2D
from sklearn.metrics import confusion_matrix, f1_score
from keras.models import Model
In [80]:
imagenet=inception_v3.InceptionV3(weights='imagenet',include_top=False)
imagenet_new=imagenet.output
new_model = models.Sequential()
new_model.add(imagenet)
new_model.add(GlobalAveragePooling2D())
new_model.add(Dense(1024,activation='relu'))
new_model.add(Dense(1024,activation='relu')) #dense layer 2
new_model.add(Dense(512,activation='relu')) #dense layer 3
new_model.add(Dense(1,activation='sigmoid')) #final layer
In [81]:
for i,layer in enumerate(imagenet.layers):
print(i,layer.name, layer.trainable)
In [82]:
#don't train the image net (or you will wait too long which is the first layer)
for layer in new_model.layers[:1]:
layer.trainable=False
In [83]:
new_model.compile(optimizer='Adam',loss='binary_crossentropy',metrics=['accuracy'])
# step_size_train=train_generator.n//train_generator.batch_size
new_model.fit(X_train,
y_train,
epochs=10,
batch_size=50,
validation_data=(X_val, y_val))
Out[83]:
In [84]:
new_model.evaluate(x_te, y_te)
Out[84]:
In [85]:
predictions_transfer = new_model.predict(x_te)
predictions_transfer = np.around(predictions_transfer)
In [86]:
f1_score(y_te, predictions_transfer, average='macro')
Out[86]:
In [87]:
plt.figure()
plot_confusion_matrix(confusion_matrix(y_te, predictions_transfer), classes=['Elliptical', 'Spiral'])
plt.title('Confusion Matrix for Inception')
plt.savefig('Inception_ConfusionMatrix.png')
In [112]:
fpr_inception, tpr_inception, thresholds = roc_curve(y_te, new_model.predict(x_te))
plt.plot(fpr_inception,tpr_inception)
plt.plot([0, 1], [0, 1],'r--')
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
plt.title('ROC Curve for Inception')
plt.savefig('Inception_ROC.png')
# plt.show()
In [ ]:
Simple Neural Network Baseline Model¶
In [89]:
from keras.models import Sequential, Input
from keras.layers import Dense
model = Sequential()
In [90]:
model.add(layers.Flatten(input_shape=(106, 106, 3)))
layer_1 = model.add(Dense(units=128, activation='relu'))
layer_2 = model.add(Dense(units=1, activation='sigmoid'))
In [91]:
model.compile(loss='binary_crossentropy',
optimizer='sgd',
metrics=['acc'])
In [92]:
history = model.fit(X_train, y_train,
epochs=50,
validation_data=(X_val, y_val),
batch_size=500)
In [93]:
# history.history['val_acc']
In [94]:
predictions_mlp = model.predict(x_te)
predictions_mlp = np.around(predictions_mlp)
In [132]:
plt.figure()
plot_confusion_matrix(confusion_matrix(y_te, predictions_mlp), classes=['Elliptical', 'Spiral'])
plt.title('Confusion Matrix for Neural Network')
plt.savefig('Baseline_MLP_ConfusionMatrix.png')
In [133]:
fpr_neural, tpr_neural, thresholds = roc_curve(y_te, model.predict(x_te))
plt.plot(fpr_neural,tpr_neural)
plt.plot([0, 1], [0, 1],'r--')
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
plt.title('ROC Curve for Neural Network')
plt.savefig('BaselineMLP_ROC.png')
In [97]:
fpr
Out[97]:
In [98]:
tpr
Out[98]:
In [99]:
plt.figure()
plt.plot(range(len(history.history['val_acc'])), history.history['val_acc'])
plt.xlabel('Epoch')
plt.ylabel('Validation Accuracy')
Out[99]:
In [100]:
loss_and_metrics = model.evaluate(x_te, y_te)
In [101]:
loss_and_metrics
Out[101]:
In [102]:
model.metrics_names
Out[102]:
In [103]:
model.save('simple_nn_baseline.h5')
In [113]:
# fpr_cnn, tpr_cnn, thresholds = roc_curve(y_te, cnn.predict(x_te)) # CNN -- already have
fpr_decision, tpr_decision, thresholds = roc_curve(y_te, [x[1] for x in y_predict_proba])
# fpr, tpr, thresholds = roc_curve(y_te, new_model.predict(x_te)) # inception
fpr_mlp, tpr_mlp, thresholds = roc_curve(y_te, model.predict(x_te)) # simple MLP
plt.plot(fpr_cnn, tpr_cnn, label='CNN')
plt.plot(fpr_inception, tpr_inception, label='Inception')
plt.plot(fpr_decision, tpr_decision, label="Decision Tree")
plt.plot(fpr_mlp, tpr_mlp, label="Baseline MLP")
plt.plot([0, 1], [0, 1],linestyle='--',c='black')
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
plt.title('ROC Curves')
plt.legend()
plt.savefig('all_models_roc.png')
In [ ]: