Merge pull request #4 from schroedk/develop

Develop
schroedk · Oct 25, 2020 · 0298ab6 · 0298ab6
2 parents 956dfaa + b6488cd
commit 0298ab6
Show file tree

Hide file tree

Showing 8 changed files with 157 additions and 78 deletions.
diff --git a/.github/workflows/tox.yaml b/.github/workflows/tox.yaml
@@ -1,10 +1,10 @@
-name: Merge master, run tests and build docu
+name: Run tests and build docu
 
 on:
   push:
-    branches: [master]
+    branches: [main]
   pull_request:
-    branches: [master]
+    branches: [main]
 
 jobs:
   build:
@@ -14,10 +14,6 @@ jobs:
       - uses: actions/checkout@v2
         with:
           fetch-depth: 0
-      - name: Merge master into current branch
-        run: |
-          git fetch origin master:master --update-head-ok
-          git merge master
       - name: Setup Python 3.8
         uses: actions/setup-python@v1
         with:

diff --git a/README.md b/README.md
@@ -7,7 +7,7 @@ This repository contains python code to implement a basic variant of the Harmoni
 > Available online: https://jmlr.org/papers/volume19/17-244/17-244.pdf
 
 ## Version history
-* Version 0.0.1, 10/21/2020
+* Version 0.0.1, 10/25/2020
 
 ## Author
-Kristof Schröder
+Kristof Schröder
diff --git a/docs/conf.py b/docs/conf.py
@@ -27,8 +27,19 @@
 
 # Add any Sphinx extension module names here, as strings. They can be extensions
 # coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
-extensions = ['sphinx.ext.napoleon', 'sphinx.ext.autodoc', 'sphinx.ext.doctest', 'sphinx.ext.linkcode', 'sphinx_rtd_theme']
+extensions = ['sphinx.ext.napoleon',
+              'sphinx.ext.doctest', 'sphinx.ext.intersphinx','sphinx.ext.linkcode', 'sphinx_rtd_theme', 'sphinx.ext.autodoc', 'sphinx.ext.autodoc.typehints'
+              ]
 
+autodoc_typehints = 'description'
+
+intersphinx_mapping = {
+    'numpy': ('https://numpy.org/doc/stable/', None),
+    'scipy': ('https://docs.scipy.org/doc/scipy/reference', None),
+}
+
+napoleon_use_param = True
+typehints_fully_qualified = True
 
 # adding links to source files
 # see https://www.sphinx-doc.org/en/master/usage/extensions/linkcode.html#module-sphinx.ext.linkcode
@@ -40,7 +51,8 @@ def linkcode_resolve(domain, info):
 
     path, linkExtension = getPathAndLinkExtension(info['module'])
     objectName = info['fullname']
-    if "." in objectName:  # don't add source link to methods within classes (we might want to change that in the future)
+    # don't add source link to methods within classes (we might want to change that in the future)
+    if "." in objectName:
         return None
     lineno = findLineFromObjectName(path, objectName)
 
@@ -67,7 +79,8 @@ def getPathAndLinkExtension(module: str):
         linkExtension = f"src/{filename}/__init__.py"
         return os.path.join(sourcePathPrefix, "__init__.py"), linkExtension
     else:
-        raise Exception(f"{sourcePathPrefix} is neither a module nor a package with init - did you fortet to add an __init__.py?")
+        raise Exception(
+            f"{sourcePathPrefix} is neither a module nor a package with init - did you fortet to add an __init__.py?")
 
 
 def findLineFromObjectName(sourceFile, objectName):
@@ -81,6 +94,7 @@ def findLineFromObjectName(sourceFile, objectName):
     else:
         return desiredNode.lineno
 
+
 autodoc_default_options = {
     'exclude-members': 'log, DuplicateColumnNamesException',
     'member-order': 'bysource',

diff --git a/docs/requirements.txt b/docs/requirements.txt
@@ -1,2 +1,3 @@
 numpy
-scipy
+scipy
+sphinx_rtd_theme
diff --git a/scripts/test_run.py → scripts/simple_example.py b/scripts/test_run.py → scripts/simple_example.py
@@ -3,6 +3,16 @@
 from hmirls.operators import SamplingMatrixOperator
 from hmirls.problem import Problem
 
+"""
+Simple example for recovering a low-rank matrix via minimization
+
+    .. math::
+
+        \\min_{x \\in \\mathbb{C}^{d_1 \\times d_2}} \\operatorname{rank}(x), \\text{s.t.} \\Phi(x) = y,
+
+
+"""
+
 u = np.array([1.0, 10.0, -2.0, 0.1]).reshape((-1, 1))
 v = np.array([1.0, 2.0, 3.0, 4.0]).reshape((-1, 1))
 X = np.matmul(u, v.transpose())

diff --git a/src/hmirls/operators.py b/src/hmirls/operators.py
@@ -1,8 +1,8 @@
 from enum import Enum
-from typing import Tuple, List, Union, Callable
+from typing import Tuple, List, Union, Callable, Optional
 import numpy as np
 from scipy.sparse import csr_matrix, kron, eye
-from scipy.sparse.linalg import LinearOperator as ScipyLinearOperator, aslinearoperator
+from scipy.sparse.linalg import LinearOperator, aslinearoperator
 
 
 class MatrixOperatorCompatibility:
@@ -64,26 +64,42 @@ def _same_indexing_order(
 
 class MatrixOperator:
     class IndexingOrder(Enum):
+        """
+        Enum for order argument for :func:`numpy.reshape` and :func:`numpy.flatten`
+        """
         ROW_MAJOR = "C"
         COLUMN_MAJOR = "F"
 
     def __init__(
         self,
-        flattened_operator: Union[ScipyLinearOperator, np.ndarray, csr_matrix],
+        flattened_operator: Union[LinearOperator, np.ndarray, csr_matrix],
         input_shape: Tuple[int, int],
         output_shape: Tuple[int, int],
-        representing_matrix: Union[np.ndarray, csr_matrix] = None,
-        order=IndexingOrder,
+        order: IndexingOrder,
+        representing_matrix: Optional[Union[np.ndarray, csr_matrix]] = None,
     ):
         """
-        # ToDO More detailed description, doc tests
-        Takes care on flattening and reshaping.
-        Wraps a scipy linear operator or numpy array.
+        Represents an operator
+
+        .. math::
+
+            \\Phi: \\mathbb{C}^{d_1 \\times d_2} \\rightarrow \\mathbb{C}^{m_1 \\times m_2}.
+
+        Takes care on flattening and reshaping. Wraps a scipy linear operator or numpy array.
+
+        :param flattened_operator: Represents the flattened map
+
+            .. math::
+
+                \\varphi: \\mathbb{C}^{d_1 \\cdot d_2} \\rightarrow \\mathbb{C}^{m_1 \\cdot m_2},
+
+            such that :math:`\\varphi(\\operatorname{vec}(x)) = \\operatorname{vec}(\\Phi(x))`.
+
+        :param input_shape: :math:`(d_1,d_2)`
+        :param output_shape: :math:`(m_1,m_2)`
+        :param order: order for vectorization :math:`\\operatorname{vec}`
+        :param representing_matrix: matrix representation of flattened operator
 
-        :type output_shape:
-        :param flattened_operator:
-        :param input_shape:
-        :param order:
         """
         self._output_shape = output_shape
         self._input_shape = input_shape
@@ -109,8 +125,8 @@ def __add__(self, other: "MatrixOperator"):
                     self._flattened_operator + other._flattened_operator,
                     self._input_shape,
                     self._output_shape,
-                    representing_matrix,
                     self._order,
+                    representing_matrix,
                 )
             raise ValueError(f"Operators are not compatible")
         return NotImplemented
@@ -126,11 +142,11 @@ def dot(self, x):
 
         Parameters
         ----------
-        x : array_like or MatrixOperator or scalar
+        x : array_like of shape :math:`(d_1, d_2)` or MatrixOperator compatible for composition or scalar
 
         Returns
         -------
-        Ax : array or MatrixOperator that represents
+            :math:`\\Phi(x)` : array or MatrixOperator that represents
             the result of applying this linear operator on x.
 
         """
@@ -146,18 +162,18 @@ def dot(self, x):
                     )
                 return MatrixOperator(
                     self.flattened_operator.dot(x.flattened_operator),
-                    input_shape=x.input_shape,
-                    output_shape=self.output_shape,
-                    representing_matrix=representing_matrix,
-                    order=self.order,
+                    x.input_shape,
+                    self.output_shape,
+                    self.order,
+                    representing_matrix,
                 )
         elif np.isscalar(x):
             return MatrixOperator._from_callable(self, lambda op: x * op)
         else:
             if x.shape == self.input_shape:
-                x_flattened = x.flatten(self.order)
+                x_flattened = x.flatten(order=self.order.value)
                 return self.flattened_operator(x_flattened).reshape(
-                    self.output_shape, order=self.order
+                    self.output_shape, order=self.order.value
                 )
             else:
                 raise ValueError(
@@ -207,8 +223,8 @@ def _from_callable(
         cls,
         matrix_operator: "MatrixOperator",
         fcn: Callable[
-            [Union[ScipyLinearOperator, np.ndarray]],
-            Union[ScipyLinearOperator, np.ndarray],
+            [Union[LinearOperator, np.ndarray]],
+            Union[LinearOperator, np.ndarray],
         ],
     ):
         """
@@ -220,13 +236,12 @@ def _from_callable(
         representing_matrix = None
         if matrix_operator.representing_matrix is not None:
             representing_matrix = fcn(matrix_operator.representing_matrix)
-        return cls(
-            fcn(matrix_operator.flattened_operator),
-            input_shape=matrix_operator.input_shape,
-            output_shape=matrix_operator.output_shape,
-            representing_matrix=representing_matrix,
-            order=matrix_operator.order,
-        )
+        return cls(fcn(matrix_operator.flattened_operator),
+                   matrix_operator.input_shape,
+                   matrix_operator.output_shape,
+                   matrix_operator.order,
+                   representing_matrix,
+                   )
 
     @property
     def output_shape(self):
@@ -257,7 +272,7 @@ def adjoint(self):
 
         Returns
         -------
-        A_H : :class:`~MatrixOperator`
+        :math:`\\Phi^{\\star}` : :class:`~MatrixOperator`
             Hermitian adjoint of self.
         """
         return self._adjoint()
@@ -272,7 +287,7 @@ def transpose(self):
 
         Returns
         -------
-        A_T : :class:`~MatrixOperator`
+        :math:`\\Phi^{\\star}` : :class:`~MatrixOperator`
              adjoint of self.
         """
 
@@ -284,31 +299,22 @@ def _adjoint(self):
         """Default implementation of _adjoint"""
         representing_matrix = None
         if self.representing_matrix is not None:
-            representing_matrix = self._representing_matrix.H
+            representing_matrix = self._representing_matrix.conjugate().transpose()
         return MatrixOperator(
             self._flattened_operator.H,
             self._output_shape,
             self._input_shape,
-            representing_matrix,
             self.order,
+            representing_matrix,
         )
 
     def _transpose(self):
         """ Default implementation of _transpose"""
-        representing_matrix = None
-        if self.representing_matrix is not None:
-            representing_matrix = self._representing_matrix.T
-        return MatrixOperator(
-            self._flattened_operator.T,
-            self._output_shape,
-            self._input_shape,
-            representing_matrix,
-            self.order,
-        )
+        return MatrixOperator._from_callable(self, lambda op: op.T)
 
 
-class SamplingOperator(ScipyLinearOperator):
-    def __init__(self, indices: List[int], input_dimension: int):
+class SamplingOperator(LinearOperator):
+    def __init__(self, indices: Union[List[int], np.ndarray], input_dimension: int):
         """
 
         :type input_dimension:
@@ -349,11 +355,12 @@ def _construct_sampling_matrix(self):
 
     @classmethod
     def from_matrix_indices(
-        cls, row_indices: List[int], column_indices: List[int], shape, order="F"
+        cls, row_indices: List[int], column_indices: List[int], shape: Tuple[int, int],
+            order: MatrixOperator.IndexingOrder = MatrixOperator.IndexingOrder.COLUMN_MAJOR
     ):
         return cls(
-            np.ravel_multi_index((row_indices, column_indices), shape, order=order),
-            input_dimension=np.prod(shape),
+            np.ravel_multi_index((row_indices, column_indices), shape, order=order.value),
+            input_dimension=np.product(shape),
         )
 
 
@@ -363,17 +370,17 @@ def __init__(
         row_indices: List[int],
         column_indices: List[int],
         shape: Tuple[int, int],
-        order="F",
+        order: MatrixOperator.IndexingOrder = MatrixOperator.IndexingOrder.COLUMN_MAJOR,
     ):
         flattened_operator = SamplingOperator.from_matrix_indices(
             row_indices, column_indices, shape, order
         )
-        super(SamplingMatrixOperator, self).__init__(
+        super().__init__(
             flattened_operator,
             shape,
             (len(row_indices), 1),
-            flattened_operator.sampling_matrix,
             order,
+            flattened_operator.sampling_matrix,
         )
 
 
@@ -401,6 +408,7 @@ def __init__(
             flattened_operator,
             shape,
             shape,
-            representing_matrix=weight_matrix,
-            order=order,
+            order,
+            weight_matrix,
+
         )