Merge pull request #24 from hippylib/doc_build

Sphinx documentation fixed, added for algorithms, modeling and utils

docs/conf.py

@@ -18,7 +18,7 @@
 sys.path.insert(
     0,
     os.path.join(
-        os.path.abspath(os.path.dirname(os.path.dirname(__file__))), "hippylibX"
+        os.path.abspath(os.path.dirname(os.path.dirname(__file__))), "hippylibx"
     ),
 )
 sys.path.insert(0, os.path.abspath(".."))
@@ -32,6 +32,14 @@
     "scipy",
     "numpy",
 ]
+
+autodoc_default_options = {
+    "members": True,
+    "private-members": True,
+    "undoc-members": True,
+    "show-inheritance": True,
+}
+
 autodoc_default_flags = ["members", "private-members", "undoc-members"]
 autoclass_content = "both"
 
@@ -68,7 +76,6 @@
     "sphinx.ext.todo",
     "sphinx.ext.mathjax",
     "sphinx.ext.viewcode",
-    # "m2r",
     "sphinx_mdinclude",
 ]
 

hippylibX/algorithms/NewtonCG.py

@@ -310,6 +310,9 @@ def _solve_ls(self, x: list) -> list:
         return x
 
     def _solve_tr(self, x):
+        """
+        TR method has not been implemented in v0.1.0.
+        """
         rel_tol = self.parameters["rel_tolerance"]
         abs_tol = self.parameters["abs_tolerance"]
         max_iter = self.parameters["max_iter"]

hippylibX/algorithms/cgsolverSteihaug.py

@@ -46,18 +46,22 @@ class CGSolverSteihaug:
     """
     Solve the linear system :math:`A x = b` using preconditioned conjugate gradient ( :math:`B` preconditioner)
     and the Steihaug stopping criterion:
+
     - reason of termination 0: we reached the maximum number of iterations (no convergence)
     - reason of termination 1: we reduced the residual up to the given tolerance (convergence)
     - reason of termination 2: we reached a negative direction (premature termination due to not spd matrix)
     - reason of termination 3: we reached the boundary of the trust region
 
     The stopping criterion is based on either
+
     - the absolute preconditioned residual norm check: :math:`|| r^* ||_{B^{-1}} < atol`
     - the relative preconditioned residual norm check: :math:`|| r^* ||_{B^{-1}}/|| r^0 ||_{B^{-1}} < rtol,`
+
     where :math:`r^* = b - Ax^*` is the residual at convergence and :math:`r^0 = b - Ax^0` is the initial residual.
 
     The operator :code:`A` is set using the method :code:`set_operator(A)`.
     :code:`A` must provide the following two methods:
+
     - :code:`A.mult(x,y)`: `y = Ax`
     - :code:`A.init_vector(x, dim)`: initialize the vector `x` so that it is compatible with the range `(dim = 0)` or
       the domain `(dim = 1)` of :code:`A`.

hippylibX/algorithms/multivector.py

@@ -15,6 +15,10 @@
 
 
 class MultiVector:
+    """
+    Multivector object is created as a list of PETSc.petsc4py.Vec objects.
+    """
+
     def __init__(self, example_vec: petsc4py.PETSc.Vec, nvec: int):
         self.nvec = nvec
         self.data = []
@@ -25,10 +29,16 @@ def __init__(self, example_vec: petsc4py.PETSc.Vec, nvec: int):
     def createFromVec(
         cls, example_vec: petsc4py.PETSc.Vec, nvec: int
     ) -> Type["MultiVector"]:
+        """
+        Create multivector from sample petsc4py vector whose parallel distribution is to be replicated.
+        """
         return cls(example_vec, nvec)
 
     @classmethod
     def createFromMultiVec(cls, mv: Type["MultiVector"]) -> Type["MultiVector"]:
+        """
+        Create multivector from another MultiVector whose parallel distribution is to be replicated.
+        """
         mv_copy = cls(mv[0], mv.nvec)
         for i in range(mv_copy.nvec):
             mv_copy.data[i] = mv.data[i].duplicate(mv.data[i].getArray())
@@ -43,6 +53,9 @@ def __getitem__(self, k: int) -> petsc4py.PETSc.Vec:
         return self.data[k]
 
     def scale(self, alpha: Union[float, np.ndarray]) -> None:
+        """
+        Scale each value in the Multivector - either by a single float value or a numpy array of float values.
+        """
         if isinstance(alpha, float):
             for d in self.data:
                 d.scale(alpha)
@@ -51,6 +64,9 @@ def scale(self, alpha: Union[float, np.ndarray]) -> None:
                 d.scale(alpha[i])
 
     def dot(self, v: Union[petsc4py.PETSc.Vec, Type["MultiVector"]]) -> np.array:
+        """
+        Perform dot product of a MultiVector object and petsc4py Vec object, store result in numpy array.
+        """
         if isinstance(v, petsc4py.PETSc.Vec):
             return_values = np.zeros(self.nvec)
             for i in range(self.nvec):
@@ -65,10 +81,16 @@ def dot(self, v: Union[petsc4py.PETSc.Vec, Type["MultiVector"]]) -> np.array:
         return return_values
 
     def reduce(self, y: petsc4py.PETSc.Vec, alpha: np.array) -> None:
+        """
+        Reduction of petsc4py Vec using values in a numpy array stored in each data element of MultiVector object.
+        """
         for i in range(self.nvec):
             y.axpy(alpha[i], self[i])
 
     def axpy(self, alpha: Union[float, np.array], Y: Type["MultiVector"]) -> None:
+        """
+        Reduction of MultiVector object with a float or values in a numpy array stored in another MultiVector object.
+        """
         if isinstance(alpha, float):
             for i in range(self.nvec):
                 self[i].axpy(alpha, Y[i])
@@ -77,17 +99,55 @@ def axpy(self, alpha: Union[float, np.array], Y: Type["MultiVector"]) -> None:
                 self[i].axpy(alpha[i], Y[i])
 
     def norm(self, norm_type: petsc4py.PETSc.NormType) -> np.array:
+        """
+        Return numpy array containing norm of each data element of MultiVector.
+        """
         norm_vals = np.zeros(self.nvec)
         for i in range(self.nvec):
             norm_vals[i] = self[i].norm(norm_type)
         return norm_vals
 
     def Borthogonalize(self, B: petsc4py.PETSc.Mat):
+        """
+        Returns :math:`QR` decomposition of self. :math:`Q` and :math:`R` satisfy the following relations in exact arithmetic
+
+        .. math::
+            QR \\,= \\,Z, && (1),
+
+            Q^*BQ\\, = \\, I, && (2),
+
+            Q^*BZ \\, = \\,R, && (3),
+
+            ZR^{-1} \\, = \\, Q, && (4).
+
+        Returns:
+
+            :code:`Bq` of type :code:`MultiVector` -> :code:`B`:math:`^{-1}`-orthogonal vectors
+            :code:`r` of type :code:`ndarray` -> The :math:`r` of the QR decomposition.
+
+        .. note:: :code:`self` is overwritten by :math:`Q`.
+        """
+
         return self._mgs_stable(B)
 
     def _mgs_stable(
         self, B: petsc4py.PETSc.Mat
     ) -> tuple[Type["MultiVector"], np.array]:
+        """ 
+        Returns :math:`QR` decomposition of self, which satisfies conditions (1)--(4).
+        Uses Modified Gram-Schmidt with re-orthogonalization (Rutishauser variant)
+        for computing the :math:`B`-orthogonal :math:`QR` factorization.
+        
+        References:
+            1. `A.K. Saibaba, J. Lee and P.K. Kitanidis, Randomized algorithms for Generalized \
+            Hermitian Eigenvalue Problems with application to computing \
+            Karhunen-Loe've expansion http://arxiv.org/abs/1307.6885`
+            2. `W. Gander, Algorithms for the QR decomposition. Res. Rep, 80(02), 1980`
+        
+        https://github.com/arvindks/kle
+        
+        """
+
         n = self.nvec
         Bq = MultiVector(self[0], n)
         r = np.zeros((n, n), dtype="d")

hippylibX/algorithms/randomizedEigensolver.py

@@ -23,6 +23,24 @@ def doublePassG(
     k: int,
     s=1,
 ) -> tuple[np.array, MultiVector]:
+    """
+    The double pass algorithm for the GHEP as presented in [2].
+
+    Inputs:
+
+    - :code:`A`: the operator for which we need to estimate the dominant generalized eigenpairs.
+    - :code:`B`: the right-hand side operator.
+    - :code:`Binv`: the inverse of the right-hand side operator.
+    - :code:`Omega`: a random gassian matrix with :math:`m \\geq k` columns.
+    - :code:`k`: the number of eigenpairs to extract.
+    - :code:`s`: the number of power iterations for selecting the subspace.
+
+    Outputs:
+
+    - :code:`d`: the estimate of the :math:`k` dominant eigenvalues of :math:`A`.
+    - :code:`U`: the estimate of the :math:`k` dominant eigenvectors of :math:`A,\\, U^T B U = I_k`.
+    """
+
     nvec = Omega.nvec
     assert nvec >= k
     Ybar = MultiVector.createFromVec(Omega[0], nvec)

hippylibX/modeling/Regularization.py

@@ -16,6 +16,10 @@
 
 # functional handler for prior:
 class H1TikhonvFunctional:
+    """
+    The functional handler for the Variational Regularization prior.
+    """
+
     def __init__(self, gamma, delta, m0=None):
         self.gamma = gamma  # These are dlx Constant, Expression, or Function
         self.delta = delta

hippylibX/modeling/misfit.py

@@ -15,6 +15,15 @@
 
 
 class NonGaussianContinuousMisfit(object):
+    """
+    Abstract class to model the misfit component of the cost functional.
+    In the following :code:`x` will denote the variable :code:`[u, m, p]`, denoting respectively
+    the state :code:`u`, the parameter :code:`m`, and the adjoint variable :code:`p`.
+
+    The methods in the class misfit will usually access the state u and possibly the
+    parameter :code:`m`. The adjoint variables will never be accessed.
+    """
+
     def __init__(self, Vh: list, form, bc0=[]):
         self.Vh = Vh
         self.form = form
@@ -30,6 +39,10 @@ def __init__(self, Vh: list, form, bc0=[]):
         self.xfun = [dlx.fem.Function(Vhi) for Vhi in Vh]
 
     def cost(self, x: list) -> float:
+        """
+        Given x evaluate the cost functional.
+        Only the state u and (possibly) the parameter m are accessed.
+        """
         hpx.updateFromVector(self.xfun[hpx.STATE], x[hpx.STATE])
         u_fun = self.xfun[hpx.STATE]
 
@@ -41,6 +54,10 @@ def cost(self, x: list) -> float:
         return self.Vh[hpx.STATE].mesh.comm.allreduce(glb_cost_proc, op=MPI.SUM)
 
     def grad(self, i: int, x: list, out: dlx.la.Vector) -> None:
+        """
+        Given the state and the paramter in :code:`x`, compute the partial gradient of the misfit
+        functional in with respect to the state (:code:`i == STATE`) or with respect to the parameter (:code:`i == PARAMETER`).
+        """
         hpx.updateFromVector(self.xfun[hpx.STATE], x[hpx.STATE])
         u_fun = self.xfun[hpx.STATE]
 
@@ -69,6 +86,9 @@ def setLinearizationPoint(self, x: list, gauss_newton_approx=False) -> None:
         self.gauss_newton_approx = gauss_newton_approx
 
     def apply_ij(self, i: int, j: int, dir: dlx.la.Vector, out: dlx.la.Vector) -> None:
+        """
+        Apply the second variation :math:`\delta_{ij}` (:code:`i,j = STATE,PARAMETER`) of the cost in direction :code:`dir`.
+        """
         form = self.form(*self.x_lin_fun)
         if j == hpx.STATE:
             dlx.fem.set_bc(dir.array, self.bc0)

hippylibX/modeling/model.py

@@ -315,7 +315,7 @@ def Rsolver(self) -> Any:
         operator :math:`R`.
 
         The solver object should implement the method :code:`Rsolver.solve(z,r)` such that
-        :math:`Rz \approx r`.
+        :math: `Rz approx r`
         """
         return self.prior.Rsolver
 

hippylibX/modeling/prior.py

@@ -81,6 +81,12 @@ def generate_vector(self) -> petsc4py.PETSc.Vec:
 
 
 class SqrtPrecisionPDE_Prior:
+    """
+    This class implement a prior model with covariance matrix
+    :math:`C = A^{-1} M A^-1`,
+    where A is the finite element matrix arising from discretization of sqrt_precision_varf_handler
+    """
+
     def __init__(
         self,
         Vh: dlx.fem.FunctionSpace,
@@ -90,16 +96,12 @@ def __init__(
         """
         Construct the prior model.
         Input:
+
         - :code:`Vh`:              the finite element space for the parameter
         - :code:sqrt_precision_varf_handler: the PDE representation of the  sqrt of the covariance operator
         - :code:`mean`:            the prior mean
         """
 
-        """
-        This class implement a prior model with covariance matrix
-        :math:`C = A^{-1} M A^-1`,
-        where A is the finite element matrix arising from discretization of sqrt_precision_varf_handler
-        """
         self.dx = ufl.Measure("dx", metadata={"quadrature_degree": 4})
         self.ds = ufl.Measure("ds", metadata={"quadrature_degree": 4})
 
@@ -202,7 +204,8 @@ def R(self) -> petsc4py.PETSc.Mat:
 
     def generate_parameter(self, dim: int) -> dlx.la.Vector:
         """
-        Inizialize a vector :code:`x` to be compatible with the range/domain of :math:`R`.
+        Initialize a vector :code:`x` to be compatible with the range/domain of :math:`R`.
+
         If :code:`dim == "noise"` inizialize :code:`x` to be compatible with the size of
         white noise used for sampling.
         """

hippylibX/modeling/reducedHessian.py

@@ -23,6 +23,7 @@ class ReducedHessian:
     """
     This class implements matrix free application of the reduced Hessian operator.
     The constructor takes the following parameters:
+
     - :code:`model`:               the object which contains the description of the problem.
     - :code:`misfit_only`:         a boolean flag that describes whenever the full Hessian or only the misfit component of the Hessian is used.
 

hippylibX/utils/projection.py

@@ -13,12 +13,17 @@
 
 
 def projection(v, target_func, bcs=[]):
-    # reference:
-    # https://github.com/michalhabera/dolfiny/blob/master/dolfiny/projection.py
+    """
+    Return projection of given expression :code:`v` onto the finite element
+    space of a function :code:`target_func`.
 
-    # v -> expression to project
-    # target_func -> function that contains the projection
+    reference:
+    https://github.com/michalhabera/dolfiny/blob/master/dolfiny/projection.py
 
+    Inputs:
+        - :code:`v`:  expression to project
+        - :code:`target_func` : function that contains the projection
+    """
     V = target_func.function_space
     dx = ufl.dx(V.mesh)
     w = ufl.TestFunction(V)