nonlinear-optimization (empty) → 0.1
raw patch · 8 files changed
+3553/−0 lines, 8 filesdep +basedep +primitivedep +vectorsetup-changed
Dependencies added: base, primitive, vector
Files
- CG_DESCENT-C-3.0/README +80/−0
- CG_DESCENT-C-3.0/cg_descent.c +1752/−0
- CG_DESCENT-C-3.0/cg_descent.h +148/−0
- CG_DESCENT-C-3.0/cg_user.h +164/−0
- LICENSE +674/−0
- Numeric/Optimization/Algorithms/HagerZhang05.hsc +673/−0
- Setup.lhs +3/−0
- nonlinear-optimization.cabal +59/−0
+ CG_DESCENT-C-3.0/README view
@@ -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"
+ CG_DESCENT-C-3.0/cg_descent.c view
@@ -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+*/
+ CG_DESCENT-C-3.0/cg_descent.h view
@@ -0,0 +1,148 @@+#include <math.h>+#include <limits.h>+#include <float.h>+#include <string.h>+#include <ctype.h>+#include <stdlib.h>+#include <stdio.h>++#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+) ;
+ CG_DESCENT-C-3.0/cg_user.h view
@@ -0,0 +1,164 @@+#include <limits.h>+#include <string.h>+#include <stdlib.h>+#include <stdio.h>++#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+) ;
+ LICENSE view
@@ -0,0 +1,674 @@+ GNU GENERAL PUBLIC LICENSE+ Version 3, 29 June 2007++ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>+ 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.++ <one line to give the program's name and a brief idea of what it does.>+ Copyright (C) <year> <name of author>++ 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 <http://www.gnu.org/licenses/>.++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:++ <program> Copyright (C) <year> <name of author>+ 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+<http://www.gnu.org/licenses/>.++ 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+<http://www.gnu.org/philosophy/why-not-lgpl.html>.
+ Numeric/Optimization/Algorithms/HagerZhang05.hsc view
@@ -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] <http://www.math.ufl.edu/~hager/papers/CG/CG_DESCENT-C-3.0.tar.gz>+--+--------------------------------------------------------------------------+++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)
+ Setup.lhs view
@@ -0,0 +1,3 @@+#!/usr/bin/runhaskell+> import Distribution.Simple+> main = defaultMain
+ nonlinear-optimization.cabal view
@@ -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 <felipe.lessa@gmail.com>+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+