# Clean environment and load libraries
rm(list = ls())  # Clear environment to avoid conflicts


# Source files
library(MixStable)

# Generate synthetic data
# Generate smaller synthetic dataset for CRAN
data <- generate_synthetic_data(n = 200)  # instead of large default

# Skip plotting on CRAN
if (!identical(Sys.getenv("NOT_CRAN"), "true")) {
  cat("Skipping heavy plots on CRAN\n")
  quit("no")
}

# Estimate bandwidth (if needed)
h_n <- tryCatch({
  stats::bw.SJ(data)
}, error = function(e) {
  message("Bandwidth estimation failed: ", e$message)
  1.0
})

# Run EM methods with error handling
safe_run_method <- function(method_func, data, method_name) {
  cat("Running", method_name, "...\n")
  tryCatch({
    result <- method_func(data)
    cat("✓", method_name, "completed successfully\n")
    return(result)
  }, error = function(e) {
    cat("✗", method_name, "failed:", e$message, "\n")
    return(NULL)
  })
}

methods <- list()
methods[["EM_Estimated"]] <- safe_run_method(em_fit_alpha_stable_mixture, data, "EM_Estimated")
methods[["EM_ECF_recursive_ecf"]] <- safe_run_method(em_estimate_stable_recursive_ecf, data, "EM_ECF_recursive_ecf")
methods[["EM_ECF_Kernel"]] <- safe_run_method(em_estimate_stable_kernel_ecf, data, "EM_ECF_Kernel")
methods[["EM_ECF_weighted_ols"]] <- safe_run_method(em_estimate_stable_weighted_ols, data, "EM_ECF_weighted_ols")
methods[["EM_ECF_from_cdf"]] <- safe_run_method(em_estimate_stable_from_cdf, data, "EM_ECF_from_cdf")

# Remove failed methods
methods <- methods[!sapply(methods, is.null)]

if (length(methods) == 0) {
  stop("All methods failed. Check your data and functions.")
}

to_list <- function(p) {
  if (is.list(p)) return(p)
  names(p) <- c("alpha", "beta", "gamma", "delta")
  as.list(p)
}

# Plotting
x_vals <- seq(min(data), max(data), length.out = 1000)

for (name in names(methods)) {
  result <- methods[[name]]

  # Skip if result is NULL or doesn't have required structure
  if (is.null(result) || is.null(result$params1) || is.null(result$params2)) {
    cat("Skipping", name, "- invalid result structure\n")
    next
  }

  tryCatch({
    param1 <- to_list(result$params1)
    param2 <- to_list(result$params2)
    weight <- result$weight

    # Validate weight
    if (is.null(weight) || !is.numeric(weight) || weight < 0 || weight > 1) {
      cat("Skipping", name, "- invalid weight\n")
      next
    }

    p1_vals <- unpack_params(param1)
    p2_vals <- unpack_params(param2)

    # Use explicit parameter names to avoid confusion
    pdf_comp1 <- r_stable_pdf(x_vals,
                              alpha = param1$alpha,
                              beta = param1$beta,
                              scale = param1$gamma,
                              location = param1$delta)

    pdf_comp2 <- r_stable_pdf(x_vals,
                              alpha = param2$alpha,
                              beta = param2$beta,
                              scale = param2$gamma,
                              location = param2$delta)

    mixture_pdf <- weight * pdf_comp1 + (1 - weight) * pdf_comp2

    cat(sprintf("\n🔍 Plotting %s\n", name))
    cat(sprintf("Component 1: α=%.3f, β=%.3f, γ=%.3f, δ=%.3f\n",
                param1$alpha, param1$beta, param1$gamma, param1$delta))
    cat(sprintf("Component 2: α=%.3f, β=%.3f, γ=%.3f, δ=%.3f\n",
                param2$alpha, param2$beta, param2$gamma, param2$delta))
    cat(sprintf("Weight: %.3f\n", weight))

    plot_title <- sprintf("%s vs Data", gsub("_", " ", name))
    png_filename <- sprintf("mixture_alpha_stable_%s.png", tolower(name))

    png(png_filename, width = 900, height = 500)
    hist(data, breaks = 40, freq = FALSE, col = rgb(0.7, 0.7, 0.7, 0.4),
         main = plot_title, xlab = "x", ylab = "Density")
    lines(x_vals, mixture_pdf, col = "blue", lwd = 2)
    legend("topright",
           legend = c("Data", "Estimated Mixture"),
           col = c("gray", "blue"),
           lty = c(NA, 1),
           pch = c(15, NA),
           bty = "n")
    dev.off()

    cat("✓ Plot saved as", png_filename, "\n")

  }, error = function(e) {
    cat("Error plotting", name, ":", e$message, "\n")
  })
}