%Ekta Gujral,Tianxiong Yang,Ravdeep Pasricha and Evangelos E. Papalexakis
%Department of Computer Science and Engineering, University of California Riverside 
%"OCTen:Online Compression-based Tensor Decomposition", Submitted in PKDD 2018
function [A,B,C,U,V,summaries]=OCTen_initilization(X,ops)
    if ~isfield(ops, 'p'), ops.p = 10; end %number of compressed cube 
    if ~isfield(ops, 'Q'), ops.Q = 10; end  %size of compressed dimension 
    if ~isfield(ops, 'shared'), ops.shared = round(ops.Q/2); end % number of shared projector columns
    if ~isfield(ops, 'R'), ops.R = 2; end % number of factors
    sz=size(X);
    I=sz(1);J=sz(2);K=sz(3);
    U=randn(I,ops.Q,ops.p); V=randn(J,ops.Q,ops.p);W=randn(K,ops.Q,ops.p);
     for pp=2:ops.p
        U(:,1:ops.shared,pp)=U(:,1:ops.shared,1); 
        V(:,1:ops.shared,pp)=V(:,1:ops.shared,1); 
        W(:,1:ops.shared,pp)=W(:,1:ops.shared,1); 
    end
      % Takes up more memory but faster:
      As=zeros(ops.Q,ops.R,ops.p); % holds all P reduced mode loadings for the A mode
      Bs=zeros(ops.Q,ops.R,ops.p); % likewise for B mode 
      Cs=zeros(ops.Q,ops.R,ops.p); % likewise for C mode 
 parfor pp=1:ops.p 
    %     fprintf('Now computing compressed cube number %d\n',pp);
        G1 = reshape(double(X),I,J*K);
        G2=reshape(U(:,:,pp).'*G1,ops.Q*J,K).'; 
        G3=reshape(W(:,:,pp).'*G2,ops.Q*ops.Q,J).'; 
        G4=reshape(V(:,:,pp).'*G3,ops.Q*ops.Q,ops.Q).'; 
        summaries{pp}=reshape(G4,ops.Q,ops.Q,ops.Q); 
        Xs= summaries{pp};
        [As(:,:,pp),Bs(:,:,pp),Cs(:,:,pp),~]=comfac(Xs, ops.R,2);
 end
    clear G1; clear G2;clear G3;clear G4;clear X; % not needed anymore
 AsTop=As(1:ops.shared,:,:); 
 BsTop=Bs(1:ops.shared,:,:); 
 CsTop=Cs(1:ops.shared,:,:); 
  for pp=1:ops.p
    for f=1:ops.R 
        [~,mloc]=max(abs(AsTop(:,f,pp)));
        AsTop(:,f,pp)=AsTop(:,f,pp)/AsTop(mloc,f,pp);
        [~,mloc]=max(abs(BsTop(:,f,pp)));
        BsTop(:,f,pp)=BsTop(:,f,pp)/BsTop(mloc,f,pp);
        [~,mloc]=max(abs(CsTop(:,f,pp)));
       CsTop(:,f,pp)=CsTop(:,f,pp)/CsTop(mloc,f,pp);
        
    end
  end
    
  for pp=2:ops.p 
    % permute columns of each AtildeTop(:,:,p) to match AtildeTop(:,:,1)
  %  fprintf('Now computing permutation for partial result number %d \n',p);
    [AsTop(:,:,pp),PermpA] = perm2match(AsTop(:,:,pp),AsTop(:,:,1));
    [BsTop(:,:,pp),PermpB] = perm2match(BsTop(:,:,pp),BsTop(:,:,1));
    [CsTop(:,:,pp),PermpC] = perm2match(CsTop(:,:,pp),CsTop(:,:,1));
    
    
    % Now permute columns of Atilde(:,:,p) to match the perm in Atilde(:,:,1)
    As(:,:,pp) = As(:,:,pp)*PermpA;
    Bs(:,:,pp) = Bs(:,:,pp)*PermpB;
    Cs(:,:,pp) = Cs(:,:,pp)*PermpC;
    
        
    % NOW equalize matched column scales by looking at Atilde(1:S,:,:) which are common: we pick the first row for this here: 
    LambdapA = diag(As(1,:,1)./As(1,:,pp));
    LambdapB = diag(Bs(1,:,1)./Bs(1,:,pp));
    LambdapC = diag(Cs(1,:,1)./Cs(1,:,pp));
    
    % apply to the full Atilde(:,:,p):
    As(:,:,pp) = As(:,:,pp)*LambdapA;
    Bs(:,:,pp) = Bs(:,:,pp)*LambdapB;
    Cs(:,:,pp) = Cs(:,:,pp)*LambdapC;
 
  end  
    clear AsTop; clear BsTop;clear CsTop; clear LambdapA; clear LambdapB; clear LambdapC; % no need now

   Q=ops.Q; 
 allAs = zeros(Q+(Q-ops.shared)*(ops.p-1),ops.R);
allUTA=zeros(Q+(Q-ops.shared)*(ops.p-1),I);
allAs(1:Q,:)=As(:,:,1);
allUTA(1:Q,:)=U(:,:,1)';
for pp=2:ops.p 
    allAs(Q+(Q-ops.shared)*(pp -2)+1:Q+(Q-ops.shared)*(pp-1),:)=As(ops.shared+1:Q,:,pp);
    allUTA(Q+(Q-ops.shared)*(pp -2)+1:Q+(Q-ops.shared)*(pp-1),:)=U(:,ops.shared+1:Q,pp)';
end


allBs = zeros(Q+(Q-ops.shared)*(ops.p-1),ops.R);
allUTB=zeros(Q+(Q-ops.shared)*(ops.p-1),J);
allBs(1:Q,:)=Bs(:,:,1);
allUTB(1:Q,:)=V(:,:,1)';
for pp=2:ops.p 
    allBs(Q+(Q-ops.shared)*(pp -2)+1:Q+(Q-ops.shared)*(pp-1),:)=Bs(ops.shared+1:Q,:,pp);
    allUTB(Q+(Q-ops.shared)*(pp -2)+1:Q+(Q-ops.shared)*(pp-1),:)=V(:,ops.shared+1:Q,pp)';
end

allCs = zeros(Q+(Q-ops.shared)*(ops.p-1),ops.R);
allUTC=zeros(Q+(Q-ops.shared)*(ops.p-1),K);
% allUTC=zeros(Q+(Q-ops.shared)*(ops.p-1),K);
allCs(1:Q,:)=Cs(:,:,1);
allUTC(1:Q,:)=W(:,:,1)';
for pp=2:ops.p 
    allCs(Q+(Q-ops.shared)*(pp-2)+1:Q+(Q-ops.shared)*(pp-1),:)=Cs(ops.shared+1:Q,:,pp);
    allUTC(Q+(Q-ops.shared)*(pp-2)+1:Q+(Q-ops.shared)*(pp-1),:)=W(:,ops.shared+1:Q,pp)';
end
    
  %solve for the full size decomposition using least squares
    A=pinv(allUTA)*allAs; % solve the master equation; equiv. to A=pinv(allUT)*allAtilde;
    B=pinv(allUTB)*allBs;
    C=pinv(allUTC)*allCs;
     
end
