%% CS294A/CS294W Linear Decoder Exercise

%  Instructions

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

%

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

%  linear decoder exericse. For this exercise, you will only need to modify

%  the code in sparseAutoencoderLinearCost.m. You will not need to modify

%  any code in this file.

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

%% STEP : Initialization

%  Here we initialize some parameters used for the exercise.

imageChannels = ;     % number of channels (rgb, so )

patchDim   = ;          % patch dimension

numPatches = ;   % number of patches

visibleSize = patchDim * patchDim * imageChannels;  % number of input units

outputSize  = visibleSize;   % number of output units

hiddenSize  = ;           % number of hidden units 

sparsityParam = 0.035; % desired average activation of the hidden units.

lambda = 3e-;         % weight decay parameter

beta = ;              % weight of sparsity penalty term       

epsilon = 0.1;           % epsilon for ZCA whitening

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

%% STEP : Create and modify sparseAutoencoderLinearCost.m to use a linear decoder,

%          and check gradients

%  You should copy sparseAutoencoderCost.m from your earlier exercise

%  and rename it to sparseAutoencoderLinearCost.m.

%  Then you need to rename the function from sparseAutoencoderCost to

%  sparseAutoencoderLinearCost, and modify it so that the sparse autoencoder

%  uses a linear decoder instead. Once that is done, you should check

% your gradients to verify that they are correct.

% NOTE: Modify sparseAutoencoderCost first!

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

% dummy patches

debugHiddenSize = ;

debugvisibleSize = ;

patches = rand([ ]);

theta = initializeParameters(debugHiddenSize, debugvisibleSize); 

[cost, grad] = sparseAutoencoderLinearCost(theta, debugvisibleSize, debugHiddenSize, ...

                                           lambda, sparsityParam, beta, ...

                                           patches);

% Check gradients

numGrad = computeNumericalGradient( @(x) sparseAutoencoderLinearCost(x, debugvisibleSize, debugHiddenSize, ...

                                                  lambda, sparsityParam, beta, ...

                                                  patches), 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');

% NOTE: Once your gradients check out, you should run step  again to

%       reinitialize the parameters

%}

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

%% STEP : 从pathes中学习特征 Learn features on small patches

%  In this step, you will use your sparse autoencoder (which now uses a

%  linear decoder) to learn features on small patches sampled from related

%  images.

%% STEP 2a: 加载数据 Load patches

%  In this step, we load 100k patches sampled from the STL10 dataset and

%  visualize them. Note that these patches have been scaled to [,]

load stlSampledPatches.mat  %怎么就就这个变量加到pathes上了呢？因为它里面自己定义了变量patches的值！

figure;

displayColorNetwork(patches(:, :)); 

%% STEP 2b: 预处理 Apply preprocessing

%  In this sub-step, we preprocess the sampled patches, in particular,

%  ZCA whitening them.

%

%  In a later exercise on convolution and pooling, you will need to replicate

%  exactly the preprocessing steps you apply to these patches before

%  using the autoencoder to learn features on them. Hence, we will save the

%  ZCA whitening and mean image matrices together with the learned features

%  later on.

% Subtract mean patch (hence zeroing the mean of the patches)

meanPatch = mean(patches, );  %为什么是对每行求平均，以前是对每列即每个样本求平均呀？因为以前是灰度图，现在是彩色图，如果现在对每列平均就是对三个通道求平均，这肯定不行

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

% Apply ZCA whitening

sigma = patches * patches' / numPatches; %协方差矩阵

[u, s, v] = svd(sigma);

ZCAWhite = u * diag( ./ sqrt(diag(s) + epsilon)) * u';

patches = ZCAWhite * patches;

figure;

displayColorNetwork(patches(:, :));

%% STEP 2c: Learn features

%  You will now use your sparse autoencoder (with linear decoder) to learn

%  features on the preprocessed patches. This should take around  minutes.

theta = initializeParameters(hiddenSize, visibleSize);

% Use minFunc to minimize the function

addpath minFunc/

options = struct;

options.Method = 'lbfgs';

options.maxIter = ;

options.display = 'on';

[optTheta, cost] = minFunc( @(p) sparseAutoencoderLinearCost(p, ...

                                   visibleSize, hiddenSize, ...

                                   lambda, sparsityParam, ...

                                   beta, patches), ...

                              theta, options);

% Save the learned features and the preprocessing matrices for use in

% the later exercise on convolution and pooling

fprintf('Saving learned features and preprocessing matrices...\n');

save('STL10Features.mat', 'optTheta', 'ZCAWhite', 'meanPatch');

fprintf('Saved\n');

%% STEP 2d: Visualize learned features

W = reshape(optTheta(:visibleSize * hiddenSize), hiddenSize, visibleSize);

b = optTheta(*hiddenSize*visibleSize+:*hiddenSize*visibleSize+hiddenSize);

figure;

displayColorNetwork( (W*ZCAWhite)');

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

                                                            lambda, sparsityParam, beta, data)

%计算线性解码器代价函数及其梯度

% visibleSize:输入层神经单元节点数

% hiddenSize:隐藏层神经单元节点数

% lambda: 权重衰减系数

% sparsityParam: 稀疏性参数

% beta: 稀疏惩罚项的权重

% data: 训练集

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

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

% Instructions:

%   Copy sparseAutoencoderCost in sparseAutoencoderCost.m from your

%   earlier exercise onto this file, renaming the function to

%   sparseAutoencoderLinearCost, and changing the autoencoder to use a

%   linear decoder.

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

% The input theta is a vector because minFunc only deal with vectors. In

% this step, we will convert theta to matrix format such that they follow

% the notation in the lecture notes.

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

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

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

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

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

m = size(data, ); % 样本数量

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

%  Instructions: Compute the loss for the Sparse Autoencoder and gradients

%                W1grad, W2grad, b1grad, b2grad

%

%  Hint: ) data(:,i) is the i-th example

%        ) your computation of loss and gradients should match the size

%        above for loss, W1grad, W2grad, b1grad, b2grad

% z2 = W1 * x + b1

% a2 = f(z2)

% z3 = W2 * a2 + b2

% h_Wb = a3 = f(z3)

z2 = W1 * data + repmat(b1, [, m]);

a2 = sigmoid(z2);

z3 = W2 * a2 + repmat(b2, [, m]);

a3 = z3;

rhohats = mean(a2,);

rho = sparsityParam;

KLsum = sum(rho * log(rho ./ rhohats) + (-rho) * log((-rho) ./ (-rhohats)));

squares = (a3 - data).^;

squared_err_J = (/) * (/m) * sum(squares(:));              %均方差项

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

sparsity_J = beta * KLsum;                                    %惩罚项

cost = squared_err_J + weight_decay_J + sparsity_J;%损失函数值

% delta3 = -(data - a3) .* fprime(z3);

% but fprime(z3) = a3 * (-a3)

delta3 = -(data - a3);

beta_term = beta * (- rho ./ rhohats + (-rho) ./ (-rhohats));

delta2 = ((W2' * delta3) + repmat(beta_term, [1,m]) ) .* a2 .* (1-a2);

W2grad = (/m) * delta3 * a2' + lambda * W2;   % W2梯度

b2grad = (/m) * sum(delta3, );               % b2梯度

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

b1grad = (/m) * sum(delta2, );               % b1梯度

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

% Convert weights and bias gradients to a compressed form

% This step will concatenate and flatten all your gradients to a vector

% which can be used in the optimization method.

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

end

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

% We are giving you the sigmoid function, you may find this function

% useful in your computation of the loss and the gradients.

function sigm = sigmoid(x)

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

end

function displayColorNetwork(A)

% display receptive field(s) or basis vector(s) for image patches

%

% A         the basis, with patches as column vectors

% In case the midpoint is not set at , we shift it dynamically

if min(A(:)) >=

    A = A - mean(A(:)); % 0均值化

end

cols = round(sqrt(size(A, )));% 每行大图像中小图像块的个数

channel_size = size(A,) / ;

dim = sqrt(channel_size);   % 小图像块内每行或列像素点个数

dimp = dim+;

rows = ceil(size(A,)/cols);   % 每列大图像中小图像块的个数

B = A(:channel_size,:);                   % R通道像素值

C = A(channel_size+:channel_size*,:);    % G通道像素值

D = A(*channel_size+:channel_size*,:);  % B通道像素值

B=B./(ones(size(B,),)*max(abs(B)));% 归一化

C=C./(ones(size(C,),)*max(abs(C)));

D=D./(ones(size(D,),)*max(abs(D)));

% Initialization of the image

I = ones(dim*rows+rows-,dim*cols+cols-,);

%Transfer features to this image matrix

for i=:rows-

  for j=:cols-

    if i*cols+j+ > size(B, )

        break

    end

    % This sets the patch

    I(i*dimp+:i*dimp+dim,j*dimp+:j*dimp+dim,) = ...

         reshape(B(:,i*cols+j+),[dim dim]);

    I(i*dimp+:i*dimp+dim,j*dimp+:j*dimp+dim,) = ...

         reshape(C(:,i*cols+j+),[dim dim]);

    I(i*dimp+:i*dimp+dim,j*dimp+:j*dimp+dim,) = ...

         reshape(D(:,i*cols+j+),[dim dim]);

  end

end

I = I + ; % 使I的范围从[-,]变为[，]

I = I / ; % 使I的范围从[，]变为[, ]

imagesc(I);

axis equal  % 等比坐标轴：设置屏幕高宽比，使得每个坐标轴的具有均匀的刻度间隔

axis off    % 关闭所有的坐标轴标签、刻度、背景

end