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
|
### ============================================================================
### Check forcing function data set, event inputs and time-lag input
### ============================================================================
checkforcings <- function (forcings, times, dllname,
initforc, verbose, fcontrol = list()) {
## Check the names of the initialiser function
if (is.null(initforc))
stop(paste("initforc should be loaded if there are forcing functions ",initforc))
if (inherits (initforc, "CFunc")) {
ModelForc <- body(initforc)[[2]]
} else if (is.loaded(initforc, PACKAGE = dllname, type = "") ||
is.loaded(initforc, PACKAGE = dllname, type = "Fortran")) {
ModelForc <- getNativeSymbolInfo(initforc, PACKAGE = dllname)$address
} else
stop(paste("initforc should be loaded if there are forcing functions ",initforc))
## Check the type of the forcing function data series
if (is.data.frame(forcings)) forcings <- list(a=forcings)
if (! is.list(forcings)) forcings <- list(a=forcings)
nf <- length(forcings)
#1 check if each forcing function consists of a 2-columned matrix
for (i in 1:nf) {
if (ncol(forcings[[i]]) != 2)
stop("forcing function data sets should consist of two-column matrix")
}
## Check the control elements (see optim code)
con <- list(method="linear", rule = 2, f = 0, ties = "ordered")
nmsC <- names(con)
con[(namc <- names(fcontrol))] <- fcontrol
if (length(noNms <- namc[!namc %in% nmsC]) > 0)
warning("unknown names in fcontrol: ", paste(noNms, collapse = ", "))
method <- pmatch(con$method, c("linear", "constant"))
if (is.na(method))
stop("invalid interpolation method for forcing functions")
# 1 if linear, 2 if constant...
## Check the timespan of the forcing function data series
# time span of forcing function data sets should embrace simulation time...
# although extrapolation is allowed if con$rule = 2 (the default)
r_t <- range(times)
for (i in 1:nf) {
r_f <- range(forcings[[i]][,1]) # time range of this forcing function
if (r_f[1] > r_t[1]) {
if (con$rule == 2) {
mint <- c(r_t[1],forcings[[i]][1,2] )
forcings[[i]] <- rbind(mint,forcings[[i]])
if(verbose)
warning(paste("extrapolating forcing function data sets to first timepoint",i))
} else
stop(paste("extrapolating forcing function data sets to first timepoint",i))
}
nr <- nrow(forcings[[i]])
if (r_f[2] < r_t[2]) {
if (con$rule == 2) {
maxt <- c(r_t[2],forcings[[i]][nr,2] )
forcings[[i]] <- rbind(forcings[[i]],maxt)
if(verbose)
warning(paste("extrapolating forcing function data sets to last timepoint",i))
} else
stop(paste("extrapolating forcing function data sets to last timepoint",i))
}
}
## Check what needs to be done in case the time series is not "ordered"
if (!identical(con$ties, "ordered")) { # see approx code
for (i in 1:nf) {
x <- forcings[[i]][,1]
nx <- length(x)
if (length(ux <- unique(x)) < nx) { # there are non-unique values
y <- forcings[[i]][,2]
ties <- con$tiesn
if (missing(ties))
warning("collapsing to unique 'x' values")
y <- as.vector(tapply(y, x, ties))
x <- sort(ux)
forcings[[i]] <- cbind(x, y)
} else { # values are unique, but need sorting
y <- forcings[[i]][,2]
o <- order(x)
x <- x[o]
y <- y[o]
forcings[[i]] <- cbind(x,y)
}
} # i
}
## In case the interpolation is of type "constant" and f not equal to 0
## convert y-series, so that always the left value is taken
if (method == 2 & con$f != 0) {
for (i in 1:nf) {
y <- forcings[[i]][,2]
YY <- c(y,y[length(y)])[-1]
forcings[[i]][,2] <- (1-con$f)*y + con$f*YY
}
}
## all forcings in one vector; adding index to start/end
fmat <- tmat <- NULL
imat <- rep(1,nf+1)
for (i in 1:nf) {
# Karline: check for NA in forcing series and remove those
ii <- apply(forcings[[i]],1,function(x)any(is.na(x)))
if (sum(ii) > 0) forcings[[i]] <- forcings[[i]][!ii,]
tmat <- c(tmat, forcings[[i]][,1])
fmat <- c(fmat, forcings[[i]][,2])
imat[i+1]<-imat[i]+nrow(forcings[[i]])
}
storage.mode(tmat) <- storage.mode(fmat) <- "double"
storage.mode(imat) <- "integer"
# DIRTY trick not to inflate the number of arguments:
# add method (linear/constant) to imat
return(list(tmat = tmat, fmat = fmat, imat = c(imat, method),
ModelForc = ModelForc))
}
### ============================================================================
### Check timelags data set - also passes "dllname" now (not yet used)
### ============================================================================
checklags <- function (lags, dllname) {
if (!is.null(lags)) {
lags$islag = 1L
if (is.null(lags$mxhist))
lags$mxhist <- 1e4
if (lags$mxhist <1)
lags$mxhist <- 1e4
lags$mxhist<-as.integer(lags$mxhist)
if (is.null(lags$interpol)) # 1= hermitian, 2 = higher order interpolation
lags$interpol <- 1
lags$interpol<-as.integer(lags$interpol)
lags$isfun <- 0L
} else
lags$islag <- 0L
return(lags)
}
|
