Merge branch 'master' of github.com:prody/ProDy into devel-5

prody/__init__.py

@@ -60,13 +60,11 @@ def turnonDepracationWarnings(action='always'):
 
 from . import utilities
 from .utilities import *
-from .utilities import PackageLogger, PackageSettings
+from .utilities import LOGGER, PackageSettings
 from .utilities import getPackagePath, joinRepr, tabulate
 __all__.extend(utilities.__all__)
 __all__.append('utilities')
 
-LOGGER = PackageLogger('.prody')
-
 SETTINGS = PackageSettings('prody', logger=LOGGER)
 SETTINGS.load()
 

prody/apps/apptools.py

@@ -1,4 +1,4 @@
-"""This module define tools for application development."""
+"""This module defines tools for application development."""
 
 from copy import copy
 try:

prody/atomic/functions.py

@@ -339,6 +339,8 @@ def extendAtoms(nodes, atoms, is3d=False):
     #LOGGER.progress('Extending atoms...', n_nodes, '_prody_extendAtoms_extend')
     for i, node in enumerate(i_nodes):
         #LOGGER.update(i, label='_prody_extendAtoms')
+        if node is None:
+            continue
         res = get(node.getChid() or None, node.getResnum(),
                   node.getIcode() or None, node.getSegname() or None)
         if res is None:

prody/chromatin/cluster.py

@@ -27,7 +27,7 @@ def KMeans(V, **kwargs):
     for details.
 
     :arg V: row-normalized eigenvectors for the purpose of clustering.
-    :type V: :class:`numpy.ndarray`
+    :type V: :class:`~numpy.ndarray`
 
     :arg n_clusters: specifies the number of clusters. 
     :type n_clusters: int
@@ -54,7 +54,7 @@ def Hierarchy(V, **kwargs):
     explaination of the arguments. Here lists arguments that are different from those of scipy.
 
     :arg V: row-normalized eigenvectors for the purpose of clustering.
-    :type V: :class:`numpy.ndarray`
+    :type V: :class:`~numpy.ndarray`
 
     :arg inconsistent_percentile: if the clustering *criterion* for :func:`scipy.cluster.hierarchy.fcluster`
     is ``inconsistent`` and threshold *t* is not given (default), then the function will use the percentile specified 
@@ -207,7 +207,7 @@ def GaussianMixture(V, **kwargs):
     for details.
 
     :arg V: row-normalized eigenvectors for the purpose of clustering.
-    :type V: :class:`numpy.ndarray`
+    :type V: :class:`~numpy.ndarray`
 
     :arg n_clusters: specifies the number of clusters. 
     :type n_clusters: int
@@ -236,7 +236,7 @@ def BayesianGaussianMixture(V, **kwargs):
     for details.
 
     :arg V: row-normalized eigenvectors for the purpose of clustering.
-    :type V: :class:`numpy.ndarray`
+    :type V: :class:`~numpy.ndarray`
 
     :arg n_clusters: specifies the number of clusters. 
     :type n_clusters: int
@@ -264,7 +264,7 @@ def showLinkage(V, **kwargs):
     """Shows the dendrogram of hierarchical clustering on *V*. See :func:`scipy.cluster.hierarchy.dendrogram` for details.
 
     :arg V: row-normalized eigenvectors for the purpose of clustering.
-    :type V: :class:`numpy.ndarray`
+    :type V: :class:`~numpy.ndarray`
 
     """
 

prody/dynamics/anm.py

@@ -8,11 +8,11 @@
 from prody import LOGGER
 from prody.atomic import Atomic, AtomGroup
 from prody.proteins import parsePDB
-from prody.utilities import checkCoords
+from prody.utilities import checkCoords, solveEig
 from prody.kdtree import KDTree
 
 from .nma import NMA, MaskedNMA
-from .gnm import GNMBase, solveEig, checkENMParameters
+from .gnm import GNMBase, checkENMParameters
 
 __all__ = ['ANM', 'MaskedANM', 'calcANM']
 

prody/dynamics/bbenm.py

@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
-"""This module defines a class and a function for rotating translating blocks
-(RTB) calculations."""
+"""This module defines a class and a function for using the tensorial network model, 
+which includes bond-bending and twist elasticities."""
 
 import numpy as np
 

prody/dynamics/functions.py

@@ -175,7 +175,7 @@ def loadModel(filename, **kwargs):
             if attr in ('type', '_name', '_title'):
                 continue
             elif attr in ('_trace', '_cutoff', '_gamma'):
-                dict_[attr] = float(attr_dict[attr])
+                dict_[attr] = attr_dict[attr][()]
             elif attr in ('_dof', '_n_atoms', '_n_modes'):
                 dict_[attr] = int(attr_dict[attr])
             elif attr in ('masked', ):

prody/dynamics/gnm.py

@@ -11,126 +11,12 @@
 from prody.atomic import Atomic, AtomGroup
 from prody.proteins import parsePDB
 from prody.kdtree import KDTree
-from prody.utilities import importLA, checkCoords, div0
+from prody.utilities import importLA, checkCoords, solveEig, ZERO
 
 from .nma import NMA, MaskedNMA
 from .gamma import Gamma
 
-__all__ = ['GNM', 'MaskedGNM', 'solveEig', 'calcGNM']
-
-ZERO = 1e-6
-
-
-def solveEig(M, n_modes=None, zeros=False, turbo=True, is3d=False):
-    linalg = importLA()
-    dof = M.shape[0]
-
-    expct_n_zeros = 6 if is3d else 1
-
-    if n_modes is None:
-        eigvals = None
-        n_modes = dof
-    else:
-        if n_modes >= dof:
-            eigvals = None
-            n_modes = dof
-        else:
-            eigvals = (0, n_modes+expct_n_zeros-1)
-
-    def _eigh(M, eigvals=None, turbo=True):
-        if linalg.__package__.startswith('scipy'):
-            from scipy.sparse import issparse
-
-            if eigvals:
-                turbo = False
-            if not issparse(M):
-                values, vectors = linalg.eigh(M, turbo=turbo, eigvals=eigvals)
-            else:
-                try:
-                    from scipy.sparse import linalg as scipy_sparse_la
-                except ImportError:
-                    raise ImportError('failed to import scipy.sparse.linalg, '
-                                      'which is required for sparse matrix '
-                                      'decomposition')
-                if eigvals:
-                    j = eigvals[0]
-                    k = eigvals[-1] + 1
-                else:
-                    j = 0
-                    k = dof
-
-                if k >= dof:
-                    k -= 1
-                    LOGGER.warning('Cannot calculate all eigenvalues for sparse matrices, thus '
-                                   'the last eigenvalue is omitted. See scipy.sparse.linalg.eigsh '
-                                   'for more information')
-                values, vectors = scipy_sparse_la.eigsh(M, k=k, which='SA')
-                values = values[j:k]
-                vectors = vectors[:, j:k]
-        else:
-            if n_modes is not None:
-                LOGGER.info('Scipy is not found, all modes were calculated.')
-            else:
-                n_modes = dof
-            values, vectors = linalg.eigh(M)
-        return values, vectors
-
-    def _calc_n_zero_modes(M):
-        from scipy.sparse import issparse
-
-        if not issparse(M):
-            w = linalg.eigvalsh(M)
-        else:
-            try:
-                from scipy.sparse import linalg as scipy_sparse_la
-            except ImportError:
-                raise ImportError('failed to import scipy.sparse.linalg, '
-                                    'which is required for sparse matrix '
-                                    'decomposition')
-            w, _ = scipy_sparse_la.eigsh(M, k=dof-1, which='SA')
-        n_zeros = sum(w < ZERO)
-        return n_zeros
-
-    values, vectors = _eigh(M, eigvals, turbo)
-    n_zeros = sum(values < ZERO)
-
-    if n_zeros < n_modes + expct_n_zeros:
-        if n_zeros < expct_n_zeros:
-            LOGGER.warning('Fewer than %d (%d) zero eigenvalues were calculated.'%(expct_n_zeros, n_zeros))
-        elif n_zeros > expct_n_zeros:
-            LOGGER.warning('More than %d (%d) zero eigenvalues were calculated.'%(expct_n_zeros, n_zeros))
-    else:
-        LOGGER.warning('More than %d zero eigenvalues were detected.'%expct_n_zeros)
-
-    if not zeros:
-        if n_zeros > expct_n_zeros:
-            if n_zeros == n_modes + expct_n_zeros and n_modes < dof:
-                LOGGER.debug('Determing the number of zero eigenvalues...')
-                # find the actual number of zero modes
-                n_zeros = _calc_n_zero_modes(M)
-                LOGGER.debug('%d zero eigenvalues detected.'%n_zeros)
-            LOGGER.debug('Solving for additional eigenvalues...')
-
-            if n_modes < dof:
-                start = min(n_modes+expct_n_zeros, dof-1); end = min(n_modes+n_zeros-1, dof-1)
-                values_, vectors_ = _eigh(M, eigvals=(start, end))
-                values = np.concatenate((values, values_))
-                vectors = np.hstack((vectors, vectors_))
-
-        # final_n_modes may exceed len(eigvals) - no need to fix for the sake of the simplicity of the code
-        final_n_modes = n_zeros + n_modes
-        eigvals = values[n_zeros:final_n_modes]
-        eigvecs = vectors[:, n_zeros:final_n_modes]
-        vars = 1 / eigvals
-    else:
-        eigvals = values[:n_modes]
-        eigvecs = vectors[:, :n_modes]
-        vars = div0(1, values)
-        vars[:n_zeros] = 0.
-        vars = vars[:n_modes]
-
-    return eigvals, eigvecs, vars
-
+__all__ = ['GNM', 'MaskedGNM', 'calcGNM']
 
 class GNMBase(NMA):
 

prody/dynamics/pca.py

@@ -10,7 +10,7 @@
 from prody.atomic import Atomic
 from prody.ensemble import Ensemble, PDBEnsemble
 from prody.trajectory import TrajBase
-from prody.utilities import importLA
+from prody.utilities import importLA, solveEig, ZERO
 
 from .nma import NMA
 
@@ -36,7 +36,13 @@ def setCovariance(self, covariance, is3d=True):
             raise TypeError('covariance must be an ndarray')
         elif not (covariance.ndim == 2 and
                   covariance.shape[0] == covariance.shape[1]):
-            raise TypeError('covariance must be square matrix')
+            raise ValueError('covariance must be square matrix')
+        elif covariance.dtype != float:
+            try:
+                covariance = covariance.astype(float)
+            except:
+                raise ValueError('covariance.dtype must be float')
+
         self._reset()
 
         self._is3d = is3d
@@ -190,39 +196,19 @@ def calcModes(self, n_modes=20, turbo=True):
             calculate modes, default is **True**
         :type turbo: bool"""
 
-        linalg = importLA()
         if self._cov is None:
             raise ValueError('covariance matrix is not built or set')
         start = time.time()
-        dof = self._dof
         self._clear()
         if str(n_modes).lower() == 'all':
             n_modes = None
-        if linalg.__package__.startswith('scipy'):
-            if n_modes is None:
-                eigvals = None
-                n_modes = dof
-            else:
-                n_modes = int(n_modes)
-                if n_modes >= self._dof:
-                    eigvals = None
-                    n_modes = dof
-                else:
-                    eigvals = (dof - n_modes, dof - 1)
-            values, vectors = linalg.eigh(self._cov, turbo=turbo,
-                                          eigvals=eigvals)
-        else:
-            if n_modes is not None:
-                LOGGER.info('Scipy is not found, all modes are calculated.')
-            values, vectors = linalg.eigh(self._cov)
-        # Order by descending SV
-        revert = list(range(len(values)-1, -1, -1))
-        values = values[revert]
-        vectors = vectors[:, revert]
-        which = values > 1e-8
+
+        values, vectors, _ = solveEig(self._cov, n_modes=n_modes, zeros=True, 
+                                      turbo=turbo, warn_zeros=False, reverse=True)
+        which = values > ZERO
         self._eigvals = values[which]
         self._array = vectors[:, which]
-        self._vars = self._eigvals
+        self._vars = values
         self._n_modes = len(self._eigvals)
         LOGGER.debug('{0} modes were calculated in {1:.2f}s.'
                      .format(self._n_modes, time.time()-start))

prody/dynamics/plotting.py

@@ -823,8 +823,12 @@ def showContactMap(enm, **kwargs):
 def showOverlap(mode, modes, *args, **kwargs):
     """Show overlap :func:`~matplotlib.pyplot.bar`.
 
-    :arg mode: a single mode/vector
-    :type mode: :class:`.Mode`, :class:`.Vector`
+    :arg mode: a single mode/vector or multiple modes.
+        If multiple modes are provided, then the overlaps are calculated 
+        by going through them one by one, i.e. mode i from this set is 
+        compared with mode i from the other set.
+    :type mode: :class:`.Mode`, :class:`.Vector`, :class:`.ModeSet`, 
+        :class:`.ANM`, :class:`.GNM`, :class:`.PCA`
 
     :arg modes: multiple modes
     :type modes: :class:`.ModeSet`, :class:`.ANM`, :class:`.GNM`, :class:`.PCA`
@@ -836,13 +840,15 @@ def showOverlap(mode, modes, *args, **kwargs):
     if SETTINGS['auto_show']:
         plt.figure()
 
-    if not isinstance(mode, (Mode, Vector)):
-        raise TypeError('mode must be Mode or Vector, not {0}'
+    if not isinstance(mode, (Mode, Vector, NMA, ModeSet)):
+        raise TypeError('mode must be Mode, Vector, NMA or ModeSet, not {0}'
                         .format(type(mode)))
+
     if not isinstance(modes, (NMA, ModeSet)):
         raise TypeError('modes must be NMA or ModeSet, not {0}'
                         .format(type(modes)))
-    overlap = abs(calcOverlap(mode, modes))
+
+    overlap = abs(calcOverlap(mode, modes, diag=True))
     if isinstance(modes, NMA):
         arange = np.arange(len(modes)) + 1
     else:

prody/dynamics/rtb.py

@@ -26,17 +26,11 @@ def __call__(self, i=1):
 class RTB(ANMBase):
 
     """Class for Rotations and Translations of Blocks (RTB) method ([FT00]_).
-    Optional arguments permit imposing constraints along Z-direction as in
-    *imANM* method described in [TL12]_.
 
     .. [FT00] Tama F, Gadea FJ, Marques O, Sanejouand YH. Building-block
        approach for determining low-frequency normal modes of macromolecules.
        *Proteins* **2000** 41:1-7.
 
-    .. [TL12] Lezon TR, Bahar I, Constraints Imposed by the Membrane
-       Selectively Guide the Alternating Access Dynamics of the Glutamate
-       Transporter GltPh
-
     """
 
     def __init__(self, name='Unknown'):