brahmap.core.compute_GLS_maps_from_PTS

Computes the Generalized Least Squares (GLS) maps using a pre-instantiated ProcessTimeSamples instance.

Parameters:

Name	Type	Description	Default
processed_samples	ProcessTimeSamples	The pre-processed time samples object containing pointing and map-making metadata	required
time_ordered_data	NDArray[number]	The 1D vector representing the time-ordered data (TOD) streams	required
inv_noise_cov_operator	DTypeNoiseCov | None	The inverse noise covariance linear operator (\(N^{-1}\)), by default None. If None, the identity matrix will be used as the inverse noise covariance.	None
gls_parameters	GLSParameters	The parameter configuration dictating the map-making behavior, by default GLSParameters()	GLSParameters()
x0	NDArray[number] | None	Initial guess for the GLS solution in the form of interleaved maps (e.g. \([I_1, Q_1, U_1, I_2, Q_2, U_2, \dots]\)), by default None	None

Returns:

Type	Description
GLSResult	The dataclass containing the final output from the GLS map-maker

Source code in brahmap/core/GLS.py

def compute_GLS_maps_from_PTS(
    processed_samples: ProcessTimeSamples,
    time_ordered_data: npt.NDArray[np.number],
    inv_noise_cov_operator: DTypeNoiseCov | None = None,
    gls_parameters: GLSParameters = GLSParameters(),
    x0: npt.NDArray[np.number] | None = None,
) -> GLSResult:
    r"""Computes the Generalized Least Squares (GLS) maps using a
    pre-instantiated `ProcessTimeSamples` instance.

    Parameters
    ----------
    processed_samples : ProcessTimeSamples
        The pre-processed time samples object containing pointing and
        map-making metadata
    time_ordered_data : npt.NDArray[np.number]
        The 1D vector representing the time-ordered data (TOD) streams
    inv_noise_cov_operator : DTypeNoiseCov | None, optional
        The inverse noise covariance linear operator ($N^{-1}$), by
        default `None`. If `None`, the identity matrix will be used as the
        inverse noise covariance.
    gls_parameters : GLSParameters, optional
        The parameter configuration dictating the map-making behavior, by
        default `GLSParameters()`
    x0 : npt.NDArray[np.number] | None, optional
        Initial guess for the GLS solution in the form of interleaved
        maps (e.g. $[I_1, Q_1, U_1, I_2, Q_2, U_2, \dots]$), by default `None`

    Returns
    -------
    GLSResult
        The dataclass containing the final output from the GLS map-maker
    """
    time_ordered_data = np.asarray(time_ordered_data)
    if processed_samples.nsamples != len(time_ordered_data):
        raise ValueError(
            f"Size of `pointings` must be equal to the size of `time_ordered_data` "
            f"array:\nlen(pointings) = {processed_samples.nsamples}\n"
            f"len(time_ordered_data) = {len(time_ordered_data)}"
        )

    try:
        time_ordered_data = time_ordered_data.astype(
            dtype=processed_samples.dtype_float, casting="safe", copy=False
        )
    except TypeError:
        raise TypeError(
            f"The `time_ordered_data` array has higher dtype than "
            f"`processed_samples.dtype_float={processed_samples.dtype_float}`. "
            f"Please compute `processed_samples` again with "
            f"`dtype_float={time_ordered_data.dtype}`"
        )

    if inv_noise_cov_operator is None:
        inv_noise_cov_operator = InvNoiseCovLO_Diagonal(
            size=processed_samples.nsamples, dtype=processed_samples.dtype_float
        )
    else:
        if inv_noise_cov_operator.shape[0] != processed_samples.nsamples:
            raise ValueError(
                f"The shape of `inv_noise_cov_operator` must be same as "
                f"`(len(time_ordered_data), len(time_ordered_data))`:\n"
                f"len(time_ordered_data) = {len(time_ordered_data)}\n"
                f"inv_noise_cov_operator.shape = ({inv_noise_cov_operator.shape}, "
                f"{inv_noise_cov_operator.shape})"
            )

    pointing_operator = PointingLO(
        processed_samples=processed_samples, solver_type=gls_parameters.solver_type
    )

    blockdiagprecond_operator = BlockDiagonalPreconditionerLO(
        processed_samples=processed_samples, solver_type=gls_parameters.solver_type
    )

    b = pointing_operator.T * inv_noise_cov_operator * time_ordered_data

    num_iterations = 0
    if gls_parameters.use_iterative_solver:

        def callback_function(x, r, norm_residual) -> None:
            nonlocal num_iterations
            num_iterations += 1
            if gls_parameters.callback_function is not None:
                gls_parameters.callback_function(x, r, norm_residual)

        A = pointing_operator.T * inv_noise_cov_operator * pointing_operator

        map_vector, pcg_status = cg(
            A=A,  # type: ignore
            b=b,  # type: ignore
            x0=x0,
            atol=gls_parameters.isolver_threshold,
            maxiter=gls_parameters.isolver_max_iterations,
            M=blockdiagprecond_operator,
            callback=callback_function,
            parallel=False,
        )
    else:
        pcg_status = 0
        map_vector = blockdiagprecond_operator * b

    output_maps = separate_map_vectors(
        map_vector=map_vector,  # type: ignore
        processed_samples=processed_samples,
    )

    if gls_parameters.return_hit_map:
        hit_map = processed_samples.get_hit_counts()
    else:
        hit_map = None

    if pcg_status != 0:
        convergence_status = False
    else:
        convergence_status = True

    gls_result = GLSResult(
        solver_type=processed_samples.solver_type,
        npix=processed_samples.npix,
        new_npix=processed_samples.new_npix,
        GLS_maps=output_maps,
        hit_map=hit_map,
        convergence_status=convergence_status,
        num_iterations=num_iterations,
        GLSParameters=gls_parameters,
    )

    return gls_result