Nothing
R/Data2LD.opt.R
In Data2LD: Functional Data Analysis with Linear Differential Equations
Defines functions
stepit
stepchkData2LD
Data2LD.opt
Documented in
Data2LD.opt
Data2LD.opt <- function(yList, XbasisList, modelList, rhoMat,
convcrit=1e-4, iterlim=20, dbglev=1,
parMap=diag(rep(1,nparam))) {
# Data2LD.opt optimizes a parameter vector theta defining a
# linear differential operator object used to smooth a set of data.
#
# Arguments:
# YLIST an array containing values of curves
# If the array is a matrix, rows must correspond to argument
# values and columns to replications, and it will be assumed
# that there is only one variable per observation.
# If Y is a three-dimensional array, the first dimension
# corresponds to argument values, the second to replications,
# and the third to variables within replications.
# If Y is a vector, only one replicate and variable are
# assumed.
# BASISLIST A functional data object or a BASIS object. If so, the
# smoothing parameter LAMBDA is set to 0.
# MODELLIST A list of length NVAR. Each list contains a
# struct object with members:
# XList list of length number of homogeneous terms
# Each list contains a struct object with members:
# WfdPar a fdPar object for the coefficient
# variable the index of the variable
# derivative the order of its derivative
# npar if coefficient estimated, its location
# in the composite vector
# FList list of length number of forcing terms
# Each list contains a struct object with members:
# AfdPar an fdPar object for the coefficient
# Ufd an fd object for the forcing function
# npar if coefficient estimated, its location
# in the composite vector
# order the highest order of derivative
# name a tag for the variable
# nallXterm the number of homogeneous terms
# nallFterm the number of forcing functions
# RHOMAT A matrix with the number of rows equal to the number of values of the
# smoothing parameter per variable to be used, and number of columns
# equal to the number of variables.
# Each entry in the matrix is a value \rho in the interval [0,1)
# so that 0 <= \rho < 1.
# For each variable ivar in the system of equations and each optimization
# iopt, the data are weighted by \rho(iopt,ivar) and the
# roughness penalty by 1-\rho(iopt,ivar).
# It is expected that the values of \rho within each column will be
# in ascending order, and the estimated parameters for each row
# are passed along as initial values to be used for the optimization
# defined by the values of \rho in the next row.
# CONVCRIT One convergence criterion, or a vector of two criteria.
# The first criterion is applied to the function change,
# and the second is applied to the gradient norm.
# ITERLIM Maximum number of iterations allowed.
# DBGLEV An integer controlling amount of output per iteration.
# Defaults to 1, which prints summary results for
# each iteration.
# ACTIVE A vector of indices of the parameters to be optimized. If the index
# of a parameter does not appear in this vector, it is left fixed
# at its initial value even if its specification in coefList is that
# it be estimated. If the parameters are constrained by K linear
# constraints, ACTIVE contains integers from 1 to NPAR - K, where
# NPAR is the number of parameters in the model. In other words,
# ACTIVE specifies the constrained parameters to be estimated.
# PARMAP A rectangular matrix with number of rows equal to the
# number of parameters to be estimated as defined in
# BWTLIST, and number of columns equal to the number of
# parameters less the number of linear constraints on the
# estimated parameters. The columns of PARMAP must be
# orthnormal so that t(PARMAP) %*% PARMAP is an identity matrix.
# t(PARMAP) %*% THETA maps unconstrained parameters and the
# corresponding gradient into constrained parameter space.
# PARMAP %*% THETA maps constrained parameters and the
# corresponding gradient into unconstrained parameter space.
# PARMAP will usually be set up using the full QR
# decomposition of a linear constraint coefficient matrix t(A)
# where the constraints are of the form A P <- B, A and B
# being known matrices. An example of such a constraint
# that arises often is one where two estimated coefficients
# are constrained to be equal. For example, if a vvariable
# X involved in an equation the form a(x - x.0), where x.0
# is a fixed set point or defined target level for variable
# X, then this would be set up as a.1 x + a.2 x.0, where
# coefficients a.1 and a.2 are constrained to be equal in
# magnitude but opposite in sign, or a.1 + a.2 <- 0.
# Returns: A named list object with these fields:
# THETA.OPT The optimal parameter values.
# BWTLIST.OPT The corresponding optimal coefficients for the
# homogeneous terms in the differential equation
# AWTLIST.OPT The corresponding optimal coefficients for the
# forcing terms in the differential equation
# Last modified 23 April 2020
# ------------------------------------------------------------------------
# Check input parameters
# ------------------------------------------------------------------------
nvar <- length(yList)
nparam <- 0
for (ivar in 1:nvar) {
modelListi <- modelList[[ivar]]
if (modelListi$nallXterm > 0) {
for (iw in 1:modelListi$nallXterm) {
modelListiw <- modelListi$XList[[iw]]
nparam <- nparam + length(modelListiw$index)
}
}
if (modelListi$nallFterm > 0) {
for (jforce in 1:modelListi$nallFterm) {
modelListj <- modelListi$FList[[jforce]]
nparam = nparam + length(modelListj$index)
}
}
}
# Get the vector of parameters
theta <- modelList2Vec(modelList, nparam)
ntheta <- length(theta)
thetaCon <- t(parMap) %*% theta
nparCon <- length(thetaCon)
climit <- matrix(c(-1000,1000),2,1) %*% matrix(1,1,nparCon)
active <- 1:nparCon
# check rhoMat and get nopt
if (!is.matrix(rhoMat)) {
rhoMat = matrix(rhoMat,length(rhoMat),1)
}
rhodim <- dim(rhoMat)
if (rhodim[1] != nvar) {
stop(paste("The first dimension of RHOMAT is not equal to",
"the number of variables"))
}
nopt <- rhodim[2]
# Lists to contain results if more than a single set of rho's are involved
thetastore <- matrix(0,ntheta,nopt)
dfstore <- matrix(0,nopt,1)
gcvstore <- matrix(0,nopt,1)
modelList.optList <- vector("list", nopt)
modelListnew <- modelList
# ------------------------------------------------------------------------
# loop through rho vectors
# ------------------------------------------------------------------------
for (iopt in 1:nopt) {
rhoVeci <- rhoMat[,iopt]
# compute initial criterion value, gradient and hessian
Data2LDList <- Data2LD(yList, XbasisList, modelListnew, rhoVeci)
fvec <- Data2LDList$MSE
grad <- Data2LDList$DpMSE
gradCon <- t(parMap) %*% grad
fvecsum <- sum(fvec)
norm <- sqrt(mean(gradCon^2))
# evaluate the initial update vector for correcting the initial theta
deltac <- -gradCon
# initialize iteration status arrays
iternum <- 0
status <- c(iternum, fvec, norm)
if (dbglev >= 1) {
cat("\nIter. Criterion Grad Length")
cat("\n")
# # cat(format(c(iternum,status[2:3]),digits=4))
# cat(iternum)
# cat(" ")
# cat(format(status[2],nsmall=4,digits=4))
# cat(" ")
# cat(format(status[3],nsmall=4,digits=4))
# cat("\n")
cat(iternum)
cat(" ")
cat(round(status[2],6))
cat(" ")
cat(round(status[3],6))
} else {
cat(".")
}
cat("\n")
iterhist <- matrix(0,iterlim+1,length(status))
iterhist[1,] <- status
if (iterlim == 0) {
return ()
}
# ------- Begin main iterations -----------
MAXSTEPITER <- 5
MAXSTEP <- 1000
trial <- 0.1
reset <- 0
linemat <- matrix(0,3,5)
thetaoldCon <- thetaCon
fvecsumold <- fvecsum
gradoldCon <- gradCon
dbgwrd <- dbglev >= 2
# --------------- beginning of optimization loop -----------
for (iter in 1:iterlim) {
iternum <- iternum + 1
# set logical parameters
dblwrd <- c(FALSE,FALSE)
limwrd <- c(FALSE,FALSE)
stpwrd <- FALSE
ind <- 0
ips <- 0
# compute slope
linemat[2,1] <- sum(deltac*gradoldCon)
# normalize search direction vector
sdg <- sqrt(sum(deltac^2))
deltac <- deltac/sdg
linemat[2,1] <- linemat[2,1]/sdg
# initialize line search vectors
linemat[,1:4] <- matrix(c(0, linemat[2,1], fvecsum),3,1) %*% matrix(1,1,4)
stepiter <- 0
if (dbglev >= 2) {
cat(paste(" ", stepiter, " "))
cat(format(round(t(linemat[,1]),6)))
cat("\n")
}
# return with error condition if initial slope is nonnegative
if (linemat[2,1] >= 0) {
if (dbglev >= 2) {
cat("Initial slope ")
cat(format(round(linemat[2,1])))
cat(" is nonnegative.\n")
cat("\n")
}
deltac <- -gradCon
linemat[2,1] <- -sum(gradCon^2)
sdg <- sqrt(sum(deltac^2))
deltac <- deltac/sdg
linemat[,1:4] <- matrix(c(0, linemat[2,1], fvecsum),3,1) %*% matrix(1,1,4)
linestore <- matrix(linemat[,1],3,1)
if (dbglev > 2) {
cat(paste(" ", stepiter, " "))
cat(format(round(t(linemat[,1]),6)))
cat("\n")
}
}
# return successfully if initial slope is very small
if (linemat[2,1] >= -min(c(1e-3,convcrit))) {
if (dbglev >= 2) {
cat("Initial slope too small.\n")
status <- c(iternum, fvecsum, norm)
cat("\n")
}
modelListnew <- modelList
break
}
# first step set to trial
if (is.nan(trial)) stop("is.nan(trial)")
linemat[1,5] <- trial
# ------------ begin line search iteration loop ----------
thetanewCon <- thetaCon
for (stepiter in 1:MAXSTEPITER) {
# check the step size and modify if limits exceeded
stepval <- linemat[1,5]
result <- stepchkData2LD(linemat[1,5], thetaCon, deltac, limwrd, ind,
climit, active, dbgwrd)
linemat[1,5] <- result[[1]]
ind <- result[[2]]
limwrd <- result[[3]]
# break if limit hit twice in a row
if (ind == 1) {
thetanew <- theta
fvecsumnew <- fvecsum
break
}
# break if current step size too small
if (linemat[1,5] <= 1e-7) {
if (dbglev >= 2) {
cat("\nStepsize too small: ")
cat(linemat[1,5])
cat("\n")
}
thetanew <- theta
fvecsumnew <- fvecsum
break
}
# update parameter vector
thetanewCon <- thetaCon + linemat[1,5]*deltac
# --------- update function, gradient and hessian -----------
thetanew <- parMap %*% thetanewCon
modelListnew <- modelVec2List(modelList, thetanew)
Data2LDList <- Data2LD(yList, XbasisList, modelListnew,
rhoVeci, summary=FALSE)
fvecnew <- Data2LDList$MSE
gradnew <- Data2LDList$DpMSE
hessmatnew <- Data2LDList$D2ppMSE
gradnewCon <- t(parMap) %*% gradnew
hessmatnewCon <- t(parMap) %*% hessmatnew %*% parMap
fvecsumnew <- sum(fvecnew)
# -------------------------------------------------------------
linemat[3,5] <- fvecsumnew
# compute new directional derivative
linemat[2,5] <- sum(deltac*gradnewCon)
if (dbglev >= 2) {
cat(paste(" ", stepiter, " "))
cat(format(round(t(linemat[,5]),6)))
cat("\n")
}
# compute next line search step, also testing for convergence
result <- stepit(linemat, ips, dblwrd, MAXSTEP)
linemat <- result[[1]]
ips <- result[[2]]
ind <- result[[3]]
dblwrd <- result[[4]]
trial <- linemat[1,5]
# ind == 0 implies convergence
if (ind == 0 || ind == 5) {
break
}
}
# ------------ end line search iteration loop ----------
thetaCon <- thetanewCon
theta <- thetanew
fvecsum <- fvecsumnew
gradCon <- gradnewCon
# test for function value made worse
if (fvecsum > fvecsumold) {
# Function value worse warn and terminate
if (dbglev >= 2) {
cat("\n")
cat("Criterion increased: ")
cat(format(round(c(fvecsumold, fvecsum),4)))
cat("\n")
}
# reset parameters and fit
thetaCon <- thetaoldCon
fvecsum <- fvecsumold
deltac <- gradCon
if (reset == 1) {
# This is the second time in a row that this
# has happened quit
if (dbglev >= 2) {
cat("Reset twice, terminating.\n")
cat("Convergence not attained.\n")
}
# return current status of optimization
break
} else {
reset <- 1
}
} else {
# function value has not increased, check for convergence
RMSgrad <- sqrt(mean(gradCon^2))
if (length(convcrit) > 1) {
convtest <- fvecsumold - fvecsum < convcrit[1] && RMSgrad < convcrit[2]
} else {
convtest <- fvecsumold - fvecsum < convcrit
}
if (convtest) {
norm <- sqrt(mean(gradCon^2))
status <- c(iternum, fvecsum, norm)
if (dbglev >= 1) {
cat(iternum)
cat(" ")
cat(round(status[2],6))
cat(" ")
cat(round(status[3],6))
cat("\n")
cat("Convergence reached.\n")
}
# return current status of optimization
break
}
# update old parameter vectors and fit structure
thetaoldCon <- thetaCon
fvecsumold <- fvecsum
gradoldCon <- gradCon
hessmatCon <- t(parMap) %*% hessmatnew %*% parMap
# update the line search direction vector
deltac <- -solve(hessmatCon,gradCon)
reset <- 0
}
norm <- sqrt(mean(gradCon^2))
status <- c(iternum, fvecsum, norm)
iterhist[iter+1,] <- status
if (dbglev >= 1) {
cat(iternum)
cat(" ")
cat(round(status[2],6))
cat(" ")
cat(round(status[3],6))
cat("\n")
}
if (dbglev >= 1 && iter == iterlim) {
cat(
"\nMaximum iterations reached but convergence not attained.\n")
}
}
# --------------- end of optimization loop -----------
cat("\n")
# evaluate the solution at the final optimal value of theta
Data2LDList <- Data2LD(yList, XbasisList, modelListnew, rhoVeci, summary=TRUE)
# store the optimal values for this level of the \rho's the
# storage objects
theta <- modelList2Vec(modelListnew, nparam)
thetastore[,iopt] <- theta
dfstore[iopt] <- Data2LDList$df
gcvstore[iopt] <- Data2LDList$gcv
modelList.optList[[iopt]] <- modelListnew
}
# ------------------- end of rho loop ------------------
return(list(thetastore=thetastore, dfstore=dfstore, gcvstore=gcvstore,
modelList.optList=modelList.optList))
}
# ----------------------------------------------------------------------------
stepchkData2LD = function(stepval, cvec, deltac, limwrd, ind,
climit=50*matrix(rep(1,2),2,1) %*% matrix(1,1,npar),
active=1:npar, dbgwrd=FALSE) {
#STEPCHK checks the step along a line for producing parameters within the
# limits specified by BOT and TOP
# LIMWRD Logical variable permitting detection that parameter
# was on the boundary two steps in a row
# Last modified 28 August 2017
npar = length(deltac)
bot = climit[1,]
top = climit[2,]
# if step is too small, return with flag ind = 1
if (stepval < 1e-7) {
ind = 1
return(list(stepval=stepval,ind=ind,limwrd=limwrd))
}
# ensure that step does not go beyond lower limit on parameters
stepvali = stepval*deltac
if (any(stepvali[active] < bot[active]-cvec[active])) {
index = active[stepvali[active] < bot[active]-cvec[active] &
deltac[active] != 0]
stepvalnew = min((bot[index]-cvec[index])/deltac[index])
stepval = stepvalnew
# check whether lower limit has been reached twice in a row
if (limwrd[1]) {
ind = 1
return(list(stepval=stepval,ind=ind,limwrd=limwrd))
} else {
limwrd[1] = TRUE
}
} else {
limwrd[1] = FALSE
}
if (stepval < 1e-7) {
ind = 1
return(list(stepval=stepval,ind=ind,limwrd=limwrd))
}
# ensure that step does not go beyond upper limit on parameters
stepvali = stepval*deltac
if (any(stepvali[active] > top[active]-cvec[active])) {
index = active[stepvali[active] > top[active]-cvec[active] &
deltac[active] != 0]
stepvalnew = min((top[index]-cvec[index])/deltac[index])
stepval = stepvalnew
# check whether upper limit has been reached twice in a row
if (limwrd[2]) {
ind = 1
return(list(stepval=stepval,ind=ind,limwrd=limwrd))
} else {
limwrd[2] = TRUE
}
} else {
limwrd[2] = FALSE
}
return(list(stepval=stepval,ind=ind,limwrd=limwrd))
}
# ----------------------------------------------------------------------------
stepit <- function(linemat, ips, dblwrd, MAXSTEP) {
#STEPIT computes next step size in line search algorithm
# Arguments:
# LINEMAT: Row 1 contains step values
# Row 2 contains slope values
# Row 3 contains function values
# IPS: If 1, previous slope was positive
# DBLWRD: Vector of length 2: dblwrd[1] T means step halved
# dblwrd[2] T means step doubled
# MAXSTEP: maximum size of step
# Last modified 4 March 2014 by Jim Ramsay
# Wolfe condition 1
test1.1 = linemat[3,5] <= linemat[3,1] + linemat[1,5]*linemat[2,1]/20
# Wolfe condition 2
test1.2 = abs(linemat[2,5]) <= abs(linemat[2,1])/10
# disp([test1.1, test1.2])
test1 = test1.1 && test1.2
# test1 = test1.2
test2 = linemat[3,5] > linemat[3,1]
test3 = linemat[2,5] > 0
if ((test1 || !test3) && test2) {
# ************************************************************
# function is worse and either slope is satisfory or negative
ips = 0 # step is halved
if (dblwrd[2]) {
ind = 5
return(list(linemat = linemat, ips = ips, ind = ind, dblwrd = dblwrd))
}
linemat[1,5] = min(c(linemat[1,5]/2, MAXSTEP))
linemat[,2] = linemat[,1]
linemat[,3] = linemat[,1]
dblwrd = c(1, 0)
ind = 2
return(list(linemat = linemat, ips = ips, ind = ind, dblwrd = dblwrd))
}
# *********************************************************
if (test1) {
# test1 means successful convergence
ind = 0
return(list(linemat = linemat, ips = ips, ind = ind, dblwrd = dblwrd))
}
# **********************************************************
if (test3) {
# Current slope is positive
ips = 1
linemat[,4] = linemat[,5]
deltaf = linemat[3,3] - linemat[3,5]
z = (3/(linemat[1,5] - linemat[1,3]))*deltaf + linemat[2,3] + linemat[2,5]
w = z * z - linemat[2,3] * linemat[2,5]
if (abs(linemat[2,3] + linemat[2,5] + 2 * z) >= 1e-05 && w > 0) {
if (w < 0) stop("w < 0")
w = sqrt(w)
denom = linemat[2,5] - linemat[2,3] + 2 * w
if (denom == 0) stop("denom == 0")
linemat[1,5] = linemat[1,3] + (1 - ((linemat[2,5] + w - z)/
(denom))) * (linemat[1,5] - linemat[1,3])
if (is.nan(linemat[1,5])) stop("first test3: is.nan(linemat[1,5])")
} else {
# linear interpolation necessary
aerror = linemat[1,3]
if (linemat[1,5] > linemat[1,3]) {
aerror = linemat[1,5]
}
linemat[1,5] = linemat[1,3] - linemat[2,3] * ((linemat[1,5] - linemat[1,3])/
(linemat[2,5] - linemat[2,3]))
if (linemat[1,5] > 2 * aerror) {
linemat[1,5] = 2 * aerror
}
}
linemat[1,5] = min(c(linemat[1,5], MAXSTEP))
dblwrd = c(0,0)
ind = 2
return(list(linemat = linemat, ips = ips, ind = ind, dblwrd = dblwrd))
}
# *************************************************************
# Current slope is negative or zero
linemat[,2] = linemat[,3]
linemat[,3] = linemat[,5]
if (ips == 1) {
# *****************************************************
# previous slope is positive
deltaf = linemat[3,5] - linemat[3,4]
z = c(3/(linemat[1,4] - linemat[1,5])) * deltaf +
linemat[2,5] + linemat[2,4]
w = z * z - linemat[2,5] * linemat[2,4]
if (abs(linemat[2,5] + linemat[2,4] + 2 * z) >= 1e-05 && w > 0) {
w = sqrt(w)
linemat[1,5] = linemat[1,5] + (1 - ((linemat[2,4] + w - z)/
(linemat[2,4] - linemat[2,5] +
2 * w))) * (linemat[1,4] - linemat[1,5])
} else {
aerror = linemat[1,5]
if (linemat[1,4] > linemat[1,5]) {
aerror = linemat[1,4]
}
linemat[1,5] = linemat[1,5] - linemat[2,5] *
((linemat[1,4] - linemat[1,5])/
(linemat[2,4] - linemat[2,5]))
if (linemat[1,5] > 2 * aerror) {
linemat[1,5] = 2 * aerror
}
}
linemat[1,5] = min(c(linemat[1,5], MAXSTEP))
dblwrd = c(0,0)
ind = 2
return(list(linemat = linemat, ips = ips, ind = ind, dblwrd = dblwrd))
}
# ******************************************************
if ((linemat[2,3] - linemat[2,2]) * (linemat[1,3] - linemat[1,2]) > 0) {
# previous slope is negative
z = c(3/(linemat[1,3] - linemat[1,2])) * (linemat[3,2] - linemat[3,3]) +
linemat[2,2] + linemat[2,3]
w = z * z - linemat[2,2] * linemat[2,3]
if (abs(linemat[2,2] + linemat[2,3] + 2 * z) >= 1e-05 && w > 0) {
w = sqrt(w)
linemat[1,5] = linemat[1,2] +
(1 - ((linemat[2,3] + w - z)/(linemat[2,3] - linemat[2,2] +
2 * w))) * (linemat[1,3] - linemat[1,2])
} else {
linemat[1,5] = linemat[1,2] - linemat[2,2] *
((linemat[1,3] - linemat[1,2])/
(linemat[2,3] - linemat[2,2]))
}
linemat[1,5] = min(c(linemat[1,5], MAXSTEP))
dblwrd = c(0,0)
ind = 2
return(list(linemat = linemat, ips = ips, ind = ind, dblwrd = dblwrd))
} else {
# previous slope also negative but not as much
if (dblwrd[1]) {
ind = 5
return(
list(linemat = linemat, ips = ips, ind = ind, dblwrd = dblwrd))
} else {
linemat[1,5] = 2 * linemat[1,5]
linemat[1,5] = min(c(linemat[1,5], MAXSTEP))
dblwrd = c(0,1)
ind = 2
return(
list(linemat = linemat, ips = ips, ind = ind, dblwrd = dblwrd))
}
}
ind = 2
}
Try the Data2LD package in your browser
Any scripts or data that you put into this service are public.
Data2LD documentation built on Aug. 6, 2020, 1:08 a.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 stepit stepchkData2LD Data2LD.opt
Documented in Data2LD.opt
Data2LD.opt <- function(yList, XbasisList, modelList, rhoMat,
convcrit=1e-4, iterlim=20, dbglev=1,
parMap=diag(rep(1,nparam))) {
# Data2LD.opt optimizes a parameter vector theta defining a
# linear differential operator object used to smooth a set of data.
#
# Arguments:
# YLIST an array containing values of curves
# If the array is a matrix, rows must correspond to argument
# values and columns to replications, and it will be assumed
# that there is only one variable per observation.
# If Y is a three-dimensional array, the first dimension
# corresponds to argument values, the second to replications,
# and the third to variables within replications.
# If Y is a vector, only one replicate and variable are
# assumed.
# BASISLIST A functional data object or a BASIS object. If so, the
# smoothing parameter LAMBDA is set to 0.
# MODELLIST A list of length NVAR. Each list contains a
# struct object with members:
# XList list of length number of homogeneous terms
# Each list contains a struct object with members:
# WfdPar a fdPar object for the coefficient
# variable the index of the variable
# derivative the order of its derivative
# npar if coefficient estimated, its location
# in the composite vector
# FList list of length number of forcing terms
# Each list contains a struct object with members:
# AfdPar an fdPar object for the coefficient
# Ufd an fd object for the forcing function
# npar if coefficient estimated, its location
# in the composite vector
# order the highest order of derivative
# name a tag for the variable
# nallXterm the number of homogeneous terms
# nallFterm the number of forcing functions
# RHOMAT A matrix with the number of rows equal to the number of values of the
# smoothing parameter per variable to be used, and number of columns
# equal to the number of variables.
# Each entry in the matrix is a value \rho in the interval [0,1)
# so that 0 <= \rho < 1.
# For each variable ivar in the system of equations and each optimization
# iopt, the data are weighted by \rho(iopt,ivar) and the
# roughness penalty by 1-\rho(iopt,ivar).
# It is expected that the values of \rho within each column will be
# in ascending order, and the estimated parameters for each row
# are passed along as initial values to be used for the optimization
# defined by the values of \rho in the next row.
# CONVCRIT One convergence criterion, or a vector of two criteria.
# The first criterion is applied to the function change,
# and the second is applied to the gradient norm.
# ITERLIM Maximum number of iterations allowed.
# DBGLEV An integer controlling amount of output per iteration.
# Defaults to 1, which prints summary results for
# each iteration.
# ACTIVE A vector of indices of the parameters to be optimized. If the index
# of a parameter does not appear in this vector, it is left fixed
# at its initial value even if its specification in coefList is that
# it be estimated. If the parameters are constrained by K linear
# constraints, ACTIVE contains integers from 1 to NPAR - K, where
# NPAR is the number of parameters in the model. In other words,
# ACTIVE specifies the constrained parameters to be estimated.
# PARMAP A rectangular matrix with number of rows equal to the
# number of parameters to be estimated as defined in
# BWTLIST, and number of columns equal to the number of
# parameters less the number of linear constraints on the
# estimated parameters. The columns of PARMAP must be
# orthnormal so that t(PARMAP) %*% PARMAP is an identity matrix.
# t(PARMAP) %*% THETA maps unconstrained parameters and the
# corresponding gradient into constrained parameter space.
# PARMAP %*% THETA maps constrained parameters and the
# corresponding gradient into unconstrained parameter space.
# PARMAP will usually be set up using the full QR
# decomposition of a linear constraint coefficient matrix t(A)
# where the constraints are of the form A P <- B, A and B
# being known matrices. An example of such a constraint
# that arises often is one where two estimated coefficients
# are constrained to be equal. For example, if a vvariable
# X involved in an equation the form a(x - x.0), where x.0
# is a fixed set point or defined target level for variable
# X, then this would be set up as a.1 x + a.2 x.0, where
# coefficients a.1 and a.2 are constrained to be equal in
# magnitude but opposite in sign, or a.1 + a.2 <- 0.
# Returns: A named list object with these fields:
# THETA.OPT The optimal parameter values.
# BWTLIST.OPT The corresponding optimal coefficients for the
# homogeneous terms in the differential equation
# AWTLIST.OPT The corresponding optimal coefficients for the
# forcing terms in the differential equation
# Last modified 23 April 2020
# ------------------------------------------------------------------------
# Check input parameters
# ------------------------------------------------------------------------
nvar <- length(yList)
nparam <- 0
for (ivar in 1:nvar) {
modelListi <- modelList[[ivar]]
if (modelListi$nallXterm > 0) {
for (iw in 1:modelListi$nallXterm) {
modelListiw <- modelListi$XList[[iw]]
nparam <- nparam + length(modelListiw$index)
}
}
if (modelListi$nallFterm > 0) {
for (jforce in 1:modelListi$nallFterm) {
modelListj <- modelListi$FList[[jforce]]
nparam = nparam + length(modelListj$index)
}
}
}
# Get the vector of parameters
theta <- modelList2Vec(modelList, nparam)
ntheta <- length(theta)
thetaCon <- t(parMap) %*% theta
nparCon <- length(thetaCon)
climit <- matrix(c(-1000,1000),2,1) %*% matrix(1,1,nparCon)
active <- 1:nparCon
# check rhoMat and get nopt
if (!is.matrix(rhoMat)) {
rhoMat = matrix(rhoMat,length(rhoMat),1)
}
rhodim <- dim(rhoMat)
if (rhodim[1] != nvar) {
stop(paste("The first dimension of RHOMAT is not equal to",
"the number of variables"))
}
nopt <- rhodim[2]
# Lists to contain results if more than a single set of rho's are involved
thetastore <- matrix(0,ntheta,nopt)
dfstore <- matrix(0,nopt,1)
gcvstore <- matrix(0,nopt,1)
modelList.optList <- vector("list", nopt)
modelListnew <- modelList
# ------------------------------------------------------------------------
# loop through rho vectors
# ------------------------------------------------------------------------
for (iopt in 1:nopt) {
rhoVeci <- rhoMat[,iopt]
# compute initial criterion value, gradient and hessian
Data2LDList <- Data2LD(yList, XbasisList, modelListnew, rhoVeci)
fvec <- Data2LDList$MSE
grad <- Data2LDList$DpMSE
gradCon <- t(parMap) %*% grad
fvecsum <- sum(fvec)
norm <- sqrt(mean(gradCon^2))
# evaluate the initial update vector for correcting the initial theta
deltac <- -gradCon
# initialize iteration status arrays
iternum <- 0
status <- c(iternum, fvec, norm)
if (dbglev >= 1) {
cat("\nIter. Criterion Grad Length")
cat("\n")
# # cat(format(c(iternum,status[2:3]),digits=4))
# cat(iternum)
# cat(" ")
# cat(format(status[2],nsmall=4,digits=4))
# cat(" ")
# cat(format(status[3],nsmall=4,digits=4))
# cat("\n")
cat(iternum)
cat(" ")
cat(round(status[2],6))
cat(" ")
cat(round(status[3],6))
} else {
cat(".")
}
cat("\n")
iterhist <- matrix(0,iterlim+1,length(status))
iterhist[1,] <- status
if (iterlim == 0) {
return ()
}
# ------- Begin main iterations -----------
MAXSTEPITER <- 5
MAXSTEP <- 1000
trial <- 0.1
reset <- 0
linemat <- matrix(0,3,5)
thetaoldCon <- thetaCon
fvecsumold <- fvecsum
gradoldCon <- gradCon
dbgwrd <- dbglev >= 2
# --------------- beginning of optimization loop -----------
for (iter in 1:iterlim) {
iternum <- iternum + 1
# set logical parameters
dblwrd <- c(FALSE,FALSE)
limwrd <- c(FALSE,FALSE)
stpwrd <- FALSE
ind <- 0
ips <- 0
# compute slope
linemat[2,1] <- sum(deltac*gradoldCon)
# normalize search direction vector
sdg <- sqrt(sum(deltac^2))
deltac <- deltac/sdg
linemat[2,1] <- linemat[2,1]/sdg
# initialize line search vectors
linemat[,1:4] <- matrix(c(0, linemat[2,1], fvecsum),3,1) %*% matrix(1,1,4)
stepiter <- 0
if (dbglev >= 2) {
cat(paste(" ", stepiter, " "))
cat(format(round(t(linemat[,1]),6)))
cat("\n")
}
# return with error condition if initial slope is nonnegative
if (linemat[2,1] >= 0) {
if (dbglev >= 2) {
cat("Initial slope ")
cat(format(round(linemat[2,1])))
cat(" is nonnegative.\n")
cat("\n")
}
deltac <- -gradCon
linemat[2,1] <- -sum(gradCon^2)
sdg <- sqrt(sum(deltac^2))
deltac <- deltac/sdg
linemat[,1:4] <- matrix(c(0, linemat[2,1], fvecsum),3,1) %*% matrix(1,1,4)
linestore <- matrix(linemat[,1],3,1)
if (dbglev > 2) {
cat(paste(" ", stepiter, " "))
cat(format(round(t(linemat[,1]),6)))
cat("\n")
}
}
# return successfully if initial slope is very small
if (linemat[2,1] >= -min(c(1e-3,convcrit))) {
if (dbglev >= 2) {
cat("Initial slope too small.\n")
status <- c(iternum, fvecsum, norm)
cat("\n")
}
modelListnew <- modelList
break
}
# first step set to trial
if (is.nan(trial)) stop("is.nan(trial)")
linemat[1,5] <- trial
# ------------ begin line search iteration loop ----------
thetanewCon <- thetaCon
for (stepiter in 1:MAXSTEPITER) {
# check the step size and modify if limits exceeded
stepval <- linemat[1,5]
result <- stepchkData2LD(linemat[1,5], thetaCon, deltac, limwrd, ind,
climit, active, dbgwrd)
linemat[1,5] <- result[[1]]
ind <- result[[2]]
limwrd <- result[[3]]
# break if limit hit twice in a row
if (ind == 1) {
thetanew <- theta
fvecsumnew <- fvecsum
break
}
# break if current step size too small
if (linemat[1,5] <= 1e-7) {
if (dbglev >= 2) {
cat("\nStepsize too small: ")
cat(linemat[1,5])
cat("\n")
}
thetanew <- theta
fvecsumnew <- fvecsum
break
}
# update parameter vector
thetanewCon <- thetaCon + linemat[1,5]*deltac
# --------- update function, gradient and hessian -----------
thetanew <- parMap %*% thetanewCon
modelListnew <- modelVec2List(modelList, thetanew)
Data2LDList <- Data2LD(yList, XbasisList, modelListnew,
rhoVeci, summary=FALSE)
fvecnew <- Data2LDList$MSE
gradnew <- Data2LDList$DpMSE
hessmatnew <- Data2LDList$D2ppMSE
gradnewCon <- t(parMap) %*% gradnew
hessmatnewCon <- t(parMap) %*% hessmatnew %*% parMap
fvecsumnew <- sum(fvecnew)
# -------------------------------------------------------------
linemat[3,5] <- fvecsumnew
# compute new directional derivative
linemat[2,5] <- sum(deltac*gradnewCon)
if (dbglev >= 2) {
cat(paste(" ", stepiter, " "))
cat(format(round(t(linemat[,5]),6)))
cat("\n")
}
# compute next line search step, also testing for convergence
result <- stepit(linemat, ips, dblwrd, MAXSTEP)
linemat <- result[[1]]
ips <- result[[2]]
ind <- result[[3]]
dblwrd <- result[[4]]
trial <- linemat[1,5]
# ind == 0 implies convergence
if (ind == 0 || ind == 5) {
break
}
}
# ------------ end line search iteration loop ----------
thetaCon <- thetanewCon
theta <- thetanew
fvecsum <- fvecsumnew
gradCon <- gradnewCon
# test for function value made worse
if (fvecsum > fvecsumold) {
# Function value worse warn and terminate
if (dbglev >= 2) {
cat("\n")
cat("Criterion increased: ")
cat(format(round(c(fvecsumold, fvecsum),4)))
cat("\n")
}
# reset parameters and fit
thetaCon <- thetaoldCon
fvecsum <- fvecsumold
deltac <- gradCon
if (reset == 1) {
# This is the second time in a row that this
# has happened quit
if (dbglev >= 2) {
cat("Reset twice, terminating.\n")
cat("Convergence not attained.\n")
}
# return current status of optimization
break
} else {
reset <- 1
}
} else {
# function value has not increased, check for convergence
RMSgrad <- sqrt(mean(gradCon^2))
if (length(convcrit) > 1) {
convtest <- fvecsumold - fvecsum < convcrit[1] && RMSgrad < convcrit[2]
} else {
convtest <- fvecsumold - fvecsum < convcrit
}
if (convtest) {
norm <- sqrt(mean(gradCon^2))
status <- c(iternum, fvecsum, norm)
if (dbglev >= 1) {
cat(iternum)
cat(" ")
cat(round(status[2],6))
cat(" ")
cat(round(status[3],6))
cat("\n")
cat("Convergence reached.\n")
}
# return current status of optimization
break
}
# update old parameter vectors and fit structure
thetaoldCon <- thetaCon
fvecsumold <- fvecsum
gradoldCon <- gradCon
hessmatCon <- t(parMap) %*% hessmatnew %*% parMap
# update the line search direction vector
deltac <- -solve(hessmatCon,gradCon)
reset <- 0
}
norm <- sqrt(mean(gradCon^2))
status <- c(iternum, fvecsum, norm)
iterhist[iter+1,] <- status
if (dbglev >= 1) {
cat(iternum)
cat(" ")
cat(round(status[2],6))
cat(" ")
cat(round(status[3],6))
cat("\n")
}
if (dbglev >= 1 && iter == iterlim) {
cat(
"\nMaximum iterations reached but convergence not attained.\n")
}
}
# --------------- end of optimization loop -----------
cat("\n")
# evaluate the solution at the final optimal value of theta
Data2LDList <- Data2LD(yList, XbasisList, modelListnew, rhoVeci, summary=TRUE)
# store the optimal values for this level of the \rho's the
# storage objects
theta <- modelList2Vec(modelListnew, nparam)
thetastore[,iopt] <- theta
dfstore[iopt] <- Data2LDList$df
gcvstore[iopt] <- Data2LDList$gcv
modelList.optList[[iopt]] <- modelListnew
}
# ------------------- end of rho loop ------------------
return(list(thetastore=thetastore, dfstore=dfstore, gcvstore=gcvstore,
modelList.optList=modelList.optList))
}
# ----------------------------------------------------------------------------
stepchkData2LD = function(stepval, cvec, deltac, limwrd, ind,
climit=50*matrix(rep(1,2),2,1) %*% matrix(1,1,npar),
active=1:npar, dbgwrd=FALSE) {
#STEPCHK checks the step along a line for producing parameters within the
# limits specified by BOT and TOP
# LIMWRD Logical variable permitting detection that parameter
# was on the boundary two steps in a row
# Last modified 28 August 2017
npar = length(deltac)
bot = climit[1,]
top = climit[2,]
# if step is too small, return with flag ind = 1
if (stepval < 1e-7) {
ind = 1
return(list(stepval=stepval,ind=ind,limwrd=limwrd))
}
# ensure that step does not go beyond lower limit on parameters
stepvali = stepval*deltac
if (any(stepvali[active] < bot[active]-cvec[active])) {
index = active[stepvali[active] < bot[active]-cvec[active] &
deltac[active] != 0]
stepvalnew = min((bot[index]-cvec[index])/deltac[index])
stepval = stepvalnew
# check whether lower limit has been reached twice in a row
if (limwrd[1]) {
ind = 1
return(list(stepval=stepval,ind=ind,limwrd=limwrd))
} else {
limwrd[1] = TRUE
}
} else {
limwrd[1] = FALSE
}
if (stepval < 1e-7) {
ind = 1
return(list(stepval=stepval,ind=ind,limwrd=limwrd))
}
# ensure that step does not go beyond upper limit on parameters
stepvali = stepval*deltac
if (any(stepvali[active] > top[active]-cvec[active])) {
index = active[stepvali[active] > top[active]-cvec[active] &
deltac[active] != 0]
stepvalnew = min((top[index]-cvec[index])/deltac[index])
stepval = stepvalnew
# check whether upper limit has been reached twice in a row
if (limwrd[2]) {
ind = 1
return(list(stepval=stepval,ind=ind,limwrd=limwrd))
} else {
limwrd[2] = TRUE
}
} else {
limwrd[2] = FALSE
}
return(list(stepval=stepval,ind=ind,limwrd=limwrd))
}
# ----------------------------------------------------------------------------
stepit <- function(linemat, ips, dblwrd, MAXSTEP) {
#STEPIT computes next step size in line search algorithm
# Arguments:
# LINEMAT: Row 1 contains step values
# Row 2 contains slope values
# Row 3 contains function values
# IPS: If 1, previous slope was positive
# DBLWRD: Vector of length 2: dblwrd[1] T means step halved
# dblwrd[2] T means step doubled
# MAXSTEP: maximum size of step
# Last modified 4 March 2014 by Jim Ramsay
# Wolfe condition 1
test1.1 = linemat[3,5] <= linemat[3,1] + linemat[1,5]*linemat[2,1]/20
# Wolfe condition 2
test1.2 = abs(linemat[2,5]) <= abs(linemat[2,1])/10
# disp([test1.1, test1.2])
test1 = test1.1 && test1.2
# test1 = test1.2
test2 = linemat[3,5] > linemat[3,1]
test3 = linemat[2,5] > 0
if ((test1 || !test3) && test2) {
# ************************************************************
# function is worse and either slope is satisfory or negative
ips = 0 # step is halved
if (dblwrd[2]) {
ind = 5
return(list(linemat = linemat, ips = ips, ind = ind, dblwrd = dblwrd))
}
linemat[1,5] = min(c(linemat[1,5]/2, MAXSTEP))
linemat[,2] = linemat[,1]
linemat[,3] = linemat[,1]
dblwrd = c(1, 0)
ind = 2
return(list(linemat = linemat, ips = ips, ind = ind, dblwrd = dblwrd))
}
# *********************************************************
if (test1) {
# test1 means successful convergence
ind = 0
return(list(linemat = linemat, ips = ips, ind = ind, dblwrd = dblwrd))
}
# **********************************************************
if (test3) {
# Current slope is positive
ips = 1
linemat[,4] = linemat[,5]
deltaf = linemat[3,3] - linemat[3,5]
z = (3/(linemat[1,5] - linemat[1,3]))*deltaf + linemat[2,3] + linemat[2,5]
w = z * z - linemat[2,3] * linemat[2,5]
if (abs(linemat[2,3] + linemat[2,5] + 2 * z) >= 1e-05 && w > 0) {
if (w < 0) stop("w < 0")
w = sqrt(w)
denom = linemat[2,5] - linemat[2,3] + 2 * w
if (denom == 0) stop("denom == 0")
linemat[1,5] = linemat[1,3] + (1 - ((linemat[2,5] + w - z)/
(denom))) * (linemat[1,5] - linemat[1,3])
if (is.nan(linemat[1,5])) stop("first test3: is.nan(linemat[1,5])")
} else {
# linear interpolation necessary
aerror = linemat[1,3]
if (linemat[1,5] > linemat[1,3]) {
aerror = linemat[1,5]
}
linemat[1,5] = linemat[1,3] - linemat[2,3] * ((linemat[1,5] - linemat[1,3])/
(linemat[2,5] - linemat[2,3]))
if (linemat[1,5] > 2 * aerror) {
linemat[1,5] = 2 * aerror
}
}
linemat[1,5] = min(c(linemat[1,5], MAXSTEP))
dblwrd = c(0,0)
ind = 2
return(list(linemat = linemat, ips = ips, ind = ind, dblwrd = dblwrd))
}
# *************************************************************
# Current slope is negative or zero
linemat[,2] = linemat[,3]
linemat[,3] = linemat[,5]
if (ips == 1) {
# *****************************************************
# previous slope is positive
deltaf = linemat[3,5] - linemat[3,4]
z = c(3/(linemat[1,4] - linemat[1,5])) * deltaf +
linemat[2,5] + linemat[2,4]
w = z * z - linemat[2,5] * linemat[2,4]
if (abs(linemat[2,5] + linemat[2,4] + 2 * z) >= 1e-05 && w > 0) {
w = sqrt(w)
linemat[1,5] = linemat[1,5] + (1 - ((linemat[2,4] + w - z)/
(linemat[2,4] - linemat[2,5] +
2 * w))) * (linemat[1,4] - linemat[1,5])
} else {
aerror = linemat[1,5]
if (linemat[1,4] > linemat[1,5]) {
aerror = linemat[1,4]
}
linemat[1,5] = linemat[1,5] - linemat[2,5] *
((linemat[1,4] - linemat[1,5])/
(linemat[2,4] - linemat[2,5]))
if (linemat[1,5] > 2 * aerror) {
linemat[1,5] = 2 * aerror
}
}
linemat[1,5] = min(c(linemat[1,5], MAXSTEP))
dblwrd = c(0,0)
ind = 2
return(list(linemat = linemat, ips = ips, ind = ind, dblwrd = dblwrd))
}
# ******************************************************
if ((linemat[2,3] - linemat[2,2]) * (linemat[1,3] - linemat[1,2]) > 0) {
# previous slope is negative
z = c(3/(linemat[1,3] - linemat[1,2])) * (linemat[3,2] - linemat[3,3]) +
linemat[2,2] + linemat[2,3]
w = z * z - linemat[2,2] * linemat[2,3]
if (abs(linemat[2,2] + linemat[2,3] + 2 * z) >= 1e-05 && w > 0) {
w = sqrt(w)
linemat[1,5] = linemat[1,2] +
(1 - ((linemat[2,3] + w - z)/(linemat[2,3] - linemat[2,2] +
2 * w))) * (linemat[1,3] - linemat[1,2])
} else {
linemat[1,5] = linemat[1,2] - linemat[2,2] *
((linemat[1,3] - linemat[1,2])/
(linemat[2,3] - linemat[2,2]))
}
linemat[1,5] = min(c(linemat[1,5], MAXSTEP))
dblwrd = c(0,0)
ind = 2
return(list(linemat = linemat, ips = ips, ind = ind, dblwrd = dblwrd))
} else {
# previous slope also negative but not as much
if (dblwrd[1]) {
ind = 5
return(
list(linemat = linemat, ips = ips, ind = ind, dblwrd = dblwrd))
} else {
linemat[1,5] = 2 * linemat[1,5]
linemat[1,5] = min(c(linemat[1,5], MAXSTEP))
dblwrd = c(0,1)
ind = 2
return(
list(linemat = linemat, ips = ips, ind = ind, dblwrd = dblwrd))
}
}
ind = 2
}
Try the Data2LD 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.