#Loading require package.
require(VineCopula)
require(lmom)

#' Calculate LMF and decompose into mean/variability/dependence components
#'
#' @param tas_da 3D array [lon, lat, time] monthly temperature (tas)
#' @param pr_da  3D array [lon, lat, time] monthly precipitation (pr)
#' @param i lon index
#' @param j lat index
#' @param threshold percentile threshold in (0,1), e.g., 0.9
#' @param window_yr moving window length (years), e.g., 30
#' @return list with LMF and components

cal_lmf <- function(tas_da, pr_da, i, j, threshold, window_yr) {
  
  # --- basic checks
  stopifnot(threshold > 0, threshold < 1)
  stopifnot(window_yr >= 2)
  
  nt <- length(tas_da[i, j, ])
  if (nt != length(pr_da[i, j, ])) stop("tas and pr time lengths do not match.")
  if (nt %% 12 != 0) stop("Time length must be a multiple of 12 (monthly data).")
  
  yr <- nt / 12
  
  # --- reshape to [month, year]
  tas_ijk <- matrix(tas_da[i, j, ], nrow = 12, ncol = yr)
  pr_ijk  <- matrix(pr_da [i, j, ], nrow = 12, ncol = yr)
  
  # --- monthly anomalies: remove climatological monthly mean
  tas_ijk <- tas_ijk - rowMeans(tas_ijk, na.rm = TRUE)
  pr_ijk  <- pr_ijk  - rowMeans(pr_ijk , na.rm = TRUE)
  
  # --- rolling windows over years (each window is 12*window_yr months)
  len <- yr - window_yr + 1
  if (len < 2) stop("window_yr too large for the available record length.")
  
  loca_sel <- embed(1:yr, window_yr)  # each row is a window of years (reversed order)
  
  # --- helper for robust sd check
  ok_series <- function(x) {
    x <- as.vector(x)
    x <- x[is.finite(x)]
    length(x) >= 5 && stats::sd(x) > 0
  }
  
  # --- pre-allocate as numeric vectors (faster and cleaner)
  p_fut   <- rep(NA_real_, nrow(loca_sel) - 1) #that is p_fut   
  p_pr_mv <- rep(NA_real_, nrow(loca_sel) - 1) #that is p_fut^(m-v-pr)
  p_ta_mv <- rep(NA_real_, nrow(loca_sel) - 1) #that is p_fut^(m-v-Tas)
  p_pr_v  <- rep(NA_real_, nrow(loca_sel) - 1) #that is p_fut^(v-pr) 
  p_ta_v  <- rep(NA_real_, nrow(loca_sel) - 1) #that is p_fut^(v-Tas)
  p_d2    <- rep(NA_real_, nrow(loca_sel) - 1) #that is p_fut^D  
  p_d1    <- rep(NA_real_, nrow(loca_sel) - 1) #that is p_his   
  # --- main loop: kk=2..nrow(loca_sel); store at idx=kk-1
  
  for (kk in 2:nrow(loca_sel)) {#kk=20
    
    idx <- kk - 1
    
    pr_fut  <- matrix(pr_ijk [, loca_sel[kk, ]], ncol=1)
    tas_fut <- matrix(tas_ijk[, loca_sel[kk, ]], ncol=1)
    pr_his  <- matrix(pr_ijk [, loca_sel[1 , ]], ncol=1)
    tas_his <- matrix(tas_ijk[, loca_sel[1 , ]], ncol=1)
    
    # require usable variance
    if (!ok_series(pr_fut) || !ok_series(pr_his) || !ok_series(tas_fut) || !ok_series(tas_his)) {
      next
    }
    
    # --- require detrend() exists
    if (!exists("detrend", mode = "function")) {
      stop("Function detrend() not found. Please define it or load the package providing it.")
    }
    
    # --- fit Frank copula on pseudo-observations
    m_fut <- cbind(pobs(detrend(pr_fut)),  pobs(detrend(tas_fut)))
    m_his <- cbind(pobs(detrend(pr_his)),  pobs(detrend(tas_his)))
    
    pa_fut <- BiCopEst(m_fut[, 1], m_fut[, 2], family = 5)
    pa_his <- BiCopEst(m_his[, 1], m_his[, 2], family = 5)
    
    cop_fut <- BiCop(family = 5, par = pa_fut$par, par2 = pa_fut$par2)
    cop_his <- BiCop(family = 5, par = pa_his$par, par2 = pa_his$par2)
    
    # --- map historical threshold in physical space to future CDF space
    # temperature
    ta_thre <- cdfnor(
      quanor(threshold, para = c(mean(tas_his, na.rm = TRUE),sd(detrend(tas_his), na.rm = TRUE))),
      para = c(mean(tas_fut, na.rm = TRUE), sd(detrend(tas_fut), na.rm = TRUE))
    )
    # precipitation
    pr_thre <- cdfnor(
      quanor(threshold, para = c(mean(pr_his, na.rm = TRUE), sd(detrend(pr_his), na.rm = TRUE))),
      para = c(mean(pr_fut, na.rm = TRUE), sd(detrend(pr_fut), na.rm = TRUE))
    )
    
    # --- mean+var adjusted (keep dependence = future copula)
    p_pr_mv[idx] <- 1 - pr_thre  - threshold + BiCopCDF(pr_thre,  threshold, cop_fut)
    p_ta_mv[idx] <- 1 - threshold - ta_thre  + BiCopCDF(threshold, ta_thre,  cop_fut)
    
    # --- variability-only (set mean = 0 in both periods)
    ta_thre_det <- cdfnor(
      quanor(threshold, para = c(0, sd(detrend(tas_his), na.rm = TRUE))),
      para = c(0, sd(detrend(tas_fut), na.rm = TRUE))
    )
    pr_thre_det <- cdfnor(
      quanor(threshold, para = c(0, sd(detrend(pr_his), na.rm = TRUE))),
      para = c(0, sd(detrend(pr_fut), na.rm = TRUE))
    )
    
    p_pr_v[idx] <- 1 - pr_thre_det - threshold + BiCopCDF(pr_thre_det, threshold, cop_fut)
    p_ta_v[idx] <- 1 - threshold   - ta_thre_det + BiCopCDF(threshold, ta_thre_det, cop_fut)
    
    # --- dependence-only (fixed thresholds in uniform space)
    p_d2[idx] <- 1 - threshold - threshold + BiCopCDF(threshold, threshold, cop_fut)
    p_d1[idx] <- 1 - threshold - threshold + BiCopCDF(threshold, threshold, cop_his)
    
    # --- full future probability (thresholds mapped into future)
    p_fut[idx] <- 1 - pr_thre - ta_thre + BiCopCDF(pr_thre, ta_thre, cop_fut)
  }
  
  # --- LMF + decomposition (naming aligned)
  lmf      <- p_fut / p_d1
  lmf_pr_m <- p_pr_mv / p_pr_v
  lmf_ta_m <- p_ta_mv / p_ta_v
  lmf_pr_v <- p_pr_v / p_d2
  lmf_ta_v <- p_ta_v / p_d2
  lmf_dep  <- p_d2 / p_d1
  

  return(list(
    lmf      = lmf,
    lmf_pr_m = lmf_pr_m,
    lmf_ta_m = lmf_ta_m,
    lmf_pr_v = lmf_pr_v,
    lmf_ta_v = lmf_ta_v,
    lmf_dep  = lmf_dep
  ))
}

#----------A specific example
#' pr_da = readRDS("D:/17LUCC_Compound/CMIP6_data_use_for_supercomputer/01cmip6_data/1pct_2degree/pr/pr_Amon_ACCESS-ESM1-5_1pctCO2_r1i1p1f1_gn_010101-025012.nc_2degree")[,,1:1680]*86400*30
#' tas_da = readRDS("D:/17LUCC_Compound/CMIP6_data_use_for_supercomputer/01cmip6_data/1pct_2degree/tas/tas_Amon_ACCESS-ESM1-5_1pctCO2_r1i1p1f1_gn_010101-025012.nc_2degree")[,,1:1680]
#' i = 120
#' j = 50
#' threshold = 0.9
#' window_yr = 30
#' result_i = cal_lmf(tas_da, pr_da, i, j, threshold, window_yr)


#---------check function
#' ddi = readRDS("D:/17LUCC_Compound/Results/lmf_1pct_0.90/hot-wet/lmf_far_monthly_tas_pr_140yr_detrend__1pctco2_ACCESS-ESM1-5_0.9_30")
#' len = dim(ddi)[3]/12
#' 
#' p_dep_pr_ta_sd_mean = ddi[,,(len*05+1):(len*06)]#p_dep_sm_ta_sd_mean 
#' p_dep_pr_sd_mean    = ddi[,,(len*06+1):(len*07)]#p_dep_sm_sd_mean    
#' p_dep_ta_sd_mean    = ddi[,,(len*07+1):(len*08)]#p_dep_ta_sd_mean    
#' p_dep_pr_sd         = ddi[,,(len*08+1):(len*09)]#p_dep_sm_sd         
#' p_dep_ta_sd         = ddi[,,(len*09+1):(len*10)]#p_dep_ta_sd         
#' p_dep2              = ddi[,,(len*10+1):(len*11)]#p_dep2              
#' p_dep1              = ddi[,,(len*11+1):(len*12)]#p_dep1 
#' lmf_tot  = p_dep_pr_ta_sd_mean/p_dep1
#' 
#' plot(result_i$lmf, lmf_tot[i,j,])
#' abline(0,1, col='blue')