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;  }  /**   * 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() );	  }}