% Copyright (c) 2017, Domenico L. Gatti
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without 
% modification, are permitted provided that the following conditions are 
% met:
% 
%     * Redistributions of source code must retain the above copyright 
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright 
%       notice, this list of conditions and the following disclaimer in 
%       the documentation and/or other materials provided with the 
%       distribution
%       
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
% IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
%% General dependencies
% We always start from the CODE directory and we add to the path
% subdirectories containing various tools described in the book chapters.
addpath(genpath('../GENERAL_SCRIPTS_FUNCTIONS'));
addpath(genpath('../DATABASE'));

% In order to run the code of this chapter you will first need to install
% the following toolboxes in your TOOLBOXES directory:
% Tensor Toolbox/Poblano Toolbox: (http://www.sandia.gov/~tgkolda/TensorToolbox/index-2.6.html)
% Tensorlab Toolbox (http://www.tensorlab.net/)
% and add them to your path using 'addpath'. For example:
% addpath(genpath('../TOOLBOXES/tensorlab_2016-03-28'))
% addpath(genpath('../TOOLBOXES/tensor_toolbox'))
% addpath(genpath('../TOOLBOXES/poblano_toolbox'))

%% CHAPTER 21: Non-negative matrix factorization and tensor decomposition

% Non-negative matrix factorization + Sparse component analysis
clear, clc

X = randi(20,50,10)
openvar A, openvar D, openvar Bt

% NNMF

[ A,D,Bt,X_hat,niter,sse,sse_diff,sparseness ] = ...
    nnmf_sca(X,5);
sse

% Compare to MATLAB nnmf
[A2,Bt2] = nnmf(X,5);
X_hat2 = A2*Bt2;
R2 = X - X_hat2;
sse2 = R2(:)'*R2(:) 


[ A,D,Bt,X_hat,niter,sse,sse_diff,sparseness ] = ...
    nnmf_sca(X,5,'diag','nobalance');
sse

[ A,D,Bt,X_hat,niter,sse,sse_diff,sparseness ] = ...
    nnmf_sca(X,5,'diag','balance');
sse


% SCA
[ A,D,Bt,X_hat,niter,sse,sse_diff,sparseness ] = ...
    nnmf_sca(X,5,'ident','nobalance','sparse');
sse

[ A,D,Bt,X_hat,niter,sse,sse_diff,tot_sparseness,col_sparseness ] = ...
    nnmf_sca(X,5,'ident','nobalance','sparse',0.3);
sse

[ A,D,Bt,X_hat,niter,sse,sse_diff,tot_sparseness,col_sparseness ] = ...
    nnmf_sca(X,5,'ident','nobalance','sparse',0.3,'bycols');
sse

% NNMF + SCA
[ A,D,Bt,X_hat,niter,sse,sse_diff,tot_sparseness,col_sparseness ] = ...
    nnmf_sca(X,5,'ident','nobalance','both',0.3,'bycols');
sse

[ A,D,Bt,X_hat,niter,sse,sse_diff,tot_sparseness,col_sparseness ] = ...
    nnmf_sca(X,5,'diag','balance','both',0.3,'random');
sse


% Example with more variables than observations
clear, clc

X = randi(20,30,50)
openvar A, openvar D, openvar Bt

[ A,D,Bt,X_hat,niter,sse,sse_diff,col_sparseness,tot_sparseness ] = ...
    nnmf_sca(X,5);
sse

% Compare to MATLAB nnmf
[A2,Bt2] = nnmf(X,5);
X_hat2 = A2*Bt2;
R2 = X - X_hat2;
sse2 = R2(:)'*R2(:) 


%% Example 1 
clear, clc, close all
load fisheriris

[ A,D,Bt,X_hat,niter,sse ] = nnmf_sca(meas,2);
sse
[A2,Bt2] = nnmf(meas,2);
X = meas;
X_hat2 = A2*Bt2;
R2 = X - X_hat2;
sse2 = R2(:)'*R2(:) 

biplot(Bt','varlabels',{'sl','sw','pl','pw'},'scores',A);
xlabel('Column 1')
ylabel('Column 2')

%% Example 2
clear, clc, close all

load moore
X = moore(:,1:5);

% Use nnmf_sca
[ A,D,Bt,X_hat,niter,sse ] = nnmf_sca(X,2);
sse

figure
biplot(Bt','scores',A,'varlabels',{'','','v3','','v5'});
xlabel('Column 1')
ylabel('Column 2')

% Use MATLAB nnmf
[W,H] = nnmf(X,2);
X_hat2 = W*H;
R2 = X - X_hat2;
sse2 = R2(:)'*R2(:) 

figure
biplot(H','scores',W,'varlabels',{'','','v3','','v5'});
xlabel('Column 1')
ylabel('Column 2')

%% Tensors
addpath(genpath('../TOOLBOXES/tensorlab_2016-03-28'))
addpath(genpath('../TOOLBOXES/tensor_toolbox'))
clear, clc, close all

% Generate a random tensor:
T = randn(10,20,30);
S = randn(10,20,30);
V = T + S; % Elementwise addition
W = T.*S; % Elementwise multiplication

% Generate an incomplete tensor:
T = NaN(9,9,9);
T(1,2,3) = -1;
T(4,5,6) = -2;
T(7,8,9) = -3;
T = fmt(T)
Tf = ful(T)

% Sandia Toolbox
Ts = tensor(Tf)

% Generate a sparse tensor:
T = zeros(5,5,5);
T(1:31:end) = 1;

Ts = sptensor(T)    % Sandia
T = fmt(T)          % Tensorlab
Tf = ful(T)
size(T)

Ts = tensor(Tf); % Sandia

% or:
clear T
T.val = [-1 -2 -3 ]; 
T.sub = [1 2 3 ; 4 5 6 ; 7 8 9 ]; 
T.size = [9 9 9]; 
T.sparse = true; 
T = fmt(T);
ful(T)

% Matricization and tensorization
clear T
T = randn(3,5,7,9);
T = fmt(T)
M = tens2mat(T,[1 3],[4 2]);
size(T)
size(M)

% Sandia Toolbox
Ms = tenmat(T,[1 3],[4 2]);

% mode-n matricization
T = zeros(3,3,3,3);
for i = 1:3
    for j = 1:3
        for k = 1:3
            for l = 1:3
                T(i,j,k,l) = l+10*k+100*j+1000*i;
            end
        end
    end
end

T = fmt(T)
Ts = tensor(T)

M = tens2mat(T,[1 2],[3 4]); 
size_tens = size(T)
size(M)
M_nrows = prod(size_tens([1 2]))
M_ncols = prod(size_tens([3 4]))

M1 = tens2mat(T,1);
% M1 = tens2mat(T,1,[2 3 4]);
M2 = tens2mat(T,2);
% M2 = tens2mat(T,2,[1 3 4]);
M3 = tens2mat(T,3);
% M3 = tens2mat(T,3,[1 2 4]);
M4 = tens2mat(T,4);
% M4 = tens2mat(T,4,[1 2 3]);
size(M4)

% Sandia Toolbox
Ms = tenmat(T,[1 3],[4 2]);
%
M1s = tenmat(T,1);
% M1s = tenmat(T,1,[2 3 4]);
M2s = tenmat(T,2);
% M2s = tenmat(T,2,[1 3 4]);
M3s = tenmat(T,3);
% M3s = tenmat(T,3,[1 2 4]);
M4s = tenmat(T,4);
% M4s = tenmat(T,4,[1 2 3]);
size(M4s)

% Same as M using MATLAB reshape
Mm = reshape(T,9,9)

% Same as M1 using MATLAB reshape
M1m = reshape(T,3,27)

% Vectorization
V4 = M4(:)

% Folding
T3 = mat2tens(M3,size(T),3,[1 2 4]);    % tensorlab
T3s = tensor(M3s);                      % Sandia

% Tensors multiplications
% Tensor x vector: a single tensor x vector product produces a matrix
T = randn(3,5,7)
Tm = tens2mat(T,3)
size_tens = size(Tm)
u = randn(7,1)
S = tmprod(T,u',3)
size(S)
S = reshape(u'*Tm,3,5)

% Sandia
Ts = tensor(T);
Ss = ttv(Ts,u,3);

Tm = tens2mat(T,2)
size_tens = size(Tm)
v = randn(5,1)
S = tmprod(T,v',2)
size(S)
S = reshape(v'*Tm,3,1,7)
S = reshape(v'*Tm,3,7)

% Sandia
Ss = ttv(Ts,v,2)
Ss = double(Ss)

% Two consecutive tensor x vector products produce a vector:
s = tmprod(T,{u',v'},[3 2])
S1 = tmprod(T,u',3)
s = tmprod(S1,v',2)

Tm = tens2mat(T,3)
S1 = reshape(u'*Tm,3,5)
s = reshape(v'*S1',3,1)

% Sandia
ss = ttv(Ts,{u,v},[3,2])
ss = double(ss)

% Tensor x matrix:
T = randn(3,5,7)
Tm = tens2mat(T,3)
size_tens = size(Tm)
U = randn(3,7)
S = tmprod(T,U,3)
size(S)
S = reshape((U*Tm)',3,5,3)

% Sandia
Ts = tensor(T)
Ss = ttm(Ts,U,3)


Tm = tens2mat(T,2)
size_tens = size(Tm)
V = randn(4,5)
S = tmprod(T,V,2)
size(S)
S = reshape((V*Tm)',3,4,7)

% Sandia
Ss = ttm(Ts,V,2)

% Two consecutive tensor x matrix products produce a tensor:
S = tmprod(T,{U,V},[3 2])
S1 = tmprod(T,U,3)
S2 = tmprod(S1,V,2)

Tm = tens2mat(T,3)
S1 = reshape((U*Tm)',3,5,3)
S1m = tens2mat(S1,2)
S2 = reshape((V*S1m)',3,4,3)

% Sandia
Ss = ttm(Ts,{U,V},[3 2])

% Kronecker product
A = randn(5,3)
B = randn(3,4)
C = randn(2,6)
AB = kron(A,B)
size(AB)
D = kron(A,B,C)
D = kron(AB,C)
D = kron(kron(A,B),C)
size(D)

a = randn(5,1)
b = randn(3,1)
C = a*b'
c = kron(a,b)

% Contraction via matricization
T = randn(5,2,7)
u = {randn(7,1)}

Tm = tens2mat(T,3);
size_tens = size(Tm);
u_v = u{1};
S = tmprod(T,u_v',3)
size(S)
S = reshape(u_v'*Tm,5,2)

% Tensorlab
Tu_c = contract(T,u,3)
size(Tu_c)

% Sandia
Ts = tensor(T);
Us = tensor(u{1},7);
Ss = ttt(Ts,Us,3,1)
size(Ss)

% TxU
T = randn(5,2,7)
Tm = tens2mat(T,3);
size(Tm);
U = randn(7,12,3)
Um = tens2mat(U,[2 3],1)
size(Um)
size(Tm)
S = reshape((Um*Tm)',5,2,12,3)

% Sandia
Ts = tensor(T)
Us = tensor(U)
Ss = ttt(Ts,Us,3,1)
size(Ss)

% Outer tensor product
A = randn(5,2,7)
B = randn(3,6,4)
ABo = outprod(A,B)
size(ABo)

% Sandia
As = tensor(A);               
Bs = tensor(B);               
ABos = ttt(As,Bs)                
size(ABos)                      


a = randn(7,1)
b = randn(5,1)
c = randn(12,1)
X = outprod(a,b,c)
size(X)

% Sandia
Xs = ktensor({a,b,c})
size(Xs)

% Inner tensor product, Frobenius norm
A = randn(5,2,7)
B = randn(5,2,7)
Ai = inprod(A,A)
frob(A)
frob(A,'squared')
Bi = inprod(B,B)
frob(B)
frob(B,'squared')
ABi = inprod(A,B)
ABi = A(:)'*B(:)

% Sandia
ABis = ttt(tensor(A),tensor(B),[1:3],[1:3])
ABis = ttt(tensor(A),tensor(B),[1:3])

% Kronecker tensor product
T = randn(5,2,7,3)
S = randn(4,3,6,2)
U = kron(T,S)
size(U)
% here we do some additional steps to obtain the correct kronecker product
usize = size(T).*size(S)
U = mat2tens(U,usize,1,[2 3 4])
size(U)

% Khatri-Rao product
A = randn(5,3)
B = randn(4,3)
C = randn(20,3)
AB = [kron(A(:,1),B(:,1)) kron(A(:,2),B(:,2)) kron(A(:,3),B(:,3))]
AB = kr(A,B)
size(AB)
D = kr(kr(A,B),C)
D = kr(A,B,C)
size(D)

% Sandia
ABs = khatrirao(A,B)
size(ABs)
Ds = khatrirao(A,B,C)
size(Ds)


% CPD
% generate pseudorandom factor matrices U(n)
size_tens = [7 8 9]; R = 4;
U = cpd_rnd(size_tens,R);

% By default, cpd_rnd generates each factor U{n} as randn(size_tens(n),R).
% Given a cell array of factor matrices U = {U{1},U{2},...}, its associated
% full tensor T can be computed with

T = cpdgen(U);
T1 = outprod(U{1}(:,1),U{2}(:,1),U{3}(:,1));
T2 = outprod(U{1}(:,2),U{2}(:,2),U{3}(:,2));
T3 = outprod(U{1}(:,3),U{2}(:,3),U{3}(:,3));
T4 = outprod(U{1}(:,4),U{2}(:,4),U{3}(:,4));
T_sum = T1 + T2 + T3 + T4 ;

M = U{1}*kr(U(end:-1:2)).';
size_tens = cellfun('size',U,1);
T_kr = mat2tens(M,size_tens,1);

% Alternating least square (ALS) requiring initialization
U_hat = cpd_als(T,cpd_rnd(size_tens,R))
% Higher order algorithm including ALS and NLS, that does its own internal
% initialization
U_hat = cpd(T,R)
U{1}
U_hat{1}
T_hat = cpdgen(U_hat)
frob(T)
frob(T_hat)

% Residual:
T_diff = T - T_hat
T_diff = cpdres(T,U_hat)

relerr = frob(cpdres(T,U_hat))/frob(T); 
relerr = frobcpdres(T, U_hat)/frob(T);

% Permutation and scaling:
[relerr,P,D,U_hat_ps] = cpderr(U,U_hat);
U{1}
U_hat_ps{1}
U_hat{1}*P*D{1}

% Scaling factors
lambda = diag(D{1}*D{2}*D{3})
U_hat_p = {U_hat{1}*P U_hat{2}*P U_hat{3}*P}

T1_hat_p = outprod(U_hat_p{1}(:,1),U_hat_p{2}(:,1),U_hat_p{3}(:,1));
T2_hat_p = outprod(U_hat_p{1}(:,2),U_hat_p{2}(:,2),U_hat_p{3}(:,2));
T3_hat_p = outprod(U_hat_p{1}(:,3),U_hat_p{2}(:,3),U_hat_p{3}(:,3));
T4_hat_p = outprod(U_hat_p{1}(:,4),U_hat_p{2}(:,4),U_hat_p{3}(:,4));
T_hat_p_sum = T1_hat_p*lambda(1) + T2_hat_p*lambda(2) ...
      + T3_hat_p*lambda(3) + T4_hat_p*lambda(4);
relerr = frob(T - T_hat_p_sum)/frob(T);

% Tensor rank
rankest(T)
R = rankest(T)


% Sandia
Ts = tensor(T)
[Ps1,~,out1] = cp_als(Ts,4)
fit1 = 1 - norm(Ts-full(Ps1))/norm(Ts)
% Sandia + Poblano Toolbox
[Ps2,~,out2] = cp_opt(Ts,4)
fit2 = 1 - norm(Ts-full(Ps2))/norm(Ts)
% Sandia + Poblano Toolbox for sparse tensors
P = (double(Ts)~=0)
[Ps3,~,out3] = cp_wopt(Ts,tensor(P),4)
fit3 = 1 - norm(Ts-full(Ps3))/norm(Ts)

% Rank
for i = 1:6
    [Ps,] = cp_opt(Ts,i);
    fit(i) = 1 - norm(Ts-full(Ps))/norm(Ts);
end
plot(fit,'-r')
ylim([0.2 1.2]);ylabel('fit')
xlim([1 6]);xlabel('rank');grid on
     
%% Tucker decomposition
clear,clc,close all
% We can use the Sandia Tensor Toolbox to create synthetic data (a Tucker
% tensor):
synth_data = create_problem('Type', 'Tucker', 'Size', [5 4 3], 'Num_Factors', [3 3 2]);
Ss = synth_data.Soln
Xs = synth_data.Data

% If we know the number of factors we can use them directly in the
% decomposition:
Ts = tucker_als(Xs,[3 3 2])
fit = 1 - norm(Xs-full(Ts))/norm(Xs)

% Alternatively we can scan all possible combinations to identify an
% optimal decomposition:
Xsize = size(Xs);
fit_mat = zeros(Xsize)
for i = 1:Xsize(1)
    for j = 1:Xsize(2)
        for k = 1:Xsize(3)
            Ts = tucker_als(Xs,[i,j,k]);
            fit_mat(i,j,k) = 1 - norm(Xs-full(Ts))/norm(Xs);
        end
    end
end
fit_mat

[Xaxis,Yaxis] = meshgrid(1:5,1:4)
Tucker_fit = figure
set(gcf,'Unit','Normalized','Position',[0.2 0.2 0.4 0.8]);
subplot(3,1,1)
surf(Xaxis,Yaxis,fit_mat(:,:,1)')
title('R_3 = 1')
xlabel('R_1');ylabel('R_2')
subplot(3,1,2)
surf(Xaxis,Yaxis,fit_mat(:,:,2)')
title('R_3 = 2')
xlabel('R_1');ylabel('R_2')
subplot(3,1,3)
surf(Xaxis,Yaxis,fit_mat(:,:,3)')
title('R_3 = 3')
xlabel('R_1');ylabel('R_2')

% Of course a perfect fit is obtained taking all the possible components:
Ts = tucker_als(Xs,Xsize);
fit = 1 - norm(Xs-full(Ts))/norm(Xs)

% Horizontal (fiber) slices
sv1 = zeros(Xsize(1),1);
sv2 = zeros(Xsize(2),1);
sv3 = zeros(Xsize(3),1);

for i = 1:Xsize(1)
    sv1(i) = norm(Ts.core(i,:,:));
end
for i = 1:Xsize(2)
    sv2(i) = norm(Ts.core(:,i,:));
end
for i = 1:Xsize(3)
    sv3(i) = norm(Ts.core(:,:,i));
end
sv1
sv2
sv3

mlr = mlrankest(double(Xs))

%%
Ts = tucker_als(Xs,[2 2 2]);
fit = 1 - norm(Xs-full(Ts))/norm(Xs)
fTs = full(Ts);
        
%% MLSVD (HOSVD)
clear, clc, close all
S = randn(2,2,3);
U = {rand(4,2), rand(5,2), rand(6,3)};
T = lmlragen(U,S);
size(T)
[Ue,Se,sve] = mlsvd(T)

% Subspaces
subspace(Ue{1}(:,1:2),U{1})
subspace(Ue{2}(:,1:2),U{2})
subspace(Ue{3}(:,1:3),U{3})

% Orthogonal complements
U{1}'*Ue{1}(:, 3:end)
U{2}'*Ue{2}(:, 3:end)
U{3}'*Ue{3}(:, 4:end)

% Multilinear singular values
nrm = zeros(size(Se,1),1);
for i = 1:size(Se,1)
nrm(i) = frob(Se(i,:,:));
end
nrm % mode-1 singular values
sve{1}

% Multilinear rank
for n = 1:3
subplot(1,3,n);
semilogy(sve{n},'.-');
xlim([1 length(sve{n})])
grid on
end

% or
mlR = mlrankest(T)

%% Low multilinear rank approximation
clear, clc, close all
size_tens = [10 20 30];
T = randn(size_tens);
[U,S,sv] = mlsvd(T);
size(S),size(U{1}),size(U{2}),size(U{3})

% Orthogonality of the S core
M = tens2mat(S, 2);
size(M)
M*M'

% As expected all the off-diagonal entries are zero, and the diagonal
% contains the squared multilinear singular values for the second mode.
[sv{2} sv{2}.^2]

% Obtain back the original tensor from the multilinear product
T_hat = tmprod(S,U,1:3);

% Multilinear singular values
nrm = zeros(size(S,1),1);
for i = 1:size(S,1)
nrm(i) = frob(S(i,:,:));
end
nrm % mode-1 singular values
sv{1}

% Multilinear rank
for n = 1:3
subplot(1,3,n);
semilogy(sv{n},'.-');
xlim([1 length(sv{n})])
grid on
end

% or
mlR = mlrankest(T_hat)

%% Truncated MLSVD 
clc, close all
size_core = [4 5 6];
[U_tr,S_tr,sv_tr] = mlsvd(T,size_core);
T_hat_tr = tmprod(S_tr,U_tr,1:3);

% However, S_tr is not orthogonal, so we rerun the mlsvd:
[U_tr,S_tr,sv_tr] = mlsvd(T_hat_tr);

% Multilinear singular values
nrm = zeros(size(S_tr,1),1);
for i = 1:size(S_tr,1)
nrm(i) = frob(S_tr(i,:,:));
end
nrm % mode-1 singular values
sv_tr{1}

% Multilinear rank
for n = 1:3
subplot(1,3,n);
semilogy(sv_tr{n},'.-');
xlim([1 length(sv_tr{n})])
grid on
end

mlR_tr = mlrankest(T_hat_tr)

%% LMLRA
clc, close all
[U_lmlra,S_lmlra] = lmlra(T,size_core);
T_hat_lmlra = tmprod(S_lmlra,U_lmlra,1:3);

% Multilinear singular values
nrm1 = zeros(size(S_lmlra,1),1);
for i = 1:size(S_lmlra,1)
nrm1(i) = frob(S_lmlra(i,:,:));
end
nrm1 % mode-1 singular values
nrm2 = zeros(size(S_lmlra,2),1);
for i = 1:size(S_lmlra,2)
nrm2(i) = frob(S_lmlra(:,i,:));
end
nrm2 % mode-2 singular values
nrm3 = zeros(size(S_lmlra,3),1);
for i = 1:size(S_lmlra,3)
nrm3(i) = frob(S_lmlra(:,:,i));
end
nrm3 % mode-3 singular values
sv_lmlra = {nrm1 nrm2 nrm3};

% Multilinear rank
for n = 1:3
subplot(1,3,n);
semilogy(sv_tr{n},'.-');
xlim([1 length(sv_tr{n})])
grid on
end

mlR_lmlra = mlrankest(T_hat_lmlra)

%% Block Term Decomposition
clear, clc, close all
% Multilinear rank (Lr,Lr,1)
% The command ll1_rnd generates pseudorandom factors for a decomposition in
% multilinear rank-(Lr,Lr,1). For example, consider a tensor of size
% 10�11�12 that can be written as a sum of three terms with multilinear
% ranks (2,2,1), (3,3,1) and (4,4,1), hence L = [2 3 4]:

size_tens = [10 11 12];
L = [2 3 4];
% Default format (btd)
U = ll1_rnd(size_tens, L);
U
U{1}
U{2}
U{3}
U{1}{4}

% CPD format
U = ll1_rnd(size_tens, L,'OutputFormat', 'cpd');
% To expand the factorized representation U of the decomposition in
% multilinear rank-(Lr,Lr,1) terms to a full tensor, ll1gen can be used:
T = ll1gen(U,L)

% or for BTD format
U = ll1_rnd(size_tens, L,'OutputFormat', 'btd');
T = ll1gen(U)

% The easiest way to compute a decomposition in multilinear rank-(Lr,Lr,1)
% terms of a dense, sparse, incomplete or structured tensor ? with
% L=[L1,L2,?,LR] is to use the ll1 method:
[Uhat,output] = ll1(T, L);

% The residual between a tensor T and its decomposition in multilinear
% rank-(Lr,Lr,1) terms given by Uhat can be computed with

res = ll1res(T, Uhat);    % Uhat in BTD format
res = ll1res(T, Uhat, L); % Uhat in CPD or BTD format

% If the tensor is dense, then the result is equivalent to ll1gen(Uhat)-T
% or ll1gen(Uhat,L)-T. The relative error between the tensor and its CPD
% approximation can then be computed as

relerr = frob(ll1res(T,Uhat))/frob(T);   % Uhat in BTD format
relerr = frob(ll1res(T,Uhat,L))/frob(T); % Uhat in CPD or BTD format


% Multilinear rank (Lr,Mr,Nr)
size_tens = [17 19 21];
size_core = {[3 5 7],[6 3 5],[4 3 4]};

% Here we generate the BTD terms, each containing factor matrices U(r,n)
% and core tensors S(r), stored as a cell of cells U.
U = btd_rnd(size_tens,size_core);

% Given a cell array of block terms U = {U{1},U{2},...}, its associated
% full tensor T can be computed with
T = btdgen(U);

% In contrast to the CPD, the LMLRA and the decomposition in multilinear
% rank-(Lr,Lr,1) terms, no specialized high-level algorithm is available
% for a general BTD. This means that the user has to generate his or her
% own initialization and select a specific (family of) algorithm(s) such as
% btd_nls or btd_minf to compute the BTD given the initialization.
% 
% To compute a BTD of a dense, sparse, incomplete or structured tensor T
% given an initialization U0, the command btd_nls(T,U0) can be used. For
% example:
U0 = noisy(U,20);
[Uhat,output] = btd_nls(T,U0);

% The residual between a tensor ? and its BTD defined by the factors �
% (r,n) and the core tensors ?? (r)(stored in Uhat) can be computed with
res = btdgen(Uhat)-T
res = btdres(T,Uhat);

% The relative error between the tensor and its BTD can be computed as
relerr = frob(btdres(T,Uhat))/frob(T); % or
relerr = frobbtdres(T,Uhat)/frob(T); % also works for structured tensors


