Nothing
cnvtst <- function (alpha, pnorm, xnorm,
dif, ftest, gnorm, gtp, f, flast, g, ipivot, accrcy) {
## Note: trace is integer (changed 180329)
## ---------------------------------------------------------
## test for convergence
## ---------------------------------------------------------
## set up
## ---------------------------------------------------------
conv <- 0;
eps <- .Machine$double.eps
toleps <- sqrt(accrcy) + sqrt(eps);
rtleps <- accrcy + eps;
imax <- 0;
ltest <- (flast - f <= -0.5*gtp);
## ---------------------------------------------------------
## if anti-zigzag test satisfied, test multipliers;
## if appropriate, modify active set
## ---------------------------------------------------------
if ( ! ltest) {
ind <- which( (ipivot != 0) & (ipivot != 2))
if ( length(ind) > 0 ) { ## how to ensure ind not empty
t <- -sum(ipivot[ind]*g[ind])
cmax <- min(t) # ?? why [cmax, imax] = min(t) ??
imax <- cmax
if (cmax >= 0) { imax <- 0 }
}
}
if (imax != 0) {
ipivot[ ind[imax] ] <- 0;
flast <- f;
} else {
conv = ( ( (alpha*pnorm < toleps*(1 + xnorm) )
&& (abs(dif) < rtleps*ftest)
&& (gnorm < accrcy^(1/3)*ftest) ) ||
( gnorm < .01*sqrt(accrcy)*ftest) )
}
result <- list(conv=conv, flast1=flast, ipivot1=ipivot)
}
crash <- function (x, low, up) {
##---------------------------------------------------------
## this initializes the constraint information, and
## ensures that the initial point satisfies
## low <= x <= up.
## the constraints are checked for consistency.
##---------------------------------------------------------
ierror <- 0
if (any(low > up)) { ierror = - max(which(low > up)) }
# above is check on error in bounds specification
xnew <- pmax (low, x) # force params into bounds
xnew <- pmin (up, xnew) # No diagnostic! ??
# we output revised parameters
n <- length(x)
ind <- which(low == x)
ipivot <- rep(0, n) ## zeros(size(x)) in Matlab
if (length(ind) > 0) { ipivot[ind] <- -1 }
ind <- which(x == up)
if (length(ind) > 0) { ipivot[ind] <- 1 }
ind <- which(low == up)
if (length(ind) > 0) { ipivot[ind] <- 2 }
list(ipivot=ipivot, ierror=ierror, xnew=xnew)
}
gtims <- function (v, x, g, accrcy, xnorm, sfun, ...) {
##---------------------------------------------------------
## compute the product of the Hessian times the vector v;
## store result in the vector gv
## (finite-difference version)
##---------------------------------------------------------
## cat("v, x, g:")
## print(v)
## print(x)
## print(g)
delta <- sqrt(accrcy)*(1 + xnorm)/sqrt(sum(v^2))
hg <- x + delta*v
tresult<-sfun(hg, ...)
## cat("tresult: ")
## print(tresult)
gv <- attr(tresult, "gradient")
gv <- (gv - g)/delta
gv
}
initpc <- function(d, upd1, ireset) {
##---------------------------------------------------------
## initialize the diagonal preconditioner (d -- a vector!)
## ---------------------------------------------------------
## global hyk sk yk sr yr yksk yrsr ?? do we have these
# global vectors hyk sk yk sr yr & scalars yksk yrsr
## ---------------------------------------------------------
#% cat("In initpc -- sk and hyk if upd1 false: ", upd1,"\n")
#% print(envjn$sk)
if (upd1) {
td <- d
} else {
if (ireset) {
envjn$hyk <- d * envjn$sk # vector * vector by element??
sds <- as.numeric(crossprod(envjn$sk, envjn$hyk)) # matrix multiply
if (all(envjn$hyk == 0) && trace > 1) { cat("INITPC: envjn$hyk = 0 \n") }
if (sds == 0 && trace > 1) { cat("INITPC: sds = 0 \n") }
td <- d - d*d*envjn$sk*envjn$sk/sds + envjn$yk*envjn$yk/envjn$yksk
# by element
} else {
envjn$hyk <- d * envjn$sr # vector * vector
sds <- as.numeric(crossprod(envjn$sr, envjn$hyk))
srds <- as.numeric(crossprod(envjn$sk, envjn$hyk))
yrsk <- as.numeric(crossprod(envjn$yr, envjn$sk))
envjn$hyk <- d*envjn$sk - envjn$hyk*srds/sds + envjn$yr*yrsk/envjn$yrsr
td <- d - d*d*envjn$sr*envjn$sr/sds+envjn$yr*envjn$yr/envjn$yrsr
sds <- as.numeric(crossprod(envjn$sk, envjn$hyk))
td <- td - envjn$hyk*envjn$hyk/sds + envjn$yk*envjn$yk/envjn$yksk
}
}
#% cat("td:")
#% print(td)
ans<-list(td=td)
}
lin1 <- function(p, x, f, alpha, g, sfun, ...){
## ---------------------------------------------------------
## line search (naive)
## ---------------------------------------------------------
## set up
## ---------------------------------------------------------
# cat("lin1: alpha=",alpha," p:\n")
# print(p)
if (is.null(alpha)) alpha <- 0
ierror <- 3
xnew <- x
fnew <- f
gnew <- g
maxit <- 15
if (alpha == 0) {
ierror <- 0
maxit <- 1
}
alpha1 <- alpha
## ---------------------------------------------------------
## line search
## ---------------------------------------------------------
for (itcnt in 1:maxit) {
xt <- x + alpha1*p
fg <- sfun(xt, ...) # Note: added dots 140902
ft<-fg
gt<- attr(fg,"gradient") # may simplify later
if (ft < f) {
ierror <- 0
xnew <- xt
fnew <- ft
gnew <- gt
## cat("about to break in lin1\n")
break
}
alpha1 <- alpha1 / 2
}
if (ierror == 3) { alpha1 <- 0 }
nf1 <- itcnt
## never used in SGN code
## if (nargout == 7) { ## ?? what is nargout?
## dfdp <- as.numeric(crossprod(gt, p) )
## varargout{1} <- dfdp ## ?? what is varargout
## end
result<-list(xnew=xnew, fnew=fnew, gnew=gnew, nf1=nf1,
ierror=ierror, alpha1=alpha1)
}
lmqnbc <- function (x, sfun, lower, upper, maxit, maxfun, stepmx, accrcy, trace, ...) {
## ---------------------------------------------------------
## This is a bounds-constrained truncated-newton method.
## The truncated-newton method is preconditioned by a
## limited-memory quasi-newton method (computed by
## this routine) with a diagonal scaling (routine ndia3).
## For further details, see routine tnbc.
## ---------------------------------------------------------
## global hyk sk yk sr yr yksk yrsr
## ---------------------------------------------------------
## check that initial x is feasible and that the bounds
## are consistent
## ---------------------------------------------------------
n <- length(x)
# JN: Define globals here
gtn<-list(yrsr=0, yksk=0, yr = rep(0, n), yk = rep(0, n),
sr = rep(0, n), sk = rep(0, n),
hg=rep(0,n), hyk=rep(0,n), hyr=rep(0,n) )
envjn<<-list2env(gtn) # Note this important tool for globals in Rtnmin
# end globals
## [ipivot, ierror, x] = crash(x, lower, upper)
if (trace > 1) cat("lmqnbc -- crout:")
crout<-crash(x, lower, upper)
if (trace) print(crout)
ierror <- crout$ierror
ipivot <- crout$ipivot
x<-crout$xnew # in case x changed by bounds
f <- 0
g <- rep(0,n)
if (ierror != 0) {
stop('LMQNBC: terminating (no feasible point)')
## return
}
## ---------------------------------------------------------
## initialize variables, parameters, and constants
## ---------------------------------------------------------
options(digits=5)
if (trace > 1) cat("stempmx, accrcy, maxfun:",stepmx, accrcy, maxfun, "\n")
if (trace > 0) cat(' it nf cg f |g|\n')
eps <- .Machine$double.eps
upd1 <- TRUE
ncg <- 0
conv <- FALSE
xnorm <- max(abs(x))
ierror <- 0
if ( (stepmx < sqrt(accrcy)) || (maxfun < 1) ) {
ierror <- -1
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror,
nfngr=ncg)
if (trace > 1) {
cat("Exiting lmqnbc - almgn:")
print(almqn)
}
return(almqn) # return
}
## ---------------------------------------------------------
## compute initial function value and related information
## ---------------------------------------------------------
if (trace > 1) cat("Try initial fn\n")
fg <- sfun(x, ...)
nf <- 1
nit <- 0
g<- attr(fg,"gradient")
f<-fg
flast <- f
# if (is.null(g) ) { ## 160922 change
if (! is.numeric(g) ) { # try fix 180328
gnorm <- 1.0/eps
} else { gnorm <- max(abs(g)) } ## norm(g,'inf')
## ---------------------------------------------------------
## Test if Lagrange multipliers are non-negative.
## Because the constraints are only bounds, the Lagrange
## multipliers are components of the gradient.
## Then form the projected gradient.
## ---------------------------------------------------------
ind <- which((ipivot != 2) &
(as.numeric(crossprod(ipivot,g)) >0 ) )
if (length(ind) > 0) {
ipivot[ind] <- rep(0, length(ind))
}
g <- ztime (g, ipivot)
# if (is.null(g) ) { ## 160922 change
if (! is.numeric(g) ) { # try fix 180328
gnorm <- 1.0/eps
} else { gnorm <- max(abs(g)) } ## norm(g,'inf')
if (trace > 0) cat(nit,"\t", nf,"\t", ncg,"\t", f," ", gnorm,"\n")
if (trace > 1) {
print(x)
print(g)
}
## ---------------------------------------------------------
## check if the initial point is a local minimum.
## ---------------------------------------------------------
ftest <- 1 + abs(f)
if (gnorm < .01*sqrt(eps)*ftest) {
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror,
nfngr=ncg)
return(almqn)
}
if (trace > 1) cat("before set initial, flast=",flast,"\n")
## ---------------------------------------------------------
## set initial values to other parameters
## ---------------------------------------------------------
icycle <- n-1
ireset <- 0
bounds <- TRUE
difnew <- 0
epsred <- .05
fkeep <- f
d <- rep(1,n)
## ---------------------------------------------------------
## ..........main iterative loop..........
## ---------------------------------------------------------
## compute the search direction
## ---------------------------------------------------------
argvec <- c(accrcy, gnorm, xnorm)
##[p, gtp, ncg1, d] <- ...
## modlnp (d, x, g, maxit, upd1, ireset, bounds, ipivot, ## argvec, sfun)
mres <- modlnp (d, x, g, maxit, upd1, ireset,
bounds, ipivot, argvec, sfun, ...)
ncg1 <- mres$ncg1
gtp <- mres$gtp
d <- mres$dnew
p <- mres$p
## cat("p:")
## print(p)
## tmp<-readline("cont.")
ncg <- ncg + ncg1
while (!conv) {
oldg <- g
pnorm <- max(abs(p)) # norm2(p, 'inf')
oldf <- f
## ---------------------------------------------------------
## line search
## ---------------------------------------------------------
pe <- pnorm + eps
spe <- stpmax (stepmx, pe, x, p, ipivot, lower, upper)
## cat("spe=",spe,"\n")
alpha <- step1 (f, gtp, spe)
alpha0 <- alpha
## cat("alpha0=",alpha0,"\n")
## [x_new, f_new, g_new, nf1, ierror, alpha] <- lin1 (p, x,
## f, alpha0, g, sfun)
reslin <- lin1 (p, x, f, alpha0, g, sfun, ...)
ierror <- reslin$ierror
alpha <- reslin$alpha1
## ---------------------------------------------------------
if ((alpha == 0) && (alpha0 != 0) || (ierror == 3)){
cat('Error in Line Search\n')
cat(' ierror = ', ierror, "\n")
cat(' alpha = ',alpha, "\n")
cat(' alpha0 = ', alpha0, "\n")
cat(' gtp = ', gtp, "\n")
## ############################
cat(' |g| = ', norm2(g), "\n")
cat(' |p| = ', norm2(p), "\n")
tmp <- readline('Hit any key to continue')
## ############################
}
## #######################
x <- reslin$xnew # need fixup
f <- reslin$fnew
g <- reslin$gnew
nf1 <- reslin$nf1
nf <- nf + nf1
nit <- nit + 1
## ---------------------------------------------------------
## update active set, if appropriate
## ---------------------------------------------------------
newcon <- FALSE
## cat("Check active set - alpha, spe, f:",alpha, spe,f,"\n")
if (abs(alpha-spe) <= 10*eps) {
newcon <- TRUE
ierror <- 0
## cat("flast:",flast,"\n")
# [ipivot, flast] <- modz (x, p, ipivot, lower, upper,
# flast, f, alpha)
## cat("ipivot before and after update by modz:\n")
## print(ipivot)
modzres<-modz(x, p, ipivot, lower, upper,
flast, f, alpha)
ipivot <- modzres$ipivot1
flast<-modzres$flast1
## print(ipivot)
}
if (ierror == 3) {
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror,
nfngr=ncg)
if (trace > 0) {
cat("Error updating active set -- almqn:")
print(almqn)
}
return(almqn)
}
## ---------------------------------------------------------
## stop if more than maxfun evaluations have been made
## ---------------------------------------------------------
if (nf > maxfun) {
ierror <- 2
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror,
nfngr=ncg)
if (trace > 0) {
cat("Too many function evaluations -- almqn:")
print(almqn)
}
return(almqn)
}
## ---------------------------------------------------------
## set up for convergence and resetting tests
## ---------------------------------------------------------
difold <- difnew
difnew <- oldf - f # scalar
if (icycle == 1) {
if (difnew > 2*difold) { epsred <- 2*epsred }
if (difnew < .5*difold) { epsred <- .5*epsred }
}
gv <- ztime (g, ipivot)
gnorm <- max(abs(gv))
ftest <- 1 + abs(f)
xnorm <- max(abs(x))
############### DISPLAY ##############
if (trace > 0) cat(nit,"\t", nf,"\t", ncg,"\t", f," ", gnorm,"\n")
## print(x)
## print(g)
## ---------------------------------------------------------
## test for convergence
## ---------------------------------------------------------
## [conv, flast, ipivot] <- cnvtst (alpha, pnorm, xnorm, ...
## difnew, ftest, gnorm, gtp, f, flast, g, ...
## ipivot, accrcy)
## cat("before cnvtst, flast=",flast,"\n")
ctres <- cnvtst (alpha, pnorm, xnorm,
difnew, ftest, gnorm, gtp, f, flast, g,
ipivot, accrcy)
conv <- ctres$conv
flast <- ctres$flast1
ipivot <- ctres$ipivot1
## cat("after cnvtst - conv, flast, ipivot:", conv, flast,"\n")
## print(ipivot)
if (conv) {
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror,
nfngr=ncg)
return(almqn)
}
g <- ztime (g, ipivot)
## ---------------------------------------------------------
## modify data for LMQN preconditioner
## ---------------------------------------------------------
if (! newcon) { ## 0 is FALSE and 1 is TRUE supposedly
envjn$yk <- g - oldg
envjn$sk <- alpha*p
envjn$yksk <- as.numeric(crossprod(envjn$yk, envjn$sk))
## cat("yksk:",envjn$yksk,"\n")
ireset <- ( (icycle == n-1) |
(difnew < epsred*(fkeep-f)) )
if (! ireset) {
envjn$yrsr <-
as.numeric(crossprod(envjn$yr,envjn$sr))
ireset <- (envjn$yrsr <= 0)
}
upd1 <- (envjn$yksk <= 0)
}
## cat("newcon, upd1:", newcon, upd1,"\n")
## ---------------------------------------------------------
## compute the search direction
## ---------------------------------------------------------
argvec <- c(accrcy, gnorm, xnorm)
## [p, gtp, ncg1, d] <- ...
## modlnp (d, x, g, maxit, upd1, ireset, bounds,
## ipivot, argvec, sfun)
mres <- modlnp (d, x, g, maxit, upd1, ireset,
bounds, ipivot, argvec, sfun, ...)
ncg1 <- mres$ncg1
gtp <- mres$gtp
d <- mres$dnew
p <- mres$p
## cat("New p:")
## print(p)
## tmp<-readline("cont.")
ncg <- ncg + ncg1
## ---------------------------------------------------------
## update LMQN preconditioner
## ---------------------------------------------------------
if (! newcon) {
if (ireset) {
envjn$sr <- envjn$sk
envjn$yr <- envjn$yk
fkeep <- f
icycle <- 1
} else {
envjn$sr <- envjn$sr + envjn$sk
envjn$yr <- envjn$yr + envjn$yk
icycle <- icycle + 1
}
}
} # end while !conv
}
lmqn <- function (x, sfun, maxit, maxfun, stepmx, accrcy, trace, ...) {
## ---------------------------------------------------------
## truncated-newton method for unconstrained minimization
## (customized version)
## ---------------------------------------------------------
# global vectors hyk sk yk sr yr & scalars yksk yrsr
## ---------------------------------------------------------
## set up
## ---------------------------------------------------------
## format compact
## format short e
n<-length(x)
# JN: Define globals here. Is it necessary to set up.
gtn<-list(yrsr=0, yksk=0, yr = rep(0, n), yk = rep(0, n), sr = rep(0, n), sk = rep(0, n))
envjn<<-list2env(gtn)
# end globals
eps <- .Machine$double.eps
upd1 <- 1
ncg <- 0
xnorm <- max(abs(x)) # norm(x,'inf')
ierror <- 0
if (stepmx < sqrt(accrcy) || maxfun < 1) {
ierror <- -1
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror,
nfngr=ncg)
return(almqn)
}
## ---------------------------------------------------------
## compute initial function value and related information
## ---------------------------------------------------------
fg <- sfun(x, ...)
#% print(fg)
g<- attr(fg, "gradient")
#% print(g)
if (is.null(g)) stop("Must have gradient defined for Rtnmin")
f<-fg
# if (is.null(g) ) { ## 160922 change
if (! is.numeric(g) ) { # try fix 180328
gnorm <- 1.0/eps
} else { gnorm <- max(abs(g)) } ## norm(g,'inf')
nf <- 1
nit <- 0
if (trace) cat("Itn ",nit," ",nf," ",ncg, " ",f, " ", gnorm,"\n")
## ---------------------------------------------------------
## check for small gradient at the starting point.
## ---------------------------------------------------------
ftest <- 1 + abs(f)
if (gnorm < .01*sqrt(eps)*ftest) {
ierror <- 0
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror, nfngr=ncg)
return(almqn)
}
## ---------------------------------------------------------
## set initial values to other parameters
## ---------------------------------------------------------
n <- length(x)
icycle <- n-1
toleps <- sqrt(accrcy) + sqrt(eps)
rtleps <- accrcy + eps
difnew <- 0
epsred <- .05
fkeep <- f
conv <- FALSE
ireset <- 0
ipivot <- 0
if (trace > 1) cat("end initial values\n")
## ---------------------------------------------------------
## initialize diagonal preconditioner to the identity
## ---------------------------------------------------------
d <- rep(1,n) # as a vector
## ---------------------------------------------------------
## ..........main iterative loop..........
## ---------------------------------------------------------
## compute search direction
## ---------------------------------------------------------
argvec <- c(accrcy, gnorm, xnorm)
# cat("call modlnp\n")
mres <- modlnp (d, x, g, maxit, upd1, ireset, bounds=FALSE, ipivot, argvec, sfun, ...)
p <- mres$p
# cat("p from first call to modlnp\n")
# print(p)
# tmp<-readline("cont.")
gtp <- mres$gtp
ncg1<-mres$ncg1
d <- mres$dnew
#% cat("d from modlnp:")
#% print(d)
ncg <- ncg + ncg1
while ( ! conv) {
#% cat("top while ")
#% tmp <- readline("Top of iteration")
oldg <- g
pnorm <- max(abs(p)) # norm(p,'inf')
oldf <- f
## ---------------------------------------------------------
## line search
## ---------------------------------------------------------
if (trace > 1) cat("start linesearch\n")
pe <- pnorm + eps
spe <- stepmx/pe
if (trace > 1) cat("gtp, spe:", gtp, spe,"\n")
alpha0 <- step1 (f, gtp, spe)
reslin <- lin1 (p, x, f, alpha0, g, sfun, ...)
## [x, f, g, nf1, ierror, alpha] <-
x <- reslin$xnew # need fixup
f <- reslin$fnew
g <- reslin$gnew
nf1 <- reslin$nf1
ierror <- reslin$ierror
alpha <- reslin$alpha1
# cat("after lin1, alpha=",alpha,"\n")
# tmp<-readline("cont.")
nf <- nf + nf1
## ---------------------------------------------------------
nit <- nit + 1
# if (is.null(g) ) { ## 160922 change
if (! is.numeric(g) ) { # try fix 180328
gnorm <- 1.0/eps
} else { gnorm <- max(abs(g)) } ## norm(g,'inf')
### Display info
if (trace > 0) cat("Itn ",nit," ",nf," ",ncg, " ",f, " ", gnorm,"\n")
if (ierror == 3) {
if (length(ncg) == 0) { ncg <- 0 } # ?? is.null(ncg)??
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror, nfngr=ncg)
return(almqn)
}
## ---------------------------------------------------------
## stop if more than maxfun evalutations have been made
## ---------------------------------------------------------
if (nf >= maxfun) {
ierror <- 2
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror, nfngr=ncg)
return(almqn)
}
## ---------------------------------------------------------
## set up for convergence and resetting tests
## ---------------------------------------------------------
ftest <- 1 + abs(f)
xnorm <- max(abs(x)) # norm(x,'inf')
difold <- difnew
difnew <- oldf - f
envjn$yk <- g - oldg
envjn$sk <- alpha*p
if (icycle == 1) {
if (difnew > 2*difold) { epsred <- 2*epsred }
if (difnew < 0.5*difold) { epsred <- 0.5*epsred }
}
## ---------------------------------------------------------
## convergence test
## ---------------------------------------------------------
conv <- ( ( (alpha*pnorm < toleps*(1 + xnorm)) &&
(abs(difnew) < rtleps*ftest) &&
(gnorm < accrcy^(1/3)*ftest) )||
( gnorm < .01*sqrt(accrcy)*ftest ) )
if (conv) {
ierror <- 0
xstar <- x
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror, nfngr=ncg)
return(almqn)
}
## ---------------------------------------------------------
## update lmqn preconditioner
## ---------------------------------------------------------
envjn$yksk <- as.numeric(crossprod(envjn$yk, envjn$sk) )
ireset <- ((icycle == n-1) || (difnew < epsred*(fkeep-f)) )
if ( ! ireset) {
envjn$yrsr <- as.numeric(crossprod(envjn$yr, envjn$sr))
ireset <- (envjn$yrsr <= 0)
}
upd1 <- (envjn$yksk <= 0)
## ---------------------------------------------------------
## compute search direction
## ---------------------------------------------------------
argvec <- c(accrcy, gnorm, xnorm)
## cat("ireset =", ireset," upd1=",upd1,"\n")
#% cat("New d to modlnp:")
#% print(d)
## [p, gtp, ncg1, d] <- ...
mres <- modlnp (d, x, g, maxit, upd1, ireset, 0, ipivot, argvec, sfun, ...)
p <- mres$p
gtp <- mres$gtp
ncg1<-mres$ncg1
d <- mres$dnew
#% cat("returned dnew, envjn$john: ", envjn$john)
#% print(d)
ncg <- ncg + ncg1
## ---------------------------------------------------------
## store information for lmqn preconditioner
## ---------------------------------------------------------
if (ireset) {
envjn$sr <- envjn$sk
envjn$yr <- envjn$yk
fkeep <- f
icycle <- 1
} else {
envjn$sr <- envjn$sr + envjn$sk
envjn$yr <- envjn$yr + envjn$yk
icycle <- icycle + 1
}
} # end while
## [xstar, f, g, ierror] = ..
almqn<-list(xstar=xstar, f=f, g=g, ierror=ierror, nfngr=ncg)
} # end lmqn
modlnp <- function(d, x, g, maxit, upd1, ireset, bounds,
ipivot, argvec, sfun, ...) {
##---------------------------------------------------------
## this routine performs a preconditioned conjugate-gradient
## iteration to solve the Newton equations for a search
## direction for a truncated-newton algorithm.
## When the value of the quadratic model is sufficiently
## reduced, the iteration is terminated.
##---------------------------------------------------------
## parameters
##
## p - computed search direction
## g - current gradient
## gv,gz1,v - scratch vectors
## r - residual
## d - diagonal preconditoning matrix
## feval - value of quadratic function
##------------------------------------------------------------
## initialization
##------------------------------------------------------------
if (is.null(d)) stop("Null d")
## print(x)
## print(g)
## cat(bounds,"\n")
accrcy <- argvec[1]
gnorm <- argvec[2]
xnorm <- argvec[3]
if (maxit == 0) {
p <- -g
gtp <- as.numeric(crossprod(p, g))
ncg1 <- 1
dnew <- d
if (sqrt(sum(p^2))==0) {
## cat("MODLNP 01: |p| = 0\n")
## pause(1)
}
#% cat("modlnp - dout01:")
#% print(dnew)
result<-list(p=p, gtp=gtp, ncg1=ncg1, dnew=dnew)
return(result)
}
first <- 1
tol <- 1e-6 ######### was 1.d-12 #########
qold <- 0
ainit <- initpc (d, upd1, ireset)
dnew<-ainit$td
#% cat("after initpc dnew:")
#% print(dnew)
r <- -g
v <- rep(0, length(r))
p <- v
gtp <- as.numeric(crossprod(p, g))
rho <- rep(0, maxit+1)
beta <- rep(0, maxit)
v.gv <- rep(0, maxit)
rho[[1]] <- as.numeric(crossprod(r))
##------------------------------------------------------------
## main iteration (conjugate-gradient iterations for Ax = b)
##------------------------------------------------------------
ind <- 0
ncg1 <- 0
for (k in 1:maxit) {
ncg1 <- ncg1 + 1
if (bounds) { r <- ztime(r, ipivot) }
amsolve <- msolve (r, upd1, ireset, first, d)
zk<-amsolve$y
if (bounds) { zk <- ztime (zk, ipivot) }
#% cat("r:")
#% print(r)
#% cat("zk:")
#% print(zk)
rz <- as.numeric(crossprod(r,zk))
if (rz/gnorm < tol) {
ind <- 80
if (sqrt(sum(p^2))==0) {
p <- -g
gtp <- as.numeric(crossprod(p, g))
}
#% cat("modlnp - dout - ind80:")
#% print(dnew)
result<-list(p=p, gtp=gtp, ncg1=ncg1, dnew=dnew)
return(result)
}
if (k > 1) {
beta[k] <- rz/rzold
} else {
beta[k] <- 0
}
#% cat("beta[",k,"] =",beta[k],"\n")
v <- zk + beta[k]*v
if (bounds) { v <- ztime( v, ipivot) }
## cat("about to call gtims\n")
gv <- gtims(v, x, g, accrcy, xnorm, sfun, ...)
#% cat("After gtims: gv=")
#% print(gv)
## cat(bounds,"\n")
if (bounds) { gv <- ztime (gv, ipivot) }
#% cat("gv, v:")
## cat("gv=", gv, "\n")
#% print(gv)
#% print(v)
v.gv[[k]] <- as.numeric(crossprod(v, gv)) ## ?? v.gv has underscore!! ??
#% cat("v.gv[[",k,"]]=",v.gv[[k]],"\n")
if (v.gv[[k]]/gnorm < tol) {
ind <- 50
if (sqrt(sum(p^2))==0) {
p <- -g ## Mod SGN 140912
if (bounds) {
p <- ztime(p, ipivot)
}
gtp <- crossprod(p, g)
## cat("MODLNP 03: |p| = 0 \n")
## pause(1)
}
#% cat("modlnp - dout03:")
#% print(dnew)
result<-list(p=p, gtp=gtp, ncg1=ncg1, dnew=dnew)
return(result)
}
dnew <- ndia3(dnew, v, gv, r, v.gv[[k]]) # NB need to return something but it doesn't get used
#% cat("after ndia3 below dout03, dnew:")
#% print(dnew)
##------------------------------------------------------------
## compute current solution and related vectors
##------------------------------------------------------------
alpha <- rz / v.gv[[k]]
p <- p + alpha* v
r <- r - alpha*gv
rho[[k+1]] <- as.numeric(crossprod(r))
##------------------------------------------------------------
## test for convergence
##------------------------------------------------------------
gtp <- as.numeric(crossprod(p,g))
pr <- as.numeric(crossprod(r,p))
q <- (gtp + pr) / 2
qtest <- k * (1 - qold/q)
if (qtest <= 0.5) {
if (sqrt(sum(p^2))==0) {
## cat("MODLNP 04: |p| = 0\n")
## pause(1)
}
#% cat("modlnp - dout04:")
#% print(dnew)
result<-list(p=p, gtp=gtp, ncg1=ncg1, dnew=dnew)
return(result)
}
qold <- q
##------------------------------------------------------------
## perform cautionary test
##------------------------------------------------------------
if (gtp > 0) {
ind <- 40
if (sqrt(sum(p^2))==0) {
## cat("MODLNP 05: |p| = 0 \n")
## pause(1)
}
#% cat("modlnp - dout05:")
#% print(dnew)
result<-list(p=p, gtp=gtp, ncg1=ncg1, dnew=dnew)
return(result)
}
rzold <- rz
} ## end loop over k
k <- k-1
##------------------------------------------------------------
## terminate algorithm
##------------------------------------------------------------
if (ind == 40) {
p <- p - alpha*v
}
##------------------------------------------------------------
if (ind == 50 && k <= 1) {
amsolve <- msolve (g, upd1, ireset, first, d)
p<-amsolve$y
p <- -p
if (bounds) { p <- ztime (p, ipivot) }
}
##------------------------------------------------------------
if (ind == 80 && k <= 1) {
p <- -g
if (bounds) { p <- ztime (p, ipivot) }
}
##------------------------------------------------------------
## store new diagonal preconditioner
##------------------------------------------------------------
gtp <- as.numeric(crossprod(p, g))
ncg1 <- k + 1
if (sqrt(sum(p^2))==0) {
## cat("MODLNP 06: |p| = 0 \n")
## pause(1)
}
## cat("modlnp - dout06:")
#% print(dnew)
result<-list(p=p, gtp=gtp, ncg1=ncg1, dnew=dnew)
return(result)
}
modz <- function (x, p, ipivot, low, up, flast, f, alpha){
##---------------------------------------------------------------------
## update the constraint matrix if a new constraint is encountered
##---------------------------------------------------------------------
eps <- .Machine$double.eps
indl <- which(ipivot == 0 & p < 0)
if (length(indl) > 0) {
toll <- 10 * eps * (abs(low[indl]) + 1)
hitl <- which(x[indl]-low[indl] <= toll)
if (length(hitl) > 0) {
flast <- f
ipivot[indl[hitl]] <- -1
}
}
##---------------------------------------------------------------------
indu <- which((ipivot == 0) & (p > 0));
if (length(indu) > 0) {
tolu <- 10 * eps * (abs( up[indu]) + 1)
hitu <- which(up[indu]-x[indu] <= tolu)
if (length(hitu) > 0) {
flast <- f
ipivot[indu[hitu]] <- 1
}
}
##---------------------------------------------------------------------
flast1 <- flast
ipivot1 <- ipivot
list(flast1 = flast, ipivot1=ipivot)
}
msolve <- function(g, upd1, ireset, first, d) {
##---------------------------------------------------------
## This routine acts as a preconditioning step for the
## linear conjugate-gradient routine. It is also the
## method of computing the search direction from the
## gradient for the non-linear conjugate-gradient code.
## It represents a two-step self-scaled bfgs formula.
##---------------------------------------------------------
#% cat("msolve, envjn$john=",envjn$john,"\n")
envjn$john<-"msolve"
if (upd1) {
#% cat("upd1 is TRUE\n")
y <- g/d
} else {
gsk <- as.numeric(crossprod(g, envjn$sk))
if (ireset) {
envjn$hg <- g/d
if (first) {
envjn$hyk <- envjn$yk/d
ykhyk <- as.numeric(crossprod(envjn$yk, envjn$hyk))
}
ghyk <- as.numeric(crossprod(g, envjn$hyk))
y <- ssbfgs(envjn$sk,envjn$hg,envjn$hyk,envjn$yksk,ykhyk,gsk,ghyk)
} else {
envjn$hg <- g/d
if (first) {
envjn$hyk <- envjn$yk/d
envjn$hyr <- envjn$yr/d
envjn$yksr <- as.numeric(crossprod(envjn$yk, envjn$sr))
ykhyr <- as.numeric(crossprod(envjn$yk, envjn$hyr))
}
gsr <- as.numeric(crossprod(g, envjn$sr))
ghyr <- as.numeric(crossprod(g, envjn$hyr))
if (first) {
yrhyr <- as.numeric(crossprod(envjn$yr, envjn$hyr))
}
envjn$hg <- ssbfgs(envjn$sr,envjn$hg,envjn$hyr,envjn$yrsr,yrhyr,gsr,ghyr)
if (first) {
envjn$hyk <- ssbfgs(envjn$sr,envjn$hyk,envjn$hyr,envjn$yrsr,yrhyr,envjn$yksr,ykhyr)
}
ykhyk <- as.numeric(crossprod(envjn$hyk, envjn$yk))
ghyk <- as.numeric(crossprod(envjn$hyk, g))
y <- ssbfgs(envjn$sk,envjn$hg,envjn$hyk,envjn$yksk,ykhyk,gsk,ghyk)
}
}
amsolve<-list(y=y) ## returns y
}
ndia3<-function(e, v, gv, r, vgv){
##---------------------------------------------------------
## update the preconditioning matrix based on a diagonal
## version of the bfgs quasi-newton update.
##---------------------------------------------------------
tol <- 1e-6
vr <- as.numeric(crossprod(v,r))
#% cat("ndia3: vgv, vr, then e:",vgv, vr,"\n")
#% print(e)
if (abs(vr)>tol && abs(vgv)>tol) {
e <- e - (r*r)/vr + (gv*gv)/vgv # CAUTION! May be crossprod
ind <- which(e < tol)
if (length(ind)>0){
e[ind] <- 1
}
}
#% print(e)
e # Need to return the object?? Possibly not!
}
norm2 <- function (x) {
##---------------------------------------------------------
## compute the 2 norm of vector x
##---------------------------------------------------------
res<-sqrt(as.numeric(crossprod(x)))
}
ssbfgs <- function (s, hv, hy, ys, yhy, vs, vhy){
##---------------------------------------------------------
## self-scaled bfgs quasi-Newton update
## (used by preconditioner)
##---------------------------------------------------------
delta <- (1 + yhy/ys)*vs/ys - vhy/ys
beta <- -vs/ys
z <- hv + delta*s + beta*hy
}
step1 <- function(f, gtp, smax) {
##---------------------------------------------------------
## step1 returns the length of the initial step to be
## taken along the vector p in the next linear search.
##---------------------------------------------------------
## [fm is supposed to be an estimate of the optimal function value]
eps<-.Machine$double.eps
fm <- 0
d <- abs(f-fm)
alpha <- min(1, smax)
if ((2*d <= (-gtp)) && (d >= eps)) {
alpha = min(-2*d/gtp, smax)
}
alpha # need to ensure returned value
}
stpmax <- function(stepmx, pe, x, p, ipivot, low, up) {
##------------------------------------------------
## compute the maximum allowable step length
## (spe is the standard (unconstrained) max step)
##------------------------------------------------
## cat("stpmax: stepmx, pe:",stepmx, pe,"\n")
spe <- stepmx / pe
al <- spe
au <- spe
##------------------------------------------------
indl <- which(ipivot==0 & p < 0)
if ( length(indl)>0) {
tl <- low[indl] - x[indl]
al <- min(tl/p[indl])
}
##------------------------------------------------
indu <- which(ipivot==0 & p > 0)
if (length(indu)>0) {
tu <- up[indu] - x[indu]
au <- min(tu/p[indu])
}
##------------------------------------------------
spe <- min(c(spe, al, au))
}
ztime<-function(x,ipivot) {
## ---------------------------------------------------------
## this routine multiplies the vector x by the
## constraint matrix z
## ---------------------------------------------------------
ind <- which(ipivot!=0);
x[ind] <- 0
x1 <- x;
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.