import numpy as np
import matplotlib.pyplot as plt
from cycler import cycler

from plot_helpers import CO, bold_ticks


def fig6_9():
    # Integrate f(x) over [a, b] using composite simpson and simpson-romberg

    a = 0
    b = 1
    exact = (np.exp(3) - 1) / 3

    # interv = number of subinternals
    interv = np.array([4, 8, 16, 32, 64, 128])
    nk = len(interv)
    tol = 1e-10
    n = 2
    intervals = [n]
    S = [simp(a, b, 2)]
    errS = [abs(S[0] - exact) / exact]
    err = 10 * tol
    R = [1.0]
    errR = [-1.0]  # this is a fake value
    k = 0

    while err > tol:
        k = k + 1
        n = 2 * n
        intervals.append(n)
        S.append(simp(a, b, n))
        R.append((16 * S[k] - S[k - 1]) / 15)
        err = abs(1 - R[k] / R[k - 1])
        errR.append(abs(R[k] - exact) / exact)
        errS.append(abs(S[k] - exact) / exact)

    intervals = np.asarray(intervals)
    errS = np.asarray(errS)
    errR = np.asarray(errR)

    plt.close("all")
    # get(gcf)
    fig, ax = plt.subplots(figsize=(8.8, 3.3))
    fig.subplots_adjust(left=0.10, right=0.98, bottom=0.22, top=0.98)
    ax.set_prop_cycle(cycler(color=CO))

    ax.loglog(intervals, errS, "--ob", markersize=8, linewidth=1.6)
    ax.loglog(intervals[1 : k + 1], errR[1 : k + 1], "--*r", markersize=8, linewidth=1.6)
    ax.grid(True)
    ax.set_xlabel("n-axis")
    ax.set_ylabel("Error")
    ax.legend(
        [" Simp", " Simp-Romb"],
        loc="lower left",
        prop={"size": 14, "weight": "bold"},
    )
    ax.set_xlim(1, 2e2)
    ax.set_ylim(1e-14, 1e-1)
    ax.set_yticks([1e-14, 1e-10, 1e-6, 1e-2])
    ax.grid(False, which="minor", axis="y")
    bold_ticks(ax, 14)
    plt.show()

    # plt.savefig('/Users/mark/Desktop/rom2.eps')


def f(x):
    y = np.exp(3 * x)
    return y


### composite simpson: n = number of subintervals
def simp(a, b, n):
    xd = np.linspace(a, b, n + 1)
    h = xd[1] - xd[0]
    s = f(xd[0])
    ff2 = 0.0
    for j in range(1, n, 2):
        ff2 = f(xd[j + 1])
        s = s + 4 * f(xd[j]) + 2 * ff2
    s = h * (s - ff2) / 3
    return s


if __name__ == "__main__":
    fig6_9()
