import time

import numpy as np


ifun = 0


def simpA():
    # Integrate f(x) over [a, b] using adaptive Simpson
    # Include cpu time and number of function evaluations

    global ifun

    a = 0
    b = 1

    # n = initial number of subintervals
    n = 4

    # error tolerance
    tol = 1e-8

    # L = max number of allowed levels
    L = 25

    ### ifun = number of function evaluations
    ifun = 0

    tic = time.perf_counter()
    nf = 1
    m = 0.5 * (a + b)
    fm = f(m)
    fa = f(a)
    fb = f(b)
    s = simp(a, m, b, fa, fm, fb)
    F = np.array([[a, m, b, fa, fm, fb, s]], dtype=float)
    A = np.empty((0, 3), dtype=float)
    Lmax = L

    # this is the loop for the levels
    for il in range(2, L + 1):
        FF = np.empty((0, 7), dtype=float)

        ### use simpson on each remaining subinterval
        for is_ in range(nf):
            mL = 0.5 * (F[is_, 0] + F[is_, 1])
            fmL = f(mL)
            mR = 0.5 * (F[is_, 1] + F[is_, 2])
            fmR = f(mR)
            s2L = simp(F[is_, 0], mL, F[is_, 1], F[is_, 3], fmL, F[is_, 4])
            s2R = simp(F[is_, 1], mR, F[is_, 2], F[is_, 4], fmR, F[is_, 5])
            s2 = s2L + s2R
            Err = abs(s2 - F[is_, 6]) / 15
            tolA = tol * (F[is_, 2] - F[is_, 0]) / (b - a)

            ### determine if the accuracy is attained
            if Err < tolA:
                A = np.vstack(
                    (A, [F[is_, 0], F[is_, 2], (16 * s2 - F[is_, 6]) / 15])
                )
            else:
                FF = np.vstack(
                    (
                        FF,
                        [
                            F[is_, 0],
                            mL,
                            F[is_, 1],
                            F[is_, 3],
                            fmL,
                            F[is_, 4],
                            s2L,
                        ],
                    )
                )
                FF = np.vstack(
                    (
                        FF,
                        [
                            F[is_, 1],
                            mR,
                            F[is_, 2],
                            F[is_, 4],
                            fmR,
                            F[is_, 5],
                            s2R,
                        ],
                    )
                )

        nn = FF.shape

        ### if no remaining subintervals, then leave the for-loop
        nf = nn[0]
        if nf == 0:
            Lmax = il
            break
        F = FF

    ### if max number of levels exceeded, just use what was calcuated
    if nf > 0:
        for is_ in range(nf):
            A = np.vstack((A, [FF[is_, 0], FF[is_, 2], FF[is_, 6]]))

    ### add up values from each subinterval
    ss = A.shape
    sum = 0.0
    for is_ in range(ss[0]):
        sum = sum + A[is_, 2]

    print(f"Elapsed time: {time.perf_counter() - tic:.6f} seconds")
    print()
    print(f" {sum:12.4e}  &  {ifun:d}")
    print()

    # number of levels and subintevals used
    levels_used = Lmax
    subintervals = ss[0]
    print(f"levels_used = {levels_used}")
    print(f"subintervals = {subintervals}")

    return sum, ifun, levels_used, subintervals


def simp(a, m, b, fa, fm, fb):
    g = (b - a) * (fa + 4 * fm + fb) / 6
    return g


def f(x):
    global ifun
    a = np.sqrt(15)
    lambda_ = np.pi
    y = (
        (2 * a * np.sinh(x) - 2 * x)
        / (1 + a ** 2 + x ** 2 + 2 * a * np.cosh(x))
        - 1
    ) * np.sin(lambda_ * x)
    ifun = ifun + 1
    return y


if __name__ == "__main__":
    simpA()
