calcore.m

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

function [G]=calcore(X,DimX,Fac,Factors,Options,O,MissingExist);
%
% $ Version 1.01 $ Date 17. Sep. 1998 $ Not compiled $
%
% ----------------------------------------------------
%           Calculate the Tucker core
% ----------------------------------------------------
%
%	Copyright
%	Claus A. Andersson 1995-1997
%	Chemometrics Group, Food Technology
%	Department of Food and Dairy Science
%	Royal Veterinary and Agricultutal University
%	Rolighedsvej 30, T254
%	DK-1958 Frederiksberg
%	Denmark
%
%	Phone 	[phone]
%	Fax	   [phone]
%	E-mail	[e-mail]
%
%	
% [G]=calcore(X,DimX,Fac,Factors,Options);
% [G]=calcore(X,DimX,Fac,Factors);
%
% This algorithm applies to the general N-way case, so
% the unfolded X can have any number of dimensions. The principles of
% 'projections' and 'systematic unfolding methodology (SUM)' are used
% in this algorithm so orthogonality is required.
% This algorithm can handle missing values in X and
% also allows for TUCKER2 models using the '-1' flag in 'Fac'.
% The variable 'Factors' must contain the stringed-out factors.

format compact
format long

% Initialize system variables
if length(Fac)==1,
    Fac=Fac*ones(size(DimX));
end;

Fac_orig=Fac;
i=find(Fac==-1);
Fac(i)=zeros(1,length(i));
N=size(Fac,2);
FIdx0=zeros(1,N);
FIdx1=zeros(1,N);
if ~exist('MissingExist')
    if sum(isnan(X(:)))>0,
        MissingExist=1;
    else
        MissingExist=0;
    end;
end;
FIdx0=cumsum([1 DimX(1:N-1).*Fac(1:N-1)]);
FIdx1=cumsum([DimX.*Fac]);
if ~exist('O') | isempty(O),
    O=1;
end;


if O, %means orthogonality
    CurDimX=DimX;
    RedData=X;
    for c=1:N,
        
        if Fac_orig(c)==-1,
            kthFactor=eye(DimX(c));
            CurDimX(c)=DimX(c);
        else
            kthFactor=reshape(Factors(FIdx0(c):FIdx1(c)),DimX(c),Fac(c));
            CurDimX(c)=Fac(c);
        end;      
        if MissingExist
            RedData=missmult(kthFactor',RedData);
        else
            RedData=kthFactor'*RedData;
        end;
        
        if c~=N,
            newi=CurDimX(c+1);
            newj=prod(CurDimX)/CurDimX(c+1);
        else
            newi=CurDimX(1);
            newj=prod(CurDimX)/CurDimX(1);
        end;
        
        RedData=reshape(RedData',newi,newj);
    end;
    G=RedData;
else %oblique factors
    
    LMatTmp=1;
    if Fac_orig(1)==-1,
        LMatTmp=eye(DimX(c));
    else
        LMatTmp=reshape(Factors(FIdx0(1):FIdx1(1)),DimX(1),Fac(1));
    end;    
    
    RMatTmp=1;
    for c=2:N,
        if Fac_orig(c)==-1,
            kthFactor=eye(DimX(c));
        else
            kthFactor=reshape(Factors(FIdx0(c):FIdx1(c)),DimX(c),Fac(c));
        end;    
        RMatTmp=ckron(kthFactor',RMatTmp);
    end;
    
    if MissingExist
        RedData=missmult(pinv(LMatTmp),X);
        RedData=missmult(RedData,pinv(RMatTmp));
    else
        RedData=LMatTmp\X;
        RedData=RedData/RMatTmp;
    end;
    
    G=RedData;
    
end;    

return