#!/usr/bin/env python3 """Finite mapped spectral-collocation reference benchmark, not a QNM classifier. For the black-hole cases the code uses dimensionless radius r=r_phys/(GM), frequency omega=(GM)*omega_phys, and a=a_phys/(GM), with physical c_light=1. For the horizon-free PT control use any reference length L0: omega=L0*omega_phys, alpha=L0*alpha_phys, V0=L0^2*V0_phys. The map parameter c is NOT c_light. The collocation polynomial lives in computational u; x is the compact physical coordinate. All evaluations and derivatives reconstruct this one polynomial. The full declared grid is retained. Reciprocal, ambiguity-checked matching uses frequencies and eigenfunctions without reference frequencies; a later label selects only the finite fundamental-reference benchmark. Endpoint rows, finite residual samples, and matched matrices do not prove a retarded analytic domain or continuum spectral convergence. No Kerr angular eigensolve is run. """ import argparse import csv import hashlib import json from pathlib import Path import platform import numpy as np from scipy.interpolate import BarycentricInterpolator from scipy.fft import dct from scipy.linalg import eig import scipy ORDERS = (12, 20, 24, 28, 32, 36, 40, 44) MAPS = (0.7, 1.0, 1.4) THRESHOLDS = { "frequency_difference": 5e-6, "second_best_gap": 5e-6, "eigenfunction_mismatch": 1e-3, "pencil_residual": 1e-10, "transformed_residual": 1e-4, "unfactored_residual": 1e-4, "horizon_row_residual": 1e-4, "infinity_row_residual": 1e-4, "ingoing_phase_residual": 1e-4, "outgoing_phase_residual": 1e-4, "coefficient_tail": 1e-4, "reference_difference": 5e-6, } GAUSS_NODES, GAUSS_WEIGHTS = np.polynomial.legendre.leggauss(64) COMPARE_POINTS = (GAUSS_NODES+1)/2 COMPARE_WEIGHTS = GAUSS_WEIGHTS/2 REFERENCES = { "poschl-teller": 0.994987437106620 - 0.100000000000000j, "schwarzschild": 0.3736716844180418358 - 0.0889623156889356983j, "kerr-radial": 0.4079791087153231765 - 0.0872571447358583068j, } def chebyshev_lobatto(order): j = np.arange(order + 1) z = np.cos(np.pi * j / order) weights = np.ones(order + 1) weights[[0, -1]] = 2 weights *= (-1.0) ** j zz = z[:, None] dz = (weights[:, None] / weights[None, :]) / ( zz - zz.T + np.eye(order + 1) ) dz -= np.diag(dz.sum(axis=1)) return (z + 1) / 2, 2 * dz def mapped_derivatives(order, map_parameter): u, du = chebyshev_lobatto(order) denominator = map_parameter + (1 - map_parameter) * u x = u / denominator x_u = map_parameter / denominator**2 x_uu = -2 * map_parameter * (1 - map_parameter) / denominator**3 dx = np.diag(1 / x_u) @ du d2x = ( np.diag(1 / x_u**2) @ (du @ du) - np.diag(x_uu / x_u**3) @ du ) return u, x, dx, d2x def quadratic_pencil(c0, c1, c2, dx, d2x): matrices = [ np.diag(c2[q]) @ d2x + np.diag(c1[q]) @ dx + np.diag(c0[q]) for q in range(3) ] size = matrices[0].shape[0] zero = np.zeros_like(matrices[0]) identity = np.eye(size, dtype=complex) left = np.block([[zero, identity], [-matrices[0], -matrices[1]]]) right = np.block([[identity, zero], [zero, matrices[2]]]) values, vectors = eig(left, right) return values, vectors[:size], matrices def poschl_teller(order, map_parameter): u, x, dx, d2x = mapped_derivatives(order, map_parameter) alpha, potential = 0.2, 1.0 zero = np.zeros_like(x) c2 = [alpha**2 * x * (1 - x), zero, zero] c1 = [alpha**2 * (1 - 2 * x), 1j * alpha * (2 * x - 1), zero] c0 = [-potential + zero, 1j * alpha + zero, 1 + zero] return u, x, *quadratic_pencil(c0, c1, c2, dx, d2x) def schwarzschild(order, map_parameter): u, x, dx, d2x = mapped_derivatives(order, map_parameter) ell = 2 zero = np.zeros_like(x) c2 = [x * (x - 1) ** 2 / 4, zero, zero] c1 = [ (x - 1) * (3 * x - 1) / 4, 1j * (-2 * x**2 + 4 * x - 1), zero, ] c0 = [ -(ell**2 + ell + 3 * x - 3) / 4, -2j * (x - 1), -4 * (x - 2), ] return u, x, *quadratic_pencil(c0, c1, c2, dx, d2x) def kerr_coefficients(x): a = 0.4 angular_a = 3.767574385986149 + 0.053240584492497j d = 2 * np.sqrt(1 - a**2) zero = np.zeros_like(x) p = x * (1 - x) ** 2 b0 = (1 - x) * (3 - x + 2j * a * (1 - x) / d) b1 = 2j * ( d * x - 3 * x**2 + 4 * x - 1 - 2 * (1 - x) ** 2 / d ) c0 = (1 - x) * (1 - 2j * a / d) - angular_a c1 = ( -2 * a + 8 * a * (x - 1) / d + 1j * ((d + 2) ** 2 + 2 * (d - 2) * x) / d ) c2 = (d + 2) * (d**2 + 14 * d + 32 * (1 - x)) / (4 * d) return [c0, c1, c2], [b0, b1, zero], [p, zero, zero] def kerr_radial(order, map_parameter): u, x, dx, d2x = mapped_derivatives(order, map_parameter) c0, c1, c2 = kerr_coefficients(x) return u, x, *quadratic_pencil(c0, c1, c2, dx, d2x) BUILDERS = { "poschl-teller": poschl_teller, "schwarzschild": schwarzschild, "kerr-radial": kerr_radial, } class MappedPolynomial: """The degree-N collocation polynomial in u, NOT a new polynomial in x. Fixed Lobatto weights remove the interpolator's randomized weight construction. Off-grid and endpoint derivatives always use the inverse map and both terms of the second-derivative chain rule. """ def __init__(self, u, values, map_parameter): self.u = np.asarray(u) self.values = np.asarray(values) self.map_parameter = map_parameter weights = (-1.0)**np.arange(len(u)) weights[[0, -1]] *= 0.5 self.interpolant = BarycentricInterpolator(self.u, self.values, wi=weights) def evaluate(self, targets, derivative=0): x = np.asarray(targets) c = self.map_parameter denominator = 1-(1-c)*x u = c*x/denominator if derivative == 0: return self.interpolant(u) ux = c/denominator**2 first = self.interpolant.derivative(u, der=1) if derivative == 1: return first*ux if derivative == 2: uxx = 2*c*(1-c)/denominator**3 return self.interpolant.derivative(u, der=2)*ux**2+first*uxx raise ValueError("Only derivatives 0,1,2 are used by this finite benchmark") def coefficients(self): coefficients = dct(self.values, type=1)/(len(self.u)-1) coefficients[[0, -1]] *= 0.5 return coefficients def continuous_coefficients(case, x): x = np.asarray(x) zero = np.zeros_like(x) if case == "poschl-teller": alpha = 0.2 return ( [-1 + zero, 1j * alpha + zero, 1 + zero], [alpha**2 * (1 - 2 * x), 1j * alpha * (2 * x - 1), zero], [alpha**2 * x * (1 - x), zero, zero], ) if case == "schwarzschild": ell = 2 return ( [-(ell**2 + ell + 3 * x - 3) / 4, -2j * (x - 1), -4 * (x - 2)], [(x - 1) * (3 * x - 1) / 4, 1j * (-2 * x**2 + 4 * x - 1), zero], [x * (x - 1) ** 2 / 4, zero, zero], ) return kerr_coefficients(x) def polynomial_residual(omega, vector, matrices): residual = (matrices[0] + omega * matrices[1] + omega**2 * matrices[2]) @ vector scale = ( np.linalg.norm(matrices[0]) + abs(omega) * np.linalg.norm(matrices[1]) + abs(omega) ** 2 * np.linalg.norm(matrices[2]) ) * np.linalg.norm(vector) return float(np.linalg.norm(residual) / scale) def sampled_residual(terms): """Global sampled backward residual, defined through field nodes. Unlike max(pointwise ratios), this does not divide by a zero local denominator at an exact eigenfunction node. A zero or nonfinite global denominator is undefined and fails the gate; it is not regularized away. """ scale = float(np.max(sum(np.abs(term) for term in terms))) if not np.isfinite(scale) or scale<=0: return float("inf") return float(np.max(np.abs(sum(terms)))/scale) def transformed_residual(case, omega, polynomial): points = np.linspace(0.05, 0.95, 37) y = polynomial.evaluate(points) yp = polynomial.evaluate(points, 1) ypp = polynomial.evaluate(points, 2) c0, c1, c2 = continuous_coefficients(case, points) terms = [c2[q] * ypp + c1[q] * yp + c0[q] * y for q in range(3)] return sampled_residual([terms[0],omega*terms[1],omega**2*terms[2]]) def unfactored_terms(case, omega, polynomial): t = np.linspace(0.05, 0.95, 37) y = polynomial.evaluate(t) yp = polynomial.evaluate(t, 1) ypp = polynomial.evaluate(t, 2) if case == "poschl-teller": alpha, p = 0.2, -1j * omega / 0.4 g = p * (1 / t - 1 / (1 - t)) gp = p * (-1 / t**2 - 1 / (1 - t) ** 2) psi_x = yp + g * y psi_xx = ypp + 2 * g * yp + (g**2 + gp) * y xdot = 2 * alpha * t * (1 - t) xddot = 4 * alpha**2 * t * (1 - t) * (1 - 2 * t) terms = [ xdot**2 * psi_xx + xddot * psi_x, (omega**2 - 4 * t * (1 - t)) * y, ] elif case == "schwarzschild": r = 2 / (1 - t) f = t dx_dr = 2 / r**2 d2x_dr2 = -4 / r**3 g = 1j * omega - 2j * omega / (r - 2) + 4j * omega / r gp = 2j * omega / (r - 2) ** 2 - 4j * omega / r**2 y_r = yp * dx_dr y_rr = ypp * dx_dr**2 + yp * d2x_dr2 x_r = y_r + g * y x_rr = y_rr + 2 * g * y_r + (g**2 + gp) * y potential = f * (6 / r**2 - 6 / r**3) terms = [f**2 * x_rr, f * dx_dr * x_r, (omega**2 - potential) * y] else: a, m, s = 0.4, 1, -2 angular_a = 3.767574385986149 + 0.053240584492497j d = 2 * np.sqrt(1 - a**2) rp, rm = 1 + d / 2, 1 - d / 2 r = (rp - t * rm) / (1 - t) dx_dr = (1 - t) ** 2 / d d2x_dr2 = -2 * (1 - t) ** 3 / d**2 sigma = (2 * rp * omega - m * a) / d p = -s - 1j * sigma q = -1 - s + 2j * omega + 1j * sigma g = 1j * omega + p / (r - rp) + q / (r - rm) gp = -p / (r - rp) ** 2 - q / (r - rm) ** 2 y_r = yp * dx_dr y_rr = ypp * dx_dr**2 + yp * d2x_dr2 radial_r = y_r + g * y radial_rr = y_rr + 2 * g * y_r + (g**2 + gp) * y delta = (r - rp) * (r - rm) k = (r**2 + a**2) * omega - m * a potential = ( (k**2 + 4j * (r - 1) * k) / delta - 8j * omega * r - angular_a - a**2 * omega**2 + 2 * a * m * omega ) terms = [ delta * radial_rr, (s + 1) * 2 * (r - 1) * radial_r, potential * y, ] return terms def unfactored_residual(case, omega, polynomial): return sampled_residual(unfactored_terms(case,omega,polynomial)) def endpoint_row_residuals(case, omega, polynomial): points = np.array([0.0, 1.0]) y = polynomial.evaluate(points) yp = polynomial.evaluate(points, 1) c0, c1, _ = continuous_coefficients(case, points) derivative_term = (c1[0] + omega * c1[1] + omega**2 * c1[2]) * yp value_term = (c0[0] + omega * c0[1] + omega**2 * c0[2]) * y scale = (1 + abs(omega) + abs(omega) ** 2) * (np.abs(y) + np.abs(yp)) + 1e-300 residual = np.abs(derivative_term + value_term) / scale return float(residual[0]), float(residual[1]) def phase_residuals(case, omega, polynomial): points = np.array([1e-7, 1 - 1e-7]) y = polynomial.evaluate(points) yp = polynomial.evaluate(points, 1) if np.min(np.abs(y)) <= 1e-12*np.max(np.abs(y)): return float("inf"), float("inf") if case == "poschl-teller": alpha, p = 0.2, -1j * omega / 0.4 g = p * (1 / points - 1 / (1 - points)) dlog = 2 * alpha * points * (1 - points) * (g + yp / y) return float(abs(dlog[0] + 1j * omega)), float(abs(dlog[1] - 1j * omega)) if case == "schwarzschild": r = 2 / (1 - points) dx_dr = 2 / r**2 g = 1j * omega - 2j * omega / (r - 2) + 4j * omega / r dlog = points * (g + (yp / y) * dx_dr) return float(abs(dlog[0] + 1j * omega)), float(abs(dlog[1] - 1j * omega)) a, m, s = 0.4, 1, -2 d = 2 * np.sqrt(1 - a**2) rp, rm = 1 + d / 2, 1 - d / 2 r = (rp - points * rm) / (1 - points) dx_dr = (1 - points) ** 2 / d sigma = (2 * rp * omega - m * a) / d p = -s - 1j * sigma q = -1 - s + 2j * omega + 1j * sigma g = 1j * omega + p / (r - rp) + q / (r - rm) delta = (r - rp) * (r - rm) dr_drstar = delta / (r**2 + a**2) dlog_r = dr_drstar * (g + (yp / y) * dx_dr) omega_h = a / (rp**2 + a**2) horizon = dlog_r[0] + dr_drstar[0] * s * 2 * (r[0] - 1) / delta[0] infinity = dlog_r[1] + dr_drstar[1] * (1 + 2 * s) / r[1] return ( float(abs(horizon + 1j * (omega - m * omega_h))), float(abs(infinity - 1j * omega)), ) def solve(case, order, map_parameter): u, x, values, vectors, matrices = BUILDERS[case](order, map_parameter) rows = [] for column, omega in enumerate(values): if not np.isfinite(omega) or not (0 < omega.real < 1 and -1 < omega.imag < 0): continue vector = vectors[:, column] vector /= np.linalg.norm(vector) polynomial = MappedPolynomial(u, vector, map_parameter) horizon_row, infinity_row = endpoint_row_residuals(case, omega, polynomial) ingoing_phase, outgoing_phase = phase_residuals(case, omega, polynomial) coefficients = polynomial.coefficients() rows.append({ "case": case, "N": order, "map": map_parameter, "omega_real": omega.real, "omega_imag": omega.imag, "pencil_residual": polynomial_residual(omega, vector, matrices), "transformed_residual": transformed_residual(case, omega, polynomial), "unfactored_residual": unfactored_residual(case, omega, polynomial), "horizon_row_residual": horizon_row, "infinity_row_residual": infinity_row, "ingoing_phase_residual": ingoing_phase, "outgoing_phase_residual": outgoing_phase, "coefficient_tail": ( np.max(np.abs(coefficients[-4:])) / (np.max(np.abs(coefficients)) + 1e-300) ), "reference_difference": abs(omega - REFERENCES[case]), "_vector": vector, "_x": x, "_u": u, "_polynomial": polynomial, "_samples": polynomial.evaluate(COMPARE_POINTS), }) rows = sorted(rows, key=lambda row: (row["omega_real"], row["omega_imag"])) for index, row in enumerate(rows): row["candidate_id"] = "%s:N%d:c%.1f:j%d" % (case, order, map_parameter, index) row["matching_checks"] = [] return rows def frequency(row): return row["omega_real"] + 1j * row["omega_imag"] def eigenfunction_mismatch(first, second): left, right = first["_samples"], second["_samples"] norm_left = float(np.vdot(left, COMPARE_WEIGHTS*left).real) norm_right = float(np.vdot(right, COMPARE_WEIGHTS*right).real) if not np.isfinite(norm_left) or not np.isfinite(norm_right) or norm_left<=0 or norm_right<=0: return float("inf") overlap = abs(np.vdot(left, COMPARE_WEIGHTS*right))/np.sqrt(norm_left*norm_right) if not np.isfinite(overlap): return float("inf") return float(max(0, 1 - min(1, overlap))) def nearest_diagnostics(left_frequencies, right_frequencies): """Reference-free two-way frequency ranking and explicit ambiguity gap.""" if len(left_frequencies) == 0: return [] if len(right_frequencies) == 0: return [{"nearest_index": None, "best_distance": None, "second_distance": None, "second_best_gap": None, "reciprocal": False, "unambiguous": False} for _ in left_frequencies] distances = np.abs(np.asarray(left_frequencies)[:, None]-np.asarray(right_frequencies)[None, :]) reverse = np.argmin(distances, axis=0) result = [] for i, values in enumerate(distances): ranked = np.argsort(values, kind="stable") j = int(ranked[0]) best = float(values[j]) second = float(values[ranked[1]]) if len(ranked)>1 else None gap = second-best if second is not None else None result.append({"nearest_index": j, "best_distance": best, "second_distance": second, "second_best_gap": gap, "reciprocal": int(reverse[j])==i, "unambiguous": gap is None or gap>THRESHOLDS["second_best_gap"]}) return result def configuration(row): return (row["case"], row["N"], row["map"]) def compare_configurations(left, right, edge_kind, peer_config_left, peer_config_right): forward = nearest_diagnostics([frequency(r) for r in left], [frequency(r) for r in right]) backward = nearest_diagnostics([frequency(r) for r in right], [frequency(r) for r in left]) edges = [] for rows, peers, selected, reversed_selected, peer_config in ( (left, right, forward, backward, peer_config_right), (right, left, backward, forward, peer_config_left), ): for row, diagnostic in zip(rows, selected): record = {"edge_kind": edge_kind, "peer_configuration": list(peer_config), **diagnostic} index = diagnostic["nearest_index"] if index is None: record.update({"peer_candidate_id": None, "reverse_nearest_candidate_id": None, "reverse_best_distance": None, "reverse_second_distance": None, "reverse_second_best_gap": None, "reverse_unambiguous": False, "eigenfunction_mismatch": None, "passed": False, "failure_reasons": ["empty_peer_configuration"]}) else: peer = peers[index] reverse = reversed_selected[index] reverse_id = rows[reverse["nearest_index"]]["candidate_id"] mismatch = eigenfunction_mismatch(row, peer) reasons = [] if not diagnostic["reciprocal"]: reasons.append("nonreciprocal") if not diagnostic["unambiguous"]: reasons.append("forward_ambiguous") if not reverse["unambiguous"]: reasons.append("reverse_ambiguous") if diagnostic["best_distance"]>=THRESHOLDS["frequency_difference"]: reasons.append("frequency_difference") if not np.isfinite(mismatch) or mismatch>=THRESHOLDS["eigenfunction_mismatch"]: reasons.append("eigenfunction_mismatch") record.update({"peer_candidate_id": peer["candidate_id"], "reverse_nearest_candidate_id": reverse_id, "reverse_best_distance": reverse["best_distance"], "reverse_second_distance": reverse["second_distance"], "reverse_second_best_gap": reverse["second_best_gap"], "reverse_unambiguous": reverse["unambiguous"], "eigenfunction_mismatch": mismatch, "passed": not reasons, "failure_reasons": reasons}) if not reasons: edges.append(tuple(sorted((row["candidate_id"], peer["candidate_id"])))) row["matching_checks"].append(record) return set(edges) def components_and_acceptance(rows, cases, orders, maps): by_config = {(case, n, c): [] for case in cases for n in orders for c in maps} by_id = {row["candidate_id"]: row for row in rows} for row in rows: by_config[configuration(row)].append(row) edges = set() for case in cases: for c in maps: for n1, n2 in zip(orders[:-1], orders[1:]): a, b = (case, n1, c), (case, n2, c) edges |= compare_configurations(by_config[a], by_config[b], "resolution", a, b) for n in orders: for i, c1 in enumerate(maps): for c2 in maps[i+1:]: a, b = (case, n, c1), (case, n, c2) edges |= compare_configurations(by_config[a], by_config[b], "map", a, b) graph = {key:set() for key in by_id} for a, b in edges: graph[a].add(b) graph[b].add(a) components = [] unseen = set(by_id) while unseen: stack = [min(unseen)] visited = set() while stack: key = stack.pop() if key not in visited: visited.add(key) stack.extend(sorted(graph[key]-visited, reverse=True)) unseen -= visited members = sorted(visited) configs = [configuration(by_id[key]) for key in members] unique = len(configs)==len(set(configs)) branch_id = "branch-%04d" % len(components) components.append({"branch_id":branch_id, "members":members, "branch_unique":unique, "reference_label":None}) for key in members: by_id[key]["branch_id"] = branch_id by_id[key]["branch_unique"] = unique # Labels are assigned after reference-free matching. A single highest-N, # identity-map anchor identifies each finite reference branch, never an # independent nearest-to-reference selection in each candidate bucket. anchors = [] for case in cases: group = by_config[(case, orders[-1], 1.0)] ranked = sorted(group, key=lambda row: abs(frequency(row)-REFERENCES[case])) best = abs(frequency(ranked[0])-REFERENCES[case]) if ranked else None second = abs(frequency(ranked[1])-REFERENCES[case]) if len(ranked)>1 else None gap = second-best if second is not None else None labeled = bool(ranked and bestTHRESHOLDS["second_best_gap"]) and ranked[0]["branch_unique"]) branch_id = ranked[0]["branch_id"] if labeled else None label = case+":fundamental-reference" if labeled else None anchors.append({"case":case,"configuration":[case,orders[-1],1.0], "candidate_id":ranked[0]["candidate_id"] if ranked else None, "reference_difference":best,"second_reference_difference":second, "reference_gap":gap,"labeled":labeled,"branch_id":branch_id}) if labeled: next(component for component in components if component["branch_id"]==branch_id)["reference_label"] = label labels = {component["branch_id"]:component["reference_label"] for component in components} for row in rows: local_fields = [key for key in THRESHOLDS if key not in ("frequency_difference","second_best_gap","eigenfunction_mismatch","reference_difference")] local_failures = [key for key in local_fields if not np.isfinite(row[key]) or row[key]>=THRESHOLDS[key]] row["local_pass"] = not local_failures row["matching_pass"] = bool(row["matching_checks"]) and all(check["passed"] for check in row["matching_checks"]) row["reference_label"] = labels[row["branch_id"]] row["reference_pass"] = bool(row["reference_label"] and row["reference_difference"]minimum, (name,value,minimum) negatives[name] = {"residual":value,"minimum":minimum,"rejected":True} for c in MAPS: u, x, dx, d2x = mapped_derivatives(12,c) nodal = 1+u+u**3+u**6 polynomial = MappedPolynomial(u,nodal,c) points = np.r_[np.linspace(0.03,0.97,37),0.0,1.0] denominator = 1-(1-c)*points target_u = c*points/denominator ux, uxx = c/denominator**2,2*c*(1-c)/denominator**3 exact = [1+target_u+target_u**3+target_u**6, (1+3*target_u**2+6*target_u**5)*ux, (6*target_u+30*target_u**4)*ux**2+(1+3*target_u**2+6*target_u**5)*uxx] for order in (0,1,2): error = np.abs(polynomial.evaluate(points,order)-exact[order]) below("mapped_polynomial_c%.1f_offgrid_d%d"%(c,order),np.max(error[:-2]),1e-8) below("mapped_polynomial_c%.1f_endpoint_d%d"%(c,order),np.max(error[-2:]),1e-8) below("mapped_polynomial_c%.1f_nodal_first"%c,np.max(np.abs(dx@nodal-polynomial.evaluate(x,1))),1e-8) below("mapped_polynomial_c%.1f_nodal_second"%c,np.max(np.abs(d2x@nodal-polynomial.evaluate(x,2))),1e-8) coefficients = polynomial.coefficients() below("mapped_polynomial_c%.1f_DCT_reconstruction"%c,np.max(np.abs(np.polynomial.chebyshev.chebval(2*target_u-1,coefficients)-exact[0])),1e-12) phase_x = np.array([1e-7,1-1e-7]) phase_den = 1-(1-c)*phase_x phase_u = c*phase_x/phase_den exact_log = (1+3*phase_u**2+6*phase_u**5)*c/phase_den**2/(1+phase_u+phase_u**3+phase_u**6) below("mapped_polynomial_c%.1f_log_derivative"%c,np.max(np.abs(polynomial.evaluate(phase_x,1)/polynomial.evaluate(phase_x)-exact_log)),1e-9) if c!=1: differences = x[:,None]-x[None,:] np.fill_diagonal(differences,1) weights = 1/np.prod(differences,axis=1) old = BarycentricInterpolator(x,nodal,wi=weights/np.max(np.abs(weights))) above("wrong_physical_x_interpolation_c%.1f"%c,np.max(np.abs(old.derivative(points,der=2)-exact[2])),1e-8) above("omit_inverse_second_chain_term_c%.1f"%c,np.max(np.abs((1+3*target_u**2+6*target_u**5)*uxx)),0.1) alpha,V0 = 0.2,1.0 omega = np.sqrt(V0-alpha*alpha/4)-1.5j*alpha x = np.linspace(0.05,0.95,37) y,yp,ypp = 1-2*x,-2+0*x,0*x c0,c1,c2 = continuous_coefficients("poschl-teller",x) transformed = sum(omega**q*(c2[q]*ypp+c1[q]*yp+c0[q]*y) for q in range(3)) below("PT_n1_nonconstant_transformed",np.max(np.abs(transformed)),1e-12) p = -1j*omega/(2*alpha) factor = np.exp(p*(np.log(x)+np.log(1-x))) g,gp = p*(1/x-1/(1-x)),p*(-1/x**2-1/(1-x)**2) derivative = factor*(yp+g*y) second = factor*(ypp+2*g*yp+(g*g+gp)*y) xdot = 2*alpha*x*(1-x) xddot = 4*alpha**2*x*(1-x)*(1-2*x) raw = [xdot*xdot*second,xddot*derivative,(omega*omega-4*V0*x*(1-x))*factor*y] below("PT_n1_nonconstant_full_raw_ODE",sampled_residual(raw),1e-12) above("old_PT_first_derivative_sign_n1",np.max(np.abs(-2*alpha*alpha*(1-2*x)*yp)),0.1) for c in MAPS: u,x,values,vectors,matrices = poschl_teller(12,c) exact_vector = 1-2*x polynomial = MappedPolynomial(u,exact_vector,c) below("PT_n1_exact_grid_pencil_c%.1f"%c,polynomial_residual(omega,exact_vector,matrices),1e-8) below("PT_n1_exact_grid_raw_c%.1f"%c,unfactored_residual("poschl-teller",omega,polynomial),1e-4) if c==0.7: raw_terms = unfactored_terms("poschl-teller",omega,polynomial) old_pointwise = np.max(np.abs(sum(raw_terms))/(sum(np.abs(term) for term in raw_terms)+1e-300)) above("exact_PT_node_pointwise_relative_residual_instability",old_pointwise,1e-2) undefined = sampled_residual([np.zeros(3),np.zeros(3)]) assert not np.isfinite(undefined) negatives["zero_global_denominator_is_undefined"] = {"rejected":True,"normalization_defined":False,"normalized_residual":None} # The common physical-x weighted overlap must not depend on nodal or # overall complex normalization. Degree-one y(x) is represented exactly # at c=1; the test comparison is with its independent closed form. u,x,_,_ = mapped_derivatives(12,1.0) first = MappedPolynomial(u,1-2*x,1.0) second = MappedPolynomial(u,(2+3j)*(1-2*x),1.0) overlap_error = eigenfunction_mismatch({"_samples":first.evaluate(COMPARE_POINTS)}, {"_samples":second.evaluate(COMPARE_POINTS)}) below("physical_grid_overlap_complex_rescaling",overlap_error,1e-14) analytical = 1-2*COMPARE_POINTS below("physical_grid_nonconstant_eigenfunction_samples",np.max(np.abs(first.evaluate(COMPARE_POINTS)-analytical)),1e-12) undefined_overlap = eigenfunction_mismatch({"_samples":np.zeros_like(COMPARE_POINTS)}, {"_samples":analytical}) assert not np.isfinite(undefined_overlap) negatives["zero_eigenfunction_overlap_is_undefined"] = {"rejected":True,"normalization_defined":False,"eigenfunction_mismatch":None} ambiguous = nearest_diagnostics([0.4-0.1j],[0.400001-0.1j,0.400002-0.1j])[0] assert not ambiguous["unambiguous"] negatives["nearest_reference_does_not_resolve_ambiguity"] = {"rejected":True,"diagnostic":ambiguous} nonmutual = nearest_diagnostics([0.4-0.1j,0.400001-0.1j],[0.4000001-0.1j])[1] assert not nonmutual["reciprocal"] negatives["one_way_nearest_neighbor_is_not_mutual"] = {"rejected":True,"diagnostic":nonmutual} return {"residuals":residuals,"negative_controls":negatives} def clean(value): if isinstance(value,dict): return {key:clean(item) for key,item in value.items() if not key.startswith("_")} if isinstance(value,(list,tuple)): return [clean(item) for item in value] if isinstance(value,np.generic): return clean(value.item()) if isinstance(value,float) and not np.isfinite(value): return None return value def main(): parser = argparse.ArgumentParser(description="Finite corrected spectral reference benchmark; no network or package installation.") parser.add_argument("--cases",default=",".join(BUILDERS)) parser.add_argument("--orders",default=",".join(map(str,ORDERS))) parser.add_argument("--maps",default=",".join(map(str,MAPS))) parser.add_argument("--write-candidates",type=Path,default=Path("candidates.csv")) parser.add_argument("--write-residuals",type=Path,default=Path("residuals.csv")) parser.add_argument("--write-checks",type=Path,default=Path("spectral-qnm-checks.json")) args = parser.parse_args() cases,orders,maps = args.cases.split(","),tuple(map(int,args.orders.split(","))),tuple(map(float,args.maps.split(","))) if cases!=list(BUILDERS) or orders!=ORDERS or maps!=MAPS: parser.error("This reproducible benchmark supports exactly the declared three cases, eight orders and three maps.") requirements = Path(__file__).with_name("requirements.txt") requirements_bytes = requirements.read_bytes() controls = numerical_controls() print("Mapped-polynomial, nonconstant PT and matching controls passed.",flush=True) rows = [] configurations = [] for case in cases: for n in orders: for c in maps: selected = solve(case,n,c) rows.extend(selected) configurations.append({"case":case,"N":n,"map":c,"candidate_count":len(selected)}) print("Completed all24 configurations:",case,flush=True) branches,anchors,edges = components_and_acceptance(rows,cases,orders,maps) accepted = [row for row in rows if row["accepted"]] fields = ["case","N","map","candidate_id","omega_real","omega_imag", "pencil_residual","transformed_residual","unfactored_residual", "horizon_row_residual","infinity_row_residual","ingoing_phase_residual","outgoing_phase_residual", "coefficient_tail","reference_difference","max_resolution_difference","resolution_eigenfunction_mismatch", "max_map_difference","map_eigenfunction_mismatch","local_pass","matching_pass","branch_id","branch_unique", "reference_label","reference_pass","accepted","failure_reasons"] for output,selected in ((args.write_candidates,rows),(args.write_residuals,accepted)): with output.open("w",newline="",encoding="utf-8") as stream: writer = csv.DictWriter(stream,fieldnames=fields) writer.writeheader() for row in selected: public = clean(row) public["failure_reasons"] = "|".join(public["failure_reasons"]) writer.writerow({key:public[key] for key in fields}) report = {"schema":"spectral-qnm-corrected-v1","status":"computed finite benchmark; independent acceptance recorded separately", "versions":{"python":platform.python_version(),"numpy":np.__version__,"scipy":scipy.__version__}, "source_sha256":hashlib.sha256(Path(__file__).read_bytes()).hexdigest(), "requirements_sha256":hashlib.sha256(requirements_bytes).hexdigest(), "inputs":{"cases":cases,"orders":orders,"maps":maps,"frequency_window":{"real":[0,1],"imag":[-1,0],"strict":True}}, "thresholds":THRESHOLDS, "interpolation":{"space":"degree N in computational u","map":"x=u/[c+(1-c)u]","inverse":"u=cx/[1-(1-c)x]","second_derivative":"y_xx=y_uu*u_x^2+y_u*u_xx","lobatto_weights":"(-1)^j, endpoints halved","offgrid_points":{"start":0.05,"end":0.95,"count":37},"ODE_residual":"max_i abs(sum_a term_ai)/max_i sum_a abs(term_ai); zero or nonfinite global denominator fails","phase_points":[1e-7,1-1e-7],"phase_zero_relative_floor":1e-12,"eigenfunction_quadrature":"64-point Gauss-Legendre in physical x on[0,1]","eigenfunction_mismatch":"1-abs(sum(w conj(y1) y2))/sqrt(sum(w abs(y1)^2) sum(w abs(y2)^2))"}, "matching_contract":{"resolution_peers":"Every available adjacent N in the declared ordered list at fixed map; lowest/highest N have one neighbor, interior N two.","map_peers":"Both other maps at the same N.","ranking":"Mutual nearest complex frequency, independent of references.","ambiguity":"Second-best minus best distance must exceed5e-6 in both directions; absent second candidate has null gap and is unambiguous.","branch_uniqueness":"Each connected component of passed pair matches has at most one candidate per(case,N,map).","reference_label":"After matching only: highest N, map1 anchor within5e-6 and with reference gap>5e-6 labels the component.","acceptance":"local_pass and matching_pass and branch_unique and reference_pass","undefined_diagnostic":"JSON null and empty CSV field; its gate fails, never a finite success sentinel."}, "scope":{"Schwarzschild":"GM=1 axial s=-2 ell2 fundamental reference","Kerr":"GM=1 s=-2 ell2 m1 a=.4, fixed complex A=3.767574385986149+.053240584492497i; not a coupled Kerr solver","PT":"V0=1 alpha=.2 fundamental reference; separate exact nonconstant n1 controls","limitations":["Finite double-precision reference benchmark, not a general physical-mode classifier.","Smooth endpoint rows and finite samples do not prove a retarded analytic domain or continuum convergence.","No reference participates in numerical pair matching.","All candidate failures are retained; increased N need not improve precision."]}, "controls":controls,"configurations":configurations,"candidates":rows,"branches":branches,"reference_anchors":anchors,"matching_edges":edges, "accepted_candidate_ids":[row["candidate_id"] for row in accepted], "counts":{"configurations":len(configurations),"candidates":len(rows),"accepted":len(accepted),"branches":len(branches),"matching_edges":len(edges),"positive_controls":len(controls["residuals"]),"negative_controls":len(controls["negative_controls"])} } args.write_checks.write_text(json.dumps(clean(report),indent=2,sort_keys=True,allow_nan=False)+"\n") print(json.dumps(report["counts"],sort_keys=True)) if __name__ == "__main__": main()