#
#   function [X,Y,Sigma]=homals2(H,maxit,x00,epsilon)
#
#   HOMALS algorithm.
#   First solution (dimension). Keeping intermediate results.
#   Normalization y'*D*y=1.
#
#   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
#   Sigma   loss function
#
#           these three matrices contain all the intermediate
#           values, till convergence. The i-th column corresponds
#           to the i-th iteration.
#
function [X,Y,Sigma]=homals2(H,maxit,a00,epsilon)
#
   [n,m]=size(H);
# ------------------------------------------------------------------
#
#   Default values for parameters
#
   default_epsilon=1.0e-4;
   default_maxit=50;
   u=rand(n,1);
   default_x0=u-(sum(u)/n)*ones(n,1);
# ------------------------------------------------------------------
#
#   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.
#   D1 is the inverse of the matrix D of univariate marginals.
#

   [G,M,mm]=htoind(H);
   d=sum(G);
   D=diag(d);
   D1=diag(d.**(-1));
# ------------------------------------------------------------------
#
#   Check if maxit (=maximum number of iterations) was entered.
#   Otherwise, set a default value.
#
#
   if nargin==1,
      maxit=default_maxit;
   endif
# ------------------------------------------------------------------
#
#   Check if x0 (=initial scores vector) was entered.
#   if so, check its dimensions.
#   Otherwise, set a default value.
#
   if nargin>=3
      [n1,n2]=size(x00);
      if n1==n && n2==1,
      x0=x00;
      else
         if n1==1 && n2==n,
            x0=x00';
         else
            disp('*** HOMALS2: Wrong sized x0. Set default');
            x0=default_x0;
         endif
      endif
   else
      x0=default_x0;
   endif
# ------------------------------------------------------------------
#
#   Check if epsilon was entered.
#   Otherwise, set a default value
#
   if nargin<4
      epsilon=default_epsilon;
   endif
# ------------------------------------------------------------------
#
#   Initialize matrices for storing all the iterations.
#
   X=zeros(n,maxit);
   Y=zeros(mm,maxit);
   Sigma=zeros(1,maxit);
# ------------------------------------------------------------------
#
#   Initialize x and y
#
   x=x0;
   y=D1*G'*x;
   sigma=norm(x*y'-G,'fro')/m;
   X(:,1)=x;
   Y(:,1)=y;
# ------------------------------------------------------------------
#
#   Iterate
#
    for i=2:maxit
         y=D1*G'*x;
         ly=sqrt(y'*D*y);
         y=y/ly;
         x=G*y/m;
         sigma=norm(x*y'-G,'fro')/m;
         X(:,i)=x;
         Y(:,i)=y;
         Sigma(i)=sigma;
         if norm(Y(:,i-1)-y)<epsilon
            disp('*** HOMALS2: Convergence for y attained');
            break
         endif
    endfor
    num_iter=i;
    X=X(:,1:num_iter);
    Y=Y(:,1:num_iter);
    Sigma=Sigma(:,1:num_iter);
endfunction