library(tidyr)
library(ETAS)
library(tictoc)
#First run synthetic.R

#1. Loading data
df <- read.csv("...\\ETAS\\Italy_2016.csv", header = TRUE, stringsAsFactors = FALSE)

df
df <- df %>% separate(time2, into = c("date", "time"), sep = " ") 

#2.Creating catalog
lat = c( min(df$lat),max(df$lat))#(,27.7)20.5
long = c(min(df$long),max(df$long))#(90.2, 97.2)
start = "2005/08/24"
end = "2016/08/24"
mag = 2.7

italy.cat <- catalog(df, time.begin=start, study.start=start, study.end=end,lat.range=lat, long.range=long, mag.threshold=mag, dist.unit ="degree")
plot(italy.cat)
italy.cat

#3. Fitting
param0 <- c(1.010084e+00, 3.004448e-01, 1.055809e-02, 1.222385e+00, 1.167979e+00, 7.565286e-05, 1.822504e+00, 8.426979e-01) #2016 mag>=2.7result
nthreads <- parallel::detectCores()
nthreads
italy.fit <- etas(italy.cat, param0,  nthreads=nthreads)
italy.fit 

#4. Model validation
resid.etas(italy.fit)
italy.rates <- rates(italy.fit, dimyx=c(100, 125))

#5. Preparing data for simulation
region <- list(lat=lat, long=long)
mu_kde <- estimate_mu_density(df, region) #from synthetic.R

#6. Simulation
res <- simulate_etas_catalog(
  T = diff(italy.cat$rtperiod),
  region = region,
  mu = 2876.534/diff(italy.cat$rtperiod), K = italy.fit$param['A'], alpha = italy.fit$param['alpha'] , c = italy.fit$param['c'] , 
  p = italy.fit$param['p'],  D = italy.fit$param['D'], q = italy.fit$param['q'], gamma=italy.fit$param['gamma'], 
  M0 = mag, 
  mag_b = 2.7052, 
  dist_unit = "degree", 
  mu_kde = mu_kde, 
  time_start = start,
  #seed = 42
)

#7. Synthetic vs. real catalog comparsion
cat("Synthetic:", nrow(res$data), "\n")
cat("Real data:", nrow(italy.cat$longlat.coord), "\n")
cat( (1 / (mean(res$data$mag) - mag)))
cat( max(res$data$mag) )
head(res$data)

#8. Ploting synthetic catalog
plot(res$catalog)

#9. Saving synthetic catalog
write.csv(res$data, file = "...\\ETAS//synthetic_catalog.csv", row.names = FALSE)

#10. Generating 1000 synthetic catalogs
tic("Generowanie")
num_iterations <- 3000
for (i in 1:num_iterations) {
  #6. Simulation
  res <- simulate_etas_catalog(
    T = diff(italy.cat$rtperiod),
    region = region,
    mu = 2876.534/diff(italy.cat$rtperiod), K = italy.fit$param['A'], alpha = italy.fit$param['alpha'] , c = italy.fit$param['c'] , 
    p = italy.fit$param['p'],  D = italy.fit$param['D'], q = italy.fit$param['q'], gamma=italy.fit$param['gamma'], 
    M0 = mag, 
    mag_b = 2.7052, 
    dist_unit = "degree", 
    mu_kde = mu_kde, 
    time_start = start,
    #seed = 42
  )
  #11. Saving synthetic catalogs
  file_name <- paste0("...\\ETAS//1000//italy_merged_", i, ".csv")
  write.csv(res$data, file = file_name, row.names = FALSE)
  if (i %% 100 == 0) {
    cat(paste0("Zapisano plik numer: ", i, " (", file_name, ")\n"))
  }
}
toc()

num_value <- numeric(1000)
beta_value <- numeric(1000)
for (i in 1:num_iterations) {
  file_name1 <- paste0("...\\ETAS\\1000\\italy_merged_", i, ".csv")
  df <- read.csv(file_name1, header = TRUE, stringsAsFactors = FALSE)
  #12. Statistics
  num_value[i] <- nrow(df)
  beta_value[i] <- (1 / (mean(df[['mag']]) - mag))
  #13. Plotting synthetic catalogs
  #file_name <- paste0("...\\plots\\s", i, ".png")
  #png(filename = file_name, width = 732, height = 739)
  #plot(italy.cat)
  #dev.off()
  if (i %% 100 == 0) {
    cat(paste0("Zapisano plik numer: ", i, " (", file_name, ")\n"))
  }
}

#14. Saving statistics
write.table(
  x = beta_value,
  file = "...\\ETAS//beta_value.txt",
  row.names = FALSE,
  col.names = FALSE,
  quote = FALSE,
  sep = "\n"
)
write.table(
  x = mu_kde,
  file = "...\\ETAS\\mu_kde.txt",
  row.names = TRUE,
  col.names = TRUE,
  quote = FALSE,
  sep = ";"
)