sensitivitiesSymb: Compute sensitivity equations of a function symbolically
In dlill/cOde2ndsens: Automated C Code Generation for Use with the 'deSolve' and 'bvpSolve' Packages
Description
Usage
Arguments
Details
Value
Examples
Description
Compute sensitivity equations of a function symbolically
Usage
1
2
Arguments
f
named vector of type character, the functions
states
Character vector. Sensitivities are computed with respect to initial
values of these states
parameters
Character vector. Sensitivities are computed with respect to initial
values of these parameters
inputs
Character vector. Input functions or forcings. They are excluded from
the computation of sensitivities.
reduce
Logical. Attempts to determine vanishing sensitivities, removes their
equations and replaces their right-hand side occurences by 0.
secondSens
Logical. Determines whether the second sensitivities are to be computed.
Details
The sensitivity equations are ODEs that are derived from the original ODE f.
They describe the sensitivity of the solution curve with respect to parameters like
initial values and other parameters contained in f. These equtions are also useful
for parameter estimation by the maximum-likelihood method. For consistency with the
time-continuous setting provided by adjointSymb, the returned equations contain
attributes for the chisquare functional and its gradient.
Value
Named vector of type character with the sensitivity equations. Furthermore,
attributes "chi" (the integrand of the chisquare functional), "grad" (the integrand
of the gradient of the chisquare functional), "forcings" (Character vector of the
additional forcings being necessare to compute chi and grad) and "yini" (
The initial values of the sensitivity equations) are returned.
Examples
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
## Not run:
##############################################################################################
## Sensitivity analysis of ozone formation
##############################################################################################
library(deSolve)
# O2 + O <-> O3
f <- c(
O3 = " build_O3 * O2 * O - decay_O3 * O3",
O2 = "-build_O3 * O2 * O + decay_O3 * O3",
O = "-build_O3 * O2 * O + decay_O3 * O3"
)
# Compute sensitivity equations
f_s <- sensitivitiesSymb(f)
# Generate ODE function
func <- funC(c(f, f_s))
# Initialize times, states, parameters and forcings
times <- seq(0, 15, by = .1)
yini <- c(O3 = 0, O2 = 3, O = 2, attr(f_s, "yini"))
pars <- c(build_O3 = .1, decay_O3 = .01)
# Solve ODE
out <- odeC(y = yini, times = times, func = func, parms = pars)
# Plot solution
par(mfcol=c(2,3))
t <- out[,1]
M1 <- out[,2:4]
M2 <- out[,5:7]
M3 <- out[,8:10]
M4 <- out[,11:13]
M5 <- out[,14:16]
M6 <- out[,17:19]
matplot(t, M1, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="solution")
legend("topright", legend = c("O3", "O2", "O"), lty=1, col=1:3)
matplot(t, M2, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="d/(d O3)")
matplot(t, M3, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="d/(d O2)")
matplot(t, M4, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="d/(d O)")
matplot(t, M5, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="d/(d build_O3)")
matplot(t, M6, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="d/(d decay_O3)")
## End(Not run)
## Not run:
##############################################################################################
## Estimate parameter values from experimental data
##############################################################################################
library(deSolve)
# O2 + O <-> O3
# diff = O2 - O3
# build_O3 = const.
f <- c(
O3 = " build_O3 * O2 * O - decay_O3 * O3",
O2 = "-build_O3 * O2 * O + decay_O3 * O3",
O = "-build_O3 * O2 * O + decay_O3 * O3"
)
# Compute sensitivity equations and get attributes
f_s <- sensitivitiesSymb(f)
chi <- attr(f_s, "chi")
grad <- attr(f_s, "grad")
forcings <- attr(f_s, "forcings")
# Generate ODE function
func <- funC(f = c(f, f_s, chi, grad), forcings = forcings, fcontrol = "nospline")
# Initialize times, states, parameters
times <- seq(0, 15, by = .1)
yini <- c(O3 = 0, O2 = 2, O = 2.5)
yini_s <- attr(f_s, "yini")
yini_chi <- c(chi = 0)
yini_grad <- rep(0, length(grad)); names(yini_grad) <- names(grad)
pars <- c(build_O3 = .2, decay_O3 = .1)
# Initialize forcings (the data)
data(oxygenData)
forcData <- data.frame(time = oxygenData[,1],
name = rep(colnames(oxygenData[,-1]), each=dim(oxygenData)[1]),
value = as.vector(oxygenData[,-1]))
forc <- setForcings(func, forcData)
# Solve ODE
out <- odeC(y = c(yini, yini_s, yini_chi, yini_grad),
times = times, func = func, parms = pars, forcings = forc,
method = "lsodes")
# Plot solution
par(mfcol=c(1,2))
t <- out[,1]
M1 <- out[,2:4]
M2 <- out[,names(grad)]
tD <- oxygenData[,1]
M1D <- oxygenData[,2:4]
matplot(t, M1, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="states")
matplot(tD, M1D, type="b", lty=2, col=1:3, pch=4, add=TRUE)
legend("topright", legend = names(f), lty=1, col=1:3)
matplot(t, M2, type="l", lty=1, col=1:5, xlab="time", ylab="value", main="gradient")
legend("topleft", legend = names(grad), lty=1, col=1:5)
# Define objective function
obj <- function(p) {
out <- odeC(y = c(p[names(f)], yini_s, yini_chi, yini_grad),
times = times, func = func, parms = p[names(pars)],
forcings = forc, method="lsodes")
value <- as.vector(tail(out, 1)[,"chi"])
gradient <- as.vector(tail(out, 1)[,paste("chi", names(p), sep=".")])
hessian <- gradient%*%t(gradient)
return(list(value = value, gradient = gradient, hessian = hessian))
}
# Fit the data
myfit <- optim(par = c(yini, pars),
fn = function(p) obj(p)$value,
gr = function(p) obj(p)$gradient,
method = "L-BFGS-B",
lower=0,
upper=5)
# Model prediction for fit parameters
prediction <- odeC(y = c(myfit$par[1:3], yini_s, yini_chi, yini_grad),
times = times, func = func, parms = myfit$par[4:5],
forcings = forc, method = "lsodes")
# Plot solution
par(mfcol=c(1,2))
t <- prediction[,1]
M1 <- prediction[,2:4]
M2 <- prediction[,names(grad)]
tD <- oxygenData[,1]
M1D <- oxygenData[,2:4]
matplot(t, M1, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="states")
matplot(tD, M1D, type="b", lty=2, col=1:3, pch=4, add=TRUE)
legend("topright", legend = names(f), lty=1, col=1:3)
matplot(t, M2, type="l", lty=1, col=1:5, xlab="time", ylab="value", main="gradient")
legend("topleft", legend = names(grad), lty=1, col=1:5)
## End(Not run)
dlill/cOde2ndsens documentation built on May 30, 2019, 1:37 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.
Description Usage Arguments Details Value Examples
Description
Compute sensitivity equations of a function symbolically
Usage
1 2 |
Arguments
f |
named vector of type character, the functions |
states |
Character vector. Sensitivities are computed with respect to initial values of these states |
parameters |
Character vector. Sensitivities are computed with respect to initial values of these parameters |
inputs |
Character vector. Input functions or forcings. They are excluded from the computation of sensitivities. |
reduce |
Logical. Attempts to determine vanishing sensitivities, removes their equations and replaces their right-hand side occurences by 0. |
secondSens |
Logical. Determines whether the second sensitivities are to be computed. |
Details
The sensitivity equations are ODEs that are derived from the original ODE f. They describe the sensitivity of the solution curve with respect to parameters like initial values and other parameters contained in f. These equtions are also useful for parameter estimation by the maximum-likelihood method. For consistency with the time-continuous setting provided by adjointSymb, the returned equations contain attributes for the chisquare functional and its gradient.
Value
Named vector of type character with the sensitivity equations. Furthermore,
attributes "chi" (the integrand of the chisquare functional), "grad" (the integrand
of the gradient of the chisquare functional), "forcings" (Character vector of the
additional forcings being necessare to compute chi and grad) and "yini" (
The initial values of the sensitivity equations) are returned.
Examples
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 | ## Not run:
##############################################################################################
## Sensitivity analysis of ozone formation
##############################################################################################
library(deSolve)
# O2 + O <-> O3
f <- c(
O3 = " build_O3 * O2 * O - decay_O3 * O3",
O2 = "-build_O3 * O2 * O + decay_O3 * O3",
O = "-build_O3 * O2 * O + decay_O3 * O3"
)
# Compute sensitivity equations
f_s <- sensitivitiesSymb(f)
# Generate ODE function
func <- funC(c(f, f_s))
# Initialize times, states, parameters and forcings
times <- seq(0, 15, by = .1)
yini <- c(O3 = 0, O2 = 3, O = 2, attr(f_s, "yini"))
pars <- c(build_O3 = .1, decay_O3 = .01)
# Solve ODE
out <- odeC(y = yini, times = times, func = func, parms = pars)
# Plot solution
par(mfcol=c(2,3))
t <- out[,1]
M1 <- out[,2:4]
M2 <- out[,5:7]
M3 <- out[,8:10]
M4 <- out[,11:13]
M5 <- out[,14:16]
M6 <- out[,17:19]
matplot(t, M1, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="solution")
legend("topright", legend = c("O3", "O2", "O"), lty=1, col=1:3)
matplot(t, M2, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="d/(d O3)")
matplot(t, M3, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="d/(d O2)")
matplot(t, M4, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="d/(d O)")
matplot(t, M5, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="d/(d build_O3)")
matplot(t, M6, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="d/(d decay_O3)")
## End(Not run)
## Not run:
##############################################################################################
## Estimate parameter values from experimental data
##############################################################################################
library(deSolve)
# O2 + O <-> O3
# diff = O2 - O3
# build_O3 = const.
f <- c(
O3 = " build_O3 * O2 * O - decay_O3 * O3",
O2 = "-build_O3 * O2 * O + decay_O3 * O3",
O = "-build_O3 * O2 * O + decay_O3 * O3"
)
# Compute sensitivity equations and get attributes
f_s <- sensitivitiesSymb(f)
chi <- attr(f_s, "chi")
grad <- attr(f_s, "grad")
forcings <- attr(f_s, "forcings")
# Generate ODE function
func <- funC(f = c(f, f_s, chi, grad), forcings = forcings, fcontrol = "nospline")
# Initialize times, states, parameters
times <- seq(0, 15, by = .1)
yini <- c(O3 = 0, O2 = 2, O = 2.5)
yini_s <- attr(f_s, "yini")
yini_chi <- c(chi = 0)
yini_grad <- rep(0, length(grad)); names(yini_grad) <- names(grad)
pars <- c(build_O3 = .2, decay_O3 = .1)
# Initialize forcings (the data)
data(oxygenData)
forcData <- data.frame(time = oxygenData[,1],
name = rep(colnames(oxygenData[,-1]), each=dim(oxygenData)[1]),
value = as.vector(oxygenData[,-1]))
forc <- setForcings(func, forcData)
# Solve ODE
out <- odeC(y = c(yini, yini_s, yini_chi, yini_grad),
times = times, func = func, parms = pars, forcings = forc,
method = "lsodes")
# Plot solution
par(mfcol=c(1,2))
t <- out[,1]
M1 <- out[,2:4]
M2 <- out[,names(grad)]
tD <- oxygenData[,1]
M1D <- oxygenData[,2:4]
matplot(t, M1, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="states")
matplot(tD, M1D, type="b", lty=2, col=1:3, pch=4, add=TRUE)
legend("topright", legend = names(f), lty=1, col=1:3)
matplot(t, M2, type="l", lty=1, col=1:5, xlab="time", ylab="value", main="gradient")
legend("topleft", legend = names(grad), lty=1, col=1:5)
# Define objective function
obj <- function(p) {
out <- odeC(y = c(p[names(f)], yini_s, yini_chi, yini_grad),
times = times, func = func, parms = p[names(pars)],
forcings = forc, method="lsodes")
value <- as.vector(tail(out, 1)[,"chi"])
gradient <- as.vector(tail(out, 1)[,paste("chi", names(p), sep=".")])
hessian <- gradient%*%t(gradient)
return(list(value = value, gradient = gradient, hessian = hessian))
}
# Fit the data
myfit <- optim(par = c(yini, pars),
fn = function(p) obj(p)$value,
gr = function(p) obj(p)$gradient,
method = "L-BFGS-B",
lower=0,
upper=5)
# Model prediction for fit parameters
prediction <- odeC(y = c(myfit$par[1:3], yini_s, yini_chi, yini_grad),
times = times, func = func, parms = myfit$par[4:5],
forcings = forc, method = "lsodes")
# Plot solution
par(mfcol=c(1,2))
t <- prediction[,1]
M1 <- prediction[,2:4]
M2 <- prediction[,names(grad)]
tD <- oxygenData[,1]
M1D <- oxygenData[,2:4]
matplot(t, M1, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="states")
matplot(tD, M1D, type="b", lty=2, col=1:3, pch=4, add=TRUE)
legend("topright", legend = names(f), lty=1, col=1:3)
matplot(t, M2, type="l", lty=1, col=1:5, xlab="time", ylab="value", main="gradient")
legend("topleft", legend = names(grad), lty=1, col=1:5)
## End(Not run)
|
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.