"""
innerproduct.py — numerical companion to the glossary entry 'inner product'.

Purpose
-------
Verifies the entry's axioms, geometry, and ML tie-ins for the standard
Euclidean inner product <x, x'> = x^T x' on R^d, and for the implicit
inner products encoded by a Gaussian kernel.  Self-contained (numpy/
matplotlib only; the optional network check inside [P-cities] uses the
stdlib urllib to query OpenStreetMap and is skipped without network),
fixed seed.

Blocks
------
One block per paragraph of the entry (marked [P...]), in the order the
paragraphs appear: each block verifies numerically what its paragraph asserts.

[P-array]   The opening framing: a colour image and a sampled sensor signal
            both flatten into the numeric array that a data point carries as
            its feature vector.
[P-axioms]  The three defining properties hold for 1000 random triples and
            scalars (d = 5); the complex case satisfies conjugate symmetry,
            linearity in the first argument and positive-definiteness; and the
            E{x y} inner product on zero-expectation RVs is another instance,
            where the cosine of the angle is the correlation coefficient.
[P-norm]    The induced norm sqrt(<x, x>) equals the Euclidean norm and
            satisfies the triangle inequality; the induced metric satisfies
            the metric axioms.
[P-cos]     <x, x'> = ||x|| ||x'|| cos(theta) in R^2, and the Cauchy-Schwarz
            bound holds for 1000 random pairs.
[P-project] The closest point of a subspace: the error is orthogonal to every
            vector of the subspace, and the Pythagorean identity holds.
[P-convex]  Projection onto a closed convex set satisfies the variational
            inequality <v - vhat, u - vhat> <= 0 and is the nearest point.
[P-linreg]  Least squares: the error is orthogonal to every column of the
            feature matrix and the predictions are the orthogonal projection
            of the label vector onto its column space; an orthogonal Q leaves
            inner products, ERM values and predictions unchanged.
[P-weight]  The weighted inner product x^T A x' satisfies the axioms and
            CHANGES the similarity ranking: for A = diag(4, 1/4) the candidate
            preferred under the standard inner product loses under it.
[P-basis]   Coordinates with respect to an orthonormal basis are inner
            products; declaring an arbitrary basis orthonormal defines an
            inner product, which equals the weighted one with A = (B B^T)^-1.
[P-cities]  The geometry of the entry's city figure: five European cities as
            vectors in R^2, drawn in an orthonormal frame that puts Helsinki
            on the horizontal line, with the dashed rulers hitting it at the
            projected lengths; a neuron's pre-activation is maximized by the
            input aligned with its weights; and the pinned coordinates are
            re-checked against OpenStreetMap when the network is available.
[P-kernel]  The Gram matrix of a Gaussian kernel on 40 points is positive
            semi-definite: kernel evaluations are inner products between
            transformed feature vectors.
[P-dist]    The squared Euclidean distance expands into three inner products,
            which is why distance-based methods touch the feature vectors only
            through inner products.

Outputs
-------
innerproduct.png : matplotlib preview — scatter of <x, x'> against
                   ||x|| ||x'|| cos(theta) on the identity line
                   (checking only; the entry's figure is schematic TikZ).
"""

import numpy as np
import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt

report = []


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


rng = np.random.default_rng(0)

# ----------------------------------------------------------- [P-array]
# a colour image and a sampled sensor signal both flatten to feature vectors
img = rng.random((4, 5, 3))                        # 4 x 5 pixels, 3 channels
x_img = img.reshape(-1)                            # the numeric array of features
tt = np.linspace(0.0, 1.0, 64, endpoint=False)
x_sig = np.sin(2 * np.pi * 5 * tt)                 # 64 amplitudes of a signal
ok_arr = (x_img.shape == (4 * 5 * 3,)
          and np.isclose(x_img[0], img[0, 0, 0])
          and np.isclose(x_img[-1], img[-1, -1, -1])
          and x_sig.shape == (64,))
check("[P-array]   image pixels and signal amplitudes flatten to a feature "
      "vector", ok_arr)

# ---------------------------------------------------------- [P-axioms]
ok_sym = ok_lin = ok_pd = True
for _ in range(1000):
    x, xp, xpp = rng.standard_normal((3, 5))
    b = rng.standard_normal()
    ok_sym &= np.isclose(x @ xp, xp @ x)
    ok_lin &= np.isclose((b * x + xpp) @ xp, b * (x @ xp) + xpp @ xp)
    ok_pd &= x @ x >= 0
ok_pd &= np.isclose(np.zeros(5) @ np.zeros(5), 0.0)
check("[P-axioms] symmetry, linearity, positive-definiteness",
      ok_sym and ok_lin and ok_pd)

# the complex case <x, x'> = sum_j x_j conj(x'_j): conjugate symmetry,
# linearity in the FIRST argument, positive-definiteness (as in the entry)
def cip(x, xp):
    return np.sum(x * np.conj(xp))


ok_csym = ok_clin = ok_cpd = True
for _ in range(1000):
    z, zp, zpp = (rng.standard_normal((3, 5))
                  + 1j * rng.standard_normal((3, 5)))
    b = rng.standard_normal() + 1j * rng.standard_normal()
    ok_csym &= np.isclose(cip(z, zp), np.conj(cip(zp, z)))
    ok_clin &= np.isclose(cip(b * z + zpp, zp), b * cip(z, zp) + cip(zpp, zp))
    ok_cpd &= (np.isclose(cip(z, z).imag, 0.0) and cip(z, z).real >= 0)
ok_cpd &= np.isclose(cip(np.zeros(5), np.zeros(5)), 0.0)
check("[P-axioms] conjugate symmetry, first-argument linearity, pd over C",
      ok_csym and ok_clin and ok_cpd)

pw = np.array([0.10, 0.15, 0.20, 0.25, 0.18, 0.12])  # outcome probabilities


def Ew(v):                                         # expectation
    return pw @ v


xr, yr = rng.standard_normal((2, 6)) * 1.5
xr, yr = xr - Ew(xr), yr - Ew(yr)                  # zero-expectation RVs
ipE = Ew(xr * yr)                                  # <x, y> = E{x y}
cos_ang = ipE / np.sqrt(Ew(xr * xr) * Ew(yr * yr))  # cosine of the angle
# correlation from its defining ratio (means subtracted explicitly)
cov_xy = Ew(xr * yr) - Ew(xr) * Ew(yr)
var_x = Ew(xr * xr) - Ew(xr) ** 2
var_y = Ew(yr * yr) - Ew(yr) ** 2
corr = cov_xy / np.sqrt(var_x * var_y)             # correlation coefficient
ok_corr = np.isclose(corr, cos_ang)
# uncorrelated = orthogonal: subtract the projection to decorrelate
yo = yr - (ipE / Ew(xr * xr)) * xr
check("[P-axioms]   correlation = cosine of the angle; uncorrelated = "
      "orthogonal", ok_corr and np.isclose(Ew(xr * yo), 0.0))

# ------------------------------------------------------------ [P-norm]
ok_norm = True
for _ in range(1000):
    x, xp, xpp = rng.standard_normal((3, 5))
    ok_norm &= np.isclose(np.sqrt(x @ x), np.linalg.norm(x))
    ok_norm &= np.linalg.norm(x + xp) <= np.linalg.norm(x) \
        + np.linalg.norm(xp) + 1e-12
    # metric axioms: symmetry and triangle inequality
    ok_norm &= np.isclose(np.linalg.norm(x - xp), np.linalg.norm(xp - x))
    ok_norm &= np.linalg.norm(x - xpp) <= np.linalg.norm(x - xp) \
        + np.linalg.norm(xp - xpp) + 1e-12
check("[P-norm]   induced norm and metric", ok_norm)

# ------------------------------------------------------------- [P-cos]
ok_cos = ok_cs = True
ips, cos_prods = [], []
for _ in range(1000):
    x, xp = rng.standard_normal((2, 2))
    theta = np.arctan2(xp[1], xp[0]) - np.arctan2(x[1], x[0])
    lhs = x @ xp
    rhs = np.linalg.norm(x) * np.linalg.norm(xp) * np.cos(theta)
    ok_cos &= np.isclose(lhs, rhs)
    ok_cs &= abs(lhs) <= np.linalg.norm(x) * np.linalg.norm(xp) + 1e-12
    ips.append(lhs), cos_prods.append(rhs)
check("[P-cos]    cos-theta identity and Cauchy-Schwarz bound",
      ok_cos and ok_cs)

# --------------------------------------------------------- [P-project]
# the closest point of a subspace: orthogonality condition and Pythagoras
Us = rng.standard_normal((6, 2))                   # a 2-dim subspace of R^6
vs = rng.standard_normal(6)
coef = np.linalg.lstsq(Us, vs, rcond=None)[0]
vhat = Us @ coef                                   # the closest point
err = vs - vhat
ok_orth = np.allclose(Us.T @ err, 0.0, atol=1e-10)  # error orthogonal to U
u_other = Us @ rng.standard_normal(2)
ok_pyth = np.isclose(np.linalg.norm(vs - u_other) ** 2,
                     np.linalg.norm(err) ** 2
                     + np.linalg.norm(vhat - u_other) ** 2)
check("[P-project] closest point of a subspace: error is orthogonal to it, "
      "and Pythagoras holds", ok_orth and ok_pyth)

# ---------------------------------------------------------- [P-convex]
# projection onto a closed convex set: the variational inequality
vc = np.array([2.5, 1.0])                          # a point outside the unit ball
vhat_c = vc / np.linalg.norm(vc)                   # its projection onto the ball
pts = rng.standard_normal((4000, 2))
pts = pts[np.linalg.norm(pts, axis=1) <= 1.0]      # points of the set C
ok_var = np.all((pts - vhat_c) @ (vc - vhat_c) <= 1e-12)
ok_near = np.all(np.linalg.norm(pts - vc, axis=1)
                 >= np.linalg.norm(vhat_c - vc) - 1e-12)
check("[P-convex]  projection onto a closed convex set satisfies the "
      "variational inequality and is the nearest point", ok_var and ok_near)

# ---------------------------------------------------------- [P-linreg]
# orthogonal transformations preserve inner products, hence the linear-model
# ERM objective values and predictions (the entry's invariance statement)
m = 50
X = rng.standard_normal((m, 3))
y = X @ np.array([1.0, -2.0, 0.5]) + 0.1 * rng.standard_normal(m)
Q, _ = np.linalg.qr(rng.standard_normal((3, 3)))     # a random orthogonal Q
w = rng.standard_normal(3)
ok_ip = np.allclose((X @ Q.T) @ (Q @ w), X @ w)       # <Qx, Qw> = <x, w>
f_orig = np.mean((y - X @ w) ** 2)                    # ERM objective at w
f_rot = np.mean((y - (X @ Q.T) @ (Q @ w)) ** 2)       # after rotating x and w
w_hat = np.linalg.lstsq(X, y, rcond=None)[0]          # ERM solution, original
w_hat_rot = np.linalg.lstsq(X @ Q.T, y, rcond=None)[0]  # ERM solution, rotated
ok_pred = np.allclose((X @ Q.T) @ w_hat_rot, X @ w_hat)  # same predictions
check("[P-linreg]  orthogonal Q preserves inner products, ERM values, "
      "predictions", ok_ip and np.isclose(f_orig, f_rot) and ok_pred)

# the least-squares characterization: the error is orthogonal to every column
# of the feature matrix, so the predictions are the orthogonal projection of
# the label vector onto its column space
resid = y - X @ w_hat
ok_normal = np.allclose(X.T @ resid, 0.0, atol=1e-9)
Pcol = X @ np.linalg.pinv(X)                       # projector onto col(X)
ok_proj = np.allclose(Pcol @ y, X @ w_hat)
check("[P-linreg]  least squares: error orthogonal to every column, "
      "predictions are the projection of the label vector",
      ok_normal and ok_proj)

# ---------------------------------------------------------- [P-weight]
Aw = np.diag([4.0, 0.25])                          # positive definite weights


def ipw(x, yv):                                    # weighted inner product
    return x @ Aw @ yv


xq = np.array([1.0, 1.0])                          # query vector
c1, c2 = np.array([0.9, 1.2]), np.array([1.3, 0.4])  # two candidates
ok_axioms = (np.isclose(ipw(c1, c2), ipw(c2, c1))  # symmetry
             and np.isclose(ipw(2.0 * c1 + c2, xq),
                            2.0 * ipw(c1, xq) + ipw(c2, xq))  # bilinearity
             and ipw(xq, xq) > 0)                  # positive definite
ok_swap = (xq @ c1 > xq @ c2) and (ipw(xq, c1) < ipw(xq, c2))
check("[P-weight] weighted inner product: axioms hold, similarity "
      "ranking swaps", ok_axioms and ok_swap)

# ----------------------------------------------------------- [P-basis]
# (i) coordinates w.r.t. an orthonormal basis = inner products with it
Qb, _ = np.linalg.qr(rng.standard_normal((4, 4)))    # orthonormal basis (cols)
u4 = rng.standard_normal(4)
coords = Qb.T @ u4                                   # <u, b^(j)> for each j
ok_coord = np.allclose(Qb @ coords, u4)              # u = sum_j coeff_j b^(j)
# (ii) declaring an arbitrary basis orthonormal defines an inner product
Bb = rng.standard_normal((4, 4)) + 4.0 * np.eye(4)   # an arbitrary basis


def ipb(u, v):                                       # via coordinate vectors
    cu, cv = np.linalg.solve(Bb, u), np.linalg.solve(Bb, v)
    return cu @ cv


v4, w4 = rng.standard_normal((2, 4))
a4 = rng.standard_normal()
ok_ax = (np.isclose(ipb(u4, v4), ipb(v4, u4))        # symmetry
         and np.isclose(ipb(a4 * u4 + w4, v4),
                        a4 * ipb(u4, v4) + ipb(w4, v4))  # linearity
         and ipb(u4, u4) > 0)                        # positive definite
G = np.array([[ipb(Bb[:, i], Bb[:, j]) for j in range(4)] for i in range(4)])
ok_onb = np.allclose(G, np.eye(4))                   # B is orthonormal under ipb
Ab = np.linalg.inv(Bb @ Bb.T)                        # the weighted form
ok_wt = np.isclose(ipb(u4, v4), u4 @ Ab @ v4)        # ipb = u^T A v
check("[P-basis]  orthonormal-basis coordinates; declared basis defines "
      "an inner product", ok_coord and ok_ax and ok_onb and ok_wt)

# ---------------------------------------------------------- [P-cities]
# (latitude, longitude) in degrees, retrieved 2026-08-15 from OpenStreetMap
# via the Nominatim geocoder (https://nominatim.openstreetmap.org/search,
# format=json, limit=1, queries "<City>, <Country>"); data (c) OpenStreetMap
# contributors, ODbL.  The optional check below re-verifies these values
# against the live service whenever network access is available.
cities = {
    "Helsinki": np.array([60.1666, 24.9435]),
    "Paris": np.array([48.8535, 2.3484]),
    "Rome": np.array([41.8933, 12.4829]),
    "Madrid": np.array([40.4168, -3.7035]),
    "Reykjavik": np.array([64.1460, -21.9422]),
}
QUERIES = {"Helsinki": "Helsinki, Finland", "Paris": "Paris, France",
           "Rome": "Rome, Italy", "Madrid": "Madrid, Spain",
           "Reykjavik": "Reykjavik, Iceland"}
hel = cities["Helsinki"]
ips_city = {name: float(hel @ v) for name, v in cities.items()
            if name != "Helsinki"}
ok_rank = max(ips_city, key=ips_city.get) == "Reykjavik"
# the figure's frame: each city (x1, x2) is DRAWN at x1 u + x2 v for an
# orthonormal pair u, v chosen such that Helsinki lies on the horizontal
# line and the other cities lie above it
alpha = np.arctan2(hel[1], hel[0])
u_fr = np.array([np.cos(alpha), np.sin(alpha)])
v_fr = np.array([np.sin(alpha), -np.cos(alpha)])
ok_onf = (np.isclose(u_fr @ v_fr, 0.0)              # orthonormal pair
          and np.isclose(u_fr @ u_fr, 1.0) and np.isclose(v_fr @ v_fr, 1.0))
pos = {n: c[0] * u_fr + c[1] * v_fr for n, c in cities.items()}
pos_fig = {"Helsinki": (65.13, 0.00), "Paris": (46.03, 16.54),
           "Rome": (43.48, 4.51), "Madrid": (35.92, 18.90),
           "Reykjavik": (50.85, 44.84)}
ok_pos = all(np.allclose(pos[n], pos_fig[n], atol=0.01) for n in pos_fig)
ok_horiz = np.isclose(pos["Helsinki"][1], 0.0) and pos["Helsinki"][0] > 0
ok_above = all(pos[n][1] > 0 for n in cities if n != "Helsinki")
# orthonormality preserves inner products: drawn positions reproduce the
# lat/lon inner products with Helsinki exactly
ok_pres = all(np.isclose(float(pos[n] @ pos["Helsinki"]), ips_city[n])
              for n in ips_city)
# vertical rulers: each foot sits on the horizontal axis at the projected
# length <hel, x> / ||hel||, so left-to-right order = inner-product order
ok_feet = all(np.isclose(pos[n][0], ips_city[n] / np.linalg.norm(hel))
              for n in ips_city)
ok_order = (sorted(ips_city, key=ips_city.get)
            == sorted(ips_city, key=lambda n: pos[n][0]))
check("[P-cities] orthonormal frame puts Helsinki horizontal, others "
      "above; positions match the figure and preserve inner products",
      ok_onf and ok_pos and ok_horiz and ok_above and ok_pres
      and ok_feet and ok_order and ok_rank)
print("           " + ", ".join(
    f"{n}: <hel,x>={ips_city[n]:.0f}, drawn=({pos[n][0]:.2f},{pos[n][1]:.2f})"
    for n in ips_city))

w_t = np.array([2.0, 1.2])                        # the neuron's weight vector
angles_n = np.linspace(0.0, 2.0 * np.pi, 720, endpoint=False)
inputs = np.c_[np.cos(angles_n), np.sin(angles_n)]  # unit-norm inputs
best_in = inputs[np.argmax(inputs @ w_t)]
ang_n = np.degrees(np.arccos(np.clip(
    best_in @ (w_t / np.linalg.norm(w_t)), -1.0, 1.0)))
check("[P-cities] pre-activation <w, x> maximized by the template "
      "direction", ang_n < 0.5)

# optional: re-verify the pinned coordinates against the live authoritative
# source (OpenStreetMap Nominatim); skipped gracefully without network
try:
    import json
    import time
    import urllib.parse
    import urllib.request

    fetched = {}
    for name, q in QUERIES.items():
        url = ("https://nominatim.openstreetmap.org/search"
               "?format=json&limit=1&q=" + urllib.parse.quote(q))
        req = urllib.request.Request(url, headers={
            "User-Agent": "dictionaryappliedml-demo/1.0 (alex.jung@aalto.fi)"})
        with urllib.request.urlopen(req, timeout=15) as resp:
            hit = json.load(resp)[0]
        fetched[name] = np.array([float(hit["lat"]), float(hit["lon"])])
        time.sleep(1.1)                 # Nominatim usage policy: max 1 req/s
    ok_src = all(np.allclose(cities[n], fetched[n], atol=0.01)
                 for n in cities)       # pinned values current to 0.01 degrees
    check("[P-cities] pinned (lat, lon) values match OpenStreetMap Nominatim",
          ok_src)
except Exception as exc:                # no network: keep the demo runnable
    print(f"  [--] [P-cities] skipped (no network / service unavailable: "
          f"{type(exc).__name__})")

n_ok = sum(ok for _, ok in report)
print(f"\n{n_ok}/{len(report)} checks pass")
print("wrote pythondemos/innerproduct.png")
if n_ok != len(report):
    raise SystemExit(1)

# ---------------------------------------------------------- [P-kernel]
P = rng.standard_normal((40, 2))
d2 = ((P[:, None, :] - P[None, :, :]) ** 2).sum(-1)
K = np.exp(-d2 / 2.0)
eigmin = float(np.linalg.eigvalsh(K).min())
check("[P-kernel]   Gaussian-kernel Gram matrix is psd", eigmin > -1e-10)

# ------------------------------------------------------------ [P-dist]
ok_dist = True
for _ in range(1000):
    x, yv = rng.standard_normal((2, 5))
    ok_dist &= np.isclose(np.linalg.norm(x - yv) ** 2,
                          x @ x - 2 * (x @ yv) + yv @ yv)
check("[P-dist]   ||x - y||^2 = <x,x> - 2<x,y> + <y,y>", ok_dist)

# ---- preview figure (runs last; uses cos_prods and ips from [P-cos])
fig, ax = plt.subplots(figsize=(3.8, 3.8))
ax.plot(cos_prods, ips, "k.", ms=2)
lim = [min(ips), max(ips)]
ax.plot(lim, lim, "k--", lw=0.8)
ax.set_xlabel(r"$\|x\|\,\|x'\|\cos\theta$")
ax.set_ylabel(r"$\langle x, x'\rangle$")
ax.set_aspect("equal")
fig.tight_layout()
fig.savefig("pythondemos/innerproduct.png", dpi=110)