Nothing
R/solv-fpeq.R
In GNE: Computation of Generalized Nash Equilibria
Defines functions
extrapolFP
relaxationAlgoVH
relaxationAlgoUR
pureFP
fpeq
#############################################################################
# Copyright (c) 2012 Christophe Dutang
#
# 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.,
# 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
#
#############################################################################
### fixed-point equations in GNE
###
### R functions
###
fpeq <- function(xinit, fn, merit,
method=c("pure", "UR", "vH", "RRE", "MPE", "SqRRE", "SqMPE"),
control=list(), stepfunc, argstep, silent=TRUE, order.method=1, ...)
{
method <- match.arg(method, c("pure","UR", "vH", "RRE", "MPE", "SqRRE", "SqMPE"))
if(method %in% c("SqRRE", "SqMPE"))
{
order.method <- ifelse(order.method == 1, 2, order.method)
method <- substr(method, 3, 5)
}
if(method == "UR")
{
if(missing(fn) || missing(xinit) || missing(stepfunc) || missing(merit))
stop("missing parameters for UR method")
if(!missing(argstep))
finalstep <- function(x) stepfunc(x, argstep)
else
finalstep <- stepfunc
}
if(method == "vH")
if(missing(fn) || missing(xinit) || missing(merit))
stop("missing parameters for vH method")
noitercount <- FALSE
#inner iterations to compute fixed point function fn
inner.counts.fpfn <- c(0, 0) #call to gap and grad gap
inner.iter.fpfn <- 0 #iter number
#fixed point function
wrapfn <- function(x)
{
fx <- fn(x)
inner.counts.fpfn <<- inner.counts.fpfn + fx$counts
inner.iter.fpfn <<- inner.iter.fpfn + fx$iter
fx$par
}
if(!is.null(merit))
{
#inner iterations to compute merit function merit
inner.counts.merit <- c(0, 0) #call to gap and grad gap
inner.iter.merit <- 0 #iter number
wrapmerit <- function(x, ...)
{
vx <- merit(x, ...)
inner.counts.merit <<- inner.counts.merit + vx$counts
inner.iter.merit <<- inner.iter.merit + vx$iter
vx$value
}
}else
{
wrapmerit <- NULL
inner.counts.merit <- inner.iter.merit <- 0
}
#default control parameters
#for relaxationAlgoVH()
con <- list(sigma=9/10, beta=1/2, tol=1e-6, maxit=100, echo=0)
namc <- names(con)
con[namc <- names(control)] <- control
#for fpiter()
confpiter <- list(tol=1e-6, maxiter=100, trace=FALSE)
if(!is.null(control$tol))
confpiter$tol <- control$tol
if(!is.null(control$maxit))
confpiter$maxiter <- control$maxit
if(!is.null(control$trace))
confpiter$trace <- control$trace
#for squarem()
consquarem <- list(tol=1e-6, maxiter=100, trace=FALSE, K=order.method,
method=ifelse(order.method==1, 1*(method == "RRE")+2*(method == "MPE"), method))
if(!is.null(control$tol))
consquarem$tol <- control$tol
if(!is.null(control$maxit))
consquarem$maxiter <- control$maxit
if(!is.null(control$trace))
consquarem$trace <- control$trace
if(method == "pure" && !is.null(wrapmerit))
{
myGNE <- try( pureFP(xinit, wrapfn, wrapmerit, control=con, ...),
silent=silent)
}
if(method == "pure" && is.null(wrapmerit))
{
myGNE <- try( fpiter(xinit, wrapfn, wrapmerit, control=confpiter, ...),
silent=silent)
if(!is(myGNE,"try-error"))
{
if(!silent)
print(myGNE)
myGNE$value <- max(abs(fn(myGNE$par)$par))
myGNE$counts <- c(fn=myGNE$fpevals, merit=myGNE$objfevals)
myGNE$code <- 1*(myGNE$convergence == 0) + 4*(myGNE$fpevals[1] >= confpiter$maxiter)
myGNE$code <- myGNE$code + 2*(myGNE$convergence != 0 & myGNE$fpevals[1] <= consquarem$maxiter)
}
}
if(method == "UR")
{
myGNE <- try( relaxationAlgoUR(xinit, finalstep, wrapfn, wrapmerit,
control=con, ...), silent=silent)
}
if(method == "vH")
{
myGNE <- try( relaxationAlgoVH(xinit, wrapfn, wrapmerit,
control=con, ...), silent=silent)
}
if(method %in% c("RRE","MPE","SqRRE","SqMPE") && !is.null(wrapmerit))
{
myGNE <- try( extrapolFP(xinit, wrapfn, wrapmerit,
control=con, method=method, ...), silent=FALSE)
}
if(method %in% c("RRE","MPE","SqRRE","SqMPE") && is.null(wrapmerit))
{
myGNE <- try( squarem(xinit, wrapfn, control=consquarem, ...), silent=FALSE)
if(!is(myGNE,"try-error"))
{
if(!silent)
print(myGNE)
myGNE$value <- max(abs(fn(myGNE$par)$par))
myGNE$counts <- c(fn=myGNE$fpevals, merit=myGNE$objfevals)
myGNE$code <- 1*(myGNE$convergence == 0) + 4*(myGNE$fpevals[1] > consquarem$maxiter)
myGNE$code <- myGNE$code + 2*(myGNE$convergence != 0 & myGNE$fpevals[1] <= consquarem$maxiter)
}
}
if(!is(myGNE,"try-error"))
res <- list(par=myGNE$par, value=myGNE$value,
outer.counts=myGNE$counts, outer.iter=myGNE$counts[1],
code=myGNE$code, inner.iter=inner.iter.fpfn+inner.iter.merit,
inner.counts=inner.counts.fpfn+inner.counts.merit,
message=myGNE$message)
else
res <- list(par= NA, value=NA, outer.counts=NA, outer.iter=NA, code=100,
message=paste("Error in the fixed-point problem:", myGNE, "."),
inner.counts=NA, inner.iter=NA)
res
}
#pure fixed point iteration with a merit function
pureFP <- function(xinit, fn, merit, control, ...)
{
echo <- control$echo
maxit <- control$maxit
tol <- control$tol
k <- 1
xk_1 <- xinit
xk <- fn(xk_1)
xkp1 <- fn(xk)
merit_xk <- merit(xk, y=xkp1)
if(echo >= 1)
cat("**** k", k, "\n x_k", xk, "\n")
if(echo >= 2)
cat(" m(x_k)", merit_xk, "\n")
while( abs( merit_xk ) > tol && k < maxit)
{
xk_1 <- xk
k <- k+1
xk <- xkp1
xkp1 <- fn(xk)
merit_xk <- merit(xk, y=xkp1)
if(echo >= 1)
cat("**** k", k, "\n x_k", xk, "\n")
if(echo >= 2)
cat(" m(x_k)", merit_xk, "\n")
}
list(par = xk, value=merit_xk , counts=c(fn=k+1, merit=k+1), iter = k,
code=(k >= maxit)*4 + (abs(merit_xk) <= tol)*1 + (max(abs(xk - xk_1)/abs(xk)) <= tol)*2)
}
#Uryasev and Rubinstein (non optimized stepsize)
relaxationAlgoUR <- function(xinit, stepfunc, fn, merit, control, ...)
{
echo <- control$echo
maxit <- control$maxit
tol <- control$tol
k <- 1
xk_1 <- xinit
alphak <- stepfunc(k)
xk <- (1-alphak) * xk_1 + alphak * fn(xk_1)
if(!is.null(merit))
merit_xk <- merit(xk)
else
merit_xk <- max(abs(xk - xk_1))
if(echo >= 1)
cat("**** k", k, "\n x_k", xk, "\n")
if(echo >= 2)
cat(" m(x_k)", merit_xk, "\n")
if(echo >= 3)
cat("step size", alphak, "\n\n")
while( abs( merit_xk ) > tol && k < maxit)
{
xk_1 <- xk
k <- k+1
alphak <- stepfunc(k)
xk <- (1-alphak) * xk_1 + alphak * fn(xk_1)
if(!is.null(merit))
merit_xk <- merit(xk)
else
merit_xk <- max(abs(xk - xk_1))
if(echo >= 1)
cat("**** k", k, "\n x_k", xk, "\n")
if(echo >= 2)
cat(" m(x_k)", merit_xk, "\n")
if(echo >= 3)
cat("step size", alphak, "\n\n")
}
list(par = xk, value=merit_xk , counts=c(fn=k+1, merit=k+1), iter = k,
code=(k >= maxit)*4 + (abs(merit_xk) <= tol)*1 + (max(abs(xk - xk_1)/abs(xk)) <= tol)*2)
}
#von Heusinger
relaxationAlgoVH <- function(xinit, fn, merit, control, ...)
{
sigma <- control$sigma
beta <- control$beta
tol <- control$tol
echo <- control$echo
maxit <- control$maxit
k <- 0
xk_1 <- xinit
xk <- fn(xk_1)
merit_xk <- merit(xk)
counts <- c(fn=1, merit=1)
if(echo >= 1)
cat("**** k", k, "\n x_k", xk, "\n")
if(echo >= 2)
cat(" m(x_k)", merit_xk, "\n")
while( abs( merit_xk ) > tol && k < maxit)
{
k <- k+1
xk_1 <- xk
dk <- fn(xk) - xk
normdk <- sqrt(sum(dk^2))
merit_xk <- merit(xk)
#backtracking line search
tk <- 1
l <- 0
merit_xktkdk <- merit(xk + tk*dk)
# use remark 4.32
while( merit_xktkdk > merit_xk - sigma * tk^2 * normdk^2 )
{
tk <- tk * beta
merit_xktkdk <- merit(xk + tk*dk)
if(echo >= 3)
{
cat(l, "\t", merit_xktkdk, "\t <= ")
cat(merit_xk - sigma * tk^2 * normdk^2, "?\t")
cat("tk", tk, "\n\n")
}
l <- l+1
}
xk <- xk_1 + tk * dk
merit_xk <- merit_xktkdk
if(echo >= 1)
cat("**** k", k, "\n x_k", xk, "\n")
if(echo >= 2)
cat(" m(x_k)", merit_xk, "\n")
counts <- counts+1
counts[2] <- counts[2]+l
}
list(par = xk, value=merit_xk , counts=counts, iter = k,
code=(k >= maxit)*4 + (abs(merit_xk) <= tol)*1 + (max(abs(xk - xk_1)/abs(xk)) <= tol)*2)
}
#extrapolation method for fixed point iteration
extrapolFP <- function(xinit, fn, merit, control, method, ...)
{
echo <- control$echo
maxit <- control$maxit
tol <- control$tol
k <- 1
xk_1 <- xinit
xk <- fn(xk_1)
merit_xk <- merit(xk)
counts <- c(fn=1, merit=1)
if(echo >= 1)
cat("**** k", k, "\n x_k", xk, "\n")
if(echo >= 2)
cat(" m(x_k)", merit_xk, "\n")
while( abs( merit_xk ) > tol && k < maxit)
{
xk_1 <- xk
k <- k+1
xk <- fn(xk_1)
xkp1 <- fn(xk)
#RRE/MPE cycle of order 1, equivalent to Aitken Delta process when xk is univariate
Delta1_xk <- xk - xk_1
Delta2_xk <- xkp1 - 2*xk + xk_1
if(method == "RRE" || method == "SqRRE")
mystep <- as.vector(crossprod(Delta2_xk, Delta1_xk) / crossprod(Delta2_xk, Delta2_xk))
if(method == "MPE" || method == "SqMPE")
mystep <- as.vector(crossprod(Delta1_xk, Delta1_xk) / crossprod(Delta1_xk, Delta2_xk))
if(is.nan(mystep) || is.infinite(mystep))
{
warning("Delta square xk is too small: Use the last iterate.")
xk <- xkp1
}else if(substr(method, 1, 2) == "Sq")
{
xk <- xk_1 - 2*Delta1_xk * mystep + Delta2_xk * mystep^2
}else #simple version
{
xk <- xk_1 - Delta1_xk * mystep
}
merit_xk <- merit(xk)
if(echo >= 1)
cat("**** k", k, "\n x_k", xk, "\n")
if(echo >= 2)
cat(" m(x_k)", merit_xk, "\n")
if(echo >= 3)
cat(" Delta(xk)", Delta1_xk, "\n", "Delta^2(xk)", Delta2_xk, "\n\n")
counts <- counts + c(2,1)
}
list(par = xk, value=merit_xk , counts=counts, iter = k,
code=(k >= maxit)*4 + (abs(merit_xk) <= tol)*1 + (max(abs(xk - xk_1)/abs(xk)) <= tol)*2)
}
Try the GNE package in your browser
Any scripts or data that you put into this service are public.
GNE documentation built on March 31, 2023, 9:25 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions extrapolFP relaxationAlgoVH relaxationAlgoUR pureFP fpeq
#############################################################################
# Copyright (c) 2012 Christophe Dutang
#
# 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.,
# 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
#
#############################################################################
### fixed-point equations in GNE
###
### R functions
###
fpeq <- function(xinit, fn, merit,
method=c("pure", "UR", "vH", "RRE", "MPE", "SqRRE", "SqMPE"),
control=list(), stepfunc, argstep, silent=TRUE, order.method=1, ...)
{
method <- match.arg(method, c("pure","UR", "vH", "RRE", "MPE", "SqRRE", "SqMPE"))
if(method %in% c("SqRRE", "SqMPE"))
{
order.method <- ifelse(order.method == 1, 2, order.method)
method <- substr(method, 3, 5)
}
if(method == "UR")
{
if(missing(fn) || missing(xinit) || missing(stepfunc) || missing(merit))
stop("missing parameters for UR method")
if(!missing(argstep))
finalstep <- function(x) stepfunc(x, argstep)
else
finalstep <- stepfunc
}
if(method == "vH")
if(missing(fn) || missing(xinit) || missing(merit))
stop("missing parameters for vH method")
noitercount <- FALSE
#inner iterations to compute fixed point function fn
inner.counts.fpfn <- c(0, 0) #call to gap and grad gap
inner.iter.fpfn <- 0 #iter number
#fixed point function
wrapfn <- function(x)
{
fx <- fn(x)
inner.counts.fpfn <<- inner.counts.fpfn + fx$counts
inner.iter.fpfn <<- inner.iter.fpfn + fx$iter
fx$par
}
if(!is.null(merit))
{
#inner iterations to compute merit function merit
inner.counts.merit <- c(0, 0) #call to gap and grad gap
inner.iter.merit <- 0 #iter number
wrapmerit <- function(x, ...)
{
vx <- merit(x, ...)
inner.counts.merit <<- inner.counts.merit + vx$counts
inner.iter.merit <<- inner.iter.merit + vx$iter
vx$value
}
}else
{
wrapmerit <- NULL
inner.counts.merit <- inner.iter.merit <- 0
}
#default control parameters
#for relaxationAlgoVH()
con <- list(sigma=9/10, beta=1/2, tol=1e-6, maxit=100, echo=0)
namc <- names(con)
con[namc <- names(control)] <- control
#for fpiter()
confpiter <- list(tol=1e-6, maxiter=100, trace=FALSE)
if(!is.null(control$tol))
confpiter$tol <- control$tol
if(!is.null(control$maxit))
confpiter$maxiter <- control$maxit
if(!is.null(control$trace))
confpiter$trace <- control$trace
#for squarem()
consquarem <- list(tol=1e-6, maxiter=100, trace=FALSE, K=order.method,
method=ifelse(order.method==1, 1*(method == "RRE")+2*(method == "MPE"), method))
if(!is.null(control$tol))
consquarem$tol <- control$tol
if(!is.null(control$maxit))
consquarem$maxiter <- control$maxit
if(!is.null(control$trace))
consquarem$trace <- control$trace
if(method == "pure" && !is.null(wrapmerit))
{
myGNE <- try( pureFP(xinit, wrapfn, wrapmerit, control=con, ...),
silent=silent)
}
if(method == "pure" && is.null(wrapmerit))
{
myGNE <- try( fpiter(xinit, wrapfn, wrapmerit, control=confpiter, ...),
silent=silent)
if(!is(myGNE,"try-error"))
{
if(!silent)
print(myGNE)
myGNE$value <- max(abs(fn(myGNE$par)$par))
myGNE$counts <- c(fn=myGNE$fpevals, merit=myGNE$objfevals)
myGNE$code <- 1*(myGNE$convergence == 0) + 4*(myGNE$fpevals[1] >= confpiter$maxiter)
myGNE$code <- myGNE$code + 2*(myGNE$convergence != 0 & myGNE$fpevals[1] <= consquarem$maxiter)
}
}
if(method == "UR")
{
myGNE <- try( relaxationAlgoUR(xinit, finalstep, wrapfn, wrapmerit,
control=con, ...), silent=silent)
}
if(method == "vH")
{
myGNE <- try( relaxationAlgoVH(xinit, wrapfn, wrapmerit,
control=con, ...), silent=silent)
}
if(method %in% c("RRE","MPE","SqRRE","SqMPE") && !is.null(wrapmerit))
{
myGNE <- try( extrapolFP(xinit, wrapfn, wrapmerit,
control=con, method=method, ...), silent=FALSE)
}
if(method %in% c("RRE","MPE","SqRRE","SqMPE") && is.null(wrapmerit))
{
myGNE <- try( squarem(xinit, wrapfn, control=consquarem, ...), silent=FALSE)
if(!is(myGNE,"try-error"))
{
if(!silent)
print(myGNE)
myGNE$value <- max(abs(fn(myGNE$par)$par))
myGNE$counts <- c(fn=myGNE$fpevals, merit=myGNE$objfevals)
myGNE$code <- 1*(myGNE$convergence == 0) + 4*(myGNE$fpevals[1] > consquarem$maxiter)
myGNE$code <- myGNE$code + 2*(myGNE$convergence != 0 & myGNE$fpevals[1] <= consquarem$maxiter)
}
}
if(!is(myGNE,"try-error"))
res <- list(par=myGNE$par, value=myGNE$value,
outer.counts=myGNE$counts, outer.iter=myGNE$counts[1],
code=myGNE$code, inner.iter=inner.iter.fpfn+inner.iter.merit,
inner.counts=inner.counts.fpfn+inner.counts.merit,
message=myGNE$message)
else
res <- list(par= NA, value=NA, outer.counts=NA, outer.iter=NA, code=100,
message=paste("Error in the fixed-point problem:", myGNE, "."),
inner.counts=NA, inner.iter=NA)
res
}
#pure fixed point iteration with a merit function
pureFP <- function(xinit, fn, merit, control, ...)
{
echo <- control$echo
maxit <- control$maxit
tol <- control$tol
k <- 1
xk_1 <- xinit
xk <- fn(xk_1)
xkp1 <- fn(xk)
merit_xk <- merit(xk, y=xkp1)
if(echo >= 1)
cat("**** k", k, "\n x_k", xk, "\n")
if(echo >= 2)
cat(" m(x_k)", merit_xk, "\n")
while( abs( merit_xk ) > tol && k < maxit)
{
xk_1 <- xk
k <- k+1
xk <- xkp1
xkp1 <- fn(xk)
merit_xk <- merit(xk, y=xkp1)
if(echo >= 1)
cat("**** k", k, "\n x_k", xk, "\n")
if(echo >= 2)
cat(" m(x_k)", merit_xk, "\n")
}
list(par = xk, value=merit_xk , counts=c(fn=k+1, merit=k+1), iter = k,
code=(k >= maxit)*4 + (abs(merit_xk) <= tol)*1 + (max(abs(xk - xk_1)/abs(xk)) <= tol)*2)
}
#Uryasev and Rubinstein (non optimized stepsize)
relaxationAlgoUR <- function(xinit, stepfunc, fn, merit, control, ...)
{
echo <- control$echo
maxit <- control$maxit
tol <- control$tol
k <- 1
xk_1 <- xinit
alphak <- stepfunc(k)
xk <- (1-alphak) * xk_1 + alphak * fn(xk_1)
if(!is.null(merit))
merit_xk <- merit(xk)
else
merit_xk <- max(abs(xk - xk_1))
if(echo >= 1)
cat("**** k", k, "\n x_k", xk, "\n")
if(echo >= 2)
cat(" m(x_k)", merit_xk, "\n")
if(echo >= 3)
cat("step size", alphak, "\n\n")
while( abs( merit_xk ) > tol && k < maxit)
{
xk_1 <- xk
k <- k+1
alphak <- stepfunc(k)
xk <- (1-alphak) * xk_1 + alphak * fn(xk_1)
if(!is.null(merit))
merit_xk <- merit(xk)
else
merit_xk <- max(abs(xk - xk_1))
if(echo >= 1)
cat("**** k", k, "\n x_k", xk, "\n")
if(echo >= 2)
cat(" m(x_k)", merit_xk, "\n")
if(echo >= 3)
cat("step size", alphak, "\n\n")
}
list(par = xk, value=merit_xk , counts=c(fn=k+1, merit=k+1), iter = k,
code=(k >= maxit)*4 + (abs(merit_xk) <= tol)*1 + (max(abs(xk - xk_1)/abs(xk)) <= tol)*2)
}
#von Heusinger
relaxationAlgoVH <- function(xinit, fn, merit, control, ...)
{
sigma <- control$sigma
beta <- control$beta
tol <- control$tol
echo <- control$echo
maxit <- control$maxit
k <- 0
xk_1 <- xinit
xk <- fn(xk_1)
merit_xk <- merit(xk)
counts <- c(fn=1, merit=1)
if(echo >= 1)
cat("**** k", k, "\n x_k", xk, "\n")
if(echo >= 2)
cat(" m(x_k)", merit_xk, "\n")
while( abs( merit_xk ) > tol && k < maxit)
{
k <- k+1
xk_1 <- xk
dk <- fn(xk) - xk
normdk <- sqrt(sum(dk^2))
merit_xk <- merit(xk)
#backtracking line search
tk <- 1
l <- 0
merit_xktkdk <- merit(xk + tk*dk)
# use remark 4.32
while( merit_xktkdk > merit_xk - sigma * tk^2 * normdk^2 )
{
tk <- tk * beta
merit_xktkdk <- merit(xk + tk*dk)
if(echo >= 3)
{
cat(l, "\t", merit_xktkdk, "\t <= ")
cat(merit_xk - sigma * tk^2 * normdk^2, "?\t")
cat("tk", tk, "\n\n")
}
l <- l+1
}
xk <- xk_1 + tk * dk
merit_xk <- merit_xktkdk
if(echo >= 1)
cat("**** k", k, "\n x_k", xk, "\n")
if(echo >= 2)
cat(" m(x_k)", merit_xk, "\n")
counts <- counts+1
counts[2] <- counts[2]+l
}
list(par = xk, value=merit_xk , counts=counts, iter = k,
code=(k >= maxit)*4 + (abs(merit_xk) <= tol)*1 + (max(abs(xk - xk_1)/abs(xk)) <= tol)*2)
}
#extrapolation method for fixed point iteration
extrapolFP <- function(xinit, fn, merit, control, method, ...)
{
echo <- control$echo
maxit <- control$maxit
tol <- control$tol
k <- 1
xk_1 <- xinit
xk <- fn(xk_1)
merit_xk <- merit(xk)
counts <- c(fn=1, merit=1)
if(echo >= 1)
cat("**** k", k, "\n x_k", xk, "\n")
if(echo >= 2)
cat(" m(x_k)", merit_xk, "\n")
while( abs( merit_xk ) > tol && k < maxit)
{
xk_1 <- xk
k <- k+1
xk <- fn(xk_1)
xkp1 <- fn(xk)
#RRE/MPE cycle of order 1, equivalent to Aitken Delta process when xk is univariate
Delta1_xk <- xk - xk_1
Delta2_xk <- xkp1 - 2*xk + xk_1
if(method == "RRE" || method == "SqRRE")
mystep <- as.vector(crossprod(Delta2_xk, Delta1_xk) / crossprod(Delta2_xk, Delta2_xk))
if(method == "MPE" || method == "SqMPE")
mystep <- as.vector(crossprod(Delta1_xk, Delta1_xk) / crossprod(Delta1_xk, Delta2_xk))
if(is.nan(mystep) || is.infinite(mystep))
{
warning("Delta square xk is too small: Use the last iterate.")
xk <- xkp1
}else if(substr(method, 1, 2) == "Sq")
{
xk <- xk_1 - 2*Delta1_xk * mystep + Delta2_xk * mystep^2
}else #simple version
{
xk <- xk_1 - Delta1_xk * mystep
}
merit_xk <- merit(xk)
if(echo >= 1)
cat("**** k", k, "\n x_k", xk, "\n")
if(echo >= 2)
cat(" m(x_k)", merit_xk, "\n")
if(echo >= 3)
cat(" Delta(xk)", Delta1_xk, "\n", "Delta^2(xk)", Delta2_xk, "\n\n")
counts <- counts + c(2,1)
}
list(par = xk, value=merit_xk , counts=counts, iter = k,
code=(k >= maxit)*4 + (abs(merit_xk) <= tol)*1 + (max(abs(xk - xk_1)/abs(xk)) <= tol)*2)
}
Try the GNE package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.