"""Recompute matched circles in a toy flat cubic three-torus. Standard library only. Lengths are in units of the comoving distance to last scattering, chi_rec = 1. Temperatures are in units of the signal RMS. Usage: python3 -I recompute_matched_circles.py tutorial-inputs.json output.json The field is a synthetic Gaussian-like random field made from Gaussian bumps repeated on the lattice L Z^3. It is a teaching model, not a CMB simulation: there is no Sachs-Wolfe, Doppler or ISW physics, no beam and no mask. """ import json import math import random import sys SCRIPT_VERSION = "ac86-recompute-matched-circles-1" def normal_generator(seed): """Box-Muller from random.random(), so results do not depend on gauss().""" uniform = random.Random(seed) spare = [] def draw(): if spare: return spare.pop() radius = math.sqrt(-2.0 * math.log(1.0 - uniform.random())) angle = 2.0 * math.pi * uniform.random() spare.append(radius * math.sin(angle)) return radius * math.cos(angle) return uniform, draw class PeriodicBumpField: """T(x) = sum_j w_j exp(-|x - p_j - n L|^2 / (2 s^2)) summed over lattice images n.""" def __init__(self, box_side, bump_width, bump_count, seed): uniform, normal = normal_generator(seed) self.box_side = box_side self.bump_width = bump_width self.cutoff = 4.0 * bump_width self.cell = self.cutoff bumps = [((uniform.random() * box_side, uniform.random() * box_side, uniform.random() * box_side), normal()) for _ in range(bump_count)] # Copy every bump into the images that can reach the observable sphere # (or the periodicity probe at distance L/2, whichever is farther). reach = max(1.0, 0.5 * box_side + 0.05) + self.cutoff image_range = int(math.ceil(reach / box_side)) + 1 self.grid = {} for (x, y, z), weight in bumps: for i in range(-image_range, image_range + 1): for j in range(-image_range, image_range + 1): for k in range(-image_range, image_range + 1): px, py, pz = x + i * box_side, y + j * box_side, z + k * box_side if max(abs(px), abs(py), abs(pz)) > reach + box_side: continue if px * px + py * py + pz * pz > (reach + 0.01) ** 2: continue key = (math.floor(px / self.cell), math.floor(py / self.cell), math.floor(pz / self.cell)) self.grid.setdefault(key, []).append((px, py, pz, weight)) self.normalisation = 1.0 self.normalisation = self._rms_on_sphere(uniform, 4000) def raw(self, point): x, y, z = point cx, cy, cz = math.floor(x / self.cell), math.floor(y / self.cell), math.floor(z / self.cell) total, inverse, cutoff2 = 0.0, 1.0 / (2.0 * self.bump_width ** 2), self.cutoff ** 2 for i in (cx - 1, cx, cx + 1): for j in (cy - 1, cy, cy + 1): for k in (cz - 1, cz, cz + 1): for px, py, pz, weight in self.grid.get((i, j, k), ()): distance2 = (x - px) ** 2 + (y - py) ** 2 + (z - pz) ** 2 if distance2 > cutoff2: continue total += weight * math.exp(-distance2 * inverse) return total def _rms_on_sphere(self, uniform, count): values = [] for _ in range(count): direction = random_direction(uniform) values.append(self.raw(direction)) mean = sum(values) / count self.mean = mean return math.sqrt(sum((value - mean) ** 2 for value in values) / count) def temperature(self, point): return (self.raw(point) - self.mean) / self.normalisation def random_direction(uniform): z = 2.0 * uniform.random() - 1.0 angle = 2.0 * math.pi * uniform.random() radius = math.sqrt(max(0.0, 1.0 - z * z)) return (radius * math.cos(angle), radius * math.sin(angle), z) def normalise(vector): length = math.sqrt(sum(component * component for component in vector)) return tuple(component / length for component in vector) def perpendicular_frame(axis): helper = (1.0, 0.0, 0.0) if abs(axis[0]) < 0.9 else (0.0, 1.0, 0.0) first = normalise((axis[1] * helper[2] - axis[2] * helper[1], axis[2] * helper[0] - axis[0] * helper[2], axis[0] * helper[1] - axis[1] * helper[0])) second = (axis[1] * first[2] - axis[2] * first[1], axis[2] * first[0] - axis[0] * first[2], axis[0] * first[1] - axis[1] * first[0]) return first, second def circle_points(axis, alpha_degrees, samples, mirrored=False): """Points on the unit last-scattering sphere at angle alpha from axis. The phase phi runs counter-clockwise about +axis; mirrored=True runs it the other way. """ alpha = math.radians(alpha_degrees) first, second = perpendicular_frame(axis) points = [] for index in range(samples): phi = 2.0 * math.pi * index / samples * (-1.0 if mirrored else 1.0) points.append(tuple(math.cos(alpha) * a + math.sin(alpha) * (math.cos(phi) * b + math.sin(phi) * c) for a, b, c in zip(axis, first, second))) return points def circle_statistic(first, second, shift): """CSS 1998 statistic S = 2 / with b shifted by `shift` samples.""" count = len(first) numerator = denominator = 0.0 for index in range(count): a, b = first[index], second[(index + shift) % count] numerator += 2.0 * a * b denominator += a * a + b * b return numerator / denominator def best_over_shifts(first, second): values = [circle_statistic(first, second, shift) for shift in range(len(first))] best = max(range(len(values)), key=values.__getitem__) return values[best], best def lattice_pairs(box_side): """Clone translations g = L(i, j, k) with 0 < |g| < 2; one entry per +-g pair.""" limit = int(math.floor(2.0 / box_side)) pairs = [] for i in range(-limit, limit + 1): for j in range(-limit, limit + 1): for k in range(-limit, limit + 1): if (i, j, k) <= (0, 0, 0) and (-i, -j, -k) != (i, j, k): continue if (i, j, k) == (0, 0, 0): continue length = box_side * math.sqrt(i * i + j * j + k * k) if length < 2.0 - 1e-12: pairs.append(((i, j, k), length, math.degrees(math.acos(length / 2.0)))) pairs.sort(key=lambda pair: (pair[1], pair[0])) return pairs def sample(field, points, noise_ratio, normal): return [field.temperature(point) + noise_ratio * normal() for point in points] def summarise(values): ordered = sorted(values) count = len(ordered) mean = sum(ordered) / count deviation = math.sqrt(sum((value - mean) ** 2 for value in ordered) / (count - 1)) return {"count": count, "mean": round(mean, 4), "standardDeviation": round(deviation, 4), "percentile99": round(ordered[min(count - 1, int(math.ceil(0.99 * count)) - 1)], 4), "maximum": round(ordered[-1], 4)} def run(settings): box_side, samples = settings["boxSide"], settings["circleSamples"] field = PeriodicBumpField(box_side, settings["bumpWidth"], settings["bumpCount"], settings["seed"]) _, normal = normal_generator(settings["seed"] + 1) result = {"scriptVersion": SCRIPT_VERSION, "python": sys.version.split()[0], "units": {"length": "chi_rec (comoving distance to last scattering) = 1", "angle": "degree", "temperature": "signal RMS on the sphere = 1", "statistic": "dimensionless, -1 to 1"}} # 1. Periodicity: the field repeats after one lattice step (both points stay inside the stored images). probe = (0.5 * box_side, 0.03, -0.02) shifted = (probe[0] - box_side, probe[1], probe[2]) result["periodicityDifference"] = round(abs(field.temperature(probe) - field.temperature(shifted)), 12) # 2. Predicted pairs and their measured match. measured = [] for vector, length, alpha in lattice_pairs(box_side): axis = normalise(vector) first_points = circle_points(axis, alpha, samples) second_points = [tuple(p - box_side * v for p, v in zip(point, vector)) for point in first_points] first = sample(field, first_points, 0.0, normal) second = sample(field, second_points, 0.0, normal) best, shift = best_over_shifts(first, second) noisy_first = sample(field, first_points, settings["noiseRatio"], normal) noisy_second = sample(field, second_points, settings["noiseRatio"], normal) mirrored = sample(field, [tuple(p - box_side * v for p, v in zip(point, vector)) for point in circle_points(axis, alpha, samples, mirrored=True)], 0.0, normal) measured.append({"latticeVector": list(vector), "separation": round(length, 4), "alphaDegrees": round(alpha, 3), "statisticAtZeroShift": round(circle_statistic(first, second, 0), 4), "bestStatistic": round(best, 4), "bestShiftSamples": shift, "noisyStatisticAtZeroShift": round(circle_statistic(noisy_first, noisy_second, 0), 4), "noisyBestStatistic": round(best_over_shifts(noisy_first, noisy_second)[0], 4), "mirroredBestStatistic": round(best_over_shifts(first, mirrored)[0], 4)}) result["predictedPairs"] = measured xi2 = 1.0 / settings["noiseRatio"] ** 2 result["noiseExpectation"] = {"signalToNoise": round(1.0 / settings["noiseRatio"], 4), "css1998Equation5": round(xi2 / (1.0 + xi2), 4)} # 3. False matches: back-to-back circles about random axes, best over all shifts, # with the same noise as the noisy true pairs so the two are compared fairly. uniform = random.Random(settings["seed"] + 2) false_best, false_zero_shift = [], [] for _ in range(settings["falsePairTrials"]): axis = random_direction(uniform) alpha = uniform.choice(settings["falsePairRadiiDegrees"]) opposite = tuple(-component for component in axis) first = sample(field, circle_points(axis, alpha, samples), settings["noiseRatio"], normal) # Same frame vectors for both circles, as a pure translation would give. # Reflect each point through the plane normal to the axis: the circle about -axis. second_points = [tuple(p + 2.0 * math.cos(math.radians(alpha)) * o for p, o in zip(point, opposite)) for point in circle_points(axis, alpha, samples)] second = sample(field, second_points, settings["noiseRatio"], normal) false_best.append(best_over_shifts(first, second)[0]) false_zero_shift.append(circle_statistic(first, second, 0)) result["falsePairs"] = {"bestOverShifts": summarise(false_best), "singleShift": summarise(false_zero_shift), "independentSampleNullStandardDeviation": round(1.0 / math.sqrt(samples + 1), 4)} return result def main(): if len(sys.argv) != 3: raise SystemExit("usage: recompute_matched_circles.py INPUTS.json OUTPUT.json") with open(sys.argv[1], encoding="utf-8") as handle: inputs = json.load(handle) output = {"inputs": inputs["inputsVersion"], "runs": {}} for name, settings in inputs["runs"].items(): output["runs"][name] = run(settings) with open(sys.argv[2], "w", encoding="utf-8") as handle: json.dump(output, handle, indent=2) handle.write("\n") if __name__ == "__main__": main()