tucker.m

多维数据处理：MATLAB源程序用于处理多维数据

                  end;
               end;
               %...using the fast NIPALS
               if Mth(c)==2,
                  if MissingExist,
                     [cthFactor]=fnipals(RedData,Fac(c),reshape(Factors(FIdx0(c):FIdx1(c)),DimX(c),Fac(c)));
                  else
                     [cthFactor]=fnipals(RedData,Fac(c),reshape(Factors(FIdx0(c):FIdx1(c)),DimX(c),Fac(c)));
                  end;
               end;
               %...using simplified continuous Gram-Schmidt orthogonalization
               if Mth(c)==3,
                  if MissingExist,
                     TempMat=RedData*RedData';
                     cthFactor=reshape(Factors(FIdx0(c):FIdx1(c)),DimX(c),Fac(c));
                     for i=1:2,
                        [cthFactor]=gsm(TempMat*cthFactor);
                     end;
                  else
                     TempMat=RedData*RedData';
                     cthFactor=reshape(Factors(FIdx0(c):FIdx1(c)),DimX(c),Fac(c));
                     for i=1:2,
                        [cthFactor]=gsm(TempMat*cthFactor);
                     end;
                  end;
               end;
               %...this is void (no compression for this mode)
               if Mth(c)==4,
               end;
               %...this is void (Keep factors unchanged)
               if Mth(c)==7,
               end;
               %Update the 'Factors' with the current estimates
               if Mth(c)~=4 & Mth(c)~=7
                  Factors(FIdx0(c):FIdx1(c))=cthFactor(:)';
               end;
            end;
         end;
         
         if ~Core_const & Core_uncon==1 & Core_nonneg==0,
            G=calcore(X,DimX,Fac_orig,Factors,[],1,MissingExist);
         elseif Core_nonneg==1,
            g=T3core(X,DimX,Fac,Factors(:),0,1);
            G=reshape(g,Fac(1),prod(Fac(2:N)));
         else
            tmpM2=1;
            for k=1:N;
               if Mth(k)==4,
                  tmpM1=eye(DimX(k));
               else
                  tmpM1=reshape(Factors(FIdx0(k):FIdx1(k)),DimX(k),Fac(k));
               end;
               tmpM2=ckron(tmpM2,tmpM1);
            end
            G=G(:);
            G(fwz)=tmpM2(:,fwz)\X(:);
            enda=size(Fac,2);
            G=reshape(G,Fac(1),prod(Fac(2:enda)));
         end;
         
         SSG=sum(sum(G.^2));
         if MissingExist,
            SSG=SSG-SSMis;
         end;
         if abs(SSG-SSGOld)<Options12*SSGOld,
            Converged1=1;
         end;
         if it1>=Options61,
            itlim1=1;
            Converged1=1;
            Converged2=1;
         end;
         SSGOld=SSG;
         js=0;
         %Save on time count
         if Options101>0 & (etime(clock,t0)>Options101),
            save('temp.mat','Factors','G','DimX','Fac');
            t0=clock;
            js=1;
         end;
         %Save on iteration count
         %if (Options101<0) & (mod(it1,abs(Options101))==0),
         keval=it1/Options51;
         if (Options101<0) & ( abs( keval - floor(keval) ) <=eps),
            save('temp.mat','Factors','G','DimX','Fac');
            js=1;
         end;
         %if mod(it1,Options51)==0 | it1==1  | js==1,
         keval=it1/Options51;
         if (abs( keval - floor(keval) ) <=eps) | it1==1  | js==1, %Matlab 4.2 comp.
            Xm=nmodel(Factors,G,DimX,Fac_orig);
            
            if MissingExist,
               X(IdxIsNans)=Xm(IdxIsNans);
               SSMis=sum(sum( Xm(IdxIsNans).^2 ));
               if abs(SSMis-SSMisOld)<Options11*SSMisOld,
                  Converged2=1;
               end;
               SSMisOld=SSMis;
            else
               Converged2=1;
            end;
            ExplX=100*(1-sum(sum((X-Xm).^2))/SSX);
            pout=pout+1;
            if pout>pmore,
               if prlvl > 0,
                  fprintf('%s\n',str1);
                  fprintf('%s\n',str2);
               end;
               pout=0;
            end;
            
            if prlvl>0,
               if MissingExist,
                  fprintf(' %6i       %14.3f     %14.3f         %8.4f',it1,SSG,SSMis,ExplX);
               else
                  fprintf(' %6i        %14.3f      %8.4f',it1,SSG,ExplX);
               end;
               if js,
                  fprintf(' - saved to ''temp.mat'' \n')
               else
                  fprintf('\n')
               end;
            end;
         end;
      end; %Inner loop
   end; %Outer loop
   if prlvl>0,
      if itlim1==0,
         fprintf('   Stopped. Convergence criteria reached.\n');
      else
         fprintf('   Stopped. Iteration limits reached in model and expectation loops.\n');
      end;      
      if MissingExist,
         fprintf(' %6i       %14.3f     %14.3f         %8.4f',it1,SSG,SSMis,ExplX);
      else
         fprintf(' %6i        %14.3f      %8.4f',it1,SSG,ExplX);
      end;
   end;
   if Options101~=0,
      save('temp.mat','Factors','G','DimX','Fac');
      if prlvl>0,
         fprintf(' - saved to ''temp.mat'' \n')
      end;
   else
      if prlvl>0,
         fprintf('\n')
      end;
   end;
   
elseif MethodO==2, %Must use slower but more general schemes
   
   SSGOld=0;
   OldFactors=0*Factors;
   Converged2=0;
   it1=0;
   t0=clock;
   itlim1=0;
   while ~Converged2, 
      Converged1=0;
      while ~Converged1,
         it1=it1+1;   
         %Iterate over the modes
         for c=1:N,
            faclist1=[N:-1:1 N:-1:1];
            faclist=faclist1(N-c+2:N-c+2+(N-2));
            tmpM2=1;
            tmpM2Pinv=1;
            for k=1:N-1;
               if Mth(faclist(k))==4,
                  tmpM1=eye(DimX(faclist(k)));
               else
                  tmpM1=reshape(Factors(FIdx0(faclist(k)):FIdx1(faclist(k))),DimX(faclist(k)),Fac(faclist(k)));
               end;
               if CalcOrdinar(c)==1,
                  tmpM2=ckron(tmpM2,tmpM1');
               end
            end
            %Estimate Factors for the cth way
            %...this is void (no compression for this mode)
            if any(Mth(c)==[1 2 3]),
               tmpM4=G*tmpM2;
               if MissingExist,
                  %cthFactor=(X*tmpM4')/((tmpM4*X'*X*tmpM4')^(1/2));
                  [U S V]=svd(tmpM4*X'*X*tmpM4',0);
                  mS=min(size(S));
                  Sm=S;
                  Sm(1:mS,1:mS)=diag(1./sqrt(diag(S(1:mS,1:mS))));
                  pinvm=U*Sm*V';
                  cthFactor=X*tmpM4'*pinvm;
               else
                  %cthFactor=(X*tmpM4')/((tmpM4*X'*X*tmpM4')^(1/2));
                  [U S V]=svd(tmpM4*X'*X*tmpM4',0);
                  mS=min(size(S));
                  Sm=S;
                  Sm(1:mS,1:mS)=diag(1./sqrt(diag(S(1:mS,1:mS))));
                  pinvm=U*Sm*V';
                  cthFactor=X*tmpM4'*pinvm;
               end;
            end;
            if Mth(c)==4,
            end;
            if Mth(c)==5, %Nonneg
               tmpM3=reshape(Factors(FIdx0(c):FIdx1(c)),DimX(c),Fac(c));
               if it1==1,
                  %tmpM3(find(tmpM3<0))=0;
                  i__=find(tmpM3<0);
                  if ~isempty(i__),
                     tmpM3(i__)=zeros(1,length(i__));
                  end;
               end;
               if CalcOrdinar(c) == 1,
                  tmpM4=G*tmpM2;
                  if dbg, ss1=sum(sum( (X-tmpM3*tmpM4).^2 ));end;
                  [cthFactor,SS]=nonneg(X,tmpM4,tmpM3);
                  if dbg, ss2=sum(sum( (X-cthFactor*tmpM4).^2 ));
                     fprintf('Nonng report (Ordi) %15.8d  %15.8d\n',ss1,ss2); end;
               end;
            end;
            if Mth(c)==6, %Uncon
               cthFactor=reshape(Factors(FIdx0(c):FIdx1(c)),DimX(c),Fac(c));
               if CalcOrdinar(c) == 1,
                  M_=G*tmpM2;
                  if dbg, ss1=sum(sum( (X-cthFactor*G*tmpM2).^2 ));end;
                  cthFactor=X/M_;
                  if dbg, ss2=sum(sum( (X-cthFactor*G*tmpM2).^2 ));
                     fprintf('Uncon report (Ordi) %15.8d  %15.8d\n',ss1,ss2);end;
               end;
            end;
            if Mth(c)==7, %Not updated
            end;
            if Mth(c)==8, %Unimodality
               tmpM3=reshape(Factors(FIdx0(c):FIdx1(c)),DimX(c),Fac(c));
               if it1==1,
                  %tmpM3(find(tmpM3<0))=0;
                  i__=find(tmpM3<0);
                  if ~isempty(i__),
                     tmpM3(i__)=zeros(1,length(i__));
                  end;
               end;
               if CalcOrdinar(c) == 1,
                  tmpM4=G*tmpM2;
                  if dbg, ss1=sum(sum( (X-tmpM3*tmpM4).^2 ));end;
                  cthFactor = unimodal(tmpM4',X');
                  if dbg, ss2=sum(sum( (X-cthFactor*tmpM4).^2 ));fprintf('Unimod report (Ordi) %15.8d  %15.8d\n',ss1,ss2); end;
               end;
            end;
            
            
            %Update 'Factors' with the current estimates
            if Mth(c)~=4 & Mth(c)~=7
               cn=sum(cthFactor.^2);
               for i=1:Fac(c);
                  if it1<=10 & cn(i)<eps,
                     cthFactor(:,i)=rand(size(cthFactor,1),1);
                  end;
               end;
               if Core_const,
                  if c>1
                     cthFactor=normit(cthFactor);
                  end;
               else
                  cthFactor=normit(cthFactor);
               end;
               Factors(FIdx0(c):FIdx1(c))=cthFactor(:)';
            end;
            
            %Estimate the new core after each SUBiteration
            CoreupdatedInner=0;
            if ~Core_const
               if UpdateCore(c)==1,
                  tmpM3=reshape(Factors(FIdx0(c):FIdx1(c)),DimX(c),Fac(c));
                  if CalcOrdinar(c) == 1,
                     if dbg, ss1=sum(sum( (X-tmpM3*G*tmpM2).^2 ));end;
                     G=tmpM3\(X/tmpM2);
                     if dbg, ss2=sum(sum( (X-tmpM3*G*tmpM2).^2 ));
                        fprintf('Core report  (Ordi) %15.8d  %15.8d\n',ss1,ss2);end;
                  end;
                  CoreupdatedInner=1;
               end;
               if UpdateCore(c)==2,
                  G=zeros(Fac(c),prod(Fac(nsetdiff(1:N,c))));
                  tmpM2=1;
                  for k=[(c-1):-1:1 N:-1:c];
                     if Mth(k)==4,
                        tmpM1=eye(DimX(k));
                     else
                        tmpM1=reshape(Factors(FIdx0(k):FIdx1(k)),DimX(k),Fac(k));
                     end;
                     tmpM2=ckron(tmpM2,tmpM1);
                  end
                  G=G(:);
                  fwz=find(G_cons==1);
                  G(fwz)=tmpM2(:,fwz)\X(:);
                  G=reshape(G,Fac(c),prod(Fac(nsetdiff(1:N,c))));
                  CoreupdatedInner=1;
               end;
            end;
            
            %Reshape to the next unfolding
            if c~=N,
               newix=DimX(c+1);
               newjx=prod(DimX)/DimX(c+1);
               newig=Fac(c+1);
               newjg=prod(Fac)/Fac(c+1);
            else
               newix=DimX(1);
               newjx=prod(DimX)/DimX(1);
               newig=Fac(1);
               newjg=prod(Fac)/Fac(1);
            end;
            X=reshape(X',newix,newjx);
            G=reshape(G',newig,newjg);
            if length(G_cons(:))>1,
               G_cons=reshape(G_cons',newig,newjg);
            end;
         end;
         
         %Estimate the new core after each MAIN iteration
         if CoreupdatedInner==0,
            if ~Core_const & Core_cmplex==0,
               if Core_nonneg==1,
                  g=t3core(X,DimX,Fac,Factors,0,1);
                  G=reshape(g,Fac(1),prod(Fac(2:N)));
               elseif Core_nonneg==0,
                  G=calcore(X,DimX,Fac,Factors,[],0,MissingExist);
               end
            end;
         end;
         CoreupdatedInner=0;
         
         if ~Core_const
            SSG=sum(sum(G.^2));
            if MissingExist,
               SSG=SSG-SSMis;
            end;
            if abs(SSG-SSGOld)<Options12*SSGOld,
               Converged1=1;
            end;
            if it1>=Options61,
               Converged1=1;
               Converged2=1;
               itlim1=1;
            end;
            SSGOld=SSG;
         else
            SSF=sum(sum((Factors-OldFactors).^2));
            if SSF<Options12*sum(sum(Factors.^2));
               Converged1=1;
            end;
            if it1>=Options61,
               Converged1=1;
               Converged2=1;
               itlim1=1;
            end;
            OldFactors=Factors;
            SSG=SSF;
         end;
         
         js=0;
         if Options101>0 & (etime(clock,t0)>Options101),
            save('temp.mat','Factors','G','DimX','Fac');
            t0=clock;
            js=1;
         end;
         %if Options101<0 & mod(it1,abs(Options101))==0,
         keval=it1/Options101;
         if (Options101<0) & (abs( keval - floor(keval) ) <=eps), %Matlab 4.2 comp.
            save('temp.mat','Factors','G','DimX','Fac');
            js=1;
         end;
         %if mod(it1,Options51)==0 | it1==1 | js==1,
         keval=it1/Options51;
         if (abs( keval - floor(keval) ) <=eps) | it1==1  | js==1, %Matlab 4.2 comp.
            Xm=nmodel(Factors,G,DimX,Fac_orig);
            if MissingExist,
               X(IdxIsNans)=Xm(IdxIsNans);
               SSMis=sum(sum( Xm(IdxIsNans).^2 ));
               if abs(SSMis-SSMisOld)<Options11*SSMisOld,
                  Converged2=1;
               end;
               SSMisOld=SSMis;
            else
               Converged2=1;
            end;
            ExplX=100*(1 - sum(sum((X-Xm).^2))/SSX );
            pout=pout+1;
            if pout>pmore,
               if prlvl>0,
                  fprintf('%s\n',str1);
                  fprintf('%s\n',str2);
               end;
               pout=0;
            end;              
            if prlvl>0,
               if MissingExist,
                  fprintf(' %6i       %14.3f     %14.3f         %8.4f',it1,SSG,SSMis,ExplX);
               else
                  fprintf(' %6i        %14.3f      %8.4f',it1,SSG,ExplX);
               end;
               if js,
                  fprintf(' - saved to ''temp.mat'' \n')
               else
                  fprintf('\n')
               end;
            end;
         end;
      end; %Inner loop
   end; %Outer loop
   if prlvl>0,
      if itlim1==0,
         fprintf('   Stopped. Convergence criteria reached.\n');
      else
         fprintf('   Stopped. Iteration limits reached in model and expectation loops.\n');
      end;      
      if MissingExist,
         fprintf(' %6i       %14.3f     %14.3f         %8.4f',it1,SSG,SSMis,ExplX);
      else
         fprintf(' %6i        %14.3f      %8.4f',it1,SSG,ExplX);
      end;
   end;
   if Options101~=0,
      save('temp.mat','Factors','G','DimX','Fac');
      if prlvl>0,
         fprintf(' - saved to ''temp.mat'' \n')
      end;
   else
      if prlvl>0,
         fprintf('\n')
      end;
   end;
end; %Outer loop

if MissingExist,
   Xf=Xm(IdxIsNans);
end;

Factors=Factors';
format