#ETAS simulator by Igor Szumny 
#ETAS by A. Jalilian in use (https://CRAN.R-project.org/package=ETAS )

library(MASS)
library(ETAS)
library(spatstat)
library(gstat)
library(sf)
library(sp)

#Estimation of mu density
estimate_mu_density <- function(historical, region, n_grid = 200) {
  kde <- MASS::kde2d(
    x = historical$long,
    y = historical$lat,
    n = n_grid,
    lims = c(region$long[1], region$long[2],
             region$lat[1], region$lat[2])
  )
  kde$z <- kde$z / sum(kde$z) #normalization
  return(kde)
}

#sampling mu
sample_from_mu <- function(kde, n) {
  probs <- as.vector(kde$z)
  idx <- sample(seq_along(probs), size = n, replace = TRUE, prob = probs)
  
  nx <- length(kde$x)
  ny <- length(kde$y)
  
  ix <- ((idx - 1) %% nx) + 1
  iy <- ((idx - 1) %/% nx) + 1
  
  long <- kde$x[ix]
  lat  <- kde$y[iy]
  
  data.frame(lat, long)
}

#simulate catalog
simulate_etas_catalog <- function(
    T = 365, # time [days]
    region = list(lat = c(34, 36), long = c(-118, -116)), # recktangle (latmin, latmax), (longmin, longmax)
    mu = 0.05, K = 0.8, alpha = 1.0,
    c = 0.01, p = 1.1,
    D = 1.0, q = 1.5, gamma = 0.0,
    M0 = 2.5, mag_b = 1.0,
    dist_unit = "degree", #km OR degree
    seed = NULL,
    mu_kde = NULL,   # mu density, if NULL density is calculate by runif()
    time_start = "2005-05-11 00:00:00"
) {
  
  if(!is.null(seed)) set.seed(seed)
  
  # magnitude Gutenberg-Richter
  rmagnitud <- function(n, M0, b = mag_b) {
    u <- runif(n)
    mags <- M0 - (1 / b) * log(1 - u)
    return(mags)
  }
  
  #offspring
  generate_offspring <- function(parent_time, parent_lat, parent_long, parent_mag) {
 
    # k(m)
    lambda_off <- K * exp(alpha * (parent_mag - M0))
    if (isTRUE(!lambda_off)) return(NULL)
    n_off <- rpois(1, lambda_off)
    if (isTRUE(!n_off)) return(NULL)
    #g(t)
    u <- runif(n_off)
    if(p == 1) {
      t_off <- c * (exp(u * log(1 + (T * 10) / c)) - 1)
    } else {
      t_off <- c * ((1 - u)^(1/(1 - p)) - 1)
      t_off <- pmax(0, t_off)
    }
    
    #f(x)
    u2 <- runif(n_off)
    D_m <- D * exp(gamma * (parent_mag - M0))
    r <-  sqrt(D_m * ( (1 - u2)^(1/(1 - q)) - 1 ) )
    theta <- runif(n_off, 0, 2*pi)
    r_deg <- if(dist_unit == "km") r / 111.0 else r #NOT TESTED
    lat_off <- parent_lat + r_deg * sin(theta)
    long_off <- parent_long + r_deg * cos(theta) / cos(parent_lat * pi/180)
    
    #Guteberg-Richter
    mags <- rmagnitud(n_off, M0)
    
    #return offspring
    df <- data.frame(
      time = parent_time + t_off,
      lat = lat_off,
      long = long_off,
      mag = mags,
      parent_time = parent_time,
      stringsAsFactors = FALSE
    )
    return(df)
  }
  
  # ------------ MAIN ------------ 
  #mu(x,y,t)
  n_background <- rpois(1, mu*T) #*T
  times_bg <- sort(runif(n_background, 0, T))
  if (!is.null(mu_kde)) {
    #mu(x,y,t)
    bg_points <- sample_from_mu(mu_kde, n_background)
    lats_bg <- bg_points$lat
    longs_bg <- bg_points$long
  } else {
    #mu(t)
    lats_bg <- runif(n_background, region$lat[1], region$lat[2])
    longs_bg <- runif(n_background, region$long[1], region$long[2])
  }
  
  #Gutenerg-Richter
  mags_bg <- rmagnitud(n_background, M0)
  
  catalog_events <- data.frame(
    time = times_bg,
    lat = lats_bg,
    long = longs_bg,
    mag = mags_bg,
    parent_time = NA_real_,
    stringsAsFactors = FALSE
  )
  
  # ------------ offsprings ------------ 
  i <- 1
  while(i <= nrow(catalog_events)) {
    ev <- catalog_events[i, ]
    tmp <- try(offs <- generate_offspring(ev$time, ev$lat, ev$long, ev$mag), silent=TRUE)
    if(!is.null(offs)) {
      
      #time filter
      offs <- offs[offs$time <= T,  c("time", "lat", "long", "mag", "parent_time")]
      
      #latitude and longitude filter
      offs <- offs[
        offs$lat >= region$lat[1] & offs$lat <= region$lat[2] &
          offs$long >= region$long[1] & offs$long <= region$long[2],
      ]
      
      if(nrow(offs) > 0)
        catalog_events <- rbind(catalog_events, offs)
    }
    i <- i + 1
    
    if(nrow(catalog_events) > 1000000) {
      warning("Katalog przekroczył 1e6 zdarzeń — przerwano.")
      break
    }
    
  }
  
  catalog_events <- catalog_events[order(catalog_events$time), ]
  rownames(catalog_events) <- NULL
  
  # ------------ time ------------ 
  t0 <- as.POSIXct(time_start, tz = "UTC")
  event_posix <- t0 + as.difftime(catalog_events$time, units = "days")
  
  out_df <- data.frame(
    date = format(event_posix, "%Y-%m-%d"),
    time = format(event_posix, "%H:%M:%S"),
    lat = catalog_events$lat,
    long = catalog_events$long,
    mag = round(catalog_events$mag, 3),
    stringsAsFactors = FALSE
  )
  attr(out_df, "event_time_days") <- catalog_events$time
  attr(out_df, "parent_time") <- catalog_events$parent_time
  
  # ------------ return ------------   
  if(requireNamespace("ETAS", quietly = TRUE)) {
   #as ETAS catalog
   cat_obj <- ETAS::catalog(out_df,
                            time.begin = as.character(format(t0, "%Y-%m-%d %H:%M:%S")),
                             study.start = as.character(format(t0, "%Y-%m-%d %H:%M:%S")),
                             study.length = T)
     return(list(data = out_df, catalog = cat_obj))
  } else {
  #as list
  return(list(data = out_df, catalog = NULL))
  }
}