package org.openscience.cdk.modeling.forcefield;import javax.vecmath.GVector;import org.openscience.cdk.tools.LoggingTool;/** *  Find a decrease direction of the energy fuction from a point of the 3xN coordinates space using the conjugate gradient approach. * *@author     vlabarta *@cdk.module     forcefield * */public class ConjugateGradientMethod {	double uk_FletcherReeves = 0;	double uk_PolankRibiere = 0;	GVector conjugatedGradientDirection = null;	GVector previousConjugatedGradientDirection = null;	boolean orthogonalDirectionsProperty = true;	GVector diffgk_gkminus1 = null;	private LoggingTool logger;	/**	 *  Constructor for the ConjugateGradientMethod object	 */	public ConjugateGradientMethod() {		logger = new LoggingTool(this);	}	/**	 *  Fletcher-Reeves: uk = gk gk / gk-1 gk-1	 *	 */	public void initialize(GVector gradient) {		conjugatedGradientDirection = new GVector(gradient);		//conjugatedGradientDirection.normalize();		conjugatedGradientDirection.scale(-1);		previousConjugatedGradientDirection = new GVector(gradient.getSize());		diffgk_gkminus1 = new GVector(gradient.getSize());	}	/**	 *  Fletcher-Reeves: uk = gk gk / gk-1 gk-1	 *	 */	public void setFletcherReeves_uk(GVector gkminus1, GVector gk) {		uk_FletcherReeves = gk.dot(gk) / gkminus1.dot(gkminus1);				logger.debug("uk_FletcherReeves = " + uk_FletcherReeves);		return;	}	/**	 *  Polak-Ribiere plus: uk = Max(0, (gk - gk-1) gk / gk-1 gk-1)	 *	 */	public void setPolankRibierePlus_uk(GVector gkminus1, GVector gk) {		diffgk_gkminus1.set(gk);		diffgk_gkminus1.sub(gkminus1);		uk_PolankRibiere = Math.max(0, diffgk_gkminus1.dot(gk) / gkminus1.dot(gkminus1));		if (uk_PolankRibiere == 0) {			//logger.debug("uk_PolankRibiere == 0");		}		return;	}	/**	 *  Check if two consecutive conjugate gradient direction are mutually orthogonal.	 *	 * @param  pkminus1  	Conjugate Gradient direction at coordinates Xk-1	 * @param  pk  		Conjugate Gradient direction at coordinates Xk	 */	private void checkingOrthogonality(GVector pkminus1, GVector pk) {		//logger.debug("Math.abs(pk.dot(pkminus1)) / Math.pow(pk.norm(),2) = " + Math.abs(pk.dot(pkminus1)) / Math.pow(pk.norm(),2));		//logger.debug("Math.abs(pk.dot(pkminus1)) / Math.pow(pk.normSquared(),2) = " + Math.abs(pk.dot(pkminus1)) / Math.pow(pk.normSquared(),2));		if (Math.abs(pk.dot(pkminus1)) / Math.pow(pk.normSquared(),2) >= 0.1) {			orthogonalDirectionsProperty = false;			//logger.debug("orthogonalDirectionsProperty = false");		} else {orthogonalDirectionsProperty = true;}	}	/**	 *  Restart conjugate gradient direction: assign the gradient of xk as conjugate gradient direction.	 *	 * @param  gk  		gradient at coordinates Xk	 */	private void restartConjugateGradient(GVector gk) {		conjugatedGradientDirection.set(gk);		//conjugatedGradientDirection.normalize();		conjugatedGradientDirection.scale(-1);		//logger.debug("vectorvk : " + direction);	}	/**	 *  Set the new direction conjugated to the previous direction: vk=-gk + uk vk-1	 *	 * @param  gradient  				gradient at coordinates Xk	 * @param  previousGradient  		gradient at coordinates Xk-1	 */	private void setConjugateGradientDirection(GVector gradient, GVector previousGradient) {		setPolankRibierePlus_uk(previousGradient,gradient);		previousConjugatedGradientDirection.scale(uk_PolankRibiere);		conjugatedGradientDirection.set(gradient);		//conjugatedGradientDirection.normalize();		conjugatedGradientDirection.scale(-1);		conjugatedGradientDirection.add(previousConjugatedGradientDirection);		previousConjugatedGradientDirection.scale(1/uk_PolankRibiere);		//direction.normalize();		//logger.debug("vector direction : " + direction);	}	/**	 *  Calculate the conjugate gradient direction.	 *	 * @param  gradient  			Energy function gradient at coordinates Xk	 * @param  previousGradient  	Energy function gradient at coordinates Xk-1	*/	public void setDirection(GVector gradient, GVector previousGradient) {		previousConjugatedGradientDirection.set(conjugatedGradientDirection);		setConjugateGradientDirection(gradient, previousGradient);		checkingOrthogonality(previousConjugatedGradientDirection,conjugatedGradientDirection);		if (orthogonalDirectionsProperty == false) {restartConjugateGradient(gradient);}	}			}