function output=canalysis(N,k,Lr,Lc,Sup,Ls);
% output=canalysis(N,k,Lr,Lc,Sup,Ls);
%
% INPUT:
%
%       N                 Contingence matrix
%       k                 Number of dimensions to retain
%       Lr                Vector of row labels
%       Lc                Vector of column labels
%       Sup               Structure containing supplementary rows or columns (they has to be supplementary rows or columns of the contingence table)
%                         Sup.r    Supplementary rows
%                         Sup.c    Supplementary columns
%       Ls                Structure containing supplementary rows and column labels
%                         Ls.r     Vector of supplementary row labels
%                         Ls.c     Vector of supplementary column labels
%
% OUTPUT:
%
% 		output.N         Contigence table
% 		output.SupRow    Contigence table related to supplementary rows
% 		output.SupCol    Contigence table related to supplementary columns
% 		output.labels.r  Row labels
% 		output.labels.c  Column labels
%       output.labels.sr Row supplementary labels
%       output.labels.sc Column supplementary labels
% 		output.I         Number of rows
% 		output.J         Number of columns
% 		output.n         Grand total
% 		output.P         Correspondence matrix (proportions)
%
% 		output.r         Row masses (vector)
% 		output.Dr        Row masses (diagonal matrix)
% 		output.c         Column masses (vector); c is the centroid of row profiles 
% 		output.Dc        Column masses (diagonal matrix)
% 		output.RP        Rows profiles
% 		output.CP        Column profiles
%
% 		output.K         Maximun number of dimensions
% 		output.k         Number of retained dimensions
%
% 		output.A         Standarized residuals
%
% 		output.chis      Chi-square statistic for the contingence table
% 		output.ti        Total inertia
% 		output.CV        Cramer's V   
% 		output.eig        Principal inertias [Raw Percentage Acumulate_percentage]
%       output.aev        Average explained variance (dimensions explaining less variance should be excluded from the map)
%
% 		output.F          Principal coordinates of rows
% 		output.G          Principal coordinates of columns
% 		output.X          Standard coordinates of rows
% 		output.Y          Standard coordinates of columns
% 		output.H          Symmetrical coordinates of rows
% 		output.Z          Symmetrical coordinates of columns
%
% 		output.ARC        Absolute row contributions: Contributions of points to axes.
% 		output.ACC        Absolute column contributions: Contributions of points to axes.
% 		output.INr        Inertia of each row [raw percentage]
% 		output.INc        Inertia of each column [raw percentage]
% 		output.RRC        Relative row contributions: Contributions of axes to points; squared correlations of points with axes
% 		output.RCC        Relative column contributions: Contributions of axes to points; squared correlations of points with axes
% 		output.QLTr       Measure of quality of display for each row in the reduced k-dimensional map. The qualities are equivalent to the communalities in PCA 
% 		output.QLTc       Measure of quality of display for each column in the reduced k-dimensional map. The qualities are equivalent to the communalities in PCA 
% 		output.Br         Correlation of rows with axes. Similar to component loadings in PCA
% 		output.Bc         Correlation of columns with axes. Similar to component loadings in PCA
%
% 		output.Fsup       Principal coordinates of supplementary columns
% 		output.Xsup       Standard coordinates of supplementary columns
% 		output.Hsup       Symmetrical coordinates of supplementary columns
%
% 		output.Gsup       Principal coordinates of supplementary columns
% 		output.Ysup       Standard coordinates of supplementary columns
% 		output.Zsup       Symmetrical coordinates of supplementary columns
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Size of N
[I,J]=size(N);

output.N = N;             %contigence table
if exist('Sup','var'),    %Contigence tables related to supplementary variables
   if isfield(Sup,'r'), 
      output.SupRow = Sup.r;
  end;
   if isfield(Sup,'c'), 
      output.SupCol = Sup.c;
  end;
end;

if exist('Lr','var'),
    output.labels.r = Lr;     %row labels;
    if ischar(output.labels.r),output.labels.r=cellstr(output.labels.r);end;
else,
    [f1,f2]=size(N);
    output.labels.r=cellstr(num2str([1:f1]'));
end;

if exist('Lc','var'),
    output.labels.c = Lc;     %column labels;
    if ischar(output.labels.c),output.labels.c=cellstr(output.labels.c);end;
else,
    [f1,f2]=size(N);
    output.labels.c=cellstr(num2str([1:f2]'));    
end;

%Supplementary row and column labels
if exist('Ls','var'),
   if isfield(Ls,'r'), 
      output.labels.sr = Ls.r;
      if ischar(output.labels.sr),output.labels.sr=cellstr(output.labels.sr);end;
   end;
   if isfield(Ls,'c'), 
      output.labels.sc = Ls.c;
      if ischar(output.labels.sc),output.labels.sc=cellstr(output.labels.sc);end;
   end;
end;

output.I=I;            %number of rows
output.J=J;            %number of columns

%vectors of ones
Ir = ones(I,1);
Ic = ones(J,1);

%grand total
n=Ir'*N*Ic;
output.n=n;            %grand total

%correspondence matrix P
P = (1/n) * N;
output.P=P;            %Correspondence matrix (proportions)

%row masses
r  = P * Ic ;
Dr = diag(r);
output.r = r;          %row masses (vector)
output.Dr = Dr;        %row masses (diagonal matrix)

%column masses
% Centroid of the row profiles   r' * Dr^(-1) * P = 1' * P = c'

c  = (Ir' * P)';
Dc = diag(c);
output.c = c;          %column masses (vector); c is the centroid of row profiles 
output.Dc = Dc;        %column masses (diagonal matrix)


%Rows profiles
RP = Dr^(-1) * P;
output.RP = RP;        %rows profiles

%Columns profiles
CP = P * Dc^(-1);
output.CP = CP;        %column profiles

%Number fo dimensions
%
K = min(I-1,J-1);
output.K = K;          %maximun number of dimensions

if exist('k','var'),   %number of retained dimensions
   output.k = k;
else,
   k=2;
   output.k = k;
end;

% Standarized residuals
%
% A = Dr^(1/2) * (Dr^(-1) * P - Ir * c') * Dc^(-1/2) = Dr^(-1/2) * (P - r * c') * Dc^(-1/2);  aij = sqrt( (pij - ricj)^2 /ricj ) =(pij - ricj)/sqrt(ricj)
% Chi-square distances
%
A =  Dr^(-1/2) * (P - r * c') * Dc^(-1/2);
output.A= A;           % Standarized residuals

%
% SVD of A = U*L*V'

[U,L,V] = svd(A);
L=L(1:K,1:K);
U=U(:,1:K);
V=V(:,1:K);

% Chi-square statistic for the contingence table
%
%  chis = n * sum_i sum_j (pij - ricj)^2 / ricj
%
chis = Ir'* (A.*A) * Ic * n;
output.chis = chis;    % Chi-square statistic for the contingence table

% Total inertia
%   
% ti = sum_i sum_j (pij - ricj)^2 / ricj = chis/n
%
ti = chis/n;
output.ti = ti;        %total inertia

output.CV = sqrt(chis/(n*(min(I,J)-1)));   %Cramer's V   

% Principal inertias 
% 
% eigenvalues of A'A
% percentages of the total inertia 

api = [];tmp=0;
for i=1:K,
    tmp = tmp + L(i,i)*L(i,i)/ti;
    api = [api; tmp];
end;
output.eig = [diag(L.*L) diag(L.*L / ti) api] ;      % Principal inertias [Raw Percentage Acumulate_percentage]

output.aev = 100/K;                        %Average explained variance (dimensions explaining less variance should be excluded from the map)

% Principal coordinates of rows
% 
F = Dr^(-1/2) * U*L;
output.F = F;                               % Principal coordinates of rows

% Principal coordinates of columns
% 
G = Dc^(-1/2) * V*L;                        
output.G = G;                               % Principal coordinates of columns

% Standard coordinates of rows
%
% X = F * L^(-1) = Dr^(-1/2) * U * L * L^(-1) = Dr^(-1/2) * U  
%
X = Dr^(-1/2) * U ;
output.X = X;                               % Standard coordinates of rows

% Standard coordinates of columns
%
% X = G * L^(-1) = Dc^(-1/2) * V * L * L^(-1) = Dc^(-1/2) * V  
%
Y = Dc^(-1/2) * V ;
output.Y=Y;                                 % Standard coordinates of columns

% Simmetrical coordinates of rows
% 
H = Dr^(-1/2) * U*L^(1/2);
output.H = H;                               % Symmetrical coordinates of rows

% Simmetrical coordinates of columns
% 
Z = Dc^(-1/2) * V*L^(1/2);
output.Z=Z;                                 % Symmetrical coordinates of columns


% Contributions of rows to inertia 
% 
ARC = Dr*X.*X;
output.ARC = ARC;                          % Absolute row contributions: Contributions of points to axes.

% Contributions of columns to inertia
% 
ACC = Dc*Y.*Y;
output.ACC = ACC;                          % Absolute column contributions: Contributions of points to axes.

%Inertia of each row
%
INr = Dr*(F.*F)*ones(K,1);
%Inertia relative to their total (used in the output)
INr_p=INr / trace(L.*L);
output.INr = [INr INr_p];                   %Inertia of each row [raw percentage]

%Inertia of each column
%
INc = Dc*(G.*G)*ones(K,1);
%Inertia relative to their total (used in the output)
INc_p=INc / trace(L.*L);
output.INc = [INc INc_p];                   %Inertia of each column [raw percentage]

%Squared correlations of rows with axes
%
QCOR_r = diag(INr)^(-1)*Dr*(F.*F);
output.RRC=QCOR_r;                          %Relative row contributions: Contributions of axes to points; squared correlations of points with axes

%Squared correlations of columns with axes
%
QCOR_c = diag(INc)^(-1)*Dc*(G.*G);
output.RCC=QCOR_c;                          %Relative columns contributions: Contributions of axes to points; squared correlations of points with axes

%Measure of quality of display for each row in the reduced k-dimensional map. The qualities are equivalent to the communalities in PCA 
%
QLTr = QCOR_r(:,1:k)*ones(k,1);
output.QLTr = QLTr;                         %Measure of quality of display for each row in the reduced k-dimensional map. The qualities are equivalent to the communalities in PCA 

%Measure of quality of display for each column in the reduced k-dimensional map. The qualities are equivalent to the communalities in PCA 
%
QLTc = QCOR_c(:,1:k)*ones(k,1);
output.QLTc = QLTc;                         %Measure of quality of display for each column in the reduced k-dimensional map. The qualities are equivalent to the communalities in PCA 

%Correlation of rows with axes. Similar to component loadings in PCA
%
Br = sign(F).*sqrt(QCOR_r);
output.Br = Br;                             %Correlation of rows with axes. Similar to component loadings in PCA

%Correlation of columns with axes. Similar to component loadings in PCA
%
Bc = sign(G).*sqrt(QCOR_c);
output.Bc = Bc;                             %Correlation of columns with axes. Similar to component loadings in PCA

if exist('Sup','var'),
    %Supplementary rows
    if isfield(Sup,'r'),
        h=Sup.r/(diag(sum(Sup.r)));
        Fsup=h*Y;                              %Principal coordinates of supplementary rows
        Xsup=h*Y*L^(-1);                       %Standard coordinates of supplementary rows
        Hsup=h*Y*L^(-1/2);                     %Symmetrical coordinates of supplementary rows
        output.Fsup=Fsup;
        output.Xsup=Xsup;
        output.Hsup=Hsup;                                 
%       rrc=(Gsup(:,1:k).*Gsup(:,1:k))/trace(Gsup(:,1:k)'*Gsup(:,1:k));
    end;
    
    %Supplementary columns
    if isfield(Sup,'c'),
        h=Sup.c/(diag(sum(Sup.c)));
        Gsup=h'*X;                              %Principal coordinates of supplementary columns
        Ysup=h'*X*L^(-1);                       %Standard coordinates of supplementary columns
        Zsup=h'*X*L^(-1/2);                     %Symmetrical coordinates of supplementary columns
        output.Gsup=Gsup;
        output.Ysup=Ysup;
        output.Zsup=Zsup;                                 
%       rrc=(Gsup(:,1:k).*Gsup(:,1:k))/trace(Gsup(:,1:k)'*Gsup(:,1:k));
    end;
end

return;
