
com.actelion.research.chem.optimization.OptimizerLBFGS Maven / Gradle / Ivy
/*
* 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;
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy