diff --git a/autofit/non_linear/bijector.py b/autofit/non_linear/bijector.py index 7654a8a99..3b9eaa287 100644 --- a/autofit/non_linear/bijector.py +++ b/autofit/non_linear/bijector.py @@ -230,7 +230,8 @@ def _check_resolved(self): raise RuntimeError( f"{type(self).__name__}.from_model(model) must be called before " "forward / inverse / log_det_jacobian / bounds_forward / kinds " - "-- the per-coordinate kind arrays are not yet resolved." + "/ identity_scales -- the per-coordinate kind arrays are not " + "yet resolved." ) @property @@ -239,6 +240,34 @@ def kinds(self) -> List[str]: self._check_resolved() return [_KIND_NAMES[int(code)] for code in self._kind_code] + @property + def identity_scales(self) -> List[Optional[float]]: + """ + The per-coordinate **linear** scale, or ``None`` where there is not one. + + Entry ``i`` is ``self._scale[i]`` when coordinate ``i`` is + ``identity``-kind — the one case in which ``forward`` is exactly + ``theta / scale``, an affine map — and ``None`` for ``log`` / ``logit``, + where ``scale`` is not read at all and no single linear factor describes + the coordinate. + + It exists so a consumer can ask "is this coordinate merely rescaled, + and by how much?" without reaching into ``_kind_code`` / ``_scale``. + :class:`~autofit.non_linear.clipper.ClipperPriorBoxJoint` is the caller: + a ball is a statement about physical coordinates, and it stays a ball + under a *common* linear rescale of both its members (radius ``R`` maps + to ``R / s``) but not otherwise — so it needs the scale where there is + one and an explicit "there is not" where there is not. ``None`` rather + than ``1.0`` for the non-identity kinds, deliberately: a ``log`` + coordinate is not a coordinate scaled by one, and reporting it as such + would let a caller compose with it silently and wrongly. + """ + self._check_resolved() + return [ + float(scale) if int(code) == _IDENTITY else None + for code, scale in zip(self._kind_code, self._scale) + ] + def forward(self, theta, xp=np): """ The physical-to-transformed map, ``theta -> phi``, vectorised over the diff --git a/autofit/non_linear/clipper.py b/autofit/non_linear/clipper.py index 170155e8e..af11a7f83 100644 --- a/autofit/non_linear/clipper.py +++ b/autofit/non_linear/clipper.py @@ -97,6 +97,11 @@ :meth:`ClipperPriorBox._inset_from_model`); the plain, no-bijector ``bounds_from_model`` path is unaffected — every existing caller keeps its current (physical) inset exactly. + +A **ball** composes with the same maps under a narrower condition, resolved per +pair rather than per map: it survives a common linear rescale of both its +members and nothing else. See :class:`ClipperPriorBoxJoint` ("Composition with +a scaler/bijector") and :meth:`ClipperPriorBoxJoint.pairs_in_stepped_coordinates`. """ import logging @@ -497,6 +502,43 @@ class ClipperPriorBoxJoint(ClipperPriorBox): mask's consumer is ``MultiStartGradient``'s momentum reset: a lane pushed out of the disk carries outward momentum in *both* coordinates, and zeroing only one of them leaves the pair spiralling back out on the next step. + + Composition with a scaler/bijector + --------------------------------- + + A ball is a statement about **physical** coordinates, so a search stepping + in some other coordinates can only be given a ball projection if the ball + is still a ball there. That is a question about the *pair*, not about the + map as a whole, and it has a clean answer: + + - both members ``identity``-kind under a **common** linear scale ``s``: the + step coordinates are ``(theta_i / s, theta_j / s)``, so + ``theta_i**2 + theta_j**2 <= R**2`` is exactly + ``phi_i**2 + phi_j**2 <= (R / s)**2``. A disk of radius ``R`` in physical + coordinates **is** a disk of radius ``R / s`` in the stepped ones, and + the radial shrink below — which is scale-covariant, being a pure + multiply — projects onto it correctly. The pair is accepted, with its + radius divided; + - anything else — ``log``, ``logit``, or two ``identity`` members with + *different* scales — is not a disk. Unequal scales give an **ellipse**, + whose nearest-point projection is not a radial shrink at all (it is the + root of a quartic), and the non-linear kinds give a region with no closed + form. Projecting a circle onto those coordinates would enforce a + different, silently wrong constraint, so :meth:`project` raises instead, + naming the offending index pair and the kinds it found. + + This is deliberately resolved **per pair** rather than per map. A model + reparameterised through ``log`` on, say, an ``einstein_radius`` while its + ``ell_comps`` stay linear is the common case, and refusing it wholesale is + what left those lanes at the box corner ``|e| = 1.414`` — inside every prior + box, outside the disk, and in the flat region the clipper exists to escape. + + The projection is **not** round-tripped through the bijector (map back to + physical, project, map forward). That would be correct for every kind, but + it would also drag the ``logit`` epsilon clamps across coordinates the ball + has nothing to do with, breaking the bit-identity the interior-point path + promises, saturating gradients at the clamped edges, and adding traced ops + to every step. Dividing the radius costs nothing and is exact. """ # Divisor floor for the radial shrink. Small enough to be irrelevant for any @@ -505,6 +547,93 @@ class ClipperPriorBoxJoint(ClipperPriorBox): # this), and large enough that `r / tiny` cannot overflow float32. _TINY = 1.0e-30 + def pairs_in_stepped_coordinates(self, pairs, scale=None, bijector=None): + """ + The declared ball pairs, restated in whatever coordinates the search is + actually stepping in — or a ``ValueError`` naming the pair that cannot + be restated. + + Parameters + ---------- + pairs + ``(index_0, index_1, radius)`` triples as + :meth:`~autofit.mapper.prior_model.abstract.AbstractPriorModel. + ball_constraint_index_pairs` returns them, in **physical** + coordinates. + scale + The per-parameter linear scale a ``scaler`` is applying, or + ``None``. + bijector + A **resolved** bijector (``from_model`` already called) the search + is stepping through, or ``None``. Mutually exclusive with ``scale``. + + Returns + ------- + The same triples with each ``radius`` divided by that pair's common + linear scale — unchanged when neither a scaler nor a bijector is in + play, since the common scale is then ``1.0``. + + Raises + ------ + ValueError + When a pair is not a disk in the stepped coordinates. See the class + docstring's "Composition with a scaler/bijector" for which pairs + those are and why the projection cannot be defined for them. + + This is a **model-resolution-time** check, separable from + :meth:`project` on purpose: ``MultiStartGradient`` calls it once, as + soon as it has resolved its bijector against the model and before it + compiles a step, so a combination that cannot work dies there rather + than a minute into a multi-hour fit. + """ + if scale is None and bijector is None: + return list(pairs) + + if scale is not None: + # A raw scaler is linear by construction, so every coordinate is + # `identity`-kind; only the equality of the two scales is in doubt. + scales = [float(value) for value in np.asarray(scale, dtype=float)] + kinds = ["identity"] * len(scales) + else: + scales = bijector.identity_scales + kinds = bijector.kinds + + resolved = [] + + for index_0, index_1, radius in pairs: + scale_0 = scales[index_0] + scale_1 = scales[index_1] + + common = ( + scale_0 is not None + and scale_1 is not None + and scale_0 == scale_1 + and scale_0 > 0.0 + and np.isfinite(scale_0) + ) + + if not common: + raise ValueError( + f"{type(self).__name__} cannot project onto the ball " + f"declared on parameter indices ({index_0}, {index_1}) " + "while the search steps in scaled or transformed " + f"coordinates: that pair maps as ({kinds[index_0]}, " + f"{kinds[index_1]}) with linear scales ({scale_0}, " + f"{scale_1}), which is not a disk. A ball survives a " + "COMMON linear scale `s` (radius `R` becomes `R / s`) and " + "is composed with; unequal scales give an ellipse and a " + "log/logit coordinate a region with no closed form, and " + "projecting a circle onto either would enforce a " + "different, silently wrong constraint. Leave this pair " + "identity-mapped under one common scale (e.g. a " + "`BijectorPerPath` carrying the non-linear kinds on the " + "OTHER paths), or use ClipperPriorBox instead." + ) + + resolved.append((index_0, index_1, radius / scale_0)) + + return resolved + def project(self, vector, model, xp=np, scale=None, bijector=None): projected, clipped_mask = super().project( vector=vector, @@ -519,20 +648,15 @@ def project(self, vector, model, xp=np, scale=None, bijector=None): if not pairs: return projected, clipped_mask - if scale is not None or bijector is not None: - # A ball is a statement about PHYSICAL coordinates, and neither - # change of variables preserves it: a per-parameter `scale` maps the - # disk to an ellipse, and a `bijector` maps it to something with no - # closed form at all. Projecting the scaled coordinates onto a circle - # would enforce a different, silently wrong constraint -- so this - # says so rather than doing it. - raise ValueError( - f"{type(self).__name__} cannot project onto a ball while the " - "search steps in scaled or transformed coordinates: the ball is " - "a statement about physical parameters, and neither a scaler nor " - "a bijector maps a disk to a disk. Use ClipperPriorBox with the " - "scaler/bijector, or drop them to project onto the ball." - ) + # Restated in the coordinates `projected` is already in (the box above + # returns the STEPPED vector on both the scaler and the bijector path), + # so the radial shrink below is unchanged -- it just shrinks onto a + # radius that has been divided by the pair's common scale. + pairs = self.pairs_in_stepped_coordinates( + pairs=pairs, + scale=scale, + bijector=bijector, + ) # One multiplicative factor per parameter, defaulting to 1.0, so the # whole projection is a single elementwise multiply and every diff --git a/autofit/non_linear/search/mle/multi_start_gradient/search.py b/autofit/non_linear/search/mle/multi_start_gradient/search.py index e732c46f7..a352a4d42 100644 --- a/autofit/non_linear/search/mle/multi_start_gradient/search.py +++ b/autofit/non_linear/search/mle/multi_start_gradient/search.py @@ -160,6 +160,17 @@ def __init__( enforcement on moves where the search converges and would shift every stored multi-start benchmark, so the default flip is deliberately a separate, re-baselined change. + + ``ClipperPriorBoxJoint`` additionally projects onto the balls the + model's classes declare (the ``ell_comps`` disk being the canonical + one). It composes with ``scaler`` / ``bijector`` **per ball pair**, + not wholesale: a pair both of whose members the map merely rescales + by one common linear factor ``s`` is still a disk in the stepped + coordinates, of radius ``R / s``, and is projected onto normally; a + pair a ``log`` / ``logit`` kind or two unequal scales touch is not a + disk at all and raises. That check runs once against the model at + the start of ``_fit`` — not at construction, where there is no model + to ask — and always before the first likelihood evaluation. scaler A per-parameter step scale derived from the priors, applied as a change of variables so the rule steps in ``phi = theta / scale`` while @@ -207,6 +218,13 @@ def __init__( Default ``BijectorNone`` -- a no-op, skipped entirely rather than applied as an identity, so the default path and its compiled step are both unchanged. + + Composable with ``ClipperPriorBoxJoint`` where the map leaves each + declared ball pair identity-mapped under one common scale — so a + model wanting ``log`` on (say) an ``einstein_radius`` while its + ``ell_comps`` step linearly gets both the reparameterisation and + the disk. See ``clipper`` above for the condition and where it is + checked. resurrect Restart-on-death. When ``True``, any start whose objective goes non-finite is redrawn each step (fresh params from the start band + @@ -362,26 +380,16 @@ def __init__( "once.)" ) - # Same reasoning, one rung along: a `ClipperPriorBoxJoint` projects onto - # a ball declared in PHYSICAL coordinates, and neither change of - # variables maps a disk to a disk (a diagonal scale gives an ellipse, a - # bijector gives something with no closed form). Surfaced here rather - # than at the first step, so a multi-hour fit does not die a minute in - # on a configuration that was wrong before it started. - if isinstance(self.clipper, ClipperPriorBoxJoint) and not ( - isinstance(self.scaler, ScalerNone) - and isinstance(self.bijector, BijectorNone) - ): - raise ValueError( - f"{type(self).__name__} received a " - f"{type(self.clipper).__name__} together with a non-default " - "`scaler` or `bijector`. The joint clipper projects onto a ball " - "declared in physical parameters, and neither change of " - "variables maps a disk to a disk -- projecting in the stepped " - "coordinates would enforce a different, silently wrong " - "constraint. Use ClipperPriorBox with the scaler/bijector, or " - "drop them to project onto the ball." - ) + # A `ClipperPriorBoxJoint` alongside a scaler/bijector is NOT rejected + # here, deliberately, even though the pair above is. Whether the two + # compose is a question about the MODEL -- which coordinate pairs carry + # a ball, and how the map treats each of them -- and no model exists at + # construction. A ball survives a common linear rescale of both its + # members and nothing else, so the same search+clipper+bijector triple + # is perfectly well defined on one model and undefinable on the next. + # The check is made once in `_fit`, as soon as the bijector has been + # resolved against the model and before any step is compiled, so a bad + # combination still dies before the first likelihood evaluation. self.reset_momentum_on_clip = reset_momentum_on_clip self.record_lane_nan_history = bool(record_lane_nan_history) @@ -884,6 +892,27 @@ def _fit( if has_bijector: self.bijector.from_model(model=model) + # The joint clipper's ball is declared in PHYSICAL coordinates, and it + # is still a ball in the stepped ones only where the map rescales both + # its members by one common linear factor (see + # `ClipperPriorBoxJoint.pairs_in_stepped_coordinates`). Resolved HERE: + # this is the first moment both halves are known -- the bijector was + # resolved against the model one line above, and the model is what says + # which pairs carry a ball at all -- and it is still before `_vmapped` + # is built, so an unresolvable combination dies without a single + # likelihood evaluation rather than a minute into a multi-hour fit. + # Its return value is discarded: the step loop calls `clipper.project`, + # which re-derives the same restated pairs itself. What is wanted here + # is only the raise. + if isinstance(self.clipper, ClipperPriorBoxJoint) and ( + has_scaler or has_bijector + ): + self.clipper.pairs_in_stepped_coordinates( + pairs=model.ball_constraint_index_pairs(), + scale=scale if has_scaler else None, + bijector=self.bijector if has_bijector else None, + ) + def _to_physical(vector): if scale_jnp is not None: return vector * scale_jnp diff --git a/test_autofit/non_linear/search/mle/test_multi_start_gradient.py b/test_autofit/non_linear/search/mle/test_multi_start_gradient.py index 1f8625653..00532ca6b 100644 --- a/test_autofit/non_linear/search/mle/test_multi_start_gradient.py +++ b/test_autofit/non_linear/search/mle/test_multi_start_gradient.py @@ -1350,3 +1350,78 @@ def test__the_default_clipper_leaves_a_lane_outside_the_disk_where_it_is(): for clipper in (af.ClipperNone(), af.ClipperPriorBox(margin=0.0)): projected, _ = clipper.project(vector=params, model=model, xp=np) assert np.hypot(projected[0, 0], projected[0, 1]) == pytest.approx(1.4) + + +def _log_intensity_lane_model(): + """`_lane_model()` with a `log`-eligible `intensity`, so a per-path bijector + has something to actually transform OUTSIDE the ball pair.""" + model = _lane_model() + model.intensity = af.LogUniformPrior(lower_limit=1.0e-3, upper_limit=1.0e3) + return model + + +def test__the_joint_clipper_and_a_bijector_construct_together(): + """No longer refused at construction: whether the two compose is a question + about the MODEL, and there is no model here to ask.""" + search = af.MultiStartAdam( + clipper=af.ClipperPriorBoxJoint(margin=0.0), + bijector=af.BijectorPerPath({"intensity": "log"}), + ) + + assert isinstance(search.clipper, af.ClipperPriorBoxJoint) + assert isinstance(search.bijector, af.BijectorPerPath) + + +def test__an_identity_mapped_ball_pair_is_accepted_at_model_resolution(): + """The check `_fit` runs once, immediately after resolving the bijector + against the model and before a step is compiled. Exercised here through the + search's own clipper rather than by calling `_fit`, which is JAX-native -- + this suite stays NumPy-only, and the call below is the identical one. + + `log` sits on `intensity`; the `ell_comps` pair is untouched by the map, so + the disk survives at its declared radius and the lane keeps the projection.""" + search = af.MultiStartAdam( + clipper=af.ClipperPriorBoxJoint(margin=0.0), + bijector=af.BijectorPerPath({"intensity": "log"}), + ) + model = _log_intensity_lane_model() + bijector = search.bijector.from_model(model=model) + + assert bijector.kinds == ["identity", "identity", "log"] + + pairs = search.clipper.pairs_in_stepped_coordinates( + pairs=model.ball_constraint_index_pairs(), + bijector=bijector, + ) + + assert pairs == [(0, 1, BALL_RADIUS)] + + # And the projection the step loop then makes still pulls the corner lane in. + params = np.array([[1.4 / np.sqrt(2.0), 1.4 / np.sqrt(2.0), np.log(5.0)]]) + projected, clipped_mask = search.clipper.project( + vector=params, model=model, xp=np, scale=None, bijector=bijector + ) + + assert np.hypot(projected[0, 0], projected[0, 1]) <= BALL_RADIUS + assert clipped_mask[0].tolist() == [True, True, False] + + +def test__a_logit_on_the_ell_comps_pair_is_refused_before_the_first_step(): + """The combination that genuinely cannot work still dies at model + resolution, so a multi-hour fit does not start and then fail a minute in.""" + search = af.MultiStartAdam( + clipper=af.ClipperPriorBoxJoint(margin=0.0), + bijector=af.BijectorLogit(), + ) + model = _lane_model() + bijector = search.bijector.from_model(model=model) + + assert bijector.kinds[:2] == ["logit", "logit"] + + with pytest.raises(ValueError) as exc_info: + search.clipper.pairs_in_stepped_coordinates( + pairs=model.ball_constraint_index_pairs(), + bijector=bijector, + ) + + assert "(0, 1)" in str(exc_info.value) diff --git a/test_autofit/non_linear/test_bijector.py b/test_autofit/non_linear/test_bijector.py index 9517758c5..57b8df4d5 100644 --- a/test_autofit/non_linear/test_bijector.py +++ b/test_autofit/non_linear/test_bijector.py @@ -409,6 +409,50 @@ def test__objective_composed_through_the_bijector_is_the_SAME_objective(): ) == pytest.approx(_log_coordinate_objective(theta), rel=1.0e-9) +def test__round_tripping_a_per_path_map_is_bit_exact_where_it_is_identity(): + """ + The F5 pin. A per-path map touches the paths it names and NOTHING else, and + "nothing else" has to mean bit-identity, not agreement to some tolerance. + + It is load-bearing for + :class:`~autofit.non_linear.clipper.ClipperPriorBoxJoint`, which composes + with such a map by dividing a ball's radius rather than by round-tripping + the vector through the bijector. That shortcut is only sound if the + untouched coordinates really are untouched: `theta / 1.0 * 1.0 == theta` + exactly in IEEE 754, so an identity coordinate survives forward-then-inverse + with every bit intact, while a `log` one does not (`exp(log(x))` is + correctly rounded, not exact) and is therefore checked to a tolerance. + + The composed objective is equal to `rel=1e-12` across the whole vector, so + the log coordinate's rounding is confirmed to be the only difference and to + be negligible. + """ + model = log_uniform_model() + bijector = BijectorPerPath({"normalization": "log"}).from_model(model) + + assert bijector.kinds == ["identity", "log", "identity"] + + rng = np.random.default_rng(5) + + for _ in range(50): + theta = rng.uniform( + np.array([-10.0, 0.5, -10.0]), np.array([10.0, 500.0, 10.0]) + ) + round_tripped = bijector.inverse(bijector.forward(theta)) + + # `==`, not `approx`: the identity coordinates are bit-exact. + assert round_tripped[0] == theta[0] + assert round_tripped[2] == theta[2] + + # The log coordinate round trips through `exp(log(x))`, which is + # correctly rounded rather than exact. + assert round_tripped[1] == pytest.approx(theta[1], rel=1.0e-12) + + assert _log_coordinate_objective(round_tripped) == pytest.approx( + _log_coordinate_objective(theta), rel=1.0e-12 + ) + + def test__bijector_cannot_move_the_MAP(): """ The end-to-end statement: minimise the raw objective, and minimise the diff --git a/test_autofit/non_linear/test_clipper.py b/test_autofit/non_linear/test_clipper.py index 9cc5305ea..2c232467f 100644 --- a/test_autofit/non_linear/test_clipper.py +++ b/test_autofit/non_linear/test_clipper.py @@ -7,7 +7,13 @@ from autofit import exc from autofit.mapper.identifier import Identifier from autofit.mapper.prior.tuple_prior import TuplePrior -from autofit.non_linear.bijector import BijectorAuto +from autofit.non_linear.bijector import ( + BijectorAuto, + BijectorDiagonal, + BijectorLogit, + BijectorNone, + BijectorPerPath, +) from autofit.non_linear.clipper import ( ClipperNone, ClipperPriorBox, @@ -697,28 +703,146 @@ def test__the_prior_box_is_still_enforced(self): assert projected[0] == pytest.approx(0.5) assert mask[0] - def test__a_scaler_or_bijector_is_refused_rather_than_applied(self): - """A ball is a statement about physical parameters; a diagonal scale maps - the disk to an ellipse and a bijector to something with no closed form. - Projecting in the stepped coordinates would enforce a different, silently - wrong constraint.""" + +class TestJointBallComposesWithAMap: + """A ball is a statement about PHYSICAL parameters, but it survives a common + linear rescale of both its members: a disk of radius `R` becomes a disk of + radius `R / s`. Only a genuinely non-linear pair -- a `log`/`logit` kind, or + two identity coordinates with DIFFERENT scales (an ellipse) -- has no + projection, and only those are refused.""" + + def test__the_no_op_bijector_composes_and_is_bit_identical(self): + """`BijectorNone` is every coordinate identity-scaled by 1.0, so passing + it must be indistinguishable from passing nothing -- not merely close.""" model = _ball_model() clipper = ClipperPriorBoxJoint(margin=0.0) + vector = np.array([-1.0, -1.0, 5.0]) - with pytest.raises(ValueError): + plain, plain_mask = clipper.project(vector=vector, model=model) + composed, composed_mask = clipper.project( + vector=vector, + model=model, + bijector=BijectorNone().from_model(model=model), + ) + + assert (composed == plain).all() + assert (composed_mask == plain_mask).all() + assert np.hypot(composed[0], composed[1]) == pytest.approx(BALL_RADIUS) + + def test__a_log_on_a_path_the_ball_does_not_touch_composes(self): + """The case the blanket refusal cost: a model reparameterised through + `log` on some OTHER parameter still has a linear `ell_comps` pair, so the + disk is still a disk and the corner is still projected onto it.""" + model = _ball_model() + model.intensity = af.LogUniformPrior(lower_limit=1.0e-3, upper_limit=1.0e3) + + bijector = BijectorPerPath({"intensity": "log"}).from_model(model=model) + assert bijector.kinds == ["identity", "identity", "log"] + + projected, mask = ClipperPriorBoxJoint(margin=0.0).project( + vector=np.array([-1.0, -1.0, np.log(5.0)]), + model=model, + bijector=bijector, + ) + + # The ball pair is untouched by the map, so the physical answer stands: + # the corner projected radially onto the disk, angle preserved. + assert np.hypot(projected[0], projected[1]) == pytest.approx(BALL_RADIUS) + assert projected[0] == pytest.approx(projected[1]) + assert mask[0] + assert mask[1] + + def test__a_logit_on_the_ball_pair_itself_is_refused_and_names_the_pair(self): + """`logit` maps the disk to a region with no closed form. The refusal has + to say WHICH pair and WHAT kinds, or a user with a large model cannot act + on it.""" + model = _ball_model() + clipper = ClipperPriorBoxJoint(margin=0.0) + bijector = BijectorLogit().from_model(model=model) + + assert bijector.kinds[:2] == ["logit", "logit"] + + with pytest.raises(ValueError) as exc_info: clipper.project( vector=np.array([-1.0, -1.0, 5.0]), model=model, - scale=np.array([1.0, 1.0, 1.0]), + bijector=bijector, ) - with pytest.raises(ValueError): - clipper.project( - vector=np.array([-1.0, -1.0, 5.0]), + message = str(exc_info.value) + assert "(0, 1)" in message + assert "logit" in message + + def test__a_common_scale_projects_onto_the_divided_radius(self): + """The composition rule itself: in coordinates scaled by `s`, the disk of + radius `R` IS the disk of radius `R / s`.""" + model = _ball_model() + + projected, mask = ClipperPriorBoxJoint(margin=0.0).project( + # The stepped vector: the physical corner `(-1, -1)` divided by 2. + vector=np.array([-0.5, -0.5, 5.0]), + model=model, + scale=np.array([2.0, 2.0, 1.0]), + ) + + assert np.hypot(projected[0], projected[1]) == pytest.approx(BALL_RADIUS / 2.0) + assert mask[0] + assert mask[1] + + def test__an_unequal_scale_on_the_pair_is_refused(self): + """Two different scales map the disk to an ELLIPSE, whose nearest-point + projection is not a radial shrink at all.""" + model = _ball_model() + + with pytest.raises(ValueError) as exc_info: + ClipperPriorBoxJoint(margin=0.0).project( + vector=np.array([-0.5, -0.5, 5.0]), model=model, - bijector=BijectorAuto().from_model(model=model), + scale=np.array([2.0, 3.0, 1.0]), ) + assert "(0, 1)" in str(exc_info.value) + + def test__a_diagonal_bijector_with_equal_pair_widths_is_accepted(self): + """`BijectorDiagonal(ScalerPriorWidth())` reports every coordinate as + `identity`, and the two `ell_comps` priors have equal widths -- so the + pair shares one scale and composes.""" + model = _ball_model() + bijector = BijectorDiagonal(af.ScalerPriorWidth()).from_model(model=model) + + scales = bijector.identity_scales + assert scales[0] == scales[1] + assert scales[0] is not None + + projected, mask = ClipperPriorBoxJoint(margin=0.0).project( + vector=np.array([-1.0, -1.0, 5.0]) / scales[0], + model=model, + bijector=bijector, + ) + + assert np.hypot(projected[0], projected[1]) == pytest.approx( + BALL_RADIUS / scales[0] + ) + assert mask[0] + assert mask[1] + + def test__a_model_with_no_ball_is_unaffected_by_any_of_this(self): + """The resolution is per DECLARED pair, so a model declaring none never + reaches it -- a pipeline may configure the joint clipper globally and + still use any bijector it likes.""" + model = _model( + alpha=af.UniformPrior(lower_limit=-1.0, upper_limit=1.0), + beta=af.UniformPrior(lower_limit=-1.0, upper_limit=1.0), + ) + + projected, _ = ClipperPriorBoxJoint(margin=0.0).project( + vector=np.array([5.0, -5.0]), + model=model, + bijector=BijectorLogit().from_model(model=model), + ) + + assert np.isfinite(projected).all() + class TestJointBallIdentifier: def test__the_base_clippers_identifier_fields_are_pinned(self): @@ -739,6 +863,34 @@ def test__the_joint_clipper_is_a_different_identifier(self): assert str(box) != str(joint) + def test__composing_with_a_bijector_leaves_the_identifier_unchanged(self): + """The composition adds no constructor argument, deliberately: had it + taken one, every stored `ClipperPriorBoxJoint` result would have been + silently re-keyed by a change that does not alter what those runs did. + + Note what this does NOT say. The equality holds partly because the + search identifier does not yet tag the `bijector` at all -- a separate, + pre-existing gap that predates this change and is tracked on its own; + this test pins that the clipper's own key is untouched, not that two + runs differing only in bijector *should* share an output directory.""" + model = _ball_model() + + plain = Identifier( + [af.MultiStartAdam(clipper=ClipperPriorBoxJoint()), model, None] + ) + composed = Identifier( + [ + af.MultiStartAdam( + clipper=ClipperPriorBoxJoint(), + bijector=BijectorPerPath({"intensity": "log"}), + ), + model, + None, + ] + ) + + assert str(plain) == str(composed) + # The `jit`/`grad`/`vmap` behaviour of the radial shrink -- in particular that the # "double where" keeps `grad` finite at the origin, where `sqrt(0)` has an infinite @@ -771,17 +923,51 @@ def test__lbfgs_accepts_it_on_a_model_that_declares_no_ball(self): assert bounds.lb == pytest.approx(search.clipper.bounds_from_model(model)[0]) - def test__a_scaler_is_refused_at_construction_not_at_the_first_step(self): - """A multi-hour fit must not die a minute in on a configuration that was - wrong before it started.""" + def test__a_scaler_is_accepted_at_construction_because_no_model_exists_yet(self): + """Whether a ball composes with a map is a question about the MODEL -- + which pairs carry a ball, and how the map treats each of them -- and no + model exists at construction. The same triple is well defined on one + model and undefinable on the next, so construction cannot answer it.""" + search = af.MultiStartAdam( + clipper=ClipperPriorBoxJoint(), scaler=af.ScalerPriorWidth() + ) + + assert isinstance(search.clipper, ClipperPriorBoxJoint) + assert isinstance(search.scaler, af.ScalerPriorWidth) + + def test__a_bijector_is_accepted_at_construction_too(self): + search = af.MultiStartAdam( + clipper=ClipperPriorBoxJoint(), bijector=BijectorAuto() + ) + + assert isinstance(search.clipper, ClipperPriorBoxJoint) + assert isinstance(search.bijector, BijectorAuto) + + def test__an_unresolvable_map_is_refused_at_model_resolution(self): + """Refused before the first STEP, not at construction: a multi-hour fit + must still not die a minute in on a combination that was wrong before it + started, so the check runs once against the model at the top of `_fit`, + before anything is compiled or evaluated.""" + model = _ball_model() + clipper = ClipperPriorBoxJoint() + with pytest.raises(ValueError): - af.MultiStartAdam( - clipper=ClipperPriorBoxJoint(), scaler=af.ScalerPriorWidth() + clipper.pairs_in_stepped_coordinates( + pairs=model.ball_constraint_index_pairs(), + bijector=BijectorLogit().from_model(model=model), ) - def test__a_bijector_is_refused_at_construction_too(self): - with pytest.raises(ValueError): - af.MultiStartAdam(clipper=ClipperPriorBoxJoint(), bijector=BijectorAuto()) + def test__a_resolvable_map_passes_model_resolution(self): + model = _ball_model() + model.intensity = af.LogUniformPrior(lower_limit=1.0e-3, upper_limit=1.0e3) + clipper = ClipperPriorBoxJoint() + + pairs = clipper.pairs_in_stepped_coordinates( + pairs=model.ball_constraint_index_pairs(), + bijector=BijectorPerPath({"intensity": "log"}).from_model(model=model), + ) + + assert pairs == [(0, 1, BALL_RADIUS)] def test__the_plain_box_clipper_still_composes_with_a_scaler(self): """The refusal is scoped to the joint clipper: nothing about the existing