source: src/main/java/weka/classifiers/functions/pace/PaceMatrix.java @ 16

/*
 *    This program is free software; you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation; either version 2 of the License, or (at
 *    your option) any later version.
 *
 *    This program is distributed in the hope that it will be useful, but
 *    WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *    General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with this program; if not, write to the Free Software
 *    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.  */

/*
 *    PaceMatrix.java
 *    Copyright (C) 2002 University of Waikato, Hamilton, New Zealand
 *
 */

package weka.classifiers.functions.pace;

import weka.core.RevisionUtils;
import weka.core.matrix.DoubleVector;
import weka.core.matrix.FlexibleDecimalFormat;
import weka.core.matrix.IntVector;
import weka.core.matrix.Matrix;
import weka.core.matrix.Maths;

import java.util.Random;
import java.text.DecimalFormat;

/**
 * Class for matrix manipulation used for pace regression. <p>
 *
 * REFERENCES <p>
 * 
 * Wang, Y. (2000). "A new approach to fitting linear models in high
 * dimensional spaces." PhD Thesis. Department of Computer Science,
 * University of Waikato, New Zealand. <p>
 * 
 * Wang, Y. and Witten, I. H. (2002). "Modeling for optimal probability
 * prediction." Proceedings of ICML'2002. Sydney. <p>
 *
 * @author Yong Wang (yongwang@cs.waikato.ac.nz)
 * @version $Revision: 1.6 $
 */
public class PaceMatrix 
  extends Matrix {
  
  /** for serialization */
  static final long serialVersionUID = 2699925616857843973L;
    
  /* ------------------------
     Constructors
     * ------------------------ */
  
  /** Construct an m-by-n PACE matrix of zeros. 
      @param m    Number of rows.
      @param n    Number of colums.
  */
  public PaceMatrix( int m, int n ) {
    super( m, n );
  }

  /** Construct an m-by-n constant PACE matrix.
      @param m    Number of rows.
      @param n    Number of colums.
      @param s    Fill the matrix with this scalar value.
  */
  public PaceMatrix( int m, int n, double s ) {
    super( m, n, s );
  }
    
  /** Construct a PACE matrix from a 2-D array.
      @param A    Two-dimensional array of doubles.
      @throws  IllegalArgumentException All rows must have the same length
  */
  public PaceMatrix( double[][] A ) {
    super( A );
  }

  /** Construct a PACE matrix quickly without checking arguments.
      @param A    Two-dimensional array of doubles.
      @param m    Number of rows.
      @param n    Number of colums.
  */
  public PaceMatrix( double[][] A, int m, int n ) {
    super( A, m, n );
  }
    
  /** Construct a PaceMatrix from a one-dimensional packed array
      @param vals One-dimensional array of doubles, packed by columns (ala Fortran).
      @param m    Number of rows.
      @throws  IllegalArgumentException Array length must be a multiple of m.
  */
  public PaceMatrix( double vals[], int m ) {
    super( vals, m );
  }
    
  /** Construct a PaceMatrix with a single column from a DoubleVector 
      @param v    DoubleVector
  */
  public PaceMatrix( DoubleVector v ) {
    this( v.size(), 1 );
    setMatrix( 0, v.size()-1, 0, v );
  }
    
  /** Construct a PaceMatrix from a Matrix 
      @param X    Matrix 
  */
  public PaceMatrix( Matrix X ) {
    super( X.getRowDimension(), X.getColumnDimension() );
    A = X.getArray();
  }
    
  /* ------------------------
     Public Methods
     * ------------------------ */

  /** Set the row dimenion of the matrix
   *  @param rowDimension the row dimension
   */
  public void setRowDimension( int rowDimension ) 
  {
    m = rowDimension;
  }

  /** Set the column dimenion of the matrix
   *  @param columnDimension the column dimension
   */
  public void setColumnDimension( int columnDimension ) 
  {
    n = columnDimension;
  }

  /** 
   * Clone the PaceMatrix object.
   * 
   * @return the clone
   */
  public Object clone () {
    PaceMatrix X = new PaceMatrix(m,n);
    double[][] C = X.getArray();
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        C[i][j] = A[i][j];
      }
    }
    return (Object) X;
  }
    
  /** Add a value to an element and reset the element
   *  @param i    the row number of the element
   *  @param j    the column number of the element
   *  @param s    the double value to be added with
   */
  public void setPlus(int i, int j, double s) {
    A[i][j] += s;
  }

  /** Multiply a value with an element and reset the element
   *  @param i    the row number of the element
   *  @param j    the column number of the element
   *  @param s    the double value to be multiplied with
   */
  public void setTimes(int i, int j, double s) {
    A[i][j] *= s;
  }

  /** Set the submatrix A[i0:i1][j0:j1] with a same value 
   *  @param i0 the index of the first element of the column
   *  @param i1 the index of the last element of the column
   *  @param j0 the index of the first column
   *  @param j1 the index of the last column
   *  @param s the value to be set to
   */
  public void setMatrix( int i0, int i1, int j0, int j1, double s ) {
    try {
      for( int i = i0; i <= i1; i++ ) {
        for( int j = j0; j <= j1; j++ ) {
          A[i][j] = s;
        }
      }
    } catch( ArrayIndexOutOfBoundsException e ) {
      throw new ArrayIndexOutOfBoundsException( "Index out of bounds" );
    }
  }
  
  /** Set the submatrix A[i0:i1][j] with the values stored in a
   *  DoubleVector
   *  @param i0 the index of the first element of the column
   *  @param i1 the index of the last element of the column
   *  @param j  the index of the column
   *  @param v the vector that stores the values*/
  public void setMatrix( int i0, int i1, int j, DoubleVector v ) {
    for( int i = i0; i <= i1; i++ ) {
      A[i][j] = v.get(i-i0);
    }
  }

  /** Set the whole matrix from a 1-D array 
   *  @param v    1-D array of doubles
   *  @param columnFirst   Whether to fill the column first or the row.
   *  @throws  ArrayIndexOutOfBoundsException Submatrix indices
   */
  public void setMatrix ( double[] v, boolean columnFirst ) {
    try {
      if( v.length != m * n ) 
        throw new IllegalArgumentException("sizes not match.");
      int i, j, count = 0;
      if( columnFirst ) {
        for( i = 0; i < m; i++ ) {
          for( j = 0; j < n; j++ ) {
            A[i][j] = v[count];
            count ++;
          }
        }
      }
      else {
        for( j = 0; j < n; j++ ) {
          for( i = 0; i < m; i++ ){
            A[i][j] = v[count];
            count ++;
          }
        }
      }

    } catch( ArrayIndexOutOfBoundsException e ) {
      throw new ArrayIndexOutOfBoundsException( "Submatrix indices" );
    }
  }

  /** Returns the maximum absolute value of all elements 
      @return the maximum value
  */
  public double maxAbs () {
    double ma = Math.abs(A[0][0]);
    for (int j = 0; j < n; j++) {
      for (int i = 0; i < m; i++) {
        ma = Math.max(ma, Math.abs(A[i][j]));
      }
    }
    return ma;
  }

  /** Returns the maximum absolute value of some elements of a column,
      that is, the elements of A[i0:i1][j].
      @param i0 the index of the first element of the column
      @param i1 the index of the last element of the column
      @param j  the index of the column
      @return the maximum value */
  public double maxAbs ( int i0, int i1, int j ) {
    double m = Math.abs(A[i0][j]);
    for (int i = i0+1; i <= i1; i++) {
      m = Math.max(m, Math.abs(A[i][j]));
    }
    return m;
  }

  /** Returns the minimum absolute value of some elements of a column,
      that is, the elements of A[i0:i1][j].
      @param i0 the index of the first element of the column
      @param i1 the index of the last element of the column
      @param column the index of the column
      @return the minimum value 
  */
  public double minAbs ( int i0, int i1, int column ) {
    double m = Math.abs(A[i0][column]);
    for (int i = i0+1; i <= i1; i++) {
      m = Math.min(m, Math.abs(A[i][column]));
    }
    return m;
  }
    
  /** Check if the matrix is empty
   *   @return true if the matrix is empty
   */
  public boolean  isEmpty(){
    if(m == 0 || n == 0) return true;
    if(A == null) return true;
    return false;
  }
    
  /** Return a DoubleVector that stores a column of the matrix 
   *  @param j the index of the column
   *  @return the column
   */
  public DoubleVector  getColumn( int j ) {
    DoubleVector v = new DoubleVector( m );
    double [] a = v.getArray();
    for(int i = 0; i < m; i++)
      a[i] = A[i][j];
    return v;
  }

  /** Return a DoubleVector that stores some elements of a column of the
   *  matrix 
   *  @param i0 the index of the first element of the column
   *  @param i1 the index of the last element of the column
   *  @param j  the index of the column
   *  @return the DoubleVector
   */
  public DoubleVector  getColumn( int i0, int i1, int j ) {
    DoubleVector v = new DoubleVector( i1-i0+1 );
    double [] a = v.getArray();
    int count = 0;
    for( int i = i0; i <= i1; i++ ) {
      a[count] = A[i][j];
      count++;
    }
    return v;
  }
  
  
  /** Multiplication between a row (or part of a row) of the first matrix
   *  and a column (or part or a column) of the second matrix.
   *  @param i the index of the row in the first matrix
   *  @param j0 the index of the first column in the first matrix
   *  @param j1 the index of the last column in the first matrix
   *  @param B the second matrix
   *  @param l the index of the column in the second matrix
   *  @return the result of the multiplication
   */
  public double  times( int i, int j0, int j1, PaceMatrix B, int l ) {
    double s = 0.0;
    for(int j = j0; j <= j1; j++ ) {
      s += A[i][j] * B.A[j][l];
    }
    return s;
  }
  
  /** Decimal format for converting a matrix into a string
   *  @return the default decimal format
   */
  protected DecimalFormat []  format() {
    return format(0, m-1, 0, n-1, 7, false );
  }
  
  /** Decimal format for converting a matrix into a string
   *  @param digits the number of digits
   *  @return the decimal format
   */
  protected DecimalFormat []  format( int digits ) {
    return format(0, m-1, 0, n-1, digits, false);
  }

  /** Decimal format for converting a matrix into a string
   *  @param digits the number of digits
   *  @param trailing
   *  @return the decimal format
   */
  protected DecimalFormat []  format( int digits, boolean trailing ) {
    return format(0, m-1, 0, n-1, digits, trailing);
  }
  
  /** Decimal format for converting a matrix into a string
   *  @param i0
   *  @param i1
   *  @param j
   *  @param digits the number of digits
   *  @param trailing
   *  @return the decimal format
   */
  protected DecimalFormat  format(int i0, int i1, int j, int digits, 
                                  boolean trailing) {
    FlexibleDecimalFormat df = new FlexibleDecimalFormat(digits, trailing);
    df.grouping( true );
    for(int i = i0; i <= i1; i ++ )
      df.update( A[i][j] );
    return df;
  }
  
  /** Decimal format for converting a matrix into a string
   *  @param i0
   *  @param i1
   *  @param j0
   *  @param j1
   *  @param trailing
   *  @param digits the number of digits
   *  @return the decimal format
   */
  protected DecimalFormat []  format(int i0, int i1, int j0, int j1, 
                                     int digits, boolean trailing) {
    DecimalFormat [] f = new DecimalFormat[j1-j0+1];
    for( int j = j0; j <= j1; j++ ) {
      f[j] = format(i0, i1, j, digits, trailing);
    }
    return f;
  }
  
  /** 
   * Converts matrix to string
   * 
   * @return the matrix as string
   */ 
  public String  toString() {
    return toString( 5, false );
  }
  
  /** 
   * Converts matrix to string
   * 
   * @param digits number of digits after decimal point
   * @param trailing true if trailing zeros are padded
   * @return the matrix as string
   */ 
  public String  toString( int digits, boolean trailing ) {
    
    if( isEmpty() ) return "null matrix";
    
    StringBuffer text = new StringBuffer();
    DecimalFormat [] nf = format( digits, trailing );
    int numCols = 0;
    int count = 0;
    int width = 80;
    int lenNumber;
    
    int [] nCols = new int[n];
    int nk=0;
    for( int j = 0; j < n; j++ ) {
      lenNumber = nf[j].format( A[0][j]).length(); 
      if( count + 1 + lenNumber > width -1 ) {
        nCols[nk++]  = numCols;
        count = 0;
        numCols = 0;
      }
      count += 1 + lenNumber;
      ++numCols;
    }
    nCols[nk] = numCols;
    
    nk = 0;
    for( int k = 0; k < n; ) {
      for( int i = 0; i < m; i++ ) {
        for( int j = k; j < k + nCols[nk]; j++)
          text.append( " " + nf[j].format( A[i][j]) );
        text.append("\n");
      }
      k += nCols[nk];
      ++nk;
      text.append("\n");
    }
    
    return text.toString();
  }
  
  /** Squared sum of a column or row in a matrix
   * @param j the index of the column or row
   * @param i0 the index of the first element
   * @param i1 the index of the last element
   * @param col if true, sum over a column; otherwise, over a row
   * @return the squared sum
   */
  public double sum2( int j, int i0, int i1, boolean col ) {
    double s2 = 0;
    if( col ) {   // column 
      for( int i = i0; i <= i1; i++ ) 
        s2 += A[i][j] * A[i][j];
    }
    else {
      for( int i = i0; i <= i1; i++ ) 
        s2 += A[j][i] * A[j][i];
    }
    return s2;
  }
  
  /** Squared sum of columns or rows of a matrix
   * @param col if true, sum over columns; otherwise, over rows
   * @return the squared sum
   */
  public double[] sum2( boolean col ) {
    int l = col ? n : m;
    int p = col ? m : n;
    double [] s2 = new double[l];
    for( int i = 0; i < l; i++ ) 
      s2[i] = sum2( i, 0, p-1, col );
    return s2;
  }

  /** Constructs single Householder transformation for a column
   *
   @param j    the index of the column
   @param k    the index of the row
   @return     d and q 
  */
  public double [] h1( int j, int k ) {
    double dq[] = new double[2];
    double s2 = sum2(j, k, m-1, true);  
    dq[0] = A[k][j] >= 0 ? - Math.sqrt( s2 ) : Math.sqrt( s2 );
    A[k][j] -= dq[0];
    dq[1] = A[k][j] * dq[0];
    return dq;
  }
  
  /** Performs single Householder transformation on one column of a matrix
   *
   @param j    the index of the column 
   @param k    the index of the row
   @param q    q = - u'u/2; must be negative
   @param b    the matrix to be transformed
   @param l    the column of the matrix b
  */
  public void h2( int j, int k, double q, PaceMatrix b, int l ) {
    double s = 0, alpha;
    for( int i = k; i < m; i++ )
      s += A[i][j] * b.A[i][l];
    alpha = s / q;
    for( int i = k; i < m; i++ )
      b.A[i][l] += alpha * A[i][j];
  }
  
  /** Constructs the Givens rotation
   *  @param a 
   *  @param b
   *  @return a double array that stores the cosine and sine values
   */
  public double []  g1( double a, double b ) {
    double cs[] = new double[2];
    double r = Maths.hypot(a, b);
    if( r == 0.0 ) {
      cs[0] = 1;
      cs[1] = 0;
    }
    else {
      cs[0] = a / r;
      cs[1] = b / r;
    }
    return cs;
  }
  
  /** Performs the Givens rotation
   * @param cs a array storing the cosine and sine values
   * @param i0 the index of the row of the first element
   * @param i1 the index of the row of the second element
   * @param j the index of the column
   */
  public void  g2( double cs[], int i0, int i1, int j ){
    double w =   cs[0] * A[i0][j] + cs[1] * A[i1][j];
    A[i1][j] = - cs[1] * A[i0][j] + cs[0] * A[i1][j];
    A[i0][j] = w;
  }
  
  /** Forward ordering of columns in terms of response explanation.  On
   *  input, matrices A and b are already QR-transformed. The indices of
   *  transformed columns are stored in the pivoting vector.
   *  
   *@param b     the PaceMatrix b
   *@param pvt   the pivoting vector
   *@param k0    the first k0 columns (in pvt) of A are not to be changed
   **/
  public void forward( PaceMatrix b, IntVector pvt, int k0 ) {
    for( int j = k0; j < Math.min(pvt.size(), m); j++ ) {
      steplsqr( b, pvt, j, mostExplainingColumn(b, pvt, j), true );
    }
  }

  /** Returns the index of the column that has the largest (squared)
   *  response, when each of columns pvt[ks:] is moved to become the
   *  ks-th column. On input, A and b are both QR-transformed.
   *  
   * @param b    response
   * @param pvt  pivoting index of A
   * @param ks   columns pvt[ks:] of A are to be tested 
   * @return the index of the column
   */
  public int  mostExplainingColumn( PaceMatrix b, IntVector pvt, int ks ) {
    double val;
    int [] p = pvt.getArray();
    double ma = columnResponseExplanation( b, pvt, ks, ks );
    int jma = ks;
    for( int i = ks+1; i < pvt.size(); i++ ) {
      val = columnResponseExplanation( b, pvt, i, ks );
      if( val > ma ) {
        ma = val;
        jma = i;
      }
    }
    return jma;
  }
  
  /** Backward ordering of columns in terms of response explanation.  On
   *  input, matrices A and b are already QR-transformed. The indices of
   *  transformed columns are stored in the pivoting vector.
   * 
   *  A and b must have the same number of rows, being the (pseudo-)rank. 
   *  
   * @param b     PaceMatrix b
   * @param pvt   pivoting vector
   * @param ks    number of QR-transformed columns; psuedo-rank of A 
   * @param k0    first k0 columns in pvt[] are not to be ordered.
   */
  public void backward( PaceMatrix b, IntVector pvt, int ks, int k0 ) {
    for( int j = ks; j > k0; j-- ) {
      steplsqr( b, pvt, j, leastExplainingColumn(b, pvt, j, k0), false );
    }
  }

  /** Returns the index of the column that has the smallest (squared)
   *  response, when the column is moved to become the (ks-1)-th
   *  column. On input, A and b are both QR-transformed.
   *  
   * @param b    response
   * @param pvt  pivoting index of A
   * @param ks   psudo-rank of A
   * @param k0   A[][pvt[0:(k0-1)]] are excluded from the testing.
   * @return the index of the column
   */
  public int  leastExplainingColumn( PaceMatrix b, IntVector pvt, int ks, 
                                     int k0 ) {
    double val;
    int [] p = pvt.getArray();
    double mi = columnResponseExplanation( b, pvt, ks-1, ks );
    int jmi = ks-1;
    for( int i = k0; i < ks - 1; i++ ) {
      val = columnResponseExplanation( b, pvt, i, ks );
      if( val <= mi ) {
        mi = val;
        jmi = i;
      }
    }
    return jmi;
  }
  
  /** Returns the squared ks-th response value if the j-th column becomes
   *  the ks-th after orthogonal transformation.  A[][pvt[ks:j]] (or
   *  A[][pvt[j:ks]], if ks > j) and b[] are already QR-transformed
   *  on input and will remain unchanged on output.
   *
   *  More generally, it returns the inner product of the corresponding
   *  row vector of the response PaceMatrix. (To be implemented.)
   *
   *@param b    PaceMatrix b
   *@param pvt  pivoting vector
   *@param j    the column A[pvt[j]][] is to be moved
   *@param ks   the target column A[pvt[ks]][]
   *@return     the squared response value */
  public double  columnResponseExplanation( PaceMatrix b, IntVector pvt,
                                            int j, int ks ) {
    /*  Implementation: 
     *
     *  If j == ks - 1, returns the squared ks-th response directly.
     *
     *  If j > ks -1, returns the ks-th response after
     *  Householder-transforming the j-th column and the response.
     *
     *  If j < ks - 1, returns the ks-th response after a sequence of
     *  Givens rotations starting from the j-th row. */

    int k, l;
    double [] xxx = new double[n];
    int [] p = pvt.getArray();
    double val;
    
    if( j == ks -1 ) val = b.A[j][0];
    else if( j > ks - 1 ) {
      int jm = Math.min(n-1, j);
      DoubleVector u = getColumn(ks,jm,p[j]);
      DoubleVector v = b.getColumn(ks,jm,0);
      val = v.innerProduct(u) / u.norm2();
    }
    else {                 // ks > j
      for( k = j+1; k < ks; k++ ) // make a copy of A[j][]
        xxx[k] = A[j][p[k]];
      val = b.A[j][0];
      double [] cs;
      for( k = j+1; k < ks; k++ ) {
        cs = g1( xxx[k], A[k][p[k]] );
        for( l = k+1; l < ks; l++ ) 
          xxx[l] = - cs[1] * xxx[l] + cs[0] * A[k][p[l]];
        val = - cs[1] * val + cs[0] * b.A[k][0];
      }
    }
    return val * val;  // or inner product in later implementation???
  }

  /** 
   * QR transformation for a least squares problem<br/>
   *            A x = b<br/>
   * implicitly both A and b are transformed. pvt.size() is the psuedo-rank of 
   * A.
   *  
   * @param b    PaceMatrix b
   * @param pvt  pivoting vector
   * @param k0   the first k0 columns of A (indexed by pvt) are pre-chosen. 
   *            (But subject to rank examination.) 
   * 
   *            For example, the constant term may be reserved, in which
   *            case k0 = 1.
   **/
  public void  lsqr( PaceMatrix b, IntVector pvt, int k0 ) {
    final double TINY = 1e-15;
    int [] p = pvt.getArray();
    int ks = 0;  // psuedo-rank
    for(int j = 0; j < k0; j++ )   // k0 pre-chosen columns
      if( sum2(p[j],ks,m-1,true) > TINY ){ // large diagonal element 
        steplsqr(b, pvt, ks, j, true);
        ks++;
      }
      else {                     // collinear column
        pvt.shiftToEnd( j );
        pvt.setSize(pvt.size()-1);
        k0--;
        j--;
      }
        
    // initial QR transformation
    for(int j = k0; j < Math.min( pvt.size(), m ); j++ ) {
      if( sum2(p[j], ks, m-1, true) > TINY ) { 
        steplsqr(b, pvt, ks, j, true);
        ks++;
      }
      else {                     // collinear column
        pvt.shiftToEnd( j );
        pvt.setSize(pvt.size()-1);
        j--;
      }
    }
        
    b.m = m = ks;           // reset number of rows
    pvt.setSize( ks );
  }
    
  /** QR transformation for a least squares problem <br/>
   *            A x = b <br/>
   * implicitly both A and b are transformed. pvt.size() is the psuedo-rank of A.
   *  
   * @param b    PaceMatrix b
   * @param pvt  pivoting vector
   * @param k0   the first k0 columns of A (indexed by pvt) are pre-chosen. 
   *            (But subject to rank examination.) 
   * 
   *            For example, the constant term may be reserved, in which
   *            case k0 = 1.
   **/
  public void  lsqrSelection( PaceMatrix b, IntVector pvt, int k0 ) {
    int numObs = m;         // number of instances
    int numXs = pvt.size();

    lsqr( b, pvt, k0 );

    if( numXs > 200 || numXs > numObs ) { // too many columns.  
      forward(b, pvt, k0);
    }
    backward(b, pvt, pvt.size(), k0);
  }
    
  /** 
   * Sets all diagonal elements to be positive (or nonnegative) without
   * changing the least squares solution 
   * @param Y the response
   * @param pvt the pivoted column index
   */
  public void positiveDiagonal( PaceMatrix Y, IntVector pvt ) {
     
    int [] p = pvt.getArray();
    for( int i = 0; i < pvt.size(); i++ ) {
      if( A[i][p[i]] < 0.0 ) {
        for( int j = i; j < pvt.size(); j++ ) 
          A[i][p[j]] = - A[i][p[j]];
        Y.A[i][0] = - Y.A[i][0];
      }
    }
  }

  /** Stepwise least squares QR-decomposition of the problem
   *              A x = b
   @param b    PaceMatrix b
   @param pvt  pivoting vector
   @param ks   number of transformed columns
   @param j    pvt[j], the column to adjoin or delete
   @param adjoin   to adjoin if true; otherwise, to delete */
  public void  steplsqr( PaceMatrix b, IntVector pvt, int ks, int j, 
                         boolean adjoin ) {
    final int kp = pvt.size(); // number of columns under consideration
    int [] p = pvt.getArray();
        
    if( adjoin ) {     // adjoining 
      int pj = p[j];
      pvt.swap( ks, j );
      double dq[] = h1( pj, ks );
      int pk;
      for( int k = ks+1; k < kp; k++ ){
        pk = p[k];
        h2( pj, ks, dq[1], this, pk);
      }
      h2( pj, ks, dq[1], b, 0 ); // for matrix. ???
      A[ks][pj] = dq[0];
      for( int k = ks+1; k < m; k++ )
        A[k][pj] = 0;
    }
    else {          // removing 
      int pj = p[j];
      for( int i = j; i < ks-1; i++ ) 
        p[i] = p[i+1];
      p[ks-1] = pj;
      double [] cs;
      for( int i = j; i < ks-1; i++ ){
        cs = g1( A[i][p[i]], A[i+1][p[i]] );
        for( int l = i; l < kp; l++ ) 
          g2( cs, i, i+1, p[l] );
        for( int l = 0; l < b.n; l++ )
          b.g2( cs, i, i+1, l );
      }
    }
  }
    
  /** Solves upper-triangular equation <br/>
   *    R x = b <br/>
   *  On output, the solution is stored in b
   *  @param b the response
   *  @param pvt the pivoting vector
   *  @param kp the number of the first columns involved 
   */
  public void  rsolve( PaceMatrix b, IntVector pvt, int kp) {
    if(kp == 0) b.m = 0;
    int i, j, k;
    int [] p = pvt.getArray();
    double s;
    double [][] ba = b.getArray();
    for( k = 0; k < b.n; k++ ) {
      ba[kp-1][k] /= A[kp-1][p[kp-1]];
      for( i = kp - 2; i >= 0; i-- ){
        s = 0;
        for( j = i + 1; j < kp; j++ )
          s += A[i][p[j]] * ba[j][k];
        ba[i][k] -= s;
        ba[i][k] /= A[i][p[i]];
      }
    } 
    b.m = kp;
  }
    
  /** Returns a new matrix which binds two matrices together with rows. 
   *  @param b  the second matrix
   *  @return the combined matrix
   */
  public PaceMatrix  rbind( PaceMatrix b ){
    if( n != b.n ) 
      throw new IllegalArgumentException("unequal numbers of rows.");
    PaceMatrix c = new PaceMatrix( m + b.m, n );
    c.setMatrix( 0, m - 1, 0, n - 1, this );
    c.setMatrix( m, m + b.m - 1, 0, n - 1, b );
    return c;
  }

  /** Returns a new matrix which binds two matrices with columns.
   *  @param b the second matrix 
   *  @return the combined matrix
   */
  public PaceMatrix  cbind( PaceMatrix b ) {
    if( m != b.m ) 
      throw new IllegalArgumentException("unequal numbers of rows: " + 
                                         m + " and " + b.m);
    PaceMatrix c = new PaceMatrix(m, n + b.n);
    c.setMatrix( 0, m - 1, 0, n - 1, this );
    c.setMatrix( 0, m - 1, n, n + b.n - 1, b );
    return c;
  }

  /** Solves the nonnegative linear squares problem. That is, <p>
   *   <center>   min || A x - b||, subject to x >= 0.  </center> <p>
   * 
   *  For algorithm, refer to P161, Chapter 23 of C. L. Lawson and
   *  R. J. Hanson (1974).  "Solving Least Squares
   *  Problems". Prentice-Hall.
   *    @param b the response
   *  @param pvt vector storing pivoting column indices
   *    @return solution */
  public DoubleVector nnls( PaceMatrix b, IntVector pvt ) {
    int j, t, counter = 0, jm = -1, n = pvt.size();
    double ma, max, alpha, wj;
    int [] p = pvt.getArray();
    DoubleVector x = new DoubleVector( n );
    double [] xA = x.getArray();
    PaceMatrix z = new PaceMatrix(n, 1);
    PaceMatrix bt;
        
    // step 1 
    int kp = 0; // #variables in the positive set P
    while ( true ) {         // step 2 
      if( ++counter > 3*n )  // should never happen
        throw new RuntimeException("Does not converge");
      t = -1;
      max = 0.0;
      bt = new PaceMatrix( b.transpose() );
      for( j = kp; j <= n-1; j++ ) {   // W = A' (b - A x) 
        wj = bt.times( 0, kp, m-1, this, p[j] );
        if( wj > max ) {        // step 4
          max = wj;
          t = j;
        }
      }
            
      // step 3 
      if ( t == -1) break; // optimum achieved 
            
      // step 5 
      pvt.swap( kp, t );       // move variable from set Z to set P
      kp++;
      xA[kp-1] = 0;
      steplsqr( b, pvt, kp-1, kp-1, true );
      // step 6
      ma = 0;
      while ( ma < 1.5 ) {
        for( j = 0; j <= kp-1; j++ ) z.A[j][0] = b.A[j][0];
        rsolve(z, pvt, kp); 
        ma = 2; jm = -1;
        for( j = 0; j <= kp-1; j++ ) {  // step 7, 8 and 9
          if( z.A[j][0] <= 0.0 ) { // alpha always between 0 and 1
            alpha = xA[j] / ( xA[j] - z.A[j][0] ); 
            if( alpha < ma ) {
              ma = alpha; jm = j;
            }
          }
        }
        if( ma > 1.5 ) 
          for( j = 0; j <= kp-1; j++ ) xA[j] = z.A[j][0];  // step 7 
        else { 
          for( j = kp-1; j >= 0; j-- ) { // step 10
            // Modified to avoid round-off error (which seemingly 
            // can cause infinite loop).
            if( j == jm ) { // step 11 
              xA[j] = 0.0;
              steplsqr( b, pvt, kp, j, false );
              kp--;  // move variable from set P to set Z
            }
            else xA[j] += ma * ( z.A[j][0] - xA[j] );
          }
        }
      }
    }
    x.setSize(kp);
    pvt.setSize(kp);
    return x;
  }

  /** Solves the nonnegative least squares problem with equality
   *    constraint. That is, <p>
   *  <center> min ||A x - b||, subject to x >= 0 and c x = d. </center> <p>
   *
   *  @param b the response
   *  @param c coeficients of equality constraints
   *  @param d constants of equality constraints
   *  @param pvt vector storing pivoting column indices
   *  @return the solution
   */
  public DoubleVector nnlse( PaceMatrix b, PaceMatrix c, PaceMatrix d, 
                             IntVector pvt ) {
    double eps = 1e-10 * Math.max( c.maxAbs(), d.maxAbs() ) /
    Math.max( maxAbs(), b.maxAbs() );
        
    PaceMatrix e = c.rbind( new PaceMatrix( times(eps) ) );
    PaceMatrix f = d.rbind( new PaceMatrix( b.times(eps) ) );

    return e.nnls( f, pvt );
  }

  /** Solves the nonnegative least squares problem with equality
   *    constraint. That is, <p>
   *  <center> min ||A x - b||,  subject to x >= 0 and || x || = 1. </center>
   *  <p>
   *  @param b the response 
   *  @param pvt vector storing pivoting column indices
   *  @return the solution
   */
  public DoubleVector nnlse1( PaceMatrix b, IntVector pvt ) {
    PaceMatrix c = new PaceMatrix( 1, n, 1 );
    PaceMatrix d = new PaceMatrix( 1, b.n, 1 );
        
    return nnlse(b, c, d, pvt);
  }

  /** Generate matrix with standard-normally distributed random elements
      @param m    Number of rows.
      @param n    Number of colums.
      @return An m-by-n matrix with random elements.  */
  public static Matrix randomNormal( int m, int n ) {
    Random random = new Random();
     
    Matrix A = new Matrix(m,n);
    double[][] X = A.getArray();
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        X[i][j] = random.nextGaussian();
      }
    }
    return A;
  }
  
  /**
   * Returns the revision string.
   * 
   * @return            the revision
   */
  public String getRevision() {
    return RevisionUtils.extract("$Revision: 1.6 $");
  }

  /**
   * for testing only
   * 
   * @param args the commandline arguments - ignored
   */
  public static void  main( String args[] ) {
    System.out.println("================================================" + 
                       "===========");
    System.out.println("To test the pace estimators of linear model\n" + 
                       "coefficients.\n");

    double sd = 2;     // standard deviation of the random error term
    int n = 200;       // total number of observations
    double beta0 = 100;   // intercept
    int k1 = 20;       // number of coefficients of the first cluster
    double beta1 = 0;  // coefficient value of the first cluster
    int k2 = 20;      // number of coefficients of the second cluster
    double beta2 = 5; // coefficient value of the second cluster 
    int k = 1 + k1 + k2;

    DoubleVector beta = new DoubleVector( 1 + k1 + k2 );
    beta.set( 0, beta0 );
    beta.set( 1, k1, beta1 );
    beta.set( k1+1, k1+k2, beta2 );

    System.out.println("The data set contains " + n + 
                       " observations plus " + (k1 + k2) + 
                       " variables.\n\nThe coefficients of the true model"
                       + " are:\n\n" + beta );
        
    System.out.println("\nThe standard deviation of the error term is " + 
                       sd );
        
    System.out.println("===============================================" 
                       + "============");
                
    PaceMatrix X = new PaceMatrix( n, k1+k2+1 );
    X.setMatrix( 0, n-1, 0, 0, 1 );
    X.setMatrix( 0, n-1, 1, k1+k2, random(n, k1+k2) );
        
    PaceMatrix Y = new 
      PaceMatrix( X.times( new PaceMatrix(beta) ).
                  plusEquals( randomNormal(n,1).times(sd) ) );

    IntVector pvt = (IntVector) IntVector.seq(0, k1+k2);

    /*System.out.println( "The OLS estimate (by jama.Matrix.solve()) is:\n\n" + 
      (new PaceMatrix(X.solve(Y))).getColumn(0) );*/
        
    X.lsqrSelection( Y, pvt, 1 );
    X.positiveDiagonal( Y, pvt );

    PaceMatrix sol = (PaceMatrix) Y.clone();
    X.rsolve( sol, pvt, pvt.size() );
    DoubleVector betaHat = sol.getColumn(0).unpivoting( pvt, k );
    System.out.println( "\nThe OLS estimate (through lsqr()) is: \n\n" + 
                        betaHat );

    System.out.println( "\nQuadratic loss of the OLS estimate (||X b - X bHat||^2) = " + 
                        ( new PaceMatrix( X.times( new 
                          PaceMatrix(beta.minus(betaHat)) )))
                        .getColumn(0).sum2() );

    System.out.println("=============================================" + 
                       "==============");
    System.out.println("             *** Pace estimation *** \n");
    DoubleVector r = Y.getColumn( pvt.size(), n-1, 0);
    double sde = Math.sqrt(r.sum2() / r.size());
        
    System.out.println( "Estimated standard deviation = " + sde );

    DoubleVector aHat = Y.getColumn( 0, pvt.size()-1, 0).times( 1./sde );
    System.out.println("\naHat = \n" + aHat );
        
    System.out.println("\n========= Based on chi-square mixture ============");

    ChisqMixture d2 = new ChisqMixture();
    int method = MixtureDistribution.NNMMethod;
    DoubleVector AHat = aHat.square();
    d2.fit( AHat, method ); 
    System.out.println( "\nEstimated mixing distribution is:\n" + d2 );
        
    DoubleVector ATilde = d2.pace2( AHat );
    DoubleVector aTilde = ATilde.sqrt().times(aHat.sign());
    PaceMatrix YTilde = new 
      PaceMatrix((new PaceMatrix(aTilde)).times( sde ));
    X.rsolve( YTilde, pvt, pvt.size() );
    DoubleVector betaTilde = 
    YTilde.getColumn(0).unpivoting( pvt, k );
    System.out.println( "\nThe pace2 estimate of coefficients = \n" + 
                        betaTilde );
    System.out.println( "Quadratic loss = " + 
                        ( new PaceMatrix( X.times( new 
                          PaceMatrix(beta.minus(betaTilde)) )))
                        .getColumn(0).sum2() );
        
    ATilde = d2.pace4( AHat );
    aTilde = ATilde.sqrt().times(aHat.sign());
    YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times( sde ));
    X.rsolve( YTilde, pvt, pvt.size() );
    betaTilde = YTilde.getColumn(0).unpivoting( pvt, k );
    System.out.println( "\nThe pace4 estimate of coefficients = \n" + 
                        betaTilde );
    System.out.println( "Quadratic loss = " + 
                        ( new PaceMatrix( X.times( new 
                          PaceMatrix(beta.minus(betaTilde)) )))
                        .getColumn(0).sum2() );
        
    ATilde = d2.pace6( AHat );
    aTilde = ATilde.sqrt().times(aHat.sign());
    YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times( sde ));
    X.rsolve( YTilde, pvt, pvt.size() );
    betaTilde = YTilde.getColumn(0).unpivoting( pvt, k );
    System.out.println( "\nThe pace6 estimate of coefficients = \n" + 
                        betaTilde );
    System.out.println( "Quadratic loss = " + 
                        ( new PaceMatrix( X.times( new 
                          PaceMatrix(beta.minus(betaTilde)) )))
                        .getColumn(0).sum2() );
        
    System.out.println("\n========= Based on normal mixture ============");
        
    NormalMixture d = new NormalMixture();
    d.fit( aHat, method ); 
    System.out.println( "\nEstimated mixing distribution is:\n" + d );
        
    aTilde = d.nestedEstimate( aHat );
    YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times( sde ));
    X.rsolve( YTilde, pvt, pvt.size() );
    betaTilde = YTilde.getColumn(0).unpivoting( pvt, k );
    System.out.println( "The nested estimate of coefficients = \n" + 
                        betaTilde );
    System.out.println( "Quadratic loss = " + 
                        ( new PaceMatrix( X.times( new 
                          PaceMatrix(beta.minus(betaTilde)) )))
                        .getColumn(0).sum2() );
        
        
    aTilde = d.subsetEstimate( aHat );
    YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times( sde ));
    X.rsolve( YTilde, pvt, pvt.size() );
    betaTilde = 
    YTilde.getColumn(0).unpivoting( pvt, k );
    System.out.println( "\nThe subset estimate of coefficients = \n" + 
                        betaTilde );
    System.out.println( "Quadratic loss = " + 
                        ( new PaceMatrix( X.times( new 
                          PaceMatrix(beta.minus(betaTilde)) )))
                        .getColumn(0).sum2() );
        
    aTilde = d.empiricalBayesEstimate( aHat );
    YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times( sde ));
    X.rsolve( YTilde, pvt, pvt.size() );
    betaTilde = YTilde.getColumn(0).unpivoting( pvt, k );
    System.out.println( "\nThe empirical Bayes estimate of coefficients = \n"+
                        betaTilde );
        
    System.out.println( "Quadratic loss = " + 
                        ( new PaceMatrix( X.times( new 
                          PaceMatrix(beta.minus(betaTilde)) )))
                        .getColumn(0).sum2() );
        
  }
}