Untitled

using IterTools
using Distributions
using SpecialFunctions
using QuadGK

# تعریف مقادیر ثابت
const q = 2
const Y = 500
const D0 = 50
const Z0 = 0.25
const Z1 = 0.25
const Ci = 0.5
const C0_RSS = 20
const Cq = 4.22
const Cr = 2.11
const p0 = 0.95

# تعریف مقادیر متغیر برای i های مختلف
const p = [0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 0.99]
const D1 = [2734.77, 3113.924, 3544.304, 4000, 4531.646, 5113.924, 5620.253]
const D2 = [1077.7, 1051.086, 1025.39, 1000, 975.222, 951.055, 932.395]
const Z2 = [2.155, 2.102, 2.051, 2.0, 1.95, 1.902, 1.865]
const lambda_i = [0.3523 / 1000, 0.3436 / 1000, 0.3352 / 1000, 0.3269 / 1000, 0.3188 / 1000, 0.3109 / 1000, 0.3048 / 1000]

# تابع betainc برای محاسبه تابع بتای ناقص
function betainc(a, b, x)
    f(t) = t^(a - 1) * (1 - t)^(b - 1)
    return quadgk(f, 0, x)[1] / beta(a, b)
end

# تابع برای محاسبه آلفا (α)
function calculate_alpha(r, k, L)
    I = collect(product([0:r for _ in 1:k]...))  # تولید تمام ترکیبات ممکن
    R = [sum(r) for r in I]  # محاسبه مجموع هر ترکیب
    ucl = Int(floor((k * r / 2 + (L * sqrt(k * r * delta0(k) / 4)) + k * r) / 2))

    if ucl > k * r
        return 100000  # مقدار بازگشتی در صورت عدم رعایت شرط
    end

    alpha = 0.0
    for c in 0:ucl
        r_vec = I[R.==c]  # انتخاب ترکیباتی که مجموع آن‌ها برابر c است
        for i in 1:length(r_vec)
            prod_term = 1.0
            for j in 1:k
                Ji = r_vec[i][j]  # مقدار j_i برای ترکیب فعلی
                binom_val = binomial(r, Ji)  # ضریب دوجمله‌ای
                p_val = betainc(j + 1, k - j, 0.5)  # محاسبه I_{1/2}(j + 1, k - j + 1)
                prod_term *= binom_val * (p_val^(r - Ji)) * ((1 - p_val)^Ji)
            end
            alpha += prod_term
        end
    end
    return 1 - alpha
end

# تابع برای محاسبه بتا (β_i)
function calculate_beta_i(r, k, L, i)
    I = collect(product([0:r for _ in 1:k]...))  # تولید تمام ترکیبات ممکن
    R = [sum(r) for r in I]  # محاسبه مجموع هر ترکیب
    ucl = Int(floor((k * r / 2 + (L * sqrt(k * r * delta0(k) / 4)) + k * r) / 2))

    if ucl > k * r
        return 100000  # مقدار بازگشتی در صورت عدم رعایت شرط
    end

    beta_i = 0.0
    for c in 0:ucl
        r_vec = I[R.==c]  # انتخاب ترکیباتی که مجموع آن‌ها برابر c است
        for j in 1:length(r_vec)
            prod_term = 1.0
            for l in 1:k
                Ji = r_vec[j][l]  # مقدار j_i برای ترکیب فعلی
                binom_val = binomial(r, Ji)  # ضریب دوجمله‌ای
                p_val = betainc(l + 1, k - l, p[i])  # محاسبه I_p(l + 1, k - l + 1)
                prod_term *= binom_val * (p_val^(r - Ji)) * ((1 - p_val)^Ji)
            end
            beta_i += prod_term
        end
    end
    return beta_i
end

# تابع برای محاسبه P(I)
function calculate_PI(h, r, q, k, L)
    lambda = sum(lambda_i)
    beta = [calculate_beta_i(r, k, L, i) for i in 1:7]  # محاسبه beta_i برای تمام i‌ها
    PI = sum((1 - beta[i]) * (1 - exp(-lambda_i[i] * h^q)) * exp(-(lambda - lambda_i[i]) * h^q) /
             ((1 - exp(-lambda * h^q)) * (1 - beta[i] * exp(-(lambda - lambda_i[i]) * h^q))) for i in 1:7)
    return PI
end

function calculate_P_I_i(h, q, r, k, L, lambda, i)
    lambda_i_val = lambda_i[i]  # مقدار lambda_i برای اندیس i
    beta_i = calculate_beta_i(r, k, L, i)  # محاسبه beta_i برای اندیس i

    # صورت کسر
    numerator = (1 - beta_i) * (1 - exp(-lambda_i_val * h^q)) * exp(-(lambda - lambda_i_val) * h^q)

    # مخرج کسر
    denominator = (1 - exp(-lambda * h^q)) * (1 - beta_i * exp(-(lambda - lambda_i_val) * h^q))

    # جلوگیری از تقسیم بر صفر
    if denominator == 0.0
        denominator = 1e-10
    end

    # محاسبه P_I_i
    P_I_i = numerator / denominator
    return P_I_i
end


# تابع Fi
function Fi(m, i)
    return m * binomial(m - 1, i - 1) * beta(i, m - i + 1) * betainc(i, m - i + 1, 0.5)
end

# تابع برای محاسبه delta0
function delta0(m)
    return 1 - (4 / m) * sum((Fi(m, i) - 0.5)^2 for i in 1:m)
end


# تابع برای محاسبه E(C)
function calculate_EC(r, k, L, h, alpha)
    lambda = sum(lambda_i)
    PI = calculate_PI(h, r, q, k, L)
    beta = [calculate_beta_i(r, k, L, i) for i in 1:7]  # محاسبه beta_i برای تمام i‌ها

    term1 = D0 * (1 / lambda)^(1 / q) * gamma(1 / q + 1)
    term2 = alpha * Y * exp(-lambda * h^q) / (1 - exp(-lambda * h^q))

    term3 = 0.0
    for i in 1:7
        sum_j = 0.0
        for j in 1:100  # استفاده از 100 جمله برای تقریب سری
            numerator = (exp(-lambda * h^q))^j - (beta[i] * exp(-(lambda - lambda_i[i]) * h^q))^j
            denominator = exp(-lambda_i[i] * h^q) - beta[i]
            sum_j += j^(1 / q) * h * (1 - beta[i]) * (1 - exp(-lambda_i[i] * h^q)) * (numerator / denominator)

        end
        term3 += D1[i] / PI * sum_j
    end

    term4 = 0.0
    for i in 1:7  # جمع روی i از 1 تا 7
        inner_sum = 0.0
        for j in 1:100  # تقریب بینهایت با 100 جمله
            # جزء نمایی
            exp_term = exp(-(lambda - lambda_i[i]) * h^q)
            exp_term_j = exp_term^j

            # جزء گاما
            gamma_term_prev = gamma((q + 1) / q, lambda_i[i] * (j - 1) * h^q)
            gamma_term_curr = gamma((q + 1) / q, lambda_i[i] * j * h^q)
            gamma_diff = gamma_term_prev - gamma_term_curr

            # جزء تقسیم
            denominator = lambda_i[i]^(1 / k)

            # محاسبه جمله‌ی j
            numerator = (-beta[i]) * exp_term_j / (1 - (beta[i] * exp_term))
            term_j = numerator * gamma_diff / denominator

            inner_sum += term_j
        end

        # اضافه کردن به ترم ۴
        term4 += D1[i] / PI * inner_sum
    end


    term5 = C0_RSS  # مقدار اولیه C0_RSS
    rk_term = r * k * (k * Ci + (k - 1) * Cr + Cq)
    sum_term = 0.0
    for i in 1:7
        beta_i = calculate_beta_i(r, k, L, i)  # محاسبه beta_i برای هر i
        lambda_i_val = lambda_i[i]  # مقدار lambda_i برای هر i

        exp_term_lambda = exp(-lambda * h^q)
        term1 = exp_term_lambda / (1 - exp_term_lambda)

        exp_term_lambda_diff = exp(-(lambda - lambda_i_val) * h^q)
        term2 = 1 / (1 - beta_i * exp_term_lambda_diff)

        P_I_i = calculate_P_I_i(h, q, r, k, L, lambda, i)
        sum_term += (P_I_i / PI) * (term1 + term2)
    end

    # اضافه کردن به term5
    term5 += rk_term * sum_term


    term6 = sum(calculate_P_I_i(h, q, r, k, L, lambda, i) / PI * D2[i] for i in 1:7)

    EC = term1 + term2 + term3 - term4 + term5 + term6
    return EC
end

# تابع برای محاسبه E(T)
function calculate_ET(r, k, L, h, alpha)
    lambda = sum(lambda_i)
    PI = calculate_PI(h, r, q, k, L)
    beta = [calculate_beta_i(r, k, L, i) for i in 1:7]  # محاسبه beta_i برای تمام i‌ها

    term1 = (1 / lambda)^(1 / q) * gamma(1 / q + 1) + alpha * Z0 * (exp(-lambda * h^q) / (1 - exp(-lambda * h^q))) + Z1

    term2 = 0.0
    for j in 1:100  # استفاده از 100 جمله برای تقریب سری
        for i in 1:7
            numerator = (1 - beta[i]) * (1 - exp(-lambda_i[i] * h^q)) * ((exp(-lambda * h^q))^j - (beta[i] * exp(-(lambda - lambda_i[i]) * h^q))^j)
            denominator = PI * (exp(-lambda_i[i] * h^q) - beta[i])
            term2 += (j^(1 / q) * h) * numerator / denominator
        end
    end

    term3 = 0.0
    for i in 1:7
        sum_j = 0.0
        for j in 1:100  # استفاده از 100 جمله برای تقریب سری
            numerator = (1 - beta[i]) * (exp(-(lambda - lambda_i[i]) * h^q))^j
            denominator = 1 - beta[i] * exp(-(lambda - lambda_i[i]) * h^q)
            gamma1 = gamma((q + 1) / q, lambda_i[i] * (j - 1) * h^q)
            gamma2 = gamma((q + 1) / q, lambda_i[i] * j * h^q)
            sum_j += numerator / denominator * (gamma1 - gamma2) / (lambda_i[i]^(1 / k))
        end
        term3 += sum_j / PI
    end

    term4 = sum(Z2[i] * calculate_P_I_i(h, q, r, k, L, lambda, i) / PI for i in 1:7)

    ET = term1 + term2 - term3 + term4
    return ET
end

# تابع برای محاسبه E(A)
function calculate_EA(r, k, L, h)
    alpha = calculate_alpha(r, k, L)
    EC = calculate_EC(r, k, L, h, alpha)
    ET = calculate_ET(r, k, L, h, alpha)
    EA = EC / ET
    return EA
end

# تابع برای محاسبه E(A) با مقادیر دستی
function calculate_EA_manual(h, L, r, k)
    if calculate_PI(h, r, q, k, L) >= p0
        EA = calculate_EA(r, k, L, h)
        alpha = calculate_alpha(r, k, L)
        beta_values = [calculate_beta_i(r, k, L, i) for i in 1:7]  # تغییر نام به beta_values
        return EA, alpha, beta_values
    else
        return NaN, NaN, NaN
    end
end

# مثال استفاده از تابع calculate_EA_manual
h_manual = 0.1
L_manual = 1
r_manual = 2
k_manual = 2

#EA, alpha, beta_values = calculate_EA_manual(h_manual, L_manual, r_manual, k_manual)
#println("E(A) با مقادیر دستی: ", EA)
#println("alpha با مقادیر دستی: ", alpha)
#println("beta_i با مقادیر دستی: ", beta_values)

# تابع برای پیدا کردن مقادیر بهینه و چاپ آن‌ها
function find_optimal_EA()
    best_EA = Inf  # مقدار اولیه برای مینیمم (بی‌نهایت)
    optimal_params = nothing  # برای ذخیره مقادیر بهینه
    optimal_alpha = nothing  # ذخیره مقدار بهینه آلفا
    optimal_beta = nothing   # ذخیره مقادیر بهینه بتا

    # محدوده مقادیر برای h, L, r, k
    h_values = 0.1:0.1:2.0  # مقادیر h
    L_values = 1.0:2.0      # مقادیر L
    r_values = 2:5          # مقادیر r
    k_values = 1:6          # مقادیر k

    for h in h_values
        for L in L_values
            for r in r_values
                for k in k_values
                    # محاسبه E(A) با مقادیر فعلی
                    EA, alpha, beta_values = calculate_EA_manual(h, L, r, k)

                    # بررسی اگر مقدار جدید کمتر از بهترین مقدار باشد
                    if !isnan(EA) && EA < best_EA
                        best_EA = EA
                        optimal_params = (h=h, L=L, r=r, k=k)  # ذخیره مقادیر بهینه
                        optimal_alpha = alpha                  # ذخیره آلفای بهینه
                        optimal_beta = beta_values             # ذخیره مقادیر بتای بهینه
                    end
                end
            end
        end
    end

    # چاپ مقادیر بهینه و مینیمم E(A)
    println("Minimum E(A): ", best_EA)
    println("Optimal Parameters: h=", optimal_params.h, ", L=", optimal_params.L,
        ", r=", optimal_params.r, ", k=", optimal_params.k)
    println("Optimal Alpha (α): ", optimal_alpha)
    println("Optimal Beta (β): ", optimal_beta)
end

# فراخوانی تابع برای پیدا کردن مقادیر بهینه
find_optimal_EA()





function debug_calculate_EA_to_file(h_values, L_values, r_values, k_values, output_file)
    # باز کردن فایل برای نوشتن خروجی
    open(output_file, "w") do io
        println(io, "Starting debug for E(A)...")

        for h in h_values
            for L in L_values
                for r in r_values
                    for k in k_values
                        # محاسبه E(A) با مقادیر فعلی
                        EA, alpha, beta_values = calculate_EA_manual(h, L, r, k)

                        # نوشتن اطلاعات در فایل
                        println(io, "\nChecking Parameters: h=", h, ", L=", L, ", r=", r, ", k=", k)
                        println(io, "Alpha: ", alpha)
                        println(io, "Beta Values: ", beta_values)

                        # بررسی E(A)، E(C) و E(T)
                        EC = calculate_EC(r, k, L, h, alpha)
                        ET = calculate_ET(r, k, L, h, alpha)
                        println(io, "E(C): ", EC, ", E(T): ", ET, ", E(A): ", EA)

                        # بررسی مقادیر نامعتبر
                        if isnan(EA) || isnan(EC) || isnan(ET)
                            println(io, "Warning: NaN detected in calculations.")
                        elseif EA < 0 || EC < 0 || ET < 0
                            println(io, "Warning: Negative values detected!")
                            println(io, "Details: h=", h, ", L=", L, ", r=", r, ", k=", k)
                        end
                    end
                end
            end
        end

        println(io, "Debugging complete!")
    end
end

# تعریف محدوده مقادیر برای h, L, r, k
#h_values = 0.1:0.1:2.0  # مقادیر h
#L_values = 1.0:2.0         # مقادیر L
#r_values = 2:5          # مقادیر r
#k_values = 1:6          # مقادیر k

# نام فایل خروجی
#output_file = "debug_output.txt"

# فراخوانی تابع دیباگ و ذخیره در فایل
#debug_calculate_EA_to_file(h_values, L_values, r_values, k_values, output_file)

#println("Debugging complete. Output saved to file: ", output_file)