1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
|
bmstep <- function(par, srchdirn, lower=NULL, upper=NULL, bdmsk=NULL, trace=0) {
## Find maximum steplength from par along srchdirn given bounds and masks
# ?? do we want step length as multiple of srchdirn?, or actual max step length
# 20101031 -- issue of bounds not working correctly
# - -Inf seems to upset bounds
# - no upper bounds gives troubles (applies to Rcgmin too!)
#
# Input:
# par = a vector containing the starting point
# srchdirn = the direction of search (a vector)
# 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)
# trace = control of output: 0 for none (default), >0 for output
##
# Output:
# A double giving the maximum steplength. Not bigger than maxstep.
#
########## length of vectors #########
n<-length(par)
############# bounds and masks ################
# check if there are bounds
if(is.null(lower) || ! any(is.finite(lower))) nolower<-TRUE else nolower<-FALSE
if(is.null(upper) || ! any(is.finite(upper))) noupper<-TRUE else noupper<-FALSE
# Next line NOT same as in bmchk(). Leave out bdmsk.
if(nolower && noupper) bounds<-FALSE else bounds<-TRUE
if (is.null(bdmsk)) bdmsk<-rep(1,n) # make sure we have values
if (any(bdmsk==0)) {
if (trace > 2) cat("Masks present -- adjusting search direction.\n")
srchdirn[which(bdmsk==0)]<-0 # adjust search direction for masked elements
}
if(nolower) lower<-rep(-Inf,n)
if(noupper) upper<-rep(Inf,n)
######## find maximum step (may be Inf) #############
# distance to bounds
d2lo<-par-lower
d2up<-upper-par
if (trace>0) {
cat("Distances to bounds, lower then upper\n")
print(d2lo)
print(d2up)
}
sslo<-rep(0,n)
ssup<-sslo
# Now want to get ssup -- stepsize to upper bound along directions where srchdirn>0
# Hard way, by loop
for (i in 1:n) {
if (bdmsk[i]==1) { # free parameter
sdi<-srchdirn[i]
if (sdi>0) ssup[i]<-d2up[i]/sdi
if (sdi<0) sslo[i]<- -d2lo[i]/sdi
# sdi==0, no changes
}
}
# another approach 20111022
# suppressWarnings(ssup2<-(bdmsk*(d2up/srchdirn)))
# suppressWarnings(sslo2<-((-1)*bdmsk*(d2lo/srchdirn)))
# cat("sslo2 & ssup2\n")
# print(sslo2)
# print(ssup2)
# cat("after adjustment\n")
# ssup2[which(srchdirn<=0)]<-0
# sslo2[which(srchdirn>=0)]<-0
# print(sslo2)
# print(ssup2)
# ss<-c(sslo2, ssup2)
# ss<-ss[which(ss>0)]
# ms2<-min(ss)
# cat("ms2=",ms2,"\n")
if (trace>0) {
cat("steplengths, lower then upper\n")
print(sslo)
print(ssup)
}
sslo<-sslo[which(sslo>0)]
ssup<-ssup[which(ssup>0)]
if (trace>0) {
cat("steplengths, truncated, lower then upper\n")
if (length(sslo)>0) print(sslo) else cat("sslo NULL\n")
if (length(ssup)>0) print(ssup) else cat("ssup NULL\n")
}
if (is.null(sslo) && is.null(ssup)) {# Not needed, min will return Inf
maxstep<-Inf
} else {
maxstep<-min(sslo,ssup)
}
} ## end of bmstep.R
|
