-
Notifications
You must be signed in to change notification settings - Fork 6
/
generateTPfullgrids.R
166 lines (129 loc) · 6.59 KB
/
generateTPfullgrids.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
#' Using a trained emulator, generate new full fields.
#'
#' This function takes in a trained emulator - a structure of class
#' \code{fldgen}.
#' This structure contains everything the emulator has learned about the model,
#' and is used to generate new fields of residuals. Also taking in a list of
#' generated residual fields and a global average yield, a global gridded
#' mean field is constructed accoring to the input reconstruction_function.
#' The mean field and residual fields are added to return a list of different
#' realizations of full fields.
#'
#'
#' @param emulator A trained \code{fldgen} temperature precipitation joint
#' emulator.
#' @param residgrids A list of new residual fields, each entry in the list is a
#' new realization, a matrix that is [Nyears x 2 * Ngrid]; the first 1:Ngrid
#' cols are the temperature residuals and columns (Ngrid + 1):(2*Ngrid) are the
#' precipitation residuals.
#' @param tgav A vector (or N x 1 matrix) of global mean temperatures, used to
#' calculate the mean warming response field.
#' @param tvarunconvert_fcn The function to undo any transformation done to the
#' input training data in \code{trainTP} to correct the support. This should be
#' the inverse function of the tvarconvert_fcn argument to \code{trainTP}. This
#' is stored in a \code{trainTP} \code{emulator} under
#' \code{emulator$griddattaT$tvarconvert_fcn}.
#' Defaults to NULL as temperature typically doesn't need to be transformed to
#' correct the support. WARNING: currently rely on the user to define the
#' correct inverse function here, though we do some checks and throw warnings
#' if it looks like there may be an issue.
#' @param pvarunconvert_fcn The function to undo any transformation done to the
#' input training data in \code{trainTP} to correct the support. This should be
#' the inverse function of the pvarconvert_fcn argument to \code{trainTP}. This
#' is stored in a \code{trainTP} \code{emulator} under
#' \code{emulator$griddattaP$pvarconvert_fcn}.
#' Defaults to exp() as precipitation is usually log-transformed in ordere to
#' correct the support. WARNING: currently rely on the user to define the
#' correct inverse function here, though we do some checks and throw warnings
#' if it looks like there may be an issue.
#' @param reconstruction_function A function for constructing a mean field from
#' trained pattern scaling result + a vector of global annual mean temperatures.
#' @return A list of:
#' 1) fullgrids = A list of new full fields, each entry in
#' the list is a new realization, a matrix that is [Nyears x 2 * Ngrid]; the
#' first 1:Ngrid cols are the temperature field and columns
#' (Ngrid + 1):(2*Ngrid) are the precipitation field.
#' 2) meanfieldT = the reconstructed, pattern scaled temperature mean field.
#' 3) meanfieldP = the reconstructed, pattern scaled precipitation mean field.
#' @export
generate.TP.fullgrids <- function(emulator, residgrids, tgav,
tvarunconvert_fcn = NULL, pvarunconvert_fcn = exp,
reconstruction_function = pscl_apply){
# dbl check that convert and unconvert functions are actual inverses on some test data.
# throw a warning if this is not the case but it is up to the user to track down.
if(!is.null(emulator$griddataT$tvarconvert_fcn)) {
if(is.null(tvarunconvert_fcn)) {
stop('Temperature transform function was applied, but no inverse function supplied.')
}
if(!chkinverse(emulator$griddataT$tvarconvert_fcn, tvarunconvert_fcn,
range=c(250,310))) {
warning('Your functions for transforming and inverse transforming the support of T may not be inverses of each other 1.')
}
}
if(!is.null(emulator$griddataP$pvarconvert_fcn)) {
if(is.null(pvarunconvert_fcn)) {
stop('Precipitation transform function was applied, but no inverse function supplied.')
}
if(!chkinverse(emulator$griddataP$pvarconvert_fcn, pvarunconvert_fcn,
range=c(1e-9, 1e-5), logspace=TRUE)) {
warning('Your functions for transforming and inverse transforming the support of P may not be inverses of each other 1.')
}
}
## Convert the tgav input to a column vector if necessary
assertthat::assert_that(is.numeric(tgav))
if(!is.matrix(tgav) || ncol(tgav) != 1) {
tgav <- matrix(tgav, ncol=1)
}
# save a copy of the number of grids cells for each variable
Ngrid <- ncol(residgrids[[1]])/2
meanfieldT <- reconstruction_function(emulator$meanfldT, tgav)
meanfieldP <- reconstruction_function(emulator$meanfldP, tgav)
# Add the meanfield values to the residual grids
lapply(residgrids, function(matrix, gridcells = Ngrid){
# Separate the tas and pr data from one another.
tas <- matrix[ , 1:Ngrid]
pr <- matrix[ , (Ngrid + 1):(2 * Ngrid)]
# Add the meanfield to the data
tas[ , 1:Ngrid] <- tas[ , 1:Ngrid] + meanfieldT
pr[ , 1:Ngrid] <- pr[ , 1:Ngrid] + meanfieldP
# convert from (-inf, inf) support to natural support.
if( !is.null(emulator$griddataT$tvarconvert_fcn)){
tas <- tvarunconvert_fcn(tas)
} else{
message('Generated T full fields not being transformed to a different support. Up to user to know if this is desirable.')
}
if(!is.null(emulator$griddataP$pvarconvert_fcn)){
pr <- pvarunconvert_fcn(pr)
} else{
message('Generated P full fields not being transformed to a different support. Up to user to know if this is desirable.')
}
# Return output
return(list(tas = tas, pr = pr))
}) ->
fullgrids
return(list(fullgrids = fullgrids, meanfieldT = meanfieldT, meanfieldP =
meanfieldP))
}
#' Check if two functions are inverses of one another over a given range of inputs.
#'
#' This is a quick and dirty check; it will definitely fail if the input range
#' is not in the domain of both functions.
#'
#' @param f1 First function to test
#' @param f2 Second function to test
#' @param range Range of test values
#' @param ntest Number of values to test in the range
#' @param logspace Flag indicating whether to space test values logarithmically.
#' @return logical
#' @keywords internal
chkinverse <- function(f1, f2, range=c(1.0, 5.0), ntest=10, logspace=FALSE)
{
if(logspace)
range <- log(range)
xvals <- seq(range[1], range[2], length.out = ntest)
if(logspace)
xvals <- exp(xvals)
xt1 <- f1(f2(xvals))
xt2 <- f2(f1(xvals))
isTRUE(all.equal(xt1, xvals)) && isTRUE(all.equal(xt2, xvals))
}