com.actelion.research.chem.phesa.Lnsrch Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of openchemlib Show documentation
Open Source Chemistry Library
There is a newer version: 2024.12.1
package com.actelion.research.chem.phesa;

import com.actelion.research.chem.phesa.EvaluableOverlap;
import java.util.Arrays;





/**
 * taken from DD_chem3d
 *  
 */
class Lnsrch{

	
	/**
	 * Minimize the function according to the line search algorithm. (ie find lamda so that newpos = pos+lambda*dir minimizes the energy)
	 * This function expects the initial FGValue, the direction of search and the initial move
	 * 
	 *  http://www.mcs.anl.gov/~anitescu/CLASSES/2012/LECTURES/S310-2012-lect4.pdf
	 *  
	 * @param f0
	 * @param function
	 * @param grad
	 * @param dir
	 * @param fMove
	 * @return
	 */
	public static final Object[] minimizeEnergyAroundDirection(final Evaluable function, double f0, final double[] grad, final double[] dir, final double fMove) {
		final double CAPPA = .9;
		final double STPMIN = 1e-6;
		final double STPMAX = .1;
		double fA=0, fB=0, fC=0;
		double slopeA=0, slopeB=0, slopeC=0;
		double cube = 0;
		//Compute length of Gradient
		final double len = grad.length;
		final double sNorm = OptimizerLBFGS.getNorm(dir);

		//Normalize the search vector and find the projected gradient
		double slope = 0;
		for(int i=0; iSTPMAX) step = STPMAX;
		else if(step=fA || slopeB*slopeA<0) {
						break;
					}
					
					//Adapt step
					if(slopeB>slopeA) {
						double parab = (fA - fB) / (slopeB - slopeA);
						if(parab>2*step) step = 2 * step;
						else if(parab<2*step) step = step / 2;
						else step = parab;
					} else {
						step*=2;
					}
					if(step>STPMAX) step = STPMAX;
					else if(stepstep) cube = step/2;
					if(cube<0 || cube>step) {
						return new Object[]{fC, grad, Boolean.FALSE}; //Interpolation error, stop here
					}

					
					lambda -= cube;	
					
					//Get new function and gradient
					move(function, dir, lambda, initial,v);
					fC = function.getFGValue(grad);
					slopeC = 0; for(int i=0; i=STPMIN) {
							if( (slopeA*slopeB<0 && slopeA*slopeC<0) || (slopeA*slopeB>=0 && (slopeA*slopeC<0 || fA B
								fB = fC;
								slopeB = slopeC;
								step -= cube;
							} else {
								//C becomes the left limit -> A
								fA = fC;
								slopeA = slopeC;
								step = cube;
								lambda += cube;					
							}
						} else {
	//						if((fC>fA && fC>fB) || cubStepf0) {
					move(function, dir, lambda, initial,v);
					f0 = function.getFGValue(grad);
					//System.err.println("ERR fL>f0");
					return new Object[]{f0, grad, Boolean.FALSE};
				}
				f0 = fL;
				if(sgL>0) {
					lambda = 0;
					for(int i=0; iinitialF)
				move(function, dir, 0, initial,v);
			return new Object[]{f0, grad, Boolean.FALSE};

		} catch(Exception e) {
			e.printStackTrace();
			function.setState(initial);
			f0 = function.getFGValue(grad);
			return new Object[]{new Double(f0), grad, Boolean.FALSE};				
		}
	}	

	private final static void move(Evaluable eval, double[] dir, double lambda, double[] transformOld, double[] transform) {
		for(int i=0;i