'''Train a Siamese MLP on pairs of digits from the MNIST dataset.
It follows Hadsell-et-al.'06 [1] by computing the Euclidean distance on the
output of the shared network and by optimizing the contrastive loss (see paper
for mode details).
[1] "Dimensionality Reduction by Learning an Invariant Mapping"
    http://yann.lecun.com/exdb/publis/pdf/hadsell-chopra-lecun-06.pdf
Gets to 99.5% test accuracy after 20 epochs.
3 seconds per epoch on a Titan X GPU
'''

import numpy as np
import random

# Keras imports
from keras.datasets import mnist
from keras.models import Sequential, Model
from keras.layers import Dense, Dropout, Input, Lambda
from keras.optimizers import RMSprop
from keras import backend as K

def euclidean_distance(vects):
    x, y = vects
    return K.sqrt(K.maximum(K.sum(K.square(x - y), axis=1, keepdims=True), K.epsilon()))

def eucl_dist_output_shape(shapes):
    shape1, shape2 = shapes
    return (shape1[0], 1)

def contrastive_loss(y_true, y_pred):
    '''Contrastive loss from Hadsell-et-al.'06
    http://yann.lecun.com/exdb/publis/pdf/hadsell-chopra-lecun-06.pdf
    '''
    margin = 1
    return K.mean(y_true * K.square(y_pred) +
                  (1 - y_true) * K.square(K.maximum(margin - y_pred, 0)))

def create_pairs(x, digit_indices):
    '''Positive and negative pair creation.
    Alternates between positive and negative pairs.
    '''
    pairs = []
    labels = []
    n = min([len(digit_indices[d]) for d in range(10)]) - 1
    for d in range(10):
        for i in range(n):
            z1, z2 = digit_indices[d][i], digit_indices[d][i + 1]
            pairs += [[x[z1], x[z2]]]
            inc = random.randrange(1, 10)
            dn = (d + inc) % 10
            z1, z2 = digit_indices[d][i], digit_indices[dn][i]
            pairs += [[x[z1], x[z2]]]
            labels += [1, 0]
    return np.array(pairs), np.array(labels)

def create_base_network(input_dim):
    '''Base network to be shared (eq. to feature extraction).
    '''
    seq = Sequential()
    seq.add(Dense(128, input_shape=(input_dim,), activation='relu'))
    seq.add(Dropout(0.1))
    seq.add(Dense(128, activation='relu'))
    seq.add(Dropout(0.1))
    seq.add(Dense(128, activation='relu'))
    return seq

def compute_accuracy(predictions, labels):
    '''Compute classification accuracy with a fixed threshold on distances.
    '''
    return labels[predictions.ravel() < 0.5].mean()

# the data, shuffled and split between train and test sets
(x_train, y_train), (x_test, y_test) = mnist.load_data()

x_train.shape, y_train.shape

x_test.shape, y_test.shape

x_train = x_train.reshape(60000, 784)
x_test = x_test.reshape(10000, 784)
x_train = x_train.astype('float32')
x_test = x_test.astype('float32')

"""
Explanation: Breakdown of network to make sense of components.
+ It was interesting to see the manually constructed loss function
+ multi input network
+ shared weights
Need to clean up the notebook a bit
End of explanation
"""


# normalize
x_train /= 255
x_test /= 255
# flatten images
input_dim = 784
epochs = 20

# create training+test positive and negative pairs
digit_indices = [np.where(y_train == i)[0] for i in range(10)]
tr_pairs, tr_y = create_pairs(x_train, digit_indices)

digit_indices = [np.where(y_test == i)[0] for i in range(10)]
te_pairs, te_y = create_pairs(x_test, digit_indices)

tr_pairs.shape

"""
Explanation: So basically here we are just flattening the images dimensions
End of explanation
"""


tr_pairs[:,0].shape

tr_pairs[:,1].shape

tr_y.shape

len(digit_indices)

np.where(y_test == 2)

"""
Explanation: so this dimensions are basically:
+ 108400  is the number of items in the list
+ 2 each entry in the list has 2 parts (which are the pairs)
+ 784 is the dimension of the flatten image
scratch that. There is no list anymore since the function calls for np.array() before returning
End of explanation
"""


np.where(y_test == 2)[0]

"""
Explanation: np.where returns a list, so you need the [0] to get a array only
End of explanation
"""


digit_indices[0][2], digit_indices[0][2+1]

z1, z2 = digit_indices[0][2], digit_indices[0][2+1]

"""
Explanation: so the code above creates an array for each of the numbers, each array has the indices for each of the numbers

create_pairs function
python
def create_pairs(x, digit_indices):
    '''Positive and negative pair creation.
    Alternates between positive and negative pairs.
    '''
    pairs = []
    labels = []
    n = min([len(digit_indices[d]) for d in range(10)]) - 1
    for d in range(10):
        for i in range(n):
            z1, z2 = digit_indices[d][i], digit_indices[d][i + 1]
            pairs += [[x[z1], x[z2]]]
            inc = random.randrange(1, 10)
            dn = (d + inc) % 10
            z1, z2 = digit_indices[d][i], digit_indices[dn][i]
            pairs += [[x[z1], x[z2]]]
            labels += [1, 0]
    return np.array(pairs), np.array(labels)
End of explanation
"""


x_test[z1][:10]

x_test[z2][:10]

"""
Explanation: so we are getting the indices of two intances of the same number
End of explanation
"""


pairs += [[x_test[z1], x_test[z2]]]

len(pairs)

"""
Explanation: then we are concatenating them
End of explanation
"""


inc = random.randrange(1, 10)

inc

d = 4

dn = (d + inc) % 10

dn

"""
Explanation: so we keep building the pairs
End of explanation
"""


# this adds the negative pair
# z1, z2 = digit_indices[d][i], digit_indices[dn][i]
# pairs += [[x[z1], x[z2]]]

"""
Explanation: here we are building a pair with a true and false instance thats why the inc on dn
End of explanation
"""


# labels += [1, 0]

pairs = []
labels = []

n = min([len(digit_indices[d]) for d in range(10)]) - 1
# how many indices are the per number
# get the min of all the indices

n

 for d in range(10):
        print(len(digit_indices[d]))

# then it substract 1 - I assume its because the range is non-inclusive

n = 10
ran = 5

l = []

l += [1, 0]

l

for d in range(ran):
    for i in range(n):
        z1, z2 = digit_indices[d][i], digit_indices[d][i + 1]
        print(z1.shape, z2.shape)
        pairs += [[x_train[z1], x_train[z2]]]
        print(pairs)
        inc = random.randrange(1, 10)
        dn = (d + inc) % 10
        z1, z2 = digit_indices[d][i], digit_indices[dn][i]
        pairs += [[x_train[z1], x_train[z2]]]
        labels += [1, 0]
np.array(pairs), np.array(labels)

# network definition
base_network = create_base_network(input_dim)

input_a = Input(shape=(input_dim,))
input_b = Input(shape=(input_dim,))

base_network.summary()

# because we re-use the same instance `base_network`,
# the weights of the network will be shared across the two branches
processed_a = base_network(input_a)
processed_b = base_network(input_b)

distance = Lambda(euclidean_distance,
                  output_shape=eucl_dist_output_shape)([processed_a, processed_b])

model = Model([input_a, input_b], distance)

model.summary()

# loss is specified manually
# optimizer is RMSprop

tset = [tr_pairs[:, 0], tr_pairs[:, 1]]

len(tset)

len(tr_y)

# train
rms = RMSprop()
model.compile(loss=contrastive_loss, optimizer=rms)
model.fit([tr_pairs[:, 0], tr_pairs[:, 1]], tr_y,
          batch_size=128,
          epochs=epochs,
          validation_data=([te_pairs[:, 0], te_pairs[:, 1]], te_y))

# compute final accuracy on training and test sets
pred = model.predict([tr_pairs[:, 0], tr_pairs[:, 1]])
tr_acc = compute_accuracy(pred, tr_y)
pred = model.predict([te_pairs[:, 0], te_pairs[:, 1]])
te_acc = compute_accuracy(pred, te_y)

print('* Accuracy on training set: %0.2f%%' % (100 * tr_acc))
print('* Accuracy on test set: %0.2f%%' % (100 * te_acc))

"""
Explanation: the label is 1 for the first example and zero for the fake that we just made
End of explanation
"""


