"""Literal Python conversion of fig10_23.m."""

from __future__ import annotations

from pathlib import Path

import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
import numpy as np
import shapefile

from subaxis import subaxis


BASE_DIR = Path(__file__).resolve().parent


def _mapshow(ax: plt.Axes, shape: shapefile.Shape, FaceColor: list[float], LineWidth: float) -> None:
    """Matplotlib/pyshp equivalent of MATLAB mapshow for polygon shapes."""
    points = np.asarray(shape.points)
    parts = list(shape.parts) + [len(points)]
    for part_index in range(len(parts) - 1):
        polygon_points = points[parts[part_index] : parts[part_index + 1]]
        if len(polygon_points) > 0:
            patch = Polygon(
                polygon_points,
                closed=True,
                facecolor=FaceColor,
                edgecolor="black",
                linewidth=LineWidth,
            )
            ax.add_patch(patch)


def fig10_23() -> None:
    # y(t+dt)=Ax(t)

    # measure dt in days
    dt = 7

    # COVID data
    # https://usafacts.org/visualizations/coronavirus-covid-19-spread-map/
    # open Download Data/Known Cases and save as a text file covid_data
    T = np.loadtxt(BASE_DIR / "covid_data3.txt")
    Nd, Md = T.shape

    # NY rows 1865 to 1926
    # April 1, 2020 = column 75

    # T[1864:1926,74:80]
    #
    # input("Press Enter to continue...")

    NY = T[1864:1926, 74:Md].T
    N, M = NY.shape
    print("N, M =", N, M)

    n = N - dt
    m = M
    XD = NY[0:n, :]
    YD = NY[dt : n + dt, :]

    # return

    idata = 1
    if idata == 1:
        # use normalized data
        X = np.zeros((n, m))
        Y = np.zeros((n, m))
        meansX = np.mean(XD, axis=0)
        meansY = np.mean(YD, axis=0)
        for j in range(m):
            X[:, j] = XD[:, j] - meansX[j] * np.ones(n)
            Y[:, j] = YD[:, j] - meansY[j] * np.ones(n)
        XM = np.zeros(m)
        YM = np.zeros(m)
        for j in range(m):
            XM[j] = np.max(np.abs(X[:, j]))
            YM[j] = np.max(np.abs(Y[:, j]))
            X[:, j] = X[:, j] / XM[j]
            Y[:, j] = Y[:, j] / YM[j]
    elif idata == 2:
        # use unnormalized data
        X = XD
        Y = YD
        meansX = np.zeros(m)
        meansY = meansX
        XM = np.ones(m)
        YM = np.ones(m)
    else:
        raise ValueError("idata must equal 1 or 2")

    G = X.T @ X
    GI = np.linalg.inv(G)
    P = Y.T @ X @ GI
    SSx = np.diag(1 / XM)
    Sy = np.diag(YM)
    A = Sy @ P @ SSx
    b = meansY.T - A @ meansX.T
    # np.linalg.norm(b,np.inf)

    # Albany = 1
    # Erie = 15
    # Rensselaer = 42
    # ulster = 56
    # saratoga=46, king=24, westchester=60
    r = 1
    print("r =", r)
    r_index = r - 1

    print("P(r,r) =", P[r_index, r_index])
    ss = np.sort(P[r_index, :])
    print("ss =", ss)
    for ia in range(1, len(ss) + 1):
        if P[r_index, ia - 1] == ss[0]:
            print([ia, P[r_index, ia - 1]])
        elif P[r_index, ia - 1] == ss[len(ss) - 1]:
            print([ia, P[r_index, ia - 1]])

    # NY state GIS data
    # http://gis.ny.gov/gisdata/inventories/details.cfm?DSID=927
    # Download NYS Civil Boundaries (Shape)
    # data = shapefile.Reader(BASE_DIR / 'Counties.shp').shapes()
    shoreline_path = BASE_DIR / "Counties_Shoreline.shp"
    if not shoreline_path.exists():
        raise FileNotFoundError(
            "fig10_23.py requires Counties_Shoreline.shp and its companion "
            ".shx/.dbf files in the same folder. These GIS files were referenced "
            "by the MATLAB code but were not included with the uploaded files."
        )
    data = shapefile.Reader(str(shoreline_path)).shapes()
    # mapshow(data)

    data2 = list(data)
    data2[22] = data[61]
    data2[23:26] = data[22:25]
    data2[26] = data[59]
    data2[27:32] = data[25:30]
    data2[32] = data[60]
    data2[33:44] = data[30:41]
    data2[44] = data[46]
    data2[45:50] = data[41:46]
    data2[50:62] = data[47:59]

    print("data2[30,:] =", data2[29])
    print("data2[52,:] =", data2[51])

    fig1 = plt.figure(1, figsize=(7.2, 6.0))
    fig1.clf()
    # MATLAB figure position was [3 925 560 420]; figsize is scaled to fit the plot window.
    ax1 = subaxis(1, 1, 1, 1, "MT", -0.01, "MB", -0.02, "MR", -0.1, "ML", -0.1, "P", 0.04)

    pm = np.min(P[r_index, :])
    pM = np.max(P[r_index, :])
    print("pm =", pm)
    print("pM =", pM)
    for ic in range(1, 63):
        # print(ic)
        # print(data2[ic-1])
        if P[r_index, ic - 1] > 0:
            gg = 1 - P[r_index, ic - 1] / pM
            bb = gg
            if ic == r:
                _mapshow(ax1, data2[ic - 1], FaceColor=[1, gg, bb], LineWidth=4)
            else:
                _mapshow(ax1, data2[ic - 1], FaceColor=[1, gg, bb], LineWidth=1)
        else:
            rr = 1 - P[r_index, ic - 1] / pm
            gg = rr
            if ic == r:
                _mapshow(ax1, data2[ic - 1], FaceColor=[rr, gg, 1], LineWidth=4)
            else:
                _mapshow(ax1, data2[ic - 1], FaceColor=[rr, gg, 1], LineWidth=1)
        # input("Press Enter to continue...")

    ax1.autoscale_view()
    ax1.set_aspect("equal", adjustable="box")
    ax1.axis("off")
    ax1.set_autoscale_on(False)
    xl = ax1.get_xlim()
    yl = ax1.get_ylim()
    print("xl =", xl)
    print("yl =", yl)
    print("xl[0] =", xl[0])
    say = f"d = {dt:d}"
    ax1.text(
        xl[0] + 0.1 * (xl[1] - xl[0]),
        yl[0] + 0.9 * (yl[1] - yl[0]),
        say,
        fontsize=28,  # MATLAB FontSize 30 reduced by 2
        fontweight="bold",
    )

    plt.show()


if __name__ == "__main__":
    fig10_23()
