diff --git a/CG_DESCENT-C-3.0/README b/CG_DESCENT-C-3.0/README new file mode 100644 --- /dev/null +++ b/CG_DESCENT-C-3.0/README @@ -0,0 +1,80 @@ +cg_descent is a conjugate gradient algorithm for solving +an unconstrained minimization problem of the form: + + min f (x) + +The algorithm is developed in the following papers +(see www.math.ufl.edu/~hager/papers/CG): + +[1] W. W. Hager and H. Zhang, A new conjugate gradient method + with guaranteed descent and an efficient line search, + SIAM Journal on Optimization, 16 (2005), 170-192. + +[2] W. W. Hager and H. Zhang, Algorithm 851: CG_DESCENT, + A conjugate gradient method with guaranteed descent, + ACM Transactions on Mathematical Software, 32 (2006), 113-137. + +[3] W. W. Hager and H. Zhang, A survey of nonlinear conjugate + gradient methods, Pacific Journal of Optimization, + 2 (2006), pp. 35-58. + +This directory provides a C implementation of cg_descent and +the codes needed to run cg_descent in the CUTEr testing +environment (http://hsl.rl.ac.uk/cuter-www/). The program +which calls cg_descent should include the header file cg_user.h. +Examples showing how to call cg_descent are given in driver1.c +through driver5.c. The user must provide routines to evaluate +the objective function and its gradient. Performance is often +improved if the user also provides a routine to simultaneously +evaluate the objective function and its gradient (see drive1.c). +In the simplest case, cg_descent is invoked with a statement +of the form: + +cg_descent (x, n, NULL, NULL, tol, myvalue, mygrad, NULL, NULL) ; + +where x is a pointer to an array which contains the starting +guess on input and the solution on output, n is the problem +dimension, tol is the computing tolerance (max norm of the +gradient), myvalue is a routine to evaluate the user's +objective function, and mygrad is a routine to evaluate +the gradient of the user's objective function. The 4 NULL +arguments could be replaced by the following (in order): +a structure to store execution statistics, a structure containing +algorithm parameters, a pointer to a routine which evaluates the +objective function and its gradient, and a pointer to a work +array. It the work array is not provided, then the code +allocates and frees memory. If the routine to simultaneously evaluate +objective function and its gradient is not provided, then the +code will use myvalue and mygrad to compute the value and +gradient independently. When the algorithm parameters are not +provided, then the default parameter values will be used +(see cg_default for their values). + +We also provide codes to interface cg_descent to the CUTEr +testing environment. The procedure for incorporating cg_descent +in CUTEr is as follows: + +1. Put cg_user.h into $CUTER/common/include +2. Put cg_descentma.c into $CUTER/common/src/tools +3. In $CUTER/common/src/tools, "gcc -lm -O3 -c cg_descentma.c" +4. "cp cg_descentma.o $MYCUTER/double/bin" +5. In the directory where you put cg_descent, type "make" and then + "cp cg_descent.o $MYCUTER/double/bin" +6. "cp cg_descent.pro $CUTER/build/prototypes" + "cp sdcg_descent.pro $CUTER/build/prototypes" +7. "cd $MYCUTER/bin" +8. type the following command twice: + +sed -f $MYCUTER/double/config/script.sed $CUTER/build/prototypes/pack.pro > pack + +where "pack" is first "cg_descent" and then "sdcg_descent" + +9. "chmod a+x cg_descent" and "chmod a+x sdcg_descent" + +You can run a problem by cd'ing to the directory where the sif files +are stored and typing, for example, "sdcg_descent BRYBND" + +NOTE: to run valgrind with the code, edit the program "runpackage" +found in "$MYCUTER/bin" as follows: +near the end of the program, change "$EXEC/${PAC}min" to +"valgrind $EXEC/${PAC}min" diff --git a/CG_DESCENT-C-3.0/cg_descent.c b/CG_DESCENT-C-3.0/cg_descent.c new file mode 100644 --- /dev/null +++ b/CG_DESCENT-C-3.0/cg_descent.c @@ -0,0 +1,1752 @@ +/* ========================================================================= + ============================ CG_DESCENT ================================= + ========================================================================= + ________________________________________________________________ + | A conjugate gradient method with guaranteed descent | + | C-code Version 1.1 (October 6, 2005) | + | Version 1.2 (November 14, 2005) | + | Version 2.0 (September 23, 2007) | + | Version 3.0 (May 18, 2008) | + | William W. Hager and Hongchao Zhang | + | hager@math.ufl.edu hzhang@math.ufl.edu | + | Department of Mathematics | + | University of Florida | + | Gainesville, Florida 32611 USA | + | 352-392-0281 x 244 | + | | + | Copyright by William W. Hager | + | | + | http://www.math.ufl.edu/~hager/papers/CG | + |________________________________________________________________| + ________________________________________________________________ + |This program is free software; you can redistribute it and/or | + |modify it under the terms of the GNU General Public License as | + |published by the Free Software Foundation; either version 2 of | + |the License, or (at your option) any later version. | + |This program is distributed in the hope that it will be useful, | + |but WITHOUT ANY WARRANTY; without even the implied warranty of | + |MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | + |GNU General Public License for more details. | + | | + |You should have received a copy of the GNU General Public | + |License along with this program; if not, write to the Free | + |Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, | + |MA 02110-1301 USA | + |________________________________________________________________|*/ + +#include "cg_user.h" +#include "cg_descent.h" +int cg_descent /* return: + -2 (function value became nan) + -1 (starting function value is nan) + 0 (convergence tolerance satisfied) + 1 (change in func <= feps*|f|) + 2 (total iterations exceeded maxit) + 3 (slope always negative in line search) + 4 (number secant iterations exceed nsecant) + 5 (search direction not a descent direction) + 6 (line search fails in initial interval) + 7 (line search fails during bisection) + 8 (line search fails during interval update) + 9 (debugger is on and the function value increases) + 10 (out of memory) */ +( + double *x, /* input: starting guess, output: the solution */ + INT n, /* problem dimension */ + cg_stats *Stat, /* structure with statistics (can be NULL) */ + cg_parameter *UParm, /* user parameters, NULL = use default parameters */ + double grad_tol, /* StopRule = 1: |g|_infty <= max (grad_tol, + StopFac*initial |g|_infty) [default] + StopRule = 0: |g|_infty <= grad_tol(1+|f|) */ + double (*value) (double *, INT), /* f = value (x, n) */ + void (*grad) (double *, double *, INT), /* grad (g, x, n) */ + double (*valgrad) (double *, double *, INT), /* f = valgrad (g, x, n), + NULL = compute value & gradient using value & grad */ + double *Work /* either size 4n work array or NULL */ +) +{ + INT n5, iter, maxit, nrestart, i ; + int status, StopRule ; + double delta2, eta_sq, Qk, Ck, + f, ftemp, gnorm, xnorm, gnorm2, dnorm2, denom, + t, t1, t2, t3, t4, t5, dphi, dphi0, alpha, talpha, + yk, ykyk, ykgk, dkyk, yk1, yk2, yk3, yk4, yk5, beta, tol, + *d, *g, *xtemp, *gtemp, *work ; + cg_parameter *Parm, ParmStruc ; + cg_com Com ; + +/* initialize the parameters */ + + if ( UParm == NULL ) + { + Parm = &ParmStruc ; + cg_default (Parm) ; + } + else Parm = UParm ; + Com.Parm = Parm ; + + if ( Parm->PrintParms ) cg_printParms (Parm) ; + + /* allocate work arrays */ + if ( Work == NULL ) work = malloc (4*n*sizeof (double)) ; + else work = Work ; + if ( work == NULL ) + { + printf ("Insufficient memory for specified problem dimension %e\n", + (double) n) ; + status = 10 ; + return (status) ; + } + Com.x = x ; + Com.d = d = work ; + g = d+n ; + Com.xtemp = xtemp = g+n ; + Com.gtemp = gtemp = xtemp+n ; + Com.n = n ; /* problem dimension */ + Com.nf = (INT) 0 ; /* number of function evaluations */ + Com.ng = (INT) 0 ; /* number of gradient evaluations */ + Com.AWolfe = Parm->AWolfe ; /* do not touch user's AWolfe */ + Com.cg_value = value ; + Com.cg_grad = grad ; + Com.cg_valgrad = valgrad ; + StopRule = Parm->StopRule ; + + /* the conjugate gradient algorithm is restarted every nrestart iteration */ + nrestart = (INT) (((double) n)*Parm->restart_fac) ; + + /* abort when number of iterations reaches maxit */ + if ( Parm->maxit_fac == INF ) maxit = INT_INF ; + else maxit = (INT) (((double) n)*Parm->maxit_fac) ; + + f = ZERO ; + n5 = n % 5 ; + + Ck = ZERO ; + Qk = ZERO ; + +/* initial function and gradient evaluations, initial direction */ + + f = cg_fg (g, x, &Com) ; + Com.f0 = f + f ; + xnorm = ZERO ; + for (i = 0; i < n5; i++) if ( xnorm < fabs (x [i]) ) xnorm = fabs (x [i]) ; + for (; i < n; i += 5) + { + if ( xnorm < fabs (x [i] ) ) xnorm = fabs (x [i] ) ; + if ( xnorm < fabs (x [i+1]) ) xnorm = fabs (x [i+1]) ; + if ( xnorm < fabs (x [i+2]) ) xnorm = fabs (x [i+2]) ; + if ( xnorm < fabs (x [i+3]) ) xnorm = fabs (x [i+3]) ; + if ( xnorm < fabs (x [i+4]) ) xnorm = fabs (x [i+4]) ; + } + gnorm = ZERO ; + gnorm2 = ZERO ; + for (i = 0; i < n5; i++) + { + t = g [i] ; + d [i] = -t ; + gnorm2 += t*t ; + if ( gnorm < fabs (t) ) gnorm = fabs (t) ; + } + for (; i < n;) + { + t1 = g [i] ; + d [i] = -t1 ; + if ( gnorm < fabs (t1) ) gnorm = fabs (t1) ; + i++ ; + + t2 = g [i] ; + d [i] = -t2 ; + if ( gnorm < fabs (t2) ) gnorm = fabs (t2) ; + i++ ; + + t3 = g [i] ; + d [i] = -t3 ; + if ( gnorm < fabs (t3) ) gnorm = fabs (t3) ; + i++ ; + + t4 = g [i] ; + d [i] = -t4 ; + if ( gnorm < fabs (t4) ) gnorm = fabs (t4) ; + i++ ; + + t5 = g [i] ; + d [i] = -t5 ; + if ( gnorm < fabs (t5) ) gnorm = fabs (t5) ; + i++ ; + + gnorm2 += t1*t1 + t2*t2 + t3*t3 + t4*t4 + t5*t5 ; + } + /* check that starting function value is nan */ + if ( f != f ) + { + status = -1 ; + goto Exit ; + } + + if ( Parm->StopRule ) tol = MAX (gnorm*Parm->StopFac, grad_tol) ; + else tol = grad_tol ; + + if ( Parm->PrintLevel >= 1 ) + { + printf ("iter: %5i f = %14.6e gnorm = %14.6e AWolfe = %2i\n", + (int) 0, f, gnorm, Com.AWolfe) ; + } + + if ( cg_tol (f, gnorm, StopRule, tol) ) + { + status = 0 ; + goto Exit ; + } + + dphi0 = -gnorm2 ; + delta2 = 2*Parm->delta - ONE ; + eta_sq = Parm->eta*Parm->eta ; + alpha = Parm->step ; + if ( alpha == 0. ) + { + alpha = Parm->psi0*xnorm/gnorm ; + if ( xnorm == ZERO ) + { + if ( f != ZERO ) alpha = Parm->psi0*fabs (f)/gnorm2 ; + else alpha = ONE ; + } + } + +/* start the conjugate gradient iteration + alpha starts as old step, ends as final step for current iteration + f is function value for alpha = 0 + QuadOK = TRUE means that a quadratic step was taken */ + + for (iter = 1; iter <= maxit; iter++) + { + Com.QuadOK = FALSE ; + alpha = Parm->psi2*alpha ; + if ( Parm->QuadStep ) + { + if ( f != ZERO ) t = fabs ((f-Com.f0)/f) ; + else t = ONE ; + if ( t > Parm->QuadCutOff ) /* take provisional step talpha */ + { + talpha = Parm->psi1*alpha ; + cg_step (xtemp, x, d, talpha, n) ; + ftemp = cg_f (xtemp, &Com) ; /* provisional function value */ + + /* check if function value is nan */ + if ( ftemp != ftemp ) /* reduce stepsize */ + { + for (i = 0; i < Parm->nexpand; i++) + { + talpha /= Parm->rho ; + cg_step (xtemp, x, d, talpha, n) ; + ftemp = cg_f (xtemp, &Com) ; + if ( ftemp == ftemp ) break ; + } + if ( i == Parm->nexpand ) + { + status = -2 ; + goto Exit ; + } + } + + if ( ftemp < f ) /* check if quadstep > 0 */ + { + denom = 2.*(((ftemp-f)/talpha)-dphi0) ; + if ( denom > ZERO ) /* try a quadratic fit step */ + { + Com.QuadOK = TRUE ; + alpha = -dphi0*talpha/denom ; + } + } + } + } + Com.f0 = f ; /* f0 saved as prior value */ + + if ( Parm->PrintLevel >= 1 ) + { + printf ("QuadOK: %2i initial a: %14.6e f0: %14.6e dphi: %14.6e\n", + Com.QuadOK, alpha, Com.f0, dphi0) ; + } + +/* parameters in Wolfe and approximate Wolfe conditions, and in update */ + + Qk = Parm->Qdecay*Qk + ONE ; + Ck = Ck + (fabs (f) - Ck)/Qk ; /* average cost magnitude */ + + if ( Parm->PertRule ) Com.fpert = f + Parm->eps*Ck ; + else Com.fpert = f + Parm->eps ; + + Com.wolfe_hi = Parm->delta*dphi0 ; + Com.wolfe_lo = Parm->sigma*dphi0 ; + Com.awolfe_hi = delta2*dphi0 ; + Com.alpha = alpha ; /* either prior step or quadratic fit step */ + Com.f = f ; + + if ( Com.AWolfe ) status = cg_line (dphi0, &Com) ; /* approx. Wolfe */ + else status = cg_lineW (dphi0, &Com) ;/* ordinary Wolfe */ + if ( (status > 0) && !Com.AWolfe )/*try approximate Wolfe line search*/ + { + if ( Parm->PrintLevel >= 1 ) + { + printf ("\nWOLFE LINE SEARCH FAILS\n") ; + } + Com.AWolfe = TRUE ; + status = cg_line (dphi0, &Com) ; + } + + alpha = Com.alpha ; + f = Com.f ; + dphi = Com.df ; + + if ( status ) goto Exit ; + +/*Test for convergence to within machine epsilon + [set feps to zero to remove this test] */ + + if ( -alpha*dphi0 <= Parm->feps*fabs (f) ) + { + status = 1 ; + goto Exit ; + } + +/* compute beta, yk2, gnorm, gnorm2, dnorm2, update x and g */ + + if ( iter % nrestart != 0 ) + { + cg_copy (x, xtemp, n) ; + dnorm2 = ZERO ; + for (i = 0; i < n5; i++) dnorm2 = dnorm2 + d [i]*d [i] ; + for (; i < n; i += 5) + { + dnorm2 = dnorm2 + d [i]*d [i] + d [i+1]*d [i+1] + + d [i+2]*d [i+2] + + d [i+3]*d [i+3] + + d [i+4]*d [i+4] ; + } + gnorm = ZERO ; + ykyk = ZERO ; + ykgk = ZERO ; + for (i = 0; i < n5; i++) + { + t = gtemp [i] ; + if ( gnorm < fabs (t) ) gnorm = fabs (t) ; + yk = t - g [i] ; + g [i] = t ; + ykgk += yk*t ; + ykyk += yk*yk ; + } + for (; i < n; ) + { + t1 = gtemp [i] ; + yk1 = t1 - g [i] ; + g [i] = t1 ; + if ( gnorm < fabs (t1) ) gnorm = fabs (t1) ; + i++ ; + + t2 = gtemp [i] ; + yk2 = t2 - g [i] ; + g [i] = t2 ; + if ( gnorm < fabs (t2) ) gnorm = fabs (t2) ; + i++ ; + + t3 = gtemp [i] ; + yk3 = t3 - g [i] ; + g [i] = t3 ; + if ( gnorm < fabs (t3) ) gnorm = fabs (t3) ; + i++ ; + + t4 = gtemp [i] ; + yk4 = t4 - g [i] ; + g [i] = t4 ; + if ( gnorm < fabs (t4) ) gnorm = fabs (t4) ; + i++ ; + + t5 = gtemp [i] ; + yk5 = t5 - g [i] ; + g [i] = t5 ; + if ( gnorm < fabs (t5) ) gnorm = fabs (t5) ; + + i++ ; + ykyk += yk1*yk1 + yk2*yk2 + yk3*yk3 + yk4*yk4 + yk5*yk5 ; + ykgk += yk1*t1 + yk2*t2 + yk3*t3 + yk4*t4 + yk5*t5 ; + } + + if ( cg_tol (f, gnorm, StopRule, tol) ) + { + status = 0 ; + goto Exit ; + } + dkyk = dphi - dphi0 ; + beta = (ykgk - 2.*dphi*ykyk/dkyk)/dkyk ; +/* + faster: initialize dnorm2 = gnorm2 at start, then + dnorm2 = gnorm2 + beta**2*dnorm2 - 2.*beta*dphi + gnorm2 = ||g_{k+1}||^2 + dnorm2 = ||d_{k+1}||^2 + dpi = g_{k+1}' d_k */ + + t = -ONE/sqrt (dnorm2*MIN (eta_sq, gnorm2)) ; + beta = MAX (beta, t) ; + +/* update search direction d = -g + beta*dold */ + + gnorm2 = ZERO ; + for (i = 0; i < n5; i++) + { + t = g [i] ; + d [i] = -t + beta*d [i] ; + gnorm2 += t*t ; + } + for (; i < n; ) + { + t1 = g [i] ; + d [i] = -t1 + beta*d [i] ; + i++ ; + + t2 = g [i] ; + d [i] = -t2 + beta*d [i] ; + i++ ; + + t3 = g [i] ; + d [i] = -t3 + beta*d [i] ; + i++ ; + + t4 = g [i] ; + d [i] = -t4 + beta*d [i] ; + i++ ; + + t5 = g [i] ; + d [i] = -t5 + beta*d [i] ; + i++ ; + + gnorm2 += t1*t1 + t2*t2 + t3*t3 + t4*t4 + t5*t5 ; + } + dphi0 = -gnorm2 + beta*dphi ; + if ( Parm->debug ) /* Check the dphi0 = d'g */ + { + t = ZERO ; + for (i = 0; i < n; i++) t = t + d [i]*g [i] ; + if ( fabs(t-dphi0) > Parm->debugtol*fabs(dphi0) ) + { + printf("Warning, dphi0 != d'g!\n"); + printf("dphi0:%14.6e, d'g:%14.6e\n",dphi0, t) ; + } + } + } + else + { + /* search direction d = -g */ + if ( Parm->PrintLevel >= 1 ) printf ("RESTART CG\n") ; + gnorm = ZERO ; + gnorm2 = ZERO ; + cg_copy (x, xtemp, n) ; + for (i = 0; i < n5; i++) + { + t = gtemp [i] ; + g [i] = t ; + d [i] = -t ; + if ( gnorm < fabs (t) ) gnorm = fabs (t) ; + gnorm2 += t*t ; + } + for (; i < n; ) + { + t1 = gtemp [i] ; + g [i] = t1 ; + d [i] = -t1 ; + if ( gnorm < fabs (t1) ) gnorm = fabs (t1) ; + i++ ; + + t2 = gtemp [i] ; + g [i] = t2 ; + d [i] = -t2 ; + if ( gnorm < fabs (t2) ) gnorm = fabs (t2) ; + i++ ; + + t3 = gtemp [i] ; + g [i] = t3 ; + d [i] = -t3 ; + if ( gnorm < fabs (t3) ) gnorm = fabs (t3) ; + i++ ; + + t4 = gtemp [i] ; + g [i] = t4 ; + d [i] = -t4 ; + if ( gnorm < fabs (t4) ) gnorm = fabs (t4) ; + i++ ; + + t5 = gtemp [i] ; + g [i] = t5 ; + d [i] = -t5 ; + if ( gnorm < fabs (t5) ) gnorm = fabs (t5) ; + i++ ; + gnorm2 += t1*t1 + t2*t2 + t3*t3 + t4*t4 + t5*t5 ; + } + if ( cg_tol (f, gnorm, StopRule, tol) ) + { + status = 0 ; + goto Exit ; + } + dphi0 = -gnorm2 ; + } + if ( !Com.AWolfe ) + { + if ( fabs (f-Com.f0) < Parm->AWolfeFac*Ck ) Com.AWolfe = TRUE ; + } + + if ( Parm->PrintLevel >= 1 ) + { + printf ("\niter: %5i f = %14.6e gnorm = %14.6e AWolfe = %2i\n", + (int) iter, f, gnorm, Com.AWolfe) ; + } + + if ( Parm->debug ) + { + if ( f > Com.f0 + Parm->debugtol*Ck ) + { + status = 9 ; + goto Exit ; + } + } + + if ( dphi0 > ZERO ) + { + status = 5 ; + goto Exit ; + } + } + status = 2 ; + +Exit: + if ( Stat != NULL ) + { + Stat->f = f ; + Stat->gnorm = gnorm ; + Stat->nfunc = Com.nf ; + Stat->ngrad = Com.ng ; + Stat->iter = iter ; + } + if ( status > 2 ) + { + gnorm = ZERO ; + for (i = 0; i < n; i++) + { + x [i] = xtemp [i] ; + g [i] = gtemp [i] ; + t = fabs (g [i]) ; + gnorm = MAX (gnorm, t) ; + } + if ( Stat != NULL ) Stat->gnorm = gnorm ; + } + if ( Parm->PrintFinal || Parm->PrintLevel >= 1 ) + { + const char mess1 [] = "Possible causes of this error message:" ; + const char mess2 [] = " - your tolerance may be too strict: " + "grad_tol = " ; + const char mess3 [] = "Line search fails" ; + const char mess4 [] = " - your gradient routine has an error" ; + const char mess5 [] = " - the parameter epsilon in cg_descent_c.parm " + "is too small" ; + printf ("\nTermination status: %i\n", status) ; + if ( status == -2 ) + { + printf ("At iteration %10.0f function value became nan\n", + (double) iter) ; + } + else if ( status == -1 ) + { + printf ("Objective function value is nan at starting point\n") ; + } + else if ( status == 0 ) + { + printf ("Convergence tolerance for gradient satisfied\n") ; + } + else if ( status == 1 ) + { + printf ("Terminating since change in function value " + "<= feps*|f|\n") ; + } + else if ( status == 2 ) + { + printf ("Number of iterations exceed specified limit\n") ; + printf ("Iterations: %10.0f maxit: %10.0f\n", + (double) iter, (double) maxit) ; + printf ("%s\n", mess1) ; + printf ("%s %e\n", mess2, grad_tol) ; + } + else if ( status == 3 ) + { + printf ("Slope always negative in line search\n") ; + printf ("%s\n", mess1) ; + printf (" - your cost function has an error\n") ; + printf ("%s\n", mess4) ; + } + else if ( status == 4 ) + { + printf ("Line search fails, too many secant steps\n") ; + printf ("%s\n", mess1) ; + printf ("%s %e\n", mess2, grad_tol) ; + } + else if ( status == 5 ) + { + printf ("Search direction not a descent direction\n") ; + } + else if ( status == 6 ) /* line search fails */ + { + printf ("%s\n", mess3) ; + printf ("%s\n", mess1) ; + printf ("%s %e\n", mess2, grad_tol) ; + printf ("%s\n", mess4) ; + printf ("%s\n", mess5) ; + } + else if ( status == 7 ) /* line search fails */ + { + printf ("%s\n", mess3) ; + printf ("%s\n", mess1) ; + printf ("%s %e\n", mess2, grad_tol) ; + } + else if ( status == 8 ) /* line search fails */ + { + printf ("%s\n", mess3) ; + printf ("%s\n", mess1) ; + printf ("%s %e\n", mess2, grad_tol) ; + printf ("%s\n", mess4) ; + printf ("%s\n", mess5) ; + } + else if ( status == 9 ) + { + printf ("Debugger is on, function value does not improve\n") ; + printf ("new value: %25.16e old value: %25.16e\n", f, Com.f0) ; + } + else if ( status == 10 ) + { + printf ("Insufficient memory\n") ; + } + + printf ("maximum norm for gradient: %13.6e\n", gnorm) ; + printf ("function value: %13.6e\n\n", f) ; + printf ("cg iterations: %10.0f\n", (double) iter) ; + printf ("function evaluations: %10.0f\n", (double) Com.nf) ; + printf ("gradient evaluations: %10.0f\n", (double) Com.ng) ; + printf ("===================================\n\n") ; + } + if ( Work == NULL ) free (work) ; + return (status) ; +} + +/* ========================================================================= + === cg_default ========================================================== + ========================================================================= + Set default conjugate gradient parameter values. If the parameter argument + of cg_descent is NULL, this routine is called by cg_descent automatically. + If the user wishes to set parameter values, then the cg_parameter structure + should be allocated in the main program. The user could call cg_default + to initialize the structure, and then individual elements in the structure + could be changed, before passing the structure to cg_descent. + =========================================================================*/ +void cg_default +( + cg_parameter *Parm +) +{ + /* T => print final function value + F => no printout of final function value */ + Parm->PrintFinal = TRUE ; + + /* Level 0 = no printing, ... , Level 3 = maximum printing */ + Parm->PrintLevel = 0 ; + + /* T => print parameters values + F => do not display parmeter values */ + Parm->PrintParms = FALSE ; + + /* T => use approximate Wolfe line search + F => use ordinary Wolfe line search, switch to approximate Wolfe when + |f_k+1-f_k| < AWolfeFac*C_k, C_k = average size of cost */ + Parm->AWolfe = FALSE ; + Parm->AWolfeFac = 1.e-3 ; + + /* factor in [0, 1] used to compute average cost magnitude C_k as follows: + Q_k = 1 + (Qdecay)Q_k-1, Q_0 = 0, C_k = C_k-1 + (|f_k| - C_k-1)/Q_k */ + Parm->Qdecay = .7 ; + + /* Stop Rules: + T => ||grad||_infty <= max(grad_tol, initial |grad|_infty*StopFact) + F => ||grad||_infty <= grad_tol*(1 + |f_k|) */ + Parm->StopRule = TRUE ; + Parm->StopFac = 0.e-12 ; + + /* T => estimated error in function value is eps*Ck, + F => estimated error in function value is eps */ + Parm->PertRule = TRUE ; + Parm->eps = 1.e-6 ; + + /* T => attempt quadratic interpolation in line search when + |f_k+1 - f_k|/f_k <= QuadCutoff + F => no quadratic interpolation step */ + Parm->QuadStep = TRUE ; + Parm->QuadCutOff = 1.e-12 ; + + /* T => check that f_k+1 - f_k <= debugtol*C_k + F => no checking of function values */ + Parm->debug = FALSE ; + Parm->debugtol = 1.e-10 ; + + /* if step is nonzero, it is the initial step of the initial line search */ + Parm->step = ZERO ; + + /* abort cg after maxit_fac*n iterations */ + Parm->maxit_fac = INF ; + + /* maximum number of times the bracketing interval grows or shrinks + in the line search is nexpand */ + Parm->nexpand = (int) 50 ; + + /* maximum number of secant iterations in line search is nsecant */ + Parm->nsecant = (int) 50 ; + + /* conjugate gradient method restarts after (n*restart_fac) iterations */ + Parm->restart_fac = ONE ; + + /* stop when -alpha*dphi0 (estimated change in function value) <= feps*|f|*/ + Parm->feps = ZERO ; + + /* after encountering nan, growth factor when searching for + a bracketing interval */ + Parm->nan_rho = 1.3 ; + + /* Wolfe line search parameter, range [0, .5] + phi (a) - phi (0) <= delta phi'(0) */ + Parm->delta = .1 ; + + /* Wolfe line search parameter, range [delta, 1] + phi' (a) >= sigma phi' (0) */ + Parm->sigma = .9 ; + + /* decay factor for bracket interval width in line search, range (0, 1) */ + Parm->gamma = .66 ; + + /* growth factor in search for initial bracket interval */ + Parm->rho = 5. ; + + /* conjugate gradient parameter beta_k must be >= eta*||d_k||_2 */ + Parm->eta = .01 ; + + /* starting guess for line search = + psi0 ||x_0||_infty over ||g_0||_infty if x_0 != 0 + psi0 |f(x_0)|/||g_0||_2 otherwise */ + Parm->psi0 = .01 ; /* factor used in starting guess for iteration 1 */ + + /* for a QuadStep, function evalutated at psi1*previous step */ + Parm->psi1 = .1 ; + + /* when starting a new cg iteration, our initial guess for the line + search stepsize is psi2*previous step */ + Parm->psi2 = 2. ; +} + +/* ========================================================================= + ==== cg_Wolfe =========================================================== + ========================================================================= + Check whether the Wolfe or the approximate Wolfe conditions are satisfied + ========================================================================= */ +int cg_Wolfe +( + double alpha, /* stepsize */ + double f, /* function value associated with stepsize alpha */ + double dphi, /* derivative value associated with stepsize alpha */ + cg_com *Com /* cg com */ +) +{ + if ( dphi >= Com->wolfe_lo ) + { + +/* test original Wolfe conditions */ + + if ( f - Com->f0 <= alpha*Com->wolfe_hi ) + { + if ( Com->Parm->PrintLevel >= 2 ) + { + printf ("wolfe f: %14.6e f0: %14.6e dphi: %14.6e\n", + f, Com->f0, dphi) ; + } + return (1) ; + } +/* test approximate Wolfe conditions */ + else if ( Com->AWolfe ) + { + if ( (f <= Com->fpert) && (dphi <= Com->awolfe_hi) ) + { + if ( Com->Parm->PrintLevel >= 2 ) + { + printf ("f: %14.6e fpert: %14.6e dphi: %14.6e awolf_hi: " + "%14.6e\n", f, Com->fpert, dphi, Com->awolfe_hi) ; + } + return (1) ; + } + } + } + return (0) ; +} + +/* ========================================================================= + ==== cg_f =============================================================== + Evaluate the function + =========================================================================*/ +double cg_f +( + double *x, + cg_com *Com +) +{ + double f ; + f = Com->cg_value (x, Com->n) ; + Com->nf++ ; + return (f) ; +} + +/* ========================================================================= + ==== cg_g =============================================================== + Evaluate the gradient + =========================================================================*/ +void cg_g +( + double *g, + double *x, + cg_com *Com +) +{ + Com->cg_grad (g, x, Com->n) ; + Com->ng++ ; +} + +/* ========================================================================= + ==== cg_fg ============================================================== + Evaluate the function and gradient + =========================================================================*/ +double cg_fg +( + double *g, + double *x, + cg_com *Com +) +{ + double f ; + if ( Com->cg_valgrad != NULL ) f = Com->cg_valgrad (g, x, Com->n) ; + else + { + Com->cg_grad (g, x, Com->n) ; + f = Com->cg_value (x, Com->n) ; + } + Com->nf++ ; + Com->ng++ ; + return (f) ; +} + +/* ========================================================================= + ==== cg_tol ============================================================= + ========================================================================= + Check for convergence + ========================================================================= */ +int cg_tol +( + double f, /* function value associated with stepsize */ + double gnorm, /* gradient sup-norm */ + int StopRule, /* T => |grad|_infty <=max (tol, |grad|_infty*StopFact) + F => |grad|_infty <= tol*(1+|f|)) */ + double tol /* tolerance */ +) +{ + if ( StopRule ) + { + if ( gnorm <= tol ) return (1) ; + } + else if ( gnorm <= tol*(ONE + fabs (f)) ) return (1) ; + return (0) ; +} + +/* ========================================================================= + ==== cg_dot ============================================================= + ========================================================================= + Compute dot product of x and y, vectors of length n + ========================================================================= */ +double cg_dot +( + double *x, /* first vector */ + double *y, /* second vector */ + INT n /* length of vectors */ +) +{ + INT i, n5 ; + double t ; + t = 0. ; + n5 = n % 5 ; + for (i = 0; i < n5; i++) t += x [i]*y [i] ; + for (; i < n; i += 5) + { + t += x [i]*y[i] + x [i+1]*y [i+1] + x [i+2]*y [i+2] + + x [i+3]*y [i+3] + x [i+4]*y [i+4] ; + } + return (t) ; +} + +/* ========================================================================= + === cg_copy ============================================================= + ========================================================================= + Copy vector x into vector y + ========================================================================= */ +void cg_copy +( + double *y, /* output of copy */ + double *x, /* input of copy */ + int n /* length of vectors */ +) +{ + int j, n10 ; + n10 = n % 10 ; + for (j = 0; j < n10; j++) y [j] = x [j] ; + for (; j < n; j += 10) + { + y [j] = x [j] ; + y [j+1] = x [j+1] ; + y [j+2] = x [j+2] ; + y [j+3] = x [j+3] ; + y [j+4] = x [j+4] ; + y [j+5] = x [j+5] ; + y [j+6] = x [j+6] ; + y [j+7] = x [j+7] ; + y [j+8] = x [j+8] ; + y [j+9] = x [j+9] ; + } +} + +/* ========================================================================= + ==== cg_step ============================================================ + ========================================================================= + Compute xtemp = x + alpha d + ========================================================================= */ +void cg_step +( + double *xtemp, /*output vector */ + double *x, /* initial vector */ + double *d, /* search direction */ + double alpha, /* stepsize */ + INT n /* length of the vectors */ +) +{ + INT n5, i ; + n5 = n % 5 ; + for (i = 0; i < n5; i++) xtemp [i] = x[i] + alpha*d[i] ; + for (; i < n; i += 5) + { + xtemp [i] = x [i] + alpha*d [i] ; + xtemp [i+1] = x [i+1] + alpha*d [i+1] ; + xtemp [i+2] = x [i+2] + alpha*d [i+2] ; + xtemp [i+3] = x [i+3] + alpha*d [i+3] ; + xtemp [i+4] = x [i+4] + alpha*d [i+4] ; + } +} + +/* ========================================================================= + ==== cg_line ============================================================ + ========================================================================= + Approximate Wolfe line search routine + ========================================================================= */ +int cg_line +( + double dphi0, /* function derivative at starting point (alpha = 0) */ + cg_com *Com /* cg com structure */ +) +{ + INT n, iter ; + int i, nsecant, nshrink, ngrow, status ; + double a, dphia, b, dphib, c, alpha, phi, dphi, + a0, da0, b0, db0, width, fquad, rho, *x, *xtemp, *d, *gtemp ; + cg_parameter *Parm ; + + Parm = Com->Parm ; + if ( Parm->PrintLevel >= 1 ) printf ("Approximate Wolfe line search\n") ; + alpha = Com->alpha ; + phi = Com->f ; + n = Com->n ; + x = Com->x ; /* current iterate */ + xtemp = Com->xtemp ; /* x + alpha*d */ + d = Com->d ; /* current search direction */ + gtemp = Com->gtemp ; /* gradient at x + alpha*d */ + rho = Parm->rho ; + cg_step (xtemp, x, d, alpha, n) ; + cg_g (gtemp, xtemp, Com) ; + dphi = cg_dot (gtemp, d, n) ; + + /* check if gradient is nan; if so, reduce stepsize */ + if ( dphi != dphi ) + { + for (i = 0; i < Parm->nexpand; i++) + { + alpha /= rho ; + cg_step (xtemp, x, d, alpha, n) ; + cg_g (gtemp, xtemp, Com) ; + dphi = cg_dot (gtemp, d, n) ; + if ( dphi == dphi ) break ; + } + if ( i == Parm->nexpand ) + { + status = -2 ; + goto Exit ; + } + rho = Parm->nan_rho ; + } + +/*Find initial interval [a,b] such that dphia < 0, dphib >= 0, + and phia <= phi0 + feps*fabs (phi0) */ + + a = ZERO ; + dphia = dphi0 ; + ngrow = 0 ; + nshrink = 0 ; + while ( dphi < ZERO ) + { + phi = cg_f (xtemp, Com) ; + +/* if quadstep in effect and quadratic conditions hold, check wolfe condition*/ + + if ( Com->QuadOK ) + { + if ( ngrow == 0 ) fquad = MIN (phi, Com->f0) ; + if ( phi <= fquad ) + { + if ( Parm->PrintLevel >= 2 ) + { + printf ("alpha: %14.6e phi: %14.6e fquad: %14.6e\n", + alpha, phi, fquad) ; + } + if ( cg_Wolfe (alpha, phi, dphi, Com) ) + { + status = 0 ; + goto Exit ; + } + } + } + if ( phi > Com->fpert ) + { + /* contraction phase, only break at termination or Secant step */ + b = alpha ; + while ( TRUE ) + { + alpha = .5*(a+b) ; + nshrink++ ; + if ( nshrink > Parm->nexpand ) + { + status = 6 ; + goto Exit ; + } + cg_step (xtemp, x, d, alpha, n) ; + cg_g (gtemp, xtemp, Com) ; + dphi = cg_dot (gtemp, d, n) ; + if ( dphi >= ZERO ) goto Secant ; + phi = cg_f (xtemp, Com) ; + if ( Parm->PrintLevel >= 2 ) + { + printf ("contract, a: %14.6e b: %14.6e alpha: %14.6e phi: " + "%14.6e dphi: %14.6e\n", a, b, alpha, phi, dphi) ; + } + if ( Com->QuadOK && (phi <= fquad) ) + { + if ( cg_Wolfe (alpha, phi, dphi, Com) ) + { + status = 0 ; + goto Exit ; + } + } + if ( phi <= Com->fpert ) + { + a = alpha ; + dphia = dphi ; + } + else + { + b = alpha ; + } + } + } + +/* expansion phase */ + + a = alpha ; + dphia = dphi ; + ngrow++ ; + if ( ngrow > Parm->nexpand ) + { + status = 3 ; + goto Exit ; + } + alpha = rho*alpha ; + cg_step (xtemp, x, d, alpha, n) ; + cg_g (gtemp, xtemp, Com) ; + dphi = cg_dot (gtemp, d, n) ; + if ( Parm->PrintLevel >= 2 ) + { + printf ("expand, a: %14.6e alpha: %14.6e phi: " + "%14.6e dphi: %14.6e\n", a, alpha, phi, dphi) ; + } + } + +Secant: + b = alpha ; + dphib = dphi ; + if ( Com->QuadOK ) + { + phi = cg_f (xtemp, Com) ; + if ( ngrow + nshrink == 0 ) fquad = MIN (phi, Com->f0) ; + if ( phi <= fquad ) + { + if ( cg_Wolfe (alpha, phi, dphi, Com) ) + { + status = 0 ; + goto Exit ; + } + } + } + nsecant = Parm->nsecant ; + for (iter = 1; iter <= nsecant; iter++) + { + if ( Parm->PrintLevel >= 2 ) + { + printf ("secant, a: %14.6e b: %14.6e da: %14.6e db: %14.6e\n", + a, b, dphia, dphib) ; + } + width = Parm->gamma*(b - a) ; + if ( -dphia <= dphib ) alpha = a - (a-b)*(dphia/(dphia-dphib)) ; + else alpha = b - (a-b)*(dphib/(dphia-dphib)) ; + c = alpha ; + a0 = a ; + b0 = b ; + da0 = dphia ; + db0 = dphib ; + status = cg_update (&a, &dphia, &b, &dphib, &alpha, &phi, &dphi, Com) ; + if ( status >= 0 ) goto Exit ; + else if ( status == -2 ) + { + if ( c == a ) + { + if ( dphi > da0 ) alpha = c - (c-a0)*(dphi/(dphi-da0)) ; + else alpha = a ; + } + else + { + if ( dphi < db0 ) alpha = c - (c-b0)*(dphi/(dphi-db0)) ; + else alpha = b ; + } + if ( (alpha > a) && (alpha < b) ) + { + if ( Parm->PrintLevel >= 2 ) printf ("2nd secant\n") ; + status = cg_update (&a, &dphia, &b, &dphib, &alpha, &phi, + &dphi, Com) ; + if ( status >= 0 ) goto Exit ; + } + } + +/* bisection iteration */ + + if ( b-a >= width ) + { + alpha = .5*(b+a) ; + if ( Parm->PrintLevel >= 2 ) printf ("bisection\n") ; + status = cg_update (&a, &dphia, &b, &dphib, &alpha, &phi, + &dphi, Com) ; + if ( status >= 0 ) goto Exit ; + } + else if ( b <= a ) + { + status = 7 ; + goto Exit ; + } + } + status = 4 ; + +Exit: + Com->alpha = alpha ; + Com->f = phi ; + Com->df = dphi ; + return (status) ; +} + +/* ========================================================================= + ==== cg_lineW =========================================================== + ========================================================================= + Ordinary Wolfe line search routine. + This routine is identical to cg_line except that the function + psi [a] = phi [a] - phi [0] - a*delta*dphi [0] is minimized instead of + the function phi + ========================================================================= */ +int cg_lineW +( + double dphi0, /* function derivative at starting point (alpha = 0) */ + cg_com *Com /* cg com structure */ +) +{ + INT n, iter ; + int i, nsecant, nshrink, ngrow, status ; + double a, dpsia, b, dpsib, c, alpha, phi, dphi, + a0, da0, b0, db0, width, fquad, rho, psi, dpsi, + *x, *xtemp, *d, *gtemp ; + cg_parameter *Parm ; + + Parm = Com->Parm ; + if ( Parm->PrintLevel >= 1 ) printf ("Wolfe line search\n") ; + alpha = Com->alpha ; + phi = Com->f ; + dphi = Com->df ; + n = Com->n ; + x = Com->x ; /* current iterate */ + xtemp = Com->xtemp ; /* x + alpha*d */ + d = Com->d ; /* current search direction */ + gtemp = Com->gtemp ; /* gradient at x + alpha*d */ + rho = Parm->rho ; + cg_step (xtemp, x, d, alpha, n) ; + cg_g (gtemp, xtemp, Com) ; + dphi = cg_dot (gtemp, d, n) ; + + /* check if gradient is nan; if so, reduce stepsize */ + if ( dphi != dphi ) + { + for (i = 0; i < Parm->nexpand; i++) + { + alpha /= rho ; + cg_step (xtemp, x, d, alpha, n) ; + cg_g (gtemp, xtemp, Com) ; + dphi = cg_dot (gtemp, d, n) ; + if ( dphi == dphi ) break ; + } + if ( i == Parm->nexpand ) + { + status = -2 ; + goto Exit ; + } + rho = Parm->nan_rho ; + } + dpsi = dphi - Com->wolfe_hi ; + +/*Find initial interval [a,b] such that dphia < 0, dphib >= 0, + and phia <= phi0 + feps*fabs (phi0) */ + + a = ZERO ; + dpsia = dphi0 - Com->wolfe_hi ; + ngrow = 0 ; + nshrink = 0 ; + while ( dpsi < ZERO ) + { + phi = cg_f (xtemp, Com) ; + psi = phi - alpha*Com->wolfe_hi ; + +/* if quadstep in effect and quadratic conditions hold, check Wolfe condition*/ + + if ( Com->QuadOK ) + { + if ( ngrow == 0 ) fquad = MIN (phi, Com->f0) ; + if ( phi <= fquad ) + { + if ( Parm->PrintLevel >= 2 ) + { + printf ("alpha: %14.6e phi: %14.6e fquad: %14.6e\n", + alpha, phi, fquad) ; + } + if ( cg_Wolfe (alpha, phi, dphi, Com) ) + { + status = 0 ; + goto Exit ; + } + } + } + if ( psi <= Com->fpert ) + { + a = alpha ; + dpsia = dphi ; + } + else + { + /* contraction phase, only break at termination or Secant step */ + b = alpha ; + while ( TRUE ) + { + alpha = .5*(a+b) ; + nshrink++ ; + if ( nshrink > Parm->nexpand ) + { + status = 6 ; + goto Exit ; + } + cg_step (xtemp, x, d, alpha, n) ; + cg_g (gtemp, xtemp, Com) ; + dphi = cg_dot (gtemp, d, n) ; + dpsi = dphi - Com->wolfe_hi ; + if ( dpsi >= ZERO ) goto Secant ; + phi = cg_f (xtemp, Com) ; + psi = phi - alpha*Com->wolfe_hi ; + if ( Parm->PrintLevel >= 2 ) + { + printf ("contract, a: %14.6e b: %14.6e alpha: %14.6e phi: " + "%14.6e dphi: %14.6e\n", a, b, alpha, phi, dphi) ; + } + if ( Com->QuadOK && (phi <= fquad) ) + { + if ( cg_Wolfe (alpha, phi, dphi, Com) ) + { + status = 0 ; + goto Exit ; + } + } + if ( psi <= Com->fpert ) + { + a = alpha ; + dpsia = dpsi ; + } + else + { + b = alpha ; + } + } + } + +/* expansion phase */ + + ngrow++ ; + if ( ngrow > Parm->nexpand ) + { + status = 3 ; + goto Exit ; + } + alpha *= rho ; + cg_step (xtemp, x, d, alpha, n) ; + cg_g (gtemp, xtemp, Com) ; + dphi = cg_dot (gtemp, d, n) ; + dpsi = dphi - Com->wolfe_hi ; + if ( Parm->PrintLevel >= 2 ) + { + printf ("expand, a: %14.6e alpha: %14.6e phi: " + "%14.6e dphi: %14.6e\n", a, alpha, phi, dphi) ; + } + } + +Secant: + b = alpha ; + dpsib = dpsi ; + if ( Com->QuadOK ) + { + phi = cg_f (xtemp, Com) ; + if ( ngrow + nshrink == 0 ) fquad = MIN (phi, Com->f0) ; + if ( phi <= fquad ) + { + if ( cg_Wolfe (alpha, phi, dphi, Com) ) + { + status = 0 ; + goto Exit ; + } + } + } + nsecant = Parm->nsecant ; + for (iter = 1; iter <= nsecant; iter++) + { + if ( Parm->PrintLevel >= 2 ) + { + printf ("secant, a: %14.6e b: %14.6e da: %14.6e db: %14.6e\n", + a, b, dpsia, dpsib) ; + } + width = Parm->gamma*(b - a) ; + if ( -dpsia <= dpsib ) alpha = a - (a-b)*(dpsia/(dpsia-dpsib)) ; + else alpha = b - (a-b)*(dpsib/(dpsia-dpsib)) ; + c = alpha ; + a0 = a ; + b0 = b ; + da0 = dpsia ; + db0 = dpsib ; + status = cg_updateW (&a, &dpsia, &b, &dpsib, &alpha, &phi, &dphi, + &dpsi, Com) ; + if ( status >= 0 ) goto Exit ; + else if ( status == -2 ) + { + if ( c == a ) + { + if ( dpsi > da0 ) alpha = c - (c-a0)*(dpsi/(dpsi-da0)) ; + else alpha = a ; + } + else + { + if ( dpsi < db0 ) alpha = c - (c-b0)*(dpsi/(dpsi-db0)) ; + else alpha = b ; + } + if ( (alpha > a) && (alpha < b) ) + { + if ( Parm->PrintLevel >= 2 ) printf ("2nd secant\n") ; + status = cg_updateW (&a, &dpsia, &b, &dpsib, &alpha, &phi, + &dphi, &dpsi, Com) ; + if ( status >= 0 ) goto Exit ; + } + } + +/* bisection iteration */ + + if ( b-a >= width ) + { + alpha = .5*(b+a) ; + if ( Parm->PrintLevel >= 2 ) printf ("bisection\n") ; + status = cg_updateW (&a, &dpsia, &b, &dpsib, &alpha, &phi, &dphi, + &dpsi, Com) ; + if ( status >= 0 ) goto Exit ; + } + else if ( b <= a ) + { + status = 7 ; + goto Exit ; + } + } + status = 4 ; + +Exit: + Com->alpha = alpha ; + Com->f = phi ; + Com->df = dphi ; + return (status) ; +} + +/* ========================================================================= + ==== cg_update ========================================================== + ========================================================================= + update returns: 8 if too many iterations + 0 if Wolfe condition is satisfied + -1 if interval is updated and a search is done + -2 if the interval updated successfully + ========================================================================= */ +int cg_update +( + double *a , /* left side of bracketing interval */ + double *dphia , /* derivative at a */ + double *b , /* right side of bracketing interval */ + double *dphib , /* derivative at b */ + double *alpha , /* trial step (between a and b) */ + double *phi , /* function value at alpha (returned) */ + double *dphi , /* function derivative at alpha (returned) */ + cg_com *Com /* cg com structure */ +) +{ + INT n ; + int nshrink, status ; + double *x, *xtemp, *d, *gtemp ; + cg_parameter *Parm ; + + Parm = Com->Parm ; + n = Com->n ; + x = Com->x ; /* current iterate */ + xtemp = Com->xtemp ; /* x + alpha*d */ + d = Com->d ; /* current search direction */ + gtemp = Com->gtemp ; /* gradient at x + alpha*d */ + cg_step (xtemp, x, d, *alpha, n) ; + *phi = cg_fg (gtemp, xtemp, Com) ; + *dphi = cg_dot (gtemp, d, n) ; + if ( Parm->PrintLevel >= 2 ) + { + printf ("update alpha: %14.6e phi: %14.6e dphi: %14.6e\n", + *alpha, *phi, *dphi) ; + } + if ( cg_Wolfe (*alpha, *phi, *dphi, Com) ) + { + status = 0 ; + goto Exit2 ; + } + status = -2 ; + if ( *dphi >= ZERO ) + { + *b = *alpha ; + *dphib = *dphi ; + goto Exit2 ; + } + else + { + if ( *phi <= Com->fpert ) + { + *a = *alpha ; + *dphia = *dphi ; + goto Exit2 ; + } + } + nshrink = 0 ; + *b = *alpha ; + while ( TRUE ) + { + *alpha = .5*(*a + *b) ; + nshrink++ ; + if ( nshrink > Parm->nexpand ) + { + status = 8 ; + goto Exit2 ; + } + cg_step (xtemp, x, d, *alpha, n) ; + *phi = cg_fg (gtemp, xtemp, Com) ; + *dphi = cg_dot (gtemp, d, n) ; + if ( Parm->PrintLevel >= 2 ) + { + printf ("contract, a: %14.6e alpha: %14.6e " + "phi: %14.6e dphi: %14.6e\n", *a, *alpha, *phi, *dphi) ; + } + if ( cg_Wolfe (*alpha, *phi, *dphi, Com) ) + { + status = 0 ; + goto Exit2 ; + } + if ( *dphi >= ZERO ) + { + *b = *alpha ; + *dphib = *dphi ; + goto Exit1 ; + } + if ( *phi <= Com->fpert ) + { + if ( Parm->PrintLevel >= 2 ) + { + printf ("update a: %14.6e dphia: %14.6e\n", *alpha, *dphi) ; + } + *a = *alpha ; + *dphia = *dphi ; + } + else *b = *alpha ; + } +Exit1: + status = -1 ; +Exit2: + if ( Parm->PrintLevel >= 2 ) + { + printf ("UP a: %14.6e b: %14.6e da: %14.6e db: %14.6e status: %i\n", + *a, *b, *dphia, *dphib, status) ; + } + return (status) ; +} + +/* ========================================================================= + ==== cg_updateW ========================================================= + ========================================================================= + This routine is identical to cg_update except that the function + psi [a] = phi [a] - phi [0] - a*delta*dphi [0] is minimized instead of + the function phi. The return int has the following meaning: + 8 if too many iterations + 0 if Wolfe condition is satisfied + -1 if interval is updated and a search is done + -2 if the interval updated successfully + ========================================================================= */ +int cg_updateW +( + double *a , /* left side of bracketing interval */ + double *dpsia , /* derivative at a */ + double *b , /* right side of bracketing interval */ + double *dpsib , /* derivative at b */ + double *alpha , /* trial step (between a and b) */ + double *phi , /* function value at alpha (returned) */ + double *dphi , /* derivative of phi at alpha (returned) */ + double *dpsi , /* derivative of psi at alpha (returned) */ + cg_com *Com /* cg com structure */ +) +{ + INT n ; + int nshrink, status ; + double psi, *x, *xtemp, *d, *gtemp ; + cg_parameter *Parm ; + + Parm = Com->Parm ; + n = Com->n ; + x = Com->x ; /* current iterate */ + xtemp = Com->xtemp ; /* x + alpha*d */ + d = Com->d ; /* current search direction */ + gtemp = Com->gtemp ; /* gradient at x + alpha*d */ + cg_step (xtemp, x, d, *alpha, n) ; + *phi = cg_fg (gtemp, xtemp, Com) ; + psi = *phi - *alpha*Com->wolfe_hi ; + *dphi = cg_dot (gtemp, d, n) ; + *dpsi = *dphi - Com->wolfe_hi ; + if ( Parm->PrintLevel >= 2 ) + { + printf ("update alpha: %14.6e psi: %14.6e dpsi: %14.6e\n", + *alpha, psi, *dpsi) ; + } + if ( cg_Wolfe (*alpha, *phi, *dphi, Com) ) + { + status = 0 ; + goto Exit2 ; + } + status = -2 ; + if ( *dpsi >= ZERO ) + { + *b = *alpha ; + *dpsib = *dpsi ; + goto Exit2 ; + } + else + { + if ( psi <= Com->fpert ) + { + *a = *alpha ; + *dpsia = *dpsi ; + goto Exit2 ; + } + } + nshrink = 0 ; + *b = *alpha ; + while ( TRUE ) + { + *alpha = .5*(*a + *b) ; + nshrink++ ; + if ( nshrink > Parm->nexpand ) + { + status = 8 ; + goto Exit2 ; + } + cg_step (xtemp, x, d, *alpha, n) ; + *phi = cg_fg (gtemp, xtemp, Com) ; + *dphi = cg_dot (gtemp, d, n) ; + *dpsi = *dphi - Com->wolfe_hi ; + psi = *phi - *alpha*Com->wolfe_hi ; + if ( Parm->PrintLevel >= 2 ) + { + printf ("contract, a: %14.6e alpha: %14.6e " + "phi: %14.6e dphi: %14.6e\n", *a, *alpha, *phi, *dphi) ; + } + if ( cg_Wolfe (*alpha, *phi, *dphi, Com) ) + { + status = 0 ; + goto Exit2 ; + } + if ( *dpsi >= ZERO ) + { + *b = *alpha ; + *dpsib = *dpsi ; + goto Exit1 ; + } + if ( psi <= Com->fpert ) + { + if ( Parm->PrintLevel >= 2 ) + { + printf ("update a: %14.6e dpsia: %14.6e\n", *alpha, *dpsi) ; + } + *a = *alpha ; + *dpsia = *dpsi ; + } + else *b = *alpha ; + } +Exit1: + status = -1 ; +Exit2: + if ( Parm->PrintLevel >= 2 ) + { + printf ("UP a: %14.6e b: %14.6e da: %14.6e db: %14.6e status: %i\n", + *a, *b, *dpsia, *dpsib, status) ; + } + return (status) ; +} + +/* ========================================================================= + ==== cg_printParms ====================================================== + ========================================================================= + Print the contents of the cg_parameter structure + ========================================================================= */ +void cg_printParms +( + cg_parameter *Parm +) +{ + printf ("PARAMETERS:\n") ; + printf ("\n") ; + printf ("Wolfe line search parameter ..................... delta: %e\n", + Parm->delta) ; + printf ("Wolfe line search parameter ..................... sigma: %e\n", + Parm->sigma) ; + printf ("decay factor for bracketing interval ............ gamma: %e\n", + Parm->gamma) ; + printf ("growth factor for bracket interval ................ rho: %e\n", + Parm->rho) ; + printf ("growth factor for bracket interval after nan .. nan_rho: %e\n", + Parm->nan_rho) ; + printf ("truncation factor for cg beta ..................... eta: %e\n", + Parm->eta) ; + printf ("perturbation parameter for function value ......... eps: %e\n", + Parm->eps) ; + printf ("factor for computing average cost .............. Qdecay: %e\n", + Parm->Qdecay) ; + printf ("relative change in cost to stop quadstep ... QuadCufOff: %e\n", + Parm->QuadCutOff) ; + printf ("factor multiplying gradient in stop condition . StopFac: %e\n", + Parm->StopFac) ; + printf ("cost change factor, approx Wolfe transition . AWolfeFac: %e\n", + Parm->AWolfeFac) ; + printf ("restart cg every restart_fac*n iterations . restart_fac: %e\n", + Parm->restart_fac) ; + printf ("stop when cost change <= feps*|f| ................. eps: %e\n", + Parm->eps) ; + printf ("starting guess parameter in first iteration ...... psi0: %e\n", + Parm->psi0) ; + printf ("starting step in first iteration if nonzero ...... step: %e\n", + Parm->step) ; + printf ("factor multiply starting guess in quad step ...... psi1: %e\n", + Parm->psi1) ; + printf ("initial guess factor for general iteration ....... psi2: %e\n", + Parm->psi2) ; + printf ("max iterations is n*maxit_fac ............... maxit_fac: %e\n", + Parm->maxit_fac) ; + printf ("max expansions in line search ................. nexpand: %i\n", + Parm->nexpand) ; + printf ("max secant iterations in line search .......... nsecant: %i\n", + Parm->nsecant) ; + printf ("print level (0 = none, 2 = maximum) ........ PrintLevel: %i\n", + Parm->PrintLevel) ; + printf ("Logical parameters:\n") ; + if ( Parm->PertRule ) + printf (" Error estimate for function value is eps\n") ; + else + printf (" Error estimate for function value is eps*Ck\n") ; + if ( Parm->QuadStep ) + printf (" Use quadratic interpolation step\n") ; + else + printf (" No quadratic interpolation step\n") ; + if ( Parm->PrintFinal ) + printf (" Print final cost and statistics\n") ; + else + printf (" Do not print final cost and statistics\n") ; + if ( Parm->PrintParms ) + printf (" Print the parameter structure\n") ; + else + printf (" Do not print parameter structure\n") ; + if ( Parm->AWolfe) + printf (" Approximate Wolfe line search\n") ; + else + printf (" Wolfe line search") ; + if ( Parm->AWolfeFac > 0. ) + printf (" ... switching to approximate Wolfe\n") ; + else + printf ("\n") ; + if ( Parm->StopRule ) + printf (" Stopping condition uses initial grad tolerance\n") ; + else + printf (" Stopping condition weighted by absolute cost\n") ; + if ( Parm->debug) + printf (" Check for decay of cost, debugger is on\n") ; + else + printf (" Do not check for decay of cost, debugger is off\n") ; +} + +/* +Version 1.2 Change: + The variable dpsi needs to be included in the argument list for + subroutine cg_updateW (update of a Wolfe line search) + +Version 2.0 Changes: + The user interface was redesigned. The parameters no longer need to + be read from a file. For compatibility with earlier versions of the + code, we include the routine cg_readParms to read parameters. + In the simplest case, the user can use NULL for the + parameter argument of cg_descent, and the code sets the default + parameter values. If the user wishes to modify the parameters, call + cg_default in the main program to initialize a cg_parameter + structure. Individual elements of the structure could be modified. + The header file cg_user.h contains the structures and prototypes + that the user may need to reference or modify, while cg_descent.h + contains header elements that only cg_descent will access. Note + that the arguments of cg_descent have changed. + +Version 3.0 Changes: + Major overhaul +*/ diff --git a/CG_DESCENT-C-3.0/cg_descent.h b/CG_DESCENT-C-3.0/cg_descent.h new file mode 100644 --- /dev/null +++ b/CG_DESCENT-C-3.0/cg_descent.h @@ -0,0 +1,148 @@ +#include +#include +#include +#include +#include +#include +#include + +#define ZERO ((double) 0) +#define ONE ((double) 1) +#define MAX(a,b) (((a) > (b)) ? (a) : (b)) +#define MIN(a,b) (((a) < (b)) ? (a) : (b)) + +typedef struct cg_com_struct /* common variables */ +{ + /* parameters computed by the code */ + INT n ; /* problem dimension, saved for reference */ + INT nf ; /* number of function evaluations */ + INT ng ; /* number of gradient evaluations */ + int QuadOK ; /* T (quadratic step successful) */ + double alpha ; /* stepsize along search direction */ + double f ; /* function value for step alpha */ + double df ; /* function derivative for step alpha */ + double fpert ; /* perturbation is eps*Ck if PertRule is T */ + double f0 ; /* old function value */ + double Ck ; /* average cost as given by the rule: + Qk = Qdecay*Qk + 1, Ck += (fabs (f) - Ck)/Qk */ + double wolfe_hi ; /* upper bound for slope in Wolfe test */ + double wolfe_lo ; /* lower bound for slope in Wolfe test */ + double awolfe_hi ; /* upper bound for slope, approximate Wolfe test */ + int AWolfe ; /* F (use Wolfe line search) + T (use approximate Wolfe line search) + do not change user's AWolfe, this value can be + changed based on AWolfeFac */ + double *x ; /* current iterate */ + double *xtemp ; /* x + alpha*d */ + double *d ; /* current search direction */ + double *gtemp ; /* gradient at x + alpha*d */ + double (*cg_value) (double *, INT) ; /* f = cg_value (x, n) */ + void (*cg_grad) (double *, double *, INT) ; /* cg_grad (g, x, n) */ + double (*cg_valgrad) (double *, double *, INT) ; /* f = cg_valgrad (g,x,n)*/ + cg_parameter *Parm ; /* user parameters */ +} cg_com ; + +/* prototypes */ + +int cg_Wolfe +( + double alpha, /* stepsize */ + double f, /* function value associated with stepsize alpha */ + double dphi, /* derivative value associated with stepsize alpha */ + cg_com *Com /* cg com */ +) ; + +double cg_f +( + double *x, + cg_com *Com +) ; + +void cg_g +( + double *g, + double *x, + cg_com *Com +) ; + +double cg_fg +( + double *g, + double *x, + cg_com *Com +) ; + + +int cg_tol +( + double f, /* function value associated with stepsize */ + double gnorm, /* gradient sup-norm */ + int StopRule, /* T => |grad|_infty <=max (tol, |grad|_infty*StopFact) + F => |grad|_infty <= tol*(1+|f|)) */ + double tol /* tolerance */ +) ; + +double cg_dot +( + double *x, /* first vector */ + double *y, /* second vector */ + INT n /* length of vectors */ +) ; + +void cg_copy +( + double *y, /* output of copy */ + double *x, /* input of copy */ + int n /* length of vectors */ +) ; + +void cg_step +( + double *xtemp, /*output vector */ + double *x, /* initial vector */ + double *d, /* search direction */ + double alpha, /* stepsize */ + INT n /* length of the vectors */ +) ; + +int cg_line +( + double dphi0, /* function derivative at starting point (alpha = 0) */ + cg_com *Com /* cg com structure */ +) ; + +int cg_lineW +( + double dphi0, /* function derivative at starting point (alpha = 0) */ + cg_com *Com /* cg com structure */ +) ; + +int cg_update +( + double *a, /* left side of bracketing interval */ + double *dphia, /* derivative at a */ + double *b, /* right side of bracketing interval */ + double *dphib, /* derivative at b */ + double *alpha, /* trial step (between a and b) */ + double *phi, /* function value at alpha (returned) */ + double *dphi, /* function derivative at alpha (returned) */ + cg_com *Com /* cg com structure */ +) ; + +int cg_updateW +( + double *a, /* left side of bracketing interval */ + double *dpsia, /* derivative at a */ + double *b, /* right side of bracketing interval */ + double *dpsib, /* derivative at b */ + double *alpha, /* trial step (between a and b) */ + double *phi, /* function value at alpha (returned) */ + double *dphi, /* derivative of phi at alpha (returned) */ + double *dpsi, /* derivative of psi at alpha (returned) */ + cg_com *Com /* cg com structure */ +) ; + +void cg_printParms +( + cg_parameter *Parm +) ; diff --git a/CG_DESCENT-C-3.0/cg_user.h b/CG_DESCENT-C-3.0/cg_user.h new file mode 100644 --- /dev/null +++ b/CG_DESCENT-C-3.0/cg_user.h @@ -0,0 +1,164 @@ +#include +#include +#include +#include + +#define INT long int +#define INT_INF LONG_MAX +#define INF DBL_MAX + +#ifndef FALSE +#define FALSE 0 +#endif + +#ifndef TRUE +#define TRUE 1 +#endif + +#ifndef NULL +#define NULL 0 +#endif + +/*============================================================================ + user controlled parameters for the conjugate gradient algorithm + (default values in cg_default) */ +typedef struct cg_parameter_struct /* user controlled parameters */ +{ + /* parameters values that the user may wish to modify */ +/*----------------------------------------------------------------------------*/ + /* T => print final statistics + F => no printout of statistics */ + int PrintFinal ; + + /* Level 0 = no printing), ... , Level 3 = maximum printing */ + int PrintLevel ; + + /* T => print parameters values + F => do not display parmeter values */ + int PrintParms ; + + /* T => use approximate Wolfe line search + F => use ordinary Wolfe line search, switch to approximate Wolfe when + |f_k+1-f_k| < AWolfeFac*C_k, C_k = average size of cost */ + int AWolfe ; + double AWolfeFac ; + + /* factor in [0, 1] used to compute average cost magnitude C_k as follows: + Q_k = 1 + (Qdecay)Q_k-1, Q_0 = 0, C_k = C_k-1 + (|f_k| - C_k-1)/Q_k */ + double Qdecay ; + + /* Stop Rules: + T => ||proj_grad||_infty <= max(grad_tol,initial ||grad||_infty*StopFact) + F => ||proj_grad||_infty <= grad_tol*(1 + |f_k|) */ + int StopRule ; + double StopFac ; + + /* T => estimated error in function value is eps*Ck, + F => estimated error in function value is eps */ + int PertRule ; + double eps ; + + /* T => attempt quadratic interpolation in line search when + |f_k+1 - f_k|/f_k <= QuadCutoff + F => no quadratic interpolation step */ + int QuadStep ; + double QuadCutOff ; + + /* T => check that f_k+1 - f_k <= debugtol*C_k + F => no checking of function values */ + int debug ; + double debugtol ; + + /* if step is nonzero, it is the initial step of the initial line search */ + double step ; + + /* abort cg after maxit_fac*n iterations */ + double maxit_fac ; + + /* maximum number of times the bracketing interval grows or shrinks + in the line search is nexpand */ + int nexpand ; + + /* maximum number of secant iterations in line search is nsecant */ + int nsecant ; + + /* conjugate gradient method restarts after (n*restart_fac) iterations */ + double restart_fac ; + + /* stop when -alpha*dphi0 (estimated change in function value) <= feps*|f|*/ + double feps ; + + /* after encountering nan, growth factor when searching for + a bracketing interval */ + double nan_rho ; + +/*============================================================================ + technical parameters which the user probably should not touch */ + double delta ; /* Wolfe line search parameter */ + double sigma ; /* Wolfe line search parameter */ + double gamma ; /* decay factor for bracket interval width */ + double rho ; /* growth factor when searching for initial + bracketing interval */ + double eta ; /* lower bound for the conjugate gradient update + parameter beta_k is eta*||d||_2 */ + double psi0 ; /* factor used in starting guess for iteration 1 */ + double psi1 ; /* in performing a QuadStep, we evaluate the + function at psi1*previous step */ + double psi2 ; /* when starting a new cg iteration, our initial + guess for the line search stepsize is + psi2*previous step */ +} cg_parameter ; + +typedef struct cg_stats_struct /* statistics returned to user */ +{ + double f ; /*function value at solution */ + double gnorm ; /* max abs component of gradient */ + INT iter ; /* number of iterations */ + INT nfunc ; /* number of function evaluations */ + INT ngrad ; /* number of gradient evaluations */ +} cg_stats ; + +/* prototypes */ + +int cg_descent /* return: + -2 (function value became nan) + -1 (starting function value is nan) + 0 (convergence tolerance satisfied) + 1 (change in func <= feps*|f|) + 2 (total iterations exceeded maxit) + 3 (slope always negative in line search) + 4 (number secant iterations exceed nsecant) + 5 (search direction not a descent direction) + 6 (line search fails in initial interval) + 7 (line search fails during bisection) + 8 (line search fails during interval update) + 9 (debugger is on and the function value increases) + 10 (out of memory) */ +( + double *x, /* input: starting guess, output: the solution */ + INT n, /* problem dimension */ + cg_stats *Stats, /* structure with statistics (see cg_descent.h) */ + cg_parameter *UParm, /* user parameters, NULL = use default parameters */ + double grad_tol, /* StopRule = 1: |g|_infty <= max (grad_tol, + StopFac*initial |g|_infty) [default] + StopRule = 0: |g|_infty <= grad_tol(1+|f|) */ + double (*value) (double *, INT), /* f = value (x, n) */ + void (*grad) (double *, double *, INT), /* grad (g, x, n) */ + double (*valgrad) (double *, double *, INT), /* f = valgrad (g,x,n)*/ + double *Work /* either size 4n work array or NULL */ +) ; + +void cg_default /* set default parameter values */ +( + cg_parameter *Parm +) ; + +int cg_readParms /* return: + 0 (parameter file was read) + -1 (parameter file not found) + -2 (missing entry in parameter file) + -3 (comment in parameter file too long) */ +( + char *filename , + cg_parameter *Parm +) ; diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,674 @@ + GNU GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The GNU General Public License is a free, copyleft license for +software and other kinds of works. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +the GNU General Public License is intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. We, the Free Software Foundation, use the +GNU General Public License for most of our software; it applies also to +any other work released this way by its authors. You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things. + + To protect your rights, we need to prevent others from denying you +these rights or asking you to surrender the rights. Therefore, you have +certain responsibilities if you distribute copies of the software, or if +you modify it: responsibilities to respect the freedom of others. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must pass on to the recipients the same +freedoms that you received. You must make sure that they, too, receive +or can get the source code. And you must show them these terms so they +know their rights. + + Developers that use the GNU GPL protect your rights with two steps: +(1) assert copyright on the software, and (2) offer you this License +giving you legal permission to copy, distribute and/or modify it. + + For the developers' and authors' protection, the GPL clearly explains +that there is no warranty for this free software. For both users' and +authors' sake, the GPL requires that modified versions be marked as +changed, so that their problems will not be attributed erroneously to +authors of previous versions. + + Some devices are designed to deny users access to install or run +modified versions of the software inside them, although the manufacturer +can do so. This is fundamentally incompatible with the aim of +protecting users' freedom to change the software. The systematic +pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we +have designed this version of the GPL to prohibit the practice for those +products. If such problems arise substantially in other domains, we +stand ready to extend this provision to those domains in future versions +of the GPL, as needed to protect the freedom of users. + + Finally, every program is threatened constantly by software patents. +States should not allow patents to restrict development and use of +software on general-purpose computers, but in those that do, we wish to +avoid the special danger that patents applied to a free program could +make it effectively proprietary. To prevent this, the GPL assures that +patents cannot be used to render the program non-free. + + The precise terms and conditions for copying, distribution and +modification follow. + + TERMS AND CONDITIONS + + 0. Definitions. + + "This License" refers to version 3 of the GNU General Public License. + + "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + + "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations. + + To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work. + + A "covered work" means either the unmodified Program or a work based +on the Program. + + To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well. + + To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying. + + An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + + 1. Source Code. + + The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work. + + A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language. + + The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it. + + The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work. + + The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source. + + The Corresponding Source for a work in source code form is that +same work. + + 2. Basic Permissions. + + All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law. + + You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you. + + Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary. + + 3. Protecting Users' Legal Rights From Anti-Circumvention Law. + + No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures. + + When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures. + + 4. Conveying Verbatim Copies. + + You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program. + + You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee. + + 5. Conveying Modified Source Versions. + + You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions: + + a) The work must carry prominent notices stating that you modified + it, and giving a relevant date. + + b) The work must carry prominent notices stating that it is + released under this License and any conditions added under section + 7. This requirement modifies the requirement in section 4 to + "keep intact all notices". + + c) You must license the entire work, as a whole, under this + License to anyone who comes into possession of a copy. This + License will therefore apply, along with any applicable section 7 + additional terms, to the whole of the work, and all its parts, + regardless of how they are packaged. This License gives no + permission to license the work in any other way, but it does not + invalidate such permission if you have separately received it. + + d) If the work has interactive user interfaces, each must display + Appropriate Legal Notices; however, if the Program has interactive + interfaces that do not display Appropriate Legal Notices, your + work need not make them do so. + + A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate. + + 6. Conveying Non-Source Forms. + + You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways: + + a) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by the + Corresponding Source fixed on a durable physical medium + customarily used for software interchange. + + b) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by a + written offer, valid for at least three years and valid for as + long as you offer spare parts or customer support for that product + model, to give anyone who possesses the object code either (1) a + copy of the Corresponding Source for all the software in the + product that is covered by this License, on a durable physical + medium customarily used for software interchange, for a price no + more than your reasonable cost of physically performing this + conveying of source, or (2) access to copy the + Corresponding Source from a network server at no charge. + + c) Convey individual copies of the object code with a copy of the + written offer to provide the Corresponding Source. This + alternative is allowed only occasionally and noncommercially, and + only if you received the object code with such an offer, in accord + with subsection 6b. + + d) Convey the object code by offering access from a designated + place (gratis or for a charge), and offer equivalent access to the + Corresponding Source in the same way through the same place at no + further charge. You need not require recipients to copy the + Corresponding Source along with the object code. If the place to + copy the object code is a network server, the Corresponding Source + may be on a different server (operated by you or a third party) + that supports equivalent copying facilities, provided you maintain + clear directions next to the object code saying where to find the + Corresponding Source. Regardless of what server hosts the + Corresponding Source, you remain obligated to ensure that it is + available for as long as needed to satisfy these requirements. + + e) Convey the object code using peer-to-peer transmission, provided + you inform other peers where the object code and Corresponding + Source of the work are being offered to the general public at no + charge under subsection 6d. + + A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work. + + A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + + "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made. + + If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM). + + The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network. + + Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying. + + 7. Additional Terms. + + "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + + When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission. + + Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms: + + a) Disclaiming warranty or limiting liability differently from the + terms of sections 15 and 16 of this License; or + + b) Requiring preservation of specified reasonable legal notices or + author attributions in that material or in the Appropriate Legal + Notices displayed by works containing it; or + + c) Prohibiting misrepresentation of the origin of that material, or + requiring that modified versions of such material be marked in + reasonable ways as different from the original version; or + + d) Limiting the use for publicity purposes of names of licensors or + authors of the material; or + + e) Declining to grant rights under trademark law for use of some + trade names, trademarks, or service marks; or + + f) Requiring indemnification of licensors and authors of that + material by anyone who conveys the material (or modified versions of + it) with contractual assumptions of liability to the recipient, for + any liability that these contractual assumptions directly impose on + those licensors and authors. + + All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying. + + If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms. + + Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way. + + 8. Termination. + + You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11). + + However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + + Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + + Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10. + + 9. Acceptance Not Required for Having Copies. + + You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so. + + 10. Automatic Licensing of Downstream Recipients. + + Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License. + + An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts. + + You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it. + + 11. Patents. + + A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version". + + A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License. + + Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version. + + In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party. + + If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid. + + If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it. + + A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007. + + Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law. + + 12. No Surrender of Others' Freedom. + + If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program. + + 13. Use with the GNU Affero General Public License. + + Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU Affero General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the special requirements of the GNU Affero General Public License, +section 13, concerning interaction through a network will apply to the +combination as such. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + + Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU General Public License, you may choose any version ever published +by the Free Software Foundation. + + If the Program specifies that a proxy can decide which future +versions of the GNU General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program. + + Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version. + + 15. Disclaimer of Warranty. + + THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. Limitation of Liability. + + IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES. + + 17. Interpretation of Sections 15 and 16. + + If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + +Also add information on how to contact you by electronic and paper mail. + + If the program does terminal interaction, make it output a short +notice like this when it starts in an interactive mode: + + Copyright (C) + This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, your program's commands +might be different; for a GUI interface, you would use an "about box". + + You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU GPL, see +. + + The GNU General Public License does not permit incorporating your program +into proprietary programs. If your program is a subroutine library, you +may consider it more useful to permit linking proprietary applications with +the library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. But first, please read +. diff --git a/Numeric/Optimization/Algorithms/HagerZhang05.hsc b/Numeric/Optimization/Algorithms/HagerZhang05.hsc new file mode 100644 --- /dev/null +++ b/Numeric/Optimization/Algorithms/HagerZhang05.hsc @@ -0,0 +1,673 @@ +--------------------------------------------------------------------------- +-- | Module : Numeric.Statistics.Dirichlet.Mixture +-- Copyright : (c) 2009 Felipe Lessa +-- License : GPL +-- +-- Maintainer : felipe.lessa@gmail.com +-- Stability : experimental +-- Portability : portable +-- +-- This module implements the algorithms described by Hager and +-- Zhang [1]. We use bindings to @CG_DESCENT@ library by the same +-- authors, version 3.0 from 18/05/2008 [2]. The library code is +-- also licensed under the terms of the GPL. +-- +-- * [1] Hager, W. W. and Zhang, H. /A new conjugate gradient/ +-- /method with guaranteed descent and an efficient line/ +-- /search./ Society of Industrial and Applied Mathematics +-- Journal on Optimization, 16 (2005), 170-192. +-- +-- * [2] +-- +-------------------------------------------------------------------------- + + +module Math.Optimization.Algorithms.HagerZhang05 + (-- * Main function + -- $mainFunction + optimize + -- ** User-defined function types + ,Function(..) + ,Gradient(..) + ,Combined(..) + -- ** Kinds of function types + ,Simple + ,Mutable + -- * Result and statistics + ,Result(..) + ,Statistics(..) + -- * Options + ,defaultParameters + ,Parameters(..) + ,Verbose(..) + ,LineSearch(..) + ,StopRules(..) + ,EstimateError(..) + -- * Technical parameters + ,TechParameters(..) + ) where + +import qualified Data.Vector.Generic as G +import qualified Data.Vector.Generic.Mutable as GM +import qualified Data.Vector.Storable as S +import qualified Data.Vector.Storable.Mutable as SM +import Control.Exception (bracket) +import Control.Monad.Primitive (PrimMonad(..)) +import Foreign +import Foreign.C +#include "cg_user.h" + +-- $mainFunction +-- Please pay close attention to the types of @Vector@s and +-- @MVetor@s being used below. They may come from +-- "Data.Vector.Generic"/"Data.Vector.Generic.Mutable" or from +-- "Data.Vector.Storable"/"Data.Vector.Storable.Mutable". The +-- rule of thumb is that input pure vectors are @Generic@ and +-- everything else is @Storable@. + + +-- | Run the @CG_DESCENT@ optimizer and try to minimize the +-- function. +optimize :: (G.Vector v Double) + => Parameters -- ^ How should we optimize. + -> Double -- ^ @grad_tol@, see 'stopRules'. + -> v Double -- ^ Initial guess. + -> Function t1 -- ^ Function to be minimized. + -> Gradient t2 -- ^ Gradient of the function. + -> Maybe (Combined t3) -- ^ (Optional) Combined function computing + -- both the function and its gradient. + -> IO (S.Vector Double, Result, Statistics) +optimize params grad_tol initial f g c = do + -- Mutable vector used for initial guess and final solution. + let n = G.length initial + x <- GM.unstream $ G.stream initial + + -- Convert user-provided functions. + let mf = mutableF f + mg = mutableG g + mc = maybe (combine mf mg) mutableC c + cf = prepareF mf + cg = prepareG mg + cc = prepareC mc + + -- Allocate everything. + (ret, stats) <- + SM.unsafeWith x $ \x_ptr -> + alloca $ \stats_ptr -> + alloca $ \param_ptr -> + bracket (mkCFunction cf) freeHaskellFunPtr $ \cf_ptr -> + bracket (mkCGradient cg) freeHaskellFunPtr $ \cg_ptr -> + bracket (mkCCombined cc) freeHaskellFunPtr $ \cc_ptr -> + allocaArray (4*n) $ \work_ptr -> do + -- Go to C land. + poke param_ptr params + ret <- cg_descent x_ptr (fromIntegral n) + stats_ptr param_ptr grad_tol + cf_ptr cg_ptr cc_ptr work_ptr + stats <- peek stats_ptr + return (intToResult ret, stats) + + -- Retrive solution and return. + x' <- G.unsafeFreeze x + return $ ret `seq` (x', ret, stats) + +type CFunction = Ptr Double -> CInt -> IO Double +type CGradient = Ptr Double -> Ptr Double -> CInt -> IO () +type CCombined = Ptr Double -> Ptr Double -> CInt -> IO Double +foreign import ccall safe "cg_user.h" + cg_descent :: Ptr Double + -> CInt + -> Ptr Statistics + -> Ptr Parameters + -> Double + -> FunPtr CFunction + -> FunPtr CGradient + -> FunPtr CCombined + -> Ptr Double + -> IO CInt +foreign import ccall "wrapper" mkCFunction :: CFunction -> IO (FunPtr CFunction) +foreign import ccall "wrapper" mkCGradient :: CGradient -> IO (FunPtr CGradient) +foreign import ccall "wrapper" mkCCombined :: CCombined -> IO (FunPtr CCombined) + + +-- | Phantom type for simple pure functions. +data Simple +-- | Phantom type for functions using mutable data. +data Mutable + +-- | Function calculating the value of the objective function @f@ +-- at a point @x@. +data Function t where + VFunction :: G.Vector v Double + => (v Double -> Double) + -> Function Simple + MFunction :: (forall m. PrimMonad m + => SM.MVector (PrimState m) Double + -> m Double) + -> Function Mutable + +mutableF :: Function t -> Function Mutable +mutableF (VFunction f) = MFunction f' + where + f' mx = do + -- Copy the input to an immutable vector. + let s = GM.length mx + mz <- GM.new s + let go i | i > s = return () + | otherwise = GM.unsafeRead mx i >>= + GM.unsafeWrite mz i >> go (i+1) + go 0 + z <- G.unsafeFreeze mz + -- Run the user function. + return (f z) +mutableF (MFunction f) = MFunction f + +prepareF :: Function Mutable -> CFunction +prepareF (MFunction f) = + \x_ptr n -> do + let n' = fromIntegral n + x_fptr <- newForeignPtr_ x_ptr + f (SM.unsafeFromForeignPtr x_fptr 0 n') +prepareF _ = error "HagerZhang05.prepareF: never here" + + + + + +-- | Function calculating the value of the gradient of the +-- objective function @f@ at a point @x@. +-- +-- The 'MGradient' constructor uses a function receiving as +-- parameters the point @x@ being evaluated (should not be +-- modified) and the vector where the gradient should be written. +data Gradient t where + VGradient :: G.Vector v Double + => (v Double -> v Double) + -> Gradient Simple + MGradient :: (forall m. PrimMonad m + => SM.MVector (PrimState m) Double + -> SM.MVector (PrimState m) Double + -> m ()) + -> Gradient Mutable +mutableG :: Gradient t -> Gradient Mutable +mutableG (VGradient f) = MGradient f' + where + f' mx mret = do + -- Copy the input to an immutable vector. + let s = GM.length mx + mz <- GM.new s + let go i | i > s = return () + | otherwise = GM.unsafeRead mx i >>= + GM.unsafeWrite mz i >> go (i+1) + go 0 + z <- G.unsafeFreeze mz + -- Run the user function. + let !r = f z + -- Copy the output to an immutable vector + let s' = min s (G.length r) + go' i | i > s' = return () + | otherwise = let !x = G.unsafeIndex r i + in GM.unsafeWrite mret i x >> go (i+1) + go' 0 +mutableG (MGradient f) = MGradient f + +prepareG :: Gradient Mutable -> CGradient +prepareG (MGradient f) = + \ret_ptr x_ptr n -> do + let n' = fromIntegral n + x_fptr <- newForeignPtr_ x_ptr + ret_fptr <- newForeignPtr_ ret_ptr + f (SM.unsafeFromForeignPtr x_fptr 0 n') + (SM.unsafeFromForeignPtr ret_fptr 0 n') +prepareG _ = error "HagerZhang05.prepareG: never here" + + + + + + + + + +-- | Function calculating the both the value of the objective +-- function @f@ and its gradient at a point @x@. +data Combined t where + VCombined :: G.Vector v Double + => (v Double -> (Double, v Double)) + -> Combined Simple + MCombined :: (forall m. PrimMonad m + => SM.MVector (PrimState m) Double + -> SM.MVector (PrimState m) Double + -> m Double) + -> Combined Mutable +mutableC :: Combined t -> Combined Mutable +mutableC (VCombined f) = MCombined f' + where + f' mx mret = do + -- Copy the input to an immutable vector. + let s = GM.length mx + mz <- GM.new s + let go i | i > s = return () + | otherwise = GM.unsafeRead mx i >>= + GM.unsafeWrite mz i >> go (i+1) + go 0 + z <- G.unsafeFreeze mz + -- Run the user function. + let !(v,r) = f z + -- Copy the output to an immutable vector + let s' = min s (G.length r) + go' i | i > s' = return () + | otherwise = let !x = G.unsafeIndex r i + in GM.unsafeWrite mret i x >> go (i+1) + go' 0 + -- Return the value + return v +mutableC (MCombined f) = MCombined f + +prepareC :: Combined Mutable -> CCombined +prepareC (MCombined f) = + \ret_ptr x_ptr n -> do + let n' = fromIntegral n + x_fptr <- newForeignPtr_ x_ptr + ret_fptr <- newForeignPtr_ ret_ptr + f (SM.unsafeFromForeignPtr x_fptr 0 n') + (SM.unsafeFromForeignPtr ret_fptr 0 n') +prepareC _ = error "HagerZhang05.prepareC: never here" + +-- | Combine two separated functions into a single, combined one. +-- This is always a win for us since we save one jump from C to +-- Haskell land. +combine :: Function Mutable -> Gradient Mutable -> Combined Mutable +combine (MFunction f) (MGradient g) = + MCombined $ \mx mret -> g mx mret >> f mx +combine _ _ = error "HagerZhang05.combine: never here" + + + + +data Result = + ToleranceStatisfied + -- ^ Convergence tolerance was satisfied. + | FunctionChange + -- ^ Change in function value was less than @funcEpsilon * + -- |f|@. + | MaxTotalIter + -- ^ Total iterations exceeded @maxItersFac * n@. + | NegativeSlope + -- ^ Slope was always negative in line search. + | MaxSecantIter + -- ^ Number of secant iterations exceed nsecant. + | NotDescent + -- ^ Search direction not a descent direction. + | LineSearchFailsInitial + -- ^ Line search fails in initial interval. + | LineSearchFailsBisection + -- ^ Line search fails during bisection. + | LineSearchFailsUpdate + -- ^ Line search fails during interval update. + | DebugTol + -- ^ Debug tolerance was on and the test failed (see 'debugTol'). + | OutOfMemory + -- ^ Couldn't allocate enought temporary memory. + | FunctionValueNaN + -- ^ Function value became @NaN@. + | StartFunctionValueNaN + -- ^ Initial function value was @NaN@. + deriving (Eq, Ord, Show, Read, Enum) + +intToResult :: CInt -> Result +intToResult (-2) = FunctionValueNaN +intToResult (-1) = StartFunctionValueNaN +intToResult 0 = ToleranceStatisfied +intToResult 1 = FunctionChange +intToResult 2 = MaxTotalIter +intToResult 3 = NegativeSlope +intToResult 4 = MaxSecantIter +intToResult 5 = NotDescent +intToResult 6 = LineSearchFailsInitial +intToResult 7 = LineSearchFailsBisection +intToResult 8 = LineSearchFailsUpdate +intToResult 9 = DebugTol +intToResult 10 = error $ "HagerZhang05.intToResult: out of memory?! how?!" +intToResult x = error $ "HagerZhang05.intToResult: unknown value " ++ show x + +-- | Statistics given after the process finishes. +data Statistics = Statistics { + finalValue :: Double + -- ^ Value of the function at the solution. + ,gradNorm :: Double + -- ^ Maximum absolute component of the gradient at the + -- solution. + ,totalIters :: CInt + -- ^ Total number of iterations. + ,funcEvals :: CInt + -- ^ Total number of function evaluations. + ,gradEvals :: CInt + -- ^ Total number of gradient evaluations. + } deriving (Eq, Ord, Show, Read) + +instance Storable Statistics where + sizeOf _ = #{size cg_stats} + alignment _ = alignment (undefined :: Double) + peek ptr = do + v_finalValue <- #{peek cg_stats, f} ptr + v_gradNorm <- #{peek cg_stats, gnorm} ptr + v_totalIters <- #{peek cg_stats, iter} ptr + v_funcEvals <- #{peek cg_stats, nfunc} ptr + v_gradEvals <- #{peek cg_stats, ngrad} ptr + return Statistics {finalValue = v_finalValue + ,gradNorm = v_gradNorm + ,totalIters = v_totalIters + ,funcEvals = v_funcEvals + ,gradEvals = v_gradEvals} + poke ptr s = do + #{poke cg_stats, f} ptr (finalValue s) + #{poke cg_stats, gnorm} ptr (gradNorm s) + #{poke cg_stats, iter} ptr (totalIters s) + #{poke cg_stats, nfunc} ptr (funcEvals s) + #{poke cg_stats, ngrad} ptr (gradEvals s) + + + +-- | Default parameters. See the documentation for 'Parameters' +-- and 'TechParameters' to see what are the defaults. +defaultParameters :: Parameters +defaultParameters = + unsafePerformIO $ do + alloca $ \ptr -> do + cg_default ptr + peek ptr +{-# NOINLINE defaultParameters #-} +foreign import ccall unsafe "cg_user.h" + cg_default :: Ptr Parameters -> IO () + + +-- | Parameters given to the optimizer. +data Parameters = Parameters { + printFinal :: Bool + -- ^ Print final statistics to @stdout@. Defaults to @True@. + + ,printParams :: Bool + -- ^ Print parameters to @stdout@ before starting. Defaults to @False@ + + ,verbose :: Verbose + -- ^ How verbose we should be while computing. Everything is + -- printed to @stdout@. Defaults to 'Quiet'. + + ,lineSearch :: LineSearch + -- ^ What kind of line search should be used. Defaults to + -- @AutoSwitch 1e-3@. + + ,qdecay :: Double + -- ^ Factor in @[0, 1]@ used to compute average cost + -- magnitude @C_k@ as follows: + -- + -- > Q_k = 1 + (qdecay)Q_{k-1}, Q_0 = 0 + -- > C_k = C_{k-1} + (|f_k| - C_{k-1})/Q_k + -- + -- Defaults to @0.7@. + + ,stopRules :: StopRules + -- ^ Stop rules that define when the iterations should end. + -- Defaults to @DefaultStopRule 0@. + + ,estimateError :: EstimateError + -- ^ How to calculate the estimated error in the function + -- value. Defaults to @RelativeEpsilon 1e-6@. + + ,quadraticStep :: Maybe Double + -- ^ When to attempt quadratic interpolation in line search. + -- If @Nothing@ then never try a quadratic interpolation + -- step. If @Just cutoff@, then attemp quadratic + -- interpolation in line search when @|f_{k+1} - f_k| / f_k + -- <= cutoff@. Defaults to @Just 1e-12@. + + ,debugTol :: Maybe Double + -- ^ If @Just tol@, then always check that @f_{k+1} - f_k <= + -- tol * C_k@. Otherwise, if @Nothing@ then no checking of + -- function values is done. Defaults to @Nothing@. + + ,initialStep :: Maybe Double + -- ^ If @Just step@, then use @step@ as the initial step of + -- the line search. Otherwise, if @Nothing@ then the initial + -- step is programatically calculated. Defaults to + -- @Nothing@. + + ,maxItersFac :: Double + -- ^ Defines the maximum number of iterations. The process + -- is aborted when @maxItersFac * n@ iterations are done, where + -- @n@ is the number of dimensions. Defaults to infinity. + + ,nexpand :: CInt + -- ^ Maximum number of times the bracketing interval grows or + -- shrinks in the line search. Defaults to @50@. + + ,nsecant :: CInt + -- ^ Maximum number of secant iterations in line search. + -- Defaults to @50@. + + ,restartFac :: Double + -- ^ Restart the conjugate gradient method after @restartFac + -- * n@ iterations. Defaults to @1@. + + ,funcEpsilon :: Double + -- ^ Stop when @-alpha * dphi0@, the estimated change in + -- function value, is less than @funcEpsilon * |f|@. + -- Defaults to @0@. + + ,nanRho :: Double + -- ^ After encountering @NaN@ while calculating the step + -- length, growth factor when searching for a bracketing + -- interval. Defaults to @1.3@. + + ,techParameters :: TechParameters + -- ^ Technical parameters which you probably should not + -- touch. + } deriving (Eq, Ord, Show, Read) + +instance Storable Parameters where + sizeOf _ = #{size cg_parameter} + alignment _ = alignment (undefined :: Double) + peek ptr = do + v_printFinal <- #{peek cg_parameter, PrintFinal} ptr + v_printParams <- #{peek cg_parameter, PrintParms} ptr + v_verbose <- #{peek cg_parameter, PrintLevel} ptr + v_awolfe <- #{peek cg_parameter, AWolfe} ptr + v_awolfefac <- #{peek cg_parameter, AWolfeFac} ptr + v_qdecay <- #{peek cg_parameter, Qdecay} ptr + v_stopRule <- #{peek cg_parameter, StopRule} ptr + v_stopRuleFac <- #{peek cg_parameter, StopFac} ptr + v_estimateError <- #{peek cg_parameter, PertRule} ptr + v_estimateEps <- #{peek cg_parameter, eps} ptr + v_quadraticStep <- #{peek cg_parameter, QuadStep} ptr + v_quadraticCut <- #{peek cg_parameter, QuadCutOff} ptr + v_debug <- #{peek cg_parameter, debug} ptr + v_debugTol <- #{peek cg_parameter, debugtol} ptr + v_initialStep <- #{peek cg_parameter, step} ptr + v_maxItersFac <- #{peek cg_parameter, maxit_fac} ptr + v_nexpand <- #{peek cg_parameter, nexpand} ptr + v_nsecant <- #{peek cg_parameter, nsecant} ptr + v_restartFac <- #{peek cg_parameter, restart_fac} ptr + v_funcEpsilon <- #{peek cg_parameter, feps} ptr + v_nanRho <- #{peek cg_parameter, nan_rho} ptr + + v_delta <- #{peek cg_parameter, delta} ptr + v_sigma <- #{peek cg_parameter, sigma} ptr + v_gamma <- #{peek cg_parameter, gamma} ptr + v_rho <- #{peek cg_parameter, rho} ptr + v_eta <- #{peek cg_parameter, eta} ptr + v_psi0 <- #{peek cg_parameter, psi0} ptr + v_psi1 <- #{peek cg_parameter, psi1} ptr + v_psi2 <- #{peek cg_parameter, psi2} ptr + + let tech = TechParameters {techDelta = v_delta + ,techSigma = v_sigma + ,techGamma = v_gamma + ,techRho = v_rho + ,techEta = v_eta + ,techPsi0 = v_psi0 + ,techPsi1 = v_psi1 + ,techPsi2 = v_psi2} + + let b :: CInt -> Bool; b = (/= 0) + + return Parameters {printFinal = b v_printFinal + ,printParams = b v_printParams + ,verbose = case v_verbose :: CInt of + 0 -> Quiet + 1 -> Verbose + _ -> VeryVerbose + ,lineSearch = if b v_awolfe + then ApproximateWolfe + else AutoSwitch v_awolfefac + ,qdecay = v_qdecay + ,stopRules = if b v_stopRule + then DefaultStopRule v_stopRuleFac + else AlternativeStopRule + ,estimateError = if b v_estimateError + then RelativeEpsilon v_estimateEps + else AbsoluteEpsilon v_estimateEps + ,quadraticStep = if b v_quadraticStep + then Just v_quadraticCut + else Nothing + ,debugTol = if b v_debug + then Just v_debugTol + else Nothing + ,initialStep = case v_initialStep of + 0 -> Nothing + x -> Just x + ,maxItersFac = v_maxItersFac + ,nexpand = v_nexpand + ,nsecant = v_nsecant + ,restartFac = v_restartFac + ,funcEpsilon = v_funcEpsilon + ,nanRho = v_nanRho + ,techParameters = tech} + poke ptr p = do + let i b = if b p then 1 else (0 :: CInt) + m b = maybe (0 :: CInt) (const 1) (b p) + #{poke cg_parameter, PrintFinal} ptr (i printFinal) + #{poke cg_parameter, PrintParms} ptr (i printParams) + #{poke cg_parameter, PrintLevel} ptr (case verbose p of + Quiet -> 0 :: CInt + Verbose -> 1 + VeryVerbose -> 3) + let (awolfe, awolfefac) = case lineSearch p of + ApproximateWolfe -> (1, 0) + AutoSwitch x -> (0, x) + #{poke cg_parameter, AWolfe} ptr (awolfe :: CInt) + #{poke cg_parameter, AWolfeFac} ptr awolfefac + #{poke cg_parameter, Qdecay} ptr (qdecay p) + let (stopRule, stopRuleFac) = case stopRules p of + DefaultStopRule x -> (1, x) + AlternativeStopRule -> (0, 0) + #{poke cg_parameter, StopRule} ptr (stopRule :: CInt) + #{poke cg_parameter, StopFac} ptr stopRuleFac + let (pertRule, eps) = case estimateError p of + RelativeEpsilon x -> (1,x) + AbsoluteEpsilon x -> (0,x) + #{poke cg_parameter, PertRule} ptr (pertRule :: CInt) + #{poke cg_parameter, eps} ptr eps + #{poke cg_parameter, QuadStep} ptr (m quadraticStep) + #{poke cg_parameter, QuadCutOff} ptr (maybe 0 id $ quadraticStep p) + #{poke cg_parameter, debug} ptr (m debugTol) + #{poke cg_parameter, debugtol} ptr (maybe 0 id $ debugTol p) + #{poke cg_parameter, step} ptr (maybe 0 id $ initialStep p) + #{poke cg_parameter, maxit_fac} ptr (maxItersFac p) + #{poke cg_parameter, nexpand} ptr (nexpand p) + #{poke cg_parameter, nsecant} ptr (nsecant p) + #{poke cg_parameter, restart_fac} ptr (restartFac p) + #{poke cg_parameter, feps} ptr (funcEpsilon p) + #{poke cg_parameter, nan_rho} ptr (nanRho p) + + #{poke cg_parameter, delta} ptr (techDelta $ techParameters p) + #{poke cg_parameter, sigma} ptr (techSigma $ techParameters p) + #{poke cg_parameter, gamma} ptr (techGamma $ techParameters p) + #{poke cg_parameter, rho} ptr (techRho $ techParameters p) + #{poke cg_parameter, eta} ptr (techEta $ techParameters p) + #{poke cg_parameter, psi0} ptr (techPsi0 $ techParameters p) + #{poke cg_parameter, psi1} ptr (techPsi1 $ techParameters p) + #{poke cg_parameter, psi2} ptr (techPsi2 $ techParameters p) + + + + +-- | Technical parameters which you probably should not touch. +-- You should read the papers of @CG_DESCENT@ to understand how +-- you can tune these parameters. +data TechParameters = TechParameters { + techDelta :: Double + -- ^ Wolfe line search parameter. Defaults to @0.1@. + ,techSigma :: Double + -- ^ Wolfe line search parameter. Defaults to @0.9@. + ,techGamma :: Double + -- ^ Decay factor for bracket interval width. Defaults to + -- @0.66@. + ,techRho :: Double + -- ^ Growth factor when searching for initial bracketing + -- interval. Defaults to @5@. + ,techEta :: Double + -- ^ Lower bound for the conjugate gradient update parameter + -- @beta_k@ is @techEta * ||d||_2@. Defaults to @0.01@. + ,techPsi0 :: Double + -- ^ Factor used in starting guess for iteration 1. Defaults + -- to @0.01@. + ,techPsi1 :: Double + -- ^ In performing a QuadStep, we evaluate the function at + -- @psi1 * previous step@. Defaults to @0.1@. + ,techPsi2 :: Double + -- ^ When starting a new CG iteration, our initial guess for + -- the line search stepsize is @psi2 * previous step@. + -- Defaults to @2@. + } deriving (Eq, Ord, Show, Read) + + + +-- | How verbose we should be. +data Verbose = + Quiet + -- ^ Do not output anything to @stdout@, which most of the + -- time is good. + | Verbose + -- ^ Print what work is being done on each iteraction. + | VeryVerbose + -- ^ Print information about every step, may be useful for + -- troubleshooting. + deriving (Eq, Ord, Show, Read, Enum) + +-- | Line search methods that may be used. +data LineSearch = + ApproximateWolfe + -- ^ Use approximate Wolfe line search. + | AutoSwitch Double + -- ^ Use ordinary Wolfe line search, switch to approximate + -- Wolfe when + -- + -- > |f_{k+1} - f_k| < AWolfeFac * C_k + -- + -- where @C_k@ is the average size of cost and + -- @AWolfeFac@ is the parameter to this constructor. + deriving (Eq, Ord, Show, Read) + +-- | Stop rules used to decided when to stop iterating. +data StopRules = + DefaultStopRule Double + -- ^ @DefaultStopRule stop_fac@ stops when + -- + -- > |g_k|_infty <= max(grad_tol, |g_0|_infty * stop_fac) + -- + -- where @|g_i|_infty@ is the maximum absolute component of + -- the gradient at the @i@-th step. + | AlternativeStopRule + -- ^ @AlternativeStopRule@ stops when + -- + -- > |g_k|_infty <= grad_tol * (1 + |f_k|) + deriving (Eq, Ord, Show, Read) + +-- | How to calculate the estimated error in the function value. +data EstimateError = + AbsoluteEpsilon Double + -- ^ @AbsoluteEpsilon eps@ estimates the error as @eps@. + | RelativeEpsilon Double + -- ^ @RelativeEpsilon eps@ estimates the error as @eps * C_k@. + deriving (Eq, Ord, Show, Read) diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#!/usr/bin/runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/nonlinear-optimization.cabal b/nonlinear-optimization.cabal new file mode 100644 --- /dev/null +++ b/nonlinear-optimization.cabal @@ -0,0 +1,59 @@ +Cabal-Version: >= 1.2.3 +Build-Type: Simple +Tested-With: GHC +Category: Math +Name: nonlinear-optimization +Version: 0.1 +Stability: experimental +License: GPL +License-File: LICENSE +Copyright: (c) 2010 Felipe A. Lessa and William W. Hager +Author: Felipe A. Lessa (Haskell code), + William W. Hager and Hongchao Zhang (CM_DESCENT code). +Maintainer: Felipe A. Lessa +Synopsis: Various iterative algorithms for optimization of nonlinear functions. +Description: + This library implements numerical algorithms to optimize + nonlinear functions. Optimization means that we try to find + a minimum of the function. Currently all algorithms + guarantee only that local minima will be found, not global + ones. + . + Almost any continuosly differentiable function @f : R^n -> R@ + may be optimized by this library. Any further restrictions + are listed in the modules that need them. + . + We use the @vector@ package to represent vectors and + matrices, although it would be possible to use something like + @hmatrix@ easily. + . + Currently only CM_DESCENT method is implemented. +Extra-Source-Files: + CG_DESCENT-C-3.0/cg_descent.c, + CG_DESCENT-C-3.0/cg_descent.h, + CG_DESCENT-C-3.0/cg_user.h, + CG_DESCENT-C-3.0/README + +Library + Build-Depends: + base >= 3 && < 5, vector >= 0.5 && < 0.6, + primitive >= 0.2 && < 0.3 + Exposed-Modules: + Numeric.Optimization.Algorithms.HagerZhang05 + Include-Dirs: + CG_DESCENT-C-3.0 + C-Sources: + CG_DESCENT-C-3.0/cg_descent.c + Includes: + cg_user.h + Extensions: + CPP, + ForeignFunctionInterface, + EmptyDataDecls, + GADTs, + Rank2Types, + FlexibleContexts + Build-Tools: hsc2hs + Extra-Libraries: m + GHC-Options: -Wall +