.packageName <- "Rcgmin"
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/Rtmpu99x13/R.INSTALL10d63af1/Rcgmin/R/Rcgmin.R"
Rcgmin <- function( par, fn, gr=NULL, lower=NULL, upper=NULL, bdmsk=NULL, control=list(), ...) {
## An R version of the conjugate gradient minimization
## using the Dai-Yuan ideas
#  This version is for unconstrained functions.
#
# Input:
#  par  = a vector containing the starting point
#  fn = objective function (assumed to be sufficeintly differentiable)
#  gr = gradient of objective function
#  lower = vector of lower bounds on parameters
#  upper = vector of upper bounds on parameters
#    Note: free parameters outside bounds will be adjusted to bounds.
#  bdmsk = control vector for bounds and masks. Parameters for which bdmsk are 1
#         are unconstrained or "free", those with bdmsk 0 are masked i.e., fixed.
#         For historical reasons, we use the same array as an indicator that a
#         parameter is at a lower bound (-3) or upper bound (-1)
#  control = list of control parameters
#           maxit = a limit on the number of iterations (default 500)
#           maximize = TRUE to maximize the function (default FALSE)
#           trace = 0 (default) for no output, 
#                  >0 for output (bigger => more output)
#
# Output:
#    A list with components: 
#
#     par: The best set of parameters found.
#
#   value: The value of 'fn' corresponding to 'par'.
#
#  counts: A two-element integer vector giving the number of calls to
#          'fn' and 'gr' respectively. This excludes those calls needed
#          to compute the Hessian, if requested, and any calls to 'fn'
#          to compute a finite-difference approximation to the gradient.
#
# convergence: An integer code. '0' indicates successful convergence.
#          Error codes are
#          '0' converged 
#          '1' indicates that the function evaluation count 'maxfeval'
#               was reached.
#          '2' indicates initial point is infeasible
#
# message: A character string giving any additional information returned
#          by the optimizer, or 'NULL'.
#
#  Author:  John C Nash
#  Date:  April 2, 2009; revised July 28, 2009
#################################################################
  # control defaults -- idea from spg
  ctrl<-list(maxit=500, maximize=FALSE, trace=0, eps=1.0e-7, usenumDeriv=FALSE)
  namc <- names(control)
  if (! all(namc %in% names(ctrl)) )
     stop("unknown names in control: ", namc[!(namc %in% names(ctrl))])     
  ctrl[namc] <- control
  maxit <- ctrl$maxit  # limit on function evaluations
  maximize <- ctrl$maximize  # TRUE to maximize the function
  trace<- ctrl$trace         # 0 for no output, >0 for output (bigger => more output)
  if (trace>2) cat("trace = ",trace,"\n")
  eps<-ctrl$eps
  fargs <- list(...)         # the ... arguments that are extra function / gradient data
  grNULL <- is.null(gr)
#############################################
# local function defined only when user does not specify a gr
# Simple gr numerical approximation. Using fn, fmin and eps from calling env.  	
  if (grNULL) {
       warning("WARNING: Numerical gradients may be inappropriate for Rcgmin")
       if ( ctrl[["usenumDeriv"]] ) { # external gradient
         require(numDeriv)
         gr<-function(par, ...) {
           gg<-grad(myfn, par, ...)
           gg
         }
       } else { # using local gradient
        gr <- function(par, ...) {
        fbase<-myfn(par,...) # ensure we have right value, may not be necessary
    	df <- rep(NA,length(par))
        teps<-eps*(abs(par)+eps)
    	for (i in 1:length(par)) {
    	  dx <- par
    	  dx[i] <- dx[i] + teps[i]
    	  df[i] <- (myfn(dx, ...) - fbase)/teps[i]
    	 }
    	df
	}
       } # end else
} # grNULL
############# end gr ########################
  myfn <- if (maximize) function(par, ...) -fn(par, ...)
                   else function(par, ...)  fn(par, ...)

  mygr <- if (maximize && !grNULL) function(par, ...) -gr(par, ...)
                    else function(par, ...)  gr(par, ...)

## Set working parameters (See CNM Alg 22)
  if (trace>0) {
	cat("Rcgmin -- J C Nash 2009 - bounds constraint version of new CG\n")
	cat("an R implementation of Alg 22 with Yuan/Dai modification\n")
  }
  bvec <- par # copy the parameter vector
  n <- length(bvec) # number of elements in par vector
  maxfeval<-round(sqrt(n+1)*maxit) # change 091219
  ig <- 0 # count gradient evaluations
  ifn <- 1 # count function evaluations (we always make 1 try below)
  stepredn <- 0.15 # Step reduction in line search
  acctol <- 0.0001 # acceptable point tolerance
  reltest <- 100.0 # relative equality test
  ceps<-.Machine$double.eps*reltest
  setstep <- 1.75 # step increase factor
  accpoint <- as.logical(FALSE) # so far do not have an acceptable point
  fail <- as.logical(FALSE) # Method hasn't yet failed on us!
  cyclimit <- min(2.5*n,10+sqrt(n)) #!! upper bound on when we restart CG cycle 
#!! getting rid of limit makes it work on negstart BUT inefficient
  # This does not appear to be in Y H Dai & Y Yuan, Annals of Operations Research 103, 33–47, 2001 aor01.pdf
  # in Alg 22 pascal, we can set this as user. Do we wish to allow that?
  tol <- n*(n*.Machine$double.eps) # # for gradient test.  Note -- integer overflow if n*n*d.eps
  fargs<-list(...) # function arguments
  if (trace > 2) {
	cat("Extra function arguments:")
	print(fargs)
  }
# set default masks if not defined
  if (is.null(bdmsk)) {
       bdmsk<-rep(1,n) # all parameters free initially
  }
# check if there are bounds
  if(is.null(lower) || ! is.finite(lower)) nolower=TRUE else nolower=FALSE
  if(is.null(upper) || ! is.finite(upper)) noupper=TRUE else noupper=FALSE
  if(nolower && noupper && all(bdmsk == 1)) bounds=FALSE else bounds=TRUE
  if (trace > 2) cat("Bounds: nolower = ",nolower,"  noupper  ",noupper," bounds = ",bounds,"\n")
######## check bounds and masks #############
## NOTE: do this inline to avoid call (??should we change this?)
   if (bounds) {
# This implementation as a loop, but try later to vectorize
     for (i in 1:n) {
       if (bdmsk[i] == 0) { # NOTE: we do not change masked parameters, even if out of bounds
           if(bvec[i]<lower[i]) { 
              if (trace>0) cat("WARNING: ",bvec[i]," = MASKED par[",i,"] < lower bound = ",lower[i],"\n")
           }
           if(bvec[i]>upper[i]) { 
              if (trace>0) cat("WARNING: ",bvec[i]," = MASKED par[",i,"] > upper bound = ",upper[i],"\n")
           }
       } else { # not masked
           if(bvec[i] <= lower[i]) { 
              if (trace>0) cat("WARNING: par[",i,"], set ",bvec[i]," to lower bound = ",lower[i],"\n")
              bvec[i]<-lower[i]
              bdmsk[i] <- -3 # active lower bound
           }
           if(bvec[i] >= upper[i]) { 
              if (trace>0) cat("WARNING: par[",i,"], set ",bvec[i]," to upper bound = ",upper[i],"\n")
              bvec[i]<-upper[i]
              bdmsk[i] <- -1 # active upper bound
           }
       } # end not masked
      } # end loop for bound/mask check
      if (trace>0) {
	cat("adjusted bvec:")
	print (bvec)
      }
    } # end if bounds
############## end bounds check #############
# Initial function value -- may NOT be at initial point specified by user.
  cat("Try function at initial point:")
  print(bvec)
  f <- try(do.call("myfn", append(list(bvec), fargs )), silent=TRUE) # Compute the function at initial point.
  cat("f=",f,"\n")
  if ( class(f) == "try-error") { 
     msg<-"Initial point is infeasible."
     if(trace > 0) cat(msg,"\n")
     ans<-list(par, NA, c(ifn, 0), 2, msg)
     names(ans)<-c("par", "value", "counts", "convergence", "message")
     return(ans)
  } 
  fmin<-f
  if (trace>0) cat("Initial fn=",f,"\n")
  if (trace>2) print(bvec)
# Start the minimization process
  keepgoing<-TRUE
  notconv <-TRUE
  msg<-"not finished" # in case we exit somehow
  oldstep<-0.8 #!! 2/3 #!! 
####################################################################
  fdiff<-NA # initially no decrease
  cycle <- 0 # !! cycle loop counter
  while (keepgoing && notconv) {  # main loop -- must remember to break out of it!!
    t<-as.vector(rep(0,n)) # zero step vector
    c<-t # zero "last" gradient
    while (keepgoing && (cycle < cyclimit)) { ## cycle loop
      cycle<-cycle+1
      if (trace>0) cat(ifn," ",ig," ",cycle," ",fmin,"  last decrease=",fdiff,"\n")
      if( trace>2) {
        print(bvec)
        cat("\n")
      }
      if (ifn > maxfeval) {
        msg<-paste("Too many function evaluations (> ",maxfeval,") ", sep='')
        if (trace > 0 ) cat(msg,"\n")
        ans<-list(par, fmin, c(ifn, ig), 1, msg) # 1 indicates not converged in function limit
        names(ans)<-c("par", "value", "counts", "convergence", "message")
        return(ans)
      }
      par<-bvec # save best parameters
      ig<-ig+1
      if (ig > maxit) {
        msg<-paste("Too many gradient evaluations (> ",maxit,") ", sep='')
        if (trace > 0 ) cat(msg,"\n")
        ans<-list(par, fmin, c(ifn, ig), 1, msg) # 1 indicates not converged in function or gradient limit
        names(ans)<-c("par", "value", "counts", "convergence", "message")
        return(ans)
      }
      g<-mygr(bvec, ...)
      if(bounds) {
        ## Bounds and masks adjustment of gradient ##
        ## first try with looping -- later try to vectorize
        if (trace>2) {
           cat("bdmsk:")
           print(bdmsk)
        }
        for (i in 1:n) {
          if( (bdmsk[i]==0) ) { # masked, so gradient component is zero
            g[i]<-0
          } else {
            if (bdmsk[i]==1) {
               if (trace>1) cat("Parameter ",i," is free\n")
            } else {
               if ( (bdmsk[i]+2)*g[i]<0.0 ) { # test for -ve gradient at upper bound, +ve at lower bound
                 g[i] <- 0. # in which case active mask or constraint and zero gradient component
               } else {
                 bdmsk[i] <- 1 # freeing parameter i
                 if (trace>1) cat("freeing parameter ",i,"\n")
               }
            }
          }
        } # end masking loop on i
        if (trace>2) {
          cat("bdmsk adj:\n")
          print(bdmsk)
          cat("proj-g:\n")
          print(g)
        }
      } # end if bounds
      ## end bounds and masks adjustment of gradient
      g1<-sum(g*(g-c)) # gradient * grad-difference
      g2<-sum(t*(g-c)) # oldsearch * grad-difference
      gradsqr<-sum(g*g)
      if (trace > 1) {
          cat("Gradsqr = ", gradsqr," g1, g2 ", g1, " ", g2," fmin=",fmin,"\n")
          #!! tmp<-readline()
      }
      #!! gtc<-sum(g*c) #!! to check below -- but remove later
      #!! cat("t2=",sum(t*t),"\n")
      #!! if (gtc > 0.2*gradsqr) cat("Powell restart FLAGGED g*c=",gtc,"\n") #!!
      c<-g # save last gradient
      g3<-1.0 # !! Default to 1 to ensure it is defined -- t==0 on first cycle
      if (gradsqr > tol*(abs(fmin)+reltest) ) {
         if (g2 > 0.0) {
            betaDY<-gradsqr/g2
            betaHS<-g1/g2
            g3 <-max( 0.0, min(betaHS, betaDY) ) # g3 is our new 'beta' !! Dai/Yuan 2001, (4.2)
         } 
      } else {
	 msg<-paste("Very small gradient -- gradsqr =",gradsqr,sep=' ')
         if (trace > 0) cat(msg,"\n")
         keepgoing<-FALSE # done loops -- should we break ????
         break # to leave inner loop
      }    
      if(trace > 2) cat("Betak = g3 = ",g3,"\n")
      if (g3 == 0.0 || cycle>=cyclimit) {  # we are resetting to gradient in this case 
          if (trace > 0) {
             if(cycle<cyclimit) cat("Yuan/Dai cycle reset\n")
             else cat("Cycle limit reached -- reset\n")          
          }
          fdiff<-NA
          cycle<-0
          break #!!
          #!! oldstep<-1 # !!
          #!! don't reset stepsize ## oldstep<-1 #!! reset
          #!! break # to quit inner loop
      } else {
        # drop through if not Yuan/Dai cycle reset
        t <- t * g3 - g # t starts at zero, later is step vector
        gradproj <- sum(t*g) # gradient projection
        if(trace > 1) cat("Gradproj =",gradproj,"\n")
        if(bounds) {
          ## Adjust search direction for masks
          if (trace>2) {
            cat("t:\n")
            print(t)
          }
          t[which(bdmsk<=0)]<-0 # apply mask constraint
          if (trace>2) {
            cat("adj-t:\n")
            print(t)
          }
          ## end adjust search direction for masks
        } # end if bounds
	# ?? Why do we not check gradproj size??
        ########################################################
        ####                  Line search                   ####
        OKpoint<-FALSE
        if (trace>2) cat("Start linesearch with oldstep=",oldstep,"\n")
        steplength <- oldstep*1.5 #!! try a bit bigger
        f <- fmin
        changed<-TRUE # Need to set so loop will start
        while ((f >= fmin) && changed) {
          if (bounds) {
            # Box constraint -- adjust step length
            for (i in 1:n) { # loop on parameters -- vectorize??
               if( (bdmsk[i]==1) && (t[i] != 0.0) ) { 
                  # only concerned with free parameters and non-zero search dimension
                  if(t[i] < 0.0) { # going down. Look at lower bound
                     trystep<-(lower[i]-par[i])/t[i] # t[i] < 0 so this is positive
                  } else { # going up, check upper bound
                     trystep<-(upper[i]-par[i])/t[i] # t[i] > 0 so this is positive
                  }
                  if (trace>2) cat("steplength, trystep:", steplength, trystep,"\n")
                  steplength<-min(steplength, trystep) # reduce as necessary
               } # end steplength reduction
            } # end loop on i to reduce step length
            if (trace>1) cat("reset steplegth=",steplength,"\n")
            # end box constraint adjustment of step length
          } # end if bounds
          bvec<-par+steplength*t
          changed <- ( !  identical((bvec+reltest), (par+reltest)) )
          if( changed ) { # compute newstep, if possible
            f <- myfn(bvec, ...) # Because we need the value for linesearch, don't use try()
            #  instead preferring to fail out, which will hopefully be unlikely.
            ifn<-ifn+1
            if (is.na(f) || (! is.finite(f))) { 
		warning("Rcgmin - undefined function")
		f<-.Machine$double.xmax 
	    }
            if (f < fmin) {
              f1 <- f # Hold onto value 
            } else {
              steplength<-steplength*stepredn
              if (trace > 0) cat("*")
            }
          }
        }  # end while
        changed1<-changed # Change in parameters occured in step reduction
        if (changed1) {
          ## ?? should we check for reduction? or is this done in if (newstep >0) ?
          newstep <- 2.0*(f - fmin - gradproj * steplength) # JN 081219 change
          if (newstep > 0) {
              newstep = -(gradproj * steplength * steplength / newstep) 
          }
          if (bounds) {
            # Box constraint -- adjust step length
            for (i in 1:n) { # loop on parameters -- vectorize??
               if( (bdmsk[i]==1) && (t[i] != 0.0) ) { # only concerned with free parameters and non-zero search dimension
                  if(t[i] < 0.0) { # going down. Look at lower bound
                     trystep<-(lower[i]-par[i])/t[i] # t[i] < 0 so this is positive
                  } else { # going up, check upper bound
                     trystep<-(upper[i]-par[i])/t[i] # t[i] > 0 so this is positive
                  }
                  if (trace>2) cat("newstep, trystep:", newstep, trystep,"\n")
                  newstep<-min(newstep, trystep) # reduce as necessary
               } # end newstep reduction
            } # end loop on i to reduce step length
            if (trace > 2) cat("reset newstep=",newstep,"\n")
            # end box constraint adjustment of step length
          } # end if bounds
          bvec<-par+newstep*t
          changed <- ( ! identical((bvec+reltest), (par+reltest)) )
          if (changed) {
              f <- myfn(bvec, ...)
              ifn<-ifn+1
          }
          if(trace > 2) cat("fmin, f1, f: ",fmin, f1, f,"\n")
          if (f < min(fmin, f1)) {
             # success
             OKpoint<-TRUE
             accpoint <- ( f <= fmin + gradproj * newstep * acctol)
             fdiff<-(fmin - f) # check decrease
             fmin<-f
             oldstep<-newstep # !! save it
          } else { 
            if (f1 < fmin) {
                bvec<-par+steplength*t # reset best point
                accpoint <- ( f1 <= fmin + gradproj * steplength * acctol) 
                OKpoint<-TRUE # Because f1 < fmin
                fdiff<-(fmin - f1) # check decrease
                fmin<-f1
                oldstep<-steplength #!! save it
            } else {
                # no reduction
                fdiff<-NA
                accpoint<-FALSE
            } # f1<?fmin
          } # f < min(f1, fmin)
          if(trace > 1) cat("accpoint = ",accpoint," OKpoint = ",OKpoint,"\n")
          if (! accpoint) {
              msg<-"No acceptable point -- exit loop"
              if (trace > 0) cat("\n",msg,"\n")
              keepgoing<-FALSE
              break #!!
          }
        } # changed1
        else { # not changed on step redn
          if (cycle == 1) { 
            msg<-" Converged -- no progress on new CG cycle"
            if (trace > 0) cat("\n",msg,"\n")
            keekpgoing<-FALSE
            notconv<-FALSE # mark as converged
            break #!!
          } 
        } # end else
      } # end of test on Yuan/Dai condition
      #### End line search ####
      if (bounds) {
        ## Reactivate constraints?? -- should check for infinite bounds
        for (i in 1:n) {
          if (bdmsk[i]==1) { # only interested in free parameters
            if(is.finite(lower[i])) {# JN091020 -- need to use abs in case bounds negative
              if ( (bvec[i]-lower[i]) < ceps*(abs(lower[i])+1.0) ) { # are we near or lower than lower bd
                if (trace>2) cat("(re)activate lower bd ",i," at ",lower[i],"\n")
                bdmsk[i] <- -3
              } # end lower bd reactivate
            }
            if(is.finite(upper[i])){ # JN091020 -- need to use abs in case bounds negative
              if ( (upper[i]-bvec[i]) < ceps*(abs(upper[i])+1.0) ) { # are we near or above upper bd
                if (trace>2) cat("(re)activate upper bd ",i," at ",upper[i],"\n")
                bdmsk[i] <- -1
              } # end lower bd reactivate
            }
          } # end test on free params
        } # end reactivate constraints
      } # end if bounds
    } # end of inner loop (cycle)
    if (oldstep < acctol) { oldstep <- acctol } 
#!!    oldstep = setstep * newstep # change to newstep for steplength
    if (oldstep > 1.0) { oldstep <- 1.0 }
    if (trace > 1) cat("End inner loop, cycle =",cycle,"\n")
  } # end of outer loop
  msg<-"Rcgmin seems to have converged"
  if (trace>0) cat(msg,"\n")
 #  par: The best set of parameters found.
 #  value: The value of 'fn' corresponding to 'par'.
 #  counts: number of calls to 'fn' and 'gr' (2 elements)
 #  convergence: An integer code. '0' indicates successful convergence.
 #  message: A character string or 'NULL'.
  ans<-list(par, fmin, c(ifn, ig), 0, msg)
  names(ans)<-c("par", "value", "counts", "convergence", "message")
  return(ans)
} ## end of Rcgmin

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/Rtmpu99x13/R.INSTALL10d63af1/Rcgmin/R/optansout.R"
optansout <- function(ansdf, filename) {
##### OPEN ISSUES: (any date order)

##### IMPLEMENTED: (reverse date order)

#  A funtion to display and print to file (if present) the output of optimx
	if (! is.null(filename)) { sink(filename) } # may need to check paths
#	if (! exists(filename)) { sink(filename) } # may need to check paths
	tpar<-ansdf$par
	tdf<-ansdf
        ltvec<-length(tpar[[1]])
	for (i in 1:length(tpar)) {
	    tvec<-tpar[[i]]
	    ltvec<-length(tvec)
	    if (ltvec > 5) {
		tvec<-tvec[1:5]
	    }
	    tpar[[i]]<-tvec
	}
	tdf$par<-tpar
	pardf<-tdf[1]
	if (ltvec > 5) {
		names(pardf)[1]<-"first.5.par"
	} else {
		names(pardf)[1]<-"par"
	}
	print(pardf)	
	tdf$par<-NULL
	print(tdf)
	if (! is.null(filename)) { sink() } # to turn it off
#	if (! exists(filename)) { sink() } # to turn it off
	rm(tdf)
	rm(tpar)
	return(TRUE)
} ## end of optansout

