Using custom cost functions

calculate_penalized_aic <- function() {
    
    # --- 1. Setup and Helper Functions -----------------------------------------
    
    # Helper for piecewise-linear interpolation
    # It calculates a y-value for a given x based on a set of (xp, yp) points.
    piecewise_linear_interp <- function(x, xp, yp) {
        # Ensure inputs are valid
        stopifnot(length(xp) == length(yp), length(xp) >= 2)
        
        # Return 0 for NA inputs and clamp values outside the interpolation range
        if (is.na(x)) return(0)
        if (x <= xp[1]) return(yp[1])
        if (x >= tail(xp, 1)) return(tail(yp, 1))
        
        # Perform linear interpolation
        i <- findInterval(x, xp)
        x0 <- xp[i]; x1 <- xp[i + 1]
        y0 <- yp[i]; y1 <- yp[i + 1]
        
        y0 + (y1 - y0) * (x - x0) / (x1 - x0)
    }
    
    # --- 2. Penalty Calculation Definitions ------------------------------------
    # These functions define the penalty curves for different metrics.
    
    # Penalty for THETA RSE (%): 10%->0, 30%->20, 75%->100, 1000%->100000
    pen_theta_rse <- function(rse_pct) {
        piecewise_linear_interp(rse_pct,
                                xp = c(10, 30, 75, 1000),
                                yp = c(0, 20, 100, 100000))
    }
    
    # Penalty for OMEGA RSE (%): 20%->0, 40%->30, 65%->100, 1000%->100000
    pen_omega_rse <- function(rse_pct) {
        piecewise_linear_interp(rse_pct,
                                xp = c(20, 40, 65, 1000),
                                yp = c(0, 30, 100, 100000))
    }
    
    # Penalty for OMEGA value (on the standard deviation scale)
    pen_omega_value <- function(omega_sd_val) {
        # Convert SD to variance for interpolation, as thresholds are often
        # considered on the variance scale.
        omega_var_val <- omega_sd_val^2
        xp <- c(0.0000, 0.0010, 0.0100, 1.0000, 1.6500, 2.3000, 3.0000, 3.5000)
        yp <- c(1e6,    1e3,    0,      0,      100,    500,    1e4,    1e6)
        piecewise_linear_interp(omega_var_val, xp, yp)
    }
    
    # Penalty for OMEGA shrinkage (%): 40%->0, 100%->200
    pen_omega_shrink <- function(shr_pct) {
        piecewise_linear_interp(shr_pct, xp = c(40, 100), yp = c(0, 200))
    }
    
    # Penalty for Condition Number: 500->0, 1000->100, 10000->1000, 1e6->1e7
    pen_condition_number <- function(cond_num) {
        piecewise_linear_interp(cond_num,
                                xp = c(500, 1000, 10000, 1e6),
                                yp = c(0, 100, 1000, 1e7))
    }
    
    # --- 3. Extract Metrics from Monolix Project -------------------------------
    
    # Project directory for result files
    proj_dir <- getProjectSettings()$directory
    
    # AIC (from importance sampling)
    AIC <- getEstimatedLogLikelihood()$importanceSampling[["AIC"]]
    if (is.null(AIC) || is.na(AIC)) {
        warning("AIC from Importance Sampling not found. Using AIC from Linearization.")
        AIC <- getEstimatedLogLikelihood()$linearization[["AIC"]]
        if (is.null(AIC)) stop("Could not retrieve AIC value.")
    }
    
    # Population parameter estimates and standard errors
    pop_params <- getEstimatedPopulationParameters()
    std_errs <- getEstimatedStandardErrors()$stochasticApproximation
    
    # Filter for THETA (fixed effects) and OMEGA (random effect variances)
    is_theta <- grepl("(_pop|beta_)", names(pop_params)) & !grepl("omega", names(pop_params))
    is_omega <- grepl("^omega_", names(pop_params)) & !grepl("corr|cov", names(pop_params))
    
    theta_rse <- std_errs[is_theta, "rse"]
    omega_rse <- std_errs[is_omega, "rse"]
    omega_values <- pop_params[is_omega]
    
    # Shrinkage
    shrinkage_file <- file.path(proj_dir, "individualParameters", "shrinkage.txt")
    shrinkage <- if (file.exists(shrinkage_file)) {
        read.csv(shrinkage_file)$shrinkage_mode
    } else {
        warning("Shrinkage file not found. Skipping shrinkage penalty.")
        numeric(0)
    }
    
    # Condition number of the Fisher Information Matrix (FIM)
    fim_file <- file.path(proj_dir, "FisherInformation", "covarianceEstimatesSA.txt")
    condition_number <- if (file.exists(fim_file)) {
        fim <- as.matrix(read.csv(fim_file, header = FALSE, row.names = 1))
        eigen_values <- eigen(fim, only.values = TRUE)$values
        max(eigen_values) / min(eigen_values)
    } else {
        warning("FIM file not found. Skipping condition number penalty.")
        NA
    }
    
    # --- 4. Compute Penalties -------------------------------------------------
    
    # `sum(sapply(x, fun))` gracefully handles empty vectors (sum = 0)
    p_theta_rse <- sum(sapply(theta_rse, pen_theta_rse))
    p_omega_rse <- sum(sapply(omega_rse, pen_omega_rse))
    p_omega_val <- sum(sapply(omega_values, pen_omega_value))
    p_shrinkage <- sum(sapply(shrinkage, pen_omega_shrink))
    p_cond_num  <- pen_condition_number(condition_number)
    
    # --- 5. Aggregate and Return Results --------------------------------------
    
    penalty_breakdown <- c(
        THETA_RSE        = p_theta_rse,
        OMEGA_RSE        = p_omega_rse,
        OMEGA_value      = p_omega_val,
        OMEGA_shrinkage  = p_shrinkage,
        Condition_number = p_cond_num
    )
    
    total_penalty <- sum(penalty_breakdown, na.rm = TRUE)
    penalized_AIC <- AIC + total_penalty
    
    penalized_AIC
}