Switch PersistenceLength to Results class (#3289)

Fixes #3291

## Work done in this PR
* Switch PersistenceLength to Results class

package/CHANGELOG

@@ -175,6 +175,9 @@ Enhancements
     checking if it can be used in parallel analysis. (Issue #2996, PR #2950)
 
 Changes
+  * `analysis.polymer.PersistenceLength` class now stores `lb`,
+    `lp` and `fit` using the `analysis.base.Results` class 
+    (Issues #3289, #3291)
   * `analysis.hole2.HoleAnalysis` now stores ``sphpdbs``, ``outfiles``
     and ``profiles`` in the `analysis.base.Results` class (Issues #3261, #3269)
   * `helix_analysis.HELANAL` now uses the `analysis.base.Results` class to
@@ -247,6 +250,13 @@ Changes
   * Added OpenMM coordinate and topology converters (Issue #2863, PR #2917)
 
 Deprecations
+  * The `analysis.polymer.PersistenceLength.lb`,
+    `analysis.polymer.PersistenceLength.lp` and 
+    `analysis.polymer.PersistenceLength.fit` attributes are now deprecated in
+    favour of `analysis.polymer.PersistenceLength.results.lb`,
+    `analysis.polymer.PersistenceLength.results.lp` and
+    `analysis.polymer.PersistenceLength.results.fit` respectively. They will
+    be removed in 3.0.0 (Issues #3289, #3291)
   * The ``sphpdbs``, ``outfiles`` and ``profiles`` attributes of
     `analysis.hole2.HoleAnalysis` are now deprecated in favour of
     ``results.sphpdbs``, ``results.outfiles`` and

package/MDAnalysis/analysis/polymer.py

@@ -31,34 +31,6 @@
 :Copyright: GNU Public License v3
 
 This module contains various commonly used tools in analysing polymers.
-
-
-Persistence length worked example
----------------------------------
-
-This example shows how to find the persistence length of a polymer
-using MDAnalysis.
-::
-
-  from MDAnalysis.tests.datafiles import TRZ_psf, TRZ
-  import MDAnalysis as mda
-  from MDAnalysis.analysis import polymer
-  u = mda.Universe(TRZ_psf, TRZ)
-  # this system is a pure polymer melt of polyamide,
-  # so we can select the chains by using the .fragments attribute
-  chains = u.atoms.fragments
-  # select only the backbone atoms for each chain
-  backbones = [chain.select_atoms('not name O* H*') for chain in chains]
-  # sort the chains, removing any non-backbone atoms
-  sorted_backbones = [polymer.sort_backbone(bb) for bb in backbones]
-  persistence_length = polymer.PersistenceLength(sorted_backbones)
-  # Run the analysis, this will average over all polymer chains
-  # and all timesteps in trajectory
-  persistence_length = persistence_length.run()
-  print('The persistence length is: {}'.format(persistence_length.lp))
-  # always check the visualisation of this:
-  persistence_length.plot()
-
 """
 import numpy as np
 import scipy.optimize
@@ -160,34 +132,102 @@ class PersistenceLength(AnalysisBase):
     Parameters
     ----------
     atomgroups : iterable
-       List of AtomGroups.  Each should represent a single
+       List of AtomGroups. Each should represent a single
        polymer chain, ordered in the correct order.
-    verbose : bool (optional)
+    verbose : bool, optional
        Show detailed progress of the calculation if set to ``True``.
 
     Attributes
     ----------
     results.bond_autocorrelation : numpy.ndarray
        the measured bond autocorrelation
-    lb : float
+    results.lb : float
        the average bond length
-    lp : float
+
+       .. versionadded:: 2.0.0
+
+    lb : float
+
+       Alias to the :attr:`results.lb`.
+
+       .. deprecated:: 2.0.0
+            Will be removed in MDAnalysis 3.0.0. Please use
+            :attr:`results.lb` instead.
+
+    results.x : numpy.ndarray
+        length of the decorrelation predicted by *lp*
+
+        .. versionadded:: 2.0.0
+
+    results.lp : float
        calculated persistence length
-    fit : numpy.ndarray
+
+       .. versionadded:: 2.0.0
+
+    lp : float
+
+       Alias to the :attr:`results.lp`.
+
+       .. deprecated:: 2.0.0
+            Will be removed in MDAnalysis 3.0.0. Please use
+            :attr:`results.lp` instead.
+
+    results.fit : numpy.ndarray
        the modelled backbone decorrelation predicted by *lp*
 
+       .. versionadded:: 2.0.0
+
+    fit : float
+
+       Alias to the :attr:`results.fit`.
+
+       .. deprecated:: 2.0.0
+            Will be removed in MDAnalysis 3.0.0. Please use
+            :attr:`results.fit` instead.
+
     See Also
     --------
     :func:`sort_backbone`
        for producing the sorted AtomGroup required for input.
 
+    Example
+    -------
+    .. code-block:: python
+
+        from MDAnalysis.tests.datafiles import TRZ_psf, TRZ
+        import MDAnalysis as mda
+        from MDAnalysis.analysis import polymer
+        u = mda.Universe(TRZ_psf, TRZ)
+
+        # this system is a pure polymer melt of polyamide,
+        # so we can select the chains by using the .fragments attribute
+        chains = u.atoms.fragments
+
+        # select only the backbone atoms for each chain
+        backbones = [chain.select_atoms('not name O* H*') for chain in chains]
+
+        # sort the chains, removing any non-backbone atoms
+        sorted_backbones = [polymer.sort_backbone(bb) for bb in backbones]
+        persistence_length = polymer.PersistenceLength(sorted_backbones)
+
+        # Run the analysis, this will average over all polymer chains
+        # and all timesteps in trajectory
+        persistence_length = persistence_length.run()
+        print(f'The persistence length is: {persistence_length.results.lp}')
+
+        # always check the visualisation of this:
+        persistence_length.plot()
+
+
     .. versionadded:: 0.13.0
     .. versionchanged:: 0.20.0
        The run method now automatically performs the exponential fit
     .. versionchanged:: 1.0.0
        Deprecated :meth:`PersistenceLength.perform_fit` has now been removed.
     .. versionchanged:: 2.0.0
-       Former ``results`` are now stored as ``results.bond_autocorrelation``
+       Former ``results`` are now stored as ``results.bond_autocorrelation``.
+       :attr:`lb`, :attr:`lp`, :attr:`fit` are now stored in a
+       :class:`MDAnalysis.analysis.base.Results` instance.
     """
     def __init__(self, atomgroups, **kwargs):
         super(PersistenceLength, self).__init__(
@@ -219,6 +259,30 @@ def _single_frame(self):
             for i in range(n-1):
                 self._results[:(n-1)-i] += inner_pr[i, i:]
 
+    @property
+    def lb(self):
+        wmsg = ("The `lb` attribute was deprecated in "
+                "MDAnalysis 2.0.0 and will be removed in MDAnalysis 3.0.0. "
+                "Please use `results.variance` instead.")
+        warnings.warn(wmsg, DeprecationWarning)
+        return self.results.lb
+
+    @property
+    def lp(self):
+        wmsg = ("The `lp` attribute was deprecated in "
+                "MDAnalysis 2.0.0 and will be removed in MDAnalysis 3.0.0. "
+                "Please use `results.variance` instead.")
+        warnings.warn(wmsg, DeprecationWarning)
+        return self.results.lp
+
+    @property
+    def fit(self):
+        wmsg = ("The `fit` attribute was deprecated in "
+                "MDAnalysis 2.0.0 and will be removed in MDAnalysis 3.0.0. "
+                "Please use `results.variance` instead.")
+        warnings.warn(wmsg, DeprecationWarning)
+        return self.results.fit
+
     def _conclude(self):
         n = len(self._atomgroups[0])
 
@@ -237,20 +301,21 @@ def _calc_bond_length(self):
             pos = ag.positions
             b = calc_bonds(pos[:-1], pos[1:]).mean()
             bs.append(b)
-        self.lb = np.mean(bs)
+        self.results.lb = np.mean(bs)
 
     def _perform_fit(self):
         """Fit the results to an exponential decay"""
         try:
             self.results.bond_autocorrelation
         except AttributeError:
             raise NoDataError("Use the run method first") from None
-        self.x = np.arange(len(self.results.bond_autocorrelation)) * self.lb
+        self.results.x = self.results.lb *\
+                            np.arange(len(self.results.bond_autocorrelation))
 
-        self.lp = fit_exponential_decay(self.x,
-                                        self.results.bond_autocorrelation)
+        self.results.lp = fit_exponential_decay(self.results.x,
+                                                self.results.bond_autocorrelation)
 
-        self.fit = np.exp(-self.x/self.lp)
+        self.results.fit = np.exp(-self.results.x/self.results.lp)
 
     def plot(self, ax=None):
         """Visualise the results and fit
@@ -267,16 +332,16 @@ def plot(self, ax=None):
         import matplotlib.pyplot as plt
         if ax is None:
             fig, ax = plt.subplots()
-        ax.plot(self.x,
+        ax.plot(self.results.x,
                 self.results.bond_autocorrelation,
                 'ro',
                 label='Result')
-        ax.plot(self.x,
-                self.fit,
+        ax.plot(self.results.x,
+                self.results.fit,
                 label='Fit')
         ax.set_xlabel(r'x')
         ax.set_ylabel(r'$C(x)$')
-        ax.set_xlim(0.0, 40 * self.lb)
+        ax.set_xlim(0.0, 40 * self.results.lb)
 
         ax.legend(loc='best')
 

testsuite/MDAnalysisTests/analysis/test_persistencelength.py

@@ -65,11 +65,11 @@ def test_run(self, p_run):
         assert len(p_run.results.bond_autocorrelation) == 280
 
     def test_lb(self, p_run):
-        assert_almost_equal(p_run.lb, 1.485, 3)
+        assert_almost_equal(p_run.results.lb, 1.485, 3)
 
     def test_fit(self, p_run):
-        assert_almost_equal(p_run.lp, 6.504, 3)
-        assert len(p_run.fit) == len(p_run.results.bond_autocorrelation)
+        assert_almost_equal(p_run.results.lp, 6.504, 3)
+        assert len(p_run.results.fit) == len(p_run.results.bond_autocorrelation)
 
     def test_raise_NoDataError(self, p):
         #Ensure that a NoDataError is raised if perform_fit()
@@ -96,6 +96,13 @@ def test_current_axes(self, p_run):
         ax2 = p_run.plot(ax=None)
         assert ax2 is not ax
 
+    @pytest.mark.parametrize("attr", ("lb", "lp", "fit"))
+    def test(self, p, attr):
+        p_run = p.run(step=3)
+        wmsg = f"The `{attr}` attribute was deprecated in MDAnalysis 2.0.0"
+        with pytest.warns(DeprecationWarning, match=wmsg):
+            getattr(p_run, attr) is p_run.results[attr]
+
 
 class TestFitExponential(object):
     x = np.linspace(0, 250, 251)