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

from plot_helpers import CO, bold_ticks


def fig6_17(nk_midpoint=20, nk_simpson=14):
    # Integrate f(x) over [a, b] using midpoint and simpson's rule
    # with compensated summation
    #

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

    # a = 0; b = 3*np.pi
    # exact = 3*np.pi

    # midpoint
    print("\n Error Using Midpint Rule")
    print("\n n  \t I_M \t\t I_M + CS")
    nk = nk_midpoint
    n = 120
    intervM = np.zeros(nk, dtype=int)
    MM = np.zeros(nk)
    M = np.zeros(nk)
    errM = np.zeros(nk)
    errMM = np.zeros(nk)

    for k in range(nk):
        intervM[k] = 2 * n
        n = int(intervM[k])
        xd = np.linspace(a, b, n + 1)
        h = xd[1] - xd[0]
        sumMM = f(xd[0] + 0.5 * h)
        err = 0.0
        sumM = f(xd[0] + 0.5 * h)
        for j in range(1, n):
            z = f(xd[j] + 0.5 * h) + err
            q = sumMM
            sumMM = q + z
            err = q - sumMM + z
            sumM = sumM + f(xd[j] + 0.5 * h)
        MM[k] = h * sumMM
        M[k] = h * sumM
        errM[k] = abs(exact - M[k])
        errMM[k] = abs(exact - MM[k])
        print(f"{intervM[k]:d} \t {errM[k]:7.1e} \t {errMM[k]:7.1e}")
    print()

    # simpson
    print("\n Error Using Simpsons Rule")
    print("\n n  \t I_S \t\t I_S + CS")
    nk = nk_simpson
    n = 20
    # nk = 10
    # n = 10
    intervS = np.zeros(nk, dtype=int)
    S = np.zeros(nk)
    SS = np.zeros(nk)
    errS = np.zeros(nk)
    errSS = np.zeros(nk)

    for k in range(nk):
        intervS[k] = 2 * n
        n = int(intervS[k])
        xd = np.linspace(a, b, n + 1)
        h = xd[1] - xd[0]
        sumS = f(xd[0])
        err = 0.0
        sumSS = f(xd[0])
        ff2 = 0.0
        for j in range(1, n, 2):
            ff2 = f(xd[j + 1])
            z = 4 * f(xd[j]) + 2 * ff2 + err
            q = sumSS
            sumSS = q + z
            err = q - sumSS + z
            sumS = sumS + 4 * f(xd[j]) + 2 * ff2
        S[k] = h * (sumS - ff2) / 3
        SS[k] = h * (sumSS - ff2) / 3
        errS[k] = abs(exact - S[k])
        errSS[k] = abs(exact - SS[k])
        print(f"{intervS[k]:d} \t {errS[k]:7.1e} \t {errSS[k]:7.1e}")
    print()

    # plot error using loglog plot
    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(intervM, errM, "-or", markersize=8, linewidth=1.6)
    ax.loglog(intervM, errMM, "-*b", markersize=8, linewidth=1.4)
    ax.loglog(intervS, errS, "--or", markersize=8, linewidth=1.6)
    ax.loglog(intervS, errSS, "--*b", markersize=8, linewidth=1.4)
    ax.grid(True)
    ax.set_xlabel("n-axis")
    ax.set_ylabel("Error")
    ax.set_xlim(1, 1e8)
    ax.set_ylim(1e-16, 1e-4)
    ax.set_yticks([1e-16, 1e-12, 1e-8, 1e-4])
    ax.legend(
        [" Original", " Comp Sum"],
        loc="upper right",
        prop={"size": 12, "weight": "bold"},
    )

    bold_ticks(ax, 14)
    plt.show()

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


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


# def df(x):
#     y = np.exp(x)
#     return y


if __name__ == "__main__":
    fig6_17()
