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
|
#
# fields is a package for analysis of spatial data written for
# the R software environment.
# Copyright (C) 2024 Colorado School of Mines
# 1500 Illinois St., Golden, CO 80401
# Contact: Douglas Nychka, douglasnychka@gmail.com,
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with the R software environment if not, write to the Free Software
# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
# or see http://www.r-project.org/Licenses/GPL-2
##END HEADER
mKrigFastPredict <- function(object, gridList, ynew = NULL,
derivative = 0, Z = NULL, drop.Z = FALSE,
NNSize=5, setupObject= NULL,
giveWarnings=TRUE)
{
#NOTE: covariance model is specified by the arguments in object$args
#cov.args <- c( object$args, list(...) )
# For convenience the Z covariates are already assumed to be
# in the unrolled form. But this may be awkward if this
# function is called directly
# See the code in predictSurface.mKrig for details. E.g. unrollZGrid
if (derivative != 0) {
stop("Derivatives not supported with fast prediction method")
}
if( ncol(object$c.coef)>1 ){
stop("Replicated fields currently not supported for fast predict.")
}
names( gridList)<- c("x","y")
np<- NNSize
xObs<- object$x
nx<- length(gridList$x )
ny<- length(gridList$y )
if (!is.null(ynew)) {
coef.hold <- mKrig.coef(object, ynew,
collapseFixedEffect=TRUE)
c.coef <- coef.hold$c.coef
beta <- coef.hold$beta
}
else {
c.coef <- object$c.coef
beta <- object$beta
}
# fixed part of the model this a polynomial of degree m-1
# Tmatrix <- fields.mkpoly(xnew, m=object$m)
# only do this if nt>0, i.e. there is a fixed part.
#
if (!drop.Z & (object$nZ > 0) & (derivative >0) ) {
stop("derivative not supported with Z covariate included
use drop.Z = FALSE to omit Z ")
}
if( object$nt>0){
xnew<- make.surface.grid( gridList)
if (derivative == 0) {
if (drop.Z | object$nZ == 0) {
# just evaluate polynomial and not the Z covariate
temp1 <- fields.mkpoly(xnew, m = object$m) %*%
beta[object$ind.drift, ]
}
else {
if( nrow( xnew) != nrow(as.matrix(Z)) ){
stop(paste("number of rows of covariate Z",
nrow(as.matrix(Z)),
" is not the same as the number of locations",
nrow( xnew) )
)
}
temp0 <- cbind(fields.mkpoly(xnew, m = object$m),as.matrix(Z))
temp1 <- temp0 %*% beta
}
}
else {
temp1 <- fields.derivative.poly(xnew, m = object$m,
beta[object$ind.drift,])
}
}
# add nonparametric part. Covariance basis functions
# times coefficients.
# syntax is the name of the function and then a list with
# all the arguments. This allows for different covariance functions
# that have been passed as their name.
# enlarge the grid if needed so that obs have np grid point points on all margins.
if( (min(xObs[,1]) < gridList$x[1]) | (max(xObs[,1]) > gridList$x[nx] ) ) {
stop( "x obs locations can not be outside the grid ")
}
if( (min(xObs[,2]) < gridList$y[1]) | (max(xObs[,2]) > gridList$y[ny]) ){
stop( "y obs locations can not be outside the grid ")
}
# adjust grid if needed to include a margin of NNSize+1 grid points beyond xObs
# these are the slightly larger grids by adding margins.
# also create sparse matrices.
if( is.null(setupObject) ){
setupObject<- mKrigFastPredictSetup(object,
gridList = gridList,
NNSize = NNSize,
giveWarnings = giveWarnings)
}
gridListNew<- setupObject$marginInfo$gridListNew
nxNew<- length(gridListNew$x )
nyNew<- length(gridListNew$y )
# indX and indY are the subset of indices that match gridList
# ( gridListNew contains the grids in gridList)
indX<- setupObject$marginInfo$indX
indY<- setupObject$marginInfo$indY
if( ((indX[2]- indX[1] + 1)!= nx)| ((indY[2]- indY[1] + 1)!= ny)) {
cat(" indX, nx")
print( c(indX, nx) )
cat(" indY, ny")
print( c(indY, ny) )
stop("mismatch between subset of larger grid and gridList passed")
}
c.coefWghts<- colSums( diag.spam( c(object$c.coef) ) %*%
setupObject$offGridObject$B )
c.coefWghts<- matrix( c.coefWghts, nxNew, nyNew )
# fast multiplication of covariances on the grid with
# the coefficients on the grid
temp2<- stationary.image.cov( Y=c.coefWghts, cov.obj=setupObject$cov.obj)
# cut down the size of temp2 trimming off margins using to approximate
# exact covariance kernel.
temp2<- temp2[ indX[1]:indX[2], indY[1]:indY[2] ]
# add fixed part and spatial parts together and coerce to matrix
if( object$nt>0){
return( (matrix(temp1,nx,ny) + temp2) )
}
else{
# return only the spatial part because the fixed part is absent.
return( temp2)
}
}
|
