content: /*
   *    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.
   */
  
  /*
   *    Optimization.java
   *    Copyright (C) 2003 University of Waikato, Hamilton, New Zealand
   *
   */
  
  package weka.core;
  
  import weka.core.TechnicalInformation.Field;
  import weka.core.TechnicalInformation.Type;
  
  /**
   * Implementation of Active-sets method with BFGS update to solve optimization
   * problem with only bounds constraints in multi-dimensions.  In this
   * implementation we consider both the lower and higher bound constraints. <p/>
   *
   * Here is the sketch of our searching strategy, and the detailed description
   * of the algorithm can be found in the Appendix of Xin Xu's MSc thesis:<p/>
   * Initialize everything, incl. initial value, direction, etc.<p/>
   * LOOP (main algorithm):<br/>
   * 
   * 1. Perform the line search using the directions for free variables<br/>
   * 1.1  Check all the bounds that are not "active" (i.e. binding variables)
   *      and compute the feasible step length to the bound for each of them<br/>
   * 1.2  Pick up the least feasible step length, say \alpha, and set it as 
   *      the upper bound of the current step length, i.e.
   *      0&lt;\lambda&lt;=\alpha<br/>
   * 1.3  Search for any possible step length&lt;=\alpha that can result the 
   *      "sufficient function decrease" (\alpha condition) AND "positive
   *      definite inverse Hessian" (\beta condition), if possible, using
   *      SAFEGUARDED polynomial interpolation.  This step length is "safe" and
   *      thus is used to compute the next value of the free variables .<br/>
   * 1.4  Fix the variable(s) that are newly bound to its constraint(s).<p/>     
   *
   * 2. Check whether there is convergence of all variables or their gradients.
   *    If there is, check the possibilities to release any current bindings of
   *    the fixed variables to their bounds based on the "reliable" second-order
   *    Lagarange multipliers if available.  If it's available and negative for
   *    one variable, then release it.  If not available, use first-order
   *    Lagarange multiplier to test release.  If there is any released
   *    variables, STOP the loop.  Otherwise update the inverse of Hessian matrix
   *    and gradient for the newly released variables and CONTINUE LOOP.<p/>
   *
   * 3. Use BFGS formula to update the inverse of Hessian matrix.  Note the 
   *    already-fixed variables must have zeros in the corresponding entries
   *    in the inverse Hessian.<p/>  
   *
   * 4. Compute the new (newton) search direction d=H^{-1}*g, where H^{-1} is the 
   *    inverse Hessian and g is the Jacobian.  Note that again, the already-
   *    fixed variables will have zero direction.<p/>
   *
   * ENDLOOP<p/>
   *
   * A typical usage of this class is to create your own subclass of this class
   * and provide the objective function and gradients as follows:<p/>
   * <pre>
   * class MyOpt extends Optimization {
   *   // Provide the objective function
   *   protected double objectiveFunction(double[] x) {
   *       // How to calculate your objective function...
   *       // ...
   *   }
   *
   *   // Provide the first derivatives
   *   protected double[] evaluateGradient(double[] x) {
   *       // How to calculate the gradient of the objective function...
   *       // ...
   *   }
   *
   *   // If possible, provide the index^{th} row of the Hessian matrix
   *   protected double[] evaluateHessian(double[] x, int index) {
   *      // How to calculate the index^th variable's second derivative
   *      // ... 
   *   }
   * }
   * </pre>
   *
   * When it's the time to use it, in some routine(s) of other class...
   * <pre>
   * MyOpt opt = new MyOpt();
   * 
   * // Set up initial variable values and bound constraints
   * double[] x = new double[numVariables];
   * // Lower and upper bounds: 1st row is lower bounds, 2nd is upper
   * double[] constraints = new double[2][numVariables];
   * ...
   *
   * // Find the minimum, 200 iterations as default
   * x = opt.findArgmin(x, constraints); 
   * while(x == null){  // 200 iterations are not enough
   *    x = opt.getVarbValues();  // Try another 200 iterations
   *    x = opt.findArgmin(x, constraints);
   * }
   *
   * // The minimal function value
   * double minFunction = opt.getMinFunction();
   * ...
   * </pre>
   * 
   * It is recommended that Hessian values be provided so that the second-order
   * Lagrangian multiplier estimate can be calcluated.  However, if it is not
   * provided, there is no need to override the <code>evaluateHessian()</code>
   * function.<p/>
   *
   * REFERENCES (see also the <code>getTechnicalInformation()</code> method):<br/>
   * The whole model algorithm is adapted from Chapter 5 and other related
   * chapters in Gill, Murray and Wright(1981) "Practical Optimization", Academic
   * Press.  and Gill and Murray(1976) "Minimization Subject to Bounds on the
   * Variables", NPL Report NAC72, while Chong and Zak(1996) "An Introduction to
   * Optimization", John Wiley &amp; Sons, Inc. provides us a brief but helpful
   * introduction to the method. <p/>
   *
   * Dennis and Schnabel(1983) "Numerical Methods for Unconstrained Optimization
   * and Nonlinear Equations", Prentice-Hall Inc. and Press et al.(1992) "Numeric
   * Recipe in C", Second Edition, Cambridge University Press. are consulted for
   * the polynomial interpolation used in the line search implementation.  <p/>
   *
   * The Hessian modification in BFGS update uses Cholesky factorization and two
   * rank-one modifications:<br/>
   * Bk+1 = Bk + (Gk*Gk')/(Gk'Dk) + (dGk*(dGk)'))/[alpha*(dGk)'*Dk]. <br/>
   * where Gk is the gradient vector, Dk is the direction vector and alpha is the
   * step rate. <br/>
   * This method is due to Gill, Golub, Murray and Saunders(1974) ``Methods for
   * Modifying Matrix Factorizations'', Mathematics of Computation, Vol.28,
   * No.126, pp 505-535. <p/>
   *
   * @author Xin Xu (xx5@cs.waikato.ac.nz)
   * @version $Revision: 5953 $ 
   * @see #getTechnicalInformation()
   */
  public abstract class Optimization
      implements TechnicalInformationHandler, RevisionHandler {
      
      protected double m_ALF = 1.0e-4;
  
      protected double m_BETA = 0.9;    
  
      protected double m_TOLX = 1.0e-6;
     
      protected double m_STPMX = 100.0;
      
      protected int m_MAXITS = 200;
      
      protected static boolean m_Debug = false;
      
      /** function value */
      protected double m_f;    
   
      /** G'*p */
      private double m_Slope;
      
      /** Test if zero step in lnsrch */
      private boolean m_IsZeroStep = false;
      
      /** Used when iteration overflow occurs */
      private double[] m_X;
      
      /** Compute machine precision */
      protected static double m_Epsilon, m_Zero; 
      static {
  	m_Epsilon=1.0;
  	while(1.0+m_Epsilon > 1.0){
  	    m_Epsilon /= 2.0;	    
  	}
  	m_Epsilon *= 2.0;
  	m_Zero = Math.sqrt(m_Epsilon);
  	if (m_Debug)
  	    System.err.print("Machine precision is "+m_Epsilon+
  			     " and zero set to "+m_Zero);
      }
      
      /**
       * Returns an instance of a TechnicalInformation object, containing 
       * detailed information about the technical background of this class,
       * e.g., paper reference or book this class is based on.
       * 
       * @return the technical information about this class
       */
      public TechnicalInformation getTechnicalInformation() {
        TechnicalInformation 	result;
        TechnicalInformation 	additional;
        
        result = new TechnicalInformation(Type.MASTERSTHESIS);
        result.setValue(Field.AUTHOR, "Xin Xu");
        result.setValue(Field.YEAR, "2003");
        result.setValue(Field.TITLE, "Statistical learning in multiple instance problem");
        result.setValue(Field.SCHOOL, "University of Waikato");
        result.setValue(Field.ADDRESS, "Hamilton, NZ");
        result.setValue(Field.NOTE, "0657.594");
  
        additional = result.add(Type.BOOK);
        additional.setValue(Field.AUTHOR, "P. E. Gill and W. Murray and M. H. Wright");
        additional.setValue(Field.YEAR, "1981");
        additional.setValue(Field.TITLE, "Practical Optimization");
        additional.setValue(Field.PUBLISHER, "Academic Press");
        additional.setValue(Field.ADDRESS, "London and New York");
        
        additional = result.add(Type.TECHREPORT);
        additional.setValue(Field.AUTHOR, "P. E. Gill and W. Murray");
        additional.setValue(Field.YEAR, "1976");
        additional.setValue(Field.TITLE, "Minimization subject to bounds on the variables");
        additional.setValue(Field.INSTITUTION, "National Physical Laboratory");
        additional.setValue(Field.NUMBER, "NAC 72");
        
        additional = result.add(Type.BOOK);
        additional.setValue(Field.AUTHOR, "E. K. P. Chong and S. H. Zak");
        additional.setValue(Field.YEAR, "1996");
        additional.setValue(Field.TITLE, "An Introduction to Optimization");
        additional.setValue(Field.PUBLISHER, "John Wiley and Sons");
        additional.setValue(Field.ADDRESS, "New York");
        
        additional = result.add(Type.BOOK);
        additional.setValue(Field.AUTHOR, "J. E. Dennis and R. B. Schnabel");
        additional.setValue(Field.YEAR, "1983");
        additional.setValue(Field.TITLE, "Numerical Methods for Unconstrained Optimization and Nonlinear Equations");
        additional.setValue(Field.PUBLISHER, "Prentice-Hall");
        
        additional = result.add(Type.BOOK);
        additional.setValue(Field.AUTHOR, "W. H. Press and B. P. Flannery and S. A. Teukolsky and W. T. Vetterling");
        additional.setValue(Field.YEAR, "1992");
        additional.setValue(Field.TITLE, "Numerical Recipes in C");
        additional.setValue(Field.PUBLISHER, "Cambridge University Press");
        additional.setValue(Field.EDITION, "Second");
        
        additional = result.add(Type.ARTICLE);
        additional.setValue(Field.AUTHOR, "P. E. Gill and G. H. Golub and W. Murray and M. A. Saunders");
        additional.setValue(Field.YEAR, "1974");
        additional.setValue(Field.TITLE, "Methods for modifying matrix factorizations");
        additional.setValue(Field.JOURNAL, "Mathematics of Computation");
        additional.setValue(Field.VOLUME, "28");
        additional.setValue(Field.NUMBER, "126");
        additional.setValue(Field.PAGES, "505-535");
        
        return result;
      }
      
      /**
       * Subclass should implement this procedure to evaluate objective
       * function to be minimized
       * 
       * @param x the variable values
       * @return the objective function value
       * @throws Exception if something goes wrong
       */
      protected abstract double objectiveFunction(double[] x) throws Exception;
  
      /**
       * Subclass should implement this procedure to evaluate gradient
       * of the objective function
       * 
       * @param x the variable values
       * @return the gradient vector
       * @throws Exception if something goes wrong
       */
      protected abstract double[] evaluateGradient(double[] x) throws Exception;
  
      /**
       * Subclass is recommended to override this procedure to evaluate second-order
       * gradient of the objective function.  If it's not provided, it returns
       * null.
       *
       * @param x the variables
       * @param index the row index in the Hessian matrix
       * @return one row (the row #index) of the Hessian matrix, null as default
       * @throws Exception if something goes wrong
       */
      protected double[] evaluateHessian(double[] x, int index) throws Exception{
  	return null;
      }
  
      /**
       * Get the minimal function value
       *
       * @return minimal function value found
       */
      public double getMinFunction() {
        return m_f;
      }
  
      /**
       * Set the maximal number of iterations in searching (Default 200)
       *
       * @param it the maximal number of iterations
       */
      public void setMaxIteration(int it) {
        m_MAXITS=it;
      }
        
      /**
       * Set whether in debug mode
       *
       * @param db use debug or not
       */
      public void setDebug(boolean db) {
        m_Debug = db;
      }
      
      /**
       * Get the variable values.  Only needed when iterations exceeds 
       * the max threshold.
       *
       * @return the current variable values
       */
      public double[] getVarbValues() {
        return m_X;
      }
      
      /**
       * Find a new point x in the direction p from a point xold at which the
       * value of the function has decreased sufficiently, the positive 
       * definiteness of B matrix (approximation of the inverse of the Hessian)
       * is preserved and no bound constraints are violated.  Details see "Numerical 
       * Methods for Unconstrained Optimization and Nonlinear Equations".
       * "Numeric Recipes in C" was also consulted.
       *
       * @param xold old x value 
       * @param gradient gradient at that point
       * @param direct direction vector
       * @param stpmax maximum step length
       * @param isFixed indicating whether a variable has been fixed
       * @param nwsBounds non-working set bounds.  Means these variables are free and
       *                  subject to the bound constraints in this step
       * @param wsBdsIndx index of variables that has working-set bounds.  Means
       *                  these variables are already fixed and no longer subject to
       *                  the constraints
       * @return new value along direction p from xold, null if no step was taken
       * @throws Exception if an error occurs
       */
      public double[] lnsrch(double[] xold, double[] gradient, 
  			   double[] direct, double stpmax,
  			   boolean[] isFixed, double[][] nwsBounds,
  			   DynamicIntArray wsBdsIndx)
  	throws Exception {
  	
  	int i, k,len=xold.length, 
  	    fixedOne=-1; // idx of variable to be fixed
  	double alam, alamin; // lambda to be found, and its lower bound
  	
  	// For convergence and bound test
  	double temp,test,alpha=Double.POSITIVE_INFINITY,fold=m_f,sum; 
  	
  	// For cubic interpolation
  	double a,alam2=0,b,disc=0,maxalam=1.0,rhs1,rhs2,tmplam;
  	
  	double[] x = new double[len]; // New variable values
  	
  	// Scale the step 
  	for (sum=0.0,i=0;i<len;i++){
  	    if(!isFixed[i]) // For fixed variables, direction = 0
  		sum += direct[i]*direct[i];
  	}	
  	sum = Math.sqrt(sum);
  	
  	if (m_Debug)
  	    System.err.println("fold:  "+Utils.doubleToString(fold,10,7)+"\n"+
  			       "sum:  "+Utils.doubleToString(sum,10,7)+"\n"+
  			       "stpmax:  "+Utils.doubleToString(stpmax,10,7));
  	if (sum > stpmax){
  	    for (i=0;i<len;i++)
  		if(!isFixed[i])
  		    direct[i] *= stpmax/sum;		
  	}
  	else
  	    maxalam = stpmax/sum;
  	
  	// Compute initial rate of decrease, g'*d 
  	m_Slope=0.0;
  	for (i=0;i<len;i++){
  	    x[i] = xold[i];
  	    if(!isFixed[i])
  		m_Slope += gradient[i]*direct[i];
  	}
  	
  	if (m_Debug)
  	    System.err.print("slope:  " + Utils.doubleToString(m_Slope,10,7)+ "\n");
  	
  	// Slope too small
  	if(Math.abs(m_Slope)<=m_Zero){
  	    if (m_Debug)
  		System.err.println("Gradient and direction orthogonal -- "+
  				   "Min. found with current fixed variables"+
  				   " (or all variables fixed). Try to release"+
  				   " some variables now.");
  	    return x;
  	}
  	
  	// Err: slope > 0
  	if(m_Slope > m_Zero){
  	    if(m_Debug)
  		for(int h=0; h<x.length; h++)
  		    System.err.println(h+": isFixed="+isFixed[h]+", x="+
  				       x[h]+", grad="+gradient[h]+", direct="+
  				       direct[h]);
  	    throw new Exception("g'*p positive! -- Try to debug from here: line 327.");
  	}
  	
  	// Compute LAMBDAmin and upper bound of lambda--alpha
  	test=0.0;
  	for(i=0;i<len;i++){	    
  	    if(!isFixed[i]){// No need for fixed variables
  		temp=Math.abs(direct[i])/Math.max(Math.abs(x[i]),1.0);
  		if (temp > test) test=temp;
  	    }
  	}
  	
  	if(test>m_Zero) // Not converge
  	    alamin = m_TOLX/test;
  	else{
  	    if (m_Debug)
  		System.err.println("Zero directions for all free variables -- "+
  				   "Min. found with current fixed variables"+
  				   " (or all variables fixed). Try to release"+
  				   " some variables now.");
  	    return x;
  	}
  		
  	// Check whether any non-working-set bounds are "binding"
  	for(i=0;i<len;i++){
  	    if(!isFixed[i]){// No need for fixed variables
  		double alpi;
  		if((direct[i]<-m_Epsilon) && !Double.isNaN(nwsBounds[0][i])){//Not feasible
  		    alpi = (nwsBounds[0][i]-xold[i])/direct[i];
  		    if(alpi <= m_Zero){ // Zero
  			if (m_Debug)
  			    System.err.println("Fix variable "+i+
  					       " to lower bound "+ nwsBounds[0][i]+
  					       " from value "+ xold[i]);
  			x[i] = nwsBounds[0][i];
  			isFixed[i]=true; // Fix this variable
  			alpha = 0.0;
  			nwsBounds[0][i]=Double.NaN; //Add cons. to working set
  			wsBdsIndx.addElement(i);
  		    }
  		    else if(alpha > alpi){ // Fix one variable in one iteration
  			alpha = alpi;
  			fixedOne = i;
  		    }			
  		}
  		else if((direct[i]>m_Epsilon) && !Double.isNaN(nwsBounds[1][i])){//Not feasible
  		    alpi = (nwsBounds[1][i]-xold[i])/direct[i];
  		    if(alpi <= m_Zero){ // Zero
  			if (m_Debug)
  			    System.err.println("Fix variable "+i+
  					       " to upper bound "+ nwsBounds[1][i]+
  					       " from value "+ xold[i]);
  			x[i] = nwsBounds[1][i];
  			isFixed[i]=true; // Fix this variable
  			alpha = 0.0;
  			nwsBounds[1][i]=Double.NaN; //Add cons. to working set
  			wsBdsIndx.addElement(i);
  		    }
  		    else if(alpha > alpi){
  			alpha = alpi;
  			fixedOne = i;
  		    }			
  		}				
  	    }
  	}	
  	
  	if (m_Debug){
  	    System.err.println("alamin: " + Utils.doubleToString(alamin,10,7));
  	    System.err.println("alpha: " + Utils.doubleToString(alpha,10,7));
  	}
  	
  	if(alpha <= m_Zero){ // Zero	   
  	    m_IsZeroStep = true;
  	    if (m_Debug)
  		System.err.println("Alpha too small, try again");
  	    return x;
  	}
  	
  	alam = alpha; // Always try full feasible newton step 
  	if(alam > 1.0)
  	    alam = 1.0;
  	
  	// Iteration of one newton step, if necessary, backtracking is done
  	double initF=fold, // Initial function value
  	    hi=alam, lo=alam, newSlope=0, fhi=m_f, flo=m_f;// Variables used for beta condition
  	double[] newGrad;  // Gradient on the new variable values
  	
  	kloop:
  	for (k=0;;k++) {
  	    if(m_Debug)
  		System.err.println("\nLine search iteration: " + k);
  	    
  	    for (i=0;i<len;i++){
  		if(!isFixed[i]){
  		    x[i] = xold[i]+alam*direct[i];  // Compute xnew
  		    if(!Double.isNaN(nwsBounds[0][i]) && (x[i]<nwsBounds[0][i])){    
  			x[i] = nwsBounds[0][i]; //Rounding error	
  		    }
  		    else if(!Double.isNaN(nwsBounds[1][i]) && (x[i]>nwsBounds[1][i])){		
  			x[i] = nwsBounds[1][i]; //Rounding error	
  		    }
  		}
  	    }
  	    
  	    m_f = objectiveFunction(x);    // Compute fnew
  	    if(Double.isNaN(m_f))
  		throw new Exception("Objective function value is NaN!");
  	
  	    while(Double.isInfinite(m_f)){ // Avoid infinity
  		if(m_Debug)
  		    System.err.println("Too large m_f.  Shrink step by half.");
  		alam *= 0.5; // Shrink by half
  		if(alam <= m_Epsilon){
  		    if(m_Debug)
  			System.err.println("Wrong starting points, change them!");
  		    return x;
  		}
  		
  		for (i=0;i<len;i++)
  		    if(!isFixed[i])
  			x[i] = xold[i]+alam*direct[i]; 
  		
  		m_f = objectiveFunction(x); 
  		if(Double.isNaN(m_f))
  		    throw new Exception("Objective function value is NaN!");
  	
  		initF = Double.POSITIVE_INFINITY;
  	    }
  	    
  	    if(m_Debug) {
  		System.err.println("obj. function: " + 
  				   Utils.doubleToString(m_f, 10, 7));
  		System.err.println("threshold: " + 
  				   Utils.doubleToString(fold+m_ALF*alam*m_Slope,10,7));
  	    }
  	    
  	    if(m_f<=fold+m_ALF*alam*m_Slope){// Alpha condition: sufficient function decrease
  		if(m_Debug)		
  		    System.err.println("Sufficient function decrease (alpha condition): "); 
  		newGrad = evaluateGradient(x);
  		for(newSlope=0.0,i=0; i<len; i++)
  		    if(!isFixed[i])
  			newSlope += newGrad[i]*direct[i];
  
  		if(newSlope >= m_BETA*m_Slope){ // Beta condition: ensure pos. defnty.	
  		    if(m_Debug)		
  			System.err.println("Increasing derivatives (beta condition): "); 	
  
  		    if((fixedOne!=-1) && (alam>=alpha)){ // Has bounds and over
  			if(direct[fixedOne] > 0){
  			    x[fixedOne] = nwsBounds[1][fixedOne]; // Avoid rounding error
  			    nwsBounds[1][fixedOne]=Double.NaN; //Add cons. to working set
  			}
  			else{
  			    x[fixedOne] = nwsBounds[0][fixedOne]; // Avoid rounding error
  			    nwsBounds[0][fixedOne]=Double.NaN; //Add cons. to working set
  			}
  			
  			if(m_Debug)
  			    System.err.println("Fix variable "
  					       +fixedOne+" to bound "+ x[fixedOne]+
  					       " from value "+ xold[fixedOne]);
  			isFixed[fixedOne]=true; // Fix the variable
  			wsBdsIndx.addElement(fixedOne);
  		    }		
  		    return x;
  		}
  		else if(k==0){ // First time: increase alam 
  		    // Search for the smallest value not complying with alpha condition
  		    double upper = Math.min(alpha,maxalam); 
  		    if(m_Debug)
  			System.err.println("Alpha condition holds, increase alpha... ");
  		    while(!((alam>=upper) || (m_f>fold+m_ALF*alam*m_Slope))){
  			lo = alam;
  			flo = m_f;
  			alam *= 2.0;
  			if(alam>=upper)  // Avoid rounding errors
  			    alam=upper;
  
  			for (i=0;i<len;i++)
  			    if(!isFixed[i])
  				x[i] = xold[i]+alam*direct[i];
  			m_f = objectiveFunction(x);
  			if(Double.isNaN(m_f))
  			    throw new Exception("Objective function value is NaN!");
  			
  			newGrad = evaluateGradient(x);
  			for(newSlope=0.0,i=0; i<len; i++)
  			    if(!isFixed[i])
  				newSlope += newGrad[i]*direct[i];
  			
  			if(newSlope >= m_BETA*m_Slope){
  			    if (m_Debug)		
  				System.err.println("Increasing derivatives (beta condition): \n"+
  						   "newSlope = "+Utils.doubleToString(newSlope,10,7));
  			    
  			    if((fixedOne!=-1) && (alam>=alpha)){ // Has bounds and over
  				if(direct[fixedOne] > 0){
  				    x[fixedOne] = nwsBounds[1][fixedOne]; // Avoid rounding error
  				    nwsBounds[1][fixedOne]=Double.NaN; //Add cons. to working set
  				}
  				else{
  				    x[fixedOne] = nwsBounds[0][fixedOne]; // Avoid rounding error
  				    nwsBounds[0][fixedOne]=Double.NaN; //Add cons. to working set
  				}
  				
  				if(m_Debug)
  				    System.err.println("Fix variable "
  						       +fixedOne+" to bound "+ x[fixedOne]+
  						       " from value "+ xold[fixedOne]);
  				isFixed[fixedOne]=true; // Fix the variable
  				wsBdsIndx.addElement(fixedOne);
  			    }		 				    
  			    return x;
  			}
  		    }
  		    hi = alam;
  		    fhi = m_f;			
  		    break kloop;
  		}
  		else{
  		    if(m_Debug)
  			System.err.println("Alpha condition holds.");
  		    hi = alam2; lo = alam; flo = m_f;
  		    break kloop;
  		}		    
  	    }        
  	    else if (alam < alamin) { // No feasible lambda found       
  		if(initF<fold){ 
  		    alam = Math.min(1.0,alpha);
  		    for (i=0;i<len;i++)
  			if(!isFixed[i])
  			    x[i] = xold[i]+alam*direct[i]; //Still take Alpha
  		    
  		    if (m_Debug)
  			System.err.println("No feasible lambda: still take"+
  					   " alpha="+alam);
  		    
  		    if((fixedOne!=-1) && (alam>=alpha)){ // Has bounds and over
  			if(direct[fixedOne] > 0){
  			    x[fixedOne] = nwsBounds[1][fixedOne]; // Avoid rounding error
  			    nwsBounds[1][fixedOne]=Double.NaN; //Add cons. to working set
  			}
  			else{
  			    x[fixedOne] = nwsBounds[0][fixedOne]; // Avoid rounding error
  			    nwsBounds[0][fixedOne]=Double.NaN; //Add cons. to working set
  			}
  			
  			if(m_Debug)
  			    System.err.println("Fix variable "
  					       +fixedOne+" to bound "+ x[fixedOne]+
  					       " from value "+ xold[fixedOne]);
  			isFixed[fixedOne]=true; // Fix the variable
  			wsBdsIndx.addElement(fixedOne);
  		    }		 		    
  		}
  		else{   // Convergence on delta(x)
  		    for(i=0;i<len;i++) 
  			x[i]=xold[i];
  		    m_f=fold;
  		    if (m_Debug)
  			System.err.println("Cannot find feasible lambda"); 
  		}
  		
  		return x; 
  	    }
  	    else { // Backtracking by polynomial interpolation
  		if(k==0){ // First time backtrack: quadratic interpolation
  		    if(!Double.isInfinite(initF))
  			initF = m_f;		    
  		    // lambda = -g'(0)/(2*g''(0))
  		    tmplam = -0.5*alam*m_Slope/((m_f-fold)/alam-m_Slope);
  		}
  		else {    // Subsequent backtrack: cubic interpolation 
  		    rhs1 = m_f-fold-alam*m_Slope;
  		    rhs2 = fhi-fold-alam2*m_Slope;
  		    a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
  		    b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
  		    if (a == 0.0) tmplam = -m_Slope/(2.0*b);
  		    else {
  			disc=b*b-3.0*a*m_Slope;
  			if (disc < 0.0) disc = 0.0;
  			double numerator = -b+Math.sqrt(disc);
  			if(numerator >= Double.MAX_VALUE){
  			    numerator = Double.MAX_VALUE;
  			    if (m_Debug)
  				System.err.print("-b+sqrt(disc) too large! Set it to MAX_VALUE.");
  			}
  			tmplam=numerator/(3.0*a);
  		    }
  		    if (m_Debug)
  			System.err.print("Cubic interpolation: \n" + 
  					 "a:   " + Utils.doubleToString(a,10,7)+ "\n" +
  					 "b:   " + Utils.doubleToString(b,10,7)+ "\n" +
  					 "disc:   " + Utils.doubleToString(disc,10,7)+ "\n" +
  					 "tmplam:   " + tmplam + "\n" +
  					 "alam:   " + Utils.doubleToString(alam,10,7)+ "\n");	
  		    if (tmplam>0.5*alam)
  			tmplam=0.5*alam;             // lambda <= 0.5*lambda_old
  		}
  	    }
  	    alam2=alam;
  	    fhi=m_f;
  	    alam=Math.max(tmplam,0.1*alam);          // lambda >= 0.1*lambda_old
  	    
  	    if(alam>alpha){
  		throw new Exception("Sth. wrong in lnsrch:"+
  				    "Lambda infeasible!(lambda="+alam+
  				    ", alpha="+alpha+", upper="+tmplam+
  				    "|"+(-alpha*m_Slope/(2.0*((m_f-fold)/alpha-m_Slope)))+
  				    ", m_f="+m_f+", fold="+fold+
  				    ", slope="+m_Slope);
  	    }	    
  	} // Endfor(k=0;;k++)
  	
  	// Quadratic interpolation between lamda values between lo and hi.
  	// If cannot find a value satisfying beta condition, use lo.
  	double ldiff = hi-lo, lincr;
  	if(m_Debug)
  	    System.err.println("Last stage of searching for beta condition (alam between "
  			       +Utils.doubleToString(lo,10,7)+" and "
  			       +Utils.doubleToString(hi,10,7)+")...\n"+
  			       "Quadratic Interpolation(QI):\n"+
  			       "Last newSlope = "+Utils.doubleToString(newSlope, 10, 7));
  	
  	while((newSlope<m_BETA*m_Slope) && (ldiff>=alamin)){
  	    lincr = -0.5*newSlope*ldiff*ldiff/(fhi-flo-newSlope*ldiff);
  
  	    if(m_Debug)
  		System.err.println("fhi = "+fhi+"\n"+
  				   "flo = "+flo+"\n"+
  				   "ldiff = "+ldiff+"\n"+
  				   "lincr (using QI) = "+lincr+"\n");
  	    
  	    if(lincr<0.2*ldiff) lincr=0.2*ldiff;
  	    alam = lo+lincr;
  	    if(alam >= hi){ // We cannot go beyond the bounds, so the best we can try is hi
  	    	alam=hi;
  		lincr=ldiff;
  	    }
  	    for (i=0;i<len;i++)
  		if(!isFixed[i])
  		    x[i] = xold[i]+alam*direct[i];
  	    m_f = objectiveFunction(x);
  	    if(Double.isNaN(m_f))
  		throw new Exception("Objective function value is NaN!");
  	
  	    if(m_f>fold+m_ALF*alam*m_Slope){ 
  		// Alpha condition fails, shrink lambda_upper
  		ldiff = lincr;
  		fhi = m_f;
  	    }	    
  	    else{ // Alpha condition holds	    
  		newGrad = evaluateGradient(x);
  		for(newSlope=0.0,i=0; i<len; i++)
  		    if(!isFixed[i])
  			newSlope += newGrad[i]*direct[i];
  		
  		if(newSlope < m_BETA*m_Slope){
  		    // Beta condition fails, shrink lambda_lower
  		    lo = alam;
  		    ldiff -= lincr;
  		    flo = m_f;
  		}
  	    }
  	} 
  	
  	if(newSlope < m_BETA*m_Slope){ // Cannot satisfy beta condition, take lo
  	    if(m_Debug)
  		System.err.println("Beta condition cannot be satisfied, take alpha condition");
  	    alam=lo;
  	    for (i=0;i<len;i++)
  		if(!isFixed[i])
  		    x[i] = xold[i]+alam*direct[i];
  	    m_f = flo;
  	}
  	else if(m_Debug)
  	    System.err.println("Both alpha and beta conditions are satisfied. alam="
  			       +Utils.doubleToString(alam,10,7));
  	
  	if((fixedOne!=-1) && (alam>=alpha)){ // Has bounds and over
  	    if(direct[fixedOne] > 0){
  		x[fixedOne] = nwsBounds[1][fixedOne]; // Avoid rounding error
  		nwsBounds[1][fixedOne]=Double.NaN; //Add cons. to working set
  	    }
  	    else{
  		x[fixedOne] = nwsBounds[0][fixedOne]; // Avoid rounding error
  		nwsBounds[0][fixedOne]=Double.NaN; //Add cons. to working set
  	    }
  	    
  	    if(m_Debug)
  		System.err.println("Fix variable "
  				   +fixedOne+" to bound "+ x[fixedOne]+
  				   " from value "+ xold[fixedOne]);
  	    isFixed[fixedOne]=true; // Fix the variable
  	    wsBdsIndx.addElement(fixedOne);
  	}
  	
  	return x;
      }
      
      /**
       * Main algorithm.  Descriptions see "Practical Optimization"
       *
       * @param initX initial point of x, assuming no value's on the bound!
       * @param constraints the bound constraints of each variable
       *                    constraints[0] is the lower bounds and 
       *                    constraints[1] is the upper bounds
       * @return the solution of x, null if number of iterations not enough
       * @throws Exception if an error occurs
       */
      public double[] findArgmin(double[] initX, double[][] constraints) 
  	throws Exception{
  	int l = initX.length;
  	
  	// Initially all variables are free, all bounds are constraints of
  	// non-working-set constraints
  	boolean[] isFixed = new boolean[l];
  	double[][] nwsBounds = new double[2][l];
  	// Record indice of fixed variables, simply for efficiency
  	DynamicIntArray wsBdsIndx = new DynamicIntArray(constraints.length); 
  	// Vectors used to record the variable indices to be freed 	
  	DynamicIntArray toFree=null, oldToFree=null;	
  
  	// Initial value of obj. function, gradient and inverse of the Hessian
  	m_f = objectiveFunction(initX);
  	if(Double.isNaN(m_f))
  	    throw new Exception("Objective function value is NaN!");
  	
  	double sum=0;
  	double[] grad=evaluateGradient(initX), oldGrad, oldX,
  	    deltaGrad=new double[l], deltaX=new double[l],
  	    direct = new double[l], x = new double[l];
  	Matrix L = new Matrix(l, l);  // Lower triangle of Cholesky factor 
  	double[] D = new double[l];   // Diagonal of Cholesky factor
  	for(int i=0; i<l; i++){
  	    L.setRow(i, new double[l]);
  	    L.setElement(i,i,1.0);
  	    D[i] = 1.0;
  	    direct[i] = -grad[i];
  	    sum += grad[i]*grad[i];
  	    x[i] = initX[i];
  	    nwsBounds[0][i] = constraints[0][i];
  	    nwsBounds[1][i] = constraints[1][i];
  	    isFixed[i] = false;
  	}	
  	double stpmax = m_STPMX*Math.max(Math.sqrt(sum), l);
  	
  	iterates:
  	for(int step=0; step < m_MAXITS; step++){
  	    if (m_Debug)
  		System.err.println("\nIteration # " + step + ":");	    
  	    
  	    // Try at most one feasible newton step, i.e. 0<lamda<=alpha
  	    oldX = x;
  	    oldGrad = grad;
  	    
  	    // Also update grad
  	    if (m_Debug)
  		System.err.println("Line search ... ");
  	    m_IsZeroStep = false;
  	    x=lnsrch(x, grad, direct, stpmax, 
  		     isFixed, nwsBounds, wsBdsIndx);
  	    if (m_Debug)
  		System.err.println("Line search finished.");
  	    
  	    if(m_IsZeroStep){ // Zero step, simply delete rows/cols of D and L
  		for(int f=0; f<wsBdsIndx.size(); f++){
  		    int idx=wsBdsIndx.elementAt(f);
  		    L.setRow(idx, new double[l]);
  		    L.setColumn(idx, new double[l]);
  		    D[idx] = 0.0;
  		}		
  		grad = evaluateGradient(x);
  		step--;
  	    }
  	    else{
  		// Check converge on x
  		boolean finish = false;
  		double test=0.0;
  		for(int h=0; h<l; h++){
  		    deltaX[h] = x[h]-oldX[h];
  		    double tmp=Math.abs(deltaX[h])/
  			Math.max(Math.abs(x[h]), 1.0);
  		    if(tmp > test) test = tmp;				    
  		}
  		if(test < m_Zero){
  		    if (m_Debug)
  			System.err.println("\nDeltaX converge: "+test);
  		    finish = true;
  		}
  		
  		// Check zero gradient	    
  		grad = evaluateGradient(x);
  		test=0.0;
  		double denom=0.0, dxSq=0.0, dgSq=0.0, newlyBounded=0.0; 
  		for(int g=0; g<l; g++){
  		    if(!isFixed[g]){ 		   
  			deltaGrad[g] = grad[g] - oldGrad[g];		  
  			// Calculate the denominators			    
  			denom += deltaX[g]*deltaGrad[g];
  			dxSq += deltaX[g]*deltaX[g];
  			dgSq += deltaGrad[g]*deltaGrad[g];
  		    }
  		    else // Only newly bounded variables will be non-zero
  			newlyBounded +=  deltaX[g]*(grad[g]-oldGrad[g]);
  		    
  		    // Note: CANNOT use projected gradient for testing 
  		    // convergence because of newly bounded variables
  		    double tmp = Math.abs(grad[g])*
  			Math.max(Math.abs(direct[g]),1.0)/
  			Math.max(Math.abs(m_f),1.0);
  		    if(tmp > test) test = tmp;	
  		}
  		
  		if(test < m_Zero){
  		    if (m_Debug)
  			System.err.println("Gradient converge: "+test);
  		    finish = true;
  		}	    
  		
  		// dg'*dx could be < 0 using inexact lnsrch
  		if(m_Debug)
  		    System.err.println("dg'*dx="+(denom+newlyBounded));	
  		// dg'*dx = 0
  		if(Math.abs(denom+newlyBounded) < m_Zero)
  		    finish = true;
  		
  		int size = wsBdsIndx.size();
  		boolean isUpdate = true;  // Whether to update BFGS formula	    
  		// Converge: check whether release any current constraints
  		if(finish){
  		    if (m_Debug)
  			System.err.println("Test any release possible ...");
  		    	
  		    if(toFree != null)
  			oldToFree = (DynamicIntArray)toFree.copy();
  		    toFree = new DynamicIntArray(wsBdsIndx.size());
  		    
  		    for(int m=size-1; m>=0; m--){
  			int index=wsBdsIndx.elementAt(m);
  			double[] hessian = evaluateHessian(x, index);			
  			double deltaL=0.0;
  			if(hessian != null){
  			    for(int mm=0; mm<hessian.length; mm++)
  				if(!isFixed[mm]) // Free variable
  				    deltaL += hessian[mm]*direct[mm];
  			}
  			
  			// First and second order Lagrangian multiplier estimate
  			// If user didn't provide Hessian, use first-order only
  			double L1, L2;
  			if(x[index] >= constraints[1][index]) // Upper bound
  			    L1 = -grad[index];
  			else if(x[index] <= constraints[0][index])// Lower bound
  			    L1 = grad[index];
  			else
  			    throw new Exception("x["+index+"] not fixed on the"+
  						" bounds where it should have been!");
  			
  			// L2 = L1 + deltaL
  			L2 = L1 + deltaL;			
  			if (m_Debug)
  			    System.err.println("Variable "+index+
  					       ": Lagrangian="+L1+"|"+L2);
  			
  			//Check validity of Lagrangian multiplier estimate
  			boolean isConverge = 
  			    (2.0*Math.abs(deltaL)) < Math.min(Math.abs(L1),
  							      Math.abs(L2));  
  			if((L1*L2>0.0) && isConverge){ //Same sign and converge: valid
  			    if(L2 < 0.0){// Negative Lagrangian: feasible
  				toFree.addElement(index);
  				wsBdsIndx.removeElementAt(m);
  				finish=false; // Not optimal, cannot finish
  			    }
  			}
  			
  			// Although hardly happen, better check it
  			// If the first-order Lagrangian multiplier estimate is wrong,
  			// avoid zigzagging
  			if((hessian==null) && (toFree != null) && toFree.equal(oldToFree)) 
  			    finish = true;           
  		    }
  		    
  		    if(finish){// Min. found
  			if (m_Debug)
  			    System.err.println("Minimum found.");
  			m_f = objectiveFunction(x);
  			if(Double.isNaN(m_f))
  			    throw new Exception("Objective function value is NaN!");	
  			return x;
  		    }
  		    
  		    // Free some variables
  		    for(int mmm=0; mmm<toFree.size(); mmm++){
  			int freeIndx=toFree.elementAt(mmm);
  			isFixed[freeIndx] = false; // Free this variable
  			if(x[freeIndx] <= constraints[0][freeIndx]){// Lower bound
  			    nwsBounds[0][freeIndx] = constraints[0][freeIndx];
  			    if (m_Debug)
  				System.err.println("Free variable "+freeIndx+
  						   " from bound "+ 
  						   nwsBounds[0][freeIndx]);
  			}
  			else{ // Upper bound
  			    nwsBounds[1][freeIndx] = constraints[1][freeIndx];
  			    if (m_Debug)
  				System.err.println("Free variable "+freeIndx+
  						   " from bound "+ 
  						   nwsBounds[1][freeIndx]);
  			}			
  			L.setElement(freeIndx, freeIndx, 1.0);
  			D[freeIndx] = 1.0;
  			isUpdate = false;			
  		    }			
  		}
  		
  		if(denom<Math.max(m_Zero*Math.sqrt(dxSq)*Math.sqrt(dgSq), m_Zero)){
  		    if (m_Debug) 
  			System.err.println("dg'*dx negative!");
  		    isUpdate = false; // Do not update		    
  		}		
  		// If Hessian will be positive definite, update it
  		if(isUpdate){
  		    
  		    // modify once: dg*dg'/(dg'*dx)	
  		    double coeff = 1.0/denom; // 1/(dg'*dx)	
  		    updateCholeskyFactor(L,D,deltaGrad,coeff,isFixed);
  		    
  		    // modify twice: g*g'/(g'*p)	
  		    coeff = 1.0/m_Slope; // 1/(g'*p)
  		    updateCholeskyFactor(L,D,oldGrad,coeff,isFixed);  		    
  		}
  	    }
  	    
  	    // Find new direction 
  	    Matrix LD = new Matrix(l,l); // L*D
  	    double[] b = new double[l];
  	    
  	    for(int k=0; k<l; k++){
  		if(!isFixed[k])  b[k] = -grad[k];
  		else             b[k] = 0.0;
  		
  		for(int j=k; j<l; j++){ // Lower triangle	
  		    if(!isFixed[j] && !isFixed[k])
  			LD.setElement(j, k, L.getElement(j,k)*D[k]);
  		}		
  	    }	    	
  	    
  	    // Solve (LD)*y = -g, where y=L'*direct
  	    double[] LDIR = solveTriangle(LD, b, true, isFixed);	    
  	    LD = null;
  	    
  	    for(int m=0; m<LDIR.length; m++){
  		if(Double.isNaN(LDIR[m]))
  		    throw new Exception("L*direct["+m+"] is NaN!"
  					+"|-g="+b[m]+"|"+isFixed[m]
  					+"|diag="+D[m]);
  	    }
  	    
  	    // Solve L'*direct = y
  	    direct = solveTriangle(L, LDIR, false, isFixed);
  	    for(int m=0; m<direct.length; m++){
  		if(Double.isNaN(direct[m]))
  		    throw new Exception("direct is NaN!");
  	    }
  	    
  	    //System.gc();
  	}
  	
  	if(m_Debug)
  	    System.err.println("Cannot find minimum"+
  			       " -- too many interations!");
  	m_X = x;
  	return null;
      }
      
      /** 
       * Solve the linear equation of TX=B where T is a triangle matrix
       * It can be solved using back/forward substitution, with O(N^2) 
       * complexity
       * @param t the matrix T
       * @param b the vector B
       * @param isLower whether T is a lower or higher triangle matrix
       * @param isZero which row(s) of T are not used when solving the equation. 
       *               If it's null or all 'false', then every row is used.
       * @return the solution of X
       */     
      public static double[] solveTriangle(Matrix t, double[] b, 
  					 boolean isLower, boolean[] isZero){
  	int n = b.length; 
  	double[] result = new double[n];
  	if(isZero == null)
  	    isZero = new boolean[n];
  	
  	if(isLower){ // lower triangle, forward-substitution
  	    int j = 0;
  	    while((j<n)&&isZero[j]){result[j]=0.0; j++;} // go to the first row
  	    
  	    if(j<n){
  		result[j] = b[j]/t.getElement(j,j);
  		
  		for(; j<n; j++){
  		    if(!isZero[j]){
  			double numerator=b[j];
  			for(int k=0; k<j; k++)
  			    numerator -= t.getElement(j,k)*result[k];
  			result[j] = numerator/t.getElement(j,j);
  		    }
  		    else 
  			result[j] = 0.0;
  		}
  	    }
  	}
  	else{ // Upper triangle, back-substitution
  	    int j=n-1;
  	    while((j>=0)&&isZero[j]){result[j]=0.0; j--;} // go to the last row
  	    
  	    if(j>=0){
  		result[j] = b[j]/t.getElement(j,j);
  		
  		for(; j>=0; j--){
  		    if(!isZero[j]){
  			double numerator=b[j];
  			for(int k=j+1; k<n; k++)
  			    numerator -= t.getElement(k,j)*result[k];
  			result[j] = numerator/t.getElement(j,j);
  		    }
  		    else 
  			result[j] = 0.0;
  		}
  	    }
  	}
  	
  	return result;
      }
  
      /**
       * One rank update of the Cholesky factorization of B matrix in BFGS updates,
       * i.e. B = LDL', and B_{new} = LDL' + coeff*(vv') where L is a unit lower triangle
       * matrix and D is a diagonal matrix, and v is a vector.<br/>
       * When coeff > 0, we use C1 algorithm, and otherwise we use C2 algorithm described
       * in ``Methods for Modifying Matrix Factorizations'' 
       *
       * @param L the unit triangle matrix L
       * @param D the diagonal matrix D
       * @param v the update vector v
       * @param coeff the coeffcient of update
       * @param isFixed which variables are not to be updated
       */    
      protected void updateCholeskyFactor(Matrix L, double[] D, 
  					double[] v, double coeff,
  					boolean[] isFixed)
  	throws Exception{
  	double t, p, b;
  	int n = v.length;
  	double[] vp =  new double[n];	
  	for (int i=0; i<v.length; i++)
  	    if(!isFixed[i])
  		vp[i]=v[i];
  	    else
  		vp[i]=0.0;
  	
  	if(coeff>0.0){
  	    t = coeff;	    
  	    for(int j=0; j<n; j++){		
  		if(isFixed[j]) continue;		
  		
  		p = vp[j];
  		double d=D[j], dbarj=d+t*p*p;
  		D[j] = dbarj;
  		
  		b = p*t/dbarj;
  		t *= d/dbarj;
  		for(int r=j+1; r<n; r++){
  		    if(!isFixed[r]){
  			double l=L.getElement(r, j);
  			vp[r] -= p*l;
  			L.setElement(r, j, l+b*vp[r]);
  		    }
  		    else
  		    	L.setElement(r, j, 0.0);
  		}
  	    }
  	}
  	else{
  	    double[] P = solveTriangle(L, v, true, isFixed);	    
  	    t = 0.0;
  	    for(int i=0; i<n; i++)
  		if(!isFixed[i])
  		    t += P[i]*P[i]/D[i];	    	
  	    
  	    double sqrt=1.0+coeff*t;
  	    sqrt = (sqrt<0.0)? 0.0 : Math.sqrt(sqrt);
  	    
  	    double alpha=coeff, sigma=coeff/(1.0+sqrt), rho, theta;
  	    
  	    for(int j=0; j<n; j++){
  		if(isFixed[j]) continue;
  		
  		double d=D[j];
  		p = P[j]*P[j]/d;
  		theta = 1.0+sigma*p;
  		t -= p; 
  		if(t<0.0) t=0.0; // Rounding error
  
  		double plus = sigma*sigma*p*t;
  		if((j<n-1) && (plus <= m_Zero)) 
  		    plus=m_Zero; // Avoid rounding error
  		rho = theta*theta + plus;		
  		D[j] = rho*d;
  		
  		if(Double.isNaN(D[j])){
  		    throw new Exception("d["+j+"] NaN! P="+P[j]+",d="+d+
  					",t="+t+",p="+p+",sigma="+sigma+
  					",sclar="+coeff);
  		}
  		
  		b=alpha*P[j]/(rho*d);
  		alpha /= rho;
  		rho = Math.sqrt(rho);
  		double sigmaOld = sigma;
  		sigma *= (1.0+rho)/(rho*(theta+rho));	 
  		if((j<n-1) && 
  		   (Double.isNaN(sigma) || Double.isInfinite(sigma)))
  		    throw new Exception("sigma NaN/Inf! rho="+rho+
  				       ",theta="+theta+",P["+j+"]="+
  				       P[j]+",p="+p+",d="+d+",t="+t+
  				       ",oldsigma="+sigmaOld);
  		
  		for(int r=j+1; r<n; r++){
  		    if(!isFixed[r]){
  			double l=L.getElement(r, j);
  			vp[r] -= P[j]*l;
  			L.setElement(r, j, l+b*vp[r]);
  		    }
  		    else
  		    	L.setElement(r, j, 0.0);
  		}
  	    }
  	}	
      }
  
    /**
     * Implements a simple dynamic array for ints.
     */
    private class DynamicIntArray
      implements RevisionHandler {
  
      /** The int array. */
      private int[] m_Objects;
  
      /** The current size; */
      private int m_Size = 0;
  
      /** The capacity increment */
      private int m_CapacityIncrement = 1;
    
      /** The capacity multiplier. */
      private int m_CapacityMultiplier = 2;
  
      /**
       * Constructs a vector with the given capacity.
       *
       * @param capacity the vector's initial capacity
       */
      public DynamicIntArray(int capacity) {
        
        m_Objects = new int[capacity];
      }
  
      /**
       * Adds an element to this vector. Increases its
       * capacity if its not large enough.
       *
       * @param element the element to add
       */
      public final void addElement(int element) {
        
        if (m_Size == m_Objects.length) {
  	int[] newObjects;
  	newObjects = new int[m_CapacityMultiplier *
  			     (m_Objects.length +
  			      m_CapacityIncrement)];
  	System.arraycopy(m_Objects, 0, newObjects, 0, m_Size);
  	m_Objects = newObjects;
        }
        m_Objects[m_Size] = element;
        m_Size++;
      }
  
      /**
       * Produces a copy of this vector.
       *
       * @return the new vector
       */
      public final Object copy() {
        
        
        DynamicIntArray copy = new DynamicIntArray(m_Objects.length);
        
        copy.m_Size = m_Size;
        copy.m_CapacityIncrement = m_CapacityIncrement;
        copy.m_CapacityMultiplier = m_CapacityMultiplier;
        System.arraycopy(m_Objects, 0, copy.m_Objects, 0, m_Size);
        return copy;
      }
  
      /**
       * Returns the element at the given position.
       *
       * @param index the element's index
       * @return the element with the given index
       */
      public final int elementAt(int index) {
        
        return m_Objects[index];
      }
      
      /**
       * Check whether the two integer vectors equal to each other
       * Two integer vectors are equal if all the elements are the 
       * same, regardless of the order of the elements
       *
       * @param b another integer vector
       * @return whether they are equal
       */ 
      private boolean equal(DynamicIntArray b){
        if((b==null) || (size()!=b.size()))
  	  return false;
        
        int size=size();
        
        // Only values matter, order does not matter
        int[] sorta=Utils.sort(m_Objects), sortb=Utils.sort(b.m_Objects);
        for(int j=0; j<size;j++)
  	if(m_Objects[sorta[j]] != b.m_Objects[sortb[j]])
  	  return false;
        
        return true;
      }
  
      /**
       * Deletes an element from this vector.
       *
       * @param index the index of the element to be deleted
       */
      public final void removeElementAt(int index) {
        
        System.arraycopy(m_Objects, index + 1, m_Objects, index, 
  		       m_Size - index - 1);
        m_Size--;
      }
  
      /**
       * Removes all components from this vector and sets its 
       * size to zero. 
       */
      public final void removeAllElements() {
        
        m_Objects = new int[m_Objects.length];
        m_Size = 0;
      }
  
      /**
       * Returns the vector's current size.
       *
       * @return the vector's current size
       */
      public final int size() {
        
        return m_Size;
      }
      
      /**
       * Returns the revision string.
       * 
       * @return		the revision
       */
      public String getRevision() {
        return RevisionUtils.extract("$Revision: 5953 $");
      }
    }
  }
  