Merge pull request #781 from nkeim/doc-gapped-msd

DOC: Gapped MSD

Author: nkeim

doc/releases/v0.7.txt

@@ -0,0 +1,12 @@
+v0.7
+------
+
+trackpy v0.7 includes an enhancement/bug fix to how MSD calculations handle gaps
+in trajectories. The msd() and imsd() functions now output NaN values (instead
+of zeros) for portions of their output that could not be measured. The emsd()
+function now properly de-weights missing data, and so its output may be
+significantly different (and more accurate) when gaps are numerous and large.
+
+Bug fixes
+~~~~~~~~~
+- MSD properly handles gaps in trajectories (@vivarose, #773)

trackpy/motion.py

@@ -26,11 +26,15 @@ def msd(traj, mpp, fps, max_lagtime=100, detail=False, pos_columns=None):
 
         If detail is True, the DataFrame also contains a column N,
         the estimated number of statistically independent measurements
-        that comprise the result at each lagtime.
+        that comprise the result at each lagtime. Where gaps make the
+        MSD undefined, N is set to 0. All N values are then rescaled to 
+        make the total unchanged, so that the particle will be weighted
+        correctly in an ensemble calculation.
 
     Notes
     -----
     Input units are pixels and frames. Output units are microns and seconds.
+    Outputs NaN for lag times that could not be measured, due to gaps.
 
     See also
     --------
@@ -81,7 +85,13 @@ def _msd_iter(pos, lagtimes):
 
 def _msd_gaps(traj, mpp, fps, max_lagtime=100, detail=False, pos_columns=None):
     """Compute the mean displacement and mean squared displacement of one
-    trajectory over a range of time intervals."""
+    trajectory over a range of time intervals.
+       
+    If 'detail' is True, N is set to 0 wherever gaps make the MSD undefined. All
+    N values are then rescaled to make the total unchanged, so that the particle
+    will be weighted correctly in an ensemble calculation. For details, see
+    the discussion in PR #773: https://github.com/soft-matter/trackpy/pull/773
+    """
     if pos_columns is None:
         pos_columns = ['x', 'y']
     result_columns = ['<{}>'.format(p) for p in pos_columns] + \
@@ -196,6 +206,7 @@ def imsd(traj, mpp, fps, max_lagtime=100, statistic='msd', pos_columns=None):
     Notes
     -----
     Input units are pixels and frames. Output units are microns and seconds.
+    Outputs NaN for lag times that could not be measured, due to gaps.
     """
     ids = []
     msds = []

trackpy/tests/test_motion.py

@@ -161,6 +161,13 @@ def setUp(self):
                        'frame': np.arange(1, N), 'particle': np.ones(N - 1)})
         self.steppers = conformity(pandas_concat([a, b]))
 
+        self.badly_gapped_walks = conformity(DataFrame({
+            'frame': [1, 2, 3, 4, 5, 6, 1, 2, 6], # create an example trajectory where particle 2 disappears for a few frames and returns.
+            'x': [4, 5, 4, 5, 3, 2, -3, -1, -2],
+            'y': [6, 7, 8, 6, 7, 6, 10, 11, 10],
+            'particle': [1, 1, 1, 1, 1, 1, 2, 2, 2]
+        }))
+
     def test_zero_emsd(self):
         N = 10
         actual = tp.emsd(self.dead_still, 1, 1)
@@ -207,6 +214,33 @@ def test_linear_emsd_gaps(self):
         actual.index = expected.index
         assert_series_equal(np.round(actual), expected)
 
+    def test_major_gap_imsd(self):
+        """Large gap that should strongly affect emsd.
+        
+        Test data and fix (PR #773) by @vivarose
+        """
+        imsd = tp.imsd(self.badly_gapped_walks, mpp=1, fps=1)
+        nans = imsd.isna()
+
+        assert nans[2][2.0]
+        assert nans[2][3.0]
+        assert sum(nans.values.flatten()) == 2
+
+    def test_major_gap_emsd(self):
+        """Large gap that should strongly affect emsd.
+        
+        Test data and fix (PR #773) by @vivarose
+        """
+        actual = tp.emsd(self.badly_gapped_walks, mpp=1, fps=1)
+        expected = Series({
+            1.0: 4.012847555129437,
+            2.0: 4.000000000000007,
+            3.0: 4.333333333333333,
+            4.0: 4.193672099712366,
+            5.0: 2.645254074784259
+            })
+        assert_almost_equal(expected.values, actual.values)
+
     def test_direction_corr(self):
         # just a smoke test
         f1, f2 = 2, 6