#
#   function [A,X,sigma,num_iter]=als4b(H,p0,maxit0,X00,epsilon0)
#
#   Normalized weights algorithm. Multiple solutions. Columns of
#   H have unequal variances. Gram-Schmidt orthonormalization of A.
#
#   Input:
#   H  = Data matrix.
#   p0 = Requested number of solutions.
#   maxit0, X00, 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]=als4b(H,p0,maxit0,X00,epsilon0)
#
   [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_X0=rand(n,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 X0 (=initial scores matrix) was entered. If so,
#   check its dimensions. 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);
      else
         disp('*** als4b: Wrong sized X0. Set default');
         X0=default_X0;
      endif
   else
      X0=default_X0;
   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 X = X0 and A = 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=X0;
   A=zeros(m,p);
   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 A.
#
      A=D1*H'*X;
# ------------------------------------------------------------------
#
#   Orthonormalization of A using Gram-Schmidt.
#
         U=grschm(A,D);
         A=U;
# ------------------------------------------------------------------
#
#   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('*** als4b: Convergence for A attained');
            A_conv=1;
         endif
# ------------------------------------------------------------------
#
#   Compute next X.
#
         X=H*A/m;
# ------------------------------------------------------------------
#
#   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('*** als4b: Convergence for sigma attained');
         endif
         sigma_conv=1;
      else
         format long;
         sigma
         format
         input("Press ENTER to continue...");
      endif
# ------------------------------------------------------------------
      if A_conv
         break
      endif
    endfor
    num_iter=i;
endfunction