"""Literal Python conversion of fig10_20.m."""

from __future__ import annotations

import matplotlib.pyplot as plt
import numpy as np

from subaxis import subaxis


def _set_bold_tick_font(ax: plt.Axes, fontsize: int) -> None:
    ax.tick_params(labelsize=fontsize)
    ax.xaxis.label.set_size(fontsize)
    ax.yaxis.label.set_size(fontsize)
    ax.xaxis.label.set_fontweight("bold")
    ax.yaxis.label.set_fontweight("bold")
    for label in ax.get_xticklabels() + ax.get_yticklabels():
        label.set_fontweight("bold")


def fig10_20() -> None:
    # modal approximation

    # random generation of A values
    # m=10;   # number of vertices (must be even)
    # A[0:m//2,0:m//2]=10+np.random.randint(1,21,(m//2,m//2))
    # A[m//2:m,m//2:m]=10+np.random.randint(1,21,(m//2,m//2))
    # A[m//2,0]=1
    # A=A/np.sum(A,axis=0)

    # book example
    m = 4
    A = np.zeros((m, m))
    d = 1 / 20
    A[0, 0] = 1 / 4 - d
    A[1, 0] = 3 / 4 - d
    A[2, 0] = 2 * d
    A[0, 1] = 1 / 3
    A[1, 1] = 2 / 3
    A[2, 2] = 1 / 4
    A[3, 2] = 3 / 4
    A[2, 3] = 1 / 3
    A[3, 3] = 2 / 3

    # np.sum(A,axis=0)
    print("A =\n", A)
    eigens = np.linalg.eigvals(A)
    print("eigens =", eigens)
    sings = np.linalg.svd(A, compute_uv=False)
    print("sings =", sings)

    N = m + 1  # number of time steps
    y = np.zeros((N, m))
    y[0, 0] = 100
    # y[0,:]=np.ones(m)
    for i in range(1, N):
        z = A @ y[i - 1, :].T
        y[i, :] = z.T
    # q=0.1*(2*np.random.rand(m)-1)
    q = 0
    yd = (1 + q) * y
    XD = yd[0 : N - 1, :]
    YD = yd[1:N, :]
    n = N - 1

    # use the pseudo-inverse
    # X=XD
    # Y=YD
    # G=X.T@X
    # GI=np.linalg.inv(G)
    # kappa[k-1]=np.linalg.cond(G,np.inf)
    # P=Y.T@X@GI

    # use the SVD
    U, singular_values, Vh = np.linalg.svd(XD, full_matrices=True)
    S = np.zeros((n, m))
    S[:m, :m] = np.diag(singular_values)
    V = Vh.T
    SS = np.zeros((n, m))
    for i in range(m):
        SS[i, i] = 1 / S[i, i]
    P = YD.T @ U @ SS @ V.T
    Ac = P
    error = np.linalg.norm(A - Ac, np.inf) / np.linalg.norm(A, np.inf)
    print("error =", error)

    eigenvalues, VP = np.linalg.eig(P)
    DP = np.diag(eigenvalues)
    I = np.argsort(np.abs(np.diag(DP)))[::-1]
    D0 = np.abs(np.diag(DP))[I]
    D = np.zeros(m, dtype=complex)
    for i in range(m):
        D[i] = DP[I[i], I[i]]
    print("D =", D)
    VP = VP[:, I]
    x = XD[0, :].T
    alpha = np.linalg.solve(VP, x)
    print("alpha =", alpha)
    err1 = np.zeros(10)
    err2 = np.zeros(10)
    err3 = np.zeros(10)
    err4 = np.zeros(10)
    for it in range(1, 11):
        x = Ac @ x
        x1 = alpha[0] * D[0] ** it * VP[:, 0]
        x2 = x1 + alpha[1] * D[1] ** it * VP[:, 1]
        x3 = x2 + alpha[2] * D[2] ** it * VP[:, 2]
        x4 = x3 + alpha[3] * D[3] ** it * VP[:, 3]
        err1[it - 1] = np.linalg.norm(x - x1, np.inf) / np.linalg.norm(x)
        err2[it - 1] = np.linalg.norm(x - x2, np.inf) / np.linalg.norm(x)
        err3[it - 1] = np.linalg.norm(x - x3, np.inf) / np.linalg.norm(x)
        err4[it - 1] = np.linalg.norm(x - x4, np.inf) / np.linalg.norm(x)
    err4_error = np.linalg.norm(err4, np.inf)
    print("err4_error =", err4_error)

    fig1 = plt.figure(1, figsize=(7.2, 3.4))
    fig1.clf()
    # MATLAB figure position was [25 1115 658 230]; figsize is scaled to fit the plot window.
    ax1 = subaxis(1, 1, 1, 1, "MT", 0.01, "MB", 0.18, "MR", -0.02, "ML", 0.06, "P", 0.04)
    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],
        ]
    )
    ax1.set_prop_cycle(color=co.tolist())
    ax1.semilogy(np.arange(1, 11), err1, "--sk", linewidth=1.5, markersize=8)
    ax1.semilogy(np.arange(1, 11), err2, "--*r", linewidth=1.5, markersize=8)
    ax1.semilogy(np.arange(1, 11), err3, "--ob", linewidth=1.5, markersize=8)
    ax1.set_frame_on(True)
    ax1.grid(True)
    ax1.axis([1, 10, 1e-10, 5])
    # ax1.set_yticks([0,0.05,0.1])
    ax1.set_xlabel("t-axis", fontweight="bold")
    ax1.set_ylabel("Relative Error", fontweight="bold")
    legend = ax1.legend(
        [" M = 1 ", " M = 2 ", " M = 3 "],
        loc="lower left",
        fontsize=12,  # MATLAB FontSize 14 reduced by 2
    )
    for text in legend.get_texts():
        text.set_fontweight("bold")
    _set_bold_tick_font(ax1, 14)  # MATLAB FontSize 16 reduced by 2

    # fig1.savefig('cars2.eps',format='eps',bbox_inches='tight')

    plt.show()


if __name__ == "__main__":
    fig10_20()
