%% Convolution Neural Network Exercise

%  Instructions

%  ------------

%

%  This file contains code that helps you get started in building a single.

%  layer convolutional nerual network. In this exercise, you will only

%  need to modify cnnCost.m and cnnminFuncSGD.m. You will not need to

%  modify this file.

%%======================================================================

%% STEP : Initialize Parameters and Load Data

%  Here we initialize some parameters used for the exercise.

% Configuration

imageDim = ;

numClasses = ;  % Number of classes (MNIST images fall into  classes)

filterDim = ;    % Filter size for conv layer,*

numFilters = ;   % Number of filters for conv layer

poolDim = ;      % Pooling dimension, (should divide imageDim-filterDim+)

% Load MNIST Train

addpath ./common/;

images = loadMNISTImages('./common/train-images-idx3-ubyte');

images = reshape(images,imageDim,imageDim,[]);

labels = loadMNISTLabels('./common/train-labels-idx1-ubyte');

labels(labels==) = ; % Remap  to 

% Initialize Parameters,theta=(,),=**++**+

theta = cnnInitParams(imageDim,filterDim,numFilters,poolDim,numClasses);

%%======================================================================

%% STEP : Implement convNet Objective

%  Implement the function cnnCost.m.

%%======================================================================

%% STEP : Gradient Check

%  Use the file computeNumericalGradient.m to check the gradient

%  calculation for your cnnCost.m function.  You may need to add the

%  appropriate path or copy the file to this directory.

DEBUG=false;  % set this to true to check gradient

%DEBUG = true;

if DEBUG

    % To speed up gradient checking, we will use a reduced network and

    % a debugging data set

    db_numFilters = ;

    db_filterDim = ;

    db_poolDim = ;

    db_images = images(:,:,:);

    db_labels = labels(:);

    db_theta = cnnInitParams(imageDim,db_filterDim,db_numFilters,...

                db_poolDim,numClasses);

    [cost grad] = cnnCost(db_theta,db_images,db_labels,numClasses,...

                                db_filterDim,db_numFilters,db_poolDim);

    % Check gradients

    numGrad = computeNumericalGradient( @(x) cnnCost(x,db_images,...

                                db_labels,numClasses,db_filterDim,...

                                db_numFilters,db_poolDim), db_theta);

    % Use this to visually compare the gradients side by side

    disp([numGrad grad]);

    diff = norm(numGrad-grad)/norm(numGrad+grad);

    % Should be small. In our implementation, these values are usually

    % less than 1e-.

    disp(diff); 

    assert(diff < 1e-,...

        'Difference too large. Check your gradient computation again');

end;

%%======================================================================

%% STEP : Learn Parameters

%  Implement minFuncSGD.m, then train the model.

% 因为是采用的mini-batch梯度下降法，所以总共对样本的循环次数次数比标准梯度下降法要少

% 很多，因为每次循环中权值已经迭代多次了

options.epochs = ;

options.minibatch = ;

options.alpha = 1e-;

options.momentum = .;

opttheta = minFuncSGD(@(x,y,z) cnnCost(x,y,z,numClasses,filterDim,...

                      numFilters,poolDim),theta,images,labels,options);

save('theta.mat','opttheta');             

%%======================================================================

%% STEP : Test

%  Test the performance of the trained model using the MNIST test set. Your

%  accuracy should be above % after  epochs of training

testImages = loadMNISTImages('./common/t10k-images-idx3-ubyte');

testImages = reshape(testImages,imageDim,imageDim,[]);

testLabels = loadMNISTLabels('./common/t10k-labels-idx1-ubyte');

testLabels(testLabels==) = ; % Remap  to 

[~,cost,preds]=cnnCost(opttheta,testImages,testLabels,numClasses,...

                filterDim,numFilters,poolDim,true);

acc = sum(preds==testLabels)/length(preds);

% Accuracy should be around 97.4% after  epochs

fprintf('Accuracy is %f\n',acc);

function convolvedFeatures = cnnConvolve(filterDim, numFilters, images, W, b)

%cnnConvolve Returns the convolution of the features given by W and b with

%the given images

%

% Parameters:

%  filterDim - filter (feature) dimension

%  numFilters - number of feature maps

%  images - large images to convolve with, matrix in the form

%           images(r, c, image number)

%  W, b - W, b for features from the sparse autoencoder,传进来的w已经是矩阵的形式

%         W is of shape (filterDim,filterDim,numFilters)

%         b is of shape (numFilters,)

%

% Returns:

%  convolvedFeatures - matrix of convolved features in the form

%                      convolvedFeatures(imageRow, imageCol, featureNum, imageNum)

numImages = size(images, );

imageDim = size(images, );

convDim = imageDim - filterDim + ;

convolvedFeatures = zeros(convDim, convDim, numFilters, numImages);

% Instructions:

%   Convolve every filter with every image here to produce the

%   (imageDim - filterDim + ) x (imageDim - filterDim + ) x numFeatures x numImages

%   matrix convolvedFeatures, such that

%   convolvedFeatures(imageRow, imageCol, featureNum, imageNum) is the

%   value of the convolved featureNum feature for the imageNum image over

%   the region (imageRow, imageCol) to (imageRow + filterDim - , imageCol + filterDim - )

%

% Expected running times:

%   Convolving with  images should take less than  seconds

%   Convolving with  images should take around  minutes

%   (So to save time when testing, you should convolve with less images, as

%   described earlier)

for imageNum = :numImages

  for filterNum = :numFilters

    % convolution of image with feature matrix

    convolvedImage = zeros(convDim, convDim);

    % Obtain the feature (filterDim x filterDim) needed during the convolution

    %%% YOUR CODE HERE %%%

    filter = W(:,:,filterNum);

    bc = b(filterNum);

    % Flip the feature matrix because of the definition of convolution, as explained later

    filter = rot90(squeeze(filter),);

    % Obtain the image

    im = squeeze(images(:, :, imageNum));

    % Convolve "filter" with "im", adding the result to convolvedImage

    % be sure to do a 'valid' convolution

    %%% YOUR CODE HERE %%%

    convolvedImage = conv2(im, filter, 'valid');

    % Add the bias unit

    % Then, apply the sigmoid function to get the hidden activation

    %%% YOUR CODE HERE %%%

    convolvedImage = sigmoid(convolvedImage+bc);

    convolvedFeatures(:, :, filterNum, imageNum) = convolvedImage;

  end

end

end

function sigm = sigmoid(x)

    sigm = ./(+exp(-x));

end

function pooledFeatures = cnnPool(poolDim, convolvedFeatures)

%cnnPool Pools the given convolved features

%

% Parameters:

%  poolDim - dimension of pooling region

%  convolvedFeatures - convolved features to pool (as given by cnnConvolve)

%                      convolvedFeatures(imageRow, imageCol, featureNum, imageNum)

%

% Returns:

%  pooledFeatures - matrix of pooled features in the form

%                   pooledFeatures(poolRow, poolCol, featureNum, imageNum)

%     

numImages = size(convolvedFeatures, );

numFilters = size(convolvedFeatures, );

convolvedDim = size(convolvedFeatures, );

pooledFeatures = zeros(convolvedDim / poolDim, ...

        convolvedDim / poolDim, numFilters, numImages);

% Instructions:

%   Now pool the convolved features in regions of poolDim x poolDim,

%   to obtain the

%   (convolvedDim/poolDim) x (convolvedDim/poolDim) x numFeatures x numImages

%   matrix pooledFeatures, such that

%   pooledFeatures(poolRow, poolCol, featureNum, imageNum) is the

%   value of the featureNum feature for the imageNum image pooled over the

%   corresponding (poolRow, poolCol) pooling region.

%

%   Use mean pooling here.

%%% YOUR CODE HERE %%%

    % convolvedFeatures(imageRow, imageCol, featureNum, imageNum)

    % pooledFeatures(poolRow, poolCol, featureNum, imageNum)

    for numImage = :numImages

        for numFeature = :numFilters

            for poolRow = :convolvedDim / poolDim

                offsetRow = +(poolRow-)*poolDim;

                for poolCol = :convolvedDim / poolDim

                    offsetCol = +(poolCol-)*poolDim;

                    patch = convolvedFeatures(offsetRow:offsetRow+poolDim-, ...

                        offsetCol:offsetCol+poolDim-,numFeature,numImage); %取出一个patch

                    pooledFeatures(poolRow,poolCol,numFeature,numImage) = mean(patch(:));

                end

            end

        end

    end

end

function [cost, grad, preds] = cnnCost(theta,images,labels,numClasses,...

                                filterDim,numFilters,poolDim,pred)

% Calcualte cost and gradient for a single layer convolutional

% neural network followed by a softmax layer with cross entropy

% objective.

%

% Parameters:

%  theta      -  unrolled parameter vector

%  images     -  stores images in imageDim x imageDim x numImges

%                array

%  numClasses -  number of classes to predict

%  filterDim  -  dimension of convolutional filter

%  numFilters -  number of convolutional filters

%  poolDim    -  dimension of pooling area

%  pred       -  boolean only forward propagate and return

%                predictions

%

%

% Returns:

%  cost       -  cross entropy cost

%  grad       -  gradient with respect to theta (if pred==False)

%  preds      -  list of predictions for each example (if pred==True)

if ~exist('pred','var')

    pred = false;

end;

imageDim = size(images,); % height/width of image

numImages = size(images,); % number of images

lambda = 3e-; % weight decay parameter     

%% Reshape parameters and setup gradient matrices

% Wc is filterDim x filterDim x numFilters parameter matrix

% bc is the corresponding bias

% Wd is numClasses x hiddenSize parameter matrix where hiddenSize

% is the number of output units from the convolutional layer

% bd is corresponding bias

[Wc, Wd, bc, bd] = cnnParamsToStack(theta,imageDim,filterDim,numFilters,...

                        poolDim,numClasses); %the theta vector cosists wc,wd,bc,bd in order

% Same sizes as Wc,Wd,bc,bd. Used to hold gradient w.r.t above params.

Wc_grad = zeros(size(Wc));

Wd_grad = zeros(size(Wd));

bc_grad = zeros(size(bc));

bd_grad = zeros(size(bd));

%%======================================================================

%% STEP 1a: Forward Propagation

%  In this step you will forward propagate the input through the

%  convolutional and subsampling (mean pooling) layers.  You will then use

%  the responses from the convolution and pooling layer as the input to a

%  standard softmax layer.

%% Convolutional Layer

%  For each image and each filter, convolve the image with the filter, add

%  the bias and apply the sigmoid nonlinearity.  Then subsample the

%  convolved activations with mean pooling.  Store the results of the

%  convolution in activations and the results of the pooling in

%  activationsPooled.  You will need to save the convolved activations for

%  backpropagation.

convDim = imageDim-filterDim+; % dimension of convolved output

outputDim = (convDim)/poolDim; % dimension of subsampled output

% convDim x convDim x numFilters x numImages tensor for storing activations

activations = zeros(convDim,convDim,numFilters,numImages);

% outputDim x outputDim x numFilters x numImages tensor for storing

% subsampled activations

activationsPooled = zeros(outputDim,outputDim,numFilters,numImages);

%%% YOUR CODE HERE %%%

convolvedFeatures = cnnConvolve(filterDim, numFilters, images, Wc, bc); %前向传播,已经经过了激发函数

activationsPooled = cnnPool(poolDim, convolvedFeatures);

% Reshape activations into -d matrix, hiddenSize x numImages,

% for Softmax layer

activationsPooled = reshape(activationsPooled,[],numImages);

%% Softmax Layer

%  Forward propagate the pooled activations calculated above into a

%  standard softmax layer. For your convenience we have reshaped

%  activationPooled into a hiddenSize x numImages matrix.  Store the

%  results in probs.

% numClasses x numImages for storing probability that each image belongs to

% each class.

probs = zeros(numClasses,numImages);

%%% YOUR CODE HERE %%%

%Wd=(numClasses,hiddenSize),probs的每一列代表一个输出

M = Wd*activationsPooled+repmat(bd,[,numImages]);

M = bsxfun(@minus,M,max(M,[],));

M = exp(M);

probs = bsxfun(@rdivide, M, sum(M)); %why rdivide?

%%======================================================================

%% STEP 1b: Calculate Cost

%  In this step you will use the labels given as input and the probs

%  calculate above to evaluate the cross entropy objective.  Store your

%  results in cost.

cost = ; % save objective into cost

%%% YOUR CODE HERE %%%

% cost = -/numImages*labels(:)'*log(probs(:));

% 首先需要把labels弄成one-hot编码

groundTruth = full(sparse(labels, :numImages, ));

cost = -./numImages*groundTruth(:)'*log(probs(:))+(lambda/2.)*(sum(Wd(:).^2)+sum(Wc(:).^2)); %加入一个惩罚项

% cost = -./numImages*groundTruth(:)'*log(probs(:));

% Makes predictions given probs and returns without backproagating errors.

if pred

    [~,preds] = max(probs,[],);

    preds = preds';

    grad = ;

    return;

end;

%% 将c步和d步合成在一起了

%======================================================================

% STEP 1c: Backpropagation

%  Backpropagate errors through the softmax and convolutional/subsampling

%  layers.  Store the errors for the next step to calculate the gradient.

%  Backpropagating the error w.r.t the softmax layer is as usual.  To

%  backpropagate through the pooling layer, you will need to upsample the

%  error with respect to the pooling layer for each filter and each image.

%  Use the kron function and a matrix of ones to do this upsampling

%  quickly.

% STEP 1d: Gradient Calculation

%  After backpropagating the errors above, we can use them to calculate the

%  gradient with respect to all the parameters.  The gradient w.r.t the

%  softmax layer is calculated as usual.  To calculate the gradient w.r.t.

%  a filter in the convolutional layer, convolve the backpropagated error

%  for that filter with each image and aggregate over images.

%%% YOUR CODE HERE %%%

%%% YOUR CODE HERE %%%

% 网络结构: images--> convolvedFeatures--> activationsPooled--> probs

% Wd = (numClasses,hiddenSize)

% bd = (hiddenSize,)

% Wc = (filterDim,filterDim,numFilters)

% bc = (numFilters,)

% activationsPooled = zeros(outputDim,outputDim,numFilters,numImages);

% convolvedFeatures = (convDim,convDim,numFilters,numImages)

% images(imageDim,imageDim,numImges)

delta_d = -(groundTruth-probs); % softmax layer's preactivation,每一个样本都对应有自己每层的误差敏感性。

Wd_grad = (./numImages)*delta_d*activationsPooled'+lambda*Wd;

bd_grad = (./numImages)*sum(delta_d,); %注意这里是要求和

delta_s = Wd'*delta_d; %the pooling/sample layer's preactivation

delta_s = reshape(delta_s,outputDim,outputDim,numFilters,numImages);

%进行unsampling操作，由于kron函数只能对二维矩阵操作，所以还得分开弄

%delta_c = convolvedFeatures.*(-convolvedFeatures).*(./poolDim^)*kron(delta_s, ones(poolDim));

delta_c = zeros(convDim,convDim,numFilters,numImages);

for i=:numImages

    for j=:numFilters

        delta_c(:,:,j,i) = (./poolDim^)*kron(squeeze(delta_s(:,:,j,i)), ones(poolDim));

    end

end

delta_c = convolvedFeatures.*(-convolvedFeatures).*delta_c;

% Wc_grad = convn(images,rot90(delta_c,,'valid'))+ lambda*Wc;

for i=:numFilters

    Wc_i = zeros(filterDim,filterDim);

    for j=:numImages

        Wc_i = Wc_i+conv2(squeeze(images(:,:,j)),rot90(squeeze(delta_c(:,:,i,j)),),'valid');

    end

   % Wc_i = convn(images,rot180(squeeze(delta_c(:,:,i,:))),'valid');

    % add penalize

    Wc_grad(:,:,i) = (./numImages)*Wc_i+lambda*Wc(:,:,i);

    bc_i = delta_c(:,:,i,:);

    bc_i = bc_i(:);

    bc_grad(i) = sum(bc_i)/numImages;

end

%% Unroll gradient into grad vector for minFunc

grad = [Wc_grad(:) ; Wd_grad(:) ; bc_grad(:) ; bd_grad(:)];

end

function X = rot180(X)

    X = flipdim(flipdim(X, ), );

end

function [opttheta] = minFuncSGD(funObj,theta,data,labels,...

                        options)

% Runs stochastic gradient descent with momentum to optimize the

% parameters for the given objective.

%

% Parameters:

%  funObj     -  function handle which accepts as input theta,

%                data, labels and returns cost and gradient w.r.t

%                to theta.

%  theta      -  unrolled parameter vector

%  data       -  stores data in m x n x numExamples tensor

%  labels     -  corresponding labels in numExamples x  vector

%  options    -  struct to store specific options for optimization

%

% Returns:

%  opttheta   -  optimized parameter vector

%

% Options (* required)

%  epochs*     - number of epochs through data

%  alpha*      - initial learning rate

%  minibatch*  - size of minibatch

%  momentum    - momentum constant, defualts to 0.9

%%======================================================================

%% Setup

assert(all(isfield(options,{'epochs','alpha','minibatch'})),...

        'Some options not defined');

if ~isfield(options,'momentum')

    options.momentum = 0.9;

end;

epochs = options.epochs;

alpha = options.alpha;

minibatch = options.minibatch;

m = length(labels); % training set size

% Setup for momentum

mom = 0.5;

momIncrease = ;

velocity = zeros(size(theta));

%%======================================================================

%% SGD loop

it = ;

for e = :epochs

    % randomly permute indices of data for quick minibatch sampling

    rp = randperm(m);

    for s=:minibatch:(m-minibatch+)

        it = it + ;

        % increase momentum after momIncrease iterations

        if it == momIncrease

            mom = options.momentum;

        end;

        % get next randomly selected minibatch

        mb_data = data(:,:,rp(s:s+minibatch-)); % 取出当前的mini-batch的训练样本

        mb_labels = labels(rp(s:s+minibatch-));

        % evaluate the objective function on the next minibatch

        [cost grad] = funObj(theta,mb_data,mb_labels);

        % Instructions: Add in the weighted velocity vector to the

        % gradient evaluated above scaled by the learning rate.

        % Then update the current weights theta according to the

        % sgd update rule

        %%% YOUR CODE HERE %%%

        velocity = mom*velocity+alpha*grad; % 见ufldl教程Optimization: Stochastic Gradient Descent

        theta = theta-velocity;

        fprintf('Epoch %d: Cost on iteration %d is %f\n',e,it,cost);

    end;

    % aneal learning rate by factor of two after each epoch

    alpha = alpha/2.0;

end;

opttheta = theta;

end

function numgrad = computeNumericalGradient(J, theta)

% numgrad = computeNumericalGradient(J, theta)

% theta: a vector of parameters

% J: a function that outputs a real-number. Calling y = J(theta) will return the

% function value at theta. 

% Initialize numgrad with zeros

numgrad = zeros(size(theta));

%% ---------- YOUR CODE HERE --------------------------------------

% Instructions:

% Implement numerical gradient checking, and return the result in numgrad.

% (See Section 2.3 of the lecture notes.)

% You should write code so that numgrad(i) is (the numerical approximation to) the

% partial derivative of J with respect to the i-th input argument, evaluated at theta.

% I.e., numgrad(i) should be the (approximately) the partial derivative of J with

% respect to theta(i).

%

% Hint: You will probably want to compute the elements of numgrad one at a time. 

epsilon = 1e-;

for i =:length(numgrad)

    oldT = theta(i);

    theta(i)=oldT+epsilon;

    pos = J(theta);

    theta(i)=oldT-epsilon;

    neg = J(theta);

    numgrad(i) = (pos-neg)/(*epsilon);

    theta(i)=oldT;

    if mod(i,)==

       fprintf('Done with %d\n',i);

    end;

end;

%% ---------------------------------------------------------------

end

function theta = cnnInitParams(imageDim,filterDim,numFilters,...

                                poolDim,numClasses)

% Initialize parameters for a single layer convolutional neural

% network followed by a softmax layer.

%

% Parameters:

%  imageDim   -  height/width of image

%  filterDim  -  dimension of convolutional filter

%  numFilters -  number of convolutional filters

%  poolDim    -  dimension of pooling area

%  numClasses -  number of classes to predict

%

%

% Returns:

%  theta      -  unrolled parameter vector with initialized weights

%% Initialize parameters randomly based on layer sizes.

assert(filterDim < imageDim,'filterDim must be less that imageDim');

Wc = 1e-*randn(filterDim,filterDim,numFilters);

outDim = imageDim - filterDim + ; % dimension of convolved image

% assume outDim is multiple of poolDim

assert(mod(outDim,poolDim)==,...

       'poolDim must divide imageDim - filterDim + 1');

outDim = outDim/poolDim;

hiddenSize = outDim^*numFilters;

% we'll choose weights uniformly from the interval [-r, r]

r  = sqrt() / sqrt(numClasses+hiddenSize+);

Wd = rand(numClasses, hiddenSize) *  * r - r;

bc = zeros(numFilters, ); %初始化为0

bd = zeros(numClasses, );

% Convert weights and bias gradients to the vector form.

% This step will "unroll" (flatten and concatenate together) all

% your parameters into a vector, which can then be used with minFunc.

theta = [Wc(:) ; Wd(:) ; bc(:) ; bd(:)];

end

function [Wc, Wd, bc, bd] = cnnParamsToStack(theta,imageDim,filterDim,...

                                 numFilters,poolDim,numClasses)

% Converts unrolled parameters for a single layer convolutional neural

% network followed by a softmax layer into structured weight

% tensors/matrices and corresponding biases

%

% Parameters:

%  theta      -  unrolled parameter vectore

%  imageDim   -  height/width of image

%  filterDim  -  dimension of convolutional filter

%  numFilters -  number of convolutional filters

%  poolDim    -  dimension of pooling area

%  numClasses -  number of classes to predict

%

%

% Returns:

%  Wc      -  filterDim x filterDim x numFilters parameter matrix

%  Wd      -  numClasses x hiddenSize parameter matrix, hiddenSize is

%             calculated as numFilters*((imageDim-filterDim+)/poolDim)^

%  bc      -  bias for convolution layer of size numFilters x

%  bd      -  bias for dense layer of size hiddenSize x 

outDim = (imageDim - filterDim + )/poolDim;

hiddenSize = outDim^*numFilters;

%% Reshape theta

indS = ;

indE = filterDim^*numFilters;

Wc = reshape(theta(indS:indE),filterDim,filterDim,numFilters);

indS = indE+;

indE = indE+hiddenSize*numClasses;

Wd = reshape(theta(indS:indE),numClasses,hiddenSize);

indS = indE+;

indE = indE+numFilters;

bc = theta(indS:indE);

bd = theta(indE+:end);

end

tmp = conv2(convolvedFeatures(:,:,numFeature,numImage), ones(poolDim),'valid'); 
pooledFeatures(:,:,numFeature,numImage) =./(poolDim^)*tmp(:poolDim:end,:poolDim:end);