#
#   function [X,Y,Psi]=multica1(H,p0)
#
#   Multiple correspondence analysis.
#   Normalized scores algorithm: X'*X=n*eye.
#
#   Input:
#   H       Data matrix. Columns contain the values of integer-coded
#           categorical variables.
#           Minimum value = 1. Maximum value = number of categories.
#   p0      Requested number of dimensions
#
#   Output:
#   X        scores for individuals
#   Y        category quantifications
#   Psi      discrimination measures
#   sigma    loss function
#
#
function [X,Y,Psi]=multica1(H,p0)
#
   [n,m]=size(H);
# ------------------------------------------------------------------
#
#   Default values for parameters
#
   default_p = min([n,m])-1;
# ------------------------------------------------------------------
#
#   Check if p0 (= requested number of solutions) was entered.
#   Otherwise, set a default value.
#
   if nargin>=2
     p=min([p0,default_p]);
   else
     p=default_p;
   endif
# ------------------------------------------------------------------
#
#   Compute the indicator matrix G. The function htoind also
#   returns a row vector M which contains the number of columns
#   of G for each variable, and an integer mm = sum(M) = total
#   number of columns of G.
#   C is the matrix of bivariate marginals (Burt's matrix),
#   D is the matrix of univariate marginals,
#   D1   = D^(-1),
#   Dm12 = D^(-1/2).
#
   [G,M,mm]=htoind(H);
   C=G'*G;
   d=sum(G);
   if any(d)==0
      disp("*** MultiCA: found null marginal");
      return
   endif
   D=diag(d);
   D1=diag(d.^(-1));
   Dm12=diag(d.^(-1/2));
   u=ones(mm,1);
   [W,Psi]=eigsort(Dm12*(C-D*u*u'*D/n)*Dm12);
   [Psi,pc,ac,npos]=percent(Psi);
   if p<=npos
      Psi=Psi(1:p);
      pc=pc(1:p);
      ac=ac(1:p);
      W=W(:,1:p);
   else
      disp("*** MultiCA: requested number of solutions too large");
      p=npos;
   endif
   SqrtPsi =diag(Psi.^(1/2));
   Psim1=diag(Psi.^(-1));
   Y=sqrt(n)*Dm12*W*SqrtPsi;
   X=G*Y*Psim1;
endfunction