%% automatic detection of heart
clc;
clear all;
close all;
addpath('functions');
%% STEP 0:
% parameters
patchsize = 32;
visibleSize = patchsize*patchsize; % number of input units
hiddenSizeL1 = 100; % number of hidden units
hiddenSizeL2=100;
sparsityParam1 = 0.01; % desired average activation of the hidden units.
lambda = 3e-3; % weight decay parameter
beta = 3; % weight of sparsity penalty term
outputSize=visibleSize; % number of output units
%%======================================================================
%% STEP 1: laod training inputs and labels from mat files
load matFiles/training_data;
train_input=sampleIMAGES(t_I,patchsize);
train_labels=sampleIMAGES(t_yROI,patchsize);
%% train sparse Auto Encoder 1
% Randomly initialize the parameters
saeTheta1 = initializeParameters(hiddenSizeL1, 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';
[sae1OptTheta, cost] = minFunc( @(p) sparseAutoencoderCost(p, ...
visibleSize, hiddenSizeL1, ...
lambda, sparsityParam1, ...
beta, train_input), ...
saeTheta1, options);
%% STEP 5: Visualization of AE1
W1 = reshape(sae1OptTheta(1:hiddenSizeL1*visibleSize), hiddenSizeL1, visibleSize);
display_network(W1', 12);
%% compute activations from layer 1
[sae1Features] = feedForwardAutoencoder(sae1OptTheta, hiddenSizeL1, ...
visibleSize, train_input);
%% train sparse Auto Encoder 2
% Randomly initialize the parameters
sae2Theta = initializeParameters(hiddenSizeL2, hiddenSizeL1);
sparsityParam2=1e-1;
lambda2 = 3e-3; % weight decay parameter
beta2 = 3; % weight of sparsity penalty term
[sae2OptTheta, costL2] = minFunc( @(p) sparseAutoencoderCost(p, ...
hiddenSizeL1, hiddenSizeL2, ...
lambda2, sparsityParam2, ...
beta2, sae1Features), ...
sae2Theta, options);
W2 = reshape(sae2OptTheta(1:hiddenSizeL2*hiddenSizeL1), hiddenSizeL2, hiddenSizeL1);
%display_network(W2', 12);
%% compute activation from layer 2
[sae2Features] = feedForwardAutoencoder(sae2OptTheta, hiddenSizeL2, ...
hiddenSizeL1, sae1Features);
%% train multi outputs logstic regression
lambda_mr=1e-4;
options_mr.maxIter = 100;
trainLabels=train_labels;
mrModel = mrTrain(hiddenSizeL2, outputSize, lambda_mr, ...
sae2Features, trainLabels, options_mr);
saeMultRegOptTheta = mrModel.optTheta(:);
%% fine tuning
% Initialize the stack using the parameters learned
stack = cell(2,1);
inputSize=visibleSize;
stack{1}.w = reshape(sae1OptTheta(1:hiddenSizeL1*inputSize), ...
hiddenSizeL1, inputSize);
stack{1}.b = sae1OptTheta(2*hiddenSizeL1*inputSize+1:2*hiddenSizeL1*inputSize+hiddenSizeL1);
stack{2}.w = reshape(sae2OptTheta(1:hiddenSizeL2*hiddenSizeL1), ...
hiddenSizeL2, hiddenSizeL1);
stack{2}.b = sae2OptTheta(2*hiddenSizeL2*hiddenSizeL1+1:2*hiddenSizeL2*hiddenSizeL1+hiddenSizeL2);
% Initialize the parameters for the deep model
[stackparams, netconfig] = stack2params(stack);
stackedAETheta = [ saeMultRegOptTheta ; stackparams ];
[stackedAEOptTheta, loss] = minFunc( @(x) stackedAECost(x, ...
inputSize, hiddenSizeL2, outputSize, netconfig, ...
lambda, train_input, train_labels), ...
stackedAETheta, options);
%% save results
% save ROI_V32_H100_rho_0p01
%% test
% load test data
% load matFiles/validation_data;
% [xIsize,yIsize,zIsize]=size(t_I);
% test_input=sampleIMAGES(t_I,patchsize);
%
% [pyroi] = stackedAEPredict(stackedAEOptTheta, inputSize, hiddenSizeL2, ...
% outputSize, netconfig, test_input);
%
% yroi_h=reshape(pyroi,patchsize,patchsize,[]);
% %%
% for k=1:size(yroi_h,3)
% y1=yroi_h(:,:,k)';
% y1=imresize(y1,xIsize/patchsize);
% subplot(4,6,k)
% imshow(y1)
% I1=t_I(:,:,k);
% ymask(:,:,k)=y1.*I1;
% end
% disImgs(ymask,19);
