| """Research implementation of exact-evidence certificates for opaque spheres. |
| |
| The mathematical guarantee is in real arithmetic. Float64 plus conservative |
| tolerance is tested here, not a formally verified interval-arithmetic kernel. |
| Only renderer-authoritative geometry and fixed receiver/emitter identities |
| are accepted. Unreported geometry changes invalidate the contract. |
| """ |
| import numpy as np |
|
|
|
|
| def visibility_certificate(scene, points, lights, guard=1e-9): |
| """One segment query returns visibility and a conservative clearance. |
| |
| Distance to the trimmed segment is evaluated for every sphere. This costs |
| more arithmetic than a Boolean early-out shadow query; timings include it. |
| """ |
| p, l = np.broadcast_arrays(np.asarray(points, float), np.asarray(lights, float)) |
| if p.shape[-1] != 3 or not np.isfinite(p).all() or not np.isfinite(l).all(): |
| raise ValueError("Finite 3D segment endpoints required") |
| d = l-p |
| length2 = (d*d).sum(-1) |
| if (length2 <= 0).any() or guard < 0: |
| raise ValueError("Positive segment length and nonnegative guard required") |
| clearance = np.full(length2.shape, np.inf) |
| for sphere in scene.spheres: |
| a = np.clip(((sphere[:3]-p)*d).sum(-1)/length2, 1e-6, 1-1e-6) |
| signed = np.linalg.norm(p+a[..., None]*d-sphere[:3], axis=-1)-sphere[3] |
| clearance = np.minimum(clearance, signed) |
| v = (clearance > 0).astype(float) |
| margin = np.maximum(np.abs(clearance)-guard, 0) |
| ambiguous = np.abs(clearance) <= guard |
| if np.any(ambiguous): |
| |
| |
| |
| shape = v.shape |
| flat = v.reshape(-1).copy(); mask = ambiguous.reshape(-1) |
| flat[mask] = scene.visibility(p.reshape(-1,3)[mask], l.reshape(-1,3)[mask]) |
| v = flat.reshape(shape) |
| return v, np.nextafter(margin, 0) |
|
|
|
|
| class CertificateMemory: |
| """Canonical evidence plus a cumulative geometry-motion budget. |
| |
| Clearance shrinks by at most max_s(||delta center_s||+|delta radius_s|) |
| per update. The accumulated bound is conservative, including on return paths. |
| Light endpoints are fixed. New receiver positions can query old certificates |
| with an additional Euclidean displacement bound, but cannot commit to old IDs. |
| """ |
| def __init__(self, points, lights, spheres, namespace="scene", mode="margin"): |
| self.points = np.array(points, float, copy=True) |
| self.lights = np.array(lights, float, copy=True) |
| self.geometry = np.array(spheres, float, copy=True) |
| if (self.points.ndim != 2 or self.points.shape[1] != 3 or len(self.points) < 1 |
| or self.lights.ndim != 2 or self.lights.shape[1] != 3 or len(self.lights) < 1 |
| or self.geometry.ndim != 2 or self.geometry.shape[1] != 4 or len(self.geometry) < 1 |
| or not all(np.isfinite(a).all() for a in (self.points, self.lights, self.geometry)) |
| or (self.geometry[:, 3] <= 0).any() or not namespace |
| or mode not in ("margin", "epoch", "unsafe")): |
| raise ValueError("Invalid canonical certificate domain") |
| self.namespace, self.mode = str(namespace), mode |
| shape = (len(self.points), len(self.lights)) |
| self.values = np.full(shape, np.nan, np.float32) |
| self.margins = np.zeros(shape, np.float64) |
| self.stamps = np.zeros(shape, np.float64) |
| self.budget = 0.0 |
|
|
| @property |
| def nbytes(self): |
| return sum(a.nbytes for a in (self.points, self.lights, self.geometry, |
| self.values, self.margins, self.stamps)) |
|
|
| def _ids(self, ids): |
| x = np.asarray(ids) |
| if x.ndim != 1 or not np.issubdtype(x.dtype, np.integer) or (x < 0).any() or (x >= len(self.points)).any(): |
| raise ValueError("Invalid canonical receiver IDs") |
| return x |
|
|
| def begin_geometry(self, spheres, namespace=None): |
| g = np.array(spheres, float, copy=True) |
| if (g.shape != self.geometry.shape or not np.isfinite(g).all() |
| or (g[:, 3] <= 0).any()): |
| raise ValueError("Topology change requires a new canonical memory") |
| if namespace is not None and str(namespace) != self.namespace: |
| raise ValueError("Namespace mismatch: create fresh memory") |
| displacement = np.max(np.linalg.norm(g[:, :3]-self.geometry[:, :3], axis=1) |
| + np.abs(g[:, 3]-self.geometry[:, 3])) |
| if displacement > 0: |
| self.budget = float(np.nextafter(self.budget+displacement, np.inf)) |
| if self.mode == "epoch": |
| self.values.fill(np.nan) |
| self.geometry = g |
|
|
| def lookup(self, ids, query_points=None, extra_motion=0.0): |
| ids = self._ids(ids) |
| delta = np.zeros(len(ids)) |
| if query_points is not None: |
| p = np.asarray(query_points, float) |
| if p.shape != self.points[ids].shape or not np.isfinite(p).all(): |
| raise ValueError("Query points must match canonical anchor IDs") |
| delta = np.linalg.norm(p-self.points[ids], axis=1) |
| if not np.isfinite(extra_motion) or extra_motion < 0: |
| raise ValueError("Nonnegative future-motion bound required") |
| movement = self.budget-self.stamps[ids]+delta[:, None]+extra_motion |
| known = np.isfinite(self.values[ids]) |
| if self.mode != "unsafe": |
| known &= (self.margins[ids] > movement) | (movement == 0) |
| return np.nan_to_num(self.values[ids], nan=0).astype(float), known |
|
|
| def commit(self, ids, indices, visibility, margins): |
| ids = self._ids(ids) |
| j, v, m = np.asarray(indices), np.asarray(visibility), np.asarray(margins) |
| if (j.ndim != 1 or j.shape != ids.shape or not np.issubdtype(j.dtype, np.integer) |
| or v.shape != j.shape or m.shape != j.shape or (j < 0).any() |
| or (j >= len(self.lights)).any() or not np.isfinite(v).all() |
| or ((v != 0) & (v != 1)).any() or not np.isfinite(m).all() or (m < 0).any()): |
| raise ValueError("Expected one exact observation and margin per receiver") |
| self.values[ids, j] = v |
| self.margins[ids, j] = m |
| self.stamps[ids, j] = self.budget |
|
|
| def save(self, path): |
| np.savez_compressed(path, points=self.points, lights=self.lights, geometry=self.geometry, |
| values=self.values, margins=self.margins, stamps=self.stamps, |
| budget=np.float64(self.budget), namespace=np.array(self.namespace), mode=np.array(self.mode)) |
|
|
| @classmethod |
| def load(cls, path, namespace): |
| with np.load(path, allow_pickle=False) as d: |
| if str(d['namespace']) != str(namespace): |
| raise ValueError("Namespace mismatch") |
| obj = cls(d['points'], d['lights'], d['geometry'], namespace, str(d['mode'])) |
| shape = obj.values.shape |
| v, m, s = d['values'], d['margins'], d['stamps'] |
| budget = float(d['budget']) |
| if (v.shape != shape or m.shape != shape or s.shape != shape or v.dtype != np.float32 |
| or not np.isfinite(m).all() or not np.isfinite(s).all() or (m < 0).any() |
| or not np.isfinite(budget) or budget < 0 or (s < 0).any() or (s > budget).any() |
| or np.isinf(v).any() or ((v[np.isfinite(v)] != 0) & (v[np.isfinite(v)] != 1)).any()): |
| raise ValueError("Malformed certificate checkpoint") |
| obj.values, obj.margins, obj.stamps = v.copy(), m.copy(), s.copy() |
| obj.budget = budget |
| return obj |
|
|
|
|
| def enclosure(bound, values, known): |
| """Deterministic linear RGB interval, conditional on valid certificates.""" |
| b, v, k = np.asarray(bound), np.asarray(values), np.asarray(known, bool) |
| if b.ndim != 3 or v.shape != b.shape[:2] or k.shape != v.shape or not np.isfinite(b).all() or (b < 0).any(): |
| raise ValueError("Expected nonnegative contribution bounds and known mask") |
| lo = (b*np.where(k, v, 0)[..., None]).sum(1) |
| return lo, lo+(b*(~k)[..., None]).sum(1) |
|
|
