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


# Source files
library(MixStable)

# Initialize parameter history lists
L_alpha <- c(1.0, 1.1)
L_beta  <- c(0.0, 0.1)
L_delta <- c(1.0, 1.1)
L_omega <- c(0.0, 0.1)

# Ensure X1 is defined, e.g.:
X1 <- stabledist::rstable(1200, 1.2, 0.5, 1.0, 3.0, pm = 1)

# === Update α ===
alpha_new  <- tail(L_alpha, 1)
alpha_prev <- head(tail(L_alpha, 2), 1)
beta       <- tail(L_beta, 1)
delta      <- tail(L_delta, 1)
omega      <- tail(L_omega, 1)

# === Update α ===
f_alpha_new <- tryCatch(normalized_grad_alpha(alpha_new, beta, delta, omega, X1), error = function(e) NA)
f_alpha_prev <- tryCatch(normalized_grad_alpha(alpha_prev, beta, delta, omega, X1), error = function(e) NA)

if (is.finite(f_alpha_new) && is.finite(f_alpha_prev)) {
  alpha_updated <- false_position_update(
    alpha_prev, alpha_new, f_alpha_prev, f_alpha_new,
    function(a) tryCatch(normalized_grad_alpha(a, beta, delta, omega, X1), error = function(e) NA)
  )
  L_alpha <- c(L_alpha, min(max(alpha_updated, 0.2), 2.0))
} else {
  warning("Skipping alpha update due to invalid gradient.")
  L_alpha <- c(L_alpha, alpha_new)
}

# === Update β ===
beta_new <- tail(L_beta, 1)
beta_prev <- head(tail(L_beta, 2), 1)

f_beta_new <- tryCatch(normalized_objective_beta(beta_new, X1, L_alpha, L_delta, L_omega), error = function(e) NA)
f_beta_prev <- tryCatch(normalized_objective_beta(beta_prev, X1, L_alpha, L_delta, L_omega), error = function(e) NA)

if (is.finite(f_beta_new) && is.finite(f_beta_prev)) {
  beta_updated <- false_position_update(
    beta_prev, beta_new, f_beta_prev, f_beta_new,
    function(b) tryCatch(normalized_objective_beta(b, X1, L_alpha, L_delta, L_omega), error = function(e) NA)
  )
  L_beta <- c(L_beta, min(max(beta_updated, -1.0), 1.0))
} else {
  warning("Skipping beta update due to invalid gradient.")
  L_beta <- c(L_beta, beta_new)
}

# === Update δ ===
delta_new <- tail(L_delta, 1)
delta_prev <- head(tail(L_delta, 2), 1)

f_delta_new <- tryCatch(normalized_objective_delta(delta_new, X1, L_alpha, L_beta, L_omega), error = function(e) NA)
f_delta_prev <- tryCatch(normalized_objective_delta(delta_prev, X1, L_alpha, L_beta, L_omega), error = function(e) NA)

if (is.finite(f_delta_new) && is.finite(f_delta_prev)) {
  delta_updated <- false_position_update(
    delta_prev, delta_new, f_delta_prev, f_delta_new,
    function(d) tryCatch(normalized_objective_delta(d, X1, L_alpha, L_beta, L_omega), error = function(e) NA)
  )
  L_delta <- c(L_delta, max(delta_updated, 1e-3))
} else {
  warning("Skipping delta update due to invalid gradient.")
  L_delta <- c(L_delta, delta_new)
}

# === Update ω ===
omega_new <- tail(L_omega, 1)
omega_prev <- head(tail(L_omega, 2), 1)

f_omega_new <- tryCatch(normalized_objective_omega(omega_new, X1, L_alpha, L_beta, L_delta), error = function(e) NA)
f_omega_prev <- tryCatch(normalized_objective_omega(omega_prev, X1, L_alpha, L_beta, L_delta), error = function(e) NA)

if (is.finite(f_omega_new) && is.finite(f_omega_prev)) {
  omega_updated <- false_position_update(
    omega_prev, omega_new, f_omega_prev, f_omega_new,
    function(o) tryCatch(normalized_objective_omega(o, X1, L_alpha, L_beta, L_delta), error = function(e) NA)
  )
  L_omega <- c(L_omega, omega_updated)
} else {
  warning("Skipping omega update due to invalid gradient.")
  L_omega <- c(L_omega, omega_new)
}


# === Display Final Estimates ===
cat("✅ Final Parameter Estimates:\n")
cat(sprintf("alpha = %.4f\n", tail(L_alpha, 1)))
cat(sprintf("beta  = %.4f\n", tail(L_beta, 1)))
cat(sprintf("delta = %.4f\n", tail(L_delta, 1)))
cat(sprintf("omega = %.4f\n", tail(L_omega, 1)))