function [A,B,C,FIT,IT]=comfac(X,Fac,Options,DoComp,CompIt,Init);

%COMFAC Algorithm for fitting the complex-valued PARAFAC model
% 
% See e.g. Rasmus Bro, N. D. Sidiropoulos, and G. B. Giannakis. A Fast 
%          Least Squares Algorithm for Separating Trilinear Mixtures. 
%          Proc. ICA99 ?Int. Workshop on Independent Component Analysis 
%          and Blind Signal Separation. Jan. 11?5, Aussois, France:289-294, 1999.
%
%	        N. D. Sidiropoulos, G. B. Giannakis, and Rasmus Bro.. Blind PARAFAC 
%          Receivers for DS-CDMA Systems. IEEE Transactions on Signal Processing March, 2000.
%
% The algorithm works by first compressing the data using a Tucker3 models. Subsequently the 
% PARAFAC model is fitted to the compressed array, either initialized with DTLD (~ESPRIT) or 
% with PARAFAC-ALS. The solution is de-compressed to the original space, and a few safe-guard
% PARAFAC-ALS steps are performed. Further optimization of this algorithm is possible, e.g.
% if the data are very large, if the profiles are very correlated or if the noise is huge. This
% has not been pursued here, in order to have a generally applicable algorithm
% 
% INPUT
% X       : I x J x K data three-way array 
% Fac     : Number of factors in PARAFAC model
%
% OPTIONAL INPUTS
% Options : A 1x2 vector
% (1)     : Mode to compress to dimension two in
%           DTLD if smallest mode-dimension is more
%           than two. Default is the smallest dimension
% (2)     : Number of extra components in Tucker3
%           compression model. Default if no given 
%           is one
% 
% ADVANCED OPTIONS
% DoComp  : 0      => No compression
%           1      => Compression
% CompIt  : CompIt => Max number of iterations in Compression
% Init    : 0      => GRAM/DTLD
%           1      => Ten x PARAFAC with 10 iterations started from RandOrth
%
% 
% I/O
% [A,B,C,FIT,IT]=comfac(X,Fac,Options,DoComp,CompIt,Init);
% 
% Or short : [A,B,C]=comfac(X,Fac);
%
%
% Copyright, 1998 - 
% This M-file and the code in it belongs to the holder of the
% copyrights and is made public under the following constraints:
% It must not be changed or modified and code cannot be added.
% The file must be regarded as read-only. In case of doubt, 
% contact the holder of the copyrights.
%
% Copyright 1998
% Rasmus Bro
% KVL, Denmark, rb@kvl.dk
%      & 
% Nikos Sidiropoulos
% Univ. Minnesota, nikos@ece.umn.edu
% 

% 31/3/00 RB fixed bug when dimensions not suitable for DTLD

crit = 1e-5; % Criterion used for fitting with Alternating LS algorithm

if length(size(X))~=3
   error(' The data must be held in a three-way array')
end
   
DimX = size(X);
X = reshape(X,DimX(1),prod(DimX(2:3)));

dbg=1;

DeafultExtraCompInTucker=2;
if exist('DoComp')~=1|isempty(DoComp)
   DoComp = 1;
end
if exist('CompIt')~=1|isempty(CompIt)
   CompIt = 3;
end
if exist('Init')~=1|isempty(Init)
   Init   = 1;
end
if length(find(DimX>Fac))<2 % GRAM/ESPRIT cannot be used
   ReplaceTLDwithRand = 1;
else
   ReplaceTLDwithRand = 0;
end   

if dbg
%   home
   disp(' ')
   disp([' Fitting ',num2str(Fac),'-comp. PARAFAC model using COMFAC ...'])
   disp([' Array size: ',num2str(DimX(1)),' x ',num2str(DimX(2)),' x ',num2str(DimX(3))])
end

if nargin<3|isempty(Options)
  Options=[0 DeafultExtraCompInTucker];
end
if Options(1)==0
  [out,CompToTwo] = min(DimX);
else
  CompToTwo = Options(1);
end
if length(Options)<2
  Options(2)=DeafultExtraCompInTucker; % Number of additional components in compression as compared to model
end
W = Options(2)+Fac;
W = [W W W];
for i=1:3
  if W(i)>DimX(i)
    W(i)=DimX(i);
  end
end

% Display things
if DoComp & dbg
  disp([' Array compressed to size ',num2str(W(1)),' x ',num2str(W(2)),' x ',num2str(W(3))])
end

% Compression

if DoComp
   [Ut,Vt,Zt,Gt,fit]=tucker3(X,DimX,W,CompIt);
else
   Gt = X;
   W = DimX;  
   Ut=speye(DimX(1));
   Vt=speye(DimX(2));
   Zt=speye(DimX(3));
end

% Initialize PARAFAC using DTLD
if Init == 0
   if ~ReplaceTLDwithRand
      disp(' Initializing using direct trilinear decomposition')
      [Adtld,Bdtld,Cdtld]=cdtld(Gt,W,Fac,CompToTwo);
   else
      disp(' Initializing using random values because tld cannot be used due to the sizes')
      [Adtld,Bdtld,Cdtld,fit]=cparafac(Gt,W,Fac,crit,0,0,0,10);
   end
   
elseif Init == 1
   disp(' Initializing using best of 11 initial small runs')
   if ~ReplaceTLDwithRand
      [Adtld,Bdtld,Cdtld]=cdtld(Gt,W,Fac,CompToTwo);
   else
      [Adtld,Bdtld,Cdtld,fit]=cparafac(Gt,W,Fac,crit,0,0,0,10);
   end
  fitout=sum(sum(abs(  Gt-Adtld*ppp(Bdtld,Cdtld).').^2));
  for rep=1:10
     [adtld,bdtld,cdtldvar,fit]=cparafac(Gt,W,Fac,crit,0,0,0,10);
     if fit<fitout
       Adtld=adtld;Bdtld=bdtld;Cdtld=cdtldvar;fitout=fit;
     end
  end
end

% Fit PARAFAC model in compressed space
if DoComp
  disp(' Fitting PARAFAC in compressed space')
  [Acomp,Bcomp,Ccomp,fit2,it2]=cparafac(Gt,W,Fac,crit,Adtld,Bdtld,Cdtld);
  % Decompress to original domain
  disp(' Transforming interim solution to original space')
  Ainit = Ut*Acomp;
  Binit = Vt*Bcomp;
  Cinit = Zt*Ccomp;
else
  Ainit = Adtld;
  Binit = Bdtld;
  Cinit = Cdtld;
  it2=0;
  fit2=NaN;
end

% Fit PARAFAC model in raw space
disp(' Fitting PARAFAC in original space')
[A,B,C,FIT,IT]=cparafac(X,DimX,Fac,crit,Ainit,Binit,Cinit);
disp(' Algorithm converged')


function [A,B,C,fit,it]=cparafac(X,DimX,Fac,crit,A,B,C,maxit,DoLineSearch);

% Complex PARAFAC-ALS
% Fits the PARAFAC model Xk = A*Dk*B.' + E
% where Dk is a diagonal matrix holding the k'th
% row of C.
%
% Uses on-the-fly projection-compression to speed up 
% the computations. This requires that the first mode 
% is the largest to be effective
% 
% INPUT
% X       : Data
% DimX    : Dimension of X
% Fac     : Number of factors
% OPTIONAL INPUT
% crit    : Convergence criterion (default 1e-6)
% A,B,C   : Initial parameter values
%
% I/O
% [A,B,C,fit,it]=parafac(X,DimX,Fac,crit,A,B,C);
%
% Copyright 1998
% Rasmus Bro
% KVL, Denmark, rb@kvl.dk


% Initialization
if nargin<8
  maxit   = 2500;      % Maximal number of iterations
end
showfit = pi;         % Show fit every 'showfit'th iteration (set to pi to avoid)

if nargin<4
  crit=1e-6;
end

if crit==0
  crit=1e-6;
end

I = DimX(1);
J = DimX(2);
K = DimX(3);

InitWithRandom=0;
if nargin<7
   InitWithRandom=1;
end
if nargin>6 & size(A,1)~=I
  InitWithRandom=1;
end

if InitWithRandom

  if I<Fac
    A = rand(I,Fac);
  else
    A = orth(rand(I,Fac));
  end
  if J<Fac
    B = rand(J,Fac);
  else
    B = orth(rand(J,Fac));
  end
  if K<Fac
    C = rand(K,Fac);
  else
    C = orth(rand(K,Fac));
  end
end

SumSqX = sum(sum(abs(X).^2));
fit    = sum(sum(abs(X-A*ppp(B,C).')));
fit0   = fit;
fitold = 2*fit;
it     = 0;
Delta  = 5;

while abs((fit-fitold)/fitold)>crit & it<maxit & fit>10*eps
   it=it+1;
   fitold=fit;

   % Do line-search
   if rem(it+2,2)==-1
      [A,B,C,Delta]=linesrch(X,DimX,A,B,C,Ao,Bo,Co,Delta);
   end

   Ao=A;Bo=B;Co=C;
   % Update A
   Xbc=0;

   for k=1:K
     Xbc = Xbc + X(:,(k-1)*J+1:k*J)*conj(B*diag(C(k,:)));
   end
   A = Xbc*inv((B'*B).*(C'*C)).';

   % Project X down on orth(A) - saves time if first mode is large
   [Qa,Ra]=qr(A,0);
   x=Qa'*X;

   % Update B
   Xac=0;
   for k=1:K
     Xac = Xac + x(:,(k-1)*J+1:k*J).'*conj(Ra*diag(C(k,:)));
   end
   B = Xac*inv((Ra'*Ra).*(C'*C)).';
   
   % Update C
   ab=inv((Ra'*Ra).*(B'*B));
   for k=1:K
     C(k,:) = (ab*diag(Ra'* x(:,(k-1)*J+1:k*J)*conj(B))).';
   end

   % Calculating fit. Using orthogonalization instead
   %fit=0;for k=1:K,residual=X(:,(k-1)*J+1:k*J)-A*diag(C(k,:))*B.';fit=fit+sum(sum((abs(residual).^2)));end
   [Qb,Rb]=qr(B,0);
   [Z,Rc]=qr(C,0);
   fit=SumSqX-sum(sum(abs(Ra*ppp(Rb,Rc).').^2));

   if rem(it,showfit)==0
     fprintf(' %12.10f       %g        %3.4f \n',fit,it,100*(1-fit/fit0));
   end
end

%fprintf(' %12.10f       %g        %3.4f \n',fit,it,100*(1-fit/fit0));


function [A,B,C]=cgram(X1,X2,F);

% cGRAM - Complex Generalized Rank Annihilation Method
% Fits the PARAFAC model directly for the case of a 
% three-way array with only two frontal slabs.
% For noise-free trilinear data the algorithm is exact.
% 
% INPUTS:
% X1    : I x J matrix of data from observation one
% X2    : I x J matrix of data from observation two
% Fac   : Number of factors
% 
% OUTPUTS:
% A     : Components in the row mode (I x F)
% B     : Components in the column mode (J x F)
% C     : Weights for each slab; C(1,:) are the component 
%         weights for first slab such that the approximation
%         of X1 is equivalent to X1 = A*diag(C(1,:))*B.'
%
% Copyright 1998
% Rasmus Bro, KVL, DK
% rb@kvl.dk


  DontShowOutput = 1;

  [U,s,V]=svd(X1+X2);
  U=U(:,1:F);
  V=V(:,1:F);

%  S2=S1*b
%  b*k = k*l =>
%  S1*b*k = S1*k*l =>
%  S2*k = S1*k*l =>
%  inv(S1*k)*S2*k = l = diagonal

  S1=U'*X1*V;
  S2=U'*X2*V;
  [v,d]=eig(S1\S2);

  ddd=d;
  d=diag(d);out=abs(d)>eps;v=v(:,out);d=d(out); %only significant terms
  [dd,out]=sort(abs(d));out=flipud(out);
  d=d(out);d=diag(d);v(:,out);v=v/norm(v);% sort them

  A = U*S1*v;
  B=V/v';
  B=(B.')';
  C=(pinv(ppp(A,B))*[X1(:) X2(:)]).';

  if ~DontShowOutput
    fit = sum(sum(abs([X1 X2] - [A*diag(C(1,:))*B.' A*diag(C(2,:))*B.']).^2));
    disp([' GRAM fitted data with a sum-squared error of ',num2str(fit)])
 end
 
 
 
 function [A,B,C,fit]=cdtld(X,DimX,F,SmallMode);

% DIRECT TRILINEAR DECOMPOSITION
% calculate the parameters of the three-
% way PARAFAC model directly. The model
% is not the least-squares but will be close
% to for precise data with little model-error
%
% This implementation works with an optimal
% compression using least-squares Tucker3 fitting
% to generate two pseudo-observation matrices that
% maximally span the variation of all samples. per
% default the mode of smallest dimension is compressed
% to two samples, while the remaining modes are 
% compressed to dimension F.
% 
% For large arrays it is fastest to have the smallest
% dimension in the first mode
%
% INPUT
% [A,B,C]=dtld(X,DimX,F);
% X is the I x J x K array unfolded to an I x JK matrix
% DimX = [I J K]
% F is the number of factors to fit
% An optional parameter may be given to enforce which
% mode is to be compressed to dimension two
%
% Copyright 1998
% Rasmus Bro, KVL
% rb@kvl.dk


DontShowOutput = 1;

%rearrange X so smallest dimension is in first mode
if nargin<4
  [a,SmallMode] = min(DimX);
  X = nshape(X,DimX,SmallMode);
  DimX = DimX([SmallMode 1:SmallMode-1 SmallMode+1:3]);
  Fac   = [2 F F];
else
  X = nshape(X,DimX,SmallMode);
  DimX = DimX([SmallMode 1:SmallMode-1 SmallMode+1:3]);
  Fac   = [2 F F];
end

if DimX(1) < 2
  error(' The smallest dimension must be > 1')
end
if any(DimX(2:3)-Fac(2:3)<0)
  error(' This algorithm requires that two modes are of dimension not less the number of components')
end

% Compress data into a 2 x F x F array. Only 10 iterations are used since exact SL fit is insignificant; only obtaining good truncated bases is important
[At,Bt,Ct,G]=tucker3(X,DimX,Fac,10);

% Fit GRAM to compressed data
[Bg,Cg,Ag]=cgram(reshape(G(1,:),F,F),reshape(G(2,:),F,F),F);

% De-compress data and find A
BB = Bt*Bg;