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

from plot_helpers import CO, bold_ticks


def fig6_12():
    # Integrate f(x) over [a, b] using simp-rom, and gauss

    ### pick your example
    # exp(3x)
    a = 0
    b = 1
    exact = (np.exp(3) - 1) / 3

    # tanh(100*x-30)
    # a = 0; b = 1
    # r1 = 1 + np.exp(-140); r2 = 1 + np.exp(-60)
    # exact = 0.4 + np.log(r1/r2)

    # 1/(25*x^2 + 1)
    # a = -1; b = 1
    # exact = (2*np.arctan(5))/5

    # composite simpson-romberg
    n = 2
    i = 0
    S0 = simp(a, b, n)
    fnS = [3]
    R0 = S0
    errS = [abs(R0 - exact) / exact]
    while n < 16:
        n = 2 * n
        i = i + 1
        fnS.append(n + 1)
        S = simp(a, b, n)
        R = (16 * S - S0) / 15
        errS.append(abs(R - exact) / exact)
        S0 = S
        R0 = R

    # gauss
    G = np.zeros(16)
    fnG = np.arange(1, 17)
    errG = np.zeros(16)
    for k in range(1, 17):
        G[k - 1] = gauss(a, b, k)
        errG[k - 1] = abs(exact - G[k - 1]) / exact

    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.semilogy(fnS, errS, "--*r", markersize=8, linewidth=1.6)
    ax.grid(True)
    ax.set_xlim(3, 17)
    ax.set_ylim(1e-16, 1)
    ax.set_xticks([3, 5, 7, 9, 11, 13, 15, 17, 19, 23])
    ax.set_yticks([1e-16, 1e-12, 1e-8, 1e-4, 1])
    ax.set_xlabel("Number of Points Used")
    ax.set_ylabel("Error")
    ax.semilogy(fnG, errG, "--ob", markersize=8, linewidth=1.6)

    ax.legend(
        [" Simp-Romb", " Gaussian"],
        loc="upper right",
        prop={"size": 14, "weight": "bold"},
    )

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

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


def f(x):
    y = np.exp(3 * x)
    # y = np.tanh(100*x - 30)
    # y = 1/(25*x**2 + 1)
    return y


### gaussian quadrature: n = number of points
def gauss(a, b, n):
    beta = 0.5 / np.sqrt(1 - (2 * np.arange(1, n)) ** (-2.0))
    A = np.diag(beta, 1) + np.diag(beta, -1)
    x, V = np.linalg.eigh(A)
    x = a + 0.5 * (b - a) * (1 + x)
    w = 2 * V[0, :] ** 2
    sum = 0.0
    for i in range(len(x)):
        sum = sum + w[i] * f(x[i])
    s = sum * (b - a) / 2
    return s


### 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_12()
