function x = normalize(x, mu, sigma)
    x=bsxfun(@minus,x,mu);
        x=bsxfun(@rdivide,x,sigma);
end

function nn=nnsetup(architecture)

    nn.size = architecture      // 即nn.size为[784 100 10]
    
    nn.n = numel(nn.size)       //即nn.n=nemel([784 100 10])，返回数组中元素个数，nn.n为3
    
    nn.activation_function = ‘tan_opt’ //activation函数使用tanh函数，也可使用'sigm'即sigmoid
    
    nn.learningRate = 2     // 学习率为2，使用sigmoid函数和非归一化输入时，值需要小些
    
    nn.momentum = 0.5
    
    nn.scaling_learningRate             //Scaling factor for the learning rate (each epoch)
    nn.weightPenaltyL2 = 0            //L2 regularization
    nn.nonSparsityPenalty = 0       //Non sparsity penalty
    nn.sparsityTarget = 0.05          //Sparsity target
    nn.inputZeroMaskedFraction = 0      //Used for Denoising AutoEncoders
    nn.dropoutFraction = 0            //Dropout level
    nn.testing = 0                        //Internal variable. nntest sets this to one.
    nn.output = 'sigm'                // output unit 'sigm' (=logistic), 'softmax' and 'linear'

    for i=2: nn.n
    
        nn.W{i-1}=(rand(nn.size(i), nn.size(i - 1)+1) - 0.5) * 2 * 4 * sqrt(6 / (nn.size(i) + nn.size(i - 1)));     //{}用于cell型数组的定义和引用，可以用nn.W{1}和nn.W{1}表示第一层连接参数和第二层连接参数；
        //在此W<1x2 cell>为cell型数组，其中只存储<100x785>与<10x101 double>这两个信息，即W{1}为<100x785 double>，即<100x(784+1)>，其中存储了参数具体值；W{2}为<10x101 double>，即<10x(101+1)>，其中存储了参数具体值。
       //（Theta矩阵即W矩阵，Theta{j} <S{j+1} x S{j}+1>，表示每行对应一个样本的输入的权值计算S{j}+1个，要输出下一层的S{j+1}个。因此W{1} 为<100 x 784+1>，W{2}为<10 x 100+1>）

        nn.vW{i - 1} = zeros(size(nn.W{i - 1}));        //vW的大小与W相同，只不过其值全为0
        //  W{1} 为<100 x 785>，W{2}为<10x101>
        
        nn.p{i} = zeros(1, nn.size(i));                 //average activations，用于sparsity

    end

function [nn, L] = nntrain(nn, train_x, train_y, opts, val_x, val_y)

// NNTRAIN训练神经网络，输入x，输出y，设置opts.numepochs epochs与minibatches of size opts.batchsize。返回updated activations, error, weights与biases的nn网络（结构体），以及各个traininig nimibatch的均方差之和，L（结构体）

    assert(isfloat(train_x), 'train_x must be a float');    //判断train_x是否为float型，如果不是，返回error 'train_x must be a float'     

    assert(nargin == 4 || nargin == 6,'number ofinput arguments must be 4 or 6')    //判断nntrain函数输入的变量，nargin为number of input arguments

    loss.train.e  = [];          
    loss.train.e_frac = [];
    loss.val.e = [];
    loss.val.e_frac = []; 

    opts.validation = 0;        //对验证集的判断
    
    if nargin == 6                  //验证集的使用，在有验证集的情况下，设置validation为1
       opts.validation = 1;  
    end

    fhandle = [];
    if isfield(opts,'plot') && opts.plot == 1 //是否需要画图表示
        fhandle = figure();    
    end

    m = size(train_x, 1);                       //m为train_x的行数，即为样本个数

    batchsize = opts.batchsize;             //每次计算的batchsize为100

    numepochs = opts.numepochs;         //epoch次数为1

    numbatches = m / batchsize;             //样本中batch的个数，为600个batch

    assert(rem(numbatches, 1) == 0, 'numbatches must be a integer');         //rem函数的求整除的函数，判断batch个数必须为整数

    L = zeros(numepochs*numbatches,1);           //每个batch计算存储一次平均loss值，即计算存储了600个batch的loss值，共计算numepoch次，即1次。L为<600x1 double>
    
    n=1;
    for i = 1 : numepochs           //计算迭代次数
        tic                             //用来保存当前时间
        
        kk = randperm(m);       //随机打乱一个数字序列，即kk为<1x60000 double>，其中顺序随机
        for l = 1 : numbatches          //计算次数为batch的个数
            batch_x = train_x( kk( ( l-1 ) * batchsize + 1 :  l * batchsize ), : );        
                // 当l=1时，batch_x = train_x( kk( 1 : 100), :) ，即kk(1: 100)为kk<1x60000>前100个随机数，train_x( , :)为自动填充，使结果满足train_x的列数，则batch_x=train_x(50, :)， batch_x=train_x(46, :) ..... batch_x=train_x(854, :)，总共有100此这样的batch_x赋值，则batch_x为<100x784 double>
                //另外，(l-1)* batch +1 : l* batchsize表示在第l个batch中（从0开始计算），从该btach起始处计算到结尾处，每次计算batchsize个样本，总共l个numbatch
                
        if(nn.inputZeroMaskedFraction ~= 0)
           batch_x = batch_x.*(rand(size(batch_x))>nn.inputZeroMaskedFraction);
        end
        //在auto-encoder中使用，此处无用
        
        batch_y = train_y (kk( (l - 1) * batchsize + 1 : l * batchsize), :);            //batch_y<100x10>
        
        nn = nnff( nn, batch_x, batch_y );
        // batch_x <100 x 784>，行向量为单样本特征； batch_y <100 x 10> ，行向量为单样本标签
        
        nn = nnbp( nn );
        
        nn = nnapplygrads( nn );
        
        L( n ) = nn.L;
        
        n = n+1;
        
    end
      
    t=toc               //使用toc记录程序完成时间

    if opt.validation == 1
        loss = nneval( nn, loss , train_x , train_y , val_x , val_y);
        str_perf = sprintf('; Full-batch train mse = %f, val mse = %f', loss.train.e(end), loss.val.e(end));
     //如果使用验证集，则使用上式计算验证集误差
     
    else
       loss = nneval( nn, loss , train_x , train_y);
       str_perf = sprintf('; Full-batch train err = %f', loss.train.e(end));
       //如果不存在验证集，使用上式计算误差
     end

        if ishandle( fhandle )
            nnupdatefigures( nn, fhandle , loss , opts , i);
        end
        //对句柄判断，更新图表

        disp(['epoch ' num2str(i) '/' num2str(opts.numepochs) '. Took ' num2str(t) ' seconds' '. Mini-batch mean squared error on training set is ' num2str(mean(L((n-numbatches):(n-1)))) str_perf]);
        //输出训练信息

        nn.learningRate = nn.learningRate * nn.scaling_learningRate;
        //更新learning rate

end
end

function nn = nnff ( nn, x, y)
//NNFF执行freedforward pass，nn=nnff(nn, x, y)返回更新了activation，error和loss的网络nn
//nnff()在nntrain()中调用，为nnff( nn, batch_x, batch_y)，则x为batch_x，y为batch_y

    n = nn.n              // 网络的层数，n=3
    m = size( x, 1)     //记录x的行数，即batch_x<100x784>的行数，表示一个batch的样本个数
    
    x = [ ones(m,1) x];     //初始batch_x<100x784>，在batch_x第一列处添加全为1的一列，<100 x 1+784>，其中第一列全为1，即样本x0=1
    nn.a{1} = x;            //将batch_x<100 x 785>存入nn.a{1}，即其存放样本值
    // （在《Machine Learning》中，nn.a为列向量，且一次只计算一次样本；在这里nn.a{1}为<100 x 785> 每一行向量代表一样本，总共batch_x中100个样本）
    
    for i = 2 : n-1                 //在这循环中，只计算i=2，即第二层参数
            swith nn.activation_function        //根据nn的activation function选择计算公式
                    case 'sigm'
                         nn.a{i} = sigm( nn.a{i - 1} * nn.W{i - 1}' );
                         //根据上一层的weights与bias，计算下一层的输入值；nn.W{1}为<100x785>，nn.W{2}为<10x101>，则nn.W{1}'为<785x100>，nn.W{2}为<101x10>
                         //nn.a{1}为<100x785>，nn.a{2}为<100x100>，这里只计算a{2}
                         //（《Machine Learning》中，nn.a{j}为<S{j} +1 x 1>，算上bias；但在这里一个样本表示为<1 x S{j} +1>，总共100个样本，因此nn.a{1}为<100 x 785>，nn.a{2}为<100x100>）
                        // 在《Machine Learning》中，a{j+1}<S{j+1} x 1> = Theta{j}'<S{j+1} x S{j}+1> * a{j} <S{j}+1 x 1>表示一个样本的计算；在这里一个样本的计算为a{j+1}<1 x S{j+1} > = a{j}<1 x S{j}+1> * nn.W{j}'<S{j}+1 x S{j+1}>，a{j}为样本输入形式为行向量，W{j}的列向量为输入样本的权值计算，得到结果行向量a{j+1}<1 x S{j+1} > ，即a{2}<1 x 100>，第二层有100个神经元，100个样本，则 a{2}<100 x 100>
                    case 'tanh_opt'
                         nn.a{i} = tanh_opt( nn.a{i - 1} * nn.W{i - 1}' )
            end
            
            //dropout
            if( nn.dropoutFraction > 0 )
                    if ( nn.testing )
                            nn.a{i} = nn.a{i} .* ( 1 - nn.dropoutFraction )
                    else
                            nn.dropoOutMask{i} = ( rand( size( nn.a{i} ) ) > nn.dropoutFraction );      
                            //随机得到与样本大小相同的DropOutMask，根据dropoutFraction参数，取Mask中值，其中dropOutMask为二值，非0即1
                            // 此例中不使用dropout；rand( 100 , 785)，dropOutMask<100x785>，满足条件则取1，否则取0
                            nn.a{i} = nn.a{i} .* nn.dropoOutMask{i};                //根据上值取舍样本
                            
                    end
            end
            
            // 如果有sparsity则计算，calculate running exponential activations for use with sparsity
            if(nn.nonSparsityPenalty>0)
                    nn.p{i} = 0.99 * nn.p{i} + 0.01 * mean(nn.a{i}, 1);
            end
            
            nn.a{i} = [ones(m,1) nn.a{i}];          //为nn.a{i}添加bias项，由于在此例中，只添加nn.a{2}<100 x 1+100>
    end
    
    switch nn.output
            case 'sigm'
                    nn.a{n} = sigm( nn.a{n-1} * nn.W{n-1}' );       
                    //计算nn.a{3}<100x10>=nn.a{2}<100x101>*nn.W{2}'<101 x 10>，其中100行为100个样本，每个样本对应10个类别中的一个；W{2}<10 x 101>，则W{2}'<101 x 10>
                    //a{3}每行代表一个样本在输出层10个神经元的输出，为<1 x 10>总共100个样本，则nn.a{3}<100x10> 即 a{j}< m x S{j}>
            case 'linear'
                    nn.a{n} = nn.a{n-1} * nn.W{n-1}';
                    //线性计算nn.a{3}<100x10>
            case 'softmax'
                    nn.a{n} = nn.a{n-1} * nn.W{n-1}';
                    nn.a{n} = exp( bsxfun(@minus, nn.a{n}, max( nn.a{n}, [ ], 2 ) ) );
                    nn.a{n} = bsxfun( @rdivide, nn.a{n}, sum( nn.a{n}, 2 ) );
    end
    
    %计算error and loss
    nn.e = y - nn.a{n};             //batch_y<100x10>与nn.a{3<100x10>}的差值，nn.e<100x10>
    
    swith nn.output
            case {'sigm' , 'linear'}
                    nn.L = 1/2 * sum( sum( nn.e .^ 2 ) ) / m;   //sum(x)为竖向相加，求每列的和，结果为行向量
                    // sum( nn.e .^ 2)为行向量<1x100>，sum( sum( nn.e .^ 2 ))再求行向量的和，nn.L<1x1>为一值
            case 'softmax'
                    nn.L = -sum( sum( y .* log(nn.a{n} ) ) ) / m;
    end
end

function nn = nnbp( nn )
// NNBP执行backpropagation，nn = nnbp(nn)返回更新权值后的网络

     n = nn.n            //网络的层数
     sparsityError = 0;                  //稀疏参数
     switch nn.output                //根据具体output选择输出，这里output为sigm，
            case 'sigm'
                    d{n} = -nn.e .* ( nn.a{n} .* (1 - nn.a{n}));        // d{3}<100 x 10>相当于delta{3}  nn.e<100 x 10>  nn.a{3}<100 x 10>   g'(a{3}) <100 x 10>；在这里的计算与《Machine Learning》公式有差别
                    //《Machine Learning》中，error=a{3}-y，在这里nn.e=y-a{3}，因此加上负号，则-nn.e=a{3}-y；另外，《Machine Learning》中，d{3}=a{3}-y，相当于delta{3}，而在这里d{3}=(a{3}-y)<100 x 10> * a{3}<100 x 10> * (1-a{3})<100 x 10>；通常delta/{j}以列向量表示单样本误差<S{j} x 1>，而这里以行向量表示单样本误差<1 x S{j}>，即<1 x 10>，共100个样本，则d{3}=<100 x 10>
            case { 'softmax' , 'linear' }
                    d{n} = - nn.e;
    end

    //对activation函数求导
    for i = ( n-1 ) : -1 : 2
            switch nn.activation_function
                    case 'sigm'
                            d_act = nn.a{i} .* ( 1 - nn.a{i} );   //对激励函数求导g'(a{j})= a{j} * ( 1-a{j} )
                             //g'(a{2}) <100 x 101> = a{2} <100 x 101> .* (1-a{2}) <100 x 101>
                    case 'tanh_opt'
                            d_act = 1.7159 * 2/3 * (1 - 1/(1.7159)^2 * nn.a{i}.^2);
            end
        
            //对Sparsity的设置
            if( nn.nonSparsityPenalty>0 )
                    pi = repmat(nn.p{i}, size(nn.a{i}, 1), 1);
                    sparsityError = [zeros(size(nn.a{i},1),1) nn.nonSparsityPenalty * (-nn.sparsityTarget ./ pi + (1 - nn.sparsityTarget) ./ (1 - pi))];
            end
        
            //反向传播第一个导数
            if i+1 == n             //在本例中没有bias项，被移去了
                    d{i} = ( d{i+1} * nn.W{i} + sparsityError ) .* d_act;   
                    // d{2}<100 x 101> = d{3}<100 x 10> * W{2}<10 x 101> .* d_act<100 x 101>( 即g'(a{2}) ) <100 x 100>
                    //由于该例子只有3层，d{2}是倒数第二层的误差，是由最后一层误差d{3}得到的，由于a{3}无bias，a{2}的bias参与计算最后一层的输出，在此需要计算
                    //a{2}与W{2}的bias只计算下一层的权值W{2}与本层的误差d{2}，与上一层的神经元无连接，不参与计算上一层的权值W{1}与误差d{1}
                    //《Machine Learning》中，误差delta是列向量表示，delta{j}<S{j}+1 x 1>=[W{j}]T<S{j}+1 x S{j+1}> * delta{j+1}<S{j+1} x 1> .* g'(z{j}) <S{j}+1 x 1>，其中delta{j+1}去除了bias项；在这里，误差d是由行向量进行表示，d{j}<1 x S{j}+1> = d{j+1}<1 x S{j+1}> * W{j}<S{j+1} x S{j}+1>，其中d{j+1}去除了bias项；d{2}<1 x 101> = d{3}<1 x 10> * W{2}<10 x 101> .*g'(a{2})< 1 x 101>，总共100个样本，则d{2}为<100 x 101> 
            else
                    d{i} = ( d{i+1}( : , 2:end) *  nn.W{i} + sparsityError) .* d_act;
                    // d{i+1}(: , 2:end)表示除去bias项所得的梯度，以计算d{1}为例（d{1}是不需要计算的）d{1}<100 x 785> = d{2}<100 x 100> * W{1}<100 x 785> .* (a{1} <100 x 785> .* (1-a{1})<100 x 785>)，其中d{2}去掉了第一列即每个样本的bias，因为其不是由上一层神经元与权值计算所得，则不参与上一层的反向计算
                    // d{3}<100 x 10>，d{2}<100 x 101>
            end
        
            //dropout项，本例无
            if(nn.dropoutFraction>0)
                    d{i} = d{i} .* [ ones( size( d{i), 1 ),1) nn.dropOutMask{i}]
            end
    end    
    
    for i = 1 :  ( n - 1)
            if i+1 == n
                    nn.dW{i} =  = (d{i + 1}' * nn.a{i}) / size(d{i + 1}, 1);
                    // dW{2}为倒数第二层，由于最后一层无bias，则不用对delta的bias项处理
                    // 《Machine Learning》中，梯度计算为D{l}<S{l+1} x S{l}+1>=d{l+1}<S{l+1} x 1> * [a{l}]T<1 x S{l}+1>，这为一个样本的梯度计算方法
                    //  dW{i}<S{i+1} x S{i}+1> = d{i+1}' <S{i+1} x m> *  a{i} <m x S{i}+1> 经过转置后，d{i+1}与a{i}与《Machine Learning》公式相同，但样本总数为m个，因此dW{i}结果为m个样本的梯度相加和，因此需要除以样本数m，则为平均梯度，即最终结果
            else
                    nn.dW{i} = (d{i + 1}(:,2:end)' * nn.a{i}) / size(d{i + 1}, 1);
                    //对误差d{i+1}出去bias项计算平均梯度
            end
    end
end               

function nn = nnapplygrads(nn)
//NN用计算好的梯度更新 weights和 biases，返回更新好权值和偏差的网络

    for i = 1 : (nn.n - 1)
            if(nn.weightPenaltyL2>0)
                    dW = nn.dW{i} + nn.weightPenaltyL2 * [zeros(size(nn.W{i},1),1) nn.W{i}(:,2:end)];
                    //当用L2正则项时，使用上述计算
            else
                    dW = nn.dW{i};
                    //通常情况下，dW{1}和dw{2}就表示所求梯度，即 J(Theta)的求导
            end
            
            dW = nn.learningRate * dW;
            
            if(nn.momentum>0)
                    nn.vW{i} = nn.momentum*nn.vW{i} + dW;
                    // nW{1}<100 x 785> nW{2}<10 x 101>引入动量计算
                    dW = nn.vW{i};
            end
            
            nn.W{i} = nn.W{i} - dW;
            //更新weights和bias，使用theta := theta - alpha*d( J(theta) )
    end
end     

function [loss] = nneval(nn, loss, train_x, train_y, val_x, val_y)
//NNEVAL对网络进行了评价，返回更新后的loss
    
    assert(nargin == 4 || nargin == 6, 'Wrong number of arguments');           
    nn.testing = 1;
    
    nn = nnff(nn, train_x, train_y);    // 用全部训练集来跑一次网络，不训练，只输出错误率结果
    loss.train.e( end + 1) = nn.L;
    //即loss.train.e{1} = nn.L；使用end+1表示记录新的错误率，且只输出当前最新错误率
    
    if nargin == 6      //当有验证集存在时，再用验证集跑一次网络，输出错误率
            nn = nnff( nn, val_x, val_y);
            loss.val.e( end + 1) = nn.L;
    end
    
    nn.testing = 0;
    
    if strcmp(nn.output,'softmax')
            [er_train, dummy] = nntest(nn, train_x, train_y);
            loss.train.e_frac(end+1) = er_train;
            
            if nargin == 6
                    [er_val, dummy] = nntest(nn, val_x, val_y);
                    loss.val.e_frac(end+1) = er_val;
            end
    end
end

function [er, bad] = nntest(nn, x, y)
    labels = nnpredict(nn, x);
    //用nn网络对输入的test_x进行预测，得到预测值labels
    
    [dummy, expected] = max(y, [], 2);
    //dim为2时，取nn.a{end}中每行的最大值，其中dummy<10000 x 1>存储test_y每行的最大值，即数值全为1；expected<10000 x 1>储存每行最大值的位置，即确定类别
    
    bad = find(lables ~= expected);
    // 找到预测值与实际值不相同的样本，返回列向量<n x 1>其中存索引值
    
    er = numel(bad) / size(x,1);
    // er表示错误率，即错误预测数目 / 测试样本数目
end

function labels = nnpredict(nn, x)
    nn.testing = 1;
    nn = nnff(nn, x, zeros(size(x,1), nn.size(end)));
    // 用test_x与标签为0的y跑一次网络，zeros(size(x,1), nn.size(end))为zeros(10000, 10)，根据网络最后一层输出变化
    nn.testing = 0;
    // 测试完成，改变testing flag
    
    [dummy, i] = max(nn.a{end},[],2);
    //dim为2时，取nn.a{end}中每行的最大值，其中dummy<10000 x 1>存储a{3}每行最大值；i<10000 x 1>存储a{3}中每行最大值的位置，数值范围在0-9之间
    
    labels = i;    
    // 根据a{3}对各类的判断取最大值，即为其预测值
end