function fit = multislice(data_dir,sparse_file,tucker_file,I,J)
%MULTISLICE is a low RAM Tucker decomposition
%
%   Peter Turney
%   October 26, 2007
%
%   Copyright 2007, National Research Council of Canada
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see http://www.gnu.org/licenses/.
%
%%  set parameters
%
fprintf('MULTISLICE is running ...\n');
%
maxloops = 50;           % maximum number of iterations
eigopts.disp = 0;        % suppress messages from eigs()
minfitchange = 1e-4;     % minimum change in fit of tensor
%
%%  make slices of input data file
%
fprintf('   preparing slices\n');
%
mode1_dir = 'slice1';
mode2_dir = 'slice2';
mode3_dir = 'slice3';
%
slice(data_dir,sparse_file,mode1_dir,1,I);
slice(data_dir,sparse_file,mode2_dir,2,I);
slice(data_dir,sparse_file,mode3_dir,3,I);
%
%%  pseudo HO-SVD initialization
%
%   initialize B
%
M2 = zeros(I(2),I(2));
for i = 1:I(3)
    X3_slice = load_slice(data_dir,mode3_dir,i);
    M2 = M2 + (X3_slice' * X3_slice);
end
for i = 1:I(1)
    X1_slice = load_slice(data_dir,mode1_dir,i);
    M2 = M2 + (X1_slice * X1_slice');
end
[B,D] = eigs(M2*M2',J(2),'lm',eigopts);
%
%   initialize C
%
M3 = zeros(I(3),I(3));
for i = 1:I(1)
    X1_slice = load_slice(data_dir,mode1_dir,i);
    M3 = M3 + (X1_slice' * X1_slice);
end
for i = 1:I(2)
    X2_slice = load_slice(data_dir,mode2_dir,i);
    M3 = M3 + (X2_slice' * X2_slice);
end
[C,D] = eigs(M3*M3',J(3),'lm',eigopts);
%
%%  main loop
%
old_fit = 0;
%
fprintf('   entering main loop of MULTISLICE\n');
%
for loop_num = 1:maxloops
    %
    %   update A
    %
    M1 = zeros(I(1),I(1));
    for i = 1:I(2)
        X2_slice = load_slice(data_dir,mode2_dir,i);
        M1 = M1 + ((X2_slice * C) * (C' * X2_slice'));
    end
    for i = 1:I(3)
        X3_slice = load_slice(data_dir,mode3_dir,i);
        M1 = M1 + ((X3_slice * B) * (B' * X3_slice'));
    end
    [A,D] = eigs(M1*M1',J(1),'lm',eigopts);
    %
    %   update B
    %
    M2 = zeros(I(2),I(2));
    for i = 1:I(3)
        X3_slice = load_slice(data_dir,mode3_dir,i);
        M2 = M2 + ((X3_slice' * A) * (A' * X3_slice));
    end
    for i = 1:I(1)
        X1_slice = load_slice(data_dir,mode1_dir,i);
        M2 = M2 + ((X1_slice * C) * (C' * X1_slice'));
    end
    [B,D] = eigs(M2*M2',J(2),'lm',eigopts);
    %
    %   update C
    %
    M3 = zeros(I(3),I(3));
    for i = 1:I(1)
        X1_slice = load_slice(data_dir,mode1_dir,i);
        M3 = M3 + ((X1_slice' * B) * (B' * X1_slice));
    end
    for i = 1:I(2)
        X2_slice = load_slice(data_dir,mode2_dir,i);
        M3 = M3 + ((X2_slice' * A) * (A' * X2_slice));
    end
    [C,D] = eigs(M3*M3',J(3),'lm',eigopts);
    %
    %   build the core
    %
    G = zeros(I(1)*J(2)*J(3),1);
    G = reshape(G,[I(1) J(2) J(3)]);
    for i = 1:I(1)
        X1_slice = load_slice(data_dir,mode1_dir,i);
        G(i,:,:) = B' * X1_slice * C;
    end
    G = reshape(G,[I(1) (J(2)*J(3))]);
    G = A' * G;
    G = reshape(G,[J(1) J(2) J(3)]);
    %
    %   measure fit
    %
    normX = 0;
    sqerr = 0;
    for i = 1:I(1)
        X1_slice = load_slice(data_dir,mode1_dir,i);
        X1_approx = reshape(G,[J(1) (J(2)*J(3))]);
        X1_approx = A(i,:) * X1_approx;
        X1_approx = reshape(X1_approx,[J(2) J(3)]);
        X1_approx = B * X1_approx * C';
        sqerr = sqerr + norm(X1_slice-X1_approx,'fro')^2;
        normX = normX + norm(X1_slice,'fro')^2;
    end
    fit = 1 - sqrt(sqerr) / sqrt(normX);
    %
    fprintf('   loop %d: fit = %f\n', loop_num, fit);
    %
    %   stop if fit is not increasing fast enough
    %
    if ((fit - old_fit) < minfitchange)
        break;
    end
    %
    old_fit = fit;
    %
end
%
fprintf('   total loops = %d\n', loop_num);
%
%%  save tensor
%
output_file = [data_dir, '/', tucker_file];
save(output_file,'G','A','B','C');
%
fprintf('   tucker tensor is in %s\n',tucker_file);
%
fprintf('MULTISLICE is done\n');
%