#include <limits.h>
#include <float.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

/* added by nick gould to allow single (float) precision */

#ifdef REAL_32
typedef float    rp_;
#else
typedef double   rp_;
#endif

/*============================================================================
   cg_parameter is a structure containing parameters used in cg_descent
   cg_default assigns default values to these parameters */
typedef struct cg_parameter_struct /* user controlled parameters */
{
/*============================================================================
      parameters 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 LBFGS
       F => only use L-BFGS when memory >= n */
    int LBFGS ;

    /* number of vectors stored in memory */
    int memory ;

    /* SubCheck and SubSkip control the frequency with which the subspace
       condition is checked. It is checked for SubCheck*mem iterations and
       if not satisfied, then it is skipped for Subskip*mem iterations
       and Subskip is doubled. Whenever the subspace condition is statisfied,
       SubSkip is returned to its original value. */
    int SubCheck ;
    int SubSkip ;

    /* when relative distance from current gradient to subspace <= eta0,
       enter subspace if subspace dimension = mem */
    rp_ eta0 ;

    /* when relative distance from current gradient to subspace >= eta1,
       leave subspace */
    rp_ eta1 ;

    /* when relative distance from current direction to subspace <= eta2,
       always enter subspace (invariant space) */
    rp_ eta2 ;

    /* 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 ;
    rp_ 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 */
    rp_ Qdecay ;

    /* terminate after nslow iterations without strict improvement in
       either function value or gradient */
    int nslow ;

    /* Stop Rules:
       T => ||proj_grad||_infty <= max(grad_tol,initial ||grad||_infty*StopFact)
       F => ||proj_grad||_infty <= grad_tol*(1 + |f_k|) */
    int    StopRule ;
    rp_ StopFac ;

    /* T => estimated error in function value is eps*Ck,
       F => estimated error in function value is eps */
    int    PertRule ;
    rp_ eps ;

    /* factor by which eps grows when line search fails during contraction */
    rp_ egrow ;

    /* T => attempt quadratic interpolation in line search when
                |f_k+1 - f_k|/f_k <= QuadCutoff
       F => no quadratic interpolation step */
    int    QuadStep ;
    rp_ QuadCutOff ;

    /* maximum factor by which a quad step can reduce the step size */
    rp_ QuadSafe ;

    /* T => when possible, use a cubic step in the line search */
    int UseCubic ;

    /* use cubic step when |f_k+1 - f_k|/|f_k| > CubicCutOff */
    rp_ CubicCutOff ;

    /* |f| < SmallCost*starting cost => skip QuadStep and set PertRule = FALSE*/
    rp_ SmallCost ;

    /* T => check that f_k+1 - f_k <= debugtol*C_k
       F => no checking of function values */
    int    debug ;
    rp_ debugtol ;

    /* if step is nonzero, it is the initial step of the initial line search */
    rp_ step ;

    /* abort cg after maxit iterations */
    INT maxit ;

    /* maximum number of times the bracketing interval grows during expansion */
    int ntries ;

    /* maximum factor secant step increases stepsize in expansion phase */
    rp_ ExpandSafe ;

    /* factor by which secant step is amplified during expansion phase
       where minimizer is bracketed */
    rp_ SecantAmp ;

    /* factor by which rho grows during expansion phase where minimizer is
       bracketed */
    rp_ RhoGrow ;

    /* maximum number of times that eps is updated */
    int neps ;

    /* maximum number of times the bracketing interval shrinks */
    int nshrink ;

   /* maximum number of iterations in line search */
    int nline ;

    /* conjugate gradient method restarts after (n*restart_fac) iterations */
    rp_ restart_fac ;

    /* stop when -alpha*dphi0 (estimated change in function value) <= feps*|f|*/
    rp_ feps ;

    /* after encountering nan, growth factor when searching for
       a bracketing interval */
    rp_ nan_rho ;

    /* after encountering nan, decay factor for stepsize */
    rp_ nan_decay ;

/*============================================================================
       technical parameters which the user probably should not touch
  ----------------------------------------------------------------------------*/
    rp_           delta ; /* Wolfe line search parameter */
    rp_           sigma ; /* Wolfe line search parameter */
    rp_           gamma ; /* decay factor for bracket interval width */
    rp_             rho ; /* growth factor when searching for initial
                                bracketing interval */
    rp_            psi0 ; /* factor used in starting guess for iteration 1 */
    rp_          psi_lo ; /* in performing a QuadStep, we evaluate at point
                                betweeen [psi_lo, psi_hi]*psi2*previous step */
    rp_          psi_hi ;
    rp_            psi1 ; /* for approximate quadratic, use gradient at
                                psi1*psi2*previous step for initial stepsize */
    rp_            psi2 ; /* when starting a new cg iteration, our initial
                                guess for the line search stepsize is
                                psi2*previous step */
    int       AdaptiveBeta ; /* T => choose beta adaptively, F => use theta */
    rp_       BetaLower ; /* lower bound factor for beta */
    rp_           theta ; /* parameter describing the cg_descent family */
    rp_            qeps ; /* parameter in cost error for quadratic restart
                                criterion */
    rp_           qrule ; /* parameter used to decide if cost is quadratic */
    int           qrestart ; /* number of iterations the function should be
                                nearly quadratic before a restart */
} cg_parameter ;

typedef struct cg_stats_struct /* statistics returned to user */
{
    rp_               f ; /*function value at solution */
    rp_           gnorm ; /* max abs component of gradient */
    INT               iter ; /* number of iterations */
    INT            IterSub ; /* number of subspace iterations */
    INT             NumSub ; /* total number subspaces */
    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) */
(
    rp_            *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 */
    rp_      grad_tol, /* StopRule = 1: |g|_infty <= max (grad_tol,
                                           StopFac*initial |g|_infty) [default]
                             StopRule = 0: |g|_infty <= grad_tol(1+|f|) */
    rp_        (*value) (rp_ *, INT),  /* f = value (x, n) */
    void           (*grad) (rp_ *, rp_ *, INT), /* grad (g, x, n) */
    rp_      (*valgrad) (rp_ *, rp_ *, INT), /* f = valgrad (g,x,n)*/
    rp_         *Work  /* either size 4n work array or NULL */
) ;

void cg_default /* set default parameter values */
(
    cg_parameter   *Parm
) ;
