Skip to content

GitLab

Commit df6bbb23 authored by peastman's avatar peastman
Browse files

Fixed memory leaks when minimizer fails

parent 6231d53f
Hide whitespace changes
Inline Side-by-side
Showing with 262 additions and 235 deletions
libraries/lbfgs/src/lbfgs.cpp
View file @ df6bbb23
...@@ -408,210 +408,231 @@ int lbfgs( ...@@ -408,210 +408,231 @@ int lbfgs(
pf = (lbfgsfloatval_t*)vecalloc(param.past * sizeof(lbfgsfloatval_t)); pf = (lbfgsfloatval_t*)vecalloc(param.past * sizeof(lbfgsfloatval_t));
} }
/* Evaluate the function value and its gradient. */ try {
fx = cd.proc_evaluate(cd.instance, x, g, cd.n, 0); /* Evaluate the function value and its gradient. */
if (0. != param.orthantwise_c) { fx = cd.proc_evaluate(cd.instance, x, g, cd.n, 0);
/* Compute the L1 norm of the variable and add it to the object value. */ if (0. != param.orthantwise_c) {
xnorm = owlqn_x1norm(x, param.orthantwise_start, param.orthantwise_end); /* Compute the L1 norm of the variable and add it to the object value. */
fx += xnorm * param.orthantwise_c; xnorm = owlqn_x1norm(x, param.orthantwise_start, param.orthantwise_end);
owlqn_pseudo_gradient( fx += xnorm * param.orthantwise_c;
pg, x, g, n,
param.orthantwise_c, param.orthantwise_start, param.orthantwise_end
);
}
/* Store the initial value of the objective function. */
if (pf != NULL) {
pf[0] = fx;
}
/*
Compute the direction;
we assume the initial hessian matrix H_0 as the identity matrix.
*/
if (param.orthantwise_c == 0.) {
vecncpy(d, g, n);
} else {
vecncpy(d, pg, n);
}
/*
Make sure that the initial variables are not a minimizer.
*/
vec2norm(&xnorm, x, n);
if (param.orthantwise_c == 0.) {
vec2norm(&gnorm, g, n);
} else {
vec2norm(&gnorm, pg, n);
}
if (xnorm < 1.0) xnorm = 1.0;
if (gnorm / xnorm <= param.epsilon) {
ret = LBFGS_ALREADY_MINIMIZED;
goto lbfgs_exit;
}
/* Compute the initial step:
step = 1.0 / sqrt(vecdot(d, d, n))
*/
vec2norminv(&step, d, n);
k = 1;
end = 0;
for (;;) {
/* Store the current position and gradient vectors. */
veccpy(xp, x, n);
veccpy(gp, g, n);
/* Search for an optimal step. */
if (param.orthantwise_c == 0.) {
ls = linesearch(n, x, &fx, g, d, &step, xp, gp, w, &cd, &param);
} else {
ls = linesearch(n, x, &fx, g, d, &step, xp, pg, w, &cd, &param);
owlqn_pseudo_gradient( owlqn_pseudo_gradient(
pg, x, g, n, pg, x, g, n,
param.orthantwise_c, param.orthantwise_start, param.orthantwise_end param.orthantwise_c, param.orthantwise_start, param.orthantwise_end
); );
} }
if (ls < 0) {
/* Revert to the previous point. */ /* Store the initial value of the objective function. */
veccpy(x, xp, n); if (pf != NULL) {
veccpy(g, gp, n); pf[0] = fx;
ret = ls;
goto lbfgs_exit;
} }
/* Compute x and g norms. */ /*
vec2norm(&xnorm, x, n); Compute the direction;
we assume the initial hessian matrix H_0 as the identity matrix.
*/
if (param.orthantwise_c == 0.) { if (param.orthantwise_c == 0.) {
vec2norm(&gnorm, g, n); vecncpy(d, g, n);
} else { } else {
vec2norm(&gnorm, pg, n); vecncpy(d, pg, n);
}
/* Report the progress. */
if (cd.proc_progress) {
if (ret = cd.proc_progress(cd.instance, x, g, fx, xnorm, gnorm, step, cd.n, k, ls)) {
goto lbfgs_exit;
}
} }
/* /*
Convergence test. Make sure that the initial variables are not a minimizer.
The criterion is given by the following formula:
|g(x)| / \max(1, |x|) < \epsilon
*/ */
vec2norm(&xnorm, x, n);
if (param.orthantwise_c == 0.) {
vec2norm(&gnorm, g, n);
} else {
vec2norm(&gnorm, pg, n);
}
if (xnorm < 1.0) xnorm = 1.0; if (xnorm < 1.0) xnorm = 1.0;
if (gnorm / xnorm <= param.epsilon) { if (gnorm / xnorm <= param.epsilon) {
/* Convergence. */ ret = LBFGS_ALREADY_MINIMIZED;
ret = LBFGS_SUCCESS; goto lbfgs_exit;
break;
} }
/* /* Compute the initial step:
Test for stopping criterion. step = 1.0 / sqrt(vecdot(d, d, n))
The criterion is given by the following formula:
(f(past_x) - f(x)) / f(x) < \delta
*/ */
if (pf != NULL) { vec2norminv(&step, d, n);
/* We don't test the stopping criterion while k < past. */
if (param.past <= k) { k = 1;
/* Compute the relative improvement from the past. */ end = 0;
rate = (pf[k % param.past] - fx) / fx; for (;;) {
/* Store the current position and gradient vectors. */
/* The stopping criterion. */ veccpy(xp, x, n);
if (rate < param.delta) { veccpy(gp, g, n);
ret = LBFGS_STOP;
break; /* Search for an optimal step. */
if (param.orthantwise_c == 0.) {
ls = linesearch(n, x, &fx, g, d, &step, xp, gp, w, &cd, &param);
} else {
ls = linesearch(n, x, &fx, g, d, &step, xp, pg, w, &cd, &param);
owlqn_pseudo_gradient(
pg, x, g, n,
param.orthantwise_c, param.orthantwise_start, param.orthantwise_end
);
}
if (ls < 0) {
/* Revert to the previous point. */
veccpy(x, xp, n);
veccpy(g, gp, n);
ret = ls;
goto lbfgs_exit;
}
/* Compute x and g norms. */
vec2norm(&xnorm, x, n);
if (param.orthantwise_c == 0.) {
vec2norm(&gnorm, g, n);
} else {
vec2norm(&gnorm, pg, n);
}
/* Report the progress. */
if (cd.proc_progress) {
if (ret = cd.proc_progress(cd.instance, x, g, fx, xnorm, gnorm, step, cd.n, k, ls)) {
goto lbfgs_exit;
} }
} }
/* Store the current value of the objective function. */ /*
pf[k % param.past] = fx; Convergence test.
} The criterion is given by the following formula:
|g(x)| / \max(1, |x|) < \epsilon
*/
if (xnorm < 1.0) xnorm = 1.0;
if (gnorm / xnorm <= param.epsilon) {
/* Convergence. */
ret = LBFGS_SUCCESS;
break;
}
if (param.max_iterations != 0 && param.max_iterations < k+1) { /*
/* Maximum number of iterations. */ Test for stopping criterion.
ret = LBFGSERR_MAXIMUMITERATION; The criterion is given by the following formula:
break; (f(past_x) - f(x)) / f(x) < \delta
} */
if (pf != NULL) {
/* We don't test the stopping criterion while k < past. */
if (param.past <= k) {
/* Compute the relative improvement from the past. */
rate = (pf[k % param.past] - fx) / fx;
/* The stopping criterion. */
if (rate < param.delta) {
ret = LBFGS_STOP;
break;
}
}
/* /* Store the current value of the objective function. */
Update vectors s and y: pf[k % param.past] = fx;
s_{k+1} = x_{k+1} - x_{k} = \step * d_{k}. }
y_{k+1} = g_{k+1} - g_{k}.
*/
it = &lm[end];
vecdiff(it->s, x, xp, n);
vecdiff(it->y, g, gp, n);
/* if (param.max_iterations != 0 && param.max_iterations < k+1) {
Compute scalars ys and yy: /* Maximum number of iterations. */
ys = y^t \cdot s = 1 / \rho. ret = LBFGSERR_MAXIMUMITERATION;
yy = y^t \cdot y. break;
Notice that yy is used for scaling the hessian matrix H_0 (Cholesky factor). }
*/
vecdot(&ys, it->y, it->s, n);
vecdot(&yy, it->y, it->y, n);
it->ys = ys;
/* /*
Recursive formula to compute dir = -(H \cdot g). Update vectors s and y:
This is described in page 779 of: s_{k+1} = x_{k+1} - x_{k} = \step * d_{k}.
Jorge Nocedal. y_{k+1} = g_{k+1} - g_{k}.
Updating Quasi-Newton Matrices with Limited Storage. */
Mathematics of Computation, Vol. 35, No. 151, it = &lm[end];
pp. 773--782, 1980. vecdiff(it->s, x, xp, n);
*/ vecdiff(it->y, g, gp, n);
bound = (m <= k) ? m : k;
++k;
end = (end + 1) % m;
/* Compute the steepest direction. */ /*
if (param.orthantwise_c == 0.) { Compute scalars ys and yy:
/* Compute the negative of gradients. */ ys = y^t \cdot s = 1 / \rho.
vecncpy(d, g, n); yy = y^t \cdot y.
} else { Notice that yy is used for scaling the hessian matrix H_0 (Cholesky factor).
vecncpy(d, pg, n); */
} vecdot(&ys, it->y, it->s, n);
vecdot(&yy, it->y, it->y, n);
it->ys = ys;
j = end; /*
for (i = 0;i < bound;++i) { Recursive formula to compute dir = -(H \cdot g).
j = (j + m - 1) % m; /* if (--j == -1) j = m-1; */ This is described in page 779 of:
it = &lm[j]; Jorge Nocedal.
/* \alpha_{j} = \rho_{j} s^{t}_{j} \cdot q_{k+1}. */ Updating Quasi-Newton Matrices with Limited Storage.
vecdot(&it->alpha, it->s, d, n); Mathematics of Computation, Vol. 35, No. 151,
it->alpha /= it->ys; pp. 773--782, 1980.
/* q_{i} = q_{i+1} - \alpha_{i} y_{i}. */ */
vecadd(d, it->y, -it->alpha, n); bound = (m <= k) ? m : k;
} ++k;
end = (end + 1) % m;
/* Compute the steepest direction. */
if (param.orthantwise_c == 0.) {
/* Compute the negative of gradients. */
vecncpy(d, g, n);
} else {
vecncpy(d, pg, n);
}
vecscale(d, ys / yy, n); j = end;
for (i = 0;i < bound;++i) {
j = (j + m - 1) % m; /* if (--j == -1) j = m-1; */
it = &lm[j];
/* \alpha_{j} = \rho_{j} s^{t}_{j} \cdot q_{k+1}. */
vecdot(&it->alpha, it->s, d, n);
it->alpha /= it->ys;
/* q_{i} = q_{i+1} - \alpha_{i} y_{i}. */
vecadd(d, it->y, -it->alpha, n);
}
for (i = 0;i < bound;++i) { vecscale(d, ys / yy, n);
it = &lm[j];
/* \beta_{j} = \rho_{j} y^t_{j} \cdot \gamma_{i}. */
vecdot(&beta, it->y, d, n);
beta /= it->ys;
/* \gamma_{i+1} = \gamma_{i} + (\alpha_{j} - \beta_{j}) s_{j}. */
vecadd(d, it->s, it->alpha - beta, n);
j = (j + 1) % m; /* if (++j == m) j = 0; */
}
/* for (i = 0;i < bound;++i) {
Constrain the search direction for orthant-wise updates. it = &lm[j];
*/ /* \beta_{j} = \rho_{j} y^t_{j} \cdot \gamma_{i}. */
if (param.orthantwise_c != 0.) { vecdot(&beta, it->y, d, n);
for (i = param.orthantwise_start;i < param.orthantwise_end;++i) { beta /= it->ys;
if (d[i] * pg[i] >= 0) { /* \gamma_{i+1} = \gamma_{i} + (\alpha_{j} - \beta_{j}) s_{j}. */
d[i] = 0; vecadd(d, it->s, it->alpha - beta, n);
j = (j + 1) % m; /* if (++j == m) j = 0; */
}
/*
Constrain the search direction for orthant-wise updates.
*/
if (param.orthantwise_c != 0.) {
for (i = param.orthantwise_start;i < param.orthantwise_end;++i) {
if (d[i] * pg[i] >= 0) {
d[i] = 0;
}
} }
} }
}
/* /*
Now the search direction d is ready. We try step = 1 first. Now the search direction d is ready. We try step = 1 first.
*/ */
step = 1.0; step = 1.0;
}
}
catch (...) {
vecfree(pf);
/* Free memory blocks used by this function. */
if (lm != NULL) {
for (i = 0;i < m;++i) {
vecfree(lm[i].s);
vecfree(lm[i].y);
}
vecfree(lm);
}
vecfree(pg);
vecfree(w);
vecfree(d);
vecfree(gp);
vecfree(g);
vecfree(xp);
throw;
} }
lbfgs_exit: lbfgs_exit:
......
openmmapi/src/LocalEnergyMinimizer.cpp
View file @ df6bbb23
...@@ -105,83 +105,89 @@ static lbfgsfloatval_t evaluate(void *instance, const lbfgsfloatval_t *x, lbfgsf ...@@ -105,83 +105,89 @@ static lbfgsfloatval_t evaluate(void *instance, const lbfgsfloatval_t *x, lbfgsf
void LocalEnergyMinimizer::minimize(Context& context, double tolerance, int maxIterations) { void LocalEnergyMinimizer::minimize(Context& context, double tolerance, int maxIterations) {
const System& system = context.getSystem(); const System& system = context.getSystem();
int numParticles = system.getNumParticles(); int numParticles = system.getNumParticles();
lbfgsfloatval_t *x = lbfgs_malloc(numParticles*3);
if (x == NULL)
throw OpenMMException("LocalEnergyMinimizer: Failed to allocate memory");
double constraintTol = context.getIntegrator().getConstraintTolerance(); double constraintTol = context.getIntegrator().getConstraintTolerance();
double workingConstraintTol = std::max(1e-4, constraintTol); double workingConstraintTol = std::max(1e-4, constraintTol);
double k = tolerance/workingConstraintTol; double k = tolerance/workingConstraintTol;
lbfgsfloatval_t *x = lbfgs_malloc(numParticles*3);
if (x == NULL)
throw OpenMMException("LocalEnergyMinimizer: Failed to allocate memory");
try {
// Initialize the minimizer. // Initialize the minimizer.
lbfgs_parameter_t param; lbfgs_parameter_t param;
lbfgs_parameter_init(&param); lbfgs_parameter_init(&param);
if (!context.getPlatform().supportsDoublePrecision()) if (!context.getPlatform().supportsDoublePrecision())
param.xtol = 1e-7; param.xtol = 1e-7;
param.max_iterations = maxIterations; param.max_iterations = maxIterations;
param.linesearch = LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE; param.linesearch = LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE;
// Make sure the initial configuration satisfies all constraints. // Make sure the initial configuration satisfies all constraints.
context.applyConstraints(workingConstraintTol); context.applyConstraints(workingConstraintTol);
// Record the initial positions and determine a normalization constant for scaling the tolerance. // Record the initial positions and determine a normalization constant for scaling the tolerance.
vector<Vec3> initialPos = context.getState(State::Positions).getPositions(); vector<Vec3> initialPos = context.getState(State::Positions).getPositions();
double norm = 0.0; double norm = 0.0;
for (int i = 0; i < numParticles; i++) { for (int i = 0; i < numParticles; i++) {
x[3*i] = initialPos[i][0]; x[3*i] = initialPos[i][0];
x[3*i+1] = initialPos[i][1]; x[3*i+1] = initialPos[i][1];
x[3*i+2] = initialPos[i][2]; x[3*i+2] = initialPos[i][2];
norm += initialPos[i].dot(initialPos[i]); norm += initialPos[i].dot(initialPos[i]);
}
norm /= numParticles;
norm = (norm < 1 ? 1 : sqrt(norm));
param.epsilon = tolerance/norm;
// Repeatedly minimize, steadily increasing the strength of the springs until all constraints are satisfied.
double prevMaxError = 1e10;
while (true) {
// Perform the minimization.
lbfgsfloatval_t fx;
MinimizerData data(context, k);
lbfgs(numParticles*3, x, &fx, evaluate, NULL, &data, &param);
// Check whether all constraints are satisfied.
vector<Vec3> positions = context.getState(State::Positions).getPositions();
int numConstraints = system.getNumConstraints();
double maxError = 0.0;
for (int i = 0; i < numConstraints; i++) {
int particle1, particle2;
double distance;
system.getConstraintParameters(i, particle1, particle2, distance);
Vec3 delta = positions[particle2]-positions[particle1];
double r = sqrt(delta.dot(delta));
double error = fabs(r-distance);
if (error > maxError)
maxError = error;
} }
if (maxError <= workingConstraintTol) norm /= numParticles;
break; // All constraints are satisfied. norm = (norm < 1 ? 1 : sqrt(norm));
context.setPositions(initialPos); param.epsilon = tolerance/norm;
if (maxError >= prevMaxError)
break; // Further tightening the springs doesn't seem to be helping, so just give up. // Repeatedly minimize, steadily increasing the strength of the springs until all constraints are satisfied.
prevMaxError = maxError;
k *= 10; double prevMaxError = 1e10;
if (maxError > 100*workingConstraintTol) { while (true) {
// We've gotten far enough from a valid state that we might have trouble getting // Perform the minimization.
// back, so reset to the original positions.
lbfgsfloatval_t fx;
for (int i = 0; i < numParticles; i++) { MinimizerData data(context, k);
x[3*i] = initialPos[i][0]; lbfgs(numParticles*3, x, &fx, evaluate, NULL, &data, &param);
x[3*i+1] = initialPos[i][1];
x[3*i+2] = initialPos[i][2]; // Check whether all constraints are satisfied.
vector<Vec3> positions = context.getState(State::Positions).getPositions();
int numConstraints = system.getNumConstraints();
double maxError = 0.0;
for (int i = 0; i < numConstraints; i++) {
int particle1, particle2;
double distance;
system.getConstraintParameters(i, particle1, particle2, distance);
Vec3 delta = positions[particle2]-positions[particle1];
double r = sqrt(delta.dot(delta));
double error = fabs(r-distance);
if (error > maxError)
maxError = error;
}
if (maxError <= workingConstraintTol)
break; // All constraints are satisfied.
context.setPositions(initialPos);
if (maxError >= prevMaxError)
break; // Further tightening the springs doesn't seem to be helping, so just give up.
prevMaxError = maxError;
k *= 10;
if (maxError > 100*workingConstraintTol) {
// We've gotten far enough from a valid state that we might have trouble getting
// back, so reset to the original positions.
for (int i = 0; i < numParticles; i++) {
x[3*i] = initialPos[i][0];
x[3*i+1] = initialPos[i][1];
x[3*i+2] = initialPos[i][2];
}
} }
} }
} }
catch (...) {
lbfgs_free(x);
throw;
}
lbfgs_free(x); lbfgs_free(x);
// If necessary, do a final constraint projection to make sure they are satisfied // If necessary, do a final constraint projection to make sure they are satisfied
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment