-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathnumDeriv.R
243 lines (209 loc) · 8.31 KB
/
numDeriv.R
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
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
# grad case 1 and 2 are special cases of jacobian, with a scalar rather than
# vector valued function. Case 3 differs only because of the interpretation
# that the vector result is a scalar function applied to each argument, and the
# thus the result has the same length as the argument.
# The code of grad could be consolidated to use jacobian.
# There is also some duplication in genD.
############################################################################
# functions for gradient calculation
############################################################################
grad <- function (func, x, method="Richardson", side=NULL,
method.args=list(), ...) UseMethod("grad")
grad.default <- function(func, x, method="Richardson", side=NULL,
method.args=list(), ...){
# modified by Paul Gilbert from code by Xingqiao Liu.
# case 1/ scalar arg, scalar result (case 2/ or 3/ code should work)
# case 2/ vector arg, scalar result (same as special case jacobian)
# case 3/ vector arg, vector result (of same length, really 1/ applied multiple times))
f <- func(x, ...)
n <- length(x) #number of variables in argument
if (is.null(side)) side <- rep(NA, n)
else {
if(n != length(side))
stop("Non-NULL argument 'side' should have the same length as x")
if(any(1 != abs(side[!is.na(side)])))
stop("Non-NULL argument 'side' should have values NA, +1, or -1.")
}
case1or3 <- n == length(f)
if((1 != length(f)) & !case1or3)
stop("grad assumes a scalar valued function.")
if(method=="simple"){
# very simple numerical approximation
args <- list(eps=1e-4) # default
args[names(method.args)] <- method.args
side[is.na(side)] <- 1
eps <- rep(args$eps, n) * side
if(case1or3) return((func(x+eps, ...)-f)/eps)
# now case 2
df <- rep(NA,n)
for (i in 1:n) {
dx <- x
dx[i] <- dx[i] + eps[i]
df[i] <- (func(dx, ...) - f)/eps[i]
}
return(df)
}
else if(method=="complex"){ # Complex step gradient
if (any(!is.na(side))) stop("method 'complex' does not support non-NULL argument 'side'.")
eps <- .Machine$double.eps
v <- try(func(x + eps * 1i, ...))
if(inherits(v, "try-error"))
stop("function does not accept complex argument as required by method 'complex'.")
if(!is.complex(v))
stop("function does not return a complex value as required by method 'complex'.")
if(case1or3) return(Im(v)/eps)
# now case 2
h0 <- rep(0, n)
g <- rep(NA, n)
for (i in 1:n) {
h0[i] <- eps * 1i
g[i] <- Im(func(x+h0, ...))/eps
h0[i] <- 0
}
return(g)
}
else if(method=="Richardson"){
args <- list(eps=1e-4, d=0.0001, zero.tol=sqrt(.Machine$double.eps/7e-7), r=4, v=2, show.details=FALSE) # default
args[names(method.args)] <- method.args
d <- args$d
r <- args$r
v <- args$v
show.details <- args$show.details
a <- matrix(NA, r, n)
#b <- matrix(NA, (r - 1), n)
# first order derivatives are stored in the matrix a[k,i],
# where the indexing variables k for rows(1 to r), i for columns (1 to n),
# r is the number of iterations, and n is the number of variables.
h <- abs(d*x) + args$eps * (abs(x) < args$zero.tol)
pna <- (side == 1) & !is.na(side) # double these on plus side
mna <- (side == -1) & !is.na(side) # double these on minus side
for(k in 1:r) { # successively reduce h
ph <- mh <- h
ph[pna] <- 2 * ph[pna]
ph[mna] <- 0
mh[mna] <- 2 * mh[mna]
mh[pna] <- 0
if(case1or3) a[k,] <- (func(x + ph, ...) - func(x - mh, ...))/(2*h)
else for(i in 1:n) {
if((k != 1) && (abs(a[(k-1),i]) < 1e-20)) a[k,i] <- 0 #some func are unstable near zero
else a[k,i] <- (func(x + ph*(i==seq(n)), ...) -
func(x - mh*(i==seq(n)), ...))/(2*h[i])
}
if (any(is.na(a[k,]))) stop("function returns NA at ", h," distance from x.")
h <- h/v # Reduced h by 1/v.
}
if(show.details) {
cat("\n","first order approximations", "\n")
print(a, 12)
}
#------------------------------------------------------------------------
# 1 Applying Richardson Extrapolation to improve the accuracy of
# the first and second order derivatives. The algorithm as follows:
#
# -- For each column of the derivative matrix a,
# say, A1, A2, ..., Ar, by Richardson Extrapolation, to calculate a
# new sequence of approximations B1, B2, ..., Br used the formula
#
# B(i) =( A(i+1)*4^m - A(i) ) / (4^m - 1) , i=1,2,...,r-m
#
# N.B. This formula assumes v=2.
#
# -- Initially m is taken as 1 and then the process is repeated
# restarting with the latest improved values and increasing the
# value of m by one each until m equals r-1
#
# 2 Display the improved derivatives for each
# m from 1 to r-1 if the argument show.details=T.
#
# 3 Return the final improved derivative vector.
#-------------------------------------------------------------------------
for(m in 1:(r - 1)) {
a <- (a[2:(r+1-m),,drop=FALSE]*(4^m)-a[1:(r-m),,drop=FALSE])/(4^m-1)
if(show.details & m!=(r-1) ) {
cat("\n","Richarson improvement group No. ", m, "\n")
print(a[1:(r-m),,drop=FALSE], 12)
}
}
return(c(a))
} else stop("indicated method ", method, "not supported.")
}
jacobian <- function (func, x, method="Richardson", side=NULL,
method.args=list(), ...) UseMethod("jacobian")
jacobian.default <- function(func, x, method="Richardson", side=NULL,
method.args=list(), ...){
f <- func(x, ...)
n <- length(x) #number of variables.
if (is.null(side)) side <- rep(NA, n)
else {
if(n != length(side))
stop("Non-NULL argument 'side' should have the same length as x")
if(any(1 != abs(side[!is.na(side)])))
stop("Non-NULL argument 'side' should have values NA, +1, or -1.")
}
if(method=="simple"){
# very simple numerical approximation
args <- list(eps=1e-4) # default
args[names(method.args)] <- method.args
side[is.na(side)] <- 1
eps <- rep(args$eps, n) * side
df <-matrix(NA, length(f), n)
for (i in 1:n) {
dx <- x
dx[i] <- dx[i] + eps[i]
df[,i] <- (func(dx, ...) - f)/eps[i]
}
return(df)
}
else if(method=="complex"){ # Complex step gradient
if (any(!is.na(side))) stop("method 'complex' does not support non-NULL argument 'side'.")
# Complex step Jacobian
eps <- .Machine$double.eps
h0 <- rep(0, n)
h0[1] <- eps * 1i
v <- try(func(x+h0, ...))
if(inherits(v, "try-error"))
stop("function does not accept complex argument as required by method 'complex'.")
if(!is.complex(v))
stop("function does not return a complex value as required by method 'complex'.")
h0[1] <- 0
jac <- matrix(NA, length(v), n)
jac[, 1] <- Im(v)/eps
if (n == 1) return(jac)
for (i in 2:n) {
h0[i] <- eps * 1i
jac[, i] <- Im(func(x+h0, ...))/eps
h0[i] <- 0
}
return(jac)
}
else if(method=="Richardson"){
args <- list(eps=1e-4, d=0.0001, zero.tol=sqrt(.Machine$double.eps/7e-7),
r=4, v=2, show.details=FALSE) # default
args[names(method.args)] <- method.args
d <- args$d
r <- args$r
v <- args$v
a <- array(NA, c(length(f),r, n) )
h <- abs(d*x) + args$eps * (abs(x) < args$zero.tol)
pna <- (side == 1) & !is.na(side) # double these on plus side
mna <- (side == -1) & !is.na(side) # double these on minus side
for(k in 1:r) { # successively reduce h
ph <- mh <- h
ph[pna] <- 2 * ph[pna]
ph[mna] <- 0
mh[mna] <- 2 * mh[mna]
mh[pna] <- 0
for(i in 1:n) {
a[,k,i] <- (func(x + ph*(i==seq(n)), ...) -
func(x - mh*(i==seq(n)), ...))/(2*h[i])
#if((k != 1)) a[,(abs(a[,(k-1),i]) < 1e-20)] <- 0 #some func are unstable near zero
}
h <- h/v # Reduced h by 1/v.
}
for(m in 1:(r - 1)) {
a <- (a[,2:(r+1-m),,drop=FALSE]*(4^m)-a[,1:(r-m),,drop=FALSE])/(4^m-1)
}
# drop second dim of a, which is now 1 (but not other dim's even if they are 1
return(array(a, dim(a)[c(1,3)]))
} else stop("indicated method ", method, "not supported.")
}