com.actelion.research.chem.optimization.OptimizerLBFGS Maven / Gradle / Ivy

Go to download
/*
* Copyright (c) 1997 - 2016
* Actelion Pharmaceuticals Ltd.
* Gewerbestrasse 16
* CH-4123 Allschwil, Switzerland
*
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice, this
*    list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
*    this list of conditions and the following disclaimer in the documentation
*    and/or other materials provided with the distribution.
* 3. Neither the name of the the copyright holder nor the
*    names of its contributors may be used to endorse or promote products
*    derived from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
* ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
*/
package com.actelion.research.chem.optimization;

import java.util.Arrays;




/**
 *
 * taken from DD_chem3d, small changes necessary because a different output structure is needed
 * it returns not only a value (the objective function), but the transformation array for achieving the best alignment
 */
public class OptimizerLBFGS   {

	int maxIterations;	
	double minRMS;
	

	
	public OptimizerLBFGS(int maxIterations, double minRMS){
		this.maxIterations = maxIterations;
		this.minRMS = minRMS;
		
	}

	

	/**
	 * Optimization routine using the limited Broyden-Fletcher-Goldfarb-Shanno
	 * algorithm
	 */

	public synchronized double[] optimize(Evaluable eval) {
	//public synchronized double optimize(AbstractEvaluable eval) {
		//System.out.println("new opti");
		double[] initial = eval.getState();
		int N = initial.length;

		final int MSAV = N; // Math.max(1, Math.min(12, N));
		double[] alpha = new double[MSAV];
		double[] rho = new double[MSAV];
		double gamma = 1;
		int m = 0;
		int nErrors = 0;

		// evaluate the function and get the initial gradient
		double[] grad = new double[N];
		double f = eval.getFGValue(grad);
		//System.out.println(f);
		double fOld = f;
		double f0 = f;
		double gNorm;
		double fMove = 0;
		
		if (N == 0) {
			return eval.getState();
		}

			//eturn f;
		double[] oldX = initial;
		double[] oldGradient = new double[N];
		double[][] s = new double[MSAV][N];
		double[][] y = new double[MSAV][N];
		double[] h0 = new double[N];
		double[] q = new double[N];
		double[] r = new double[N];

		boolean restart = true;
		int mUse = 0;
		int iteration ;
		for (iteration = 1; iteration <= maxIterations; iteration++) {
			if (restart) {

				mUse = 0;
				f = eval.getFGValue(grad);
				//System.out.println(f);
				gamma = 1;
				fMove = .25 * getNorm(grad);
				restart = false;
			}

			gNorm = getNorm(grad);
			double RMS = gNorm / Math.sqrt(N);

			if (RMS < minRMS) {
				break;

			} else if (nErrors > 2) {
				break;
			}
			// Estimate Hessian diagonal
			m = (m + 1) % MSAV;
			Arrays.fill(h0, gamma);
			System.arraycopy(grad, 0, q, 0, N);

			int k = m;
			for (int j = 0; j < mUse; j++) {
				k = k == 0 ? MSAV - 1 : k - 1;
				alpha[k] = 0;
				for (int i = 0; i < N; i++)
					alpha[k] += s[k][i] * q[i];
				alpha[k] *= rho[k];
				for (int i = 0; i < N; i++)
					q[i] -= y[k][i] * alpha[k];
			}
			for (int i = 0; i < N; i++)
				r[i] = h0[i] * q[i];
			for (int j = 0; j < mUse; j++) {
				double beta = 0;
				for (int i = 0; i < N; i++)
					beta += y[k][i] * r[i];
				beta *= rho[k];
				for (int i = 0; i < N; i++)
					r[i] += s[k][i] * (alpha[k] - beta);
				k = (k + 1) % MSAV;
			}

			// set search direction and store current point and gradient
			for (int i = 0; i < N; i++) {
				r[i] = -r[i];
			}

			// Memorize position
			oldX = eval.getState();
			System.arraycopy(grad, 0, oldGradient, 0, N);

			// perform line search along the new conjugate direction
			Object[] res = Lnsrch.minimizeEnergyAroundDirection(eval, f, grad, r, fMove);
			f = (Double) res[0];
			grad = (double[]) res[1];
			if (res[2] == Boolean.FALSE) {
				nErrors++;
				restart = true;
			}

			// Update variables
			double ys = 0, yy = 0;
			double[] newState = eval.getState();
			for (int i = 0; i < N; i++) {
				s[m][i] = newState[i] - oldX[i];
				y[m][i] = grad[i] - oldGradient[i];
				ys += y[m][i] * s[m][i];
				yy += y[m][i] * y[m][i];
			}
			gamma = Math.abs(ys / yy);
			if (ys == 0) {
				restart = true;
				continue;
			}

			rho[m] = 1.0 / ys;
			fMove = fOld - f;
			fOld = f;
			mUse = Math.min(mUse + 1, MSAV);

		}

		if (f > f0) {
			eval.setState(initial);
			f = f0;
		}
		return eval.getState();
	}

	public final static double getRMS(double[] vector) {return Math.sqrt(getNormSq(vector) / vector.length);}	
	
	public static double getNorm(double[] vector) {
		 return Math.sqrt(getNormSq(vector));

}
	
	public final static double getNormSq(double[] vector) {
		double res = 0;
		for (int i = 0; i < vector.length; i++) res += vector[i] * vector[i];
		return res;
	}


}