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,