r28223 - in /trunk: lib/structure/statistics.py test_suite/system_tests/structure.py -- June 02, 2016 - 19:09

Author: bugman
Date: Thu Jun  2 19:09:40 2016
New Revision: 28223

URL: http://svn.gna.org/viewcvs/relax?rev=28223&view=rev
Log:
Fix for bug #24723 (https://gna.org/bugs/?24723).

This is the bug that the mean RMSD from the structure.rmsd user function is 
incorrectly
calculated - it should be a quadratic mean.

The quadratic mean and quadratic standard deviation are now correctly 
calculated, and the
structure.test_rmsd, structure.test_rmsd_molecules, and 
structure.test_rmsd_ubi system tests have
been updated for the fix.

Modified:
    trunk/lib/structure/statistics.py
    trunk/test_suite/system_tests/structure.py

Modified: trunk/lib/structure/statistics.py
URL: 
http://svn.gna.org/viewcvs/relax/trunk/lib/structure/statistics.py?rev=28223&r1=28222&r2=28223&view=diff
==============================================================================
--- trunk/lib/structure/statistics.py   (original)
+++ trunk/lib/structure/statistics.py   Thu Jun  2 19:09:40 2016
@@ -65,16 +65,16 @@
         if verbosity:
             print("Model %2s RMSD:  %s" % (i, model_rmsd[i]))
 
-    # Calculate the mean.
-    rmsd_mean = mean(model_rmsd)
+    # Calculate the quadratic mean.
+    rmsd_mean = sqrt(mean(model_rmsd**2))
 
-    # Calculate the normal non-biased standard deviation.
+    # Calculate the normal non-biased quadratic standard deviation.
     try:
-        rmsd_sd = std(model_rmsd, ddof=1)
+        rmsd_sd = sqrt(std(model_rmsd**2, ddof=1))
 
     # Handle old numpy versions not having the ddof argument.
     except TypeError:
-        rmsd_sd = std(model_rmsd) * sqrt(M / (M - 1.0))
+        rmsd_sd = sqrt(std(model_rmsd**2) * sqrt(M / (M - 1.0)))
 
     # Printout.
     if verbosity:

Modified: trunk/test_suite/system_tests/structure.py
URL: 
http://svn.gna.org/viewcvs/relax/trunk/test_suite/system_tests/structure.py?rev=28223&r1=28222&r2=28223&view=diff
==============================================================================
--- trunk/test_suite/system_tests/structure.py  (original)
+++ trunk/test_suite/system_tests/structure.py  Thu Jun  2 19:09:40 2016
@@ -21,7 +21,7 @@
 
 # Python module imports.
 from math import sqrt
-from numpy import array, average, dot, float64, sign, std, zeros
+from numpy import array, average, dot, float64, mean, sign, std, zeros
 from numpy.linalg import norm
 from os import sep
 from re import search
@@ -5070,7 +5070,7 @@
 
         # Checks.
         self.assert_(hasattr(cdp.structure, 'rmsd'))
-        self.assertAlmostEqual(cdp.structure.rmsd, 2./3*sqrt(2))
+        self.assertAlmostEqual(cdp.structure.rmsd, sqrt(4.0/3.0))
 
 
     def test_rmsd_molecules(self):
@@ -5095,7 +5095,7 @@
 
         # Checks.
         self.assert_(hasattr(cdp.structure, 'rmsd'))
-        self.assertAlmostEqual(cdp.structure.rmsd, 2./3*sqrt(2))
+        self.assertAlmostEqual(cdp.structure.rmsd, sqrt(4.0/3.0))
 
 
     def test_rmsd_spins(self):
@@ -5149,9 +5149,13 @@
         # Calculate the RMSD.
         self.interpreter.structure.rmsd()
 
+        # The per-structure RMSDs, and their quadratic average.
+        all_rmsd = array([1.06994466835, 0.411387603119, 0.647768214719, 
0.522216519591, 0.274450954939, 0.979817275482, 0.752817901842, 
1.28956426491, 1.12979370978, 0.650514765397], float64)
+        rmsd = sqrt(mean(all_rmsd**2))
+
         # Checks (the values match the VMD 1.9.1 RMSD numbers).
         self.assert_(hasattr(cdp.structure, 'rmsd'))
-        self.assertAlmostEqual(cdp.structure.rmsd, 0.77282758781333061)
+        self.assertAlmostEqual(cdp.structure.rmsd, rmsd)
 
 
     def test_sequence_alignment_central_star_nw70_blosum62(self):