import numpy as np
import matplotlib.pyplot as plt
from cycler import cycler
from scipy.interpolate import CubicSpline

from plot_helpers import CO, bold_ticks
from splineA import splineA


def fig6_10(seed=None):
    # this requires the file splineA.m

    ##### pick your example
    ### y=exp(3x)
    a = 0
    b = 1
    exact = (np.exp(3) - 1) / 3
    ### y=sin(x)
    # a = 0; b = np.pi
    # exact = 2
    ### f=x^4;
    # a = 0; b = 10
    # exact = 10**5 / 5

    rng = np.random.default_rng(seed)

    nd = 6  # number of data points
    errT = np.zeros(10)
    errS = np.zeros(10)
    errSK = np.zeros(10)

    for ic in range(10):
        # generate data
        xd = np.sort(a + (b - a) * rng.random(nd))
        xd[0] = a
        xd[nd - 1] = b
        yd = np.zeros(nd)
        for i in range(nd):
            yd[i] = f(xd[i])

        h = np.zeros(nd - 1)
        sm = np.zeros(nd - 1)
        for i in range(nd - 1):
            h[i] = xd[i + 1] - xd[i]
            sm[i] = xd[i] + 0.5 * h[i]

        # trap
        I_T = h[0] * yd[0]
        for i in range(1, nd - 1):
            I_T = I_T + (h[i] + h[i - 1]) * yd[i]
        I_T = 0.5 * (I_T + h[nd - 2] * yd[nd - 1])

        # natural
        sdN = splineA(xd, yd, sm)
        I_S = I_T / 3
        for i in range(nd - 1):
            I_S = I_S + 2 * h[i] * sdN[i] / 3

        # not-a-knot
        sdK = CubicSpline(xd, yd)(sm)
        I_SK = I_T / 3
        for i in range(nd - 1):
            I_SK = I_SK + 2 * h[i] * sdK[i] / 3

        errT[ic] = abs((exact - I_T) / exact)
        errS[ic] = abs((exact - I_S) / exact)
        errSK[ic] = abs((exact - I_SK) / exact)

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

    ax.semilogy(np.arange(1, 11), errT, "ob", markersize=9, linewidth=2)
    # ax.semilogy(np.arange(1, 11), errS, "dk", markersize=9, linewidth=2)
    ax.semilogy(np.arange(1, 11), errSK, "+r", markersize=9, linewidth=2)
    ax.grid(True)
    ax.set_xlim(1, 10)
    ax.set_ylim(1e-4, 1)
    ax.set_yticks([1e-4, 1e-2, 1])
    ax.set_xlabel("Experiment")
    ax.set_ylabel("Error")

    ax.legend(
        [" Trap", " Spline"],
        loc="center left",
        bbox_to_anchor=(1.01, 0.5),
        prop={"size": 14, "weight": "bold"},
    )
    # ax.legend([" Trap", " Natural", " Not"], loc="center left",
    #           bbox_to_anchor=(1.01, 0.5))

    bold_ticks(ax, 14)
    for spine in ax.spines.values():
        spine.set_visible(True)
    plt.show()

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


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


if __name__ == "__main__":
    fig6_10()
