R/SimClimMat.R

In EarthSystemDiagnostics/spectraluncertainty:

#' Integrate a Power Law Function
#'
#' @param a,b parameters of a power law a*nu^b
#' @param nu1 lowest positive frequency
#' @param nu2 highest positive frequency
#'
#' @return the definite integral over postive frequencies
#' @export
#'
#' @examples
#' IntegratePowerlaw(a = 0.1, b = 1, nu1 = 1/1000, nu2 = 1)
IntegratePowerlaw <- function(a, b, nu1, nu2){
  if (b < 0){
    a*(nu1^(1-b)/(b-1)-nu2^(1-b)/(b-1))
  }else if (b == 1){
    a*(log(nu2)-log(nu1))
  }else{
    a*(nu1^(1-b)/(b-1)-nu2^(1-b)/(b-1))
  }
}


#' Simulate a climate matrix for sedproxy using parameters from PSEM
#'
#' @param sim.pars
#'
#' @return A list with the simulated climate matrix and seasonal cycle
#' @export
#' @examples
#' spec.pars <- GetSpecPars("Mg_Ca")
#' clim.mat <- SimClimMat(spec.pars)
SimClimMat <- function(sim.pars, n.months = 12){

  # mean seasonal cycle over full timeseries
  mean.seas.cycle <- with(sim.pars, cos(seq(pi, 3*pi - 2*pi/n.months, length.out = n.months)))
  # subtract mean because of only 12 months
  mean.seas.cycle <- mean.seas.cycle - mean(mean.seas.cycle)
  mean.seas.cycle <- mean.seas.cycle / sd(mean.seas.cycle)
  mean.seas.cycle <- mean.seas.cycle * sqrt(VarSine(sim.pars$seas.amp))

  # calc time varying amp of seas cycle
  tax <-  with(sim.pars, seq(-pi*T*nu_a, pi*T*nu_a, length.out = T) + phi_a + pi)
  d.seas.cycle <-  with(sim.pars, 1+2*((sqrt(sig.sq_a) * (cos(tax))) / sqrt(2)))

  seas.cycle <- t(mean.seas.cycle %*% t((d.seas.cycle)))

  #Annual means
  alpha <-  (1/2)^sim.pars$s.beta * sim.pars$s.0.5

  # Get the total variance for a powerlaw timeseries with length T and a, b parameters
  total.sd <- with(sim.pars, {sqrt(2 * IntegratePowerlaw(a = clim.spec.fun.args$a,
                                                         b = clim.spec.fun.args$b,
                                                         nu1 = 1/(T*12), nu2 = 12))})

  # Generate a monthly resolution powerlaw timeseries
  # cl.1 <- with(sim.pars, {matrix(rep(FastSimPowerlaw(beta = clim.spec.fun.args$b, N = T), 12) *
  #                                total.sd, ncol = 12, byrow = FALSE)})

  cl.1 <- with(sim.pars, {matrix(FastSimPowerlaw(beta = clim.spec.fun.args$b, N = T*12) *
                                   total.sd, ncol = 12, byrow = TRUE)})

  # Add the seasonal cycle
  cl <- cl.1 + seas.cycle
  cl <- ts(cl)

  # return(list(cl = cl, cl.1 = cl.1, seas.cycle = seas.cycle,
  #             d.seas.cycle = d.seas.cycle,
  #             mean.seas.cycle = mean.seas.cycle))
  return(list(cl = cl, seas.cycle = seas.cycle))
}

#' Simulate a climate matrix for sedproxy using parameters from PSEM and an
#' empirical climate spectrum
#'
#' @param spec Spectrum of stochastic climate
#' @param sim.pars Parameters of seasonal cycle
#'
#' @return A list with the simulated climate matrix and seasonal cycle
#' @export
#'
#' @examples
#' spec.pars <- GetSpecPars("Mg_Ca")
#' f <- GetNu(T = 10000, delta_t = 1)
#' spec <- ModelSpectrum(f, latitude = 20)
#' spec$spec[is.na(spec$spec)] <- 0
#' clim.mat <- SimClimMatEmpirical(spec, spec.pars)
#' plot(clim.mat$cl[, 1])
SimClimMatEmpirical <- function(spec, sim.pars, n.months = 12){

  # mean seasonal cycle over full timeseries
  mean.seas.cycle <- with(sim.pars, cos(seq(pi, 3*pi - 2*pi/n.months, length.out = n.months)))
  # subtract mean because of only 12 months
  mean.seas.cycle <- mean.seas.cycle - mean(mean.seas.cycle)
  mean.seas.cycle <- mean.seas.cycle / sd(mean.seas.cycle)
  mean.seas.cycle <- mean.seas.cycle * sqrt(VarSine(sim.pars$seas.amp))

  # calc time varying amp of seas cycle
  tax <-  with(sim.pars, seq(-pi*T*nu_a, pi*T*nu_a, length.out = T) + phi_a + pi)
  d.seas.cycle <-  with(sim.pars, 1+2*((sqrt(sig.sq_a) * (cos(tax))) / sqrt(2)))

  seas.cycle <- t(mean.seas.cycle %*% t((d.seas.cycle)))

  #Annual means
  alpha <-  (1/2)^sim.pars$s.beta * sim.pars$s.0.5


  # Generate a monthly resolution timeseries from empirical spec

  cl.1 <- with(sim.pars, {matrix(
    rep(PaleoSpec::SimFromEmpiricalSpec(spec = spec, N = T), 12),
    ncol = 12, byrow = FALSE
    )})


  # Add the seasonal cycle
  cl <- cl.1 + seas.cycle
  cl <- ts(cl)

  # return(list(cl = cl, cl.1 = cl.1, seas.cycle = seas.cycle,
  #             d.seas.cycle = d.seas.cycle,
  #             mean.seas.cycle = mean.seas.cycle))
  return(list(cl = cl, seas.cycle = seas.cycle))
}



#' A faster version of SimPowerlaw
#'
#' @inherit PaleoSpec::SimPowerlaw
#' @importFrom fftw planFFT
#' @importFrom fftw FFT
#' @export
FastSimPowerlaw <- function(beta, N)
{
  N2 <- (3^ceiling(log(N, base = 3)))
  df  <- 1 / N2
  f <- seq(from = df, to = 1/2, by = df)
  Filter <- sqrt(1/(f^beta))
  Filter <- c(max(Filter), Filter, rev(Filter))
  x   <- rnorm(N2, 1)
  pl <- fftw::planFFT(x)
  fx  <- fftw::FFT(x, plan = pl)
  ffx <- fx * Filter
  result <- Re(fftw::FFT(ffx, plan = pl, inverse = TRUE))[1:N]
  result <- result - mean(result)
  result <- result / sd(result)
  return(result)
}