tests/testthat/test-ClimToProxyClim.R

In EarthSystemDiagnostics/sedproxy: Simulation of Sediment Archived Climate Proxy Records

library(sedproxy)
context("ClimToProxyClim")

# Fast = Slow if no mixed layer points ------
test_that("Fast = Slow if no mixed layer points", {

  clim.in <- N41.t21k.climate[nrow(N41.t21k.climate):1,] - 273.15
  clim.in <- ts(clim.in, start = -39)

  tpts <- round(seq(400, 22000, length.out = 100))

  set.seed(26052017)
  PFM.slow <- ClimToProxyClim(clim.signal = clim.in,
                              timepoints = tpts,
                              habitat.weights = N41.G.ruber.seasonality,
                              sed.acc.rate = rep(50, length(tpts)),
                              layer.width = 0,
                              sigma.meas = 0.46,
                              sigma.ind = 0,
                              n.samples = 30,
                              n.replicates = 1)
  set.seed(26052017)
  PFM.fast <- ClimToProxyClim(clim.signal = clim.in,
                              timepoints = tpts,
                              habitat.weights = N41.G.ruber.seasonality,
                              sed.acc.rate = 50,
                              layer.width = 0,
                              sigma.meas = 0.46,
                              sigma.ind = 0,
                              n.samples = 30,
                              n.replicates = 1)

  set.seed(26052017)
  PFM.slow.inf <- ClimToProxyClim(clim.signal = clim.in,
                              timepoints = tpts,
                              habitat.weights = N41.G.ruber.seasonality,
                              sed.acc.rate = rep(50, length(tpts)),
                              layer.width = 0,
                              sigma.meas = 0.46,
                              sigma.ind = 0,
                              n.samples = Inf,
                              n.replicates = 1)
  set.seed(26052017)
  PFM.fast.inf <- ClimToProxyClim(clim.signal = clim.in,
                              timepoints = tpts,
                              habitat.weights = N41.G.ruber.seasonality,
                              sed.acc.rate = 50,
                              layer.width = 0,
                              sigma.meas = 0.46,
                              sigma.ind = 0,
                              n.samples = Inf,
                              n.replicates = 1)

  expect_equal(object = data.frame(PFM.fast$simulated.proxy),
                         expected = data.frame(PFM.slow$simulated.proxy))
  expect_equal(object = data.frame(PFM.fast$everything),
               expected = data.frame(PFM.slow$everything))

  expect_equal(object = data.frame(PFM.fast.inf$simulated.proxy),
               expected = data.frame(PFM.slow.inf$simulated.proxy))
  expect_equal(object = data.frame(PFM.fast.inf$everything),
               expected = data.frame(PFM.slow.inf$everything))
})


test_that("Slow = Cached slow", {
  load("data/PFM.cache.Rdata")

  clim.in <- N41.t21k.climate[nrow(N41.t21k.climate):1,] - 273.15
  clim.in <- ts(clim.in, start = -39)

  tpts <- round(seq(4000, 20000, length.out = 10))

  set.seed(26052017)
  PFM.slow <- ClimToProxyClim(clim.signal = clim.in,
                              timepoints = tpts,
                              habitat.weights = N41.G.ruber.seasonality,
                              sed.acc.rate = rep(50, length(tpts)),
                              layer.width = 0,
                              sigma.meas = 0.46,
                              sigma.ind = 0,
                              n.samples = 30,
                              n.replicates = 1)
  #PFM.cache <- PFM.slow
  #save(PFM.cache, file = "tests/testthat/data/PFM.cache.Rdata")

  expect_equal(object = data.frame(PFM.cache$simulated.proxy),
               expected = data.frame(PFM.slow$simulated.proxy))
})

# zero mixing -----
test_that("zero mixing works", {

  clim.in <- N41.t21k.climate[nrow(N41.t21k.climate):1,] - 273.15
  #clim.in <- cbind(1:22040)
  clim.in <- ts(clim.in, start = 1)

  tpts <- round(seq(4000, 20000, length.out = 10))

  set.seed(26052017)
  PFM.slow <- ClimToProxyClim(clim.signal = clim.in,
                              timepoints = tpts,
                              #habitat.weights = rep(1/12, 12),
                              sed.acc.rate = rep(50, length(tpts)),
                              layer.width = 0,
                              bio.depth = 0,
                              sigma.meas = 0.46,
                              sigma.ind = 0,
                              n.samples = 30,
                              n.replicates = 1)

  PFM.fast <- ClimToProxyClim(clim.signal = clim.in,
                              timepoints = tpts,
                              #habitat.weights = rep(1/12, 12),
                              sed.acc.rate = 50,
                              layer.width = 0,
                              bio.depth = 0,
                              sigma.meas = 0.46,
                              sigma.ind = 0,
                              n.samples = 30,
                              n.replicates = 1)

  expect_equal(object = PFM.slow$simulated.proxy$proxy.bt,
               expected = rowMeans(clim.in)[tpts])

  expect_equal(object = PFM.fast$simulated.proxy$proxy.bt,
               expected = rowMeans(clim.in)[tpts])
})

# negative times --------
test_that("negative times work -fast", {

  set.seed(26052017)
  clim.in <- matrix(rnorm(101*12), ncol = 12)
  clim.in <- ts(clim.in, start = -49)

  tpts <- c(1, 10, 100, 200)-51

  PFM <- ClimToProxyClim(clim.signal = clim.in,
                         timepoints = tpts,
                         calibration.type = "identity",
                         sed.acc.rate = 50,
                         bio.depth = 0, layer.width = 0,
                         n.samples = 30, n.replicates = 10, n.bd = 3,
                         top.of.core = -49)

  PFM$simulated.proxy$clim.signal.ann
  rowMeans(clim.in[time(clim.in)%in%tpts,])


  expect_equal(object = PFM$simulated.proxy$clim.signal.ann,
               expected = rowMeans(clim.in[time(clim.in)%in%tpts,]))

})

test_that("negative times work - slow", {
  
  set.seed(26052017)
  clim.in <- matrix(rnorm(101*12), ncol = 12)
  clim.in <- ts(clim.in, start = -49)
  
  tpts <- c(1, 10, 100, 200)-51
  
  PFM <- ClimToProxyClim(clim.signal = clim.in,
                         timepoints = tpts,
                         calibration.type = "identity",
                         sed.acc.rate = rep(50, length(tpts)),
                         bio.depth = 0, layer.width = 0,
                         n.samples = 30, n.replicates = 10, n.bd = 3,
                         top.of.core = -49)
  
  PFM$simulated.proxy$clim.signal.ann
  rowMeans(clim.in[time(clim.in)%in%tpts,])
  
  
  expect_equal(object = PFM$simulated.proxy$clim.signal.ann,
               expected = rowMeans(clim.in[time(clim.in)%in%tpts,]))
  
})

# mean bias repwise -----
test_that("meas.bias applied repwise", {
  
  clim.in <- ts(matrix(rep(1, 12*1e04), ncol = 12))
  PFM.fast <- ClimToProxyClim(clim.in,
                         timepoints = seq(1000, 6000, 500),
                         meas.bias = 3,
                         n.replicates = 3,
                        )
  
  #PlotPFMs(PFM.fast, max.replicates = 30)
  
  PFM.sub.fast <- PFM.fast$everything[PFM.fast$everything$stage == "simulated.proxy", ]
  
  PFM.slow <- ClimToProxyClim(clim.in,
                         timepoints = seq(1000, 6000, 500),
                         sed.acc.rate = rep(50, 11),
                         meas.bias = 3,
                         n.replicates = 3,
                         n.samples = 30)
  #PlotPFMs(PFM.slow)
  
  PFM.sub.slow <- PFM.slow$everything[PFM.slow$everything$stage == "simulated.proxy", ]
  
  
  PFM.mult <- ClimToProxyClim(clim.in,
                              timepoints = seq(1000, 6000, 500),
                              sed.acc.rate = rep(50, 11),
                              meas.bias = log(3),
                              n.replicates = 3,
                              n.samples = 30,
                              noise.type = "multiplicative")
  #PlotPFMs(PFM.mult)
  
  PFM.sub.mult <- PFM.mult$everything[PFM.mult$everything$stage == "simulated.proxy", ]
  
  
  expect_equal(as.numeric(tapply(PFM.sub.fast$value, PFM.sub.fast$replicate, sd)),
               c(0,0,0))
  
  expect_equal(as.numeric(tapply(PFM.sub.slow$value, PFM.sub.slow$replicate, sd)),
               c(0,0,0))
  
  expect_equal(as.numeric(tapply(PFM.sub.mult$value, PFM.sub.mult$replicate, sd)),
               c(0,0,0))
  
  })


# mixed layer -----

test_that("mixed layer modelled", {
  
  clim.in <- cbind(1:22040)
  clim.in <- ts(clim.in, start = 1)
  
  tpts <- round(seq(10, 5000, length.out = 10))
  
  set.seed(26052017)
  PFM.slow <- ClimToProxyClim(clim.signal = clim.in,
                              timepoints = tpts,
                              #habitat.weights = rep(1/12, 12),
                              sed.acc.rate = rep(2, length(tpts)),
                              layer.width = 0,
                              bio.depth = 10,
                              sigma.meas = 0.46,
                              sigma.ind = 0,
                              n.samples = 30,
                              n.replicates = 1)
  set.seed(26052017)
  PFM.fast <- ClimToProxyClim(clim.signal = clim.in,
                              timepoints = tpts,
                              #habitat.weights = rep(1/12, 12),
                              sed.acc.rate = 2,
                              layer.width = 0,
                              bio.depth = 10,
                              sigma.meas = 0.46,
                              sigma.ind = 0,
                              n.samples = 30,
                              n.replicates = 1)
  
  
  expect_equal(object = max(diff(PFM.slow$simulated.proxy$proxy.bt)),
               expected = 0)
  
  expect_equal(object = max(diff(PFM.fast$simulated.proxy$proxy.bt)),
               expected = 0)
})

test_that("mixed layer, varying sed.acc.rate", {
  
  ## Does not work with varying S
  
  clim.in <- cbind(1:22040)
  clim.in <- ts(clim.in, start = 1)
  
  tpts <- round(seq(10, 5000, length.out = 10))
  s.rates <- c(6,6,6,2,2,2, 2,2,2,2)
  #s.rates <- rep(2.2, 10)
  
  set.seed(26052017)
  PFM.slow <- ClimToProxyClim(clim.signal = clim.in,
                              timepoints = tpts,
                              #habitat.weights = rep(1/12, 12),
                              sed.acc.rate = s.rates,
                              layer.width = 0,
                              bio.depth = 10,
                              sigma.meas = 0.46,
                              sigma.ind = 0,
                              n.samples = 30,
                              n.replicates = 1)
  
  # PlotPFMs(PFM.slow) +
  #  ggplot2::geom_vline(xintercept = 1000* 10 / min(s.rates[1:4]))
  # 
  expect_equal(object = max(diff(PFM.slow$simulated.proxy$proxy.bt)),
               expected = 0)
})

# example from paper --------
test_that("example from paper works", {
  
  load("data/PFM.ex.Rdata")
  PFM.cache <- PFM.ex
  clim.in <- N41.t21k.climate[nrow(N41.t21k.climate):1,] - 273.15
  clim.in <- ts(clim.in, start = 1)
  
  
  tpts <- N41.proxy$Published.age
  s.rates <- N41.proxy$Sed.acc.rate.cm.ka
  
  
  set.seed(26052017)
  PFM.ex <- ClimToProxyClim(clim.signal = clim.in,
                              timepoints = tpts,
                              habitat.weights = N41.G.ruber.seasonality,
                              sed.acc.rate = s.rates,
                              layer.width = 1,
                              bio.depth = 10,
                              sigma.meas = 0.46,
                              sigma.ind = 0,
                              n.samples = 30,
                              n.replicates = 1)
  
  
  #save(PFM.ex, file = "tests/testthat/data/PFM.ex.Rdata")
  
  expect_equal(object = data.frame(PFM.cache$simulated.proxy),
               expected = data.frame(PFM.ex$simulated.proxy))
  
  p <- PlotPFMs(PFM.ex) 
  testthat::expect_equal(class(p), c("gg", "ggplot"))
  
  
})