52void 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++) {
69 sum += dir[i] * dir[i];
75 for (
unsigned int i = 0; i < dim; i++) {
76 dir[i] *= maxStep / sum;
83 for (
unsigned int i = 0; i < dim; i++) {
84 slope += dir[i] * grad[i];
91 for (
unsigned int i = 0; i < dim; i++) {
92 double temp = fabs(dir[i]) / std::max(fabs(oldPt[i]), 1.0);
101 while (it < MAX_ITER_LINEAR_SEARCH) {
103 if (lambda < lambdaMin) {
108 for (
unsigned int i = 0; i < dim; i++) {
109 newPt[i] = oldPt[i] + lambda * dir[i];
111 newVal = func(newPt);
113 if (newVal - oldVal <=
FUNCTOL * lambda * slope) {
121 tmpLambda = -slope / (2.0 * (newVal - oldVal - slope));
123 double rhs1 = newVal - oldVal - lambda * slope;
124 double rhs2 = val2 - oldVal - lambda2 * slope;
125 double a = (rhs1 / (lambda * lambda) - rhs2 / (lambda2 * lambda2)) /
127 double b = (-lambda2 * rhs1 / (lambda * lambda) +
128 lambda * rhs2 / (lambda2 * lambda2)) /
131 tmpLambda = -slope / (2.0 * b);
133 double disc = b * b - 3 * a * slope;
135 tmpLambda = 0.5 * lambda;
136 }
else if (b <= 0.0) {
137 tmpLambda = (-b + sqrt(disc)) / (3.0 * a);
139 tmpLambda = -slope / (b + sqrt(disc));
142 if (tmpLambda > 0.5 * lambda) {
143 tmpLambda = 0.5 * lambda;
148 lambda = std::max(tmpLambda, 0.1 * lambda);
153 for (
unsigned int i = 0; i < dim; i++) {
184int minimize(
unsigned int dim,
double *pos,
double gradTol,
185 unsigned int &numIters,
double &funcVal, EnergyFunctor func,
186 GradientFunctor gradFunc,
unsigned int snapshotFreq,
188 unsigned int maxIts =
MAXITS) {
193 std::vector<double> grad(dim);
194 std::vector<double> dGrad(dim);
195 std::vector<double> hessDGrad(dim);
196 std::vector<double> xi(dim);
197 std::vector<double> invHessian(dim * dim, 0);
198 std::unique_ptr<double[]> newPos(
new double[dim]);
199 snapshotFreq = std::min(snapshotFreq, maxIts);
202 double fp = func(pos);
203 gradFunc(pos, grad.data());
206 for (
unsigned int i = 0; i < dim; i++) {
207 unsigned int itab = i * dim;
209 invHessian[itab + i] = 1.0;
212 sum += pos[i] * pos[i];
215 double maxStep =
MAXSTEP * std::max(sqrt(sum),
static_cast<double>(dim));
217 for (
unsigned int iter = 1; iter <= maxIts; ++iter) {
222 linearSearch(dim, pos, fp, grad.data(), xi.data(), newPos.get(), funcVal, func,
231 for (
unsigned int i = 0; i < dim; i++) {
232 xi[i] = newPos[i] - pos[i];
234 double temp = fabs(xi[i]) / std::max(fabs(pos[i]), 1.0);
243 if (snapshotVect && snapshotFreq) {
245 snapshotVect->push_back(s);
251 double gradScale = gradFunc(pos, grad.data());
255 double term = std::max(funcVal * gradScale, 1.0);
256 for (
unsigned int i = 0; i < dim; i++) {
257 double temp = fabs(grad[i]) * std::max(fabs(pos[i]), 1.0);
258 test = std::max(test, temp);
259 dGrad[i] = grad[i] - dGrad[i];
264 if (test < gradTol) {
265 if (snapshotVect && snapshotFreq) {
267 snapshotVect->push_back(s);
277 double fac = 0, fae = 0, sumDGrad = 0, sumXi = 0;
278 for (
unsigned int i = 0; i < dim; i++) {
279 double *ivh = &(invHessian[i * dim]);
280 double &hdgradi = hessDGrad[i];
281 double *dgj = dGrad.data();
283 for (
unsigned int j = 0; j < dim; ++j, ++ivh, ++dgj) {
284 hdgradi += *ivh * *dgj;
286 fac += dGrad[i] * xi[i];
287 fae += dGrad[i] * hessDGrad[i];
288 sumDGrad += dGrad[i] * dGrad[i];
289 sumXi += xi[i] * xi[i];
291 if (fac > sqrt(
EPS * sumDGrad * sumXi)) {
293 double fad = 1.0 / fae;
294 for (
unsigned int i = 0; i < dim; i++) {
295 dGrad[i] = fac * xi[i] - fad * hessDGrad[i];
298 for (
unsigned int i = 0; i < dim; i++) {
299 unsigned int itab = i * dim;
300 double pxi = fac * xi[i], hdgi = fad * hessDGrad[i],
301 dgi = fae * dGrad[i];
302 double *pxj = &(xi[i]), *hdgj = &(hessDGrad[i]), *dgj = &(dGrad[i]);
303 for (
unsigned int j = i; j < dim; ++j, ++pxj, ++hdgj, ++dgj) {
304 invHessian[itab + j] += pxi * *pxj - hdgi * *hdgj + dgi * *dgj;
305 invHessian[j * dim + i] = invHessian[itab + j];
310 for (
unsigned int i = 0; i < dim; i++) {
311 unsigned int itab = i * dim;
314 double *ivh = &(invHessian[itab]);
315 double *gj = grad.data();
316 for (
unsigned int j = 0; j < dim; ++j, ++ivh, ++gj) {
320 if (snapshotVect && snapshotFreq && !(iter % snapshotFreq)) {
322 snapshotVect->push_back(s);
323 newPos.reset(
new double[dim]);