26 const double EPS = 3e-8;
51 template <
typename EnergyFunctor>
52 void linearSearch(
unsigned int dim,
double *oldPt,
double oldVal,
double *grad,
53 double *dir,
double *newPt,
double &newVal,
54 EnergyFunctor func,
double maxStep,
int &resCode) {
60 const unsigned int MAX_ITER_LINEAR_SEARCH = 1000;
61 double sum = 0.0, slope = 0.0, test = 0.0,
lambda = 0.0;
62 double lambda2 = 0.0, lambdaMin = 0.0, tmpLambda = 0.0, val2 = 0.0;
68 for (
unsigned int i = 0; i < dim; i++) sum += dir[i] * dir[i];
73 for (
unsigned int i = 0; i < dim; i++) dir[i] *= maxStep / sum;
79 for (
unsigned int i = 0; i < dim; i++) {
80 slope += dir[i] * grad[i];
87 for (
unsigned int i = 0; i < dim; i++) {
88 double temp = fabs(dir[i]) / std::max(fabs(oldPt[i]), 1.0);
89 if (temp > test) test = temp;
92 lambdaMin = MOVETOL / test;
95 while (it < MAX_ITER_LINEAR_SEARCH) {
102 for (
unsigned int i = 0; i < dim; i++) {
103 newPt[i] = oldPt[i] +
lambda * dir[i];
105 newVal = func(newPt);
107 if (newVal - oldVal <= FUNCTOL *
lambda * slope) {
115 tmpLambda = -slope / (2.0 * (newVal - oldVal - slope));
117 double rhs1 = newVal - oldVal -
lambda * slope;
118 double rhs2 = val2 - oldVal - lambda2 * slope;
119 double a = (rhs1 / (
lambda *
lambda) - rhs2 / (lambda2 * lambda2)) /
122 lambda * rhs2 / (lambda2 * lambda2)) /
125 tmpLambda = -slope / (2.0 * b);
127 double disc = b * b - 3 * a * slope;
130 }
else if (b <= 0.0) {
131 tmpLambda = (-b + sqrt(disc)) / (3.0 * a);
133 tmpLambda = -slope / (b + sqrt(disc));
136 if (tmpLambda > 0.5 *
lambda) {
147 for (
unsigned int i = 0; i < dim; i++) {
156 delete[] hessDGrad; \ 159 delete[] invHessian; \ 186 template <
typename EnergyFunctor,
typename GradientFunctor>
187 int minimize(
unsigned int dim,
double *pos,
double gradTol,
188 unsigned int &numIters,
double &funcVal, EnergyFunctor func,
189 GradientFunctor gradFunc,
unsigned int snapshotFreq,
191 unsigned int maxIts = MAXITS) {
196 double sum, maxStep, fp;
198 double *grad, *dGrad, *hessDGrad;
202 grad =
new double[dim];
203 dGrad =
new double[dim];
204 hessDGrad =
new double[dim];
205 newPos =
new double[dim];
206 xi =
new double[dim];
207 invHessian =
new double[dim * dim];
208 snapshotFreq = std::min(snapshotFreq, maxIts);
215 memset(invHessian, 0, dim * dim *
sizeof(
double));
216 for (
unsigned int i = 0; i < dim; i++) {
217 unsigned int itab = i * dim;
219 invHessian[itab + i] = 1.0;
222 sum += pos[i] * pos[i];
225 maxStep = MAXSTEP * std::max(sqrt(sum), static_cast<double>(dim));
227 for (
unsigned int iter = 1; iter <= maxIts; iter++) {
232 linearSearch(dim, pos, fp, grad, xi, newPos, funcVal, func, maxStep,
241 for (
unsigned int i = 0; i < dim; i++) {
242 xi[i] = newPos[i] - pos[i];
244 double temp = fabs(xi[i]) / std::max(fabs(pos[i]), 1.0);
245 if (temp > test) test = temp;
251 if (snapshotVect && snapshotFreq) {
253 snapshotVect->push_back(s);
261 double gradScale = gradFunc(pos, grad);
265 double term = std::max(funcVal * gradScale, 1.0);
266 for (
unsigned int i = 0; i < dim; i++) {
267 double temp = fabs(grad[i]) * std::max(fabs(pos[i]), 1.0);
268 test = std::max(test, temp);
269 dGrad[i] = grad[i] - dGrad[i];
274 if (test < gradTol) {
275 if (snapshotVect && snapshotFreq) {
277 snapshotVect->push_back(s);
289 double fac = 0, fae = 0, sumDGrad = 0, sumXi = 0;
290 for (
unsigned int i = 0; i < dim; i++) {
292 unsigned int itab = i * dim;
294 for (
unsigned int j = 0; j < dim; j++) {
295 hessDGrad[i] += invHessian[itab + j] * dGrad[j];
299 double *ivh = &(invHessian[i * dim]);
300 double &hdgradi = hessDGrad[i];
303 for (
unsigned int j = 0; j < dim; ++j, ++ivh, ++dgj) {
304 hdgradi += *ivh * *dgj;
307 fac += dGrad[i] * xi[i];
308 fae += dGrad[i] * hessDGrad[i];
309 sumDGrad += dGrad[i] * dGrad[i];
310 sumXi += xi[i] * xi[i];
312 if (fac > sqrt(EPS * sumDGrad * sumXi)) {
314 double fad = 1.0 / fae;
315 for (
unsigned int i = 0; i < dim; i++) {
316 dGrad[i] = fac * xi[i] - fad * hessDGrad[i];
319 for (
unsigned int i = 0; i < dim; i++) {
320 unsigned int itab = i * dim;
322 for (
unsigned int j = i; j < dim; j++) {
323 invHessian[itab + j] += fac * xi[i] * xi[j] -
324 fad * hessDGrad[i] * hessDGrad[j] +
325 fae * dGrad[i] * dGrad[j];
326 invHessian[j * dim + i] = invHessian[itab + j];
329 double pxi = fac * xi[i], hdgi = fad * hessDGrad[i],
330 dgi = fae * dGrad[i];
331 double *pxj = &(xi[i]), *hdgj = &(hessDGrad[i]), *dgj = &(dGrad[i]);
332 for (
unsigned int j = i; j < dim; ++j, ++pxj, ++hdgj, ++dgj) {
333 invHessian[itab + j] += pxi * *pxj - hdgi * *hdgj + dgi * *dgj;
334 invHessian[j * dim + i] = invHessian[itab + j];
340 for (
unsigned int i = 0; i < dim; i++) {
341 unsigned int itab = i * dim;
344 for (
unsigned int j = 0; j < dim; j++) {
345 xi[i] -= invHessian[itab + j] * grad[j];
348 double &pxi = xi[i], *ivh = &(invHessian[itab]), *gj = grad;
349 for (
unsigned int j = 0; j < dim; ++j, ++ivh, ++gj) {
354 if (snapshotVect && snapshotFreq && !(iter % snapshotFreq)) {
356 snapshotVect->push_back(s);
357 newPos =
new double[dim];
380 template <
typename EnergyFunctor,
typename GradientFunctor>
381 int minimize(
unsigned int dim,
double *pos,
double gradTol,
382 unsigned int &numIters,
double &funcVal, EnergyFunctor func,
383 GradientFunctor gradFunc,
double funcTol = TOLX,
384 unsigned int maxIts = MAXITS) {
385 return minimize(dim, pos, gradTol, numIters, funcVal, func,
386 gradFunc, 0, NULL, funcTol, maxIts);
const double MAXSTEP
Default maximim step size in the minimizer.
#define CHECK_INVARIANT(expr, mess)
const double MOVETOL
Default tolerance for x changes in the minimizer.
const double lambda
scaling factor for rBO correction
const double FUNCTOL
Default tolerance for function convergence in the minimizer.
RDKIT_OPTIMIZER_EXPORT int HEAD_ONLY_LIBRARY
RDKIT_OPTIMIZER_EXPORT int REALLY_A_HEADER_ONLY_LIBRARY
void linearSearch(unsigned int dim, double *oldPt, double oldVal, double *grad, double *dir, double *newPt, double &newVal, EnergyFunctor func, double maxStep, int &resCode)
Do a Quasi-Newton minimization along a line.
std::vector< Snapshot > SnapshotVect
#define RDUNUSED_PARAM(x)
#define RDKIT_OPTIMIZER_EXPORT
const double TOLX
Default direction vector tolerance in the minimizer.
int minimize(unsigned int dim, double *pos, double gradTol, unsigned int &numIters, double &funcVal, EnergyFunctor func, GradientFunctor gradFunc, unsigned int snapshotFreq, RDKit::SnapshotVect *snapshotVect, double funcTol=TOLX, unsigned int maxIts=MAXITS)
Do a BFGS minimization of a function.
#define PRECONDITION(expr, mess)
const double EPS
Default gradient tolerance in the minimizer.
const int MAXITS
Default maximum number of iterations.