%% setup

clear all;

% load train data
data = loadMNISTImages('train-images-idx3-ubyte');
data = data(:,1:60000);

% params
global params;
% patchWidth 
params.patchWidth=28;
params.n=params.patchWidth^2;
% sparsity cost
params.lambda = 0.0005;
% number of features to extract
params.numFeatures = 100;
% epsilon used in the rica function
params.epsilon = 1e-2;   

% only if extracting patches
% data = sample_patches(data, params.patchWidth, 200000);

%% train RICA

%  Randomly initialize the parameters
randTheta = randn(params.numFeatures,params.n)*0.01;
randTheta = randTheta ./ repmat(sqrt(sum(randTheta.^2,2)), 1, size(randTheta,2)); 
randTheta = randTheta(:);

% minFunc params
paramsMinFunc.MaxIter = 1000;
paramsMinFunc.Method = 'lbfgs';
paramsMinFunc.MaxFunEvals = Inf;

% whiten the data
[data_whitened, V] = whiten(data);

% timer start
tic

% run RICA to find W
% use minFunc applied to ricaCost to extract optimum W
[W, cost, ~] = minFunc( @(theta) ricaCost(theta, data_whitened, params), randTheta, paramsMinFunc);

% timer stop
toc

% reshape to visualize
W = reshape(W, params.numFeatures, params.n);

% visualize
display_network(W');