diff --git a/brahmap/__init__.py b/brahmap/__init__.py index b8659bc..4d56167 100644 --- a/brahmap/__init__.py +++ b/brahmap/__init__.py @@ -22,6 +22,8 @@ from . import base, _extensions, core, utilities, math # noqa: E402 +from .base import TypeChangeWarning, LowerTypeCastWarning, ShapeError # noqa: E402 + from .core import ( # noqa: E402 SolverType, ProcessTimeSamples, @@ -43,10 +45,7 @@ ) from .utilities import ( # noqa: E402 - TypeChangeWarning, - LowerTypeCastWarning, modify_numpy_context, - ShapeError, ) if find_spec("litebird_sim") is not None: @@ -80,6 +79,9 @@ "MPI_RAISE_EXCEPTION", # ./base/ "base", + "TypeChangeWarning", + "LowerTypeCastWarning", + "ShapeError", # ./_extensions/ "_extensions", # ./core/ @@ -103,10 +105,7 @@ "compute_GLS_maps", # ./utilities/ "utilities", - "TypeChangeWarning", - "LowerTypeCastWarning", "modify_numpy_context", - "ShapeError", # ./math/ "math", ] diff --git a/brahmap/base/__init__.py b/brahmap/base/__init__.py index e7ed3c7..9fab217 100644 --- a/brahmap/base/__init__.py +++ b/brahmap/base/__init__.py @@ -1,3 +1,5 @@ +# Original code: +# # Copyright (c) 2008-2013, Dominique Orban # All rights reserved. # @@ -27,6 +29,18 @@ # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY # OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF # SUCH DAMAGE. +# +# +# Modified version: +# +# Copyright (c) 2023-present, Avinash Anand +# and Giuseppe Puglisi +# +# This file is part of BrahMap. +# +# Licensed under the MIT License. See the file for details. + +from .misc import TypeChangeWarning, LowerTypeCastWarning, filter_warnings, ShapeError from .linop import ( BaseLinearOperator, @@ -59,6 +73,11 @@ ) __all__ = [ + # misc.py + "TypeChangeWarning", + "LowerTypeCastWarning", + "filter_warnings", + "ShapeError", # linop.py "BaseLinearOperator", "LinearOperator", diff --git a/brahmap/base/blkop.py b/brahmap/base/blkop.py index 369e232..6b0ac07 100644 --- a/brahmap/base/blkop.py +++ b/brahmap/base/blkop.py @@ -1,3 +1,5 @@ +# Original code: +# # Copyright (c) 2008-2013, Dominique Orban # All rights reserved. # @@ -27,16 +29,26 @@ # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY # OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF # SUCH DAMAGE. +# +# +# Modified version: +# +# Copyright (c) 2023-present, Avinash Anand +# and Giuseppe Puglisi +# +# This file is part of BrahMap. +# +# Licensed under the MIT License. See the file for details. import numpy as np import itertools import warnings -from typing import List +from typing import List, Any from functools import reduce, partial from ..base import BaseLinearOperator, LinearOperator from ..base import null_log -from ..utilities import ShapeError, TypeChangeWarning +from .misc import ShapeError, TypeChangeWarning from ..mpi import MPI_RAISE_EXCEPTION from brahmap import MPI_UTILS @@ -54,9 +66,22 @@ class BlockLinearOperator(LinearOperator): need be specified, e.g., `[[A,B,C], [D,E], [F]]`, and the blocks on the diagonal must be square and symmetric. + Parameters + ---------- + blocks : List[LinearOperator] + _description_ + symmetric : bool, optional + _description_, by default False + **kwargs: Any + _description_ """ - def __init__(self, blocks, symmetric=False, **kwargs): + def __init__( + self, + blocks: List[LinearOperator], + symmetric: bool = False, + **kwargs: Any, + ): # If building a symmetric operator, fill in the blanks. # They're just references to existing objects. try: @@ -176,7 +201,21 @@ def __iter__(self): class BlockDiagonalLinearOperator(LinearOperator): - def __init__(self, block_list, **kwargs): + """Base class for a block-diagonal linear operator + + Parameters + ---------- + block_list : List[LinearOperator] + _description_ + **kwargs: Any + _description_ + """ + + def __init__( + self, + block_list: List[LinearOperator], + **kwargs: Any, + ): try: for block in block_list: __, __ = block.shape @@ -184,7 +223,7 @@ def __init__(self, block_list, **kwargs): MPI_RAISE_EXCEPTION( condition=True, exception=ValueError, - message="The `block_list` must be a flat list of linear" "operators", + message="The `block_list` must be a flat list of linearoperators", ) self.__row_size = np.asarray( @@ -317,9 +356,19 @@ class BlockHorizontalLinearOperator(BlockLinearOperator): Each block must be a linear operator. The blocks must be specified as one list, e.g., `[A, B, C]`. + Parameters + ---------- + blocks : List[LinearOperator] + _description_ + **kwargs: Any + _description_ """ - def __init__(self, blocks, **kwargs): + def __init__( + self, + blocks: List[LinearOperator], + **kwargs: Any, + ): try: for block in blocks: __ = block.shape @@ -340,9 +389,19 @@ class BlockVerticalLinearOperator(BlockLinearOperator): Each block must be a linear operator. The blocks must be specified as one list, e.g., `[A, B, C]`. + Parameters + ---------- + blocks : List[LinearOperator] + _description_ + **kwargs: Any + _description_ """ - def __init__(self, blocks, **kwargs): + def __init__( + self, + blocks: List[LinearOperator], + **kwargs: Any, + ): try: for block in blocks: __ = block.shape diff --git a/brahmap/base/linop.py b/brahmap/base/linop.py index 24d4809..71c1564 100644 --- a/brahmap/base/linop.py +++ b/brahmap/base/linop.py @@ -1,3 +1,5 @@ +# Original code: +# # Copyright (c) 2008-2013, Dominique Orban # All rights reserved. # @@ -27,14 +29,25 @@ # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY # OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF # SUCH DAMAGE. +# +# +# Modified version: +# +# Copyright (c) 2023-present, Avinash Anand +# and Giuseppe Puglisi +# +# This file is part of BrahMap. +# +# Licensed under the MIT License. See the file for details. +from typing import Callable, Optional, Tuple, Any +import numbers import numpy as np +import numpy.typing as npt import logging -from ..utilities import ShapeError - -__docformat__ = "restructuredtext" +from .misc import ShapeError # Default (null) logger. @@ -44,23 +57,39 @@ class BaseLinearOperator(object): + """Base class for defining the common interface shared by all linear + operators. + + A linear operator is a linear mapping $x \\mapsto A(x)$ such that the size + of the input vector $x$ is `nargin` and the size of the output vector is + `nargout`. The linear operator $A$ can be visualized as a matrix of shape + `(nargout, nargin)`. + + Parameters + ---------- + nargin : int + Size of the input vector $x$ + nargout : int + Size of the output vector $A(x)$ + symmetric : bool, optional + A parameter to specify whether the linear operator is symmetric, by + default `False` + dtype : npt.DTypeLike, optional + Data type of the linear operator, by default `np.float64` + **kwargs : Any + Extra keywords arguments """ - Base class defining the common interface shared by all linear operators. - A linear operator is a linear mapping x -> A(x) such that the size of the - input vector x is `nargin` and the size of the output is `nargout`. It can - be visualized as a matrix of shape (`nargout`, `nargin`). Its type is any - valid Numpy `dtype`. By default, it has `dtype` `numpy.float` but this can - be changed to, e.g., `numpy.complex` via the `dtype` keyword argument and - attribute. - - A logger may be attached to the linear operator via the `logger` keyword - argument. - - """ + # A logger may be attached to the linear operator via the `logger` keyword + # argument. def __init__( - self, nargin, nargout, symmetric: bool = False, dtype=np.float64, **kwargs + self, + nargin: int, + nargout: int, + symmetric: bool = False, + dtype: npt.DTypeLike = np.float64, + **kwargs, ): self.__nargin = nargin self.__nargout = nargout @@ -75,13 +104,13 @@ def __init__( return @property - def nargin(self): - """The size of an input vector.""" + def nargin(self) -> int: + """The size of the input vector.""" return self.__nargin @property - def nargout(self): - """The size of an output vector.""" + def nargout(self) -> int: + """The size of the output vector.""" return self.__nargout @property @@ -90,12 +119,12 @@ def symmetric(self) -> bool: return self.__symmetric @property - def shape(self): + def shape(self) -> Tuple[int, int]: """The shape of the operator.""" return self.__shape @property - def dtype(self): + def dtype(self) -> npt.DTypeLike: """The data type of the operator.""" return self.__dtype @@ -104,7 +133,7 @@ def dtype(self, dtype): self.__dtype = dtype @property - def nMatvec(self): + def nMatvec(self) -> int: """The number of products with vectors computed so far.""" return self._nMatvec @@ -112,6 +141,10 @@ def reset_counters(self): """Reset operator/vector product counter to zero.""" self._nMatvec = 0 + def dot(self, x): + """Numpy-like dot() method.""" + return self.__mul__(x) + def __call__(self, *args, **kwargs): # An alias for __mul__. return self.__mul__(*args, **kwargs) @@ -119,7 +152,7 @@ def __call__(self, *args, **kwargs): def __mul__(self, x): raise NotImplementedError("Please subclass to implement __mul__.") - def __repr__(self): + def __repr__(self) -> str: if self.symmetric: s = "Symmetric" else: @@ -129,23 +162,45 @@ def __repr__(self): s += " with shape (%d,%d)" % (self.nargout, self.nargin) return s - def dot(self, x): - """Numpy-like dot() method.""" - return self.__mul__(x) - class LinearOperator(BaseLinearOperator): """ - Generic linear operator class. - - A linear operator constructed from a `matvec` and (possibly) a - `rmatvec` function. If `symmetric` is `True`, `rmatvec` is - ignored. All other keyword arguments are passed directly to the superclass. - + A generic linear operator class. + + A linear operator constructed from a matrix-vector multiplication `matvec`, + $x \\mapsto A(x)=Ax$ and possibly with a transposed-matrix-vector + operation `rmatvec`, $x \\mapsto A(x)=A^T x$. If `symmetric` is `True`, + `rmatvec` is ignored. All other keyword arguments are passed directly to + the superclass. + + Parameters + ---------- + nargin : int + Size of the input vector $x$ + nargout : int + Size of the output vector $A(x)$ + matvec : Callable + A function that defines the matrix-vector product $x \\mapsto A(x)=Ax$ + rmatvec : Optional[Callable], optional + A function that defines the transposed-matrix-vector product + $x \\mapsto A(x)=A^T x$, by default `None` + **kwargs : Any + Extra keywords arguments """ - def __init__(self, nargin, nargout, matvec, rmatvec=None, **kwargs): - super(LinearOperator, self).__init__(nargin, nargout, **kwargs) + def __init__( + self, + nargin: int, + nargout: int, + matvec: Callable, + rmatvec: Optional[Callable] = None, + **kwargs, + ): + super(LinearOperator, self).__init__( + nargin, + nargout, + **kwargs, + ) adjoint_of = kwargs.get("adjoint_of", None) or kwargs.get("transpose_of", None) rmatvec = rmatvec or kwargs.get("matvec_transp", None) @@ -172,32 +227,30 @@ def __init__(self, nargin, nargout, matvec, rmatvec=None, **kwargs): if isinstance(adjoint_of, BaseLinearOperator): self.__H = adjoint_of else: - msg = "kwarg adjoint_of / transpose_of must be of type LinearOperator." + msg = ( + "kwarg adjoint_of / transpose_of must be of type" + " LinearOperator." + ) msg += " Got " + str(adjoint_of.__class__) raise ValueError(msg) @property def T(self): - """The transpose operator. - - .. note:: this is an alias to the adjoint operator - - """ + """The transpose operator""" return self.__H @property def H(self): - """The adjoint operator.""" + """The adjoint operator""" return self.__H def matvec(self, x): """ Matrix-vector multiplication. - The matvec property encapsulates the matvec routine specified at + The matvec property encapsulates the `matvec` routine specified at construct time, to ensure the consistency of the input and output arrays with the operator's shape. - """ x = np.asanyarray(x) M, N = self.shape @@ -207,7 +260,12 @@ def matvec(self, x): try: x = x.reshape(N) except ValueError: - msg = "input array size incompatible with operator dimensions" + msg = ( + "The size of the input array is incompatible with the " + "operator dimensions\n" + f"size of the input array: {len(x)}\n" + f"shape of the operator: {self.shape}" + ) raise ValueError(msg) y = self.__matvec(x) @@ -217,12 +275,26 @@ def matvec(self, x): try: y = y.reshape(M) except ValueError: - msg = "output array size incompatible with operator dimensions" + msg = ( + "The size of the output array is incompatible with the " + "operator dimensions\n" + f"size of the output array: {len(y)}\n" + f"shape of the operator: {self.shape}" + ) raise ValueError(msg) return y - def to_array(self): + def to_array(self) -> np.ndarray: + """Returns the dense form of the linear operator as a 2D NumPy array + + !!! Warning + + This method first allocates a NumPy array of shape `self.shape` + and data-type `self.dtype`, and then fills them with numbers. As + such it can occupy an enormous amount of memory. Don't use it + unless you understand the risk! + """ n, m = self.shape H = np.empty((n, m), dtype=self.dtype) ej = np.zeros(m, dtype=self.dtype) @@ -233,7 +305,7 @@ def to_array(self): return H def __mul_scalar(self, x): - """Product between a linear operator and a scalar.""" + # Product between a linear operator and a scalar result_type = np.result_type(self.dtype, type(x)) if x != 0: @@ -256,9 +328,14 @@ def rmatvec(y): return ZeroOperator(self.nargin, self.nargout, dtype=result_type) def __mul_linop(self, op): - """Product between two linear operators.""" + # Product between two linear operators if self.nargin != op.nargout: - raise ShapeError("Cannot multiply operators together") + msg = ( + "Cannot multiply the two operators together\n" + f"shape of the first operator: {self.shape}\n" + f"shape of the second operator: {op.shape}" + ) + raise ShapeError(msg) def matvec(x): return self(op(x)) @@ -278,31 +355,38 @@ def rmatvec(x): ) def __mul_vector(self, x): - """Product between a linear operator and a vector.""" + # Product between a linear operator and a vector self._nMatvec += 1 result_type = np.result_type(self.dtype, x.dtype) return self.matvec(x).astype(result_type, copy=False) def __mul__(self, x): - if np.isscalar(x): + # Returns a linear operator if x is a scalar or a linear operator + # Returns a vector if x is an array + if isinstance(x, numbers.Number): return self.__mul_scalar(x) elif isinstance(x, BaseLinearOperator): return self.__mul_linop(x) elif isinstance(x, np.ndarray): return self.__mul_vector(x) else: - raise ValueError("Cannot multiply") + raise ValueError("Invalid multiplier! Cannot multiply") def __rmul__(self, x): if np.isscalar(x): return self.__mul__(x) - raise ValueError("Cannot multiply") + raise ValueError("Invalid operation! Cannot multiply") def __add__(self, other): if not isinstance(other, BaseLinearOperator): - raise ValueError("Cannot add") + raise ValueError("Invalid operation! Cannot add") if self.shape != other.shape: - raise ShapeError("Cannot add") + msg = ( + "Cannot add the two operators together\n" + f"shape of the first operator: {self.shape}\n" + f"shape of the second operator: {other.shape}" + ) + raise ShapeError(msg) def matvec(x): return self(x) + other(x) @@ -326,9 +410,14 @@ def __neg__(self): def __sub__(self, other): if not isinstance(other, BaseLinearOperator): - raise ValueError("Cannot add") + raise ValueError("Invalid operation! Cannot subtract") if self.shape != other.shape: - raise ShapeError("Cannot add") + msg = ( + "Cannot subtract one operator from the other\n" + f"shape of the first operator: {self.shape}\n" + f"shape of the second operator: {other.shape}" + ) + raise ShapeError(msg) def matvec(x): return self(x) - other(x) @@ -351,7 +440,7 @@ def __truediv__(self, other): if np.isscalar(other): return self * (1 / other) else: - raise ValueError("Cannot divide") + raise ValueError("Invalid operation! Cannot divide") def __pow__(self, other): if not isinstance(other, int): @@ -368,9 +457,17 @@ def __pow__(self, other): class IdentityOperator(LinearOperator): - """Class representing the identity operator of size `nargin`.""" + """A linear operator for the identity matrix of size `nargin` + + Parameters + ---------- + nargin : int + _description_ + **kwargs: Any + _description_ + """ - def __init__(self, nargin, **kwargs): + def __init__(self, nargin: int, **kwargs: Any): if "symmetric" in kwargs: kwargs.pop("symmetric") if "matvec" in kwargs: @@ -382,16 +479,17 @@ def __init__(self, nargin, **kwargs): class DiagonalOperator(LinearOperator): + """A linear operator for a diagonal matrix + + Parameters + ---------- + diag : np.ndarray + _description_ + **kwargs: Any + _description_ """ - Class representing a diagonal operator. - A diagonal linear operator defined by its diagonal `diag` (a Numpy array.) - The type must be specified in the `diag` argument, e.g., - `np.ones(5, dtype=np.complex)` or `np.ones(5).astype(np.complex)`. - - """ - - def __init__(self, diag, **kwargs): + def __init__(self, diag: np.ndarray, **kwargs: Any): if "symmetric" in kwargs: kwargs.pop("symmetric") if "matvec" in kwargs: @@ -415,18 +513,21 @@ def __init__(self, diag, **kwargs): class MatrixLinearOperator(LinearOperator): - """ - Class representing a matrix operator. + """A linear operator for a numpy matrix A linear operator wrapping the multiplication with a matrix and its transpose (real) or conjugate transpose (complex). The operator's dtype is the same as the specified `matrix` argument. - .. versionadded:: 0.3 - + Parameters + ---------- + matrix : np.ndarray + _description_ + **kwargs: Any + _description_ """ - def __init__(self, matrix, **kwargs): + def __init__(self, matrix: np.ndarray, **kwargs: Any): if "symmetric" in kwargs: kwargs.pop("symmetric") if "matvec" in kwargs: @@ -466,9 +567,19 @@ def __init__(self, matrix, **kwargs): class ZeroOperator(LinearOperator): - """Class representing the zero operator of shape `nargout`-by-`nargin`.""" + """A linear operator for a zero matrix of shape `(nargout, nargin)` + + Parameters + ---------- + nargin : int + _description_ + nargout : int + _description_ + **kwargs: Any + _description_ + """ - def __init__(self, nargin, nargout, **kwargs): + def __init__(self, nargin: int, nargout: int, **kwargs: Any): if "matvec" in kwargs: kwargs.pop("matvec") if "rmatvec" in kwargs: @@ -495,19 +606,28 @@ def rmatvec(x): class InverseLO(LinearOperator): r""" - Construct the inverse operator of a matrix :math:`A`, as a linear operator. + Construct the inverse operator of a matrix `A`, as a linear operator. - **Parameters** - - - ``A`` : {linear operator} - the linear operator of the linear system to invert; - - ``method`` : {function } - the method to compute ``A^-1`` (see below); - - ``P`` : {linear operator } (optional) - the preconditioner for the computation of the inverse operator. + Parameters + ---------- + A : _type_ + _description_ + method : _type_, optional + _description_, by default None + preconditioner : _type_, optional + _description_, by default None """ + def __init__(self, A, method=None, preconditioner=None): + super(InverseLO, self).__init__( + nargin=A.shape[0], nargout=A.shape[1], matvec=self.mult, symmetric=True + ) + self.A = A + self.__method = method + self.__preconditioner = preconditioner + self.__converged = None + def mult(self, x): r""" It returns :math:`y=A^{-1}x` by solving the linear system :math:`Ay=x` @@ -535,15 +655,6 @@ def isconverged(self, info): else: return False - def __init__(self, A, method=None, preconditioner=None): - super(InverseLO, self).__init__( - nargin=A.shape[0], nargout=A.shape[1], matvec=self.mult, symmetric=True - ) - self.A = A - self.__method = method - self.__preconditioner = preconditioner - self.__converged = None - @property def method(self): r""" @@ -576,9 +687,9 @@ def preconditioner(self): def ReducedLinearOperator(op, row_indices, col_indices): """ - Implement reduction of a linear operator (non symmetrical). + Implements reduction of a linear operator (non symmetrical). - Reduce a linear operator by limiting its input to `col_indices` and its + Reduces a linear operator by limiting its input to `col_indices` and its output to `row_indices`. """ @@ -605,9 +716,9 @@ def rmatvec(x): def SymmetricallyReducedLinearOperator(op, indices): """ - Implement reduction of a linear operator (symmetrical). + Implements reduction of a linear operator (symmetrical). - Reduce a linear operator symmetrically by reducing boths its input and + Reduces a linear operator symmetrically by reducing boths its input and output to `indices`. """ @@ -633,7 +744,7 @@ def rmatvec(x): def aslinearoperator(A): - """Return A as a LinearOperator. + """Returns A as a LinearOperator. 'A' may be any of the following types: - linop.LinearOperator @@ -643,10 +754,7 @@ def aslinearoperator(A): - sparse matrix (e.g. csr_matrix, lil_matrix, etc.) - any object with .shape and .matvec attributes - See the :class:`LinearOperator` documentation for additonal information. - - .. versionadded:: 0.4 - + See the `LinearOperator` documentation for additonal information. """ if isinstance(A, LinearOperator): return A diff --git a/brahmap/base/misc.py b/brahmap/base/misc.py new file mode 100644 index 0000000..140b33d --- /dev/null +++ b/brahmap/base/misc.py @@ -0,0 +1,44 @@ +import warnings + + +class TypeChangeWarning(Warning): + def __init__(self, message): + self.message = message + + def __str__(self): + return repr(self.message) + + +class LowerTypeCastWarning(Warning): + def __init__(self, message): + self.message = message + + def __str__(self): + return repr(self.message) + + +def filter_warnings(wfilter): + """ + wfilter: {string} + - "ignore": never print matching warnings; + - "always": always print matching warnings + + """ + warnings.simplefilter(wfilter) + + +class ShapeError(Exception): + """ + Exception class for handling shape mismatch errors. + + Exception raised when defining a linear operator of the wrong shape or + multiplying a linear operator with a vector of the wrong shape. + + """ + + def __init__(self, value): + super(ShapeError, self).__init__() + self.value = value + + def __str__(self): + return repr(self.value) diff --git a/brahmap/base/noise_ops.py b/brahmap/base/noise_ops.py index 61e9617..800c19e 100644 --- a/brahmap/base/noise_ops.py +++ b/brahmap/base/noise_ops.py @@ -1,5 +1,5 @@ import numpy as np -from typing import Literal +from typing import Literal, List, Any from ..base import LinearOperator, BlockDiagonalLinearOperator @@ -9,13 +9,29 @@ class NoiseCovLinearOperator(LinearOperator): + """Base class for noise covariance operators + + Parameters + ---------- + nargin : int + _description_ + matvec : int + _description_ + input_type : Literal["covariance", "power_spectrum"], optional + _description_, by default "covariance" + dtype : DTypeFloat, optional + _description_, by default np.float64 + **kwargs: Any + _description_ + """ + def __init__( self, nargin: int, matvec: int, input_type: Literal["covariance", "power_spectrum"] = "covariance", dtype: DTypeFloat = np.float64, - **kwargs, + **kwargs: Any, ): MPI_RAISE_EXCEPTION( condition=(input_type not in ["covariance", "power_spectrum"]), @@ -51,13 +67,29 @@ def get_inverse(self) -> "InvNoiseCovLinearOperator": class InvNoiseCovLinearOperator(NoiseCovLinearOperator): + """Base class for inverse noise covariance operators + + Parameters + ---------- + nargin : int + _description_ + matvec : int + _description_ + input_type : Literal["covariance", "power_spectrum"], optional + _description_, by default "covariance" + dtype : DTypeFloat, optional + _description_, by default np.float64 + **kwargs: Any + _description_ + """ + def __init__( self, nargin: int, matvec: int, input_type: Literal["covariance", "power_spectrum"] = "covariance", dtype: DTypeFloat = np.float64, - **kwargs, + **kwargs: Any, ): super(InvNoiseCovLinearOperator, self).__init__( nargin, @@ -69,7 +101,21 @@ def __init__( class BaseBlockDiagNoiseCovLinearOperator(BlockDiagonalLinearOperator): - def __init__(self, block_list, **kwargs): + """Base class for block-diagonal noise covariance operator + + Parameters + ---------- + block_list : List[NoiseCovLinearOperator] + _description_ + **kwargs: Any + _description_ + """ + + def __init__( + self, + block_list: List[NoiseCovLinearOperator], + **kwargs: Any, + ): super(BaseBlockDiagNoiseCovLinearOperator, self).__init__(block_list, **kwargs) MPI_RAISE_EXCEPTION( @@ -93,7 +139,21 @@ def get_inverse(self) -> "BaseBlockDiagInvNoiseCovLinearOperator": class BaseBlockDiagInvNoiseCovLinearOperator(BaseBlockDiagNoiseCovLinearOperator): - def __init__(self, block_list, **kwargs): + """Base class for block-diagonal inverse noise covariance operator + + Parameters + ---------- + block_list : List[InvNoiseCovLinearOperator] + _description_ + **kwargs: Any + _description_ + """ + + def __init__( + self, + block_list: List[InvNoiseCovLinearOperator], + **kwargs: Any, + ): super(BaseBlockDiagInvNoiseCovLinearOperator, self).__init__( block_list, **kwargs ) diff --git a/brahmap/core/GLS.py b/brahmap/core/GLS.py index 479666e..68e7973 100644 --- a/brahmap/core/GLS.py +++ b/brahmap/core/GLS.py @@ -21,6 +21,26 @@ @dataclass class GLSParameters: + """A class to encapsulate the parameters used for GLS map-making + + Attributes + ---------- + solver_type : SolverType + _description_ + use_iterative_solver : bool + _description_ + isolver_threshold : float + _description_ + isolver_max_iterations : int + _description_ + callback_function : Callable + _description_ + return_processed_samples : bool + _description_ + return_hit_map : bool + _description_ + """ + solver_type: SolverType = SolverType.IQU use_iterative_solver: bool = True isolver_threshold: float = 1.0e-12 @@ -32,6 +52,28 @@ class GLSParameters: @dataclass class GLSResult: + """A class to store the results of the GLS map-making + + Parameters + ---------- + solver_type : SolverType + _description_ + npix : int + _description_ + new_npix : int + _description_ + GLS_maps : np.ndarray + _description_ + hit_map : np.ndarray + _description_ + convergence_status : bool + _description_ + num_iterations : int + _description_ + GLSParameters : GLSParameters + _description_ + """ + solver_type: SolverType npix: int new_npix: int @@ -45,6 +87,24 @@ class GLSResult: def separate_map_vectors( map_vector: np.ndarray, processed_samples: ProcessTimeSamples ) -> np.ndarray: + """The output maps of the GLS are in the form + [I_1, Q_1, U_1, I_2, Q_2, U_2, ...]. Following the typical conventions, + the Stokes parameters have to be separated as [I_1, I_2, ...], + [Q_1, Q_2, ...] and [U_1, U_2, ...]. This function performs this operation + ane returns the maps of different Stokes parameters separately. + + Parameters + ---------- + map_vector : np.ndarray + _description_ + processed_samples : ProcessTimeSamples + _description_ + + Returns + ------- + np.ndarray + _description_ + """ try: map_vector = np.reshape( map_vector, @@ -79,6 +139,25 @@ def compute_GLS_maps_from_PTS( inv_noise_cov_operator: Union[DTypeNoiseCov, None] = None, gls_parameters: GLSParameters = GLSParameters(), ) -> GLSResult: + """This function computes the GLS maps given an instance of + `ProcessTimeSamples`, TOD, and inverse noise covariance operator + + Parameters + ---------- + processed_samples : ProcessTimeSamples + _description_ + time_ordered_data : np.ndarray + _description_ + inv_noise_cov_operator : Union[DTypeNoiseCov, None], optional + _description_, by default None + gls_parameters : GLSParameters, optional + _description_, by default GLSParameters() + + Returns + ------- + GLSResult + _description_ + """ MPI_RAISE_EXCEPTION( condition=(processed_samples.nsamples != len(time_ordered_data)), exception=ValueError, @@ -178,6 +257,36 @@ def compute_GLS_maps( update_pointings_inplace: bool = True, gls_parameters: GLSParameters = GLSParameters(), ) -> Union[GLSResult, tuple[ProcessTimeSamples, GLSResult]]: + """The function to compute the GLS maps given pointing information and TOD + + Parameters + ---------- + npix : int + _description_ + pointings : np.ndarray + _description_ + time_ordered_data : np.ndarray + _description_ + pointings_flag : Union[np.ndarray, None], optional + _description_, by default None + pol_angles : Union[np.ndarray, None], optional + _description_, by default None + inv_noise_cov_operator : Union[DTypeNoiseCov, None], optional + _description_, by default None + threshold : float, optional + _description_, by default 1.0e-5 + dtype_float : Union[DTypeFloat, None], optional + _description_, by default None + update_pointings_inplace : bool, optional + _description_, by default True + gls_parameters : GLSParameters, optional + _description_, by default GLSParameters() + + Returns + ------- + Union[GLSResult, tuple[ProcessTimeSamples, GLSResult]] + _description_ + """ if dtype_float is None: if pol_angles is None: dtype_float = time_ordered_data.dtype diff --git a/brahmap/core/linearoperators.py b/brahmap/core/linearoperators.py index 85b0a88..71261db 100644 --- a/brahmap/core/linearoperators.py +++ b/brahmap/core/linearoperators.py @@ -6,7 +6,7 @@ from ..core import SolverType, ProcessTimeSamples -from ..utilities import TypeChangeWarning +from ..base import TypeChangeWarning from .._extensions import PointingLO_tools from .._extensions import BlkDiagPrecondLO_tools @@ -17,21 +17,16 @@ class PointingLO(LinearOperator): - r"""Derived class from the one from the :class:`LinearOperator` in :mod:`linop`. + """Derived class from the one from the :class:`LinearOperator` in :mod:`linop`. It constitutes an interface for dealing with the projection operator (pointing matrix). - Since this can be represented as a sparse matrix, it is initialized \ - by an array of observed pixels which resembles the ``(i,j)`` positions \ - of the non-null elements of the matrix,``obs_pixs``. - - **Parameters** - - - ``processed_samples``: {:class:`ProcessTimeSamples`} - the class (instantiated before :class:`PointingLO`)has already computed - the :math:`\cos 2\phi` and :math:`\sin 2\phi`, we link the ``cos2phi`` and ``sin2phi`` - attributes of this class to the :class:`ProcessTimeSamples` ones ; - + Parameters + ---------- + processed_samples : ProcessTimeSamples + _description_ + solver_type : Union[None, SolverType], optional + _description_, by default None """ def __init__( @@ -322,26 +317,18 @@ class BlockDiagonalPreconditionerLO(LinearOperator): r""" Standard preconditioner defined as: - .. math:: - - M_{BD}=( A diag(N^{-1}) A^T)^{-1} + $$M_{BD}=( P^T diag(N^{-1}) P)^{-1}$$ - where :math:`A` is the *pointing matrix* (see :class:`PointingLO`). + where $P$ is the *pointing matrix* (see `PointingLO`). Such inverse operator could be easily computed given the structure of the - matrix :math:`A`. It could be sparse in the case of Intensity only analysis (`pol=1`), - block-sparse if polarization is included (`pol=3,2`). - - - **Parameters** - - - ``n``:{int} - the size of the problem, ``npix``; - - ``CES``:{:class:`ProcessTimeSamples`} - the linear operator related to the data sample processing. Its members (`counts`, `masks`, - `sine`, `cosine`, etc... ) are needed to explicitly compute the inverse of the - :math:`n_{pix}` blocks of :math:`M_{BD}`. - - ``pol``:{int} - the size of each block of the matrix. + matrix $P$. + + Parameters + ---------- + processed_samples : ProcessTimeSamples + _description_ + solver_type : Union[None, SolverType], optional + _description_, by default None """ def __init__( diff --git a/brahmap/core/noise_ops_block_diag.py b/brahmap/core/noise_ops_block_diag.py index 0784b15..c354c0a 100644 --- a/brahmap/core/noise_ops_block_diag.py +++ b/brahmap/core/noise_ops_block_diag.py @@ -10,6 +10,24 @@ class BlockDiagNoiseCovLO(BaseBlockDiagNoiseCovLinearOperator): + """Linear operator for block-diagonal noise covariance + + Parameters + ---------- + operator : _type_ + _description_ + block_size : Union[np.ndarray, List] + _description_ + block_input : List[Union[np.ndarray, List]] + _description_ + input_type : Literal["covariance", "power_spectrum"], optional + _description_, by default "power_spectrum" + dtype : DTypeFloat, optional + _description_, by default np.float64 + extra_kwargs : Dict[str, Any], optional + _description_, by default {} + """ + def __init__( self, operator, @@ -62,6 +80,24 @@ def __build_blocks( class BlockDiagInvNoiseCovLO(BlockDiagNoiseCovLO): + """Linear operator for block-diagonal inverse noise covariance + + Parameters + ---------- + operator : _type_ + _description_ + block_size : Union[np.ndarray, List] + _description_ + block_input : List[Union[np.ndarray, List]] + _description_ + input_type : Literal["covariance", "power_spectrum"], optional + _description_, by default "power_spectrum" + dtype : DTypeFloat, optional + _description_, by default np.float64 + extra_kwargs : Dict[str, Any], optional + _description_, by default {} + """ + def __init__( self, operator, diff --git a/brahmap/core/noise_ops_circulant.py b/brahmap/core/noise_ops_circulant.py index e37644d..3a45901 100644 --- a/brahmap/core/noise_ops_circulant.py +++ b/brahmap/core/noise_ops_circulant.py @@ -2,7 +2,7 @@ import warnings from typing import List, Union, Literal -from ..utilities import TypeChangeWarning +from ..base import TypeChangeWarning from ..base import NoiseCovLinearOperator, InvNoiseCovLinearOperator from ..math import DTypeFloat from ..mpi import MPI_RAISE_EXCEPTION @@ -11,6 +11,20 @@ class NoiseCovLO_Circulant(NoiseCovLinearOperator): + """Linear operator for Circulant noise covariance + + Parameters + ---------- + size : int + _description_ + input : Union[np.ndarray, List] + _description_ + input_type : Literal["covariance", "power_spectrum"], optional + _description_, by default "power_spectrum" + dtype : DTypeFloat, optional + _description_, by default np.float64 + """ + def __init__( self, size: int, @@ -80,6 +94,20 @@ def _mult(self, vec: np.ndarray): class InvNoiseCovLO_Circulant(InvNoiseCovLinearOperator): + """Linear operator for the inverse of Circulant noise covariance + + Parameters + ---------- + size : int + _description_ + input : Union[np.ndarray, List] + _description_ + input_type : Literal["covariance", "power_spectrum"], optional + _description_, by default "power_spectrum" + dtype : DTypeFloat, optional + _description_, by default np.float64 + """ + def __init__( self, size: int, diff --git a/brahmap/core/noise_ops_diagonal.py b/brahmap/core/noise_ops_diagonal.py index a25c843..7eb3ec4 100644 --- a/brahmap/core/noise_ops_diagonal.py +++ b/brahmap/core/noise_ops_diagonal.py @@ -4,7 +4,7 @@ from typing import List, Union, Literal -from ..utilities import TypeChangeWarning +from ..base import TypeChangeWarning from ..math import DTypeFloat, linalg_tools @@ -16,6 +16,20 @@ class NoiseCovLO_Diagonal(NoiseCovLinearOperator): + """Linear operator for diagonal noise covariance + + Parameters + ---------- + size : int + _description_ + input : Union[np.ndarray, List, DTypeFloat], optional + _description_, by default 1.0 + input_type : Literal["covariance", "power_spectrum"], optional + _description_, by default "covariance" + dtype : DTypeFloat, optional + _description_, by default np.float64 + """ + def __init__( self, size: int, @@ -91,6 +105,20 @@ def _mult(self, vec: np.ndarray): class InvNoiseCovLO_Diagonal(InvNoiseCovLinearOperator): + """Linear operator for the inverse of diagonal noise covariance + + Parameters + ---------- + size : int + _description_ + input : Union[np.ndarray, List, DTypeFloat], optional + _description_, by default 1.0 + input_type : Literal["covariance", "power_spectrum"], optional + _description_, by default "covariance" + dtype : DTypeFloat, optional + _description_, by default np.float64 + """ + def __init__( self, size: int, diff --git a/brahmap/core/noise_ops_toeplitz.py b/brahmap/core/noise_ops_toeplitz.py index a533ef0..43d7124 100644 --- a/brahmap/core/noise_ops_toeplitz.py +++ b/brahmap/core/noise_ops_toeplitz.py @@ -1,8 +1,8 @@ import numpy as np import warnings -from typing import List, Union, Literal +from typing import List, Union, Literal, Callable -from ..utilities import TypeChangeWarning +from ..base import TypeChangeWarning from ..base import LinearOperator, NoiseCovLinearOperator, InvNoiseCovLinearOperator from ..math import DTypeFloat, cg from ..mpi import MPI_RAISE_EXCEPTION @@ -12,7 +12,22 @@ class NoiseCovLO_Toeplitz01(NoiseCovLinearOperator): - """The input covariance array must be at least of the size n. The input power spectrum array must be of the size 2n-2 or 2n-1.""" + """Linear operator for Toeplitz noise covariance + + The input covariance array must be at least of the size n. The input power + spectrum array must be of the size 2n-2 or 2n-1. + + Parameters + ---------- + size : int + _description_ + input : Union[np.ndarray, List] + _description_ + input_type : Literal["covariance", "power_spectrum"], optional + _description_, by default "power_spectrum" + dtype : DTypeFloat, optional + _description_, by default np.float64 + """ def __init__( self, @@ -101,7 +116,32 @@ def _mult(self, vec: np.ndarray): class InvNoiseCovLO_Toeplitz01(InvNoiseCovLinearOperator): - """The input covariance array must be at least of the size n. The input power spectrum array must be of the size 2n-2 or 2n-1.""" + """Linear operator for the inverse of Toeplitz noise covariance + + The input covariance array must be at least of the size n. The input power + spectrum array must be of the size 2n-2 or 2n-1. + + Parameters + ---------- + size : int + _description_ + input : Union[np.ndarray, List] + _description_ + input_type : Literal["covariance", "power_spectrum"], optional + _description_, by default "power_spectrum" + precond_op : Union[ LinearOperator, Literal[None, "Strang", "TChan", "RChan", "KK2"] ], optional + _description_, by default None + precond_maxiter : int, optional + _description_, by default 50 + precond_rtol : float, optional + _description_, by default 1.0e-10 + precond_atol : float, optional + _description_, by default 1.0e-10 + precond_callback : Callable, optional + _description_, by default None + dtype : DTypeFloat, optional + _description_, by default np.float64 + """ def __init__( self, @@ -111,10 +151,10 @@ def __init__( precond_op: Union[ LinearOperator, Literal[None, "Strang", "TChan", "RChan", "KK2"] ] = None, - precond_maxiter=50, - precond_rtol=1.0e-10, - precond_atol=1.0e-10, - precond_callback=None, + precond_maxiter: int = 50, + precond_rtol: float = 1.0e-10, + precond_atol: float = 1.0e-10, + precond_callback: Callable = None, dtype: DTypeFloat = np.float64, ): self.__toeplitz_op = NoiseCovLO_Toeplitz01( diff --git a/brahmap/core/process_time_samples.py b/brahmap/core/process_time_samples.py index 7a60f1c..1ccc285 100644 --- a/brahmap/core/process_time_samples.py +++ b/brahmap/core/process_time_samples.py @@ -1,6 +1,6 @@ from enum import IntEnum import numpy as np -from typing import Union +from typing import Optional from mpi4py import MPI @@ -17,6 +17,8 @@ class SolverType(IntEnum): + """Map-making level: I, QU, or IQU""" + I = 1 # noqa: E741 QU = 2 IQU = 3 @@ -24,28 +26,35 @@ class SolverType(IntEnum): class ProcessTimeSamples(object): """ - A class to store the pre-processed and pre-computed arrays that can be used later. + A class to store the pre-processed and pre-computed arrays that can be + used later. Parameters ---------- npix : int - Number of pixels on which the map-making has to be done. Equal to `healpy.nside2npix(nside)` for a healpix map of given `nside` + Number of pixels on which the map-making has to be done. Equal to + `healpy.nside2npix(nside)` for a healpix map of given `nside` pointings : np.ndarray A 1-d array of pointing indices pointings_flag : np.ndarray - A 1-d array of pointing flags. `True` means good pointing, `False` means bad pointing. + A 1-d array of pointing flags. `True` means good pointing, `False` + means bad pointing. solver_type : SolverType Map-making level: I or QU or IQU pol_angles : np.ndarray | None A 1-d array containing the orientation angles of the detectors noise_weights : np.ndarray | None - A 1-d array of noise weights, or the diagonal elements of the inverse of noise covariance matrix + A 1-d array of noise weights, or the diagonal elements of the inverse + of noise covariance matrix threshold : float The threshold to be used to flag pixels in the sky dtype_float : boh `dtype` of the floating point arrays update_pointings_inplace : bool - The class does some operations on the pointings array. Do you want to make these operations happen in-place? If yes, you will save a lot of memory. Not recommended if you are willing to use pointing arrays somewhere after doing map-making. + The class does some operations on the pointings array. Do you want to + make these operations happen in-place? If yes, you will save a lot of + memory. Not recommended if you are willing to use pointing arrays + somewhere after doing map-making. Attributes ---------- @@ -70,15 +79,17 @@ class ProcessTimeSamples(object): observed_pixels : np.ndarray Pixel indices that are considered for map-making pixel_flag : np.ndarray - A 1-d array of size `npix`. `True` indicates that the corresponding pixel index will be dropped in map-making + A 1-d array of size `npix`. `True` indicates that the corresponding + pixel index will be dropped in map-making bad_pixels : np.ndarray - A 1-d array that contains all the pixel indices that will be excluded in map-making + A 1-d array that contains all the pixel indices that will be excluded + in map-making weighted_counts : np.ndarray Weighted counts sin2phi : np.ndarray - A 1-d array of `sin(2\phi)` + A 1-d array of $sin(2\\phi)$ cos2phi : np.ndarray - A 1-d array of `cos(2\phi)` + A 1-d array of $cos(2\\phi)$ weighted_sin : np.ndarray Weighted `sin` weighted_cos : np.ndarray @@ -92,7 +103,8 @@ class ProcessTimeSamples(object): one_over_determinant : np.ndarray Inverse of determinant for each valid pixels new_npix : int - The number of pixels actually being used in map-making. Equal to `len(observed_pixels)` + The number of pixels actually being used in map-making. Equal to + `len(observed_pixels)` """ @@ -100,12 +112,12 @@ def __init__( self, npix: int, pointings: np.ndarray, - pointings_flag: Union[np.ndarray, None] = None, + pointings_flag: Optional[np.ndarray] = None, solver_type: SolverType = SolverType.IQU, - pol_angles: Union[np.ndarray, None] = None, - noise_weights: Union[np.ndarray, None] = None, + pol_angles: Optional[np.ndarray] = None, + noise_weights: Optional[np.ndarray] = None, threshold: float = 1.0e-5, - dtype_float: Union[DTypeFloat, None] = None, + dtype_float: Optional[DTypeFloat] = None, update_pointings_inplace: bool = False, ): self.__npix = npix @@ -127,7 +139,10 @@ def __init__( MPI_RAISE_EXCEPTION( condition=(len(self.pointings_flag) != self.nsamples), exception=AssertionError, - message=f"Size of `pointings_flag` must be equal to the size of `pointings` array:\nlen(pointings_flag) = {len(self.pointings_flag)}\nlen(pointings) = {self.nsamples}", + message="Size of `pointings_flag` must be equal to the size of " + "`pointings` array:\n" + f"len(pointings_flag) = {len(self.pointings_flag)}\n" + f"len(pointings) = {self.nsamples}", ) self.__solver_type = solver_type @@ -136,15 +151,24 @@ def __init__( MPI_RAISE_EXCEPTION( condition=(self.solver_type not in [1, 2, 3]), exception=ValueError, - message="Invalid `solver_type`!!!\n`solver_type` must be either SolverType.I, SolverType.QU or SolverType.IQU (equivalently 1, 2 or 3).", + message="Invalid `solver_type`!!!\n`solver_type` must be either " + "SolverType.I, SolverType.QU or SolverType.IQU " + "(equivalently 1, 2 or 3).", ) - # setting the dtype for the `float` arrays: if one or both of `noise_weights` and `pol_angles` are supplied, the `dtype_float` will be inferred from them. Otherwise, the it will be set to `np.float64` + # setting the dtype for the `float` arrays: if one or both of + # `noise_weights` and `pol_angles` are supplied, the `dtype_float` + # will be inferred from them. Otherwise, the it will be set to + # `np.float64` if dtype_float is not None: self.__dtype_float = dtype_float elif noise_weights is not None and pol_angles is not None: - # if both `noise_weights` and `pol_angles` are given, `dtype_float` will be assigned the higher `dtype` - self.__dtype_float = np.promote_types(noise_weights.dtype, pol_angles.dtype) + # if both `noise_weights` and `pol_angles` are given, + # `dtype_float` will be assigned the higher `dtype` + self.__dtype_float = np.promote_types( + noise_weights.dtype, + pol_angles.dtype, + ) elif noise_weights is not None: self.__dtype_float = noise_weights.dtype elif pol_angles is not None: @@ -158,7 +182,10 @@ def __init__( MPI_RAISE_EXCEPTION( condition=(len(noise_weights) != self.nsamples), exception=AssertionError, - message=f"Size of `noise_weights` must be equal to the size of `pointings` array:\nlen(noise_weigths) = {len(noise_weights)}\nlen(pointings) = {self.nsamples}", + message="Size of `noise_weights` must be equal to the size of " + "`pointings` array:\n" + f"len(noise_weigths) = {len(noise_weights)}\n" + f"len(pointings) = {self.nsamples}", ) try: @@ -167,14 +194,19 @@ def __init__( ) except TypeError: raise TypeError( - f"The `noise_weights` array has higher dtype than `self.dtype_float={self.dtype_float}`. Please called `ProcessTimeSamples` again with `dtype_float={noise_weights.dtype}`" + "The `noise_weights` array has higher dtype than " + f"`self.dtype_float={self.dtype_float}`. Please call " + f"`ProcessTimeSamples` again with `dtype_float={noise_weights.dtype}`" ) if self.solver_type != 1: MPI_RAISE_EXCEPTION( condition=(len(pol_angles) != self.nsamples), exception=AssertionError, - message=f"Size of `pol_angles` must be equal to the size of `pointings` array:\nlen(pol_angles) = {len(pol_angles)}\nlen(pointings) = {self.nsamples}", + message="Size of `pol_angles` must be equal to the size of " + "`pointings` array:\n" + f"len(pol_angles) = {len(pol_angles)}\n" + f"len(pointings) = {self.nsamples}", ) try: @@ -183,7 +215,9 @@ def __init__( ) except TypeError: raise TypeError( - f"The `pol_angles` array has higher dtype than `self.dtype_float={self.dtype_float}`. Please called `ProcessTimeSamples` again with `dtype_float={pol_angles.dtype}`" + "The `pol_angles` array has higher dtype than " + f"`self.dtype_float={self.dtype_float}`. Please call " + f"`ProcessTimeSamples` again with `dtype_float={pol_angles.dtype}`" ) self._compute_weights( @@ -194,7 +228,8 @@ def __init__( MPI_RAISE_EXCEPTION( condition=(self.new_npix == 0), exception=ValueError, - message="All pixels were found to be pathological. The map-making cannot be done. Please ensure that the inputs are consistent!", + message="All pixels were found to be pathological. The map-making " + "cannot be done. Please ensure that the inputs are consistent!", ) self._repixelization() @@ -269,20 +304,21 @@ def bad_pixels(self): def get_hit_counts(self): """Returns hit counts of the pixel indices""" - hit_counts_newidx = np.zeros(self.new_npix, dtype=int) - for idx in range(self.nsamples): - hit_counts_newidx[self.pointings[idx]] += self.pointings_flag[idx] + return self.hit_counts + # hit_counts_newidx = np.zeros(self.new_npix, dtype=int) + # for idx in range(self.nsamples): + # hit_counts_newidx[self.pointings[idx]] += self.pointings_flag[idx] - MPI_UTILS.comm.Allreduce(MPI.IN_PLACE, hit_counts_newidx, MPI.SUM) + # MPI_UTILS.comm.Allreduce(MPI.IN_PLACE, hit_counts_newidx, MPI.SUM) - hit_counts = np.ma.masked_array( - data=np.zeros(self.npix), - mask=np.logical_not(self.pixel_flag), - fill_value=-1.6375e30, - ) + # hit_counts = np.ma.masked_array( + # data=np.zeros(self.npix), + # mask=np.logical_not(self.pixel_flag), + # fill_value=-1.6375e30, + # ) - hit_counts[~hit_counts.mask] = hit_counts_newidx - return hit_counts + # hit_counts[~hit_counts.mask] = hit_counts_newidx + # return hit_counts def _compute_weights(self, pol_angles: np.ndarray, noise_weights: np.ndarray): self.hit_counts = np.zeros(self.npix, dtype=self.pointings.dtype) diff --git a/brahmap/lbsim/lbsim_GLS.py b/brahmap/lbsim/lbsim_GLS.py index a8c609d..449da80 100644 --- a/brahmap/lbsim/lbsim_GLS.py +++ b/brahmap/lbsim/lbsim_GLS.py @@ -14,12 +14,63 @@ @dataclass class LBSimGLSParameters(GLSParameters): + """A class to encapsulate the parameters used for GLS map-making with + `litebird_sim` data + + Parameters + ---------- + solver_type : SolverType + _description_ + use_iterative_solver : bool + _description_ + isolver_threshold : float + _description_ + isolver_max_iterations : int + _description_ + callback_function : Callable + _description_ + return_processed_samples : bool + _description_ + return_hit_map : bool + _description_ + return_processed_samples : bool + _description_ + output_coordinate_system : lbs.CoordinateSystem + _description_ + """ + return_processed_samples: bool = False output_coordinate_system: lbs.CoordinateSystem = lbs.CoordinateSystem.Galactic @dataclass class LBSimGLSResult(GLSResult): + """A class to store the results of the GLs map-making done with `litebird_sim` data + + Parameters + ---------- + solver_type : SolverType + _description_ + npix : int + _description_ + new_npix : int + _description_ + GLS_maps : np.ndarray + _description_ + hit_map : np.ndarray + _description_ + convergence_status : bool + _description_ + num_iterations : int + _description_ + GLSParameters : GLSParameters + _description_ + nside : int + _description_ + coordinate_system : lbs.CoordinateSystem + _description_ + """ + nside: int coordinate_system: lbs.CoordinateSystem @@ -30,12 +81,42 @@ def LBSim_compute_GLS_maps( pointings: Union[np.ndarray, List[np.ndarray], None] = None, hwp: Optional[lbs.HWP] = None, components: Union[str, List[str]] = "tod", - pointings_flag: Union[np.ndarray, None] = None, + pointings_flag: Optional[np.ndarray] = None, inv_noise_cov_operator: Union[DTypeNoiseCov, DTypeLBSNoiseCov, None] = None, threshold: float = 1.0e-5, - dtype_float: Union[DTypeFloat, None] = None, + dtype_float: Optional[DTypeFloat] = None, LBSim_gls_parameters: LBSimGLSParameters = LBSimGLSParameters(), ) -> Union[LBSimGLSResult, tuple[LBSimProcessTimeSamples, LBSimGLSResult]]: + """_summary_ + + Parameters + ---------- + nside : int + _description_ + observations : Union[lbs.Observation, List[lbs.Observation]] + _description_ + pointings : Union[np.ndarray, List[np.ndarray], None], optional + _description_, by default None + hwp : Optional[lbs.HWP], optional + _description_, by default None + components : Union[str, List[str]], optional + _description_, by default "tod" + pointings_flag : Optional[np.ndarray], optional + _description_, by default None + inv_noise_cov_operator : Union[DTypeNoiseCov, DTypeLBSNoiseCov, None], optional + _description_, by default None + threshold : float, optional + _description_, by default 1.0e-5 + dtype_float : Optional[DTypeFloat], optional + _description_, by default None + LBSim_gls_parameters : LBSimGLSParameters, optional + _description_, by default LBSimGLSParameters() + + Returns + ------- + Union[LBSimGLSResult, tuple[LBSimProcessTimeSamples, LBSimGLSResult]] + _description_ + """ if inv_noise_cov_operator is None: noise_weights = None else: diff --git a/brahmap/lbsim/lbsim_noise_operators.py b/brahmap/lbsim/lbsim_noise_operators.py index 06f4d97..41a740e 100644 --- a/brahmap/lbsim/lbsim_noise_operators.py +++ b/brahmap/lbsim/lbsim_noise_operators.py @@ -1,8 +1,10 @@ -from typing import List, Union, Literal, Dict, Any +from typing import List, Union, Literal, Dict, Any, Optional import numpy as np import litebird_sim as lbs +from ..base import InvNoiseCovLinearOperator + from ..core import ( InvNoiseCovLO_Diagonal, InvNoiseCovLO_Circulant, @@ -10,18 +12,32 @@ BlockDiagInvNoiseCovLO, ) +from ..math import DTypeFloat + class LBSim_InvNoiseCovLO_UnCorr(BlockDiagInvNoiseCovLO): - """The assumption is that at a given MPI process, all observations - contain same set of detectors""" + """_summary_ + + The assumption is that at a given MPI process, all observations + contain same set of detectors + + Parameters + ---------- + obs : Union[lbs.Observation, List[lbs.Observation]] + _description_ + noise_variance : Optional[dict], optional + _description_, by default None + dtype : DTypeFloat, optional + _description_, by default np.float64 + """ # Keep a note of the hard-coded factor of 1e4 def __init__( self, obs: Union[lbs.Observation, List[lbs.Observation]], - noise_variance: Union[dict, None] = None, - dtype=np.float64, + noise_variance: Optional[dict] = None, + dtype: DTypeFloat = np.float64, ): if isinstance(obs, lbs.Observation): obs_list = [obs] @@ -59,12 +75,26 @@ def __init__( class LBSim_InvNoiseCovLO_Circulant(BlockDiagInvNoiseCovLO): + """_summary_ + + Parameters + ---------- + obs : Union[lbs.Observation, List[lbs.Observation]] + _description_ + input : Union[dict, Union[np.ndarray, List]] + _description_ + input_type : Literal["covariance", "power_spectrum"], optional + _description_, by default "power_spectrum" + dtype : DTypeFloat, optional + _description_, by default np.float64 + """ + def __init__( self, obs: Union[lbs.Observation, List[lbs.Observation]], input: Union[dict, Union[np.ndarray, List]], input_type: Literal["covariance", "power_spectrum"] = "power_spectrum", - dtype=np.float64, + dtype: DTypeFloat = np.float64, ): if isinstance(obs, lbs.Observation): obs_list = [obs] @@ -136,15 +166,33 @@ def _resize_input(self, new_size, input, input_type, dtype): class LBSim_InvNoiseCovLO_Toeplitz(BlockDiagInvNoiseCovLO): - """Note that the observation length is either n or n-1.""" + """_summary_ + + Note that the observation length is either n or n-1. + + Parameters + ---------- + obs : Union[lbs.Observation, List[lbs.Observation]] + _description_ + input : Union[dict, Union[np.ndarray, List]] + _description_ + input_type : Literal["covariance", "power_spectrum"], optional + _description_, by default "power_spectrum" + operator : InvNoiseCovLinearOperator, optional + _description_, by default InvNoiseCovLO_Toeplitz01 + dtype : DTypeFloat, optional + _description_, by default np.float64 + extra_kwargs : Dict[str, Any], optional + _description_, by default {} + """ def __init__( self, obs: Union[lbs.Observation, List[lbs.Observation]], input: Union[dict, Union[np.ndarray, List]], input_type: Literal["covariance", "power_spectrum"] = "power_spectrum", - operator=InvNoiseCovLO_Toeplitz01, - dtype=np.float64, + operator: InvNoiseCovLinearOperator = InvNoiseCovLO_Toeplitz01, + dtype: DTypeFloat = np.float64, extra_kwargs: Dict[str, Any] = {}, ): if isinstance(obs, lbs.Observation): diff --git a/brahmap/lbsim/lbsim_process_time_samples.py b/brahmap/lbsim/lbsim_process_time_samples.py index a6276c2..c5d1c9b 100644 --- a/brahmap/lbsim/lbsim_process_time_samples.py +++ b/brahmap/lbsim/lbsim_process_time_samples.py @@ -5,21 +5,47 @@ import litebird_sim as lbs from ..core import SolverType, ProcessTimeSamples +from ..math import DTypeFloat class LBSimProcessTimeSamples(ProcessTimeSamples): + """A class to store the pre-processed and pre-computed arrays from `litebird_sim` observations. + + Parameters + ---------- + nside : int + Nside of the healpix map + observations : Union[lbs.Observation, List[lbs.Observation]] + An instance of the `Observation` class or a list of the same + pointings : Union[np.ndarray, List[np.ndarray], None], optional + _description_, by default None + hwp : Optional[lbs.HWP], optional + _description_, by default None + pointings_flag : Optional[np.ndarray], optional + _description_, by default None + solver_type : SolverType, optional + _description_, by default SolverType.IQU + noise_weights : Optional[np.ndarray], optional + _description_, by default None + output_coordinate_system : lbs.CoordinateSystem, optional + _description_, by default lbs.CoordinateSystem.Galactic + threshold : float, optional + _description_, by default 1.0e-5 + dtype_float : DTypeFloat, optional + _description_, by default np.float64 + """ def __init__( self, nside: int, observations: Union[lbs.Observation, List[lbs.Observation]], pointings: Union[np.ndarray, List[np.ndarray], None] = None, hwp: Optional[lbs.HWP] = None, - pointings_flag: Union[np.ndarray, None] = None, + pointings_flag: Optional[np.ndarray] = None, solver_type: SolverType = SolverType.IQU, - noise_weights: Union[np.ndarray, None] = None, + noise_weights: Optional[np.ndarray] = None, output_coordinate_system: lbs.CoordinateSystem = lbs.CoordinateSystem.Galactic, threshold: float = 1.0e-5, - dtype_float=np.float64, + dtype_float: DTypeFloat = np.float64, ): self.__nside = nside self.__coordinate_system = output_coordinate_system @@ -93,12 +119,15 @@ def __init__( @property def obs_list(self): + """List of the instances of `Observation` class""" return self.__obs_list @property def nside(self): + """Nside parameter of the healpix map""" return self.__nside @property def coordinate_system(self): + """Coordinate system used in data-processing""" return self.__coordinate_system diff --git a/brahmap/math/linalg.py b/brahmap/math/linalg.py index 95becfe..a6ab3e9 100644 --- a/brahmap/math/linalg.py +++ b/brahmap/math/linalg.py @@ -5,9 +5,23 @@ import scipy.sparse.linalg from brahmap import MPI_UTILS +from ..base import LinearOperator -def parallel_norm(x: np.ndarray): +def parallel_norm(x: np.ndarray) -> float: + """A replacement of `np.linalg.norm` to compute 2-norm of a vector + distributed among multiple MPI processes + + Parameters + ---------- + x : np.ndarray + Input array + + Returns + ------- + float + The norm of vector `x` + """ sqnorm = x.dot(x) sqnorm = MPI_UTILS.comm.allreduce(sqnorm) ret = np.sqrt(sqnorm) @@ -15,16 +29,46 @@ def parallel_norm(x: np.ndarray): def cg( - A, - b, - x0=None, - rtol=1.0e-12, - atol=1.0e-12, - maxiter=100, - M=None, - callback=None, - parallel=True, + A: LinearOperator, + b: np.ndarray, + x0: np.ndarray = None, + rtol: float = 1.0e-12, + atol: float = 1.0e-12, + maxiter: int = 100, + M: LinearOperator = None, + callback: Callable = None, + parallel: bool = True, ): + """A replacement of `scipy.sparse.linalg.cg` where `np.linalg.norm` is + replaced with `brahmap.math.parallel_norm` when the parameter `parallel` + is set `True` + + Parameters + ---------- + A : LinearOperator + _description_ + b : np.ndarray + _description_ + x0 : np.ndarray, optional + _description_, by default None + rtol : float, optional + _description_, by default 1.0e-12 + atol : float, optional + _description_, by default 1.0e-12 + maxiter : int, optional + _description_, by default 100 + M : LinearOperator, optional + _description_, by default None + callback : Callable, optional + _description_, by default None + parallel : bool, optional + _description_, by default True + + Returns + ------- + _type_ + _description_ + """ A, M, x, b, postprocess = scipy.sparse.linalg._isolve.utils.make_system( A, M, diff --git a/brahmap/utilities/__init__.py b/brahmap/utilities/__init__.py index 77025fa..4a13bad 100644 --- a/brahmap/utilities/__init__.py +++ b/brahmap/utilities/__init__.py @@ -8,10 +8,6 @@ from .tools import ( bash_colors, modify_numpy_context, - TypeChangeWarning, - LowerTypeCastWarning, - filter_warnings, - ShapeError, profile_run, output_profile, ) @@ -30,10 +26,6 @@ # tools.py "bash_colors", "modify_numpy_context", - "TypeChangeWarning", - "LowerTypeCastWarning", - "filter_warnings", - "ShapeError", "profile_run", "output_profile", ] diff --git a/brahmap/utilities/tools.py b/brahmap/utilities/tools.py index eae8b75..1b3db91 100644 --- a/brahmap/utilities/tools.py +++ b/brahmap/utilities/tools.py @@ -1,4 +1,3 @@ -import warnings import numpy as np from ..math import parallel_norm @@ -41,6 +40,8 @@ def underline(self, string): class modify_numpy_context(object): + """A context manager that replaces `np.linalg.norm` with `parallel_norm`""" + def __init__(self): self.parallel_norm = parallel_norm self.original_norm = np.linalg.norm @@ -52,52 +53,13 @@ def __exit__(self, exc_type, exc_val, exc_tb): np.linalg.norm = self.original_norm -class TypeChangeWarning(Warning): - def __init__(self, message): - self.message = message - - def __str__(self): - return repr(self.message) - - -class LowerTypeCastWarning(Warning): - def __init__(self, message): - self.message = message - - def __str__(self): - return repr(self.message) - - -def filter_warnings(wfilter): - """ - wfilter: {string} - - "ignore": never print matching warnings; - - "always": always print matching warnings - - """ - warnings.simplefilter(wfilter) - - -class ShapeError(Exception): - """ - Exception class for handling shape mismatch errors. - - Exception raised when defining a linear operator of the wrong shape or - multiplying a linear operator with a vector of the wrong shape. - - """ - - def __init__(self, value): - super(ShapeError, self).__init__() - self.value = value - - def __str__(self): - return repr(self.value) - - def profile_run(): - """ - Profile the execution with :mod:`cProfile` + """Profile the execution with module `cProfile` + + Returns + ------- + _type_ + _description_ """ import cProfile @@ -106,14 +68,12 @@ def profile_run(): def output_profile(pr): - """ - Output of the profiling with :func:`profile_run`. - - **Parameter** - - - ``pr``: - instance returned by :func:`profile_run` + """Output of the profiling with `profile_run`. + Parameters + ---------- + pr : _type_ + _description_ """ import pstats import io diff --git a/brahmap/utilities/visualizations.py b/brahmap/utilities/visualizations.py index b4cc796..9199516 100644 --- a/brahmap/utilities/visualizations.py +++ b/brahmap/utilities/visualizations.py @@ -1,10 +1,26 @@ import numpy as np import matplotlib.pyplot as plt +from ..base import LinearOperator + def plot_LinearOperator( - operator: "LinearOperator", # noqa + operator: LinearOperator, ): + """A utility function to visualize BrahMap linear operators. Make sure that `matplotlib` is installed if you want to use it. + + !!! Warning + + This method first allocates a NumPy array of shape `self.shape` + and data-type `self.dtype`, and then fills them with numbers. As + such it can occupy an enormous amount of memory. Don't use it + unless you understand the risk! + + Parameters + ---------- + operator : LinearOperator + _description_ + """ plt.figure() plt.imshow(operator.to_array()) diff --git a/docs/api_reference/base_operators/BaseBlockDiagInvNoiseCovLinearOperator.md b/docs/api_reference/base_operators/BaseBlockDiagInvNoiseCovLinearOperator.md new file mode 100644 index 0000000..2096716 --- /dev/null +++ b/docs/api_reference/base_operators/BaseBlockDiagInvNoiseCovLinearOperator.md @@ -0,0 +1,3 @@ +# `brahmap.base.BaseBlockDiagInvNoiseCovLinearOperator` + +::: brahmap.base.BaseBlockDiagInvNoiseCovLinearOperator diff --git a/docs/api_reference/base_operators/BaseBlockDiagNoiseCovLinearOperator.md b/docs/api_reference/base_operators/BaseBlockDiagNoiseCovLinearOperator.md new file mode 100644 index 0000000..50ea28a --- /dev/null +++ b/docs/api_reference/base_operators/BaseBlockDiagNoiseCovLinearOperator.md @@ -0,0 +1,3 @@ +# `brahmap.base.BaseBlockDiagNoiseCovLinearOperator` + +::: brahmap.base.BaseBlockDiagNoiseCovLinearOperator diff --git a/docs/api_reference/base_operators/BaseLinearOperator.md b/docs/api_reference/base_operators/BaseLinearOperator.md new file mode 100644 index 0000000..5ce4a48 --- /dev/null +++ b/docs/api_reference/base_operators/BaseLinearOperator.md @@ -0,0 +1,3 @@ +# `brahmap.base.BaseLinearOperator` + +::: brahmap.base.BaseLinearOperator diff --git a/docs/api_reference/base_operators/BlockDiagonalLinearOperator.md b/docs/api_reference/base_operators/BlockDiagonalLinearOperator.md new file mode 100644 index 0000000..fd4b3ff --- /dev/null +++ b/docs/api_reference/base_operators/BlockDiagonalLinearOperator.md @@ -0,0 +1,3 @@ +# `brahmap.base.BlockDiagonalLinearOperator` + +::: brahmap.base.BlockDiagonalLinearOperator diff --git a/docs/api_reference/base_operators/BlockHorizontalLinearOperator.md b/docs/api_reference/base_operators/BlockHorizontalLinearOperator.md new file mode 100644 index 0000000..b53a0cc --- /dev/null +++ b/docs/api_reference/base_operators/BlockHorizontalLinearOperator.md @@ -0,0 +1,3 @@ +# `brahmap.base.BlockHorizontalLinearOperator` + +::: brahmap.base.BlockHorizontalLinearOperator diff --git a/docs/api_reference/base_operators/BlockLinearOperator.md b/docs/api_reference/base_operators/BlockLinearOperator.md new file mode 100644 index 0000000..a498f8e --- /dev/null +++ b/docs/api_reference/base_operators/BlockLinearOperator.md @@ -0,0 +1,3 @@ +# `brahmap.base.BlockLinearOperator` + +::: brahmap.base.BlockLinearOperator diff --git a/docs/api_reference/base_operators/BlockVerticalLinearOperator.md b/docs/api_reference/base_operators/BlockVerticalLinearOperator.md new file mode 100644 index 0000000..5f8e6b7 --- /dev/null +++ b/docs/api_reference/base_operators/BlockVerticalLinearOperator.md @@ -0,0 +1,3 @@ +# `brahmap.base.BlockVerticalLinearOperator` + +::: brahmap.base.BlockVerticalLinearOperator diff --git a/docs/api_reference/base_operators/DiagonalOperator.md b/docs/api_reference/base_operators/DiagonalOperator.md new file mode 100644 index 0000000..fa1cb82 --- /dev/null +++ b/docs/api_reference/base_operators/DiagonalOperator.md @@ -0,0 +1,3 @@ +# `brahmap.base.DiagonalOperator` + +::: brahmap.base.DiagonalOperator diff --git a/docs/api_reference/base_operators/IdentityOperator.md b/docs/api_reference/base_operators/IdentityOperator.md new file mode 100644 index 0000000..219f21c --- /dev/null +++ b/docs/api_reference/base_operators/IdentityOperator.md @@ -0,0 +1,3 @@ +# `brahmap.base.IdentityOperator` + +::: brahmap.base.IdentityOperator diff --git a/docs/api_reference/base_operators/InvNoiseCovLinearOperator.md b/docs/api_reference/base_operators/InvNoiseCovLinearOperator.md new file mode 100644 index 0000000..7ac09f6 --- /dev/null +++ b/docs/api_reference/base_operators/InvNoiseCovLinearOperator.md @@ -0,0 +1,3 @@ +# `brahmap.base.InvNoiseCovLinearOperator` + +::: brahmap.base.InvNoiseCovLinearOperator diff --git a/docs/api_reference/base_operators/InverseLO.md b/docs/api_reference/base_operators/InverseLO.md new file mode 100644 index 0000000..d832dd4 --- /dev/null +++ b/docs/api_reference/base_operators/InverseLO.md @@ -0,0 +1,3 @@ +# `brahmap.base.InverseLO` + +::: brahmap.base.InverseLO diff --git a/docs/api_reference/base_operators/LinearOperator.md b/docs/api_reference/base_operators/LinearOperator.md new file mode 100644 index 0000000..4ea0be4 --- /dev/null +++ b/docs/api_reference/base_operators/LinearOperator.md @@ -0,0 +1,3 @@ +# `brahmap.base.LinearOperator` + +::: brahmap.base.LinearOperator diff --git a/docs/api_reference/base_operators/MatrixLinearOperator.md b/docs/api_reference/base_operators/MatrixLinearOperator.md new file mode 100644 index 0000000..7effa94 --- /dev/null +++ b/docs/api_reference/base_operators/MatrixLinearOperator.md @@ -0,0 +1,3 @@ +# `brahmap.base.MatrixLinearOperator` + +::: brahmap.base.MatrixLinearOperator diff --git a/docs/api_reference/base_operators/NoiseCovLinearOperator.md b/docs/api_reference/base_operators/NoiseCovLinearOperator.md new file mode 100644 index 0000000..ea452b8 --- /dev/null +++ b/docs/api_reference/base_operators/NoiseCovLinearOperator.md @@ -0,0 +1,3 @@ +# `brahmap.base.NoiseCovLinearOperator` + +::: brahmap.base.NoiseCovLinearOperator diff --git a/docs/api_reference/base_operators/ReducedLinearOperator.md b/docs/api_reference/base_operators/ReducedLinearOperator.md new file mode 100644 index 0000000..3bc5396 --- /dev/null +++ b/docs/api_reference/base_operators/ReducedLinearOperator.md @@ -0,0 +1,3 @@ +# `brahmap.base.ReducedLinearOperator` + +::: brahmap.base.ReducedLinearOperator diff --git a/docs/api_reference/base_operators/SymmetricallyReducedLinearOperator.md b/docs/api_reference/base_operators/SymmetricallyReducedLinearOperator.md new file mode 100644 index 0000000..1ebd1bd --- /dev/null +++ b/docs/api_reference/base_operators/SymmetricallyReducedLinearOperator.md @@ -0,0 +1,3 @@ +# `brahmap.base.SymmetricallyReducedLinearOperator` + +::: brahmap.base.SymmetricallyReducedLinearOperator diff --git a/docs/api_reference/base_operators/ZeroOperator.md b/docs/api_reference/base_operators/ZeroOperator.md new file mode 100644 index 0000000..e76d286 --- /dev/null +++ b/docs/api_reference/base_operators/ZeroOperator.md @@ -0,0 +1,3 @@ +# `brahmap.base.ZeroOperator` + +::: brahmap.base.ZeroOperator diff --git a/docs/api_reference/base_operators/aslinearoperator.md b/docs/api_reference/base_operators/aslinearoperator.md new file mode 100644 index 0000000..02210de --- /dev/null +++ b/docs/api_reference/base_operators/aslinearoperator.md @@ -0,0 +1,3 @@ +# `brahmap.base.aslinearoperator` + +::: brahmap.base.aslinearoperator diff --git a/docs/api_reference/base_operators/index.md b/docs/api_reference/base_operators/index.md index f05592d..6aad02f 100644 --- a/docs/api_reference/base_operators/index.md +++ b/docs/api_reference/base_operators/index.md @@ -2,8 +2,27 @@ ## `linop` sub-module -::: brahmap.base.linop +- [`BaseLinearOperator`](BaseLinearOperator.md) +- [`LinearOperator`](LinearOperator.md) +- [`IdentityOperator`](IdentityOperator.md) +- [`DiagonalOperator`](DiagonalOperator.md) +- [`MatrixLinearOperator`](MatrixLinearOperator.md) +- [`ZeroOperator`](ZeroOperator.md) +- [`InverseLO`](InverseLO.md) +- [`ReducedLinearOperator`](ReducedLinearOperator.md) +- [`SymmetricallyReducedLinearOperator`](SymmetricallyReducedLinearOperator.md) +- [`aslinearoperator`](aslinearoperator.md) ## `blkop` sub-module -::: brahmap.base.blkop +- [`BlockLinearOperator`](BlockLinearOperator.md) +- [`BlockDiagonalLinearOperator`](BlockDiagonalLinearOperator.md) +- [`BlockHorizontalLinearOperator`](BlockHorizontalLinearOperator.md) +- [`BlockVerticalLinearOperator`](BlockVerticalLinearOperator.md) + +## Base noise covariance (and inverse) operators + +- [`NoiseCovLinearOperator`](NoiseCovLinearOperator.md) +- [`InvNoiseCovLinearOperator`](InvNoiseCovLinearOperator.md) +- [`BaseBlockDiagNoiseCovLinearOperator`](BaseBlockDiagNoiseCovLinearOperator.md) +- [`BaseBlockDiagInvNoiseCovLinearOperator`](BaseBlockDiagInvNoiseCovLinearOperator.md) diff --git a/docs/api_reference/core/BlockDiagInvNoiseCovLO.md b/docs/api_reference/core/BlockDiagInvNoiseCovLO.md new file mode 100644 index 0000000..9ece8cc --- /dev/null +++ b/docs/api_reference/core/BlockDiagInvNoiseCovLO.md @@ -0,0 +1,3 @@ +# `brahmap.core.BlockDiagInvNoiseCovLO` + +::: brahmap.core.BlockDiagInvNoiseCovLO diff --git a/docs/api_reference/core/BlockDiagNoiseCovLO.md b/docs/api_reference/core/BlockDiagNoiseCovLO.md new file mode 100644 index 0000000..28fed32 --- /dev/null +++ b/docs/api_reference/core/BlockDiagNoiseCovLO.md @@ -0,0 +1,3 @@ +# `brahmap.core.BlockDiagNoiseCovLO` + +::: brahmap.core.BlockDiagNoiseCovLO diff --git a/docs/api_reference/core/InvNoiseCovLO_Circulant.md b/docs/api_reference/core/InvNoiseCovLO_Circulant.md new file mode 100644 index 0000000..fae963e --- /dev/null +++ b/docs/api_reference/core/InvNoiseCovLO_Circulant.md @@ -0,0 +1,3 @@ +# `brahmap.core.InvNoiseCovLO_Circulant` + +::: brahmap.core.InvNoiseCovLO_Circulant diff --git a/docs/api_reference/core/InvNoiseCovLO_Toeplitz01.md b/docs/api_reference/core/InvNoiseCovLO_Toeplitz01.md new file mode 100644 index 0000000..8bfdf1e --- /dev/null +++ b/docs/api_reference/core/InvNoiseCovLO_Toeplitz01.md @@ -0,0 +1,3 @@ +# `brahmap.core.InvNoiseCovLO_Toeplitz01` + +::: brahmap.core.InvNoiseCovLO_Toeplitz01 diff --git a/docs/api_reference/core/NoiseCovLO_Circulant.md b/docs/api_reference/core/NoiseCovLO_Circulant.md new file mode 100644 index 0000000..1a16466 --- /dev/null +++ b/docs/api_reference/core/NoiseCovLO_Circulant.md @@ -0,0 +1,3 @@ +# `brahmap.core.NoiseCovLO_Circulant` + +::: brahmap.core.NoiseCovLO_Circulant diff --git a/docs/api_reference/core/NoiseCovLO_Diagonal.md b/docs/api_reference/core/NoiseCovLO_Diagonal.md new file mode 100644 index 0000000..34abfe0 --- /dev/null +++ b/docs/api_reference/core/NoiseCovLO_Diagonal.md @@ -0,0 +1,3 @@ +# `brahmap.core.NoiseCovLO_Diagonal` + +::: brahmap.core.NoiseCovLO_Diagonal diff --git a/docs/api_reference/core/NoiseCovLO_Toeplitz01.md b/docs/api_reference/core/NoiseCovLO_Toeplitz01.md new file mode 100644 index 0000000..775cec6 --- /dev/null +++ b/docs/api_reference/core/NoiseCovLO_Toeplitz01.md @@ -0,0 +1,3 @@ +# `brahmap.core.NoiseCovLO_Toeplitz01` + +::: brahmap.core.NoiseCovLO_Toeplitz01 diff --git a/docs/api_reference/core/index.md b/docs/api_reference/core/index.md index c7d7853..4b34ea6 100644 --- a/docs/api_reference/core/index.md +++ b/docs/api_reference/core/index.md @@ -1,4 +1,4 @@ -# Core Interfaces +# Core Interface ## Data pre-processing @@ -12,7 +12,14 @@ ## Noise covariance (and their inverse) operators +- [`NoiseCovLO_Diagonal`](NoiseCovLO_Diagonal.md) +- [`NoiseCovLO_Circulant`](NoiseCovLO_Circulant.md) +- [`NoiseCovLO_Toeplitz01`](NoiseCovLO_Toeplitz01.md) +- [`BlockDiagNoiseCovLO`](BlockDiagNoiseCovLO.md) - [`InvNoiseCovLO_Diagonal`](InvNoiseCovLO_Diagonal.md) +- [`InvNoiseCovLO_Circulant`](InvNoiseCovLO_Circulant.md) +- [`InvNoiseCovLO_Toeplitz01`](InvNoiseCovLO_Toeplitz01.md) +- [`BlockDiagInvNoiseCovLO`](BlockDiagInvNoiseCovLO.md) ## GLS map-making functions and tools diff --git a/docs/api_reference/cpp_extensions/index.md b/docs/api_reference/cpp_extensions/index.md deleted file mode 100644 index 9735fb8..0000000 --- a/docs/api_reference/cpp_extensions/index.md +++ /dev/null @@ -1,3 +0,0 @@ -# C++ Extensions - -Under development! diff --git a/docs/api_reference/index.md b/docs/api_reference/index.md index ff1aa1c..8d7c6ea 100644 --- a/docs/api_reference/index.md +++ b/docs/api_reference/index.md @@ -1,19 +1,19 @@ # API reference -## Core API +## `litebird_sim` API -- [Core Interface](./core/index.md) - [`litebird_sim` Interface](./lbsim/index.md) + +## Core API + +- [Core Map-making Interface](./core/index.md) +- [Math Functions](./math_functions/index.md) - [Utilities](./utilities/index.md) ## Base API - [Base Operators](./base_operators/index.md) - -## C++ Extensions - -- [Math functions](./math_functions/index.md) -- [C++ extensions](./cpp_extensions/index.md) +- [Miscellaneous](./misc/index.md) ## C++ API diff --git a/docs/api_reference/lbsim/LBSim_InvNoiseCovLO_Circulant.md b/docs/api_reference/lbsim/LBSim_InvNoiseCovLO_Circulant.md new file mode 100644 index 0000000..eaaad6f --- /dev/null +++ b/docs/api_reference/lbsim/LBSim_InvNoiseCovLO_Circulant.md @@ -0,0 +1,3 @@ +# `brahmap.lbsim.LBSim_InvNoiseCovLO_Circulant` + +::: brahmap.lbsim.LBSim_InvNoiseCovLO_Circulant diff --git a/docs/api_reference/lbsim/LBSim_InvNoiseCovLO_Toeplitz.md b/docs/api_reference/lbsim/LBSim_InvNoiseCovLO_Toeplitz.md new file mode 100644 index 0000000..15b927d --- /dev/null +++ b/docs/api_reference/lbsim/LBSim_InvNoiseCovLO_Toeplitz.md @@ -0,0 +1,3 @@ +# `brahmap.lbsim.LBSim_InvNoiseCovLO_Toeplitz` + +::: brahmap.lbsim.LBSim_InvNoiseCovLO_Toeplitz diff --git a/docs/api_reference/lbsim/index.md b/docs/api_reference/lbsim/index.md index c9f645b..7b8f4fc 100644 --- a/docs/api_reference/lbsim/index.md +++ b/docs/api_reference/lbsim/index.md @@ -7,6 +7,8 @@ ## Noise covariance (and their inverse) operators - [`LBSim_InvNoiseCovLO_UnCorr`](LBSim_InvNoiseCovLO_UnCorr.md) +- [`LBSim_InvNoiseCovLO_Circulant`](LBSim_InvNoiseCovLO_Circulant.md) +- [`LBSim_InvNoiseCovLO_Toeplitz`](LBSim_InvNoiseCovLO_Toeplitz.md) ## GLS map-making functions and tools diff --git a/docs/api_reference/math_functions/cg.md b/docs/api_reference/math_functions/cg.md new file mode 100644 index 0000000..ce0ac5c --- /dev/null +++ b/docs/api_reference/math_functions/cg.md @@ -0,0 +1,3 @@ +# `brahmap.math.cg` + +::: brahmap.math.cg diff --git a/docs/api_reference/math_functions/index.md b/docs/api_reference/math_functions/index.md index d3a0ab6..a6b2e07 100644 --- a/docs/api_reference/math_functions/index.md +++ b/docs/api_reference/math_functions/index.md @@ -1,8 +1,36 @@ # Math Functions -## Linear Algebra +## Linear algebra - [`parallel_norm`](parallel_norm.md) +- [`cg`](cg.md) + +## Unary functions + +These functions are element-wise unary functions, and have common signature: + +```python +function_name( + size, # size of the input/output array + vec, # input array + result, # output array containing the result (overwritten) +) +``` + +Following functions are available: + +- `brahmap.math.sin`: sine function +- `brahmap.math.cos`: cosine function +- `brahmap.math.tan`: tangent function +- `brahmap.math.asin`: arcsine function +- `brahmap.math.acos`: arccosine function +- `brahmap.math.atan`: arctangent function +- `brahmap.math.exp`: exponential function, $e^x$ +- `brahmap.math.exp2`: exponential function with base 2, $2^x$ +- `brahmap.math.log`: natural logarithm function +- `brahmap.math.log2`: base-2 logarithm function +- `brahmap.math.sqrt`: square-root function +- `brahmap.math.cbrt`: cube-root function ## `dtype` hints diff --git a/docs/api_reference/misc/LowerTypeCastWarning.md b/docs/api_reference/misc/LowerTypeCastWarning.md new file mode 100644 index 0000000..a5ec44a --- /dev/null +++ b/docs/api_reference/misc/LowerTypeCastWarning.md @@ -0,0 +1,3 @@ +# `brahmap.base.LowerTypeCastWarning` + +::: brahmap.base.LowerTypeCastWarning diff --git a/docs/api_reference/misc/ShapeError.md b/docs/api_reference/misc/ShapeError.md new file mode 100644 index 0000000..7a47c5e --- /dev/null +++ b/docs/api_reference/misc/ShapeError.md @@ -0,0 +1,3 @@ +# `brahmap.base.ShapeError` + +::: brahmap.base.ShapeError diff --git a/docs/api_reference/misc/TypeChangeWarning.md b/docs/api_reference/misc/TypeChangeWarning.md new file mode 100644 index 0000000..827a7ba --- /dev/null +++ b/docs/api_reference/misc/TypeChangeWarning.md @@ -0,0 +1,3 @@ +# `brahmap.base.TypeChangeWarning` + +::: brahmap.base.TypeChangeWarning diff --git a/docs/api_reference/misc/filter_warnings.md b/docs/api_reference/misc/filter_warnings.md new file mode 100644 index 0000000..4671d74 --- /dev/null +++ b/docs/api_reference/misc/filter_warnings.md @@ -0,0 +1,3 @@ +# `brahmap.base.filter_warnings` + +::: brahmap.base.filter_warnings diff --git a/docs/api_reference/misc/index.md b/docs/api_reference/misc/index.md new file mode 100644 index 0000000..9545e00 --- /dev/null +++ b/docs/api_reference/misc/index.md @@ -0,0 +1,6 @@ +# Miscellaneous + +- [`TypeChangeWarning`](TypeChangeWarning.md) +- [`LowerTypeCastWarning`](LowerTypeCastWarning.md) +- [`filter_warnings`](filter_warnings.md) +- [`ShapeError`](ShapeError.md) diff --git a/docs/api_reference/utilities/LowerTypeCastWarning.md b/docs/api_reference/utilities/LowerTypeCastWarning.md deleted file mode 100644 index 6190cf1..0000000 --- a/docs/api_reference/utilities/LowerTypeCastWarning.md +++ /dev/null @@ -1,3 +0,0 @@ -# `brahmap.utilities.LowerTypeCastWarning` - -::: brahmap.utilities.LowerTypeCastWarning diff --git a/docs/api_reference/utilities/ShapeError.md b/docs/api_reference/utilities/ShapeError.md deleted file mode 100644 index 80fafe6..0000000 --- a/docs/api_reference/utilities/ShapeError.md +++ /dev/null @@ -1,3 +0,0 @@ -# `brahmap.utilities.ShapeError` - -::: brahmap.utilities.ShapeError diff --git a/docs/api_reference/utilities/TypeChangeWarning.md b/docs/api_reference/utilities/TypeChangeWarning.md deleted file mode 100644 index 6bceadd..0000000 --- a/docs/api_reference/utilities/TypeChangeWarning.md +++ /dev/null @@ -1,3 +0,0 @@ -# `brahmap.utilities.TypeChangeWarning` - -::: brahmap.utilities.TypeChangeWarning diff --git a/docs/api_reference/utilities/bash_colors.md b/docs/api_reference/utilities/bash_colors.md new file mode 100644 index 0000000..5a769b3 --- /dev/null +++ b/docs/api_reference/utilities/bash_colors.md @@ -0,0 +1,3 @@ +# `brahmap.utilities.bash_colors` + +::: brahmap.utilities.bash_colors diff --git a/docs/api_reference/utilities/index.md b/docs/api_reference/utilities/index.md index dbe8085..3e055ef 100644 --- a/docs/api_reference/utilities/index.md +++ b/docs/api_reference/utilities/index.md @@ -2,7 +2,9 @@ ## Tools +- [`bash_colors`](bash_colors.md) - [`modify_numpy_context`](modify_numpy_context.md) -- [`TypeChangeWarning`](TypeChangeWarning.md) -- [`LowerTypeCastWarning`](LowerTypeCastWarning.md) -- [`ShapeError`](ShapeError.md) + +## Visualization + +- [`plot_LinearOperator`](plot_LinearOperator.md) diff --git a/docs/api_reference/utilities/plot_LinearOperator.md b/docs/api_reference/utilities/plot_LinearOperator.md new file mode 100644 index 0000000..ce4f247 --- /dev/null +++ b/docs/api_reference/utilities/plot_LinearOperator.md @@ -0,0 +1,3 @@ +# `brahmap.utilities.plot_LinearOperator` + +::: brahmap.utilities.plot_LinearOperator diff --git a/docs/javascripts/mathjax.js b/docs/javascripts/mathjax.js new file mode 100644 index 0000000..5b34852 --- /dev/null +++ b/docs/javascripts/mathjax.js @@ -0,0 +1,19 @@ +window.MathJax = { + tex: { + inlineMath: [["\\(", "\\)"]], + displayMath: [["\\[", "\\]"]], + processEscapes: true, + processEnvironments: true + }, + options: { + ignoreHtmlClass: ".*|", + processHtmlClass: "arithmatex" + } +}; + +document$.subscribe(() => { + MathJax.startup.output.clearCache() + MathJax.typesetClear() + MathJax.texReset() + MathJax.typesetPromise() +}) diff --git a/mkdocs.yaml b/mkdocs.yaml index 6714dea..f258e3b 100644 --- a/mkdocs.yaml +++ b/mkdocs.yaml @@ -42,7 +42,7 @@ theme: - search.highlight # highlight the search terms - navigation.top - navigation.tabs - - navigation.expand + # - navigation.expand - navigation.footer - navigation.indexes - navigation.path @@ -94,8 +94,7 @@ markdown_extensions: ### Matjax support extra_javascript: - javascripts/mathjax.js - - https://polyfill.io/v3/polyfill.min.js?features=es6 - - https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js + - https://unpkg.com/mathjax@3/es5/tex-mml-chtml.js plugins: @@ -105,6 +104,13 @@ plugins: python: options: docstring_style: numpy + allow_inspection: true # import the modules and inspect them + # force_inspection: true # Force import (see if it works for extensions) + show_signature: true # show function signature + show_signature_annotations: true # show signature annotation + inherited_members: false + # summary: true + show_docstring_raises: false - search # Enable search ### Page tree @@ -125,39 +131,75 @@ nav: - API Reference: - api_reference/index.md - - Core Interface: + + - litebird_sim Interface: + - api_reference/lbsim/index.md + - LBSimProcessTimeSamples: api_reference/lbsim/LBSimProcessTimeSamples.md + - LBSim_InvNoiseCovLO_UnCorr: api_reference/lbsim/LBSim_InvNoiseCovLO_UnCorr.md + - LBSim_InvNoiseCovLO_Circulant: api_reference/lbsim/LBSim_InvNoiseCovLO_Circulant.md + - LBSim_InvNoiseCovLO_Toeplitz: api_reference/lbsim/LBSim_InvNoiseCovLO_Toeplitz.md + - LBSimGLSParameters: api_reference/lbsim/LBSimGLSParameters.md + - LBSim_compute_GLS_maps: api_reference/lbsim/LBSim_compute_GLS_maps.md + - LBSimGLSResult: api_reference/lbsim/LBSimGLSResult.md + + - Core Map-making Interface: - api_reference/core/index.md - SolverType: api_reference/core/SolverType.md - ProcessTimeSamples: api_reference/core/ProcessTimeSamples.md - PointingLO: api_reference/core/PointingLO.md - BlockDiagonalPreconditionerLO: api_reference/core/BlockDiagonalPreconditionerLO.md + - NoiseCovLO_Diagonal: api_reference/core/NoiseCovLO_Diagonal.md + - NoiseCovLO_Circulant: api_reference/core/NoiseCovLO_Circulant.md + - NoiseCovLO_Toeplitz01: api_reference/core/NoiseCovLO_Toeplitz01.md + - BlockDiagNoiseCovLO: api_reference/core/BlockDiagNoiseCovLO.md - InvNoiseCovLO_Diagonal: api_reference/core/InvNoiseCovLO_Diagonal.md + - InvNoiseCovLO_Circulant: api_reference/core/InvNoiseCovLO_Circulant.md + - InvNoiseCovLO_Toeplitz01: api_reference/core/InvNoiseCovLO_Toeplitz01.md + - BlockDiagInvNoiseCovLO: api_reference/core/BlockDiagInvNoiseCovLO.md - GLSParameters: api_reference/core/GLSParameters.md - compute_GLS_maps_from_PTS: api_reference/core/compute_GLS_maps_from_PTS.md - compute_GLS_maps: api_reference/core/compute_GLS_maps.md - separate_map_vectors: api_reference/core/separate_map_vectors.md - - GLSResult: api_reference/core/GLSResult.md - - - litebird_sim Interface: - - api_reference/lbsim/index.md - - LBSimProcessTimeSamples: api_reference/lbsim/LBSimProcessTimeSamples.md - - LBSim_InvNoiseCovLO_UnCorr: api_reference/lbsim/LBSim_InvNoiseCovLO_UnCorr.md - - LBSimGLSParameters: api_reference/lbsim/LBSimGLSParameters.md - - LBSim_compute_GLS_maps: api_reference/lbsim/LBSim_compute_GLS_maps.md - - LBSimGLSResult: api_reference/lbsim/LBSimGLSResult.md + - GLSResult: api_reference/core/GLSResult.md + - Math Functions: + - api_reference/math_functions/index.md + - parallel_norm: api_reference/math_functions/parallel_norm.md + - cg: api_reference/math_functions/cg.md + - Utilities: - api_reference/utilities/index.md + - bash_colors: api_reference/utilities/bash_colors.md - modify_numpy_context: api_reference/utilities/modify_numpy_context.md - - TypeChangeWarning: api_reference/utilities/TypeChangeWarning.md - - LowerTypeCastWarning: api_reference/utilities/LowerTypeCastWarning.md - - ShapeError: api_reference/utilities/ShapeError.md + - plot_LinearOperator: api_reference/utilities/plot_LinearOperator.md - - Base Operators: api_reference/base_operators/index.md - - Math Functions: - - api_reference/math_functions/index.md - - parallel_norm: api_reference/math_functions/parallel_norm.md - - C++ Extensions: api_reference/cpp_extensions/index.md + - Base Operators: + - api_reference/base_operators/index.md + - BaseLinearOperator: api_reference/base_operators/BaseLinearOperator.md + - LinearOperator: api_reference/base_operators/LinearOperator.md + - IdentityOperator: api_reference/base_operators/IdentityOperator.md + - DiagonalOperator: api_reference/base_operators/DiagonalOperator.md + - MatrixLinearOperator: api_reference/base_operators/MatrixLinearOperator.md + - ZeroOperator: api_reference/base_operators/ZeroOperator.md + - InverseLO: api_reference/base_operators/InverseLO.md + - ReducedLinearOperator: api_reference/base_operators/ReducedLinearOperator.md + - SymmetricallyReducedLinearOperator: api_reference/base_operators/SymmetricallyReducedLinearOperator.md + - aslinearoperator: api_reference/base_operators/aslinearoperator.md + - BlockLinearOperator: api_reference/base_operators/BlockLinearOperator.md + - BlockDiagonalLinearOperator: api_reference/base_operators/BlockDiagonalLinearOperator.md + - BlockHorizontalLinearOperator: api_reference/base_operators/BlockHorizontalLinearOperator.md + - BlockVerticalLinearOperator: api_reference/base_operators/BlockVerticalLinearOperator.md + - NoiseCovLinearOperator: api_reference/base_operators/NoiseCovLinearOperator.md + - InvNoiseCovLinearOperator: api_reference/base_operators/InvNoiseCovLinearOperator.md + - BaseBlockDiagNoiseCovLinearOperator: api_reference/base_operators/BaseBlockDiagNoiseCovLinearOperator.md + - BaseBlockDiagInvNoiseCovLinearOperator: api_reference/base_operators/BaseBlockDiagInvNoiseCovLinearOperator.md + + - Miscellaneous: + - api_reference/misc/index.md + - TypeChangeWarning: api_reference/misc/TypeChangeWarning.md + - LowerTypeCastWarning: api_reference/misc/LowerTypeCastWarning.md + - filter_warnings: api_reference/misc/filter_warnings.md + - ShapeError: api_reference/misc/ShapeError.md - Development: - development/index.md