improve performance of InterRDF_s (#5073)

* Fixes #5067
* Refactored bin index calculation in analysis.rdf.InterRDF_s._single_frame to avoid multiple redundant 
   np.histogram calls, improving performance; at least 20x for small systems, potentially >100x for realistic systems.
* fixed change entries for PR #4884 on RDF parallelization
  - corrected version changed dates for PR #4884 in analysis.rdf to 2.10.0
  - moved entry in CHANGELOG to 2.10.0
* add benchmark for InterRDF_s
* Update CHANGELOG
* update AUTHORS

---------

Co-authored-by: Oliver Beckstein <[email protected]>
Co-authored-by: Irfan Alibay <[email protected]>

Author: 3 people

benchmarks/benchmarks/analysis/rdf.py

@@ -6,13 +6,13 @@
     pass
 
 try:
-    from MDAnalysis.analysis.rdf import InterRDF
+    from MDAnalysis.analysis.rdf import InterRDF, InterRDF_s
 except:
     pass
 
 
 class SimpleRdfBench(object):
-    """Benchmarks for MDAnalysis.analysis.rdf"""
+    """Benchmarks for MDAnalysis.analysis.rdf.InterRDF"""
 
     params = ([20, 75, 200], [[0, 5], [0, 15], [0, 20]], [1, 100, 1000, 10000])
 
@@ -39,3 +39,31 @@ def time_interrdf(self, nbins, range_val, natoms):
         by MDAnalysis.analysis.rdf.InterRDF
         """
         self.rdf.run()
+
+class SimpleRdfsBench(object):
+    """Benchmarks for MDAnalysis.analysis.rdf.InterRDF_s"""
+
+    params = ([20, 75, 200], [[0, 5], [0, 15], [0, 20]], [10, 100], [1, 3, 9])
+
+    param_names = ["nbins", "range_val", "natoms", "npairs"]
+
+    def setup(self, nbins, range_val, natoms, npairs):
+
+        self.sel_str = "name OW"
+
+        self.u = MDAnalysis.Universe(TPR, XTC)
+
+        try:
+            self.sel = self.u.select_atoms(self.sel_str)[:natoms]
+        except AttributeError:
+            self.sel = self.u.selectAtoms(self.sel_str)[:natoms]
+
+        ags = [[self.sel, self.sel]] * npairs
+        self.rdf_s = InterRDF_s(self.u, ags, nbins=nbins, range=range_val)
+
+
+    def time_interrdfs(self, nbins, range_val, natoms, npairs):
+        """Benchmark a full trajectory parse
+        by MDAnalysis.analysis.rdf.InterRDF_s
+        """
+        self.rdf_s.run()

package/AUTHORS

@@ -261,6 +261,7 @@ Chronological list of authors
   - Marc Schuh
   - Sirsha Ganguly
   - Amruthesh Thirumalaiswamy
+  - Ch Zhang
 
 External code
 -------------

package/CHANGELOG

@@ -16,8 +16,8 @@ The rules for this file:
 -------------------------------------------------------------------------------
 ??/??/?? IAlibay, orbeckst, BHM-Bob, TRY-ER, Abdulrahman-PROG, pbuslaev,
          yuxuanzhuang, yuyuan871111, tanishy7777, tulga-rdn, Gareth-elliott,
-         hmacdope, tylerjereddy, cbouy, talagayev, DrDomenicoMarson, amruthesht,
-         p-j-smith
+         hmacdope, tylerjereddy, cbouy, talagayev, DrDomenicoMarson, amruthesht, 
+         p-j-smith, gitzhangch
 
 
  * 2.10.0
@@ -47,6 +47,9 @@ Fixes
     and does not result in an infinite loop (Issue #5112, PR #5113)
 
 Enhancements
+ * Improved performance of `analysis.rdf.InterRDF_s` > 20-fold  
+   by using vectorized NumPy calculations. (Issue #5067, PR #5073)
+ * Parallelize `analysis.rdf.InterRDF` and `analysis.rdf.InterRDF` (Issue #4675)   
  * Support for Python 3.14 was added (PR #5121).
  * Added support for reading and processing streamed data in `coordinates.base`
    with new `StreamFrameIteratorSliced` and `StreamReaderBase` (Issue #4827, PR #4923)
@@ -125,7 +128,6 @@ Enhancements
  * Enable parallelization for analysis.nucleicacids.NucPairDist (Issue #4670)
  * Add check and warning for empty (all zero) coordinates in RDKit converter (PR #4824)
  * Added `precision` for XYZWriter (Issue #4775, PR #4771)
- * Parallelize `analysis.rdf.InterRDF` and `analysis.rdf.InterRDF` (Issue #4675)
  * Added `precision` for XYZWriter (Issue #4775, PR #4771)
  * Add check and warning for empty (all zero) coordinates in RDKit converter
    (PR #4824)

package/MDAnalysis/analysis/rdf.py

@@ -258,7 +258,7 @@ class InterRDF(AnalysisBase):
        of the `results` attribute of
        :class:`~MDAnalysis.analysis.AnalysisBase`.
 
-    .. versionchanged:: 2.9.0
+    .. versionchanged:: 2.10.0
        Enabled **parallel execution** with the ``multiprocessing`` and ``dask``
        backends; use the new method :meth:`get_supported_backends` to see all
        supported backends.
@@ -638,7 +638,7 @@ class InterRDF_s(AnalysisBase):
        Instead of `density=True` use `norm='density'`
     .. deprecated:: 2.3.0
        The `universe` parameter is superflous.
-    .. versionchanged:: 2.9.0
+    .. versionchanged:: 2.10.0
        Enabled **parallel execution** with the ``multiprocessing`` and ``dask``
        backends; use the new method :meth:`get_supported_backends` to see all
        supported backends.
@@ -731,9 +731,16 @@ def _single_frame(self):
                 backend=self.backend,
             )
 
-            for j, (idx1, idx2) in enumerate(pairs):
-                count, _ = np.histogram(dist[j], **self.rdf_settings)
-                self.results.count[i][idx1, idx2, :] += count
+            # Fast manual implementation for distance histogram (equidistant bins)
+            nbins = self.rdf_settings["bins"]
+            rmin, rmax = self.rdf_settings["range"]
+            bin_indices = (dist - rmin) * nbins / (rmax - rmin)
+            bin_indices = bin_indices.astype(np.int64)
+            counts = (bin_indices >= 0) & (bin_indices < nbins)
+            counts = counts.astype(np.int64)
+            idx1s = pairs[:, 0]
+            idx2s = pairs[:, 1]
+            self.results.count[i][idx1s, idx2s, bin_indices] += counts
 
         if self.norm == "rdf":
             self.results.volume_cum += self._ts.volume