python/dpt/spectral_recovery.py

1"""Multi-view spectral pose recovery with constrained shared acquisition parameters.23This adapter composes the canonical material projector, spectral response,4spatial response, calibration and objective. The optimiser sees a fixed SE(3)5chart followed by named nuisance blocks. Detector-sized data stay on CUDA;6pinned control buffers move only pose, nuisance values, losses, gradients and7status flags. No observed image is synthesised or normalised by this adapter.8"""910from __future__ import annotations1112import math13from dataclasses import dataclass, field14from typing import Any, Literal1516from dpt._runtime import prepare_context, require_no_tape17from dpt.contracts import ContractError, NumericalError, binary32_scalar, finite_scalar, integer18from dpt.detector import (19    BlurSpec,20    BlurWorkspace,21    CalibrationSpec,22    DetectorWorkspace,23    blur,24    blur_transpose,25    calibrate,26    calibration_vjp,27    prepare_blur,28    prepare_detector,29)30from dpt.geometry import DetectorGeometry, RigidTransform31from dpt.material_projection import (32    MaterialProjectionWorkspace,33    material_projection_vjp,34    prepare_material_projection,35    project_material_paths,36)37from dpt.materials import Provenance38from dpt.objectives import ObjectiveSpec, ObjectiveWorkspace, evaluate_objective, prepare_objective39from dpt.projection import ProjectionSpec40from dpt.recovery import PoseChart41from dpt.registration import Evaluation, RecoveryPolicy, RecoveryResult, Vector, recover_parameters42from dpt.spectral import (43    SpectralSpec,44    SpectralWorkspace,45    prepare_spectral,46    spectral_signal,47    spectral_vjp,48)49from dpt.volumes import GridSpec5051# Public Python boundaries validate runtime inputs.52# pyright: reportUnnecessaryIsInstance=false535455def _physical_parameter(value: float, name: str, *, positive: bool = False) -> float:56    """Translate shared scalar storage checks into rejectable trial-domain errors."""57    try:58        rounded = binary32_scalar(value, f"trial {name}")59    except ContractError as error:60        raise NumericalError(str(error)) from error61    if positive and rounded <= 0:62        raise NumericalError(f"trial {name} lies outside its representable physical domain")63    return rounded646566def _scaled_product(*values: float) -> float:67    """Preserve exponent range until the final binary64 chart product is formed."""68    mantissa, exponent = 1.0, 069    for value in values:70        if not math.isfinite(value):71            raise NumericalError("non-finite physical derivative or chart scale")72        part, power = math.frexp(value)73        mantissa *= part74        exponent += power75    try:76        return math.ldexp(mantissa, exponent)77    except OverflowError as error:78        raise NumericalError("chart derivative overflows binary64") from error798081# region book:calibration-parameter-chart82@dataclass(frozen=True, slots=True)83class CalibrationBlock:84    """Named acquisition group sharing one gain/exposure/offset across its views.8586    Exactly one multiplicative scale may be fitted. Its positive chart is87    reference*exp(scale_step*z); the other scale is fixed to remove the gain /88    exposure gauge. An active offset is reference+offset_step*z in output units.89    Optional Gaussian priors act on these dimensionless chart coordinates.9091    Device parameters are rounded to binary32. The VJP differentiates the smooth92    physical chart at those represented values, not the discontinuous rounding93    map. Shared physical derivatives stay binary64 until the chart product is94    formed, so a large positive scale can preserve a small logarithmic partial.95    """9697    name: str98    gain: float = 1.099    exposure: float = 1.0100    offset: float = 0.0101    fit_scale: Literal["none", "gain", "exposure"] = "none"102    fit_offset: bool = False103    scale_step: float = 1.0104    offset_step: float = 1.0105    scale_prior_precision: float = 0.0106    offset_prior_precision: float = 0.0107108    def __post_init__(self) -> None:109        if not self.name.strip():110            raise ContractError("calibration groups need non-empty names")111        if self.fit_scale not in ("none", "gain", "exposure") or type(self.fit_offset) is not bool:112            raise ContractError("invalid calibration activity specification")113        for name in ("gain", "exposure", "scale_step", "offset_step"):114            if finite_scalar(getattr(self, name), name, minimum=0) <= 0:115                raise ContractError(f"{name} must be strictly positive in the recovery chart")116        finite_scalar(self.offset, "offset")117        for name in ("scale_prior_precision", "offset_prior_precision"):118            finite_scalar(getattr(self, name), name, minimum=0)119        if self.fit_scale == "none" and self.scale_prior_precision:120            raise ContractError("a scale prior requires an active scale")121        if not self.fit_offset and self.offset_prior_precision:122            raise ContractError("an offset prior requires an active offset")123        for name in ("gain", "exposure", "offset"):124            try:125                _physical_parameter(getattr(self, name), name, positive=name != "offset")126            except NumericalError as error:127                raise ContractError(str(error)) from error128129    @property130    def dimension(self) -> int:131        return int(self.fit_scale != "none") + int(self.fit_offset)132133    def decode(self, coordinates: Vector) -> tuple[float, float, float]:134        if len(coordinates) != self.dimension:135            raise ContractError("nuisance coordinates differ from the declared group chart")136        if not all(math.isfinite(value) for value in coordinates):137            raise NumericalError("trial nuisance coordinates are non-finite")138        gain, exposure, offset = self.gain, self.exposure, self.offset139        cursor = 0140        if self.fit_scale != "none":141            try:142                scale = math.exp(self.scale_step * coordinates[cursor])143            except OverflowError as error:144                raise NumericalError("trial logarithmic calibration scale overflowed") from error145            if self.fit_scale == "gain":146                gain *= scale147            else:148                exposure *= scale149            cursor += 1150        if self.fit_offset:151            offset += self.offset_step * coordinates[cursor]152        return (153            _physical_parameter(gain, "gain", positive=True),154            _physical_parameter(exposure, "exposure", positive=True),155            _physical_parameter(offset, "offset"),156        )157158    def gradient(self, coordinates: Vector, physical_gradient: Vector) -> Vector:159        if len(physical_gradient) != 3:160            raise ContractError("a calibration gradient needs gain, exposure and offset partials")161        gain, exposure, _ = self.decode(coordinates)162        gradient: list[float] = []163        if self.fit_scale != "none":164            index, value = (0, gain) if self.fit_scale == "gain" else (1, exposure)165            gradient.append(_scaled_product(self.scale_step, value, physical_gradient[index]))166        if self.fit_offset:167            gradient.append(_scaled_product(self.offset_step, physical_gradient[2]))168        return tuple(gradient)169170    def prior(self, coordinates: Vector) -> Evaluation:171        if len(coordinates) != self.dimension:172            raise ContractError("prior coordinates differ from calibration chart")173        precisions = () if self.fit_scale == "none" else (self.scale_prior_precision,)174        precisions += (self.offset_prior_precision,) if self.fit_offset else ()175        return Evaluation(176            math.fsum(0.5 * p * z * z for p, z in zip(precisions, coordinates, strict=True)),177            tuple(p * z for p, z in zip(precisions, coordinates, strict=True)),178        )179180181# endregion book:calibration-parameter-chart182183184@dataclass(frozen=True, slots=True)185class SpectralView:186    name: str187    geometry: DetectorGeometry188    spectral: SpectralSpec189    coefficients: Any190    weights: Any191    response: Any192    observation: Any193    calibration_group: str194    observation_provenance: Provenance195    objective: ObjectiveSpec = field(default_factory=ObjectiveSpec)196    objective_weights: Any = None197    objective_weight: float = 1.0198    spatial_response: BlurSpec | None = None199200    def __post_init__(self) -> None:201        if self.spectral.precision != "float32" or self.objective.precision != "float32":202            raise ContractError(203                "spectral pose recovery requires float32 spectral and objective precision; "204                "its calibration and optional spatial-response buffers are binary32"205            )206        if not isinstance(self.observation_provenance, Provenance):207            raise ContractError("each observation requires acquisition provenance")208        if not self.name.strip() or not self.calibration_group.strip():209            raise ContractError("view and calibration-group names are required")210        if finite_scalar(self.objective_weight, "view objective weight", minimum=0) <= 0:211            raise ContractError("view objective weight must be strictly positive")212        if self.objective.domain == "log_transmission":213            raise ContractError(214                "spectral recovery compares calibrated signal or counts, not log data"215            )216        if (217            self.spatial_response is not None218            and (self.spatial_response.height, self.spatial_response.width) != self.geometry.shape219        ):220            raise ContractError("spatial response and detector shapes differ")221        if not self.spectral.active_paths:222            raise ContractError("spectral pose recovery needs active material paths")223        if (224            self.spectral.active_weights225            or self.spectral.active_response226            or self.spectral.active_coefficients227        ):228            raise ContractError(229                "this recovery adapter keeps spectrum, response and coefficients fixed"230            )231        if not self.objective.weighted and self.objective_weights is not None:232            raise ContractError("objective weights require an explicitly weighted objective")233234235@dataclass(frozen=True, slots=True)236class SpectralPoseProblem:237    grid: GridSpec238    materials: int239    fields: Any240    fields_provenance: Provenance241    views: tuple[SpectralView, ...]242    calibration: tuple[CalibrationBlock, ...]243    samples_per_ray: int = 256244245    def __post_init__(self) -> None:246        integer(self.materials, "materials", minimum=1, maximum=32)247        integer(self.samples_per_ray, "samples_per_ray", minimum=1)248        if not isinstance(self.fields_provenance, Provenance):249            raise ContractError("fixed material fields require provenance")250        if not isinstance(self.views, tuple) or not self.views:251            raise ContractError(252                "a spectral recovery problem needs an immutable non-empty view tuple"253            )254        if not isinstance(self.calibration, tuple) or not self.calibration:255            raise ContractError("calibration blocks must be an immutable non-empty tuple")256        if len({view.name for view in self.views}) != len(self.views):257            raise ContractError("view names must be unique")258        groups = {block.name: block for block in self.calibration}259        if len(groups) != len(self.calibration):260            raise ContractError("calibration group names must be unique")261        if set(groups) != {view.calibration_group for view in self.views}:262            raise ContractError("every named calibration group must be defined and used")263        for group in groups:264            units = {265                view.spectral.output_unit for view in self.views if view.calibration_group == group266            }267            if len(units) != 1:268                raise ContractError(269                    "shared calibration groups require identical detector output units"270                )271        for view in self.views:272            if view.objective.domain == "counts" and view.spectral.output_unit != "counts":273                raise ContractError("count-domain views must declare spectral output_unit='counts'")274            if view.spectral.materials != self.materials:275                raise ContractError("all views must use the same fixed material basis")276            block = groups[view.calibration_group]277            if view.objective.kind == "poisson" and (278                block.gain != 1.0279                or block.offset != 0.0280                or block.fit_scale == "gain"281                or block.fit_offset282                or view.spatial_response is not None283            ):284                raise ContractError(285                    "Poisson counts require unit gain, zero offset and no spatial blur"286                )287288289@dataclass(slots=True)290class _PreparedView:291    view: SpectralView292    group: int293    material: MaterialProjectionWorkspace294    spectral: SpectralWorkspace295    detector: DetectorWorkspace296    objective: ObjectiveWorkspace297    blur: BlurWorkspace | None298    calibration_spec: CalibrationSpec299    mean: Any300    spread: Any301    prediction: Any302    image_seed: Any303    spread_seed: Any304    mean_seed: Any305    loss: Any306    pose_gradient: Any307    nuisance_gradients: tuple[Any, Any, Any]308309310@dataclass(frozen=True, slots=True)311class CalibrationEstimate:312    name: str313    gain: float314    exposure: float315    offset: float316317318@dataclass(frozen=True, slots=True)319class SpectralRecoveryResult:320    pose: RigidTransform321    calibration: tuple[CalibrationEstimate, ...]322    optimisation: RecoveryResult323324325class SpectralPoseEvaluator:326    """Prepared single-stream multi-view objective with manual first-order products.327328    Fields, input spectra, responses, observations and weights are immutable for329    this object's lifetime. Preparation checks their current values. Numerical330    flags are checked through each workspace after evaluation. Every path,331    image and adjoint buffer is allocated once. Concurrent/re-entrant evaluation332    and ambient tape recording are rejected before any staging buffer is changed.333    """334335    def __init__(336        self,337        problem: SpectralPoseProblem,338        chart: PoseChart,339        *,340        device: str = "cuda:0",341        stream: Any = None,342    ) -> None:343        require_no_tape()344        self.problem, self.chart = problem, chart345        self.context = prepare_context(device=device, stream=stream)346        ctx, wp = self.context, self.context.wp347        ctx.array(348            problem.fields,349            "fields",350            dtype=wp.float32,351            shape=(problem.materials * problem.grid.voxels,),352        )353        self._group_slices: tuple[slice, ...] = self._parameter_slices()354        self.dimension = 6 + sum(block.dimension for block in problem.calibration)355        self._running = False356        self.last_view_losses: tuple[float, ...] | None = None357        views = len(problem.views)358        with ctx.scope():359            self._pose_device = wp.empty(12, dtype=wp.float64, device=ctx.device)360            self._parameters_device = wp.empty(361                3 * len(problem.calibration), dtype=wp.float32, device=ctx.device362            )363            self._control_device = wp.empty(7 * views, dtype=wp.float64, device=ctx.device)364            self._nuisance_device = wp.zeros(3 * views, dtype=wp.float64, device=ctx.device)365        self._parameter_views = tuple(366            tuple(self._parameters_device[3 * i + j : 3 * i + j + 1] for j in range(3))367            for i in range(len(problem.calibration))368        )369        self._pose_host = wp.empty(12, dtype=wp.float64, device="cpu", pinned=True)370        self._parameters_host = wp.empty(371            3 * len(problem.calibration), dtype=wp.float32, device="cpu", pinned=True372        )373        self._control_host = wp.empty(7 * views, dtype=wp.float64, device="cpu", pinned=True)374        self._nuisance_host = wp.empty(3 * views, dtype=wp.float64, device="cpu", pinned=True)375        self._pose_values = self._pose_host.numpy()376        self._parameter_values = self._parameters_host.numpy()377        self._control_values = self._control_host.numpy()378        self._nuisance_values = self._nuisance_host.numpy()379        self._prepared = tuple(380            self._prepare_view(index, view) for index, view in enumerate(problem.views)381        )382        self._validate_fixed_inputs()383384    def _parameter_slices(self) -> tuple[slice, ...]:385        cursor = 6386        slices: list[slice] = []387        for block in self.problem.calibration:388            slices.append(slice(cursor, cursor + block.dimension))389            cursor += block.dimension390        return tuple(slices)391392    def _prepare_view(self, index: int, view: SpectralView) -> _PreparedView:393        ctx, wp = self.context, self.context.wp394        group = next(395            i396            for i, block in enumerate(self.problem.calibration)397            if block.name == view.calibration_group398        )399        block = self.problem.calibration[group]400        pixels, spec = view.geometry.pixels, view.spectral401        ctx.array(402            view.coefficients,403            "coefficients",404            dtype=wp.float32,405            shape=(spec.materials * spec.energies,),406        )407        ctx.array(408            view.weights,409            "weights",410            dtype=wp.float32,411            shape=(spec.energies * (1 if spec.shared_weights else pixels),),412        )413        ctx.array(414            view.response,415            "response",416            dtype=wp.float32,417            shape=(spec.energies * (1 if spec.shared_response else pixels),),418        )419        ctx.array(view.observation, "observation", dtype=wp.float32, shape=(pixels,))420        if view.objective.weighted:421            ctx.array(422                view.objective_weights, "objective_weights", dtype=wp.float32, shape=(pixels,)423            )424        with ctx.scope():425            paths = wp.empty(spec.materials * pixels, dtype=wp.float32, device=ctx.device)426            adj_paths = wp.empty_like(paths)427            mean = wp.empty(pixels, dtype=wp.float32, device=ctx.device)428            mean_seed = wp.empty_like(mean)429            spread = wp.empty_like(mean) if view.spatial_response is not None else mean430            spread_seed = wp.empty_like(mean) if view.spatial_response is not None else mean_seed431            prediction, image_seed = wp.empty_like(mean), wp.empty_like(mean)432        material = prepare_material_projection(433            self.problem.grid,434            view.geometry,435            ProjectionSpec(self.problem.samples_per_ray),436            materials=self.problem.materials,437            fields=self.problem.fields,438            out_paths=paths,439            adj_paths=adj_paths,440            device=ctx.device,441            stream=ctx.stream,442        )443        spectral = prepare_spectral(spec, max_pixels=pixels, device=ctx.device, stream=ctx.stream)444        detector = prepare_detector(max_pixels=pixels, device=ctx.device, stream=ctx.stream)445        objective = prepare_objective(446            view.objective, max_pixels=pixels, device=ctx.device, stream=ctx.stream447        )448        spatial = (449            prepare_blur(view.spatial_response, device=ctx.device, stream=ctx.stream)450            if view.spatial_response is not None451            else None452        )453        calibration_spec = CalibrationSpec(454            spec.output_unit,455            active_gain=block.fit_scale == "gain",456            active_exposure=block.fit_scale == "exposure",457            active_offset=block.fit_offset,458        )459        return _PreparedView(460            view,461            group,462            material,463            spectral,464            detector,465            objective,466            spatial,467            calibration_spec,468            mean,469            spread,470            prediction,471            image_seed,472            spread_seed,473            mean_seed,474            self._control_device[7 * index : 7 * index + 1],475            self._control_device[7 * index + 1 : 7 * index + 7],476            tuple(self._nuisance_device[3 * index + j : 3 * index + j + 1] for j in range(3)),477        )478479    def _workspaces(480        self,481    ) -> tuple[482        MaterialProjectionWorkspace483        | SpectralWorkspace484        | DetectorWorkspace485        | ObjectiveWorkspace486        | BlurWorkspace,487        ...,488    ]:489        """Use public workspace diagnostics; status encodings stay with operators."""490        return tuple(491            workspace492            for item in self._prepared493            for workspace in (494                item.material,495                item.spectral,496                item.detector,497                item.objective,498                item.blur,499            )500            if workspace is not None501        )502503    def _clear_statuses(self) -> None:504        for workspace in self._workspaces():505            workspace.clear_status()506507    def _check_statuses(self) -> None:508        for workspace in self._workspaces():509            workspace.check_status()510511    def _validate_fixed_inputs(self) -> None:512        for item in self._prepared:513            item.material.validate_inputs(stream=self.context.stream)514            item.spectral.validate_inputs(515                item.view.coefficients,516                item.view.weights,517                item.view.response,518                pixels=item.view.geometry.pixels,519                probability_response=item.view.objective.kind == "poisson",520                stream=self.context.stream,521            )522            item.objective.validate_observation(523                item.view.observation,524                weights=item.view.objective_weights,525                stream=self.context.stream,526            )527528    # region book:spectral-recovery-composition529    def __call__(self, parameters: Vector) -> Evaluation:530        require_no_tape()531        if self._running:532            raise ContractError("a spectral evaluator cannot be used concurrently or recursively")533        if len(parameters) != self.dimension:534            raise ContractError("parameter count differs from pose and calibration charts")535        pose = self.chart.pose(parameters[:6])536        decoded = tuple(537            block.decode(parameters[section])538            for block, section in zip(self.problem.calibration, self._group_slices, strict=True)539        )540        ctx, wp = self.context, self.context.wp541        self._running = True542        self.last_view_losses = None543        try:544            with ctx.scope():545                try:546                    self._clear_statuses()547                    self._pose_values[:] = pose.packed()548                    self._parameter_values[:] = tuple(value for block in decoded for value in block)549                    wp.copy(self._pose_device, self._pose_host, stream=ctx.stream)550                    wp.copy(self._parameters_device, self._parameters_host, stream=ctx.stream)551                    for item in self._prepared:552                        view = item.view553                        gain, exposure, offset = self._parameter_views[item.group]554                        project_material_paths(555                            self._pose_device,556                            workspace=item.material,557                            stream=ctx.stream,558                            validate=False,559                        )560                        spectral_signal(561                            item.material.paths,562                            view.coefficients,563                            view.weights,564                            view.response,565                            out_mean=item.mean,566                            workspace=item.spectral,567                            stream=ctx.stream,568                            validate=False,569                        )570                        if item.blur is not None:571                            blur(572                                item.mean,573                                out_signal=item.spread,574                                workspace=item.blur,575                                stream=ctx.stream,576                                validate=False,577                            )578                        calibrate(579                            item.spread,580                            gain,581                            exposure,582                            offset,583                            out_signal=item.prediction,584                            spec=item.calibration_spec,585                            workspace=item.detector,586                            stream=ctx.stream,587                            validate=False,588                        )589                        evaluate_objective(590                            item.prediction,591                            view.observation,592                            out_loss=item.loss,593                            out_seed=item.image_seed,594                            weights=view.objective_weights,595                            workspace=item.objective,596                            stream=ctx.stream,597                            validate=False,598                        )599                        calibration_vjp(600                            item.spread,601                            gain,602                            exposure,603                            offset,604                            seed=item.image_seed,605                            out_grad_mean=item.spread_seed,606                            out_grad_gain=item.nuisance_gradients[0]607                            if item.calibration_spec.active_gain608                            else None,609                            out_grad_exposure=item.nuisance_gradients[1]610                            if item.calibration_spec.active_exposure611                            else None,612                            out_grad_offset=item.nuisance_gradients[2]613                            if item.calibration_spec.active_offset614                            else None,615                            spec=item.calibration_spec,616                            workspace=item.detector,617                            stream=ctx.stream,618                            validate=False,619                        )620                        if item.blur is not None:621                            blur_transpose(622                                item.spread_seed,623                                out_grad_signal=item.mean_seed,624                                workspace=item.blur,625                                stream=ctx.stream,626                                validate=False,627                            )628                        spectral_vjp(629                            item.material.paths,630                            view.coefficients,631                            view.weights,632                            view.response,633                            seed=item.mean_seed,634                            out_grad_paths=item.material.adj_paths,635                            workspace=item.spectral,636                            stream=ctx.stream,637                            validate=False,638                        )639                        material_projection_vjp(640                            self._pose_device,641                            workspace=item.material,642                            out_pose=item.pose_gradient,643                            stream=ctx.stream,644                            validate=False,645                        )646                    wp.copy(self._control_host, self._control_device, stream=ctx.stream)647                    wp.copy(self._nuisance_host, self._nuisance_device, stream=ctx.stream)648                finally:649                    # Rejected evaluations must finish using pinned staging before650                    # the optimiser can write the next trial into the same buffers.651                    wp.synchronize_stream(ctx.stream)652            self._check_statuses()653            return self._assemble(parameters)654        finally:655            self._running = False656657    # endregion book:spectral-recovery-composition658659    def _assemble(self, parameters: Vector) -> Evaluation:660        losses = tuple(float(self._control_values[7 * i]) for i in range(len(self._prepared)))661        pose_gradient = tuple(662            math.fsum(663                item.view.objective_weight * float(self._control_values[7 * i + 1 + j])664                for i, item in enumerate(self._prepared)665            )666            for j in range(6)667        )668        gradient = list(self.chart.gradient(parameters[:6], pose_gradient))669        prior_losses: list[float] = []670        for group, (block, section) in enumerate(671            zip(self.problem.calibration, self._group_slices, strict=True)672        ):673            physical = tuple(674                math.fsum(675                    item.view.objective_weight * float(self._nuisance_values[3 * i + j])676                    for i, item in enumerate(self._prepared)677                    if item.group == group678                )679                for j in range(3)680            )681            prior = block.prior(parameters[section])682            prior_losses.append(prior.loss)683            derivative = block.gradient(parameters[section], physical)684            gradient.extend(a + b for a, b in zip(derivative, prior.gradient, strict=True))685        loss = math.fsum(686            [687                *(688                    item.view.objective_weight * value689                    for item, value in zip(self._prepared, losses, strict=True)690                ),691                *prior_losses,692            ]693        )694        if not math.isfinite(loss) or not all(math.isfinite(value) for value in gradient):695            raise NumericalError("multi-view loss or chart gradient is non-finite")696        self.last_view_losses = losses697        return Evaluation(loss, tuple(gradient))698699    def calibration_estimates(self, parameters: Vector) -> tuple[CalibrationEstimate, ...]:700        if len(parameters) != self.dimension:701            raise ContractError("parameter count differs from evaluator chart")702        return tuple(703            CalibrationEstimate(block.name, *block.decode(parameters[section]))704            for block, section in zip(self.problem.calibration, self._group_slices, strict=True)705        )706707    @property708    def predicted_images(self) -> tuple[Any, ...]:709        """Borrow the current device predictions; they may belong to a rejected trial.710711        Re-evaluate accepted parameters before explicit export. No host copy is712        made by this accessor, and callers must not mutate these borrowed buffers.713        """714        return tuple(item.prediction for item in self._prepared)715716717def recover_spectral_pose(718    evaluator: SpectralPoseEvaluator,719    *,720    initial: Vector | None = None,721    policy: RecoveryPolicy | None = None,722) -> SpectralRecoveryResult:723    """Use the same safeguarded optimiser and immutable chart as primary recovery."""724    start = (0.0,) * evaluator.dimension if initial is None else initial725    optimisation = recover_parameters(evaluator, start, policy=policy)726    return SpectralRecoveryResult(727        evaluator.chart.pose(optimisation.parameters[:6]),728        evaluator.calibration_estimates(optimisation.parameters),729        optimisation,730    )731