package org.openscience.cdk.modeling.forcefield;import javax.vecmath.GVector;//import org.openscience.cdk.tools.LoggingTool;/** *  * * @author     vlabarta *@cdk.module     forcefield *  */public class LineSearchForTheWolfeConditions {		//initial values	private GVector x = null;	private double linearFunctionInAlpha0;	private GVector dfx = null;	private GVector direction = null;	private double linearFunctionDerivativeInAlpha0;	private IPotentialFunction pf = null;	private double alphaInitialStep;		//line search algorithm	private double[] alpha = new double[3];	private double[] linearFunctionInAlpha = new double[3];	private double[] linearFunctionDerivativeInAlpha = new double[3];	private GVector[] dfInAlpha = new GVector[3];	private double[] brentStep = new double[3];	private final double c1 = 0.0001;	private double c2;		//private double linearFunctionGoldenAlpha;	private double linearFunctionAlphaInterpolation;		public boolean derivativeSmallEnough = true;	public double alphaOptimum;	public double linearFunctionInAlphaOptimum;	public GVector dfOptimum = null;		//zoom	private double alphaj;	private double linearFunctionInAlphaj;	private double linearFunctionDerivativeInAlphaj;	private GVector dfInAlphaj;	private int functionEvaluationNumber;		//energy function evaluation	private GVector xAlpha = null;	//interpolation	private double a;	private double b;		//cubic interpolation	private double alphaTemporal;	private double linearFunctionInAlphaTemporal;	private double linearFunctionDerivativeInAlphaTemporal;	private double d1;	private double d2;	private double alphaiplus1;		//private LoggingTool logger;		public LineSearchForTheWolfeConditions(IPotentialFunction pfUser, String method) {		this.pf = pfUser;		if ((method == "sdm") | (method == "cgm")) {c2 = 0.07;}		else {c2 = 0.9;} 	}		public void initialize(GVector xUser, double fxUser, GVector dfxUser, GVector directionUser, double linearFunctionDerivativeUser, double alphaInitialStepUser) {		this.x = xUser;		this.linearFunctionInAlpha0 = fxUser;		this.dfx = dfxUser;		//logger.debug("derivativeSmallEnough = " + this.derivativeSmallEnough);		this.direction = directionUser;		this.linearFunctionDerivativeInAlpha0 = linearFunctionDerivativeUser;		//logger.debug("linearFunctionDerivativeInAlpha0 = " + linearFunctionDerivativeInAlpha0);		this.alphaOptimum = 0;		this.linearFunctionInAlphaOptimum = linearFunctionInAlpha0;		dfOptimum = this.dfx;		this.alphaInitialStep = alphaInitialStepUser;		this.derivativeSmallEnough = false;		this.xAlpha = new GVector(x.getSize());	}		 /*	Line Search Algorithm for the Wolfe conditions. Jorge Nocedal and Stephen J.Wright. Numerical Optimization. 1999.	 * The algorithm has two stages. This first stage begins with a trial estimate alpha1, 	 * and keeps increasing it until it finds either an acceptable step length or an interval 	 * that brackets the desired step lengths. In the later case, the second stage is invoked 	 * by calling a function called zoom, which successively decreases the size of the interval 	 * until an acceptable step length is identified. 	 * 	 * @param alphaMax				Maximum step length	 */	public void lineSearchAlgorithm (double alphaMax) {				//logger.debug("Line search for the strong wolfe conditions");		alpha[0] = 0.0;		linearFunctionInAlpha[0] = linearFunctionInAlpha0; 		linearFunctionDerivativeInAlpha[0] = linearFunctionDerivativeInAlpha0;	//To Analyse the possibility of eliminate linearFunctionDerivativeInAlpha[0]		dfInAlpha[0] = this.dfx;				alpha[1] = this.alphaInitialStep;				//logger.debug("alpha[1] = this.alphaInitialStep = " + alpha[1]);				brentStep[0] = alpha[0];		brentStep[1] = alpha[1];				int i=1;				this.functionEvaluationNumber = 0;		if (alpha[1] > alphaMax) {			alpha[1] = alphaMax;			//logger.debug("line search algorithm error: alphaInitialStep > alphaMax");		}		//	alpha[1] = alphaMax/2;		//}		try {		do {			if (alpha[1] == 0) {				System.out.println("alpha[1] == 0");				break;			}						//logger.debug("alpha[" + i + "] = " + alpha[i]);			linearFunctionInAlpha[i] = evaluateEnergyFunction(alpha[i]);			//logger.debug("linearFunctionInAlpha[" + i + "] = " + linearFunctionInAlpha[i]);						if ((linearFunctionInAlpha[i] > linearFunctionInAlpha[0] + c1 * alpha[i] * linearFunctionDerivativeInAlpha[0]) | 					((linearFunctionInAlpha[i] >= linearFunctionInAlpha[i-1]) & (i>1))) {			//The interval alpha[i-1] and alpha[i] brackets the desired step lengths.				//logger.debug("zoom(" + alpha[i-1] + ", " + linearFunctionInAlpha[i-1] + ", " + linearFunctionDerivativeInAlpha[i-1] + ", " + dfInAlpha[i-1] + ", " + alpha[i] + ", " + linearFunctionInAlpha[i] + ")");				//dfInAlpha[i] = evaluateEnergyFunctionDerivative(alpha[i]);				//linearFunctionDerivativeInAlpha[i] = dfInAlpha[i].dot(direction);				//zoom(alpha[i-1], linearFunctionInAlpha[i-1], linearFunctionDerivativeInAlpha[i-1], dfInAlpha[i-1], alpha[i], linearFunctionInAlpha[i], linearFunctionDerivativeInAlpha[i], dfInAlpha[i]);				zoom(alpha[i-1], linearFunctionInAlpha[i-1], linearFunctionDerivativeInAlpha[i-1], dfInAlpha[i-1], alpha[i], linearFunctionInAlpha[i]);				break;			} 			//The first strong Wolfe condition is satisfied for alpha[i].			dfInAlpha[i] = evaluateEnergyFunctionDerivative(alpha[i]);			//logger.debug("dfOptimum = " + dfOptimum);			linearFunctionDerivativeInAlpha[i] = dfInAlpha[i].dot(direction);			//logger.debug("linearFunctionDerivativeInAlpha[" + i + "] = " + linearFunctionDerivativeInAlpha[i]);						if (Math.abs(linearFunctionDerivativeInAlpha[i]) <= -c2 * linearFunctionDerivativeInAlpha[0]) { //The second strong Wolfe condition is also satisfied for alpha[i]				//logger.debug("The second strong Wolfe condition is also satisfied for " + alpha[i]);				alphaOptimum = alpha[i];				linearFunctionInAlphaOptimum = linearFunctionInAlpha[i];				dfOptimum = dfInAlpha[i];				//logger.debug("alphaOptimun = " + alphaOptimum);				//logger.debug("linearFunctionInAlphaOptimun = " + linearFunctionInAlphaOptimum);				//logger.debug("dfOptimum = " + dfOptimum);				this.derivativeSmallEnough = true;				break;			}						if (linearFunctionDerivativeInAlpha[i] >= 0) {		//The interval alpha[i-1] and alpha[i] brackets the desired step lengths.				/*System.out.println("zoom(" + alpha[i-1] + ", " + linearFunctionInAlpha[i-1] + ", " + linearFunctionDerivativeInAlpha[i-1] + ", " + dfInAlpha[i-1] + ", " + 						alpha[i] + ", " + linearFunctionInAlpha[i] + ")");*/								/*zoom(alpha[i], linearFunctionInAlpha[i], linearFunctionDerivativeInAlpha[i], dfInAlpha[i], 						alpha[i-1], linearFunctionInAlpha[i-1], linearFunctionDerivativeInAlpha[i], dfInAlpha[i]);*/				zoom(alpha[i-1], linearFunctionInAlpha[i-1], linearFunctionDerivativeInAlpha[i-1], dfInAlpha[i-1], 						alpha[i], linearFunctionInAlpha[i]);				break;			}					if (alpha[i] == alphaMax) {					//logger.debug("LINE SEARCH ALGORITHM WAS TERMINATE EARLIER BECAUSE alpha[i] == alphaMax");				alphaOptimum = alpha[i];				linearFunctionInAlphaOptimum = linearFunctionInAlpha[i];				dfOptimum = dfInAlpha[i];				//logger.debug("alphaOptimun = " + alphaOptimum);				//logger.debug("linearFunctionInAlphaOptimun = " + linearFunctionInAlphaOptimum);				//logger.debug("dfOptimum = " + dfOptimum);				break;			}						functionEvaluationNumber = functionEvaluationNumber + 1;			if (functionEvaluationNumber == 10) {				//logger.debug("LINE SEARCH ALGORITHM WAS TERMINATE EARLIER BECAUSE functionEvaluationNumber == 10");				alphaOptimum = alpha[i];				linearFunctionInAlphaOptimum = linearFunctionInAlpha[i];				dfOptimum = dfInAlpha[i];				//logger.debug("alphaOptimun = " + alphaOptimum);				//logger.debug("linearFunctionInAlphaOptimun = " + linearFunctionInAlphaOptimum);				//logger.debug("dfOptimum = " + dfOptimum);				break;			}						if (i>1) {				brentStep[0] = brentStep[1];				brentStep[1] = brentStep[2];				alpha[1] = alpha[2];				linearFunctionInAlpha[1] = linearFunctionInAlpha[2];				linearFunctionDerivativeInAlpha[1] = linearFunctionDerivativeInAlpha[2];				dfInAlpha[1] = dfInAlpha[2];			} 						brentStep[2] = brentStep[1] + 1.618 * (brentStep[1]-brentStep[0]);			//logger.debug("brentStep[2] = " + brentStep[2]);			if (brentStep[2] > alphaMax) {brentStep[2] = alphaMax;}			/*linearFunctionInBrentStep = this.evaluateEnergyFunction(brentStep[2]);			linearFunctionDerivativeInBrentStep = this.evaluateEnergyFunctionDerivative(brentStep[2]).dot(direction);			*/			alpha[2] = brentStep[2];			/*alpha[2] = this.cubicInterpolation(alpha[1], linearFunctionInAlpha[1], linearFunctionDerivativeInAlpha[1], 					brentStep[2], linearFunctionInBrentStep, linearFunctionDerivativeInBrentStep, alpha[1], brentStep[2]);					*/			i=2;					} while ((alpha[2] <= alphaMax) & (alpha[1] < alpha[2]) & (functionEvaluationNumber < 10));				} catch (Exception exception) {        	System.out.println("Line search for the strong wolfe conditions: " + exception.getMessage());        	System.out.println(exception);        }			}		 /*	Each iteration of zoom generates an iterate alphaj between alphaLow and alphaHigh, 	 * and then replaces one of these endpoints by alphaj in such a way that the properties 	 * (a), (b) and (c) continue to hold.	 * (a)The interval bounded by alphaLow and alphaHigh contains step lengths that satisfy the strong Wolfe conditions.       * (b)alphaLow is, among all step lengths generated so far and satisfying the sufficient decrease condition,      * the one giving the smallest function value.     * (c)alphaHigh is chosen so that linearFunctionDerivativeInAlphaj * (alphaHigh-alphaLow) < 0     *      	 *@param  alphaLow              				Among all step lengths generated so far and satisfying the sufficient decrease condition, the one giving the smallest function value.   	 *@param  linearFunctionInAlphaLow       		Function value at alphaLow.   	 *@param  linearFunctionDerivativeInAlphaLow	Derivative value at alphaLow.   	 *@param  dfInAlphaLow              			Gradient at alphaLow.   	 *@param  alphaHigh              				AlphaHigh is chosen so that linearFunctionDerivativeInAlphaj * (alphaHigh-alphaLow) < 0   	 *@param  linearFunctionInAlphaHigh             Function value at alphaHigh.	 */	private void zoom (double alphaLow, double linearFunctionInAlphaLow, double linearFunctionDerivativeInAlphaLow, GVector dfInAlphaLow, 			double alphaHigh, double linearFunctionInAlphaHigh) {				//logger.debug("zoom");				functionEvaluationNumber = 0;				/*double a;		double b;		if (alphaLow < alphaHigh) {a = alphaLow; b = alphaHigh;}		else {a = alphaHigh; b = alphaLow;}		*/				do {			//Interpolation 						//alphaj = this.cubicInterpolation(alphaLow, linearFunctionInAlphaLow, linearFunctionDerivativeInAlphaLow, alphaHigh, linearFunctionInAlphaHigh, linearFunctionDerivativeInAlphaHigh, a, b);			/*System.out.println("interpolation(" + alphaLow + ", " + linearFunctionInAlphaLow + ", " + linearFunctionDerivativeInAlphaLow + ", "					+ alphaHigh + ", " + linearFunctionInAlphaHigh + ");");*/			alphaj = this.interpolation(alphaLow, linearFunctionInAlphaLow, linearFunctionDerivativeInAlphaLow, alphaHigh, linearFunctionInAlphaHigh);			//logger.debug("alphaj = " + alphaj);			linearFunctionInAlphaj = this.linearFunctionAlphaInterpolation;			//logger.debug("linearFunctionInAlphaj = " + linearFunctionInAlphaj);						if ((linearFunctionInAlphaj > linearFunctionInAlpha0 + c1 * alphaj * linearFunctionDerivativeInAlpha0) | //The interval 0 and alphaj brackets the desired step lengths.					(linearFunctionInAlphaj >= linearFunctionInAlphaLow)) {							//logger.debug("The minimum is between alpha1 and alphaj");				alphaHigh = alphaj;				linearFunctionInAlphaHigh = linearFunctionInAlphaj;				//dfInAlphaHigh = this.evaluateEnergyFunctionDerivative(alphaHigh); 				//linearFunctionDerivativeInAlphaHigh = dfInAlphaHigh.dot(direction);			} 			else {				dfInAlphaj = evaluateEnergyFunctionDerivative(alphaj);				linearFunctionDerivativeInAlphaj = dfInAlphaj.dot(direction);				//logger.debug("linearFunctionDerivativeInAlphaj = " + linearFunctionDerivativeInAlphaj);				if (Math.abs(linearFunctionDerivativeInAlphaj) <= -c2 * linearFunctionDerivativeInAlpha0) { //alphaj satisfied the second strong Wolfe condition.					//logger.debug("Derivative small enough : " + Math.abs(linearFunctionDerivativeInAlphaj) + " <= " + (-c2 * linearFunctionDerivativeInAlpha0));					this.derivativeSmallEnough = true;					alphaOptimum = alphaj;					linearFunctionInAlphaOptimum = linearFunctionInAlphaj;					dfOptimum = dfInAlphaj;					//logger.debug("alphaOptimun = " + alphaOptimum);					//logger.debug("linearFunctionInAlphaOptimun = " + linearFunctionInAlphaOptimum);					break;