donutDim = floor(sqrt(numFeatures));

assert(donutDim * donutDim == numFeatures,'donutDim^2 must be equal to numFeatures');

groupMatrix = zeros(numFeatures, donutDim, donutDim);

groupNum = ;

for row = :donutDim

    for col = :donutDim

        groupMatrix(groupNum, :poolDim, :poolDim) = ;

        groupNum = groupNum + ;

        groupMatrix = circshift(groupMatrix, [  -]);

    end

    groupMatrix = circshift(groupMatrix, [ -, ]);

end

%poolDim = ;

%% CS294A/CS294W Sparse Coding Exercise

%  Instructions

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

%

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

%  sparse coding exercise. In this exercise, you will need to modify

%  sparseCodingFeatureCost.m and sparseCodingWeightCost.m. You will also

%  need to modify this file, sparseCodingExercise.m slightly.

% Add the paths to your earlier exercises if necessary

% addpath /path/to/solution

addpath minFunc/;

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

%% STEP : Initialization

%  Here we initialize some parameters used for the exercise.

numPatches = ;   % number of patches

numFeatures = ;    % number of features to learn

patchDim = ;         % patch dimension

visibleSize = patchDim * patchDim; 

% dimension of the grouping region (poolDim x poolDim) for topographic sparse coding

poolDim = ;

% number of patches per batch

batchNumPatches = ; 

lambda = 5e-;  % L1-regularisation parameter (on features)

epsilon = 1e-; % L1-regularisation epsilon |x| ~ sqrt(x^ + epsilon)

gamma = 1e-;   % L2-regularisation parameter (on basis)

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

%% STEP : Sample patches

images = load('IMAGES.mat');

images = images.IMAGES;

patches = sampleIMAGES(images, patchDim, numPatches);

display_network(patches(:, :));

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

%% STEP : Implement and check sparse coding cost functions

%  Implement the two sparse coding cost functions and check your gradients.

%  The two cost functions are

%  ) sparseCodingFeatureCost (in sparseCodingFeatureCost.m) for the features

%     (used when optimizing for s, which is called featureMatrix in this exercise)

%  ) sparseCodingWeightCost (in sparseCodingWeightCost.m) for the weights

%     (used when optimizing for A, which is called weightMatrix in this exericse)

% We reduce the number of features and number of patches for debugging

debug = false;

if debug

numFeatures = ;

patches = patches(:, :);

numPatches = ;

weightMatrix = randn(visibleSize, numFeatures) * 0.005;

featureMatrix = randn(numFeatures, numPatches) * 0.005;

%% STEP 2a: Implement and test weight cost

%  Implement sparseCodingWeightCost in sparseCodingWeightCost.m and check

%  the gradient

[cost, grad] = sparseCodingWeightCost(weightMatrix, featureMatrix, visibleSize, numFeatures, patches, gamma, lambda, epsilon);

numgrad = computeNumericalGradient( @(x) sparseCodingWeightCost(x, featureMatrix, visibleSize, numFeatures, patches, gamma, lambda, epsilon), weightMatrix(:) );

% Uncomment the blow line to display the numerical and analytic gradients side by side

% disp([numgrad grad]);

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

fprintf('Weight difference: %g\n', diff);

assert(diff < 1e-, 'Weight difference too large. Check your weight cost function. ');

%% STEP 2b: Implement and test feature cost (非拓扑结构non-topographic)

%  Implement sparseCodingFeatureCost in sparseCodingFeatureCost.m and check

%  the gradient. You may wish to implement the non-topographic version of

%  the cost function first, and extend it to the topographic version later.

% Set epsilon to a larger value so checking the gradient numerically makes sense

epsilon = 1e-;

% We pass in the identity matrix as the grouping matrix, putting each

% feature in a group on its own.

groupMatrix = eye(numFeatures);

[cost, grad] = sparseCodingFeatureCost(weightMatrix, featureMatrix, visibleSize, numFeatures, patches, gamma, lambda, epsilon, groupMatrix);

numgrad = computeNumericalGradient( @(x) sparseCodingFeatureCost(weightMatrix, x, visibleSize, numFeatures, patches, gamma, lambda, epsilon, groupMatrix), featureMatrix(:) );

% Uncomment the blow line to display the numerical and analytic gradients side by side

% disp([numgrad grad]);

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

fprintf('Feature difference (non-topographic): %g\n', diff);

assert(diff < 1e-, 'Feature difference too large. Check your feature cost function. ');

%% STEP 2c: Implement and test feature cost (拓扑结构topographic)

%  Implement sparseCodingFeatureCost in sparseCodingFeatureCost.m and check

%  the gradient. This time, we will pass a random grouping matrix in to

%  check if your costs and gradients are correct for the topographic

%  version.

% Set epsilon to a larger value so checking the gradient numerically makes sense

epsilon = 1e-;

% This time we pass in a random grouping matrix to check if the grouping is

% correct.

groupMatrix = rand(, numFeatures);

[cost, grad] = sparseCodingFeatureCost(weightMatrix, featureMatrix, visibleSize, numFeatures, patches, gamma, lambda, epsilon, groupMatrix);

numgrad = computeNumericalGradient( @(x) sparseCodingFeatureCost(weightMatrix, x, visibleSize, numFeatures, patches, gamma, lambda, epsilon, groupMatrix), featureMatrix(:) );

% Uncomment the blow line to display the numerical and analytic gradients side by side

% disp([numgrad grad]);

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

fprintf('Feature difference (topographic): %g\n', diff);

assert(diff < 1e-, 'Feature difference too large. Check your feature cost function. ');

end

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

%% STEP : Iterative optimization

%  Once you have implemented the cost functions, you can now optimize for

%  the objective iteratively. The code to do the iterative optimization

%  using mini-batching and good initialization of the features has already

%  been included for you.

%

%  However, you will still need to derive and fill in the analytic solution

%  for optimizing the weight matrix given the features.

%  Derive the solution and implement it in the code below, verify the

%  gradient as described in the instructions below, and then run the

%  iterative optimization.

% Initialize options for minFunc

options.Method = 'lbfgs';

options.display = 'off';

options.verbose = ;

% Initialize matrices

weightMatrix = rand(visibleSize, numFeatures);

featureMatrix = rand(numFeatures, batchNumPatches);

% Initialize grouping matrix

assert(floor(sqrt(numFeatures)) ^ == numFeatures, 'numFeatures should be a perfect square');

donutDim = floor(sqrt(numFeatures));

assert(donutDim * donutDim == numFeatures,'donutDim^2 must be equal to numFeatures');

groupMatrix = zeros(numFeatures, donutDim, donutDim);

groupNum = ;

for row = :donutDim

    for col = :donutDim

        groupMatrix(groupNum, :poolDim, :poolDim) = ;

        groupNum = groupNum + ;

        groupMatrix = circshift(groupMatrix, [  -]);

    end

    groupMatrix = circshift(groupMatrix, [ -, ]);

end

groupMatrix = reshape(groupMatrix, numFeatures, numFeatures);

if isequal(questdlg('Initialize grouping matrix for topographic or non-topographic sparse coding?', 'Topographic/non-topographic?', 'Non-topographic', 'Topographic', 'Non-topographic'), 'Non-topographic')

    groupMatrix = eye(numFeatures);

end

% Initial batch

indices = randperm(numPatches);

indices = indices(:batchNumPatches);

batchPatches = patches(:, indices);                           

fprintf('%6s%12s%12s%12s%12s\n','Iter', 'fObj','fResidue','fSparsity','fWeight');

for iteration = :

    error = weightMatrix * featureMatrix - batchPatches;

    error = sum(error(:) .^ ) / batchNumPatches;

    fResidue = error;

    R = groupMatrix * (featureMatrix .^ );

    R = sqrt(R + epsilon);

    fSparsity = lambda * sum(R(:));    

    fWeight = gamma * sum(weightMatrix(:) .^ );

    fprintf('  %4d  %10.4f  %10.4f  %10.4f  %10.4f\n', iteration, fResidue+fSparsity+fWeight, fResidue, fSparsity, fWeight)

    % Select a new batch

    indices = randperm(numPatches);

    indices = indices(:batchNumPatches);

    batchPatches = patches(:, indices);                    

    % Reinitialize featureMatrix with respect to the new batch

    featureMatrix = weightMatrix' * batchPatches;

    normWM = sum(weightMatrix .^ )';

    featureMatrix = bsxfun(@rdivide, featureMatrix, normWM); 

    % Optimize for feature matrix

    options.maxIter = ;

    [featureMatrix, cost] = minFunc( @(x) sparseCodingFeatureCost(weightMatrix, x, visibleSize, numFeatures, batchPatches, gamma, lambda, epsilon, groupMatrix), ...

                                           featureMatrix(:), options);

    featureMatrix = reshape(featureMatrix, numFeatures, batchNumPatches);                                      

    % Optimize for weight matrix

%     weightMatrix = zeros(visibleSize, numFeatures);

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

    % Instructions:

    %   Fill in the analytic solution for weightMatrix that minimizes

    %   the weight cost here.

    %   Once that is done, use the code provided below to check that your

    %   closed form solution is correct.

    %   Once you have verified that your closed form solution is correct,

    %   you should comment out the checking code before running the

    %   optimization.

    weightMatrix = batchPatches * featureMatrix' / (featureMatrix * featureMatrix' + gamma * batchNumPatches * eye(numFeatures));

%     [cost, grad] = sparseCodingWeightCost(weightMatrix, featureMatrix, visibleSize, numFeatures, batchPatches, gamma, lambda, epsilon, groupMatrix);

%     assert(norm(grad) < 1e-, 'Weight gradient should be close to 0. Check your closed form solution for weightMatrix again.')

%     error('Weight gradient is okay. Comment out checking code before running optimization.');

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

    % Visualize learned basis

    figure();

    display_network(weightMatrix);

end

function patches = sampleIMAGES(IMAGES, patchsize, numpatches)

% sampleIMAGES

% Returns  patches for training

% load IMAGES;    % load images from disk

%

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

% numpatches = ;

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

% column per patch,  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  images

%  For instance, IMAGES(:,:,) 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(:,:,) is an image

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

%  Image

tic

image_size=size(IMAGES);

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

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

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

for num=:numpatches

        patches(:,num)=reshape(IMAGES(i(num):i(num)+patchsize-,j(num):j(num)+patchsize-,k(num)),,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 [,]

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

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

% 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 +/- standard deviations and scale to - to

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

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

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

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

patches = (patches + ) * 0.4 + 0.1;

end

function [cost, grad] = sparseCodingWeightCost(weightMatrix, featureMatrix, visibleSize, numFeatures,  patches, gamma, lambda, epsilon, groupMatrix)

%sparseCodingWeightCost - given the features in featureMatrix,

%                         computes the cost and gradient with respect to

%                         the weights, given in weightMatrix

% parameters

%   weightMatrix  - the weight matrix. weightMatrix(:, c) is the cth basis

%                   vector.

%   featureMatrix - the feature matrix. featureMatrix(:, c) is the features

%                   for the cth example

%   visibleSize   - number of pixels in the patches

%   numFeatures   - number of features

%   patches       - patches

%   gamma         - weight decay parameter (on weightMatrix)

%   lambda        - L1 sparsity weight (on featureMatrix)

%   epsilon       - L1 sparsity epsilon

%   groupMatrix   - the grouping matrix. groupMatrix(r, :) indicates the

%                   features included in the rth group. groupMatrix(r, c)

%                   is  if the cth feature is in the rth group and

%                   otherwise.

    if exist('groupMatrix', 'var')

        assert(size(groupMatrix, ) == numFeatures, 'groupMatrix has bad dimension');

    else

        groupMatrix = eye(numFeatures);

    end

    numExamples = size(patches, );

    weightMatrix = reshape(weightMatrix, visibleSize, numFeatures);

    featureMatrix = reshape(featureMatrix, numFeatures, numExamples);

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

    % Instructions:

    %   Write code to compute the cost and gradient with respect to the

    %   weights given in weightMatrix.

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

 %% 求目标的代价函数

    delta = weightMatrix*featureMatrix-patches;

    fResidue = sum(sum(delta.^))./numExamples;%重构误差

    fWeight = gamma*sum(weightMatrix(:).^);%防止基内元素值过大

    sparsityMatrix = sqrt(groupMatrix*(featureMatrix.^)+epsilon);

    fSparsity = lambda*sum(sparsityMatrix(:));% 对特征系数性的惩罚值

    cost = fResidue+fWeight+fSparsity; %目标的代价函数

%     cost = fResidue+fWeight;

    %% 求目标代价函数的偏导函数

    grad = (*weightMatrix*featureMatrix*(featureMatrix')-2*patches*featureMatrix')./numExamples+*gamma*weightMatrix;

    grad = grad(:);

end

function [cost, grad] = sparseCodingFeatureCost(weightMatrix, featureMatrix, visibleSize, numFeatures, patches, gamma, lambda, epsilon, groupMatrix)

%sparseCodingFeatureCost - given the weights in weightMatrix,

%                          computes the cost and gradient with respect to

%                          the features, given in featureMatrix

% parameters

%   weightMatrix  - the weight matrix. weightMatrix(:, c) is the cth basis

%                   vector.

%   featureMatrix - the feature matrix. featureMatrix(:, c) is the features

%                   for the cth example

%   visibleSize   - number of pixels in the patches

%   numFeatures   - number of features

%   patches       - patches

%   gamma         - weight decay parameter (on weightMatrix)

%   lambda        - L1 sparsity weight (on featureMatrix)

%   epsilon       - L1 sparsity epsilon

%   groupMatrix   - the grouping matrix. groupMatrix(r, :) indicates the

%                   features included in the rth group. groupMatrix(r, c)

%                   is  if the cth feature is in the rth group and

%                   otherwise.

    if exist('groupMatrix', 'var')

        assert(size(groupMatrix, ) == numFeatures, 'groupMatrix has bad dimension');

    else

        groupMatrix = eye(numFeatures);

    end

    numExamples = size(patches, );

    weightMatrix = reshape(weightMatrix, visibleSize, numFeatures);

    featureMatrix = reshape(featureMatrix, numFeatures, numExamples);

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

    % Instructions:

    %   Write code to compute the cost and gradient with respect to the

    %   features given in featureMatrix.

    %   You may wish to write the non-topographic version, ignoring

    %   the grouping matrix groupMatrix first, and extend the

    %   non-topographic version to the topographic version later.

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

      %% 求目标的代价函数

    delta = weightMatrix*featureMatrix-patches;

    fResidue = sum(sum(delta.^))./numExamples;%重构误差

    fWeight = gamma*sum(weightMatrix(:).^);%防止基内元素值过大

    sparsityMatrix = sqrt(groupMatrix*(featureMatrix.^)+epsilon);

    fSparsity = lambda*sum(sparsityMatrix(:)); %对特征系数性的惩罚值

    cost = fResidue++fSparsity+fWeight;%此时A为常量，可以不用

%     cost = fResidue++fSparsity;

    %% 求目标代价函数的偏导函数

    gradResidue = (-*weightMatrix'*patches+2*(weightMatrix')*weightMatrix*featureMatrix)./numExamples;

    % 非拓扑结构时：

    isTopo = ;  %isTopo = 1时，表示为拓扑结构

    if ~isTopo

        gradSparsity = lambda*(featureMatrix./sparsityMatrix);

    % 拓扑结构时

    else

        gradSparsity = lambda*groupMatrix'*(groupMatrix*(featureMatrix.^2)+epsilon).^(-0.5).*featureMatrix;%一定要小心最后一项是内积乘法

    end

    grad = gradResidue+gradSparsity;

    grad = grad(:);

end