random_state = np.random.RandomState( None )

split_ = cross_validation.train_test_split( mnist_data, mnist_labels, test_size = 0.75,
                                            random_state = random_state )
X_train, X_test, y_train, y_test = split_

print X_train.shape
print X_test.shape

"""
Explanation: <hr/>
End of explanation
"""


index_perm_ = np.random.permutation( X_train.shape[ 0 ] )
X = X_train[ index_perm_[ : 400 ] ]
fast_plot( X, u"200 random digits MNIST train dataset", n = 200 )

"""
Explanation: Draw and plot some random subset of the train dataset.
End of explanation
"""


train_digits_ = { label_ : X_train[ np.flatnonzero( y_train == label_ ) ]
                           for label_ in mnist_classes }

avg_patterns_ = { label_ : samples_.mean( axis = 0 )[ np.newaxis ]
                  for label_, samples_ in train_digits_.iteritems() }

## Construct an aligned pair of arrays
average_digits = np.concatenate( [ array_ for array_ in avg_patterns_.itervalues() ], axis = 0 )
average_labels = np.array( [ label_ for label_ in avg_patterns_.iterkeys() ], dtype = np.int )
del avg_patterns_

## Sort them
order_ = np.argsort( average_labels )
average_labels = average_labels[ order_ ]
average_digits = average_digits[ order_ ]

## Plot the patterns
fast_plot( average_digits, u"Averaged digits on train" )

"""
Explanation: Simple average pattern classifier
Generate averaged images
End of explanation
"""


## Compute the squared error ...
yp_ = - 2 * np.tensordot( X_test, average_digits, axes = [ -1, -1 ] ) \
      + average_digits.dot( average_digits.T )[:1]
## ... and based on it find the closest label.
pred_labels_ = average_labels[ yp_.argmin(axis = 1) ]

"""
Explanation: A simple minimal error calssifier based on mean patterns.
$\newcommand{\argmin}{\mathop{\mathtt{argmin}}} \newcommand{\Dcal}{\mathcal{D}}$
Basically this is a simple gaussian model:
$$ \log \mathcal{L} = - \frac{|\Dcal|}{2} \log 2\pi - \frac{|\Dcal|}{2} \log \sigma^2
                      - \frac{1}{2 \sigma^2} \sum_{x\in \Dcal} ( x - \mu_{k_x} )^2 \,. $$
Classification is done with this rule:
\begin{align}
    \hat{y}(x)
        & = \argmin_{k=1,\ldots,K} \, \sum_{n=1}^N \sum_{m=1}^M ( x_{nm} - \mu_{knm} )^2 \
        & = \argmin_{k=1,\ldots,K} \, \sum_{n=1}^N \sum_{m=1}^M - 2 x_{nm} \mu_{knm} + \mu_{knm}^2 \
        & = \argmin_{k=1,\ldots,K} \, \|\mu_k\|^2 - 2 \langle x, \mu_k \rangle \,,
\end{align}
where
$$ \langle f, g \rangle  = \sum_{n=1}^N \sum_{m=1}^M f_{nm} g_{nm}\,. $$
End of explanation
"""


confusion( y_test, pred_labels_ )

"""
Explanation: Construct the confusion matrix
End of explanation
"""


print "Achieved accuracy is %.3f%%" % ( np.mean( y_test == pred_labels_ ) * 100.0, )

"""
Explanation: Compute the accuracy
End of explanation
"""


rfc_ = ensemble.RandomForestClassifier( n_estimators = 50, n_jobs = -1,
                                        random_state = random_state ).fit( X_train, y_train )

"""
Explanation: Random Forest FTW!
Initalize and fit a random forest to the train dataset.
End of explanation
"""


pred_labels_ = rfc_.classes_[ rfc_.predict_proba( X_test ).argmax( axis = 1 ) ]

"""
Explanation: Predict on the test dataset.
End of explanation
"""


print "Achieved accuracy is %.3f%%" % ( np.mean( y_test == pred_labels_ ) * 100.0, )

"""
Explanation: Display the accuracy ...
End of explanation
"""


confusion( y_test, pred_labels_ )

"""
Explanation: ... and the confusion matrix.
End of explanation
"""


fast_plot( rfc_.feature_importances_[ np.newaxis ],
           u"Random Forest feature importances",
           n_row = 1, n_col = 1 )

"""
Explanation: Show the pixel-feature importances of the random forest.
End of explanation
"""


scl_ = preprocessing.StandardScaler( ).fit( X_train )
U, S, V = sp.linalg.svd( scl_.transform( X_train ), full_matrices = False )
order_ = S.argsort( )[::-1]
U, S, V = U[:,order_], S[ order_ ], V[order_]

"""
Explanation: PCA
Let's use PCA to learn a linear manifold from the data. Scipy's SVD returns $U$, $\Sigma$ and $V'$ (not $V$).
End of explanation
"""


var_ = S**2
leverage_ = np.cumsum( var_ ) / np.sum( var_ )
n_components = np.min( np.flatnonzero( leverage_ > 0.95 ) )

print """The least number of principal components required to""" \
    + """guarantee at least 95%% recostruction is %d.""" % ( n_components, ) 

plt.plot( leverage_ )
plt.axhline( y = 0.95, color = 'k', lw = 2 )

"""
Explanation: Compute variance within each principal direction.
End of explanation
"""


n_components = 256

"""
Explanation: Chose the number of components.
End of explanation
"""


pc_ = U[:,:n_components] * S[ :n_components ]

"""
Explanation: Extract the principal components
End of explanation
"""


X_ = scl_.inverse_transform( np.dot( pc_, V[:n_components] ) )

"""
Explanation: Embed the compnents in the original space.
End of explanation
"""


fast_plot( X_, u"Principal components of train dataset", n = 512 )

"""
Explanation: The eigenvectors are
End of explanation
"""


pc_ = np.dot( scl_.transform( X_test ), V[:n_components].T )
X_ = scl_.inverse_transform( np.dot( pc_, V[:n_components] ) )

"""
Explanation: See how well it fares on the test dataset: reconstruct ...
End of explanation
"""


fast_plot( np.abs(X_-X_test), u"Low rank reconstruction error of the test dataset", n = 32 )

"""
Explanation: ...  and display
End of explanation
"""


gnb_ = naive_bayes.GaussianNB( ).fit( X_train, y_train )

"""
Explanation: Naive Bayes
Let's use a Gaussian naive bayes classifier: $\newcommand\given{\:\vert\:}$
The Naive Bayes assumes that the features are conditionally independent, provided the class label $y$ is known.
In particular for any $i\in {1,\ldots,d}$ one has:
$$ p\bigl( x_1,\ldots,x_p \:\big\vert\: y, \Theta\bigr) = \prod_{i=1}^d p\bigl( x_i \:\big\vert\: y, \Theta\bigr) \,, $$
whence using the Bayes' rule :
\begin{align}
    p\bigl( y \:\big\vert\: x, \Theta\bigr)
            &= p\bigl( y \:\big\vert\: x, \Theta\bigr) \frac{p(y\:\big\vert\: \Theta)}{p(x\:\big\vert\: \Theta)} \
            &\propto p\bigl(y\:\big\vert\: \Theta\bigr)
                    \prod_{i=1}^d p\bigl( x_i \:\big\vert\: y, \Theta\bigr) \,,
\end{align}
where $\Theta$ is the parameter vector of the model.
Prediction is done using the maximum posterior porbability principle: $\newcommand{\argmax}{\mathop{\mathtt{argmax}}}$
$$ \hat{y}(x) = \argmax_{y} \log p\bigl(y\:\big\vert\: \Theta\bigr)
                            + \sum_{i=1}^d \log p\bigl( x_i \:\big\vert\: y, \Theta\bigr) \,. $$
Fitting the Naive Bayes classifier is done unsing the ML approach:
$$  \hat{\Theta} = \argmax_\Theta \sum_{j=1}^J \log p\bigl(y_j\:\big\vert\: \Theta\bigr)
                            + \sum_{j=1}^J \sum_{i=1}^d \log p\bigl( x_{ji} \:\big\vert\: y_j, \Theta\bigr) \,. $$
Gaussian NB classifier assumes that 
$$ p\bigl( x_i \:\big\vert\: y, \Theta\bigr) = \bigl(2\pi \sigma^2_{yi}\bigr)^{-\frac{1}{2}} \mathtt{exp}\biggl( -\frac{1}{2\sigma^2_{yi}} (x_i - \mu_{yi} ) \biggr) $$
End of explanation
"""


pred_labels_ = gnb_.classes_[ gnb_.predict_proba( X_test ).argmax( axis = 1 ) ]

"""
Explanation: Predict the labels ...
End of explanation
"""


print np.mean( pred_labels_ == y_test )

"""
Explanation: ... and report the accuracy.
End of explanation
"""


confusion( y_test, pred_labels_ )

"""
Explanation: Display the confusion matrix
End of explanation
"""


fast_plot( gnb_.theta_, u"Estimated mean values of each pixel feature given a class", n = 10 )

"""
Explanation: The estimated mean of each pixel form the following patterns:
End of explanation
"""


np.allclose( gnb_.theta_, average_digits )

"""
Explanation: They seem identical to the average patterns and indeed they are:
End of explanation
"""


fast_plot( gnb_.sigma_, u"GNB Estimated standard deviation of individual pixels", n = 10 )

"""
Explanation: Let's see the estimated variablility of each pixel.
End of explanation
"""


X_train_ = scl_.transform( X_train )
lorl2_ = linear_model.LogisticRegression( C = 1.0, penalty = "l2" ).fit( X_train_, y_train )

"""
Explanation: Logistic Regression : $l_2$
Let's employ logisitc regression to classify the digits.
End of explanation
"""


pred_labels_ = lorl2_.predict( scl_.transform( X_test ) )
print "Accuracy score is %.3f%%" % ( 100*np.mean( pred_labels_ == y_test ), )

"""
Explanation: Predict on the test.
End of explanation
"""


confusion( y_test, pred_labels_ )

"""
Explanation: Display a confusion matrix
End of explanation
"""


fast_plot( lorl2_.coef_, u"Weights of the $l_2$ logistic regression", n = 10 )

"""
Explanation: Let's visualize the weights of the estimated logistic regression.
End of explanation
"""


X_train_ = scl_.transform( X_train )
lorl1_ = linear_model.LogisticRegression( C = 1.0, penalty = "l1" ).fit( X_train_, y_train )

"""
Explanation: Test an $l_1$ regularization.
End of explanation
"""


pred_labels_ = lorl1_.predict( scl_.transform( X_test ) )
print "Accuracy score is %.3f%%" % ( 100*np.mean( pred_labels_ == y_test ), )

"""
Explanation: Report accuracy
End of explanation
"""


confusion( y_test, pred_labels_ )

"""
Explanation: confusion matrix
End of explanation
"""


fast_plot( lorl1_.coef_, u"Weights of the $l_1$ logistic regression", n = 10 )

"""
Explanation: and display the weights
End of explanation
"""


X_train_2d_ = np.reshape( X_train_, ( -1, 28, 28 ) )
fft2_2d_ = np.fft.fft2( X_train_2d_ )

fft2_2d_.shape

fast_plot( fft2_2d_.real, u"2D FFT of some MNIST digits", n = 400 )

fast_plot( fft2_2d_.imag, u"2D FFT of some MNIST digits", n = 400 )

"""
Explanation: FFT
Let's attempt convolutions.
End of explanation
"""


X_train_ = scl_.transform( X_train )
lorl1_3_ = linear_model.LogisticRegression( C = 2.0, penalty = "l1" ).fit( X_train_, y_train )


"""
Explanation: <hr/>

Try a higher coefficient for regularization.
End of explanation
"""


pred_labels_ = lorl1_3_.predict( scl_.transform( X_test ) )
print "Accuracy score is %.3f%%" % ( 100*np.mean( pred_labels_ == y_test ), )
print confusion( y_test, pred_labels_ )
fast_plot( lorl1_3_.coef_, u"Weights of the $l_1$ logistic regression", n = 10 )

"""
Explanation: Report
End of explanation
"""


g_lor_ = grid_search.GridSearchCV( linear_model.LogisticRegression( ),
                                   param_grid = {
                                        "C" : np.logspace( -2, 2, num = 5 ),
                                        "penalty" : [ "l2", "l1" ],
                                  }, cv = 7, n_jobs = -1, verbose = 50 ).fit( X_train_, y_train )

"""
Explanation: Run a Grid seach over the penalty and the degree of regularization.
End of explanation
"""


g_lor_.scores_

"""
Explanation: Let's see what the cross-validation consider the best model.
End of explanation
"""


scl_ = preprocessing.StandardScaler( ).fit( X_train )
X_train_ = scl_.transform( X_train )
pca_ = decomposition.PCA( n_components = 15 ).fit( X_train_ )

"""
Explanation: PCA + $k$-nearest neighbour classifier
Reduce dimensionality using PCA to 64.
End of explanation
"""


## Don't set the number of neighbours just yet.
PC_train = pca_.transform( X_train_ )
knn_ = neighbors.KNeighborsClassifier( n_neighbors = None, n_jobs = -1 ).fit( PC_train, y_train )

def knn_call( knn, method, k, *args, **kwargs ) :
    knn.n_neighbors = k
    return getattr( knn, method )( *args, **kwargs )

"""
Explanation: The $k$-nn implementation in Scikit-Learn, appears to prohibit chaning the number of neighbors on-the-fly.
End of explanation
"""


n_neighbors = np.asarray( np.linspace( 1, 13, num = 5 ), np.int )

PC_test = pca_.transform( scl_.transform( X_test ) )
pred_labels_ = { k : knn_call( knn_, "predict", k = k, X = PC_test )
                   for k in n_neighbors }

"""
Explanation: Predict the labels using the majority voting procedure.
End of explanation
"""


import scipy as sp
for k in n_neighbors :
    tbl = confusion( y_test, pred_labels_[ k ] )
    print "%d-nn : accuracy %0.3f%%\n" % ( k, 100*np.mean( pred_labels_[ k ] == y_test ), ), tbl, "\n"

"""
Explanation: Display accuracy of the $k$-nn classifiers and the confusion matrix.
End of explanation
"""


pipe_ = pipeline.Pipeline( [
            ( "scaler1", preprocessing.StandardScaler( ) ),
            ( "pca", decomposition.PCA( n_components = 100 ) ),
            ( "scaler2", preprocessing.StandardScaler( ) ),
            ( "lore", linear_model.LogisticRegression( C = 1.0, penalty = "l2" ) ),
        ] ).fit( X_train, y_train )

"""
Explanation: <hr/>
End of explanation
"""


pred_labels_ = pipe_.predict( X_test )

"""
Explanation: Now test the pipeline.
End of explanation
"""


print "Accuracy score is %.3f%%" % ( 100*np.mean( pred_labels_ == y_test ), )
confusion( y_test, pred_labels_ )

"""
Explanation: display the summary scores.
End of explanation
"""


import lasagne, theano.tensor as T

"""
Explanation: <hr/>
End of explanation
"""


input_var = T.tensor4( 'inputs' )
target_var = T.ivector( 'targets' ) 

## The input layer
network = lasagne.layers.InputLayer( shape = ( None, 1, 28, 28 ),
                                     input_var = input_var )

## The 2D 5x5 conv layer with ReLU and 2x2 max-pooling
## (28-7+1) // 1
network = lasagne.layers.Conv2DLayer( network, num_filters = 32, filter_size = ( 5, 5 ),
                                      nonlinearity = lasagne.nonlinearities.rectify,
                                      border_mode = "valid", W = lasagne.init.GlorotUniform( ) )
network = lasagne.layers.MaxPool2DLayer( network, pool_size = ( 2, 2 ) )

## The 2D 3x3 conv layer with ReLU and 2x2 max-pooling
network = lasagne.layers.Conv2DLayer( lasagne.layers.dropout( network, p = .2 ),
                                      num_filters = 64, filter_size = ( 3, 3 ),
                                      nonlinearity = lasagne.nonlinearities.rectify,
                                      border_mode = "valid", W = lasagne.init.GlorotUniform( ) )
network = lasagne.layers.MaxPool2DLayer( network, pool_size = ( 2, 2 ) )

## FC layer with dropout
network = lasagne.layers.DenseLayer( lasagne.layers.dropout( network, p = .5 ),
                                     num_units = 256, nonlinearity = lasagne.nonlinearities.rectify )

network = lasagne.layers.DenseLayer( lasagne.layers.dropout( network, p = .5 ),
                                     num_units = 10, nonlinearity = lasagne.nonlinearities.softmax )


prediction = lasagne.layers.get_output( network )
loss = lasagne.objectives.categorical_crossentropy( prediction, target_var ).mean( )

params = lasagne.layers.get_all_params( network, trainable = True )
updates = lasagne.updates.nesterov_momentum( loss, params, learning_rate = 0.01, momentum = 0.9 )

test_prediction = lasagne.layers.get_output( network, deterministic = True )
test_loss = lasagne.objectives.categorical_crossentropy( test_prediction, target_var ).mean( )

test_acc = T.mean( T.eq( T.argmax( test_prediction, axis = 1 ), target_var ), dtype = theano.config.floatX )

train_fn = theano.function([input_var, target_var], loss, updates=updates)
val_fn = theano.function([input_var, target_var], [test_loss, test_acc])

"""
Explanation: Implement as simple neural network with lasagne
End of explanation
"""


