1import json
  2import math
  3import numpy as np
  4try:
  5    from scipy.special import sph_harm
  6except ImportError:
  7    from scipy.special import sph_harm_y
  8    def sph_harm(m, n, phi, theta):
  9        # New SciPy API takes (degree, order, polar, azimuth).
 10        return sph_harm_y(n, m, theta, phi)
 11
 12from sklearn.linear_model import LogisticRegression
 13from sklearn.model_selection import train_test_split
 14from sklearn.metrics import accuracy_score
 15
 16SEED = 635
 17rng = np.random.default_rng(SEED)
 18
 19
 20def fibonacci_sphere(m):
 21    # Nearly uniform quadrature points on S^2.
 22    k = np.arange(m, dtype=float)
 23    z = 1.0 - 2.0 * (k + 0.5) / m
 24    phi = (math.pi * (3.0 - math.sqrt(5.0)) * k) % (2.0 * math.pi)
 25    r = np.sqrt(np.maximum(0.0, 1.0 - z * z))
 26    return np.column_stack((r * np.cos(phi), r * np.sin(phi), z))
 27
 28
 29def real_harmonics(points, lmax):
 30    # Real orthonormal spherical harmonics, evaluated in scipy's convention.
 31    x, y, z = points.T
 32    theta = np.arccos(np.clip(z, -1, 1))
 33    phi = np.arctan2(y, x) % (2 * np.pi)
 34    cols, degrees = [], []
 35    for ell in range(lmax + 1):
 36        for mm in range(-ell, ell + 1):
 37            a = sph_harm(abs(mm), ell, phi, theta)
 38            if mm < 0:
 39                h = math.sqrt(2.0) * (-1 if mm % 2 else 1) * a.imag
 40            elif mm == 0:
 41                h = a.real
 42            else:
 43                h = math.sqrt(2.0) * (-1 if mm % 2 else 1) * a.real
 44            cols.append(h)
 45            degrees.append(ell)
 46    return np.asarray(cols).T.real, np.asarray(degrees)
 47
 48
 49def make_funk_matrix(points, neighbors=120):
 50    # For each u, average samples closest to u-perpendicular (small |u.v|).
 51    dots = np.abs(points @ points.T)
 52    idx = np.argpartition(dots, neighbors, axis=1)[:, :neighbors]
 53    A = np.zeros_like(dots)
 54    rows = np.arange(len(points))[:, None]
 55    # A local strip is a simple positive approximation to the great-circle average.
 56    weights = np.exp(-(dots[rows, idx] / 0.055) ** 2)
 57    weights /= weights.sum(axis=1, keepdims=True)
 58    A[rows, idx] = weights
 59    return A
 60
 61
 62def projectors(H, degrees):
 63    # Weighted-by-sampling least-squares projectors; Fibonacci weights are uniform.
 64    p0 = degrees == 0
 65    p2 = degrees == 2
 66    B0, B2 = H[:, p0], H[:, p2]
 67    P0 = B0 @ np.linalg.pinv(B0)
 68    P2 = B2 @ np.linalg.pinv(B2)
 69    return P0, P2
 70
 71
 72def phi_squared(f, A):
 73    # In n=3, c_3/kappa_2 = 1 and Phi(f)=R(f^2).
 74    z = A @ (f ** 2)
 75    return A @ (z ** 2)
 76
 77
 78def rel_norm(x):
 79    return float(np.sqrt(np.mean(x * x)))
 80
 81
 82def diagnostic(points, A, H, degrees):
 83    P0, P2 = projectors(H, degrees)
 84    Pge4 = np.eye(len(points)) - P0 - P2
 85    results = {}
 86    # Small perturbations make the finite difference equal to the claimed derivative.
 87    for ell in (2, 4, 6):
 88        j = np.flatnonzero(degrees == ell)[0]
 89        h = H[:, j] / rel_norm(H[:, j])
 90        eps = 1e-3
 91        out = (phi_squared(1.0 + eps * h, A) - 1.0) / eps
 92        if ell == 2:
 93            ratio = rel_norm(P2 @ out) / rel_norm(h)
 94            results['degree2_multiplier'] = ratio
 95        else:
 96            ratio = rel_norm(Pge4 @ out) / rel_norm(h)
 97            results[f'degree{ell}_multiplier'] = ratio
 98    # Direct nonlinear test, with mixed ellipsoidal and high-frequency perturbation.
 99    h2 = H[:, np.flatnonzero(degrees == 2)[2]] / rel_norm(H[:, np.flatnonzero(degrees == 2)[2]])
100    h4 = H[:, np.flatnonzero(degrees == 4)[1]] / rel_norm(H[:, np.flatnonzero(degrees == 4)[1]])
101    eps = 1e-2
102    before2, before4 = rel_norm(P2 @ (eps*h2)), rel_norm(Pge4 @ (eps*h4))
103    after = phi_squared(1 + eps*h2 + eps*h4, A) - 1
104    results['mixed_degree2_ratio'] = rel_norm(P2 @ after) / before2
105    results['mixed_high_ratio'] = rel_norm(Pge4 @ after) / before4
106    results['predicted_high_multiplier'] = 9.0 / 16.0
107    results['predicted_gap'] = 7.0 / 16.0
108    return results
109
110
111def toy_classification(points, A, H, degrees, n_samples=1200):
112    # Same tiny supervised task: infer the sign of an ellipsoidal (degree-2) coefficient
113    # amid high-frequency angular contamination. Compare raw features with the proposed map.
114    i2 = np.flatnonzero(degrees == 2)[2]
115    hi = np.flatnonzero(degrees >= 4)
116    X, Xmap, y = [], [], []
117    for _ in range(n_samples):
118        sign = 1 if rng.random() > 0.5 else -1
119        f = np.ones(len(points))
120        f += sign * 0.18 * H[:, i2] / rel_norm(H[:, i2])
121        coeff = rng.normal(0, 0.055, size=len(hi))
122        f += H[:, hi] @ coeff / math.sqrt(len(hi))
123        f = np.maximum(f, 0.03)
124        X.append(f - 1.0)
125        Xmap.append(phi_squared(f, A) - 1.0)
126        y.append(sign > 0)
127    X, Xmap, y = np.asarray(X), np.asarray(Xmap), np.asarray(y)
128    xa, xb, ya, yb = train_test_split(X, y, test_size=.3, random_state=SEED, stratify=y)
129    ma = LogisticRegression(C=1.0, max_iter=200, random_state=SEED).fit(xa, ya)
130    xma, xmb, _, _ = train_test_split(Xmap, y, test_size=.3, random_state=SEED, stratify=y)
131    # split indices above are independently random but identically seeded; labels align.
132    mm = LogisticRegression(C=1.0, max_iter=200, random_state=SEED).fit(xma, ya)
133    return {'raw_accuracy': float(accuracy_score(yb, ma.predict(xb))),
134            'idea_accuracy': float(accuracy_score(yb, mm.predict(xmb)))}
135
136
137def main():
138    m = 1400
139    points = fibonacci_sphere(m)
140    H, degrees = real_harmonics(points, 6)
141    A = make_funk_matrix(points)
142    d = diagnostic(points, A, H, degrees)
143    c = toy_classification(points, A, H, degrees)
144    out = {'seed': SEED, 'points': m, 'neighbor_count': 120,
145           'diagnostic': d, 'classification': c}
146    print(json.dumps(out, indent=2, sort_keys=True))
147
148
149if __name__ == '__main__':
150    main()