#### Despite being named "test_varfield", these tests actually test the whole
#### process, end to end.
context('Field generation')
### Legacy T only
### Run the analysis and generation. The results of these calculations will be
### tested in the test case below.
tann1 <- train(system.file('extdata/tann1.nc', package='fldgen'))
griddata <- tann1$griddata
tgav <- tann1$tgav
pscl <- tann1$pscl
reof <- tann1$reof
Fxmag <- tann1$fx$mag
## Get the phase and complex fft directly from reof$x. These are only needed
## for tests, not production calculations
Fx <- mvfft(reof$x)
Fxphase <- atan2(Im(Fx), Re(Fx))
tempgrids <- list() # Empty list to hold the temperature
# realizations
length(tempgrids) <- 4
## Run with the phases of the actual time series
meanfield <- pscl_apply(pscl, tgav)
tempgrids[[1]] <- reconst_fields(reof$rotation, mkcorrts(tann1, Fxphase) , meanfield)
## Run the rest with random phases
for(i in 2:4)
## If you wanted to use only 50 PCs, like in the previous demo, you could
## use reof$rotation[,1:50] and Fxmag[,1:50] in the next call.
tempgrids[[i]] <- reconst_fields(reof$rotation, mkcorrts(tann1), meanfield)
test_that('legacy T - Values produced by global mean calculation agree with the ones produced by CDO.',
{
## matrix multiply: ntime x ngrid * ngrid x 1 = ntime x 1
tg_ours <- griddata$tas %*% griddata$tgop
tgdiff <- tg_ours - tgav
## CDO is a little more rigorous in its calculation (it doesn't use the
## small angle approximation for delta-latitude), so allow discrepancies of
## up to 1e-4 times the average magnitude.
ftol <- 1.0e-4 * mean(abs(tgav))
maxdiff <- max(abs(tgdiff))
expect_lte(maxdiff, ftol)
})
test_that('legacy T - Applying the pattern scaling and adding back the residuals recovers the original temperatures.',
{
tscl <- pscl_apply(pscl, tgav) + pscl$r
## Allow for some roundoff error
ftol <- 1.0e-6 * mean(abs(tgav))
diff <- tscl - griddata$tas
maxdiff <- max(abs(diff))
expect_lte(maxdiff, ftol)
})
test_that('legacy T - All EOFs are normalized and orthogonal to one another.',
{
## check that all eofs are normalized and orthogonal to one another. We do
## this by multiplying the matrix of basis vectors by its transpose and
## comparing to the identity matrix.
tprod <- t(reof$rotation) %*% reof$rotation
maxdiff <- max(abs(tprod - diag(ncol(tprod))))
expect_lte(maxdiff, 1.0e-8)
})
test_that('legacy T - EOFs other than PC0 have zero global average.',
{
## check that eofs other than PC0 have zero global average. (Technically
## this is guaranteed by the test above, but this test serves as independent
## confirmation.)
maxtg <- max( t(griddata$tgop) %*% reof$rotation[,-1] )
expect_lte(maxtg, 1.0e-8)
})
test_that('legacy T - PSD returned from training matches manually calculated.',
{
expect_equal(Fxmag, abs(Fx))
})
test_that('legacy T - Residual field can be recovered from basis and rotated coordinates.',
{
## We will be a little lenient on this test because there is plenty of
## opportunity for roundoff errors to accumulate in this calculation.
resid.reconst <- reconst_fields(reof$rotation, reof$x)
diff <- resid.reconst - pscl$r
maxdiff <- max(abs(diff))
ftol <- 1.0e-4 * mean(abs(pscl$r))
expect_lte(maxdiff, ftol)
})
## We will run this function with two different test data sets below.
validate_mkcorrts <- function(testts=NULL, method=1)
{
## if we weren't given any test time series, generate a pair of them
if(is.null(testts)) {
x <- 0:31
testts <- matrix(c(sin(4*pi*x/32.0), sin(5*pi*x/32.0)), ncol=2)
}
Fx <- mvfft(testts)
Fxmag <- abs(Fx)
Fxphase <- atan2(Im(Fx), Re(Fx))
phcoef <- phase_eqn_coef(Fx)
tsobj <- fldgen_object(NULL, NULL, NULL, NULL, Fxmag, phcoef, 'none')
## Running mkcorrts with these phases should reproduce the original time
## series
testout1 <- mkcorrts(tsobj, Fxphase, complexout=TRUE)
## Verify that there are no imaginary components
expect_equal(testout1, Re(testout1)*1+0i)
## Verify equality with original inputs
expect_equal(testts, Re(testout1))
## Running mkcorrts with random phases should produce statistically
## equivalent time series. It's not obvious how to test this in a
## non-tautological way. For now, we'll be satisfied with checking the mean
## and standard deviation. Note that we have to look at the absolute value
## of the mean because we can't guarantee the phase of the DC component is
## the same as the input.
testout2 <- mkcorrts(tsobj)
meansin <- abs(apply(testts, 2, mean))
meansout <- abs(apply(testout2, 2, mean))
expect_equal(meansin, meansout)
sdin <- apply(testts, 2, sd)
sdout <- apply(testout2, 2, sd)
expect_equal(sdin, sdout)
}
test_that('legacy T - Time series generated by mkcorrts method 1 have the correct properties.',
{
validate_mkcorrts()
validate_mkcorrts(reof$x)
})
test_that('legacy T - Time series generated by mkcorrts method 2 have the correct properties.',
{
validate_mkcorrts(method=2)
validate_mkcorrts(reof$x, method=2)
})
test_that('legacy T - psdest produces equivalent results to manual calculation.',
{
## do analysis with psdest
psd1 <- psdest(reof)
## repeat it manually
Fx <- mvfft(reof$x)
Fxmag <- abs(Fx)
expect_equal(sqrt(psd1$psd), Fxmag)
## A list of two identical inputs should also give back the same result.
psd2 <- psdest(list(reof, reof))
expect_equal(sqrt(psd2$psd), Fxmag)
})
###############################
### Joint TP
### Run the analysis and generation. The results of these calculations will be
### tested in the test case below.
emulator <- trainTP(c(system.file('extdata/tas_annual_esm_rcp_r2i1p1_2006-2100.nc', package='fldgen'),
system.file('extdata/pr_annual_esm_rcp_r2i1p1_2006-2100.nc', package='fldgen')),
tvarname = "tas", tlatvar='lat_2', tlonvar='lon_2',
tvarconvert_fcn = NULL,
pvarname = "pr", platvar='lat', plonvar='lon',
pvarconvert_fcn = log)
griddataT <- emulator$griddataT
griddataP <- emulator$griddataP
tgav <- emulator$tgav
meanfldT <- emulator$meanfldT
meanfldP <- emulator$meanfldP
reof <- emulator$reof
Fxmag <- emulator$fx$mag
## Get the phase and complex fft directly from reof$x. These are only needed
## for tests, not production calculations
Fx <- mvfft(reof$x)
Fxphase <- atan2(Im(Fx), Re(Fx))
## Run with one realization using the phases of the actual time series
## and the rest with random phases.
## First, the actual phases.
## Note that every entry of newgrids is in the normal space and needs to be
## converted back to the native space.
newgrids1 <- reconst_fields(reof$rotation, mkcorrts(emulator, phase = Fxphase))
Ngrid <- ncol(meanfldT$r)
residgrids1T <- unnormalize.resids(empiricalquant = emulator$tfuns$quant,
rn = newgrids1[ ,1:Ngrid])$rnew
residgrids1P <- unnormalize.resids(empiricalquant = emulator$pfuns$quant,
rn = newgrids1[ , (Ngrid+1):(2*Ngrid)])$rnew
## Run the rest with random phases and add
tmp <- generate.TP.resids(emulator, ngen = 3)
residgrids <- list()
length(residgrids) <- length(tmp) + 1
residgrids[[1]] <- cbind(residgrids1T, residgrids1P)
residgrids[[2]] <- tmp[[1]]
residgrids[[3]] <- tmp[[2]]
residgrids[[4]] <- tmp[[3]]
## Define the unconvert functions
pvarunconvert_fcn <- exp
tvarunconvert_fcn <- NULL
## use the new residuals in the native space with the mean field to reconstruct
## the full new fields
fullgrids <- generate.TP.fullgrids(emulator, residgrids,
tgav = tgav,
tvarunconvert_fcn = NULL, pvarunconvert_fcn = pvarunconvert_fcn,
reconstruction_function = pscl_apply)
fullgrids <- lapply(fullgrids$fullgrids,
function(g){
return(cbind(g[[1]], g[[2]]))
})
test_that('T Values produced by global mean calculation agree with the ones produced by CDO.',
{
## matrix multiply: ntime x ngrid * ngrid x 1 = ntime x 1
tg_ours <- griddataT$vardata %*% griddataT$globalop
tgdiff <- tg_ours - tgav
## CDO is a little more rigorous in its calculation (it doesn't use the
## small angle approximation for delta-latitude), so allow discrepancies of
## up to 1e-4 times the average magnitude.
ftol <- 1.0e-4 * mean(abs(tgav))
maxdiff <- max(abs(tgdiff))
expect_lte(maxdiff, ftol)
})
test_that('Applying the pattern scaling and adding back the residuals recovers the original temperatures.',
{
tscl <- pscl_apply(meanfldT, tgav) + meanfldT$r
## Allow for some roundoff error
ftol <- 1.0e-8 * mean(abs(tgav))
diff <- tscl - griddataT$vardata
maxdiff <- max(abs(diff))
expect_lte(maxdiff, ftol)
# and test with the actual function results
tscl2 <- fullgrids[[1]][, 1:Ngrid]
diff2 <- tscl2 - griddataT$vardata
maxdiff2 <- max(abs(diff2))
expect_lte(maxdiff2, ftol)
# finally make sure those estimates agree
expect_lte(abs(maxdiff - maxdiff2), 1e-6 * ftol)
})
test_that('Applying the pattern scaling and adding back the residuals recovers the original precipitations.',
{
prfield <- pscl_apply(meanfldP, tgav)
prscl <- pvarunconvert_fcn(prfield + meanfldP$r)
prfield <- pvarunconvert_fcn(prfield)
## Allow for some roundoff error
ftol <- 1.0e-8 * max(abs(prfield))
diff <- prscl - griddataP$vardata_raw
maxdiff <- max(abs(diff))
expect_lte(maxdiff, ftol)
# and test with the actual function results
pscl2 <- fullgrids[[1]][, (Ngrid + 1):(2*Ngrid)]
diff2 <- pscl2 - griddataP$vardata_raw
maxdiff2 <- max(abs(diff2))
expect_lte(maxdiff2, ftol)
# finally make sure those estimates are close
expect_lte(abs(maxdiff - maxdiff2), 1e-6 * ftol)
})
test_that('Full P grids produced do not feature negative precipitation.',
{
## pull off the P:
fullgridsP <- lapply(fullgrids, function(g){
g[, (Ngrid+1):(2*Ngrid)]}
)
minP <- lapply(fullgridsP, function(g){
apply(g,2,min)
})
ftol <- 0
# While precip SHOULD absolutely be >0 as
# an emergent property of the fldgen algorithm.
# we may want to give ourselves some
# numerical wiggle room in the future.
# for reference, ftol = -1e-8 is
# equivalent to saying precip just has
# to be > -0.000864 mm/day. Then we could
# add a manual adjustment to be > 0.
minminP.input <- min(minP[[1]])
minminP.gen1 <- min(minP[[2]])
minminP.gen2 <- min(minP[[3]])
minminP.gen3 <- min(minP[[4]])
expect_gte(minminP.input, ftol)
expect_gte(minminP.gen1, ftol)
expect_gte(minminP.gen2, ftol)
expect_gte(minminP.gen3, ftol)
})
test_that('All EOFs are normalized and orthogonal to one another.',
{
## check that all eofs are normalized and orthogonal to one another. We do
## this by multiplying the matrix of basis vectors by its transpose and
## comparing to the identity matrix.
tprod <- t(reof$rotation) %*% reof$rotation
maxdiff <- max(abs(tprod - diag(ncol(tprod))))
expect_lte(maxdiff, 1.0e-8)
})
test_that('EOFs other than PC0 have zero global average.',
{
## check that eofs other than PC0 have zero global average. (Technically
## this is guaranteed by the test above, but this test serves as independent
## confirmation.)
maxtg <- max( t(rep(griddataT$globalop,2)) %*% reof$rotation[,-1] )
expect_lte(maxtg, 1.0e-8)
})
test_that('Time series generated by mkcorrts method 1 have the correct properties.',
{
validate_mkcorrts()
validate_mkcorrts(reof$x)
})
test_that('Time series generated by mkcorrts method 2 have the correct properties.',
{
validate_mkcorrts(method=2)
validate_mkcorrts(reof$x, method=2)
})
test_that('PSD returned from training matches manually calculated.',
{
expect_equal(Fxmag, abs(Fx))
})
test_that('Residual field can be recovered from basis and rotated coordinates.',
{
## The reconstructed residuals from reof will be in the normal space.
## Therefore, will compare to input residual field in the normal space.
norm.input.resids <- cbind(
normalize.resids(inputresids = meanfldT$r,
empiricalcdf = emulator$tfuns$cdf)$rn,
normalize.resids(inputresids = meanfldP$r,
empiricalcdf = emulator$pfuns$cdf)$rn
)
## We will be a little lenient on this test because there is plenty of
## opportunity for roundoff errors to accumulate in this calculation.
resid.reconst <- reconst_fields(reof$rotation, reof$x)
diff <- resid.reconst - norm.input.resids
maxdiff <- max(abs(diff))
ftol <- 1.0e-6 * mean(abs(norm.input.resids))
expect_lte(maxdiff, ftol)
})
test_that('psdest produces equivalent results to manual calculation.',
{
## do analysis with psdest
psd1 <- psdest(reof)
## repeat it manually
Fx <- mvfft(reof$x)
Fxmag <- abs(Fx)
expect_equal(sqrt(psd1$psd), Fxmag)
## A list of two identical inputs should also give back the same result.
psd2 <- psdest(list(reof, reof))
expect_equal(sqrt(psd2$psd), Fxmag)
})
test_that('reducedEmulator object agrees with original emulator.',
{
reducedEmulator <- emulator_reducer(emulator)
# current approach to emulator_reducer will see
# the reducedEmulator have the same names as the
# original
expect_equal(names(reducedEmulator),
names(emulator))
# Generate new residuals using the reducedEmulator
# have to set seed to compare
set.seed(11)
resids_reduced <- generate.TP.resids(reducedEmulator, ngen = 2)
# Generate new residuals using the full emulator
# have to set seed to compare
set.seed(11)
resids_fullemulator <- generate.TP.resids(emulator, ngen = 2)
# Compare the residuals
expect_equal(resids_reduced[[1]], resids_fullemulator[[1]])
expect_equal(resids_reduced[[2]], resids_fullemulator[[2]])
# And compare the reconstructed full grids
fullgrids_reduced <- generate.TP.fullgrids(reducedEmulator,
resids_reduced,
tgav = reducedEmulator$tgav,
tvarunconvert_fcn = NULL,
pvarunconvert_fcn = exp,
reconstruction_function = pscl_apply)
fullgrids_reduced <- lapply(fullgrids_reduced$fullgrids,
function(g){
return(cbind(g[[1]], g[[2]]))
})
fullgrids_fullemulator <- generate.TP.fullgrids(emulator,
resids_fullemulator,
tgav = emulator$tgav,
tvarunconvert_fcn = NULL,
pvarunconvert_fcn = exp,
reconstruction_function = pscl_apply)
fullgrids_fullemulator <- lapply(fullgrids_fullemulator$fullgrids,
function(g){
return(cbind(g[[1]], g[[2]]))
})
expect_equal(fullgrids_reduced[[1]], fullgrids_fullemulator[[1]])
expect_equal(fullgrids_reduced[[2]], fullgrids_fullemulator[[2]])
}
)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.