"""Literal Python conversion of table10_6.m."""

from __future__ import annotations

from pathlib import Path

import numpy as np


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


def table10_6() -> None:
    # comparisons with persistence model (P=I) and P^d model

    # y=Ax(t) where y approximates x(t+d)

    # measure dt in days

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

    for it in range(1, 8):
        dt = it

        # XD = new cases; YD = new cases at t+dt
        XD = NY[1 : N - dt, :] - NY[0 : N - dt - 1, :]
        YD = NY[dt + 1 : N, :] - NY[dt : N - 1, :]
        n, m = XD.shape

        meansX = np.mean(XD, axis=0)
        meansY = np.mean(YD, axis=0)
        X = XD - meansX
        Y = YD - meansY
        XM = np.zeros(m)
        YM = np.zeros(m)
        for j in range(m):
            XM[j] = np.linalg.norm(X[:, j], np.inf)
            YM[j] = np.linalg.norm(Y[:, j], np.inf)
            X[:, j] = X[:, j] / XM[j]
            Y[:, j] = Y[:, j] / YM[j]

        # use the pseudo-inverse
        # G=X.T@X
        # GI=np.linalg.inv(G)
        # P=Y.T@X@GI

        # use the SVD
        Uq, singular_values, Vqh = np.linalg.svd(X, full_matrices=True)
        Sq = np.zeros((n, m))
        Sq[:m, :m] = np.diag(singular_values)
        Vq = Vqh.T
        SS = np.zeros((n, m))
        for i in range(m):
            SS[i, i] = 1 / Sq[i, i]
        P = Y.T @ Uq @ SS @ Vq.T

        SSx = np.diag(1 / XM)
        Sy = np.diag(YM)
        A = Sy @ P @ SSx
        b = meansY.T - A @ meansX.T

        if it == 1:
            P1 = P
            Sx1 = SSx
            Sy1 = Sy

        # Albany = 1
        # Erie = 15
        # Rensselaer = 42
        # ulster = 56
        # saratoga=46, king=24, westchester=60

        # check error for county r
        # r=15
        # PP=np.linalg.matrix_power(P1,it)
        # yp=np.zeros(n)
        # ypp=yp.copy()
        # for i in range(n):
        #     yp[i]=np.dot(P[r-1,:],X[i,:])
        #     ypp[i]=np.dot(PP[r-1,:],X[i,:])
        #
        # err1=np.linalg.norm(yp-Y[:,r-1],np.inf)
        # error in y=Px model
        # err2=np.linalg.norm(yp-Y[:,r-1],2)/n
        # error in P=I model
        # err3=np.linalg.norm(Y[:,r-1]-X[:,r-1],2)/n
        # error in P=P1^d model
        # err4=np.linalg.norm(ypp-Y[:,r-1],2)/n

        # use entire state
        PP = np.linalg.matrix_power(P1, it)
        err2 = np.linalg.norm(Y.T - P @ X.T, "fro") / n
        err3 = np.linalg.norm(Y - X, "fro") / n
        err4 = np.linalg.norm(Y.T - PP @ X.T, "fro") / n

        if it == 1:
            print(" \n  d Err = ||y-x(t+d)||     Err2 = ||y-x||/Err    Err3 = ||(P^d)x -x(t+d)||/Err \n ")
        print(f" \\, {dt:2d} \\,\\, &   {err2:6.1e}  &   {err3 / err2:6.1f}  &   {err4 / err2:6.1f} \\\\\\hline ")

    print("\n")


if __name__ == "__main__":
    table10_6()
