#
#   function [A,X,sigma,num_iter]=als4a(H,p0,maxit0,A00,epsilon0)
#
#   Normalized scores algorithm. Multiple solutions. Columns of
#   H have unequal variances. Gram-Schmidt orthonormalization of X.
#
#   Input:
#   H       Data matrix.
#   p0      Requested number of solutions
#   maxit0, A00, epsilon0, are optional arguments, to control iteration.
#
#   Output: WARNING! : X and A do NOT contain successive intermediate
#           steps, but different orthogonal solutions.
#   A = weights. X = scores. sigma = loss function.
#   num_iter = number of iterations.
#
function [A,X,sigma,num_iter]=als4a(H,p0,maxit,A00,epsilon)
#
   [n,m]=size(H);
# ------------------------------------------------------------------
#
#   Default values for parameters
#
   default_p = min([n,m]);
   default_epsilon=1.0e-10;
   default_maxit=50;
# ------------------------------------------------------------------
#
#   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
   default_A0=rand(m,p);
# ------------------------------------------------------------------
#
#   Check if maxit0 (=maximum number of iterations) was entered.
#   Otherwise, set a default value.
#
   if nargin>=3,
      maxit=maxit0;
   else
      maxit=default_maxit;
   endif
# ------------------------------------------------------------------
#
#   Check if A0 (=initial weights matrix) was entered. If so,
#   check its dimensions. Otherwise, set a default value.
#
   if nargin>=4
      [m1,m2]=size(A00);
      if m1==m && m2<=p,
         p=min([m2,p]);
         A0=A00(:,1:p);
      else
         disp('*** als4a: Wrong sized A0. Set default');
         A0=default_A0;
      endif
   else
      A0=default_A0;
   endif
# ------------------------------------------------------------------
#
#   Check if epsilon0 (=tolerance for convergence) was entered.
#   Otherwise, set a default value.
#
   if nargin>=5
      epsilon=epsilon0;
   else
      epsilon=default_epsilon;
   endif
# ------------------------------------------------------------------
#
#   Compute D (the diagonal of H'*H).
#   Initialize A = A0 and X = 0. Other initializations: Sigma is an
#   auxiliary [1,m] row vector, used to compute the loss function.
#   A_conv and sigma_conv are indicators of convergence.
#
   D=diag(diag(H'*H));
   D1=inv(D);
   X=zeros(n,p);
   A=A0;
   Sigma=zeros(1,m);
   sigma=0;
   A_conv=0;
   sigma_conv=0;
# ------------------------------------------------------------------
#
#   Iterate
#
    for i=1:maxit
      previous_sigma=sigma;
      previous_A=A;
# ------------------------------------------------------------------
#
#   Compute next X.
#
      X=H*A;
# ------------------------------------------------------------------
#
#   Check convergence for sigma.
#
      for j=1:m
         Sigma(j)=ssq(X-H(:,j)*A(j,:));
      endfor
      sigma=sum(Sigma)/m;
      if abs(sigma-previous_sigma)<epsilon
         if !sigma_conv
            disp('*** als4a: Convergence for sigma attained');
         endif
         sigma_conv=1;
      else
         format long;
         sigma
         format
         input("Press ENTER to continue...");
      endif
# ------------------------------------------------------------------
#
#   Orthonormalization of X. Either using Gram-Schmidt, or svd or qr.
#   Uncomment only one option.
#
      U=grschm(X);
#   [U,L,V]=svd(X,0);
#   [U,R]=qr(X,0);
      X=U;
# ------------------------------------------------------------------
#
#   Compute next A.
#
      A=D1*H'*X;
# ------------------------------------------------------------------
#
#   Check convergence for A. Display A if sigma already converged.
#
      if sigma_conv
         format long
         A
         format
         input("Press ENTER to continue...");
      endif
      if ssq(A-previous_A)<epsilon
         disp('*** als4a: Convergence for A attained');
         A_conv=1;
      endif
# ------------------------------------------------------------------
      if A_conv
         break
      endif
    endfor
    num_iter=i;
endfunction