%% CS294A/CS294W Programming Assignment Starter Code

% Instructions

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

%

% This file contains code that helps you get started on the

% programming assignment. You will need to complete the code in sampleIMAGES.m,

% sparseAutoencoderCost.m and computeNumericalGradient.m.

% For the purpose of completing the assignment, you do not need to

% change the code in this file.

%

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

%% STEP 0: Here we provide the relevant parameters values that will

% allow your sparse autoencoder to get good filters; you do not need to

% change the parameters below.第0步：提供可得到较好滤波器的相关参数值，不得改变以下参数

visibleSize = 8*8; % number of input units 输入层单元数

hiddenSize = 25; % number of hidden units隐藏层单元数

sparsityParam = 0.01; % desired average activation of the hidden units.稀疏值

% (This was denoted by the Greek alphabet rho, which looks like a lower-case "p",

% in the lecture notes).

lambda = 0.0001; % weight decay parameter 权重衰减系数

beta = 3; % weight of sparsity penalty term稀疏值惩罚项的权重

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

%% STEP 1: Implement sampleIMAGES 第1步：实现图片采样

%

% 实现图片采样后，函数display_network从训练集中随机显示200张

% After implementing sampleIMAGES, the display_network command should

% display a random sample of 200 patches from the dataset

patches = sampleIMAGES;

display_network(patches(:,randi(size(patches,2),200,1)),8);%从10000张中随机选择200张显示

% Obtain random parameters theta初始化参数向量theta

theta = initializeParameters(hiddenSize, visibleSize);

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

%% STEP 2: Implement sparseAutoencoderCost

%

% 在计算代价函数时，可以一次计算其所有的元素项值（均方差项、权重衰减项、惩罚项），但是一步一步地计算各元素项值，

% 然后每步完成后运行梯度检验的方法可能会更容易实现，建议按照下面的步骤来实现函数sparseAutoencoderCost：

% You can implement all of the components (squared error cost, weight decay term,

% sparsity penalty) in the cost function at once, but it may be easier to do

% it step-by-step and run gradient checking (see STEP 3) after each step. We

% suggest implementing the sparseAutoencoderCost function using the following steps:

%

% (a) Implement forward propagation in your neural network, and implement the

% squared error term of the cost function. Implement backpropagation to

% compute the derivatives. Then (using lambda=beta=0), run Gradient Checking

% to verify that the calculations corresponding to the squared error cost

% term are correct.实现神经网络中的前向传播和代价函数中的均方差项。通过反向传导计算偏导数。

% 然后运行梯度检验法来检查均方差项是否计算错误。

%

% (b) Add in the weight decay term (in both the cost function and the derivative

% calculations), then re-run Gradient Checking to verify correctness.在代价函数和偏导数计算中

% 加入权重衰减项，然后运行梯度检验法来检查其正确性。

%

% (c) Add in the sparsity penalty term, then re-run Gradient Checking to

% verify correctness.加入惩罚项，然后运行梯度检验法来检查其正确性。

%

% Feel free to change the training settings when debugging your

% code. (For example, reducing the training set size or

% number of hidden units may make your code run faster; and setting beta

% and/or lambda to zero may be helpful for debugging.) However, in your

% final submission of the visualized weights, please use parameters we

% gave in Step 0 above.

[cost, grad] = sparseAutoencoderCost(theta, visibleSize, hiddenSize, lambda, ...

sparsityParam, beta, patches);

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

%% STEP 3: Gradient Checking

%

% Hint: If you are debugging your code, performing gradient checking on smaller models

% and smaller training sets (e.g., using only 10 training examples and 1-2 hidden

% units) may speed things up.

% First, lets make sure your numerical gradient computation is correct for a

% simple function. After you have implemented computeNumericalGradient.m,

% run the following:

checkNumericalGradient();

% Now we can use it to check your cost function and derivative calculations

% for the sparse autoencoder.

numgrad = computeNumericalGradient( @(x) sparseAutoencoderCost(x, visibleSize, ...

hiddenSize, lambda, ...

sparsityParam, beta, ...

patches), theta);

% Use this to visually compare the gradients side by side

disp([numgrad grad]);

% Compare numerically computed gradients with the ones obtained from backpropagation

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

disp(diff); % Should be small. In our implementation, these values are

% usually less than 1e-9.

% When you got this working, Congratulations!!!

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

%% STEP 4: After verifying that your implementation of

% sparseAutoencoderCost is correct, You can start training your sparse

% autoencoder with minFunc (L-BFGS).

% Randomly initialize the parameters

theta = initializeParameters(hiddenSize, visibleSize);

% Use minFunc to minimize the function

addpath minFunc/

options.Method = 'lbfgs'; % Here, we use L-BFGS to optimize our cost

% function. Generally, for minFunc to work, you

% need a function pointer with two outputs: the

% function value and the gradient. In our problem,

% sparseAutoencoderCost.m satisfies this.

options.maxIter = 400;	 % Maximum number of iterations of L-BFGS to run

options.display = 'on';

[opttheta, cost] = minFunc( @(p) sparseAutoencoderCost(p, ...

visibleSize, hiddenSize, ...

lambda, sparsityParam, ...

beta, patches), ...

theta, options);

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

%% STEP 5: Visualization

W1 = reshape(opttheta(1:hiddenSize*visibleSize), hiddenSize, visibleSize);

display_network(W1', 12);

print -djpeg weights.jpg % save the visualization to a file

function patches = sampleIMAGES()

% sampleIMAGES

% Returns 10000 patches for training

load IMAGES;    % load images from disk 

patchsize = 8;  % we'll use 8x8 patches

numpatches = 10000;

% Initialize patches with zeros.  Your code will fill in this matrix--one

% column per patch, 10000 columns.

patches = zeros(patchsize*patchsize, numpatches);

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

%  Instructions: Fill in the variable called "patches" using data

%  from IMAGES.

%

%  IMAGES is a 3D array containing 10 images

%  For instance, IMAGES(:,:,6) is a 512x512 array containing the 6th image,

%  and you can type "imagesc(IMAGES(:,:,6)), colormap gray;" to visualize

%  it. (The contrast on these images look a bit off because they have

%  been preprocessed using using "whitening."  See the lecture notes for

%  more details.) As a second example, IMAGES(21:30,21:30,1) is an image

%  patch corresponding to the pixels in the block (21,21) to (30,30) of

%  Image 1

tic

image_size=size(IMAGES);

i=randi(image_size(1)-patchsize+1,1,numpatches);%生成元素值随机为大于0且小于image_size(1)-patchsize+1的1行numpatches矩阵

j=randi(image_size(2)-patchsize+1,1,numpatches);

k=randi(image_size(3),1,numpatches);

for num=1:numpatches

        patches(:,num)=reshape(IMAGES(i(num):i(num)+patchsize-1,j(num):j(num)+patchsize-1,k(num)),1,patchsize*patchsize);

end

toc

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

% For the autoencoder to work well we need to normalize the data

% Specifically, since the output of the network is bounded between [0,1]

% (due to the sigmoid activation function), we have to make sure

% the range of pixel values is also bounded between [0,1]

patches = normalizeData(patches);

end

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

function patches = normalizeData(patches)

% Squash data to [0.1, 0.9] since we use sigmoid as the activation

% function in the output layer

% Remove DC (mean of images). 把patches数组中的每个元素值都减去mean(patches)

patches = bsxfun(@minus, patches, mean(patches));

% Truncate to +/-3 standard deviations and scale to -1 to 1

pstd = 3 * std(patches(:));%把patches的标准差变为其原来的3倍

patches = max(min(patches, pstd), -pstd) / pstd;%因为根据3sigma法则，95%以上的数据都在该区域内

                                                % 这里转换后将数据变到了-1到1之间

% Rescale from [-1,1] to [0.1,0.9]

patches = (patches + 1) * 0.4 + 0.1;

end

function [cost,grad] = sparseAutoencoderCost(theta, visibleSize, hiddenSize, ...

                                             lambda, sparsityParam, beta, data)

% visibleSize:输入层神经单元节点数 the number of input units (probably 64)

% hiddenSize:隐藏层神经单元节点数 the number of hidden units (probably 25)

% lambda: 权重衰减系数weight decay parameter

% sparsityParam: 稀疏性参数The desired average activation for the hidden units (denoted in the lecture

%                           notes by the greek alphabet rho, which looks like a lower-case "p").

% beta: 稀疏惩罚项的权重weight of sparsity penalty term

% data: 训练集Our 64x10000 matrix containing the training data.  So, data(:,i) is the i-th training example.

% theta：参数向量，包含W1、W2、b1、b2

% The input theta is a vector (because minFunc expects the parameters to be a vector).

% We first convert theta to the (W1, W2, b1, b2) matrix/vector format, so that this

% follows the notation convention of the lecture notes. 

W1 = reshape(theta(1:hiddenSize*visibleSize), hiddenSize, visibleSize);

W2 = reshape(theta(hiddenSize*visibleSize+1:2*hiddenSize*visibleSize), visibleSize, hiddenSize);

b1 = theta(2*hiddenSize*visibleSize+1:2*hiddenSize*visibleSize+hiddenSize);

b2 = theta(2*hiddenSize*visibleSize+hiddenSize+1:end);

% Cost and gradient variables (your code needs to compute these values).

% Here, we initialize them to zeros.

cost = 0;

W1grad = zeros(size(W1));

W2grad = zeros(size(W2));

b1grad = zeros(size(b1));

b2grad = zeros(size(b2));

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

%  Instructions: Compute the cost/optimization objective J_sparse(W,b) for the Sparse Autoencoder,

%                and the corresponding gradients W1grad, W2grad, b1grad, b2grad.

%

% W1grad, W2grad, b1grad and b2grad should be computed using backpropagation.

% Note that W1grad has the same dimensions as W1, b1grad has the same dimensions

% as b1, etc.  Your code should set W1grad to be the partial derivative of J_sparse(W,b) with

% respect to W1.  I.e., W1grad(i,j) should be the partial derivative of J_sparse(W,b)

% with respect to the input parameter W1(i,j).  Thus, W1grad should be equal to the term

% [(1/m) \Delta W^{(1)} + \lambda W^{(1)}] in the last block of pseudo-code in Section 2.2

% of the lecture notes (and similarly for W2grad, b1grad, b2grad).

%

% Stated differently, if we were using batch gradient descent to optimize the parameters,

% the gradient descent update to W1 would be W1 := W1 - alpha * W1grad, and similarly for W2, b1, b2.

%

%1.前向传播forward propagation

data_size=size(data);

active_value2=repmat(b1,1,data_size(2));

active_value3=repmat(b2,1,data_size(2));

active_value2=sigmoid(W1*data+active_value2);                %第2层激活值

active_value3=sigmoid(W2*active_value2+active_value3);       %第3层激活值

%2.计算代价函数computing error term and cost

ave_square=sum(sum((active_value3-data).^2)./2)/data_size(2);%均方差项

weight_decay=lambda/2*(sum(sum(W1.^2))+sum(sum(W2.^2)));     %权重衰减项

p_real=sum(active_value2,2)./data_size(2); %对active_value2每行求和，再平均，得到每个隐藏单元的平均活跃度（25行1列）

p_para=repmat(sparsityParam,hiddenSize,1); %稀疏性参数（25行1列）

sparsity=beta.*sum(p_para.*log(p_para./p_real)+(1-p_para).*log((1-p_para)./(1-p_real)));%惩罚项

cost=ave_square+weight_decay+sparsity;     %代价函数

delta3=(active_value3-data).*(active_value3).*(1-active_value3); %第3层残差

average_sparsity=repmat(sum(active_value2,2)./data_size(2),1,data_size(2));%每个隐藏单元的平均活跃度（25行10000列）

default_sparsity=repmat(sparsityParam,hiddenSize,data_size(2));  %稀疏性参数（25行10000列）

sparsity_penalty=beta.*(-(default_sparsity./average_sparsity)+((1-default_sparsity)./(1-average_sparsity)));

delta2=(W2'*delta3+sparsity_penalty).*((active_value2).*(1-active_value2));%第2层残差

%3.反向传导backword propagation

W2grad=delta3*active_value2'./data_size(2)+lambda.*W2;  %W2梯度

W1grad=delta2*data'./data_size(2)+lambda.*W1;           %W1梯度

b2grad=sum(delta3,2)./data_size(2);                     %b2梯度

b1grad=sum(delta2,2)./data_size(2);                     %b1梯度

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

% After computing the cost and gradient, we will convert the gradients back

% to a vector format (suitable for minFunc).  Specifically, we will unroll

% your gradient matrices into a vector.

grad = [W1grad(:) ; W2grad(:) ; b1grad(:) ; b2grad(:)];

end

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

% Here's an implementation of the sigmoid function, which you may find useful

% in your computation of the costs and the gradients.  This inputs a (row or

% column) vector (say (z1, z2, z3)) and returns (f(z1), f(z2), f(z3)). 

function sigm = sigmoid(x)

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

end

function numgrad = computeNumericalGradient(J, theta)

% numgrad = computeNumericalGradient(J, theta)

% theta: 参数向量，包含W1、W2、b1、b2；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=0.0001;

for i=1:size(theta)

    theta_plus=theta;

    theta_minu=theta;

    theta_plus(i)=theta_plus(i)+EPSILON;

    theta_minu(i)=theta_minu(i)-EPSILON;

    numgrad(i)=(J(theta_plus)-J(theta_minu))/(2*EPSILON);

end

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

end