"""Literal Python conversion of fig10_21.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_21() -> None:
    # solve z'=Bz

    n = 6
    m = 4

    S0 = np.random.rand(m, m)
    D0 = np.diag([1, 0.5, -0.5, -1])
    B = S0 @ D0 @ np.linalg.inv(S0)
    eigenvalues_B, VB = np.linalg.eig(B)
    DB = np.diag(eigenvalues_B)

    kk = np.zeros(5)
    norm_error = np.zeros(5)
    for ik in range(1, 6):
        z = np.zeros((n, 4), dtype=complex)
        z[0, :] = [100, 0, 0, 0]
        beta = np.linalg.solve(VB, z[0, :].T)
        k = 10 ** (2 - ik)
        t = 0
        for it in range(1, n):
            t = t + k
            for im in range(m):
                z[it, :] = z[it, :] + beta[im] * np.exp(DB[im, im] * t) * VB[:, im].T

        z = np.real_if_close(z)
        XD = z[0 : n - 1, :]
        YD = z[1:n, :]

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

        BB = (np.eye(m) - np.linalg.inv(P)) / k
        kk[ik - 1] = k
        norm_error[ik - 1] = np.linalg.norm(B - BB, np.inf) / np.linalg.norm(B, np.inf)
        print(f"norm_error[{ik}] =", norm_error[ik - 1])

        eigenvalues_P, VP = np.linalg.eig(P)
        I = np.argsort(eigenvalues_P)[::-1]
        D0 = eigenvalues_P[I]
        print("(1-1./D0)/k =", (1 - 1 / D0) / k)

    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.05, "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.loglog(kk, norm_error, "--or", linewidth=1.5, markersize=8)
    ax1.set_frame_on(True)
    ax1.grid(True)
    ax1.axis([1e-3, 10, 1e-4, 300])
    ax1.set_yticks([1e-4, 1e-2, 1, 100])
    ax1.set_xlabel("k-axis", fontweight="bold")
    ax1.set_ylabel("Relative Error", fontweight="bold")
    # ax1.legend([' k = 1 ',' k = 2 ',' k = 3 '],loc='lower left',fontsize=12)
    _set_bold_tick_font(ax1, 14)  # MATLAB FontSize 16 reduced by 2

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

    plt.show()


if __name__ == "__main__":
    fig10_21()
