"""The feature matrix of one weather year, summarized by a handful of
numbers: its shape, its rank, its singular values and its condition number.

Purpose
-------
Numerical companion to the glossary entry 'featuremtx'.  The GeoSphere
Austria weather station Krems (station id 3805) records eight
measurements per day: minimum, maximum and mean air temperature,
precipitation, sunshine duration, relative humidity, air pressure and
wind speed.  This script downloads the records for 2024 from the
GeoSphere data hub (dataset klima-v2-1d) and writes them to
featuremtx_weather.csv.  A value of -1 for precipitation marks a trace
of rain too small to record and is set to 0.

Collecting the eight measurements of each day as a row gives the feature
matrix of the year.  Every measurement is scaled to zero sample mean and
unit sample variance first, because temperatures, precipitation and
pressure carry different units and the singular values of an unscaled
matrix would report the choice of units rather than the data.

The demo checks the claims the entry makes.  (1) The shape is the pair
(number of data points, number of features), and the whole year is 366
times 8 numbers while the shape is two of them.  (2) The singular values
are the square roots of the eigenvalues of X^T X, and X^T X divided by
the number of data points is the sample covariance matrix of the scaled
measurements, so the singular values carry the spectrum that PCA reads.
(3) The sum of the squared singular values is the squared Frobenius norm
of the feature matrix, which for scaled measurements is the number of
data points times the number of features.  (4) The condition number is
the ratio of the largest to the smallest singular value.  (5) The
smallest singular value is far below the largest because the mean
temperature of a day is close to the midpoint of its minimum and its
maximum, a near-dependency among three columns; dropping the mean
temperature raises the smallest singular value by more than an order of
magnitude.

Deterministic: the data are a fixed archive year and the decomposition
is a singular value decomposition.  Self-contained: numpy + matplotlib
only (stdlib urllib for the download).

Blocks
------
[B-data]      Download the 366 days with eight measurements each and
              scale every measurement to zero sample mean and unit
              sample variance.
[B-shape]     The shape of the feature matrix, and how many numbers it
              takes to state it.
[B-spectrum]  The singular values, their link to the eigenvalues of
              X^T X and to the sample covariance matrix, and the share
              of the squared Frobenius norm that the leading ones carry.
[B-cond]      The rank and the condition number, and the near-dependency
              among the three temperature columns that makes the
              smallest singular value small.
[B-fig]       The preview figure: the spectrum as a stem plot and the
              share of the squared Frobenius norm that k singular values
              carry.

Outputs
-------
featuremtx_weather.csv   : date and the eight measurements, 366 days
featuremtx_spectrum.csv  : index, singular value, cumulative share
featuremtx.png           : preview figure
"""
# /// script
# requires-python = ">=3.10"
# dependencies = ["numpy==2.3.4", "matplotlib==3.10.7"]
# ///


import json
import urllib.request
from pathlib import Path

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt                        # noqa: E402
import numpy as np                                     # noqa: E402

OUT_DIR = Path(__file__).parent

FAILED = []


def check(name, ok):
    FAILED.append(name) if not ok else None
    print(f"  [{'ok' if ok else 'FAIL'}] {name}")


# ---- [B-data] the measurements, scaled to comparable units
PARAMS = ["tlmin", "tlmax", "tl_mittel", "rr", "so_h", "rf_mittel",
          "p_mittel", "vv_mittel"]
URL = ("https://dataset.api.hub.geosphere.at/v1/station/historical/"
       f"klima-v2-1d?parameters={','.join(PARAMS)}&station_ids=3805"
       "&start=2024-01-01&end=2024-12-31")
with urllib.request.urlopen(URL, timeout=180) as resp:
    payload = json.load(resp)
params = payload["features"][0]["properties"]["parameters"]
stamps = [t[:10] for t in payload["timestamps"]]
RAW = np.stack([np.array(params[p]["data"], dtype=float) for p in PARAMS], 1)
RAW[:, 3] = np.maximum(RAW[:, 3], 0.0)       # -1 marks a trace of rain
with open(OUT_DIR / "featuremtx_weather.csv", "w") as f:
    f.write("date," + ",".join(PARAMS) + "\n")
    for day, row in zip(stamps, RAW):
        f.write(day + "," + ",".join(f"{v:g}" for v in row) + "\n")

X = (RAW - RAW.mean(0)) / RAW.std(0)         # zero mean, unit variance
print("[B-data] the feature matrix of one weather year")
check("[B-data] every measurement has zero sample mean and unit variance",
      np.allclose(X.mean(0), 0.0, atol=1e-12)
      and np.allclose(X.std(0), 1.0, atol=1e-12))

# ---- [B-shape] the shape is two numbers
m, d = X.shape
print(f"[B-shape] shape ({m}, {d}): {m} data points, {d} features")
print(f"[B-shape] {m * d} numbers in the matrix, 2 in its shape")
check("[B-shape] 366 days with eight measurements each",
      (m, d) == (366, 8))
check("[B-shape] more data points than features", m > d)

# ---- [B-spectrum] the singular values
sv = np.linalg.svd(X, compute_uv=False)
eigs = np.linalg.eigvalsh(X.T @ X)[::-1]
cum = np.cumsum(sv ** 2) / np.sum(sv ** 2)
print("[B-spectrum] singular values: "
      + ", ".join(f"{s:.2f}" for s in sv))
print(f"[B-spectrum] two of them carry {100 * cum[1]:.0f}% of the "
      f"squared Frobenius norm, four carry {100 * cum[3]:.0f}%")
check("[B-spectrum] the squared singular values are the eigenvalues "
      "of X^T X", np.allclose(sv ** 2, eigs, rtol=1e-9, atol=1e-9))
samplecov = (X.T @ X) / m
check("[B-spectrum] X^T X divided by the number of data points is the "
      "sample covariance matrix of the scaled measurements",
      np.allclose(samplecov, np.cov(X, rowvar=False, bias=True),
                  atol=1e-12))
check("[B-spectrum] the squared singular values sum to the squared "
      "Frobenius norm, which is m times d here",
      abs(np.sum(sv ** 2) - m * d) < 1e-6)
with open(OUT_DIR / "featuremtx_spectrum.csv", "w") as f:
    f.write("j,sigma,share\n")
    for j, (s, c) in enumerate(zip(sv, cum), start=1):
        f.write(f"{j},{s:.6f},{c:.6f}\n")

# ---- [B-cond] the rank and the condition number
rank = int(np.linalg.matrix_rank(X))
cond = sv[0] / sv[-1]
print(f"[B-cond] rank {rank}, condition number {cond:.1f}")
check("[B-cond] the rank is the number of features, so no measurement "
      "is an exact combination of the others", rank == d)
check("[B-cond] the condition number is the ratio of the largest to "
      "the smallest singular value",
      abs(cond - np.linalg.cond(X)) < 1e-6)
# the mean temperature is close to the midpoint of the extremes
mid = 0.5 * (RAW[:, 0] + RAW[:, 1])
gap = np.abs(RAW[:, 2] - mid)
print(f"[B-cond] mean temperature differs from the midpoint of the "
      f"extremes by {gap.mean():.2f} degrees on average")
keep = [i for i in range(d) if i != 2]            # drop mean temperature
sv_keep = np.linalg.svd(X[:, keep], compute_uv=False)
cond_keep = sv_keep[0] / sv_keep[-1]
print(f"[B-cond] dropping the mean temperature: smallest singular value "
      f"{sv[-1]:.2f} -> {sv_keep[-1]:.2f}, condition number "
      f"{cond:.1f} -> {cond_keep:.1f}")
check("[B-cond] dropping the mean temperature raises the smallest "
      "singular value by more than an order of magnitude",
      sv_keep[-1] > 10 * sv[-1])
check("[B-cond] and lowers the condition number by more than an order "
      "of magnitude", cond_keep < cond / 10)

# ---- [B-fig] the preview figure
fig, (axL, axR) = plt.subplots(1, 2, figsize=(9.4, 3.5))
axL.stem(np.arange(1, d + 1), sv, linefmt="k-", markerfmt="ko",
         basefmt=" ")
axL.set_xlabel("index $j$ of the singular value")
axL.set_ylabel("singular value $\\sigma_j$")
axL.set_title("Spectrum of the feature matrix\n(366 days, 8 scaled "
              "measurements)", fontsize=10)
axL.set_xticks(np.arange(1, d + 1))
axR.plot(np.arange(1, d + 1), 100 * cum, "k-o")
axR.axhline(90, color="gray", linestyle="--", label="90 percent")
axR.set_xlabel("number $k$ of leading singular values kept")
axR.set_ylabel("percent of the squared Frobenius norm")
axR.set_title("Share the leading singular values carry", fontsize=10)
axR.set_xticks(np.arange(1, d + 1))
axR.set_ylim(0, 105)
axR.legend(frameon=False)
fig.tight_layout()
fig.savefig(OUT_DIR / "featuremtx.png", dpi=110)
check("[B-fig] the preview figure was written",
      (OUT_DIR / "featuremtx.png").exists())

print()
if FAILED:
    print(f"{len(FAILED)} check(s) FAILED: " + "; ".join(FAILED))
    raise SystemExit(1)
print("all checks passed")