pacematrix.java

为了下东西 随便发了个 datamining 的源代码

  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 */
  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 */
  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 
   **/
  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.
   *@return      void
   */
  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.  */
  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
   *            A x = b
   *  
   *@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.
   *@return     void; implicitly both A and b are transformed. pvt.size() is 
   *            is the psuedo-rank of A.
   **/
  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
   *            A x = b
   *  
   *@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.
   *@return     void; implicitly both A and b are transformed. pvt.size() is 
   *            is the psuedo-rank of A.
   **/
  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 );