tests/testthat/test_varfield.R

In JGCRI/fieldgenr: Generate Temperature Fields with Spatial and Temporal Correlation

#### 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]])

          }
          )