fix packaging (into central file) and copyright dates.
[CommonLispStat.git] / lib / minimize.c
blob05f5028b9ce25f3d9d9711c8eab0a77711c925ae
1 /* minimize - minimization for Bayes routines in XLISP-STAT and S */
2 /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney */
3 /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz */
4 /* You may give out copies of this software; for conditions see the */
5 /* file COPYING included with this distribution. */
7 /*
8 * Nonlinear optimization modules adapted from Dennis and Schnabel,
9 * "Numerical Methods for Unconstrained Optimization and Nonlinear
10 * Equations."
13 #include <stdio.h>
14 #include "xmath.h"
15 extern char buf[200];
16 #define PRINTSTR(s) printf(s)
18 typedef struct {
19 size_t n, k;
20 void (*ffun)(), (*gfun)();
21 double f, typf, new_f;
22 double crit, new_crit;
23 double *x, *new_x, *sx, *delf, *new_delf, *qnstep, *F;
24 double **hessf, **H, **L, **new_delg;
25 double gradtol, steptol, maxstep;
26 int itncount, itnlimit, maxtaken, consecmax, retcode, termcode;
27 int use_line_search, verbose, values_supplied, change_sign;
28 double diagadd;
29 } Iteration;
32 extern double macheps();
35 extern double choldecomp();
36 static void eval_gradient(Iteration *);
37 static void eval_next_gradient(Iteration *);
40 /*
41 This is a little more advanced section on functions, but is very useful. Take this for example:
42 int applyeqn(int F(int), int max, int min) {
43 int itmp;
45 itmp = F(int) + min;
46 itmp = itmp - max;
47 return itmp;
50 What does this function do if we call it with applyeqn(square(x),
51 y, z);? What happens is that the int F(int) is a reference to the
52 function that is passed in as a parameter. Thus inside applyeqn where
53 there is a call to F, it actually is a call to square! This is very
54 useful if we have one set function, but wish to vary the input
55 according to a particular function. So if we had a different function
56 called cube we could change how we call applyeqn by calling the
57 function by applyeqn(cube(x), y, z);. */
60 /************************************************************************/
61 /** **/
62 /** Definitions and Globals **/
63 /** **/
64 /************************************************************************/
66 # define nil 0L
67 # define FALSE 0
68 # define TRUE 1
70 # define INIT_GRAD_FRAC .001
71 # define CONSEC_MAX_LIMIT 5
72 # define ALPHA .0001
73 # define MAX_STEP_FACTOR 1000
74 # define GRADTOL_POWER 1.0 / 3.0
75 # define STEPTOL_POWER 2.0 / 3.0
76 # define ITNLIMIT 100
77 # define VERBOSE 0
78 # define USE_SEARCH TRUE
80 typedef double **RMatrix, *RVector;
83 static char *termcodes[] = {"not yet terminated",
84 "gradient size is less than gradient tolerance",
85 "step size is less than step tolerance",
86 "no satisfactory step found in backtracking",
87 "iteration limit exceeded",
88 "maximum size step taken 5 iterations in a row"};
90 /************************************************************************/
91 /** **/
92 /** Utility Functions **/
93 /** **/
94 /************************************************************************/
96 static double Max(a, b)
97 double a, b;
99 return(a > b ? a : b);
102 static double Min(a, b)
103 double a, b;
105 return(a > b ? b : a);
108 /************************************************************************/
109 /** **/
110 /** Cholesky Solving Functions **/
111 /** **/
112 /************************************************************************/
115 /* solve Ly = b for y */
116 static void
117 lsolve(size_t n, double *b, double **L, double *y)
119 size_t i, j;
121 for (i = 0; i < n; i++) {
122 y[i] = b[i];
123 for (j = 0; j < i; j++) y[i] -= L[i][j] * y[j];
124 if (L[i][i] != 0) y[i] /= L[i][i];
128 /* solve (L^T)x = y for x */
129 static void
130 ltsolve(size_t n, double *y, double **L, double *x)
132 size_t i, j;
134 for (i = n - 1; i >= 0; i--) {
135 x[i] = y[i];
136 for (j = i + 1; j < n; j++) x[i] -= L[j][i] * x[j];
137 if (L[i][i] != 0) x[i] /= L[i][i];
141 /* solve (L L^T) s = -g for s */
142 static void
143 cholsolve(size_t n,
144 double *g, double **L,
145 double *s)
147 size_t i;
149 /* solve Ly = g */
150 lsolve(n, g, L, s);
152 /* solve L^Ts = y */
153 ltsolve(n, s, L, s);
155 for (i = 0; i < n; i++) s[i] = -s[i];
158 static void
159 modelhess(size_t n, double *sx, double **H, double **L, double *diagadd)
161 size_t i, j;
162 double sqrteps, maxdiag, mindiag, maxoff, maxoffl, maxposdiag, mu,
163 maxadd, maxev, minev, offrow, sdd;
165 /* scale H on both sides with sx */
166 for (i = 0; i < n; i++)
167 for (j = 0; j < n; j++) H[i][j] /= sx[i] * sx[j];
170 * find mu to add to diagonal so no diagonal elements are negative
171 * and largest diagonal dominates largest off diagonal element in H
173 sqrteps = sqrt(macheps());
174 for (maxdiag = H[0][0], i = 1; i < n; i++)
175 maxdiag = Max(maxdiag, H[i][i]);
176 for (mindiag = H[0][0], i = 1; i < n; i++)
177 mindiag = Min(mindiag, H[i][i]);
178 maxposdiag = Max(0.0, maxdiag);
180 if (mindiag <= sqrteps * maxposdiag) {
181 mu = 2 * (maxposdiag - mindiag) * sqrteps - mindiag;
182 maxdiag += mu;
184 else mu = 0.0;
186 if (n > 1) {
187 for (maxoff = fabs(H[0][1]), i = 0; i < n; i++)
188 for (j = i + 1; j < n; j++)
189 maxoff = Max(maxoff, fabs(H[i][j]));
191 else maxoff = 0.0;
193 if (maxoff * (1 + 2 * sqrteps) > maxdiag) {
194 mu += (maxoff - maxdiag) + 2 * sqrteps * maxoff;
195 maxdiag = maxoff * (1 + 2 * sqrteps);
198 if (maxdiag == 0.0) {
199 mu = 1;
200 maxdiag = 1;
203 if (mu > 0) for (i = 0; i < n; i++) H[i][i] += mu;
205 maxoffl = sqrt(Max(maxdiag, maxoff / n));
208 * compute the perturbed Cholesky decomposition
210 for (i = 0; i < n; i++)
211 for (j = 0; j < n; j++) L[i][j] = H[i][j];
212 choldecomp(L, n, maxoffl, &maxadd);
215 * add something to diagonal, if needed, to make positive definite
216 * and recompute factorization
218 if (maxadd > 0) {
219 maxev = H[0][0];
220 minev = H[0][0];
221 for (i = 0; i < n; i++) {
222 for (offrow = 0.0, j = 0; j < n; j++)
223 if (i != j) offrow += fabs(H[i][j]);
224 maxev = Max(maxev, H[i][i] + offrow);
225 minev = Min(minev, H[i][i] - offrow);
227 sdd = (maxev - minev) * sqrteps - minev;
228 sdd = Max(sdd, 0.0);
229 mu = Min(maxadd, sdd);
230 for (i = 0; i < n; i++) H[i][i] += mu;
231 for (i = 0; i < n; i++)
232 for (j = 0; j < n; j++) L[i][j] = H[i][j];
233 choldecomp(L, n, maxoffl, &maxadd);
234 *diagadd = mu;
236 else *diagadd = 0.0;
238 /* unscale H and L */
239 for (i = 0; i < n; i++)
240 for (j = 0; j < n; j++) {
241 H[i][j] *= sx[i] * sx[j];
242 L[i][j] *= sx[i];
246 /************************************************************************/
247 /** **/
248 /** Stopping Criteria **/
249 /** **/
250 /************************************************************************/
252 static double
253 gradsize(Iteration *iter, int new)
255 size_t n, i;
256 double size, term, crit, typf;
257 double *x, *delf, *sx;
259 n = iter->n + iter->k;
260 crit = iter->crit;
261 typf = iter->typf;
262 x = iter->x;
263 sx = iter->sx;
264 delf = (new) ? iter->new_delf : iter->delf;
266 for (i = 0, size = 0.0; i < n; i++) {
267 term = fabs(delf[i]) * Max(x[i], 1.0 / sx[i]) / Max(fabs(crit), typf);
268 size = Max(size, term);
270 return(size);
273 static double
274 incrsize(Iteration *iter)
276 size_t n, i;
277 double size, term;
278 double *x, *new_x, *sx;
280 new_x = iter->new_x;
281 n = iter->n + iter->k;
282 x = iter->x;
283 sx = iter->sx;
285 for (i = 0, size = 0.0; i < n; i++) {
286 term = fabs(x[i] - new_x[i]) / Max(fabs(x[i]), 1.0 / sx[i]);
287 size = Max(size, term);
289 return(size);
292 static int /* AJR:TESTME guess! */
293 stoptest0(Iteration *iter)
295 iter->consecmax = 0;
297 if (gradsize(iter, FALSE) <= INIT_GRAD_FRAC * iter->gradtol)
298 iter->termcode = 1;
299 else iter->termcode = 0;
301 return(iter->termcode);
304 static int
305 stoptest(Iteration *iter)
307 int termcode, retcode, itncount, itnlimit, maxtaken, consecmax;
308 double gradtol, steptol;
310 retcode = iter->retcode;
311 gradtol = iter->gradtol;
312 steptol = iter->steptol;
313 itncount = iter->itncount;
314 itnlimit = iter->itnlimit;
315 maxtaken = iter->maxtaken;
316 consecmax = iter->consecmax;
318 termcode = 0;
319 if (retcode == 1) termcode = 3;
320 else if (gradsize(iter, TRUE) <= gradtol) termcode = 1;
321 else if (incrsize(iter) <= steptol) termcode = 2;
322 else if (itncount >= itnlimit) termcode = 4;
323 else if (maxtaken) {
324 consecmax++;
325 if (consecmax >= CONSEC_MAX_LIMIT) termcode = 5;
327 else consecmax = 0;
329 iter->consecmax = consecmax;
330 iter->termcode = termcode;
332 return(termcode);
335 /************************************************************************/
336 /** **/
337 /** Function and Derivative Evaluation **/
338 /** **/
339 /************************************************************************/
341 static void
342 eval_funval(Iteration *iter)
344 int i;
346 (*(iter->ffun))(iter->x, &iter->f, nil, nil);
347 if (iter->k == 0) iter->crit = iter->f;
348 else {
349 eval_gradient(iter);
350 for (i = 0, iter->crit = 0.0; i < iter->n + iter->k; i++)
351 iter->crit += 0.5 * pow(fabs(iter->delf[i]), 2.0);
355 static void
356 eval_next_funval(Iteration *iter)
358 int i;
360 (*(iter->ffun))(iter->new_x, &iter->new_f, nil, nil);
361 if (iter->k == 0) iter->new_crit = iter->new_f;
362 else {
363 eval_next_gradient(iter);
364 for (i = 0, iter->new_crit = 0.0; i < iter->n + iter->k; i++)
365 iter->new_crit += 0.5 * pow(fabs(iter->new_delf[i]), 2.0);
369 static void
370 eval_gradient(Iteration *iter)
372 size_t i, j, n, k;
374 n = iter->n;
375 k = iter->k;
377 (*(iter->ffun))(iter->x, nil, iter->delf, nil);
378 if (k > 0) {
379 (*(iter->gfun))(iter->x, iter->delf + n, nil, nil);
380 (*(iter->gfun))(iter->x, nil, iter->new_delg, nil);
381 for (i = 0; i < n; i++) {
382 for (j = 0; j < k; j++)
383 iter->delf[i] += iter->x[n + j] * iter->new_delg[j][i];
384 for (j = 0; j < k; j++) {
385 iter->hessf[n + j][i] = iter->new_delg[j][i];
386 iter->hessf[i][n + j] = iter->new_delg[j][i];
392 static void
393 eval_next_gradient(Iteration *iter)
395 size_t i, j, n, k;
397 n = iter->n;
398 k = iter->k;
399 (*(iter->ffun))(iter->new_x, nil, iter->new_delf, nil);
400 if (k > 0) {
401 (*(iter->gfun))(iter->new_x, iter->new_delf + n, nil, nil);
402 (*(iter->gfun))(iter->new_x, nil, iter->new_delg, nil);
403 for (i = 0; i < n; i++) {
404 for (j = 0; j < k; j++)
405 iter->new_delf[i] += iter->new_x[n + j] * iter->new_delg[j][i];
410 static void eval_hessian(Iteration *iter)
412 size_t i, j, n, k;
414 n = iter->n;
415 k = iter->k;
416 (*(iter->ffun))(iter->x, nil, nil, iter->hessf);
417 for (i = n; i < n + k; i++)
418 for (j = n; j < n + k; j++) iter->hessf[i][j] = 0.0;
421 /************************************************************************/
422 /** **/
423 /** Backtracking Line Search **/
424 /** **/
425 /************************************************************************/
427 static void
428 linesearch(Iteration *iter)
430 size_t i, n;
431 double newtlen, maxstep, initslope, rellength, lambda, minlambda,
432 lambdatemp, lambdaprev, a, b, disc, critprev, f1, f2, a11, a12, a21, a22,
433 del;
434 double *qnstep, *delf, *x, *new_x, *sx;
436 n = iter->n + iter->k;
437 if (! iter->use_line_search) {
438 iter->maxtaken = FALSE;
439 for (i = 0; i < n; i++)
440 iter->new_x[i] = iter->x[i] + iter->qnstep[i];
441 eval_next_funval(iter);
442 iter->retcode = 0;
444 else{
445 qnstep = iter->qnstep;
446 maxstep = iter->maxstep;
447 delf = (iter->k == 0) ? iter->delf : iter->F;
448 x = iter->x;
449 new_x = iter->new_x;
450 sx = iter->sx;
452 iter->maxtaken = FALSE;
453 iter->retcode = 2;
455 for (i = 0, newtlen = 0.0; i < n; i++)
456 newtlen += pow(sx[i] * qnstep[i], 2.0);
457 newtlen = sqrt(newtlen);
459 if (newtlen > maxstep) {
460 for (i = 0; i < n; i++) qnstep[i] *= (maxstep / newtlen);
461 newtlen = maxstep;
464 for (i = 0, initslope = 0.0; i < n; i++) initslope += delf[i] * qnstep[i];
466 for (i = 0, rellength = 0.0; i < n; i++)
467 rellength = Max(rellength, fabs(qnstep[i]) / Max(fabs(x[i]), 1.0 / sx[i]));
469 minlambda = (rellength == 0.0) ? 2.0 : iter->steptol / rellength;
471 lambda = 1.0;
472 lambdaprev = 1.0; /* to make compilers happy */
473 critprev = 1.0; /* to make compilers happy */
474 while (iter->retcode > 1) {
475 for (i = 0; i < n; i++) new_x[i] = x[i] + lambda * qnstep[i];
476 eval_next_funval(iter);
477 if (iter->new_crit <= iter->crit + ALPHA * lambda * initslope) {
478 iter->retcode = 0;
479 if (lambda == 1.0 && newtlen > 0.99 * maxstep) iter->maxtaken = TRUE;
481 else if (lambda < minlambda) {
482 iter->retcode = 1;
483 iter->new_crit = iter->crit;
484 for (i = 0; i < n; i++) new_x[i] = x[i];
486 else {
487 if (lambda == 1.0) { /* first backtrack, quadratic fit */
488 lambdatemp = - initslope
489 / (2 * (iter->new_crit - iter->crit - initslope));
491 else { /* all subsequent backtracks, cubic fit */
492 del = lambda - lambdaprev;
493 f1 = iter->new_crit - iter->crit - lambda * initslope;
494 f2 = critprev - iter->crit - lambdaprev * initslope;
495 a11 = 1.0 / (lambda * lambda);
496 a12 = -1.0 / (lambdaprev * lambdaprev);
497 a21 = -lambdaprev * a11;
498 a22 = -lambda * a12;
499 a = (a11 * f1 + a12 * f2) / del;
500 b = (a21 * f1 + a22 * f2) / del;
501 disc = b * b - 3 * a * initslope;
502 if (a == 0) { /* cubic is a quadratic */
503 lambdatemp = - initslope / (2 * b);
505 else { /* legitimate cubic */
506 lambdatemp = (-b + sqrt(disc)) / (3 * a);
508 lambdatemp = Min(lambdatemp, 0.5 * lambda);
510 lambdaprev = lambda;
511 critprev = iter->new_crit;
512 lambda = Max(0.1 * lambda, lambdatemp);
513 if (iter->verbose > 0) {
514 sprintf(buf, "Backtracking: lambda = %g\n", lambda);
515 PRINTSTR(buf);
522 /************************************************************************/
523 /** **/
524 /** Status Printing Functions **/
525 /** **/
526 /************************************************************************/
528 static void
529 print_header(Iteration *iter)
531 if (iter->verbose > 0) {
532 sprintf(buf, "Iteration %d.\n", iter->itncount);
533 PRINTSTR(buf);
537 static void
538 print_status(Iteration *iter)
540 size_t i, j;
542 if (iter->verbose > 0) {
543 sprintf(buf, "Criterion value = %g\n",
544 (iter->change_sign) ? -iter->crit : iter->crit);
545 PRINTSTR(buf);
546 if (iter->verbose > 1) {
547 PRINTSTR("Location = <");
548 for (i = 0; i < iter->n + iter->k; i++) {
549 sprintf(buf, (i < iter->n + iter->k - 1) ? "%g " : "%g>\n", iter->x[i]);
550 PRINTSTR(buf);
553 if (iter->verbose > 2) {
554 PRINTSTR("Gradient = <");
555 for (i = 0; i < iter->n + iter->k; i++) {
556 sprintf(buf, (i < iter->n + iter->k - 1) ? "%g " : "%g>\n",
557 (iter->change_sign) ? -iter->delf[i] : iter->delf[i]);
558 PRINTSTR(buf);
561 if (iter->verbose > 3) {
562 PRINTSTR("Hessian:\n");
563 for (i = 0; i < iter->n + iter->k; i++) {
564 for (j = 0; j < iter->n + iter->k; j++) {
565 sprintf(buf, (j < iter->n + iter->k - 1) ? "%g " : "%g\n",
566 (iter->change_sign) ? -iter->hessf[i][j] : iter->hessf[i][j]);
567 PRINTSTR(buf);
572 if (iter->termcode != 0 && iter->verbose > 0) {
573 sprintf(buf, "Reason for termination: %s.\n", termcodes[iter->termcode]);
574 PRINTSTR(buf);
578 /************************************************************************/
579 /** **/
580 /** Iteration Driver **/
581 /** **/
582 /************************************************************************/
584 static void
585 findqnstep(Iteration *iter)
587 size_t i, j, N, l;
589 if (iter->k == 0) {
590 modelhess(iter->n, iter->sx, iter->hessf, iter->L, &iter->diagadd);
591 cholsolve(iter->n, iter->delf, iter->L, iter->qnstep);
593 else {
594 N = iter->n + iter->k;
595 for (i = 0; i < N; i++) {
596 for (l = 0, iter->F[i] = 0.0; l < N; l++)
597 iter->F[i] += iter->hessf[i][l] * iter->delf[l];
598 for (j = 0; j < N; j++)
599 for (l = 0, iter->H[i][j] = 0.0; l < N; l++)
600 iter->H[i][j] += iter->hessf[i][l] * iter->hessf[j][l];
602 modelhess(N, iter->sx, iter->H, iter->L, &iter->diagadd);
603 cholsolve(N, iter->F, iter->L, iter->qnstep);
607 static void iterupdate(Iteration *iter)
609 size_t i, j, n, k;
611 n = iter->n;
612 k = iter->k;
613 iter->f = iter->new_f;
614 iter->crit = iter->new_crit;
615 for (i = 0; i < n + k; i++) {
616 iter->x[i] = iter->new_x[i];
617 iter->delf[i] = iter->new_delf[i];
619 for (i = 0; i < k; i++) {
620 for (j = 0; j < n; j++) {
621 iter->hessf[n + i][j] = iter->new_delg[i][j];
622 iter->hessf[j][n + i] = iter->new_delg[i][j];
627 static void
628 mindriver(Iteration *iter)
630 iter->consecmax = 0;
631 iter->itncount = 0;
632 iter->termcode = 0;
633 if (! iter->values_supplied) {
634 eval_funval(iter);
635 if (iter->k == 0) eval_gradient(iter);
636 eval_hessian(iter);
639 stoptest0(iter);
640 print_header(iter);
641 print_status(iter);
642 while (iter->termcode == 0) {
643 iter->itncount++;
644 print_header(iter);
645 findqnstep(iter);
646 linesearch(iter);
647 if (iter->k == 0) eval_next_gradient(iter);
648 stoptest(iter);
649 iterupdate(iter);
650 eval_hessian(iter);
651 print_status(iter);
655 /************************************************************************/
656 /** **/
657 /** External Interface Routines **/
658 /** **/
659 /************************************************************************/
661 static Iteration myiter;
663 size_t
664 minworkspacesize(size_t n, size_t k)
666 size_t size;
668 /* x, new_x, sx, delf, new_delf, qnstep and F */
669 size = 7 * sizeof(double) * (n + k);
671 /* hessf, H and L */
672 size += 3 * (sizeof(double *) * (n + k)
673 + sizeof(double) * (n + k) * (n + k));
675 /* delg and new_delg */
676 size += 2 * (sizeof(double *) * k + sizeof(double) * n * k);
678 return(size);
681 char *
682 minresultstring(int code)
684 if (code <= 0) return("bad input data");
685 else if (code <= 5) return(termcodes[code]);
686 else return("unknown return code");
689 void
690 minsetup(size_t n, size_t k,
691 void (*ffun)(), void (*gfun)(),
692 double *x, double typf, double *typx, char *work)
693 /* int (*ffun)(), (*gfun)();*/
695 Iteration *iter = &myiter;
696 size_t i, j;
697 double nx0, ntypx;
699 n = (n > 0) ? n : 0;
700 k = (k > 0) ? k : 0;
702 iter->n = n;
703 iter->k = k;
704 iter->ffun = ffun; /* FIXME:AJR */
705 iter->gfun = gfun;
707 iter->x = (double *) work; work += sizeof(double) * (n + k);
708 iter->new_x = (double *) work; work += sizeof(double) * (n + k);
709 iter->sx = (double *) work; work += sizeof(double) * (n + k);
710 iter->delf = (double *) work; work += sizeof(double) * (n + k);
711 iter->new_delf = (double *) work; work += sizeof(double) * (n + k);
712 iter->qnstep = (double *) work; work += sizeof(double) * (n + k);
713 iter->F = (double *) work; work += sizeof(double) * (n + k);
714 for (i = 0; i < n; i++) {
715 iter->x[i] = x[i];
716 iter->sx[i] = (typx != nil && typx[i] > 0.0) ? 1.0 / typx[i] : 1.0;
718 for (i = 0; i < k; i++) {
719 iter->x[n + i] = x[n + i];
720 iter->sx[n + i] = 1.0;
723 iter->hessf = (double **) work; work += sizeof(double *) * (n + k);
724 for (i = 0; i < n + k; i++) {
725 iter->hessf[i] = (double *) work;
726 work += sizeof(double) * (n + k);
728 for (i = 0; i < n + k; i++)
729 for (j = 0; j < n + k; j++) iter->hessf[i][j] = 0.0;
730 iter->L = (double **) work; work += sizeof(double *) * (n + k);
731 for (i = 0; i < n + k; i++) {
732 iter->L[i] = (double *) work;
733 work += sizeof(double) * (n + k);
735 iter->H = (double **) work; work += sizeof(double *) * (n + k);
736 for (i = 0; i < n + k; i++) {
737 iter->H[i] = (double *) work;
738 work += sizeof(double) * (n + k);
741 iter->new_delg = (k > 0) ? (double **) work : nil;
742 work += sizeof(double *) * k;
743 for (i = 0; i < k; i++) {
744 iter->new_delg[i] = (double *) work;
745 work += sizeof(double) * n;
748 iter->typf = (typf > 0.0) ? typf : 1.0;
750 iter->verbose = VERBOSE;
751 iter->use_line_search = USE_SEARCH;
752 iter->itncount = 0;
753 iter->itnlimit = ITNLIMIT;
754 iter->gradtol = pow(macheps(), GRADTOL_POWER);
755 iter->steptol = pow(macheps(), STEPTOL_POWER);
756 for (i = 0, nx0 = 0.0, ntypx = 0.0; i < iter->n; i++) {
757 nx0 += fabs(iter->new_x[i] / iter->sx[i]);
758 ntypx += fabs(1.0 / iter->sx[i]);
760 iter->maxstep = MAX_STEP_FACTOR * Max(nx0, ntypx);
761 iter->values_supplied = FALSE;
764 void
765 minsetoptions(double gradtol, double steptol, double maxstep,
766 int itnlimit, int verbose, int use_search, int change_sign)
768 Iteration *iter = &myiter;
770 if (gradtol > 0.0) iter->gradtol = gradtol;
771 if (steptol > 0.0) iter->steptol = steptol;
772 if (maxstep > 0.0) iter->maxstep = maxstep;
773 if (itnlimit >= 0) iter->itnlimit = itnlimit;
774 if (verbose >= 0) iter->verbose = verbose;
775 iter->use_line_search = use_search;
776 iter->change_sign = change_sign;
779 void
780 minsupplyvalues(double f, double *delf, double **hessf,
781 double *g, double **delg)
783 Iteration *iter = &myiter;
784 size_t i, j, n, k;
786 n = iter->n;
787 k = iter->k;
789 iter->f = f;
790 for (i = 0; i < n; i++) {
791 iter->delf[i] = delf[i];
792 for (j = 0; j < k; j++)
793 iter->delf[i] += iter->x[n + j] * delg[j][i];
794 for (j = 0; j < n; j++) iter->hessf[i][j] = hessf[i][j];
796 for (i = 0; i < k; i++) {
797 iter->delf[n + i] = g[i];
798 for (j = 0; j < n; j++) {
799 iter->hessf[n + i][j] = delg[i][j];
800 iter->hessf[j][n + i] = delg[i][j];
803 if (k == 0) iter->crit = f;
804 else {
805 for (i = 0, iter->crit = 0.0; i < n + k; i++)
806 iter->crit += 0.5 * pow(fabs(iter->delf[i]), 2.0);
808 iter->values_supplied = TRUE;
811 void
812 minimize()
814 mindriver(&myiter);
817 void
818 minresults(double *x, double *pf, double *pcrit, double *delf, double **hessf,
819 double *g, double **delg, int *pcount, int *ptermcode, double *diagadd)
821 Iteration *iter = &myiter;
822 size_t i, j, n, k;
824 n = iter->n;
825 k = iter->k;
827 if (pf != nil) *pf = iter->f;
828 if (pcrit != nil) *pcrit = iter->crit;
829 for (i = 0; i < n; i++) {
830 if (x != nil) x[i] = iter->x[i];
831 if (delf != nil) delf[i] = iter->delf[i];
832 for (j = 0; j < n; j++) if (hessf != nil) hessf[i][j] = iter->hessf[i][j];
834 for (i = 0; i < k; i++) {
835 if (x != nil) x[n + i] = iter->x[n + i];
836 if (g != nil) g[i] = iter->delf[n + i];
837 for (j = 0; j < n; j++)
838 if (delg != nil) delg[i][j] = iter->hessf[n + i][j];
840 if (pcount != nil) *pcount = iter->itncount;
841 if (ptermcode != nil) *ptermcode = iter->termcode;
842 if (diagadd != nil) *diagadd = iter->diagadd;
846 double pdlogdet(int n, double **H)
848 int i;
849 double logdet, maxadd;
851 choldecomp(H, n, 0.0, &maxadd);
852 for (i = 0, logdet = 0.0; i < n; i++) logdet += 2.0 * log(H[i][i]);
853 return(logdet);
858 TODO
860 return amount added to make pos definite
862 scaling for constraints
864 alternate global strategies
866 callback function for verbose mode