Merge pull request #162 from ComputationalCryoEM/develop

Update functions on block diagonal matrix and basis class for polar 2D

Author: junchaoxia

CONTRIBUTORS.rst

@@ -0,0 +1,53 @@
+ASPIRE has contributions from a growing team of collaborators.
+To learn more about them please visit our `Team webpage <http://spr.math.princeton.edu/team>`_.
+
+Developers of the Python and Matlab collection of ASPIRE codes are listed below.
+
+
+|
+
+
+.. table:: Python Contributors
+
+   +------------------+--------------+---------------------------+----------------------+
+   | Name             | GithubID     | Email                     | Affiliations         |
+   +==================+==============+===========================+======================+
+   | Joakim Andén     | janden       | [email protected]             | Princeton University |
+   +------------------+--------------+---------------------------+----------------------+
+   | Vineet Bansal    | vineetbansal | [email protected]     | Princeton University |
+   +------------------+--------------+---------------------------+----------------------+
+   | Robbie Brook     | rbrook       | [email protected]          | Princeton University |
+   +------------------+--------------+---------------------------+----------------------+
+   | Ayelet Heimowitz | ayeltg       | [email protected]     | Princeton University |
+   +------------------+--------------+---------------------------+----------------------+
+   | Gabi Pragier     | pgabi        | [email protected]     | Tel Aviv University  |
+   +------------------+--------------+---------------------------+----------------------+
+   | Itay Sason       | itaysason    | [email protected]  | Tel Aviv University  |
+   +------------------+--------------+---------------------------+----------------------+
+   | Yoel Shkolnisky  | yoelsh       | [email protected]    | Tel Aviv University  |
+   +------------------+--------------+---------------------------+----------------------+
+   | Amit Singer      | amitsinger   | [email protected]       | Princeton University |
+   +------------------+--------------+---------------------------+----------------------+
+   | Garrett Wright   | garrettwrong | [email protected]    | Princeton University |
+   +------------------+--------------+---------------------------+----------------------+
+   | Junchao Xia      | junchaoxia   | [email protected] | Princeton University |
+   +------------------+--------------+---------------------------+----------------------+
+
+|
+
+
+.. table:: Additional MATLAB Contributors not listed above:
+
+
+   +-------------+----------+-------+----------------------+
+   | Name        | GithubID | Email | Affiliations         |
+   +=============+==========+=======+======================+
+   | Mar Cohen   |          |       | Tel Aviv University  |
+   +-------------+----------+-------+----------------------+
+   | Lanhui Wang |          |       | Princeton University |
+   +-------------+----------+-------+----------------------+
+   | Jane Zhao   |          |       | Princeton University |
+   +-------------+----------+-------+----------------------+
+
+
+If we've missed you or would like to update your details, please contact us and we'll be happy to do so.

azure-pipelines.yml

@@ -1,3 +1,12 @@
+schedules:
+- cron: "0 0 * * *" # cron syntax defining a schedule
+  displayName: Daily Build
+  branches:
+    include:
+    - master
+    - develop
+  always: true # Always run the Daily Build
+
 jobs:
 - job:
   displayName: ubuntu-latest
@@ -36,7 +45,7 @@ jobs:
 
 #  - bash: conda update -q -y conda
 #    displayName: Update conda
-    
+
   - bash: conda env create --name myEnvironment --file environment.yml
     displayName: Create Anaconda environment
 
@@ -96,10 +105,10 @@ jobs:
 
   - bash: conda clean -i -t -y
     displayName: Removing conda cached package tarballs.
-    
+
   - bash: conda env create --name myEnvironment --file environment.yml
     displayName: Create Anaconda environment
-    
+
   - bash: |
       source activate myEnvironment
       python setup.py install
@@ -149,7 +158,7 @@ jobs:
 
   - script: conda update -q -y conda
     displayName: Update conda
-    
+
   - script: |
       call activate myEnvironment
       python setup.py install

docs/source/authors.rst

@@ -0,0 +1,4 @@
+Authors
+=======
+
+.. include:: ../../CONTRIBUTORS.rst

docs/source/index.rst

@@ -17,6 +17,7 @@ linear and nonlinear dimensionality reduction, randomized algorithms in numerica
    installation
    quickstart
    modules
+   authors
 
 
 Indices and tables
@@ -30,4 +31,4 @@ Indices and tables
 References
 ==========
 
-.. bibliography:: references.bib
+.. bibliography:: references.bib

src/aspire/basis/__init__.py

@@ -15,7 +15,7 @@ class Basis:
     Define a base class for expanding 2D particle images and 3D structure volumes
 
     """
-    def __init__(self, size, ell_max=None):
+    def __init__(self, size, ell_max=None, dtype=np.float64):
         """
         Initialize an object for the base of basis class
 
@@ -36,6 +36,9 @@ def __init__(self, size, ell_max=None):
         self.count = 0
         self.ell_max = ell_max
         self.ndim = ndim
+        self.dtype = dtype
+        if self.dtype != np.float64:
+            raise NotImplementedError("Currently only implemented for default double (np.float64) type")
 
         self._build()
 
@@ -149,46 +152,46 @@ def mat_evaluate_t(self, X):
         """
         return mdim_mat_fun_conj(X, len(self.sz), 1, self.evaluate_t)
 
-    def expand(self, v):
+    def expand(self, x):
         """
-        Expand array in basis
+        Obtain coefficients in the basis from those in standard coordinate basis
 
-        This is a similar function to `evaluate_t` but with more accuracy by
-         using the cg optimizing of linear equation, Ax=b.
+        This is a similar function to evaluate_t but with more accuracy by using
+        the cg optimizing of linear equation, Ax=b.
 
-        If `v` is a matrix of size `basis.ct`-by-..., `B` is the change-of-basis
-        matrix of this basis, and `x` is a matrix of size `self.sz`-by-...,
-        the function calculates  v = (B' * B)^(-1) * B' * x, where the rows
-        of `B` and columns of `x` are read as vectorized arrays.
-
-        :param v: An array whose first few dimensions are to be expanded in this basis.
-            These dimensions must equal `self.sz`.
-        :return: The coefficients of `v` expanded in this basis. If more than
-            one array of size `self.sz` is found in `v`, the second and higher
-            dimensions of the return value correspond to those higher dimensions of `v`.
+        :param x: An array whose first two or three dimensions are to be expanded
+            the desired basis. These dimensions must equal `self.sz`.
+        :return : The coefficients of `v` expanded in the desired basis.
+            The first dimension of `v` is with size of `count` and the
+            second and higher dimensions of the return value correspond to
+            those higher dimensions of `x`.
 
         """
-        ensure(v.shape[:self.ndim] == self.sz, f'First {self.ndim} dimensions of v must match {self.sz}.')
+        # ensure the first dimensions with size of self.sz
+        x, sz_roll = unroll_dim(x, self.ndim + 1)
+        ensure(x.shape[:self.ndim] == self.sz,
+               f'First {self.ndim} dimensions of x must match {self.sz}.')
 
-        v, sz_roll = unroll_dim(v, self.ndim + 1)
-        b = self.evaluate_t(v)
-        operator = LinearOperator(
-            shape=(self.count, self.count),
-            matvec=lambda x: self.evaluate_t(self.evaluate(x))
-        )
+        operator = LinearOperator(shape=(self.count, self.count),
+                                  matvec=lambda v: self.evaluate_t(self.evaluate(v)))
 
         # TODO: (from MATLAB implementation) - Check that this tolerance make sense for multiple columns in v
-        tol = 10 * np.finfo(v.dtype).eps
+        tol = 10*np.finfo(x.dtype).eps
         logger.info('Expanding array in basis')
-        v, info = cg(operator, b, tol=tol)
 
-        v = v[..., np.newaxis]
+        # number of image samples
+        n_data = np.size(x, self.ndim)
+        v = np.zeros((self.count, n_data), dtype=x.dtype)
 
-        if info != 0:
-            raise RuntimeError('Unable to converge!')
+        for isample in range(0, n_data):
+            b = self.evaluate_t(x[..., isample])
+            # TODO: need check the initial condition x0 can improve the results or not.
+            v[..., isample], info = cg(operator, b, tol=tol)
+            if info != 0:
+                raise RuntimeError('Unable to converge!')
 
+        # return v coefficients with the first dimension of self.count
         v = roll_dim(v, sz_roll)
-
         return v
 
     def expand_t(self, v):

src/aspire/basis/dirac.py

@@ -1,5 +1,74 @@
+import logging
+import numpy as np
+
 from aspire.basis import Basis
+from aspire.utils.matrix import roll_dim, unroll_dim
+from aspire.utils.matlab_compat import m_flatten, m_reshape
+
+logger = logging.getLogger(__name__)
 
 
 class DiracBasis(Basis):
-    pass
+    """
+    Define a derived class for Dirac basis
+    """
+    def __init__(self, sz, mask=None):
+        """
+        Initialize an object for Dirac basis
+        :param sz: The shape of the vectors for which to define the basis.
+        :param mask: A boolean _mask of size sz indicating which coordinates
+            to include in the basis (default np.full(sz, True)).
+        """
+        if mask is None:
+            mask = np.full(sz, True)
+        self._mask = m_flatten(mask)
+
+        super().__init__(sz)
+
+    def _build(self):
+        """
+        Build the internal data structure to Dirac basis
+        """
+        logger.info('Expanding object in a Dirac basis.')
+        self.count = np.sum(self._mask)
+        self._sz_prod = self.nres ** self.ndim
+
+    def evaluate(self, v):
+        """
+        Evaluate coefficients in standard coordinate basis from those in Dirac basis
+
+        :param v: A coefficient vector (or an array of coefficient vectors) to
+            be evaluated. The first dimension must equal `self.count`.
+        :return: The evaluation of the coefficient vector(s) `v` for this basis.
+            This is an array whose first dimensions equal `self.sz` and the remaining
+            dimensions correspond to dimensions two and higher of `v`.
+        """
+        v, sz_roll = unroll_dim(v, 2)
+        x = np.zeros(shape=(self._sz_prod,) + v.shape[1:])
+        x[self._mask, ...] = v
+        x = m_reshape(x, self.sz + x.shape[1:])
+        x = roll_dim(x, sz_roll)
+
+        return x
+
+    def evaluate_t(self, x):
+        """
+        Evaluate coefficient in Dirac basis from those in standard coordinate basis
+
+        :param x: The coefficient array to be evaluated. The first dimensions
+            must equal `self.sz`.
+        :return: The evaluation of the coefficient array `v` in the dual basis
+            of `basis`. This is an array of vectors whose first dimension equals
+             `self.count` and whose remaining dimensions correspond to
+             higher dimensions of `v`.
+        """
+        x, sz_roll = unroll_dim(x, self.ndim + 1)
+        x = m_reshape(x, new_shape=(self._sz_prod,) + x.shape[self.ndim:])
+        v = np.zeros(shape=(self.count, ) + x.shape[1:])
+        v = x[self._mask, ...]
+        v = roll_dim(v, sz_roll)
+
+        return v
+
+    def expand(self, x):
+        return self.evaluate_t(x)

src/aspire/basis/ffb_2d.py

@@ -3,11 +3,8 @@
 from numpy import pi
 from scipy.special import jv
 from scipy.fftpack import ifft, fft
-from scipy.sparse.linalg import LinearOperator, cg
 
 from aspire.nfft import anufft3, nufft3
-
-from aspire.utils import ensure
 from aspire.utils.matrix import roll_dim, unroll_dim
 from aspire.utils.matlab_compat import m_reshape
 from aspire.basis.basis_utils import lgwt
@@ -246,45 +243,3 @@ def evaluate_t(self, x):
         # return v coefficients with the first dimension of self.count
         v = roll_dim(v, sz_roll)
         return v
-
-    def expand(self, x):
-        """
-        Obtain coefficients in FB basis from those in standard 2D coordinate basis
-
-        This is a similar function to evaluate_t but with more accuracy by using
-        the cg optimizing of linear equation, Ax=b.
-
-        :param x: An array whose first two dimensions are to be expanded in FB basis.
-             These dimensions must equal `self.sz`.
-        :return : The coefficients of `v` expanded in FB basis. The first dimension
-            of `v` is with size of `count` and the second and higher dimensions
-            of the return value correspond to those higher dimensions of `x`.
-
-        """
-        # ensure the first two dimensions with size of self.sz
-        x, sz_roll = unroll_dim(x, self.ndim + 1)
-        x = m_reshape(x, (self.sz[0], self.sz[1], -1))
-        ensure(x.shape[:self.ndim] == self.sz,
-               f'First {self.ndim} dimensions of x must match {self.sz}.')
-
-        operator = LinearOperator(shape=(self.count, self.count),
-                                  matvec=lambda v: self.evaluate_t(self.evaluate(v)))
-
-        # TODO: (from MATLAB implementation) - Check that this tolerance make sense for multiple columns in v
-        tol = 10*np.finfo(x.dtype).eps
-        logger.info('Expanding array in basis')
-
-        # number of image samples
-        n_data = np.size(x, self.ndim)
-        v = np.zeros((self.count, n_data), dtype=x.dtype)
-
-        for isample in range(0, n_data):
-            b = self.evaluate_t(x[..., isample])
-            # TODO: need check the initial condition x0 can improve the results or not.
-            v[..., isample], info = cg(operator, b, tol=tol)
-            if info != 0:
-                raise RuntimeError('Unable to converge!')
-
-        # return v coefficients with the first dimension of self.count
-        v = roll_dim(v, sz_roll)
-        return v