"""Literal Python conversion of fig8_10.m."""

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import minimize


xd = np.array([])
yd = np.array([])


def f(a, b, x, y):
    d = a + b * x - y
    return d


def ED(x):
    global xd, yd
    n = len(xd)
    s = 0
    for i in range(n):
        s = s + f(x[0], x[1], xd[i], yd[i]) ** 2
    s = s / (1 + x[1] ** 2)
    return s


def EM(x):
    global xd, yd
    n = len(xd)
    ee = np.zeros(n)
    for i in range(n):
        ee[i] = abs(f(x[0], x[1], xd[i], yd[i]))
    s = np.linalg.norm(ee, np.inf)
    return s


def ET(x):
    global xd, yd
    n = len(xd)
    s = 0
    for i in range(n):
        s = s + abs(f(x[0], x[1], xd[i], yd[i]))
    return s


def fig8_10():
    # includes table 8.4

    #  this requires the subaxis.m and parseArgs.m files
    #  note plots change with each run due to the randomly generated data

    global xd, yd

    plt.close("all")
    # get(gcf)
    # MATLAB Position [12 450 748 895]; the final one-panel figure below is
    # scaled to fit the plot window.

    co = np.array([
        [0, 0, 1],
        [0, 0.5, 0],
        [1, 0, 0],
        [0, 0.75, 0.75],
        [0.75, 0, 0.75],
        [0.75, 0.75, 0],
        [0.25, 0.25, 0.25],
    ])

    #  generate data for regression y=a+b*x
    n = 20
    xd = np.random.rand(n)
    qd = 0.7 * (np.random.rand(n) - 0.5)
    a = 1; b = 3
    yd = a + b * xd + 2 * (np.random.rand(n) - 0.5)

    print(np.max(yd))
    print(np.min(yd))

    #  have MATLAB find coeffs for linear fit
    # g=aa+bb*x
    ab = np.polyfit(xd, yd, 1)
    aa = ab[1]; bb = ab[0]
    s2 = np.array([aa, bb])

    # true l^2 fit
    resultD = minimize(
        ED,
        [aa, bb],
        method="Nelder-Mead",
        options={"fatol": np.finfo(float).eps, "xatol": np.finfo(float).eps},
    )
    sD = resultD.x

    # vertical l^1 fit
    resultT = minimize(
        ET,
        [aa, bb],
        method="Nelder-Mead",
        options={"fatol": np.finfo(float).eps, "xatol": np.finfo(float).eps},
    )
    sT = resultT.x

    # vertical l^inf fit
    # [sM, fval, exitflag, output] = fminsearch(@EM, a, optimset('TolFun',eps,'TolX',eps))
    resultM = minimize(
        EM,
        [aa, bb],
        method="Nelder-Mead",
        options={"fatol": np.finfo(float).eps, "xatol": np.finfo(float).eps},
    )
    sM = resultM.x

    # format for book
    print(f"\n E_2 & {s2[0]:3.2f}  & {s2[1]:3.2f}", end="")
    print(f"\n E_D & {sD[0]:3.2f}  & {sD[1]:3.2f}", end="")
    print(f"\n E_1 & {sT[0]:3.2f}  & {sT[1]:3.2f}", end="")
    print(f"\n E_inf & {sM[0]:3.2f}  & {sM[1]:3.2f}", end="")
    print("\n")

    #  evaluate interpolation fn's for use in plotting
    #  nn = number of points plotted
    nn = 30
    X = np.linspace(0, 1, nn)
    Y = np.polyval(ab, X)
    yD = np.zeros(nn)
    yT = np.zeros(nn)
    yM = np.zeros(nn)
    for ix in range(nn):
        yD[ix] = sD[1] * X[ix] + sD[0]
        yT[ix] = sT[1] * X[ix] + sT[0]
        # yll[ix]=solll[0]*X[ix]+solll[1]
        yM[ix] = sM[1] * X[ix] + sM[0]

    # get(gcf)
    # MATLAB Position [25 1115 658 230], scaled to fit the plot window.
    fig, ax = plt.subplots(figsize=(6.58, 3.0), layout="constrained")
    ax.set_prop_cycle(color=co)

    ax.plot(X, yD, "r", linewidth=1.6)
    ax.plot(X, Y, "--k", linewidth=1.6)
    # ax.plot([x1,x4],[y1,y4],'--k',linewidth=1.4)
    ax.plot(xd, yd, ".b", markersize=25)
    ax.grid(True)
    ax.set_xlabel("x-axis", fontsize=14, fontweight="bold")
    ax.set_ylabel("y-axis", fontsize=14, fontweight="bold")
    ax.tick_params(labelsize=14)
    for label in ax.get_xticklabels() + ax.get_yticklabels():
        label.set_fontweight("bold")
    ax.axis([0, 1, 0, 5])
    ax.legend([" E_D", " E_2"], loc="lower right")

    plt.show()


if __name__ == "__main__":
    fig8_10()
