#
#   function [X,Y,Eta,Psi,sigma,num_iter]=homals3a(H,p0,maxit,X00,epsilon)
#
#   HOMALS algorithm. Normalized scores: 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.
#   maxit, X00, epsilon, are optional arguments, to control iteration.
#
#   Output:
#   X        scores for individuals  Y   category quantifications
#   Eta      discrimination measures Psi eigenvalues
#   sigma    loss function
#   num_iter number of iterations
#
#
function [X,Y,Eta,Psi,sigma,num_iter]=homals3a(H,p0,maxit,X00,epsilon)
#
   [n,m]=size(H);
# ------------------------------------------------------------------
#
#   Default values for parameters
#
   default_p = min([n,m]);
   default_epsilon=1.0e-9;
   default_maxit=250;
# ------------------------------------------------------------------
#
#   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
   u=rand(n,p);
   default_X0=u-ones(n,1)*(sum(u)/n);
# ------------------------------------------------------------------
#
#   Check if maxit (=maximum number of iterations) was entered.
#   Otherwise, set a default value.
#
   if nargin<3,
      maxit=default_maxit;
   endif
# ------------------------------------------------------------------
#
#   Check if X0 (=initial scores matrix) was entered. If so,
#   check its dimensions and center it.
#   Otherwise, set a default value.
#
   if nargin>=4
      [n1,n2]=size(X00);
      if n1==n && n2<=p,
         p=min([n2,p]);
         X0=X00(:,1:p);
         X0=X0-ones(n,1)*(sum(X0)/n);
      else
         disp('*** HOMALS3a: Wrong sized X0. Set default');
         X0=default_X0;
      endif
   else
      X0=default_X0;
   endif
# ------------------------------------------------------------------
#
#   Check if epsilon was entered.
#   Otherwise, set a default value
#
   if nargin<5
      epsilon=default_epsilon;
   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.
#   D is the the matrix of univariate marginals. D1, its inverse.
#   Initialize X = X0, Y = D1*G'*X
#

   [G,M,mm]=htoind(H);
   D=diag(sum(G));
   D1=D**(-1);
   sqrtn=sqrt(n);
   X=X0;
   Y = D1*G'*X;
# ------------------------------------------------------------------
#
#   Other initializations.
#   X_conv, Y_conv and sigma_conv are indicators of convergence.
#   k is a counter for displaying intermediate results and
#   centering the scores periodically.
#
   Sigma=zeros(1,mm);
   sigma=0;
   X_conv=0;
   Y_conv=0;
   sigma_conv=0;
   k=0;
# ------------------------------------------------------------------
#
#   Iterate
#
    for i=2:maxit
         previous_sigma=sigma;
         previous_X=X;
         previous_Y=Y;
         X=G*Y;
         U=sqrtn*grschm(X);
         X=U;
         k++;
         if k>=20
            format long
            if !sigma_conv
               sigma
            endif
            X=X-ones(n,1)*(sum(X)/n)
            z=input("Press ENTER to continue...");
            format
            k=0;
         endif
         if ssq(X-previous_X)<epsilon
            disp('*** HOMALS3a: Convergence for X attained');
            X_conv=1;
         endif
#   format long
         Y = D1*G'*X;
#   format
         if ssq(Y-previous_Y)<epsilon
            disp('*** HOMALS3a: Convergence for Y attained');
            Y_conv=1;
         endif
#   z=input("Press ENTER to continue...");
         for j=1:mm
            Sigma(j)=ssq(X-G(:,j)*Y(j,:));
         endfor
         sigma=sum(Sigma)/m;
         if abs(sigma-previous_sigma)<epsilon
            if !sigma_conv
              disp('*** HOMALS3a: Convergence for sigma attained');
            endif
            sigma_conv=1;
         endif
         if X_conv
            break
         endif
    endfor
    num_iter=i;
# ------------------------------------------------------------------
#
#   Compute discrimination measures and eigenvalues
#
    Eta=zeros(m,p);
    j2=0;
    for j=1:m
      j1=j2+1;
      j2=j1+M(j)-1;
      Eta(j,:)=diag(Y(j1:j2,:)'*D(j1:j2,j1:j2)*Y(j1:j2,:))'/n;
    endfor
    Psi=sum(Eta)/m;
endfunction