R/packagefiles.R

In bayesianETAS: Bayesian Estimation of the Temporal and Spatio-Temporal ETAS Models for Earthquake Occurrences

.is_missing_initval <- function(initval) {
  length(initval) == 1 && is.na(initval)
}

.ETAS_C_UPPER <- 10
.ETAS_P_UPPER <- 10

.check_numeric_vector <- function(x, name, len = NULL) {
  if (!is.numeric(x) || any(!is.finite(x))) {
    stop(sprintf("%s must be a finite numeric vector", name), call. = FALSE)
  }
  if (!is.null(len) && length(x) != len) {
    stop(sprintf("%s must have length %s", name, len), call. = FALSE)
  }
}

.check_numeric_scalar <- function(x, name) {
  .check_numeric_vector(x, name, len = 1)
}

.check_positive_scalar <- function(x, name) {
  .check_numeric_scalar(x, name)
  if (x <= 0) {
    stop(sprintf("%s must be positive", name), call. = FALSE)
  }
}

.check_nonnegative_scalar <- function(x, name) {
  .check_numeric_scalar(x, name)
  if (x < 0) {
    stop(sprintf("%s must be non-negative", name), call. = FALSE)
  }
}

.check_model_parameters <- function(mu, K, alpha, cp, p, c_lower = 0) {
  .check_nonnegative_scalar(mu, "mu")
  .check_nonnegative_scalar(K, "K")
  .check_nonnegative_scalar(alpha, "alpha")
  .check_numeric_scalar(cp, "c")
  .check_nonnegative_scalar(c_lower, "c_lower")
  if (cp <= 0) {
    stop("c must be positive", call. = FALSE)
  }
  if (cp <= c_lower) {
    stop("c must be greater than c_lower", call. = FALSE)
  }
  if (cp >= .ETAS_C_UPPER) {
    stop("c must be less than 10", call. = FALSE)
  }
  .check_numeric_scalar(p, "p")
  if (p <= 1) {
    stop("p must be greater than 1", call. = FALSE)
  }
  if (p >= .ETAS_P_UPPER) {
    stop("p must be less than 10", call. = FALSE)
  }
}

.check_catalog <- function(ts, ms, M0, maxTime, require_strict = FALSE) {
  .check_numeric_vector(ts, "ts")
  .check_numeric_vector(ms, "ms", len = length(ts))
  .check_numeric_scalar(M0, "M0")
  .check_numeric_scalar(maxTime, "maxTime")
  if (length(ts) < 1) {
    stop("ts must contain at least one event", call. = FALSE)
  }
  if (maxTime < max(ts)) {
    stop("maxTime must be greater than or equal to max(ts)", call. = FALSE)
  }
  if (is.unsorted(ts, strictly = FALSE)) {
    stop("ts must be sorted in non-decreasing event-time order", call. = FALSE)
  }
  if (require_strict && any(diff(ts) <= 0)) {
    stop("ts must be strictly increasing after tie handling", call. = FALSE)
  }
  if (any(ms < M0)) {
    stop("all magnitudes in ms must be greater than or equal to M0", call. = FALSE)
  }
}

.check_spatial_catalog <- function(ts, ms, xs, ys, M0, maxTime, area) {
  .check_catalog(ts, ms, M0, maxTime, require_strict = TRUE)
  .check_numeric_vector(xs, "xs", len = length(ts))
  .check_numeric_vector(ys, "ys", len = length(ts))
  .check_positive_scalar(area, "area")
}

.check_burnin_sims <- function(sims, burnin) {
  .check_numeric_scalar(sims, "sims")
  .check_numeric_scalar(burnin, "burnin")
  if (sims != as.integer(sims) || sims < 1) {
    stop("sims must be a positive integer", call. = FALSE)
  }
  if (burnin != as.integer(burnin) || burnin < 0) {
    stop("burnin must be a non-negative integer", call. = FALSE)
  }
  if (burnin >= sims) {
    stop("burnin must be less than sims", call. = FALSE)
  }
}

.drop_burnin <- function(samples, burnin) {
  if (burnin > 0) {
    samples <- samples[-seq_len(burnin), , drop = FALSE]
  }
  samples
}

.validate_initval <- function(initval, len, name) {
  if (.is_missing_initval(initval)) {
    return(NULL)
  }
  .check_numeric_vector(initval, name, len = len)
  initval
}

.resolve_sigma2_prior <- function(xs, ys, sigma2_prior_shape,
                                  sigma2x_prior_rate,
                                  sigma2y_prior_rate) {
  .check_numeric_scalar(sigma2_prior_shape, "sigma2_prior_shape")
  if (sigma2_prior_shape <= 2) {
    stop("sigma2_prior_shape must be greater than 2", call. = FALSE)
  }
  if (is.null(sigma2x_prior_rate)) {
    x_range <- diff(range(xs))
    if (!is.finite(x_range) || x_range <= 0) {
      stop("sigma2x_prior_rate must be provided when xs has zero or non-finite range", call. = FALSE)
    }
    s0x <- 0.1 * x_range
    sigma2x_prior_rate <- 2 * s0x^2
  }
  if (is.null(sigma2y_prior_rate)) {
    y_range <- diff(range(ys))
    if (!is.finite(y_range) || y_range <= 0) {
      stop("sigma2y_prior_rate must be provided when ys has zero or non-finite range", call. = FALSE)
    }
    s0y <- 0.1 * y_range
    sigma2y_prior_rate <- 2 * s0y^2
  }
  .check_positive_scalar(sigma2x_prior_rate, "sigma2x_prior_rate")
  .check_positive_scalar(sigma2y_prior_rate, "sigma2y_prior_rate")
  list(shape = sigma2_prior_shape,
       rate_x = sigma2x_prior_rate,
       rate_y = sigma2y_prior_rate)
}

.strict_time_jitter <- function(ts, maxTime) {
  gaps <- diff(unique(ts))
  positive_gaps <- gaps[gaps > 0]
  max_run <- max(rle(ts)$lengths)
  if (length(positive_gaps) > 0) {
    eps <- min(positive_gaps) / (10 * (max_run + 1))
  } else {
    eps <- max(maxTime, 1) / (1e6 * (length(ts) + 1))
  }
  ts_new <- ts
  for (i in seq_along(ts_new)[-1]) {
    if (ts_new[i] <= ts_new[i - 1]) {
      ts_new[i] <- ts_new[i - 1] + eps
    }
  }
  if (max(ts_new) > maxTime) {
    stop("deterministic tie jitter would move an event beyond maxTime; increase maxTime or resolve ties before fitting", call. = FALSE)
  }
  list(ts = ts_new,
       jitter_epsilon = eps,
       jittered = TRUE,
       n_tied_events = sum(duplicated(ts)))
}

.resolve_c_lower <- function(ts, maxTime, c_lower = NULL) {
  if (!is.null(c_lower)) {
    .check_positive_scalar(c_lower, "c_lower")
    if (c_lower >= .ETAS_C_UPPER) {
      stop("c_lower must be less than 10", call. = FALSE)
    }
    return(c_lower)
  }
  positive_gaps <- diff(ts)
  positive_gaps <- positive_gaps[positive_gaps > 0]
  if (length(positive_gaps) > 0) {
    resolved <- min(positive_gaps) / 10
    if (resolved >= .ETAS_C_UPPER) {
      resolved <- .ETAS_C_UPPER / 10
    }
    return(resolved)
  }
  max(maxTime, 1) * 1e-8
}

.prepare_catalog <- function(ts, ms, M0, maxTime,
                             handle_ties = c("error", "jitter"),
                             c_lower = NULL) {
  handle_ties <- match.arg(handle_ties)
  .check_catalog(ts, ms, M0, maxTime, require_strict = FALSE)
  if (any(diff(ts) == 0)) {
    if (handle_ties == "error") {
      stop("duplicate event times detected; continuous-time ETAS requires strictly increasing event times. Use handle_ties = \"jitter\" to apply deterministic preprocessing.", call. = FALSE)
    }
    tie_info <- .strict_time_jitter(ts, maxTime)
    ts <- tie_info$ts
  } else {
    tie_info <- list(jittered = FALSE,
                     jitter_epsilon = NA_real_,
                     n_tied_events = 0)
  }
  .check_catalog(ts, ms, M0, maxTime, require_strict = TRUE)
  c_lower <- .resolve_c_lower(ts, maxTime, c_lower)
  list(ts = ts, c_lower = c_lower, tie_info = tie_info)
}

.safe_positive <- function(x, floor = 1e-8) {
  ifelse(is.finite(x) & x > floor, x, floor)
}

.logit_bounded <- function(x, lower, upper) {
  z <- (x - lower) / (upper - lower)
  z <- pmin(pmax(z, 1e-8), 1 - 1e-8)
  log(z / (1 - z))
}

.inv_logit_bounded <- function(eta, lower, upper) {
  lower + (upper - lower) / (1 + exp(-eta))
}

.temporal_to_eta <- function(params, c_lower, zeromu = FALSE) {
  params[4] <- pmin(max(params[4], c_lower * 2), .ETAS_C_UPPER - 1e-6)
  params[5] <- pmin(max(params[5], 1 + 1e-6), .ETAS_P_UPPER - 1e-6)
  if (zeromu) {
    return(log(c(.safe_positive(params[2]),
                 .safe_positive(params[3]),
                 exp(.logit_bounded(params[4], c_lower, .ETAS_C_UPPER)),
                 exp(.logit_bounded(params[5], 1, .ETAS_P_UPPER)))))
  }
  log(c(.safe_positive(params[1]),
        .safe_positive(params[2]),
        .safe_positive(params[3]),
        exp(.logit_bounded(params[4], c_lower, .ETAS_C_UPPER)),
        exp(.logit_bounded(params[5], 1, .ETAS_P_UPPER))))
}

.eta_to_temporal <- function(eta, c_lower, zeromu = FALSE) {
  if (zeromu) {
    return(c(0, exp(eta[1]), exp(eta[2]), .inv_logit_bounded(eta[3], c_lower, .ETAS_C_UPPER), .inv_logit_bounded(eta[4], 1, .ETAS_P_UPPER)))
  }
  c(exp(eta[1]), exp(eta[2]), exp(eta[3]), .inv_logit_bounded(eta[4], c_lower, .ETAS_C_UPPER), .inv_logit_bounded(eta[5], 1, .ETAS_P_UPPER))
}

.spatial_to_eta <- function(params, c_lower, zeromu = FALSE) {
  params[4] <- max(params[4], c_lower * 2)
  base <- .temporal_to_eta(params[1:5], c_lower, zeromu)
  c(base, log(.safe_positive(params[6])), log(.safe_positive(params[7])))
}

.eta_to_spatial <- function(eta, c_lower, zeromu = FALSE) {
  temporal_len <- if (zeromu) 4 else 5
  c(.eta_to_temporal(eta[seq_len(temporal_len)], c_lower, zeromu),
    exp(eta[temporal_len + 1]), exp(eta[temporal_len + 2]))
}

.temporal_start_values <- function(ts, marks, maxTime, initval, c_lower) {
  n <- length(ts)
  positive_gaps <- diff(ts)
  positive_gaps <- positive_gaps[positive_gaps > 0]
  scale_c <- if (length(positive_gaps) > 0) stats::median(positive_gaps) else maxTime / max(n, 1)
  base_c <- max(c_lower * 5, min(maxTime / 20, max(scale_c, c_lower * 10)))
  starts <- list(
    c(n / maxTime, 0.2, 0.5, base_c, 1.5),
    c(n / maxTime * 0.5, 0.05, 0.2, max(c_lower * 5, base_c / 2), 1.3),
    c(n / maxTime * 0.5, 0.1, 1.0, base_c, 1.7),
    c(n / maxTime * 0.25, 0.5, 0.5, max(c_lower * 10, base_c * 2), 2.0),
    c(n / maxTime, 0.02, 1.5, base_c, 1.4),
    c(n / maxTime * 0.1, 0.2, 1.5, max(c_lower * 5, base_c / 4), 1.8),
    c(n / maxTime, 0.5, 0.1, max(c_lower * 20, base_c * 4), 2.5)
  )
  if (!is.null(initval)) {
    starts <- c(list(initval), starts)
  }
  unique(do.call(rbind, starts))
}

.spatial_start_values <- function(ts, marks, xs, ys, maxTime, initval, c_lower, temporal_seed = NULL) {
  s0x <- max(stats::var(xs), (diff(range(xs)) * 0.1)^2, 1e-4)
  s0y <- max(stats::var(ys), (diff(range(ys)) * 0.1)^2, 1e-4)
  temporal_starts <- .temporal_start_values(ts, marks, maxTime, temporal_seed, c_lower)
  starts <- lapply(seq_len(min(nrow(temporal_starts), 5)), function(i) {
    c(temporal_starts[i, ], s0x, s0y)
  })
  starts <- c(starts,
              list(c(length(ts) / maxTime, 0.1, 0.5, max(c_lower * 10, maxTime / 100), 1.5, s0x / 2, s0y / 2),
                   c(length(ts) / maxTime * 0.5, 0.2, 1.0, max(c_lower * 10, maxTime / 200), 2.0, s0x * 2, s0y * 2)))
  if (!is.null(initval)) {
    starts <- c(list(initval), starts)
  }
  unique(do.call(rbind, starts))
}

.attach_fit_attrs <- function(fit, c_lower, tie_info) {
  attr(fit, "c_lower") <- c_lower
  attr(fit, "tie_handling") <- tie_info
  fit
}

#' Draws samples from the posterior distribution of the ETAS model
#'

#' This function implements the latent variable MCMC scheme from (Ross 2019) which draws samples from the Bayesian posterior distribution of the Epidemic Type Aftershock Sequence (ETAS)  model. Either the temporal ETAS or the spatio-temporal (with a bivariate Gaussian trigger kernel) can be sampled.

#' \deqn{}

#'The ETAS model is widely used to quantify the degree of seismic activity in a geographical region, and to forecast the occurrence of future mainshocks and aftershocks (Ross 2019). The temporal ETAS model is a point process where the probability of an earthquake occurring at time \eqn{t} depends on the previous seismicity \eqn{H_t}{Ht}, and is defined by the conditional intensity function:

#' \deqn{ \lambda(t|H_t) = \mu + \sum_{t[i] < t}  \kappa(m[i]|K,\alpha) h(t-t[i]|c,p)}{ \lambda(t|Ht) = \mu + \sum \kappa(m[i]|K,\alpha) h(t[i]|c,p)s(x-x[i],y-y[i],}

#' where

#' \deqn{\kappa(m_i|K,\alpha) = Ke^{\alpha \left( m_i-M_0 \right)}}{\kappa(m[i]|K,\alpha) = K * exp(\alpha(m[i]-M0))}

#' and

#' \deqn{ h(t_i|c,p) = \frac{(p-1)c^{p-1}}{(t-t_i+c)^p}}{ h(t[i]|c,p) = (p-1) * c^(p-1) * (t-t[i]+c)^(-p)}

#' and s() is a spatial triggering kernel. The public spatio-temporal fitting routines currently use a Gaussian spatial triggering kernel with a homogeneous spatial background over the supplied area and no spatial edge correction. s() can also be omitted to fit the temporal ETAS model without a spatial component.

#' The summation is over all previous earthquakes that occurred in the region, with the i'th such earthquake occurring at time \eqn{ t_i}{t[i]} at spatial location (\eqn{ x_i}{x[i]},,\eqn{ y_i}{y[i]}) aand having magnitude \eqn{ m_i}{m[i]}. The quantity \eqn{M_0}{M0} denotes the magnitude of completeness of the catalog, so that \eqn{m_i \geq  M_0}{m[i] \ge M0} for all i. The temporal ETAS model has 5 parameters: \eqn{\mu} controls the background rate of seismicity, \eqn{ K} and \eqn{ \alpha} determine the productivity (average number of aftershocks) of an earthquake with magnitude \eqn{m}, and \eqn{c} and \eqn{ p}  are the parameters of the Modified Omori Law (which has here been normalized to integrate to 1) and represent the speed at which the aftershock rate decays over time. Each earthquake is assumed to have a magnitude which is an independent draw from the Gutenberg-Richter law \eqn{ p(m_i) = \beta e^{\beta(m_i-M_0)}}{ p(m) = \beta * exp(\beta(m-M0)}.
#'
#' In the spatio-temporal version of the model, s() is the multivariate Gaussian with diagonal covariance parameters \eqn{\sigma^2_x} and \eqn{\sigma^2_y}. The spatial background is homogeneous over the supplied area and no spatial edge correction is implemented.

#' @param ts Vector containing the earthquake times
#' @param ms Vector containing the earthquake magnitudes
#' @param maxTime Length of the time window [0,maxTime] the catalog was observed over. If not specified, will be taken as the time of the last earthquake.
#' @param M0 Magnitude of completeness.
#' @param xs (For the spatio-temporal ETAS model only), a vector containing the x coordinates of the earthquakes
#' @param ys (For the spatio-temporal ETAS model only), a vector containing the y coordinates of the earthquakes
#' @param area (For the spatio-temporal ETAS model only), a single number containing the area of the rectangular spatial region containing the earthquakes. For example if the earthquakes are defined on [28,35]x[-119,-128] then the area would be 7x9 = 63. If this is not specified then it will be calculated based on the distance between the furthest apart earthquakes in the catalog.

#' @param sims Number of posterior samples to draw
#' @param burnin Number of burnin samples. If (eg) sims=5000 and burnin=2000 then 7000 samples will be drawn with the first 2000 discarded
 #' @param initval Initial value at which to start the estimation. If specified, should be a vector, with elements (mu, K, alpha, c, p) for the temporal model and  (mu, K, alpha, c, p, sigma2x,sigma2y) for the spatio-temporal model. If unspecified, the sampler will be initialized at the maximum likelihood estimate of the model parameters
#' @param sigma2_prior_shape Shape parameter for the inverse-gamma prior on Gaussian spatial variances, using a shape/rate convention. Must be greater than 2.
#' @param sigma2x_prior_rate Rate parameter for the inverse-gamma prior on \eqn{\sigma^2_x}. If NULL, a scale-aware default equal to \eqn{2 * (0.1 * diff(range(xs)))^2} is used.
#' @param sigma2y_prior_rate Rate parameter for the inverse-gamma prior on \eqn{\sigma^2_y}. If NULL, a scale-aware default equal to \eqn{2 * (0.1 * diff(range(ys)))^2} is used.
#' @param handle_ties How to handle duplicate event times. The default "error" stops because continuous-time ETAS requires strictly increasing event times. "jitter" applies a deterministic, order-preserving jitter before fitting.
#' @param c_lower Lower bound for the Omori c parameter. If NULL, a data-scale default equal to one tenth of the smallest positive event-time gap after tie handling is used.
#' @return A data frame containing the posterior samples. Each row is a single sample, and the columns correspond to (mu, K, alpha, c, p) for the temporal model and  (mu, K, alpha, c, p,sigma2x,sigma2y) for the spatio-temporal model

#' @examples
#'\donttest{
#' #temporal model
#' ts <- c(1, 2, 3, 5, 8)
#' ms <- c(3.0, 3.1, 3.3, 3.2, 3.0)
#' samples <- estimateETAS(ts, ms, M0=3, maxTime=10, sims=20, burnin=5,
#'     initval=c(0.1, 0.1, 0.5, 0.2, 1.5))
#' apply(samples, 2, median) #point estimate
#'
#' #spatio-temporal model
#' xs <- c(1, 2, 3, 4, 5)
#' ys <- c(1, 1.5, 2, 1.5, 1)
#' samples <- estimateETAS(ts, ms, M0=3, maxTime=10, xs=xs, ys=ys,
#'     area=25, sims=20, burnin=5,
#'     initval=c(0.1, 0.1, 0.5, 0.2, 1.5, 1, 1))
#' apply(samples, 2, median) #point estimate
#' }
#'
#' @author Gordon J Ross
#' @references Gordon J. Ross - Bayesian Estimation of the ETAS Model for Earthquake Occurrences (2019)
#' @export
estimateETAS <- function(ts, ms, M0,maxTime=NULL,xs=NULL, ys=NULL,area=NULL,sims=5000, burnin=1000, initval=NA, sigma2_prior_shape=3, sigma2x_prior_rate=NULL, sigma2y_prior_rate=NULL, handle_ties=c("error","jitter"), c_lower=NULL) {
  if (is.null(maxTime)) {maxTime <- max(ts)}
  .check_burnin_sims(sims, burnin)
  prep <- .prepare_catalog(ts, ms, M0, maxTime, handle_ties, c_lower)
  ts <- prep$ts
  c_lower <- prep$c_lower

  if (is.null(xs) && is.null(ys)) {
    out <- estimateTemporalETAS(ts, ms, maxTime, M0, numMCMCSamples=200,sims=sims, burnin=burnin, initval=initval, c_lower=c_lower)
    return(.attach_fit_attrs(out, c_lower, prep$tie_info))
  } else {
    if (is.null(ys) || is.null(xs)) {
      stop("if fitting spatial ETAS then both xs and ys must be specified", call. = FALSE)
    }
    if (is.null(area)) { area <- (max(xs)-min(xs)) * (max(ys)-min(ys))}
    out <- estimateSpatialETAS(ts, ms, xs, ys, maxTime, M0,area,numMCMCSamples=200,sims=sims, burnin=burnin, initval=initval, gaussian=TRUE, sigma2_prior_shape=sigma2_prior_shape, sigma2x_prior_rate=sigma2x_prior_rate, sigma2y_prior_rate=sigma2y_prior_rate, c_lower=c_lower)
    return(.attach_fit_attrs(out, c_lower, prep$tie_info))
  }
}



#' Estimate the parameters of the ETAS model using maximum likelihood.
#'

#'The ETAS model is widely used to quantify the degree of seismic activity in a geographical region, and to forecast the occurrence of future mainshocks and aftershocks (Ross 2019). The temporal ETAS model is a point process where the probability of an earthquake occurring at time \eqn{t} depends on the previous seismicity \eqn{H_t}{Ht}, and is defined by the conditional intensity function:

#' \deqn{ \lambda(t|H_t) = \mu + \sum_{t[i] < t}  \kappa(m[i]|K,\alpha) h(t-t[i]|c,p)}{ \lambda(t|Ht) = \mu + \sum \kappa(m[i]|K,\alpha) h(t[i]|c,p)s(x-x[i],y-y[i],}

#' where

#' \deqn{\kappa(m_i|K,\alpha) = Ke^{\alpha \left( m_i-M_0 \right)}}{\kappa(m[i]|K,\alpha) = K * exp(\alpha(m[i]-M0))}

#' and

#' \deqn{ h(t_i|c,p) = \frac{(p-1)c^{p-1}}{(t-t_i+c)^p}}{ h(t[i]|c,p) = (p-1) * c^(p-1) * (t-t[i]+c)^(-p)}

#' and s() is a spatial triggering kernel. The public spatio-temporal fitting routines currently use a Gaussian spatial triggering kernel with a homogeneous spatial background over the supplied area and no spatial edge correction. s() can also be omitted to fit the temporal ETAS model without a spatial component.

#' The summation is over all previous earthquakes that occurred in the region, with the i'th such earthquake occurring at time \eqn{ t_i}{t[i]} at spatial location (\eqn{ x_i}{x[i]},,\eqn{ y_i}{y[i]}) aand having magnitude \eqn{ m_i}{m[i]}. The quantity \eqn{M_0}{M0} denotes the magnitude of completeness of the catalog, so that \eqn{m_i \geq  M_0}{m[i] \ge M0} for all i. The temporal ETAS model has 5 parameters: \eqn{\mu} controls the background rate of seismicity, \eqn{ K} and \eqn{ \alpha} determine the productivity (average number of aftershocks) of an earthquake with magnitude \eqn{m}, and \eqn{c} and \eqn{ p}  are the parameters of the Modified Omori Law (which has here been normalized to integrate to 1) and represent the speed at which the aftershock rate decays over time. Each earthquake is assumed to have a magnitude which is an independent draw from the Gutenberg-Richter law \eqn{ p(m_i) = \beta e^{\beta(m_i-M_0)}}{ p(m) = \beta * exp(\beta(m-M0)}.
#'
#' In the spatio-temporal version of the model, s() is the multivariate Gaussian with diagonal covariance parameters \eqn{\sigma^2_x} and \eqn{\sigma^2_y}. The spatial background is homogeneous over the supplied area and no spatial edge correction is implemented.

#' This function estimates the parameters of the ETAS model (both temporal and spatio-temporal) using maximum likelihood. Parameters are optimized on constrained transformed scales from multiple starting values, but the result is still a numerical maximum of a non-convex likelihood.

#' @param ts Vector containing the earthquake times
#' @param ms Vector containing the earthquake magnitudes
#' @param maxTime Length of the time window [0,maxTime] the catalog was observed over. If not specified, will be taken as the time of the last earthquake.
#' @param M0 Magnitude of completeness.
#' @param xs (For the spatio-temporal ETAS model only), a vector containing the x coordinates of the earthquakes
#' @param ys (For the spatio-temporal ETAS model only), a vector containing the y coordinates of the earthquakes
#' @param area (For the spatio-temporal ETAS model only), a single number containing the area of the rectangular spatial region containing the earthquakes. For example if the earthquakes are defined on [28,35]x[-119,-128] then the area would be 7x9 = 63. If this is not specified then it will be calculated based on the distance between the furthest apart earthquakes in the catalog.
#' @param initval Initial value at which to start the estimation. A vector, with elements (mu, K, alpha, c, p) for temporal fitting or (mu, K, alpha, c, p, sigma2x, sigma2y) for spatio-temporal fitting.
#' @param displayOutput If TRUE, display likelihood values during model fitting.
#' @param zeromu If the ETAS model is fit to a single aftershock sequence, we may want to force mu=0 in the model to have no background rate. Set zeromu=TRUE and the ETAS model will be fit under the assumption that mu=0
#' @param handle_ties How to handle duplicate event times. The default "error" stops because continuous-time ETAS requires strictly increasing event times. "jitter" applies a deterministic, order-preserving jitter before fitting.
#' @param c_lower Lower bound for the Omori c parameter. If NULL, a data-scale default equal to one tenth of the smallest positive event-time gap after tie handling is used.
#' @return A list consisting of
#'   \item{params}{A vector containing the estimated parameters, in the order (mu,K,alpha,c,p) for the temporal model or (mu,K,alpha,c,p,sigma2x,sigma2y) for the spatio-temporal model}
#'   \item{loglik}{The corresponding loglikelihood}

#' @examples
#'\donttest{
#' #temporal model
#' ts <- c(1, 2, 3, 5, 8)
#' ms <- c(3.0, 3.1, 3.3, 3.2, 3.0)
#' maxLikelihoodETAS(ts, ms, M0=3, maxTime=10)
#'
#' #spatio-temporal model
#' xs <- c(1, 2, 3, 4, 5)
#' ys <- c(1, 1.5, 2, 1.5, 1)
#' maxLikelihoodETAS(ts, ms, M0=3, maxTime=10, xs=xs, ys=ys, area=25)
#' }
#' @author Gordon J Ross
#' @references Gordon J. Ross - Bayesian Estimation of the ETAS Model for Earthquake Occurrences (2019)
#' @export
maxLikelihoodETAS <- function(ts,ms,M0,maxTime=NULL,xs=NULL,ys=NULL,area=NULL,initval=NA,displayOutput=FALSE,zeromu=FALSE, handle_ties=c("error","jitter"), c_lower=NULL) {
  if (is.null(maxTime)) {maxTime <- max(ts)}
  prep <- .prepare_catalog(ts, ms, M0, maxTime, handle_ties, c_lower)
  ts <- prep$ts
  c_lower <- prep$c_lower

  if (is.null(xs) && is.null(ys)) {
    out <- maxLikelihoodTemporalETAS(ts, ms, maxTime, M0,initval=initval,displayOutput=displayOutput,zeromu=zeromu,c_lower=c_lower)
    return(.attach_fit_attrs(out, c_lower, prep$tie_info))
  } else {
    if (is.null(ys) || is.null(xs)) {
      stop("if fitting spatial ETAS then both xs and ys must be specified", call. = FALSE)
    }
    if (is.null(area)) { area <- (max(xs)-min(xs)) * (max(ys)-min(ys))}
    out <- maxLikelihoodETASspatial(ts,ms,xs,ys,maxTime,M0,area,initval=initval,displayOutput=displayOutput,zeromu=zeromu,c_lower=c_lower)
    return(.attach_fit_attrs(out, c_lower, prep$tie_info))
  }
}





#'Simulates synthetic data from the ETAS model
#'

#'The ETAS model is widely used to quantify the degree of seismic activity in a geographical region, and to forecast the occurrence of future mainshocks and aftershocks (Ross 2019). The temporal ETAS model is a point process where the probability of an earthquake occurring at time \eqn{t} depends on the previous seismicity \eqn{H_t}{Ht}, and is defined by the conditional intensity function:

#' \deqn{ \lambda(t|H_t) = \mu + \sum_{t[i] < t}  \kappa(m[i]|K,\alpha) h(t-t[i]|c,p)}{ \lambda(t|Ht) = \mu + \sum \kappa(m[i]|K,\alpha) h(t[i]|c,p)s(x-x[i],y-y[i],}

#' where

#' \deqn{\kappa(m_i|K,\alpha) = Ke^{\alpha \left( m_i-M_0 \right)}}{\kappa(m[i]|K,\alpha) = K * exp(\alpha(m[i]-M0))}

#' and

#' \deqn{ h(t_i|c,p) = \frac{(p-1)c^{p-1}}{(t-t_i+c)^p}}{ h(t[i]|c,p) = (p-1) * c^(p-1) * (t-t[i]+c)^(-p)}

#' and s() is a spatial triggering kernel. The public spatio-temporal routines currently use a Gaussian spatial triggering kernel with a homogeneous spatial background over the supplied area and no spatial edge correction. s() can also be omitted to fit the temporal ETAS model without a spatial component.

#' The summation is over all previous earthquakes that occurred in the region, with the i'th such earthquake occurring at time \eqn{ t_i}{t[i]} at spatial location (\eqn{ x_i}{x[i]},,\eqn{ y_i}{y[i]}) aand having magnitude \eqn{ m_i}{m[i]}. The quantity \eqn{M_0}{M0} denotes the magnitude of completeness of the catalog, so that \eqn{m_i \geq  M_0}{m[i] \ge M0} for all i. The temporal ETAS model has 5 parameters: \eqn{\mu} controls the background rate of seismicity, \eqn{ K} and \eqn{ \alpha} determine the productivity (average number of aftershocks) of an earthquake with magnitude \eqn{m}, and \eqn{c} and \eqn{ p}  are the parameters of the Modified Omori Law (which has here been normalized to integrate to 1) and represent the speed at which the aftershock rate decays over time. Each earthquake is assumed to have a magnitude which is an independent draw from the Gutenberg-Richter law \eqn{ p(m_i) = \beta e^{\beta(m_i-M_0)}}{ p(m) = \beta * exp(\beta(m-M0)}.
#'
#' In the spatio-temporal version of the model, s() is the multivariate Gaussian with diagonal covariance parameters \eqn{\sigma^2_x} and \eqn{\sigma^2_y}. Spatial simulation uses a bounded-window discard convention: proposed spatial offspring outside the requested window are discarded and discarded offspring do not generate descendants. This is not an edge-corrected simulation of an infinite-plane ETAS process observed through a bounded window.


#' This function simulates sample data from the ETAS model (either temporal or spatio-temporal) over a particular time interval [0,T], For the spatio-temporal model we also specify the range of the x axis region [xmin,xmax] and the y axis region [ymin,ymax]

#' @param mu Parameter of the ETAS model as described above.
#' @param K Parameter of the ETAS model as described above.
#' @param alpha Parameter of the ETAS model as described above.
#' @param c Parameter of the ETAS model as described above.
#' @param p Parameter of the ETAS model as described above.
#' @param beta Parameter of the Gutenberg-Richter law used to generate earthquake magnitudes.
#' @param maxTime Length of the time window [0,maxTime] the catalog was observed over. If not specified, will be taken as the time of the last earthquake.
#' @param M0 Magnitude of completeness.
#' @param sigma2x (For the spatio-temporal ETAS model only) Parameter of the spatio-temporal model as described above.
#' @param sigma2y (For the spatio-temporal ETAS model only) Parameter of the spatio-temporal model as described above.
#' @param xmin (For the spatio-temporal ETAS model only) Minimum x value of the spatial region
#' @param xmax (For the spatio-temporal ETAS model only) Maximum x value of the spatial region
#' @param ymin (For the spatio-temporal ETAS model only) Minimum y value of the spatial region
#' @param ymax (For the spatio-temporal ETAS model only) Maximum y value of the spatial region
#' @return A list consisting of
#'   \item{ts}{The simulated earthquake times}
#'   \item{ms}{The simulated earthquake magnitudes}
#'  \item{locs}{ (For the spatio-temporal ETAS model only) A two column matrix wher the i'th row is the (x,y) location of earthquake i}
#' @examples
#' set.seed(1)
#'
#' #temporal model
#' beta <- 2.4; M0 <- 3; maxTime <- 20
#' catalog <- simulateETAS(0.1, 0.05, 1.0, 0.5, 2, beta, M0, maxTime)
#'
#' #spatio-temporal model
#' sigma2x <- 0.2; sigma2y <- 0.2
#' xmin <- 0; xmax <- 10; ymin <- 0; ymax <- 10
#' catalog <- simulateETAS(0.1, 0.05, 1.0, 0.5, 2, beta, M0, maxTime,
#'    sigma2x, sigma2y, xmin, xmax, ymin, ymax)
#' @author Gordon J Ross
#' @references Gordon J. Ross - Bayesian Estimation of the ETAS Model for Earthquake Occurrences (2019)
#' @export

simulateETAS <-  function(mu,K,alpha,c,p,beta,M0,maxTime,sigma2x=NULL,sigma2y=NULL,xmin=NULL,xmax=NULL,ymin=NULL,ymax=NULL) {
  .check_model_parameters(mu, K, alpha, c, p)
  .check_positive_scalar(beta, "beta")
  .check_numeric_scalar(M0, "M0")
  .check_positive_scalar(maxTime, "maxTime")
  if (all(is.null(sigma2x),is.null(sigma2y),is.null(xmin),is.null(xmax),is.null(ymin),is.null(ymax))) {
    temp <- simulateTemporalETAS(mu,K,alpha,c,p,beta,maxTime,M0)
    temp$branching <- NULL
    return(temp)
  } else if (all(!is.null(sigma2x),!is.null(sigma2y),!is.null(xmin),!is.null(xmax),!is.null(ymin),!is.null(ymax))) {
    .check_positive_scalar(sigma2x, "sigma2x")
    .check_positive_scalar(sigma2y, "sigma2y")
    .check_numeric_scalar(xmin, "xmin")
    .check_numeric_scalar(xmax, "xmax")
    .check_numeric_scalar(ymin, "ymin")
    .check_numeric_scalar(ymax, "ymax")
    if (xmax <= xmin) {stop("xmax must be greater than xmin", call. = FALSE)}
    if (ymax <= ymin) {stop("ymax must be greater than ymin", call. = FALSE)}
    return(simulateSpatialETAS(mu,K,alpha,c,p,beta,maxTime,M0,sigma2x,sigma2y,xmin,xmax,ymin,ymax))
  } else {
    stop("either specify all of sigma2x, sigma2y, xmin, xmax, ymin, ymax to simulate from the spatial ETAS model, or specify none of them to simulate from the temporal ETAS model", call. = FALSE)
  }
}



simulateTemporalETAS <- function(mu,K,alpha,cp,p,beta,maxTime,M0) {
  if (beta <= alpha || K*beta / (beta-alpha) >= 1) {
    warning("parameter choices may lead to supercriticality. Infinite events may be generated. For a guaranteed finite catalog, require beta > alpha and K * beta / (beta - alpha) < 1", call. = FALSE)
  }

  h <- function(z) {(p-1)*cp^(p-1) * 1/( (z+cp)^p )}
  kappa  <- function(m) {
    res <- length(m); res[m<M0] <- 0; inds <- which(m>=M0)
    if (length(inds)>0) {res[inds] <- K * exp(alpha*(m[inds]-M0));}
    res
  }

  maxh <- max(h(seq(0,maxTime,length=10000)))
  rmarks <- function(n) {rexp(n,beta) + M0}

  res <- simulateHawkes(mu,kappa,h,maxTime,rmarks=rmarks,previousts=numeric(),previousmarks=numeric(),starttime=0,useCluster=TRUE)
  res$ms <- res$marks
  res$marks <- NULL
  return(res)
}

simulateSpatialETAS <-  function(mu,K,alpha,cp,p,beta,maxTime,M0,sigma2x,sigma2y,xmin,xmax,ymin,ymax) {

  #first simulate from a temporal ETAS
  pts <- simulateTemporalETAS(mu,K,alpha,cp,p,beta,maxTime,M0)

  #for each point, find a spatial location
  rimmigrant <- function() {return(c(runif(1,xmin,xmax),runif(1,ymin,ymax)))}
  rg <- function(n) {mvrnorm(n,mu=c(0,0),Sigma=matrix(c(sigma2x,0,0,sigma2y),ncol=2))}
  n <- length(pts$ts)
  branching <- pts$branching
  locs <- matrix(numeric(2*n), ncol=2)
  kept <- rep(TRUE, n)
  discarded <- 0
  for (i in seq_len(n)) {
    if (branching[i]==0) {
      locs[i,] <- rimmigrant()
    } else if (!kept[branching[i]]) {
      kept[i] <- FALSE
    } else {
      locs[i,] <- rg(1) + locs[branching[i],]
      if (locs[i,1] < xmin || locs[i,1] > xmax || locs[i,2] < ymin || locs[i,2] > ymax) {
        kept[i] <- FALSE
        discarded <- discarded + 1
      }
    }
  }

  old_to_new <- integer(n)
  old_to_new[kept] <- seq_len(sum(kept))
  branching <- branching[kept]
  positive_parents <- branching > 0
  branching[positive_parents] <- old_to_new[branching[positive_parents]]
  locs <- locs[kept,,drop=FALSE]
  locs <- as.data.frame(locs); colnames(locs) <- c('xs','ys')
  if (discarded > 0) {
    warning(sprintf("Discarded %d proposed offspring outside the spatial window; spatial simulation uses a bounded-window discard convention.", discarded), call. = FALSE)
  }
  return(list(ts=pts$ts[kept], ms=pts$ms[kept], branching=branching, locs=locs))
}










maxLikelihoodTemporalETAS <- function(ts,marks,maxTime,M0,initval=NA,displayOutput=TRUE,zeromu=FALSE,extrasims=0,c_lower=0) {
  .check_catalog(ts, marks, M0, maxTime, require_strict = TRUE)

  initval <- .validate_initval(initval, 5, "initval")
  if (!is.null(initval)) {
    .check_model_parameters(initval[1], initval[2], initval[3], initval[4], initval[5], c_lower)
  }
  starts <- .temporal_start_values(ts, marks, maxTime, initval, c_lower)

  fn <- function(eta) {
    params <- .eta_to_temporal(eta, c_lower, zeromu)
    lik <- ETASlikelihood(ts, marks, maxTime, M0, params[1], params[2], params[3], params[4], params[5], c_lower=c_lower)
    if (!is.finite(lik)) return(1e100)
    if (displayOutput==TRUE) {print(c(params, -lik))}
    -lik
  }

  diagnostics <- data.frame()
  bestpar <- NULL
  bestloglik <- -Inf
  for (i in seq_len(nrow(starts))) {
    start <- starts[i, ]
    start[4] <- max(start[4], c_lower * 2)
    eta0 <- .temporal_to_eta(start, c_lower, zeromu)
    fit <- tryCatch(optim(eta0, fn, method="BFGS", control=list(maxit=1000)),
                    error=function(e) e)
    if (inherits(fit, "error")) {
      diagnostics <- rbind(diagnostics, data.frame(start=i, convergence=NA_integer_, value=NA_real_, loglik=NA_real_, message=conditionMessage(fit), failed=TRUE))
      next
    }
    par <- .eta_to_temporal(fit$par, c_lower, zeromu)
    loglik <- ETASlikelihood(ts, marks, maxTime, M0, par[1], par[2], par[3], par[4], par[5], c_lower=c_lower)
    diagnostics <- rbind(diagnostics, data.frame(start=i, convergence=fit$convergence, value=fit$value, loglik=loglik, message=ifelse(is.null(fit$message), "", fit$message), failed=FALSE))
    if (is.finite(loglik) && loglik > bestloglik) {
      bestloglik <- loglik
      bestpar <- par
    }
  }
  if (is.null(bestpar)) {
    stop("all temporal MLE optimisation starts failed", call. = FALSE)
  }
  if (zeromu==TRUE) {bestpar[1] <- 0}
  names(bestpar) <- c("mu","K","alpha","c","p")
  final_loglik <- ETASlikelihood(ts, marks, maxTime, M0, bestpar[1], bestpar[2], bestpar[3], bestpar[4], bestpar[5], c_lower=c_lower)
  res <- list(params=bestpar, loglik=final_loglik)
  attr(res, "mle_diagnostics") <- list(num_starts=nrow(starts),
                                       starts=starts,
                                       results=diagnostics,
                                       best_loglik=final_loglik,
                                       failed_starts=sum(diagnostics$failed))
  return(res)
}

#extra sims does extra random restarts
maxLikelihoodETASspatial <- function(ts,marks,xs,ys,maxTime,M0,area,initval=NA,displayOutput=TRUE,zeromu=FALSE,extrasims=0,c_lower=0) {
  .check_spatial_catalog(ts, marks, xs, ys, M0, maxTime, area)

  initval <- .validate_initval(initval, 7, "initval")
  if (!is.null(initval)) {
    .check_model_parameters(initval[1], initval[2], initval[3], initval[4], initval[5], c_lower)
    .check_positive_scalar(initval[6], "initval[6]")
    .check_positive_scalar(initval[7], "initval[7]")
  }
  temporal_seed <- tryCatch(maxLikelihoodTemporalETAS(ts, marks, maxTime, M0, displayOutput=FALSE, zeromu=zeromu, c_lower=c_lower)$params,
                            error=function(e) NULL)
  starts <- .spatial_start_values(ts, marks, xs, ys, maxTime, initval, c_lower, temporal_seed)

  fn <- function(eta) {
    params <- .eta_to_spatial(eta, c_lower, zeromu)
    lik <- ETASspatiallikelihood(ts, marks,xs,ys, maxTime, M0, area,params[1], params[2], params[3], params[4], params[5], params[6], params[7], c_lower=c_lower)
    if (!is.finite(lik)) return(1e100)
    if (displayOutput==TRUE) {print(c(params, -lik))}
    -lik
  }

  diagnostics <- data.frame()
  bestpar <- NULL
  bestloglik <- -Inf
  for (i in seq_len(nrow(starts))) {
    start <- starts[i, ]
    start[4] <- max(start[4], c_lower * 2)
    eta0 <- .spatial_to_eta(start, c_lower, zeromu)
    fit <- tryCatch(optim(eta0, fn, method="BFGS", control=list(maxit=1000)),
                    error=function(e) e)
    if (inherits(fit, "error")) {
      diagnostics <- rbind(diagnostics, data.frame(start=i, convergence=NA_integer_, value=NA_real_, loglik=NA_real_, message=conditionMessage(fit), failed=TRUE))
      next
    }
    par <- .eta_to_spatial(fit$par, c_lower, zeromu)
    loglik <- ETASspatiallikelihood(ts, marks,xs,ys, maxTime, M0, area,par[1], par[2], par[3], par[4], par[5], par[6], par[7], c_lower=c_lower)
    diagnostics <- rbind(diagnostics, data.frame(start=i, convergence=fit$convergence, value=fit$value, loglik=loglik, message=ifelse(is.null(fit$message), "", fit$message), failed=FALSE))
    if (is.finite(loglik) && loglik > bestloglik) {
      bestloglik <- loglik
      bestpar <- par
    }
  }
  if (is.null(bestpar)) {
    stop("all spatio-temporal MLE optimisation starts failed", call. = FALSE)
  }
  if (zeromu==TRUE) {bestpar[1] <- 0}
  names(bestpar) <- c("mu","K","alpha","c","p","sigma2x","sigma2y")
  final_loglik <- ETASspatiallikelihood(ts, marks,xs,ys, maxTime, M0, area,bestpar[1], bestpar[2], bestpar[3], bestpar[4], bestpar[5], bestpar[6], bestpar[7], c_lower=c_lower)
  res <- list(params=bestpar, loglik=final_loglik)
  attr(res, "mle_diagnostics") <- list(num_starts=nrow(starts),
                                       starts=starts,
                                       results=diagnostics,
                                       best_loglik=final_loglik,
                                       failed_starts=sum(diagnostics$failed))
  return(res)
}



#check the C version of etas likelihood
ETASlikelihood <- function(ts, marks, maxTime, M0, mu, K, alpha, cp, p, c_lower=0) {
  if (mu < 0 || K < 0 || alpha < 0 || cp <= c_lower || cp >= .ETAS_C_UPPER || p <= 1 || p >= .ETAS_P_UPPER) {return(-Inf)}
  if (length(ts) != length(marks)) {stop("marks should be same length as times", call. = FALSE)}

  res <- .C("ETASlikelihoodC", ts=as.double(ts), marks=as.double(marks), n=as.integer(length(ts)), maxTime = as.double(maxTime), M0=as.double(M0), mu=as.double(mu), K=as.double(K), alpha=as.double(alpha), c=as.double(cp), p=as.double(p),loglik=as.double(0))
  return(res$loglik)
}

#check the C version of etas likelihood
ETASspatiallikelihood <- function(ts, marks, xs,ys, maxTime, M0,area,mu, K, alpha, cp, p,param1,param2,c_lower=0) {

  if (mu < 0 || K < 0 || alpha < 0 || cp <= c_lower || cp >= .ETAS_C_UPPER || p <= 1 || p >= .ETAS_P_UPPER || param1 <= 0 || param2 <= 0) {return(-Inf)}

  if (length(ts) != length(marks)) {stop("marks should be same length as times", call. = FALSE)}
  if (length(ts) != length(xs)) {stop("xs should be same length as times", call. = FALSE)}
  if (length(ts) != length(ys)) {stop("ys should be same length as times", call. = FALSE)}

  res <- .C("ETASspatiallikelihoodC", ts=as.double(ts), marks=as.double(marks), xs=as.double(xs),ys=as.double(ys),n=as.integer(length(ts)), maxTime = as.double(maxTime), M0=as.double(M0), area=as.double(area),mu=as.double(mu), K=as.double(K), alpha=as.double(alpha), c=as.double(cp), p=as.double(p),param1=as.double(param1),param2=as.double(param2),loglik=as.double(0))
  return(res$loglik)
}



estimateTemporalETAS <- function(ts, marks, maxTime, M0, numMCMCSamples=500,sims=5000, burnin=1000, initval=NA, c_lower=0) {
  .check_catalog(ts, marks, M0, maxTime, require_strict = TRUE)
  .check_burnin_sims(sims, burnin)
  initval <- .validate_initval(initval, 5, "initval")
  if (is.null(initval)) {
    initval <- maxLikelihoodTemporalETAS(ts,marks,maxTime,M0,c_lower=c_lower,displayOutput=FALSE)$params
  }
  .check_model_parameters(initval[1], initval[2], initval[3], initval[4], initval[5], c_lower)

  if (length(ts) != length(marks)) {stop("ms should be same length as ts", call. = FALSE)}

  requested_sims <- sims
  sims <- sims+burnin
  branching <- rep(0,length(ts)); mu <- initval[1]; K <- initval[2]; alpha <- initval[3]; cp <- initval[4]; p <- initval[5]

  res <- .C("estimateETASBranchingC", ts=as.double(ts), marks=as.double(marks), branching=as.integer(branching), n=as.integer(length(ts)), maxTime = as.double(maxTime), M0=as.double(M0), sims=as.integer(sims), numMCMCSamples=as.integer(numMCMCSamples), mu=as.double(mu), K=as.double(K), alpha=as.double(alpha), c=as.double(cp), p=as.double(p), c_lower=as.double(c_lower), mus=as.double(numeric(sims)), Ks=as.double(numeric(sims)), alphas=as.double(numeric(sims)), cs=as.double(numeric(sims)), ps=as.double(numeric(sims)), cp_accept=as.integer(0), cp_total=as.integer(0), kappa_accept=as.integer(0), kappa_total=as.integer(0));

  bayesbranching <- cbind(res$mus,res$Ks,res$alphas, res$cs, res$ps)
  bayesbranching <- .drop_burnin(bayesbranching, burnin)
  if (nrow(bayesbranching) != requested_sims) {stop("internal error: posterior sample count mismatch", call. = FALSE)}
  bayesbranching <- as.data.frame(bayesbranching)
  colnames(bayesbranching) <- c('mu','K','alpha','c','p')
  attr(bayesbranching, "acceptance_rates") <- list(cp=res$cp_accept/res$cp_total,
                                                   K_alpha=res$kappa_accept/res$kappa_total,
                                                   cp_accept=res$cp_accept,
                                                   cp_total=res$cp_total,
                                                   K_alpha_accept=res$kappa_accept,
                                                   K_alpha_total=res$kappa_total)
  return(bayesbranching)
}

estimateSpatialETAS<- function(ts, marks, xs, ys, maxTime, M0,area,numMCMCSamples=500,sims=5000, burnin=1000, initval=NA, gaussian=TRUE, sigma2_prior_shape=3, sigma2x_prior_rate=NULL, sigma2y_prior_rate=NULL, c_lower=0) {
  .check_spatial_catalog(ts, marks, xs, ys, M0, maxTime, area)
  .check_burnin_sims(sims, burnin)
  prior <- .resolve_sigma2_prior(xs, ys, sigma2_prior_shape, sigma2x_prior_rate, sigma2y_prior_rate)
  initval <- .validate_initval(initval, 7, "initval")
  if (is.null(initval)) {
    initval <- maxLikelihoodETASspatial(ts,marks,xs,ys,maxTime,M0,area,displayOutput = FALSE,c_lower=c_lower)$params
  }
  .check_model_parameters(initval[1], initval[2], initval[3], initval[4], initval[5], c_lower)
  .check_positive_scalar(initval[6], "initval[6]")
  .check_positive_scalar(initval[7], "initval[7]")

  if (length(ts) != length(marks)) {stop("ms should be same length as ts", call. = FALSE)}
  if (length(ts) != length(xs)) {stop("xs should be same length as ts", call. = FALSE)}
  if (length(ts) != length(ys)) {stop("ys should be same length as ts", call. = FALSE)}

  requested_sims <- sims
  sims <- sims+burnin
  branching <- rep(0,length(ts));

  mu <- initval[1]; K <- initval[2]; alpha <- initval[3]; cp <- initval[4]; p <- initval[5]; param1 <- initval[6]; param2 <- initval[7]



  res <- .C("estimateSpatialETASBranchingC", ts=as.double(ts), marks=as.double(marks), xs=as.double(xs), ys=as.double(ys), branching=as.integer(branching), n=as.integer(length(ts)), maxTime = as.double(maxTime), M0=as.double(M0), area=as.double(area),sims=as.integer(sims), numMCMCSamples=as.integer(numMCMCSamples), mu=as.double(mu), K=as.double(K), alpha=as.double(alpha), c=as.double(cp), p=as.double(p), c_lower=as.double(c_lower), param1=as.double(param1), param2=as.double(param2), sigma2_prior_shape=as.double(prior$shape), sigma2x_prior_rate=as.double(prior$rate_x), sigma2y_prior_rate=as.double(prior$rate_y), mus=as.double(numeric(sims)), Ks=as.double(numeric(sims)), alphas=as.double(numeric(sims)), cs=as.double(numeric(sims)), ps=as.double(numeric(sims)), param1s=as.double(numeric(sims)), param2s=as.double(numeric(sims)), cp_accept=as.integer(0), cp_total=as.integer(0), kappa_accept=as.integer(0), kappa_total=as.integer(0), gaussian=as.logical(gaussian))

  bayesbranching <- cbind(res$mus,res$Ks,res$alphas, res$cs, res$ps, res$param1s, res$param2s)
  bayesbranching <- .drop_burnin(bayesbranching, burnin)
  if (nrow(bayesbranching) != requested_sims) {stop("internal error: posterior sample count mismatch", call. = FALSE)}
  bayesbranching <- as.data.frame(bayesbranching)
  colnames(bayesbranching) <- c('mu','K','alpha','c','p','sigma2x','sigma2y')
  attr(bayesbranching, "sigma2_prior") <- prior
  attr(bayesbranching, "acceptance_rates") <- list(cp=res$cp_accept/res$cp_total,
                                                   K_alpha=res$kappa_accept/res$kappa_total,
                                                   cp_accept=res$cp_accept,
                                                   cp_total=res$cp_total,
                                                   K_alpha_accept=res$kappa_accept,
                                                   K_alpha_total=res$kappa_total,
                                                   sigma2="Gibbs")
  return(bayesbranching)
}





#old hawkes functions, for simulation


#simulates observatoins from on an poisson process with intensity function 'targetfn' using thinning.
#if T is set, the simulatoin is on [0,T]. Alternatively if 'n' is set, the simulation lasts until n events have been generated (with T random)
sampleNHPP <- function(targetfn, maxintensity, starttime=0,T=Inf, n=Inf) {
  if (!is.finite(maxintensity) || maxintensity <= 0) {
    return(numeric())
  }
  pts <- numeric()
  s <- starttime
  while (length(pts) < n) {
    s <- s + rexp(1,maxintensity)
    if (s > T) {break}
    if (runif(1) < targetfn(s)/maxintensity) {
      pts <- c(pts,s)
    }
  }
  return(pts)
}

#simulates from a Hawkes process using either thinning or the cluster process representation (the choice of which to use is controlled by the "useCluster" function arguments, with the default being thinning). The process can be either marked or unmarked.

#Note: at the moment, the thinning method only works if the kernel is monotonic decreasing. In this case, simulation is very fast. However for (eg) lognormal kernels that are not monotonic decreasing, it cannot be used. As such, the default is the cluster process method.

#mu = background rate. Can be constant, or a function. If a function then mu(t) should be mu evaluated at the given time point
#kappa = hawkes parameter. A constant, or a function (usually a function if the process is marked). If a function, then kappa(marks) should be kappa evaluated at the given marks value (vectorised)
#h = hawkes kernal Should be a first order function
#maxTime = the time to simulate the process until
#previousts/prevousmarks/starttime = sometimes we will want to consider the influence of historical events prior to the start of the simulation period. If so, these events should be passed in previousts, with starttime being the time at which the simulation starts. For a marked process, previousmarks is the marks.
#maxmu = maximum value of the background mu when mu is a function. If this argument is missing, the maximum will be estimated.
#maxh = maximum value of the kernel h, only used in cluster process simulation. If this argument is missing, the maximum will be estimated.
#rmarks = a function that returns n observations from the mark distribtuion
#useCluster = if TRUE then we use the cluster process simulatoin method. If FALSE then we use thinning. Note that the cluster process is typically faster, but is unable to take previous historical events into account. As such it cant be used if previousts/previousmarks/starttime is set.
simulateHawkes <-function(mu, kappa, h, maxTime, previousts=numeric(), previousmarks=numeric(), starttime=0, maxmu=NA, maxh=NA, rmarks=NA, useCluster=TRUE) {
  ###
  # check/validate arguments
  ###

  #decide whether to use cluster process of thinning for the simulation
  if (useCluster==TRUE) {
    if (length(previousts)>0 || length(previousmarks)>0 || starttime>0) {
      warning("cannot use cluster process sampling with historical events; switching to thinning instead", call. = FALSE)
      useCluster <- FALSE
    }
  }

  #if mu is a function, we need the maximum value of its intensity for thinning
  if (is.function(mu) && is.na(maxmu)) {
    maxmu <- max(mu(seq(starttime,maxTime,length=1000)))
  }

  if (useCluster==TRUE && is.na(maxh)) {
    maxh <- max(h(seq(0,maxTime,length=1000)))
  }

  #if mu is constant, turn it into a function to make the later code more consistent
  tempmu <- mu
  if (!is.function(mu)) {
    mu <- function(x) {tempmu}
    maxmu <- tempmu
  }

  if (length(previousts)>0 && max(previousts) > starttime) {
    stop("attempting to start simulation earlier than last historical event", call. = FALSE)
  }

  if (length(previousts) != length(previousmarks)) {
    stop("historical times and marks have different lengths", call. = FALSE)
  }

  ###
  # Now call the relevant function
  ###
  if (useCluster==FALSE && !is.function(rmarks)) {
    return(simulateHawkesThinningUnmarked(mu,kappa,h,maxTime,maxmu,previousts=previousts,starttime=starttime))
  } else if (useCluster==FALSE && is.function(rmarks)) {
    return(simulateHawkesThinningMarked(mu,kappa,h,maxTime,maxmu,rmarks,previousts=previousts,previousmarks=previousmarks,starttime=starttime))
  } else if (useCluster==TRUE && !is.function(rmarks)) {
    return(simulateHawkesClusterUnmarked(mu,kappa,h,maxTime,maxmu,maxh))
  } else if (useCluster==TRUE && is.function(rmarks)) {
    return(simulateHawkesClusterMarked(mu,kappa,h,maxTime,maxmu,maxh,rmarks))
  } else {
    stop("invalid arguments", call. = FALSE)
  }
}

#simulate hawkes process with no marks. Shouldnt be called directly, use the simulateHawkes() function to set up arguments.
simulateHawkesThinningUnmarked <- function(mu, kappa, h, maxTime, maxmu, previousts=numeric(),starttime=0) {
  ts <- numeric(500000) #list of events
  numprevious <- length(previousts)
  index <- 0
  if (numprevious > 0) {
    ts[1:numprevious] <- previousts
    index <- numprevious
  }

  t <- starttime #current time

  while(t < maxTime) {
    #choose homogenous poisson intensity that is greater than max hawkes value
    M <- ifelse(index==0, maxmu, maxmu + sum(kappa*h(t-ts[1:index])))
    t <- t + rexp(1,M)

    hawkeseval <- mu(t) + sum(kappa*h(t-ts[1:index]))
    if (runif(1) < hawkeseval/M) {
      index <- index+1
      ts[index] <- t
    }
  }
  ts <- ts[1:index]
  ts <- ts[ts>starttime & ts<=maxTime]
  return(list(ts=ts))
}

#simulate hawkes process with marks.  Shouldnt be called directly, use the simulateHawkes() function to set up arguments.
simulateHawkesThinningMarked <- function(mu, kappa,h, maxTime,maxmu,rmarks,previousts=numeric(),previousmarks=numeric(),starttime=0) {

  ts <- numeric(500000) #list of events
  marks <- numeric(500000) #list of marks
  numprevious <- length(previousts)
  index <- 0
  if (numprevious > 0) {
    ts[1:numprevious] <- previousts
    marks[1:numprevious] <- previousmarks
    index <- numprevious
  }

  t <- starttime #current time

  while(t < maxTime) {
    #choose homogenous poisson intensity that is greater than max hawkes value
    M <- ifelse(index==0, maxmu, maxmu + sum(kappa(marks[1:index])*h(t-ts[1:index])))
    t <- t + rexp(1,M)
    hawkeseval <- mu(t) + sum(kappa(marks[1:index])*h(t-ts[1:index]))
    if (runif(1) < hawkeseval/M) {
      index <- index+1
      ts[index] <- t
      marks[index] <- rmarks(1)
    }
  }
  ts <- ts[1:index]; marks <- marks[1:index]

  inds <- which(ts>starttime & ts<=maxTime)
  ts <- ts[inds]; marks <- marks[inds]
  return(list(ts=ts,marks=marks))
}


#simulation of hawkes using cluster process
simulateHawkesClusterUnmarked <- function(mu,kappa,h,maxTime,maxmu,maxh) {
  x <- sampleNHPP(mu,maxmu,T=maxTime) #background events
  if (length(x)==0) {return(list(ts=numeric(), branching=numeric()))}
  count <- 1
  caused <- rep(0,length(x)) #contains the index of the point which caused this one, or 0 if it is an immigrant
  while(TRUE) { #add new events onto x...
    if (count > length(x)) {break}
    pt <- x[count]

    fn <- function(t) {kappa*h(t)}

    pts <- sampleNHPP(fn, maxintensity=kappa*maxh, T=maxTime-pt)
    if (!is.null(pts)) {
      x <- c(x, pts+pt)
      caused <- c(caused, rep(pt,length(pts)))
    }
    count <- count+1
  }

  ord <- order(x)
  x <- x[ord]; caused <- caused[ord]

  if (length(caused) > 0) {
    for (i in 1:length(caused)) {
      if (caused[i]==0) {next}
      caused[i] <- which(x==caused[i])
    }
  }
  return(list(ts=x,branching=caused))
}



#simulation of hawkes using cluster process
simulateHawkesClusterMarked <- function(mu,kappa,h,maxTime,maxmu,maxh,rmarks) {
  x <- sampleNHPP(mu,maxmu,T=maxTime)
  if (length(x)==0) {return(list(ts=numeric(), marks=numeric(),branching=numeric()))}
  m <- rmarks(length(x))

  count <- 1
  caused <- rep(0,length(x))

  while(TRUE) { #add new events onto x...
    if (count > length(x)) {break}
    pt <- x[count]

    fn <- function(t) {kappa(m[count])*h(t)}
    maxintensity <- kappa(m[count])*maxh

    pts <- sampleNHPP(fn, maxintensity=maxintensity, T=maxTime-pt)
    if (!is.null(pts)) {
      x <- c(x, pts+pt)
      m <- c(m,rmarks(length(pts)))
      caused <- c(caused, rep(pt,length(pts)))
    }
    count <- count+1
  }

  ord <- order(x)
  x <- x[ord]; m <- m[ord]; caused <- caused[ord]
  for (i in 1:length(caused)) {
    if (caused[i]==0) {next}
    caused[i] <- which(x==caused[i])
  }

  return(list(ts=x,marks=m,branching=caused))
}