import numpy as np


def step(a, b, n):
    return (b - a) / n


# this is basically linspace
def interpolation_nodes(a, b, n):
    h = step(a, b, n - 1)
    return [a + i * h for i in range(n)]


def composite_trapezoidal(f, a, b, n):
    h = step(a, b, n)
    x = interpolation_nodes(a, b, n + 1)
    return h / 2 * sum([f(x[i]) + f(x[i + 1]) for i in range(n)])


def error_trapezoidal(f, a, b, n, expected):
    return expected - composite_trapezoidal(f, a, b, n)


def format_error(f, a, b, n, expected):
    return f"error for {f} from {a} to {b} with {n} subintervals is {error_trapezoidal(f, a, b, n, expected)}."


def find_subinterval_count(f, a, b, expected, error_threshold):
    n = 1
    e = lambda n: error_trapezoidal(f, a, b, n, expected)
    while np.abs(e(n)) > error_threshold:
        n += 1
    return n


def estimate_convergence_rate(f, a, b, expected):
    m = 20
    errors = [np.log(abs(error_trapezoidal(f, a, b, n, expected))) for n in range(2, m)]
    steps = [np.log(step(a, b, n)) for n in range(2, m)]
    coeffs = np.polyfit(steps, errors, 1)  # log-log regression line
    p = coeffs[0]
    C = np.exp(coeffs[1])
    return p, C


def adaptive_trapezoidal(f, a, b, tol):
    T1 = (b - a) / 2 * (f(a) + f(b))
    c = (a + b) / 2
    T2 = (b - a) / 4 * (f(a) + 2 * f(c) + f(b))

    E1 = 4 / 3 * (T2 - T1)

    if abs(E1) <= tol:
        E2 = 1 / 3 * (T2 - T1)
        return T2 + E2, 1
    else:
        left_val, left_count = adaptive_trapezoidal(f, a, c, tol * 0.7)
        right_val, right_count = adaptive_trapezoidal(f, c, b, tol * 0.7)
        return left_val + right_val, left_count + right_count


print(format_error(np.exp, 0, 1, 16, np.e - 1))
print(f"subintervals: {find_subinterval_count(np.exp, 0, 1, np.e - 1, 1e-5)}")
print(f"exp p: {estimate_convergence_rate(np.exp, 0, 1, np.e - 1)[0]:.5f}")
print(f"sqrt p: {estimate_convergence_rate(np.sqrt, 0, 1, 2/3)[0]:.5f}")
print(f"adaptive estimate exp: {adaptive_trapezoidal(np.exp, 0, 1, 1e-5)[0]:.5f}")
print(f"adaptive estimate sqrt: {adaptive_trapezoidal(np.sqrt, 0, 1, 1e-5)[0]:.5f}")
print(f"adaptive n subintervals exp: {adaptive_trapezoidal(np.exp, 0, 1, 1e-5)[1]:.0f}")
print(
    f"adaptive n subintervals sqrt: {adaptive_trapezoidal(np.sqrt, 0, 1, 1e-5)[1]:.0f}"
)