import numpy as np

# Checking current approximation
# epsilon = H^(-1)(1-R)
def inv_approx(R):
    return 0.5*(1-np.sqrt(1-(np.ones(len(R))-R)**(4/3)))

# # H(epsilon) = 1-R

def H(e):
    return -e*np.log2(e) - (1-e)*np.log2(1-e)

# epsilon = inv_approx([0.9])
# H = H(epsilon)

# print(1-H) # should be R



from scipy.optimize import minimize_scalar

def H(p):
    return -p * np.log2(p) - (1 - p) * np.log2(1 - p)

def inverse_H(arr):
    vals = []
    for val in arr:
        # Define a function that returns the absolute difference between H(p) and the desired value
        func = lambda p: abs(H(p) - val)
        # Use minimize_scalar to find the value of p that minimizes the absolute difference
        result = minimize_scalar(func, bounds=(1e-15, 1-1e-15), method='bounded')
        vals.append(result.x)
    return np.array(vals)

# Example usage:
R = np.arange(0, 1, 1/50)
H_inv = inverse_H(1-R)
H_inv[H_inv > 0.5] = 1 - H_inv[H_inv > 0.5]
H = H(H_inv)

print(H_inv[10])
print(R[10])



# #Using inverse function technique
# def binary_ent(p):
#     if p == 0 or p == 1:
#         return 0
#     return -p * np.log2(p) - (1 - p) * np.log2(1 - p)

# def inverse_binary_entropy(h):
#     from scipy.optimize import newton

#     if h == 0:
#         return 0
#     elif h == 1:
#         return 1
    
#     # Define the equation for finding the root
#     def equation(p):
#         return binary_ent(p) - h

#     # Provide a better initial guess
#     if h < 0.5:
#         initial_guess = h * 0.9
#         while equation(initial_guess) > 0:
#             initial_guess *= 0.9
#     else:
#         initial_guess = 1 - (1 - h) * 0.9
#         while equation(initial_guess) < 0:
#             initial_guess = 1 - (1 - initial_guess) * 0.9

#     # Use Newton's method with a more robust stopping criterion
#     tolerance = 1e-12
#     max_iter = 10000
#     root = newton(equation, initial_guess, tol=tolerance, maxiter=max_iter)

   
#     # Use newton's method to find the root
#     return root