|
| 1 | +"""The Flight Recorder. |
| 2 | +
|
| 3 | +A showpiece: film one brutally hard homotopy path -- the descent to a highly *singular* |
| 4 | +solution -- and read out the adaptive-precision tracker's full telemetry, step by step. |
| 5 | +
|
| 6 | +The subject is the origin (0, 0), where two rotated rose curves r = sin(m*theta) and |
| 7 | +r = sin(n*theta) meet. For (m, n) = (7, 5) that intersection has multiplicity 35: thirty-five |
| 8 | +homotopy paths pile into the same point, and the ones that get there have to fight for every |
| 9 | +digit. We ask the solver for tight final accuracy, then watch a single path's instruments as it |
| 10 | +goes: |
| 11 | +
|
| 12 | + * the endgame **spiral** -- the Cauchy endgame samples a circle around the singular endpoint; |
| 13 | + with cycle number c the solution winds c times as t -> 0 (here c = 7), a log-radial spiral; |
| 14 | + * **precision** climbing 16 -> 20 -> 30 -> ... digits as adaptive precision escalates into |
| 15 | + mpfr to keep the accuracy the tight tolerance demands; |
| 16 | + * the **condition number** blowing up by many orders of magnitude as the Jacobian degenerates; |
| 17 | + * the **step size** sawtoothing -- grown when the going is easy, cut hard (a rejected step) |
| 18 | + when it is not. |
| 19 | +
|
| 20 | +All of it is real, captured by a ``PathDataCollector`` on the path's own tracker (via a |
| 21 | +``SolutionPathCollector`` over the whole solve), then laid out as a cockpit. The point of the |
| 22 | +family is that it is *crankable*: raise (m, n) or tighten the tolerance and the singular point -- |
| 23 | +and the tracker's struggle -- gets arbitrarily worse. |
| 24 | +
|
| 25 | +Regenerated through ``tools/refresh_doc_artifacts.py``. Raster (PNG) showpiece, not a doctest. |
| 26 | +
|
| 27 | +Run standalone: python flight_recorder.py |
| 28 | +""" |
| 29 | + |
| 30 | +import math |
| 31 | +import os |
| 32 | + |
| 33 | +import numpy as np |
| 34 | + |
| 35 | +import matplotlib |
| 36 | +matplotlib.use('Agg') |
| 37 | +import matplotlib.pyplot as plt |
| 38 | +import matplotlib.colors as mcolors |
| 39 | +from matplotlib.collections import LineCollection |
| 40 | + |
| 41 | +import bertini |
| 42 | +from bertini import ZeroDimSolver, SolutionPathCollector |
| 43 | +from bertini.sympy_bridge import from_sympy |
| 44 | + |
| 45 | +_OUT = os.path.dirname(os.path.abspath(__file__)) |
| 46 | +_BG = '#070a10' |
| 47 | +_INK = '#d6e2f0' |
| 48 | +_DIM = '#7f8da0' |
| 49 | + |
| 50 | + |
| 51 | +# --- the crankable singular system -------------------------------------------------------------- |
| 52 | + |
| 53 | +def _rose(k): |
| 54 | + """The rectangular equation of the rose r = sin(k*theta): (x^2+y^2)^((k+1)/2) = Im[(x+iy)^k], |
| 55 | + a polynomial for odd k. Returned as a bertini function-tree node via the sympy bridge.""" |
| 56 | + from sympy import symbols, im, I |
| 57 | + xs, ys = symbols('x y', real=True) |
| 58 | + return from_sympy((xs**2 + ys**2)**((k + 1) // 2) - im((xs + I * ys)**k)) |
| 59 | + |
| 60 | +def system_rhodonea(m, n): |
| 61 | + """Two rose curves, the second rotated by a random angle so their only structured coincidence |
| 62 | + is the highly singular meeting at the origin. (m, n) = (7, 5) -> multiplicity 35 at (0,0). |
| 63 | + Returns (system, rotation_angle) so the geometry can be drawn to match the solved system.""" |
| 64 | + x, y = bertini.variables(list('xy')) |
| 65 | + f1, f2 = _rose(m), _rose(n) |
| 66 | + t = bertini.random_real() # random rotation (seed fixed by the caller) |
| 67 | + angle = float(t.real) |
| 68 | + rx = bertini.cos(t) * x + bertini.sin(t) * y |
| 69 | + ry = -bertini.sin(t) * x + bertini.cos(t) * y |
| 70 | + f2 = f2.subs({x: rx, y: ry}) |
| 71 | + sys = bertini.System() |
| 72 | + sys.add_variable_group(x, y) |
| 73 | + sys.add([f1, f2]) |
| 74 | + return sys, angle |
| 75 | + |
| 76 | +def rose_curve(k, rotation=0.0, n=1400): |
| 77 | + """Sample the real rose r = sin(k*theta) as (X, Y), optionally rotated to match the system.""" |
| 78 | + th = np.linspace(0, 2 * math.pi, n) |
| 79 | + r = np.sin(k * th) |
| 80 | + X, Y = r * np.cos(th), r * np.sin(th) |
| 81 | + c, s = math.cos(rotation), math.sin(rotation) |
| 82 | + return c * X - s * Y, s * X + c * Y |
| 83 | + |
| 84 | +# --- solving and recording ---------------------------------------------------------------------- |
| 85 | + |
| 86 | +def _solve(m, n, seed, endgame, final_tolerance=None): |
| 87 | + """One serial, deterministic solve of system_rhodonea(m, n) with the chosen endgame, every path |
| 88 | + collected. Returns (collector, meta, rotation).""" |
| 89 | + bertini.recording(False) # WATCH tracking, do not recall it |
| 90 | + bertini.random.set_random_seed(seed) # deterministic system + solve -> stable picture |
| 91 | + system, rotation = system_rhodonea(m, n) |
| 92 | + solver = ZeroDimSolver(system, mptype='adaptive', endgame=endgame) |
| 93 | + if final_tolerance is not None: |
| 94 | + tol = solver.get_config(bertini.nag_algorithm.TolerancesConfig) |
| 95 | + tol.final_tolerance = final_tolerance |
| 96 | + solver.set_config(tol) |
| 97 | + cfg = solver.get_config(bertini.nag_algorithm.ZeroDimConfig) |
| 98 | + cfg.num_threads = 1 # serial -> deterministic path ordering / picture |
| 99 | + solver.set_config(cfg) |
| 100 | + collector = SolutionPathCollector() |
| 101 | + solver.add_observer(collector) |
| 102 | + solver.solve() |
| 103 | + meta = {int(md.path_index): md for md in solver.solution_metadata()} |
| 104 | + return collector, meta, rotation |
| 105 | + |
| 106 | +def record_hard_path(m=7, n=5, final_tolerance=1e-24, seed=2): |
| 107 | + """The per-path cockpit: a tight-tolerance CAUCHY solve (its spiral IS the point), returning the |
| 108 | + richest singular path's telemetry (the one that escalates precision the most) plus its metadata.""" |
| 109 | + collector, meta, _ = _solve(m, n, seed, 'cauchy', final_tolerance) |
| 110 | + P = collector.series[0].DIAGNOSTIC_COLUMNS.index('precision') |
| 111 | + best = None |
| 112 | + for path in collector.series: |
| 113 | + md = meta.get(path.path_index) |
| 114 | + if not (md and md.is_singular): |
| 115 | + continue |
| 116 | + dgn = path.diagnostics() |
| 117 | + key = (int(dgn[:, P].max()), len(dgn)) # most precision, then most steps |
| 118 | + if best is None or key > best[0]: |
| 119 | + best = (key, path, md) |
| 120 | + _, path, md = best |
| 121 | + dgn = path.diagnostics() |
| 122 | + return dict(m=m, n=n, final_tolerance=final_tolerance, path_index=int(path.path_index), |
| 123 | + affine=path.points()[:, 1:] / path.points()[:, 0:1], |
| 124 | + abs_t=dgn[:, 0], condition=dgn[:, 1], precision=dgn[:, 2], stepsize=dgn[:, 3], |
| 125 | + cycle=int(md.cycle_num), multiplicity=int(md.multiplicity), |
| 126 | + precision_digits=int(md.precision_digits), accuracy_digits=int(md.accuracy_digits)) |
| 127 | + |
| 128 | +def record_convergence(m=7, n=5, seed=2): |
| 129 | + """The system-level companion: a POWER-SERIES solve, whose endgame tracks radially toward t=0 |
| 130 | + (no Cauchy loops), so every singular path is an honest single-valued real-time descent onto the |
| 131 | + singular point. Returns the rotation and each path's (abs_t, affine) trajectory.""" |
| 132 | + collector, meta, rotation = _solve(m, n, seed, 'powerseries') |
| 133 | + trajectories, multiplicity = [], 1 |
| 134 | + for path in collector.series: |
| 135 | + md = meta.get(path.path_index) |
| 136 | + if md and md.is_singular: |
| 137 | + trajectories.append((np.abs(path.times()), path.points()[:, 1:] / path.points()[:, 0:1])) |
| 138 | + multiplicity = int(md.multiplicity) |
| 139 | + return dict(m=m, n=n, rotation=rotation, multiplicity=multiplicity, trajectories=trajectories) |
| 140 | + |
| 141 | + |
| 142 | +# --- the cockpit -------------------------------------------------------------------------------- |
| 143 | + |
| 144 | +def _style_axis(ax): |
| 145 | + ax.set_facecolor(_BG) |
| 146 | + for s in ax.spines.values(): |
| 147 | + s.set_color('#26324a') |
| 148 | + ax.tick_params(colors='#8fa0bb', labelsize=8) |
| 149 | + ax.grid(True, color='#141b28', lw=0.7) |
| 150 | + ax.xaxis.label.set_color(_INK); ax.yaxis.label.set_color(_INK) |
| 151 | + |
| 152 | +def render(rec, out): |
| 153 | + plt.rcParams.update({'figure.facecolor': _BG, 'axes.facecolor': _BG, |
| 154 | + 'font.family': 'monospace'}) |
| 155 | + fig = plt.figure(figsize=(13, 8)) |
| 156 | + gs = fig.add_gridspec(3, 2, width_ratios=[1.15, 1.0], height_ratios=[1, 1, 1], |
| 157 | + hspace=0.5, wspace=0.34, left=0.06, right=0.97, top=0.86, bottom=0.09) |
| 158 | + |
| 159 | + n_steps = len(rec['abs_t']) |
| 160 | + step = np.arange(n_steps) |
| 161 | + eg = rec['abs_t'] < 0.1 # the endgame portion (small |t|) |
| 162 | + boundary = int(np.argmax(eg)) if eg.any() else n_steps |
| 163 | + |
| 164 | + # ---- the endgame spiral (left, spanning all rows) ---- |
| 165 | + axS = fig.add_subplot(gs[:, 0]); _style_axis(axS) |
| 166 | + xv = rec['affine'][:, 0] # one coordinate of the solution |
| 167 | + xe, te = xv[eg], rec['abs_t'][eg] |
| 168 | + r = np.log10(np.abs(xe)) - np.log10(np.abs(xe).min()) + 0.15 # log-radial: 0 -> singular pt |
| 169 | + disp = r * np.exp(1j * np.angle(xe)) |
| 170 | + pts = np.column_stack([disp.real, disp.imag]) |
| 171 | + segs = np.stack([pts[:-1], pts[1:]], axis=1) |
| 172 | + lc = LineCollection(segs, cmap='turbo', |
| 173 | + norm=mcolors.LogNorm(max(te.min(), 1e-30), te.max())) |
| 174 | + lc.set_array(0.5 * (te[:-1] + te[1:])); lc.set_linewidth(1.7) |
| 175 | + axS.add_collection(lc) |
| 176 | + axS.scatter([0], [0], marker='*', s=320, color='white', zorder=6) |
| 177 | + axS.scatter([0], [0], marker='*', s=1100, color='#fff0b0', alpha=0.25, zorder=5) |
| 178 | + R = np.abs(disp).max() * 1.1 |
| 179 | + axS.set_xlim(-R, R); axS.set_ylim(-R, R); axS.set_aspect(1.0) |
| 180 | + groups = rec['multiplicity'] // rec['cycle'] if rec['cycle'] else 0 |
| 181 | + axS.set_title(f"Cauchy endgame spiral into the singular point (cycle number c = {rec['cycle']})\n" |
| 182 | + f"this path winds {rec['cycle']}× as t → 0; the mult-{rec['multiplicity']} point is " |
| 183 | + f"{groups} cyclic groups of {rec['cycle']} ({groups}×{rec['cycle']} = {rec['multiplicity']})", |
| 184 | + color=_INK, fontsize=9.5, pad=6) |
| 185 | + cb = fig.colorbar(lc, ax=axS, fraction=0.035, pad=0.015) |
| 186 | + cb.set_label('|t| (log)', color=_INK, labelpad=-2) |
| 187 | + cb.ax.tick_params(colors='#8fa0bb') |
| 188 | + |
| 189 | + def gauge(ax, y, label, color, logy=False, drops=None): |
| 190 | + ax.plot(step, y, color=color, lw=1.4) |
| 191 | + if drops is not None: |
| 192 | + ax.plot(step[drops], y[drops], 'v', color='#ff5a5a', ms=4, alpha=0.8) |
| 193 | + if logy: |
| 194 | + ax.set_yscale('log') |
| 195 | + ax.axvline(boundary, color='#4de0c0', lw=1.0, ls=(0, (3, 3)), alpha=0.8) |
| 196 | + ax.set_ylabel(label); _style_axis(ax) |
| 197 | + ax.set_xlim(0, n_steps - 1) |
| 198 | + |
| 199 | + # ---- precision staircase ---- |
| 200 | + axP = fig.add_subplot(gs[0, 1]) |
| 201 | + gauge(axP, rec['precision'], 'precision\n(digits)', '#ffd24a') |
| 202 | + axP.set_ylim(min(rec['precision']) - 2, max(rec['precision']) + 6) |
| 203 | + axP.text(boundary, axP.get_ylim()[1], ' endgame', color='#4de0c0', fontsize=7, va='top') |
| 204 | + |
| 205 | + # ---- condition number ---- |
| 206 | + axC = fig.add_subplot(gs[1, 1]) |
| 207 | + gauge(axC, np.maximum(rec['condition'], 1.0), 'condition\nnumber', '#ff6fae', logy=True) |
| 208 | + |
| 209 | + # ---- step size (sawtooth), failed steps marked ---- |
| 210 | + axZ = fig.add_subplot(gs[2, 1]) |
| 211 | + drops = np.where(np.diff(rec['stepsize']) < 0)[0] + 1 # a cut step size = a rejected step |
| 212 | + gauge(axZ, np.maximum(rec['stepsize'], 1e-30), 'step size', '#5ad1ff', logy=True, drops=drops) |
| 213 | + axZ.set_xlabel('tracker step (its own clock →)') |
| 214 | + axZ.plot([], [], 'v', color='#ff5a5a', ms=5, label='step cut') |
| 215 | + axZ.legend(loc='lower left', fontsize=7, facecolor=_BG, edgecolor='#26324a', labelcolor=_INK) |
| 216 | + |
| 217 | + # ---- title / readout bar ---- |
| 218 | + fig.text(0.06, 0.955, "✈ FLIGHT RECORDER", color='#4de0c0', fontsize=16, fontweight='bold') |
| 219 | + fig.text(0.06, 0.905, |
| 220 | + f"one path to a multiplicity-{rec['multiplicity']} singular point · " |
| 221 | + f"$x^{{{rec['m']}}}$-rose ∩ $x^{{{rec['n']}}}$-rose at (0,0) · " |
| 222 | + f"final_tolerance = {rec['final_tolerance']:.0e}", |
| 223 | + color=_INK, fontsize=10) |
| 224 | + fig.text(0.97, 0.955, |
| 225 | + f"cycle {rec['cycle']} · precision {int(rec['precision'].min())}→{int(rec['precision'].max())} digits " |
| 226 | + f"· condition ×10^{int(np.log10(rec['condition'].max()))} · {n_steps} steps", |
| 227 | + color='#8fa0bb', fontsize=9, ha='right') |
| 228 | + |
| 229 | + fig.savefig(out, dpi=150, facecolor=_BG) |
| 230 | + plt.close(fig) |
| 231 | + return rec |
| 232 | + |
| 233 | + |
| 234 | +def render_setup(setup, out): |
| 235 | + """The system-level companion (whole solve, not one path): the two rose curves whose crossing |
| 236 | + at (0,0) is the singular target, and all the homotopy paths converging on it loop-free.""" |
| 237 | + plt.rcParams.update({'figure.facecolor': _BG, 'axes.facecolor': _BG, 'font.family': 'monospace'}) |
| 238 | + fig, (axR, axC) = plt.subplots(1, 2, figsize=(13, 6.6)) |
| 239 | + fig.subplots_adjust(left=0.06, right=0.95, top=0.84, bottom=0.09, wspace=0.24) |
| 240 | + |
| 241 | + # --- the geometry: two real rose curves meeting at the origin --- |
| 242 | + _style_axis(axR) |
| 243 | + Xa, Ya = rose_curve(setup['m'], 0.0) |
| 244 | + Xb, Yb = rose_curve(setup['n'], setup['rotation']) |
| 245 | + axR.plot(Xa, Ya, color='#ff6fae', lw=1.6, label=f"$r=\\sin({setup['m']}\\theta)$") |
| 246 | + axR.plot(Xb, Yb, color='#5ad1ff', lw=1.6, label=f"$r=\\sin({setup['n']}\\theta)$, rotated") |
| 247 | + axR.scatter([0], [0], s=200, facecolors='none', edgecolors='white', linewidths=1.3, zorder=6) # ring, not covering |
| 248 | + axR.set_aspect(1.0); axR.set_xlabel('x'); axR.set_ylabel('y') |
| 249 | + axR.legend(loc='upper right', fontsize=9, facecolor=_BG, edgecolor='#26324a', labelcolor=_INK) |
| 250 | + axR.set_title(f"the problem — two rose curves meet at (0,0)\n" |
| 251 | + f"a multiplicity-{setup['multiplicity']} singular intersection", color=_INK, fontsize=10, pad=6) |
| 252 | + |
| 253 | + # --- the solve: every singular path's loop-free descent, converging on (0,0) --- |
| 254 | + _style_axis(axC) |
| 255 | + cmap = plt.get_cmap('turbo') |
| 256 | + trajs = setup['trajectories'] |
| 257 | + npath = len(trajs) |
| 258 | + order = np.argsort([float(np.angle(s[1][0, 0])) for s in trajs]) # hue by start angle |
| 259 | + for rank, i in enumerate(order): |
| 260 | + abs_t, affine = trajs[i] |
| 261 | + xv = affine[:, 0] # one coordinate's complex plane |
| 262 | + pts = np.column_stack([xv.real, xv.imag]) |
| 263 | + segs = np.stack([pts[:-1], pts[1:]], axis=1) |
| 264 | + rgb = cmap((rank + 0.5) / npath) # a distinct hue per path -> 35 followable threads |
| 265 | + axC.add_collection(LineCollection(segs, colors=[rgb], linewidths=1.1, alpha=0.72)) |
| 266 | + axC.plot(xv.real[0], xv.imag[0], 'o', color=rgb, ms=6, mec='white', mew=0.7, zorder=7) # start (t≈1) |
| 267 | + axC.scatter([0], [0], s=170, facecolors='none', edgecolors='white', linewidths=1.4, zorder=8) # target ring |
| 268 | + axC.autoscale(); axC.set_aspect(1.0) |
| 269 | + axC.set_xlabel('Re(x)'); axC.set_ylabel('Im(x)') |
| 270 | + axC.set_title(f"the solve — {npath} homotopy paths converge on it\n" |
| 271 | + "real-time continuation (power-series endgame: no Cauchy loops)", color=_INK, fontsize=10, pad=6) |
| 272 | + |
| 273 | + fig.text(0.06, 0.945, "THE SINGULAR RENDEZVOUS", color='#4de0c0', fontsize=14, fontweight='bold') |
| 274 | + fig.savefig(out, dpi=150, facecolor=_BG) |
| 275 | + plt.close(fig) |
| 276 | + |
| 277 | + |
| 278 | +def main(): |
| 279 | + render_setup(record_convergence(m=7, n=5, seed=2), |
| 280 | + os.path.join(_OUT, 'flight_recorder_setup.png')) |
| 281 | + render(record_hard_path(m=7, n=5, final_tolerance=1e-24, seed=2), |
| 282 | + os.path.join(_OUT, 'flight_recorder.png')) |
| 283 | + |
| 284 | + |
| 285 | +if __name__ == '__main__': |
| 286 | + main() |
0 commit comments