Load real data


In [21]:
import pandas as pd
import numpy as np
import os
#load pandas data (train set) - to be stored locally. Dataset available online in the repository (directory "data")
in_dir = "/home/acocac/post/tmp/"
df_train = pd.read_csv(in_dir + "traindata_2classes_byXYyear_label.csv", sep=',')

In [2]:
##input data preparation
#pandas to 3d array
train_dataset = np.array(list(df_train.groupby('Mid').apply(pd.DataFrame.as_matrix)))

#settings
n_classes = 2
target = ['red','blue','nir','mir','fyear'] #variables
target_label = ['0','1'] #variables in one hot format 

index_features = map(lambda col: df_train.columns.get_loc(col), target)
index_labels = map(lambda col: df_train.columns.get_loc(col), target_label)

#data
X = train_dataset[:,:,index_features]
y = train_dataset[:,:,index_labels]

Tensorflow - Paper approach


In [3]:
import tensorflow as tf
from tensorflow.contrib import rnn as rnn_cell
import io
from util.tf_utils import tf_confusion_metrics
import inspect
import util.eval as eval
import util.plots as plots

In [4]:
class Model():
    """
    Tensorflow Graph using Recurrent LSTM layers and fully connected softmax layer for field identification
    with multispectral/temporal data acquired from satellite imagery
    Params
        tf placeholders:
            X           Input data cube of dimensions [batch_size x max_observations x n_input]
            y           Target data Tensor of dimensions [batch_size x max_observations]
            seq_lenghts Number of observations for each batch if observation < max_obs data is
                        padded with zeros [batch_size]
        input parameters:
            n_input     length of observed pixel values. [n_pixels * n_bands + n_time]
                n_pixels    number of observed pixels (default 3*3)
                n_bands     number of observed bands  (default 6)
                n_time      number of time parameters (default 1 e.g. day of year)
            n_classes   number of target classes
            batch_size  number of batches
            max_obs     maximum number of observations if seq_lengs < max_obs matrices will be padded
                        controls number of iterations in rnn layers (aka sequence length)
        network specific parameters
            n_layers    number of rnn layers (aka depth)
            learning_rate
            dropout_keep_prob
            logdir
    Marc.Russwurm@tum.de
    """

    def __init__(self, n_input=9 * 6 + 1, n_classes=20, batch_size=50, max_obs=26,
                 n_layers=2, dropout_keep_prob=.5, adam_lr=1e-3, adam_b1=0.9, adam_b2=0.999, adam_eps=1e-8,
                 fc_w_stddev=0.1, fc_b_offset=0.1, n_cell_per_input=1,rnn_cell_type="basiclstm", gpu=None):
        # save input arguments
        self.args = inspect.getargvalues(inspect.currentframe()).locals
        del self.args["self"]  # delete self argument

        self.n_classes = n_classes

        with tf.device(None):

            with tf.variable_scope('input'):
                # block of [batch_size x max_obs x n_input]
                self.X = tf.placeholder(tf.float32, [batch_size, max_obs, n_input], name="X")
                self.y_ = self.y = y_ = tf.placeholder(tf.float32, [batch_size, max_obs, n_classes], name="y")
                self.seq_lengths = seq_lengths = tf.placeholder(tf.int32, [batch_size], name="seq_lengths")
                #self.y = y = tf.reshape(self.y_, [-1, n_classes], name="y")

            with tf.name_scope('RNN'):
                self.n_rnn_cells = n_rnn_cells = n_cell_per_input * n_input

                if rnn_cell_type == "basiclstm":
                    cell = rnn_cell.BasicLSTMCell(n_rnn_cells)
                if rnn_cell_type == "lstm":
                    cell = rnn_cell.LSTMCell(n_rnn_cells)
                if rnn_cell_type == "lstm_peephole":
                    cell = rnn_cell.LSTMCell(n_rnn_cells, use_peepholes=True)
                elif rnn_cell_type == "gru":
                    cell = rnn_cell.BasicLSTMCell(n_rnn_cells)
                elif rnn_cell_type == "rnn":
                    cell = rnn_cell.BasicRNNCell(n_rnn_cells)

                # dropout Wrapper
                cell = tf.contrib.rnn.DropoutWrapper(cell=cell, output_keep_prob=dropout_keep_prob)
                self.cell = cell = rnn_cell.MultiRNNCell([cell] * n_layers)

                # tensor with class labels of dimension [batch_size x max_obs]
                # defined as Variable to carry values to next iteration (not trainable must be declared explicitly)
                self.state = state = cell.zero_state(batch_size, tf.float32)

                # rnn_outputs: block of [batch_size x max_obs x rnn_size]
                # data is padded with zeros after seq_length
                outputs, last_states = tf.nn.dynamic_rnn(cell, self.X, initial_state=state, sequence_length=seq_lengths,
                                                         time_major=False)

                self.outputs = outputs
                self.last_states = last_states

            with tf.name_scope('fc'):
                # reshape outputs to: block of [batch_size * max_obs x rnn_size]
                softmax_in = tf.reshape(outputs, [-1, n_rnn_cells])
                softmax_w = tf.Variable(tf.truncated_normal([n_rnn_cells, n_classes], stddev=fc_w_stddev), name="W_softmax")
                softmax_b = tf.Variable(tf.constant(fc_b_offset, shape=[n_classes]), name="b_softmax")

                softmax_out = tf.matmul(softmax_in, softmax_w) + softmax_b
                self.logits = logits = tf.reshape(softmax_out, [batch_size, -1, n_classes])

            with tf.name_scope('train'):
                # Define loss and optimizer

                # create mask for cross entropies incases where seq_lengths < max_max_obs
                # masking from http://stackoverflow.com/questions/34128104/tensorflow-creating-mask-of-varied-lengths

                with tf.name_scope('mask'):
                    lengths_transposed = tf.expand_dims(seq_lengths, 1)

                    range = tf.range(0, max_obs, 1)
                    range_row = tf.expand_dims(range, 0)

                    self.mask = mask = tf.less(range_row, lengths_transposed)

                self.cross_entropy_matrix = cross_entropy_matrix = tf.nn.softmax_cross_entropy_with_logits(logits=logits, labels=y_)
                self.masked_cross_entropy_matrix = masked_cross_entropy_matrix = tf.where(mask, cross_entropy_matrix,
                                                                                           tf.zeros(mask.get_shape()))
                self.cross_entropy_matrix = cross_entropy_matrix = tf.nn.softmax_cross_entropy_with_logits(logits=logits, labels=y_)

                # normalize with total number of observations
                self.cross_entropy = cross_entropy = tf.reduce_sum(cross_entropy_matrix) / tf.cast(
                    tf.reduce_sum(seq_lengths), tf.float32)
                tf.summary.scalar('cross_entropy', cross_entropy)
                # grad_train_op = tf.train.GradientDescentOptimizer(learning_rate=learning_rate).minimize(cross_entropy)
                self.train_op = tf.train.AdamOptimizer(learning_rate=adam_lr, beta1=adam_b1, beta2=adam_b2,
                                                       epsilon=adam_eps).minimize(cross_entropy)
                # tf.summary.scalar('learning_rate', learning_rate)

            with tf.name_scope('evaluation'):

                self.probabilities = probs = tf.nn.softmax(logits, name="full_probability_matrix")

                # Evaluate model
                predicted = tf.argmax(logits, 2)
                targets = tf.argmax(y_, 2)

                correct_pred = tf.equal(predicted, targets)
                masked_correct_pred = tf.logical_and(mask, correct_pred)
                self.accuracy_op = accuracy = tf.reduce_sum(tf.cast(masked_correct_pred, tf.float32)) / tf.cast(
                    tf.reduce_sum(seq_lengths), tf.float32)
                tf.summary.scalar('accuracy', accuracy)

                self.probs_list = probs_list = tf.reshape(probs, (-1, n_classes))
                predicted_list = tf.reshape(predicted, [-1])
                targets_list = tf.reshape(targets, [-1])

                mask_list = tf.reshape(mask, [-1])

                one_hot_targets = tf.one_hot(targets_list, n_classes)
                scores = tf.boolean_mask(probs_list, tf.cast(one_hot_targets, tf.bool))

                # mask of individual number of observations
                obs_list = tf.tile(tf.range(0, max_obs), [batch_size])
                obs_matrix = tf.matmul(tf.expand_dims(obs_list, 1), tf.ones([1, n_classes], dtype=tf.int32))



                probs_matrix_mask = probs_matrix_mask = tf.transpose(tf.reshape(tf.tile(mask_list, [n_classes]),[n_classes,-1]))

                self.scores = tf.boolean_mask(probs_list, probs_matrix_mask)
                self.targets = tf.boolean_mask(tf.reshape(y_, [-1,n_classes]), probs_matrix_mask)
                self.obs = tf.boolean_mask(obs_list, mask_list)

                # drop all values which are > seqlength
                #self.scores = tf.boolean_mask(scores, mask_list)
                #self.targets = tf.boolean_mask(targets_list, mask_list)
                #self.obs = tf.boolean_mask(obs_list, mask_list)

                self.confusion_matrix = confusion_matrix = tf.contrib.metrics.confusion_matrix(
                    tf.boolean_mask(targets_list, mask_list),
                    tf.boolean_mask(predicted_list, mask_list),
                    num_classes=n_classes)

                confusion_matrix = tf.cast(confusion_matrix, tf.uint8)
                confusion_matrix = tf.expand_dims(confusion_matrix, 2)
                confusion_matrix = tf.expand_dims(confusion_matrix, 0)
                tf.summary.image("confusion_matrix", confusion_matrix, max_outputs=3)

                logits_ = tf.cast(logits, tf.uint8)
                logits_ = tf.expand_dims(logits_, 3)
                tf.summary.image("logits", logits_, max_outputs=1)

                probs_ = tf.cast(probs*255, tf.uint8)
                probs_ = tf.expand_dims(probs_, 3)
                tf.summary.image("probabilities", probs_, max_outputs=1)

                targets_ = tf.cast(y_, tf.uint8)
                targets_ = tf.expand_dims(targets_, 3)
                tf.summary.image("targets", targets_, max_outputs=1)

                # tf.add_to_collection(tf.GraphKeys.SUMMARIES, cm_im_summary)

            # Merge all the summaries and write them out to /tmp/mnist_logs (by default)
            self.merge_summary_op = tf.summary.merge_all()
            self.init_op = tf.global_variables_initializer()

In [5]:
def unroll(x, y, seq_lengths):
    """
        Reshapes and masks input and output data from
        X(batchsize x n_max_obs x n_input) -> X_ (new_batchsize x n_input)
        y(batchsize x n_max_obs x n_classes) -> X_ (new_batchsize x n_classes)
        new_batch_size is variable representing batchsize * n_max_obs - invalid_observations
        with invalid observations being observations > seq_length -> means
        if at one point only 24 of maximum 26 images are available X is usually padded with zeros
        this masking removes the last two observations
    :return:
    """
    # create mask for valid times of acquisition
    batch_size, max_seqlengths, n_input = x.shape
    np.arange(0, max_seqlengths) * np.ones((batch_size, max_seqlengths))
    ones = np.ones([batch_size, max_seqlengths])
    mask = np.arange(0, max_seqlengths) * ones < (seq_lengths * ones.T).T

    new_x = x[mask]
    new_y = y[mask]

    return new_x, new_y

In [6]:
n_input = X.shape[2]
n_classes = y.shape[2]
batch_size = y.shape[0]
max_obs = X.shape[1]

seq_length = np.random.randint(22, 23, X.shape[0])
seq_length.fill(X.shape[1]) #all sequences = 23 observations

In [7]:
#1 LSTM layer with 50 cells (neurons) / neurons (n_input x n_cell_per_input)
model = Model(n_input=n_input, max_obs = max_obs, n_classes=n_classes, n_layers=1, batch_size=batch_size, adam_lr=1e-3, dropout_keep_prob=0.5, n_cell_per_input=10, rnn_cell_type = "lstm")

In [8]:
#run and write results

confusion_matrix = np.zeros((n_classes, n_classes), dtype=int)

summaryWriter = tf.summary.FileWriter("/home/acocac/post/tmp/log/real", graph=tf.get_default_graph())

with tf.Session() as sess:
    sess.run([model.init_op])

    feed = {model.X: X, model.y_: y, model.seq_lengths: seq_length}
    
    # training step
    for i in range(1, 30):
        train_op, cross_entropy, new_confusion_matrix = \
            sess.run([model.train_op,
                      model.cross_entropy,
                      model.confusion_matrix], feed_dict=feed)

        confusion_matrix += new_confusion_matrix
        print(cross_entropy)

        a,b = eval.calculate_accuracy_metrics(confusion_matrix)
        scores, targets = sess.run([model.scores, tf.reshape(model.targets, [-1])], feed_dict=feed)
        
        summary = sess.run(model.merge_summary_op, feed_dict=feed)
        summaryWriter.add_summary(summary, i)

    print("done")


0.695885
0.694988
0.693575
0.692781
0.691726
0.691208
0.689877
0.689185
0.688397
0.687387
0.686195
0.684994
0.684048
0.682736
0.681446
0.679826
0.678619
0.676765
0.674774
0.67235
0.670067
0.667228
0.663966
0.660451
0.655803
0.651314
0.645294
0.638299
0.629946
done

In [28]:
#computing performance metrix from confusion matrix
overall_accuracy, (c1, c2) = eval.calculate_accuracy_metrics(confusion_matrix)
print(overall_accuracy)


0.651592385625

Theano KERAS


In [9]:
#functions
def onehot(labels):
    Uniques,Index  = np.unique(labels,return_inverse=True)
    return np_utils.to_categorical(Index,len(Uniques))

def inverse_onehot(matrix):
    labels =[]
    for row in matrix:
        labels.append(np.argmax(row,axis=1))
    return labels

In [22]:
os.environ['THEANO_FLAGS']='mode=FAST_RUN,device=cpu'

In [24]:
import sklearn
from sklearn.model_selection import StratifiedKFold
from keras.models import Sequential
from keras.layers import Dense, Dropout
from keras.layers import LSTM
#from keras.optimizers import RMSprop
from keras.optimizers import Adam
from keras.utils import np_utils

In [19]:
LSTMRNN = {"name": 'Long Short Term Memory (LSTM)', "results" : [],"time": []}

In [10]:
#reshape label array
u = y.mean(axis=1,  keepdims=True) #mean probabilities 
a = inverse_onehot(u) #determine label
Y = np.array(a) 
Y = Y.mean(axis=1)

In [26]:
#fix random seed for reproducibility
seed = 7
np.random.seed(seed)

#define k-fold cross validation
folds = 2
kfold = StratifiedKFold(n_splits=folds, shuffle=True, random_state=seed) #ensure each fold is different and each sample is seen only once.

#LSTM and FFNN model settings
batch_size = batch_size 
n_epoch = 20 

#empty df
res = pd.DataFrame(columns=['acc', 'acc_std']) #to store results

for train, test in kfold.split(X, Y): 
    
    #input data LSTM
    train_data = X[train]
    test_data = X[test]

    model = Sequential()  
    model.add(LSTM(50, input_shape=(train_data.shape[1], train_data.shape[2])))
    model.add(Dropout(0.5))
    model.add(Dense(n_classes, activation='softmax'))
    
    model.compile(loss='categorical_crossentropy', optimizer=Adam(lr=1e-3, beta_1=0.9, beta_2=0.999, epsilon=1e-8), metrics=['acc'])

    model.fit(train_data, onehot(Y[train]), batch_size=batch_size, nb_epoch=n_epoch, verbose=0)

    scores = model.evaluate(test_data, onehot(Y[test]), verbose=0, batch_size=1000)
    print scores[0]
    
    LSTMRNN['results'].append(scores[1] * 100)
    res.loc[LSTMRNN['name'], 'acc'] = np.mean(LSTMRNN['results'])
    res.loc[LSTMRNN['name'], 'acc_std'] = np.std(LSTMRNN['results'])


0.161595404148
0.0789601057768

In [27]:
res


Out[27]:
acc acc_std
Long Short Term Memory (LSTM) 96.1636 1.15321