Merge 31ff563 into fcc6a23

openpathsampling/__init__.py

@@ -81,6 +81,9 @@
     trace_replicas_for_ensemble, condense_repeats,
     ReplicaNetworkGraph
 )
+from analysis.single_trajectory_analysis import (
+    SingleTrajectoryAnalysis 
+)
 
 from live_visualization import LiveVisualization
 

openpathsampling/analysis/single_trajectory_analysis.py

@@ -0,0 +1,66 @@
+import openpathsampling as paths
+import numpy as np
+
+
+class SingleTrajectoryAnalysis(object):
+    def __init__(self, transition, dt=None):
+        self.transition = transition
+        self.dt = dt
+        self.stateA = transition.stateA
+        self.stateB = transition.stateB
+        self.continuous_time = {self.stateA : [], self.stateB : []}
+        self.lifetime = {self.stateA : [], self.stateB : []}
+        self.flux = {self.stateA : {}, self.stateB : {}}
+
+
+    def analyze_continuous_time(self, trajectory, state):
+        ensemble = paths.AllInXEnsemble(state)
+        segments = ensemble.split(trajectory)
+        lengths = [len(seg) for seg in segments]
+        self.continuous_time[state] += lengths
+
+    def analyze_lifetime(self, trajectory, state):
+        other_state = list(set([self.stateA, self.stateB]) - set([state]))[0]
+	ensemble_BAB = paths.SequentialEnsemble([
+	    paths.AllInXEnsemble(other_state) & paths.LengthEnsemble(1),
+	    paths.PartInXEnsemble(state) & paths.AllOutXEnsemble(other_state),
+	    paths.AllInXEnsemble(other_state) & paths.LengthEnsemble(1)
+	])
+	ensemble_AB = paths.SequentialEnsemble([
+	    paths.AllInXEnsemble(state) & paths.LengthEnsemble(1),
+	    paths.OptionalEnsemble(
+		paths.PartInXEnsemble(state) &
+		paths.AllOutXEnsemble(other_state)
+	    ),
+	    paths.AllInXEnsemble(other_state) & paths.LengthEnsemble(1)
+	])
+
+
+    def analyze_flux(self, trajectory, state):
+        pass
+
+    def add_frames(self, trajectory):
+        for state in [self.stateA, self.stateB]:
+            self.analyze_continuous_time(trajectory, state)
+            self.analyze_lifetime(trajectory, state)
+            self.analyze_flux(trajectory, state)
+
+    def analyze(self, trajectory):
+        self.add_frames(trajectory)
+        return self
+
+    def continuous_time_distribution(self, state):
+        dist = np.array(self.continuous_time[state])
+        if self.dt is not None:
+            dist *= self.dt
+        return dist
+
+    def lifetime_distribution(self, state):
+        dist = np.array(self.lifetime[state])
+        if self.dt is not None:
+            dist *= self.dt
+        return dist
+
+    def summary(self):
+        # return a nice text/pandas output
+        pass

openpathsampling/tests/testtrajectory.py

@@ -77,4 +77,58 @@ def test_summarize_trajectory_volumes_with_nonevol(self):
             "A-I-B-None-B-None"
         )
 
+class testSubtrajectoryIndices(object):
+    def setup(self):
+        op = paths.CV_Function("Id", lambda snap : snap.coordinates[0][0])
+        vol1 = paths.CVRangeVolume(op, 0.1, 0.5)
+        vol2 = paths.CVRangeVolume(op, -0.1, 0.7)
+        vol3 = paths.CVRangeVolume(op, 2.0, 2.5)
+
+        self.stateA = vol1
+        self.interstitial = vol2 & ~vol1
+        self.outInterface = ~vol2 & ~vol3
+        self.stateB = vol3
+
+    def test_subtrajectory_indices(self):
+        # simplify more complicated expressions
+        stateA = self.stateA
+        stateB = self.stateB
+        pretraj = [0.20, 0.30, 0.60, 0.40, 0.65, 2.10, 2.20, 2.60, 2.10,
+                   0.80, 0.55, 0.40, 0.20]
+        # 00, 01, 02, 03, 04, 05, 06, 07, 08, 09, 10, 11, 12
+        #  A,  A,  I,  A,  I,  B,  B,  X,  B,  X,  I,  A,  A
+        trajectory = make_1d_traj(
+            coordinates=pretraj,
+            velocities=[1.0]*len(pretraj)
+        )
+        ensemble_A = paths.AllInXEnsemble(stateA)
+        ensemble_B = paths.AllInXEnsemble(stateB)
+        ensemble_ABA = paths.SequentialEnsemble([
+            paths.AllInXEnsemble(stateA) & paths.LengthEnsemble(1),
+            paths.PartInXEnsemble(stateB) & paths.AllOutXEnsemble(stateA),
+            paths.AllInXEnsemble(stateA) & paths.LengthEnsemble(1)
+        ])
+        subtrajectoriesA = ensemble_A.split(trajectory)
+        subtrajectoriesB = ensemble_B.split(trajectory)
+        subtrajectoriesABA = ensemble_ABA.split(trajectory)
 
+        # make sure we have the trajectories we expect
+        assert_equal(len(subtrajectoriesA), 3)
+        assert_equal(len(subtrajectoriesB), 2)
+        assert_equal(len(subtrajectoriesABA), 1)
+        # the following assertions check that the subtrajectories are the
+        # ones that we expect; the numbers here are linked to the indices
+        # we'll test next
+        assert_equal(subtrajectoriesA[0], trajectory[0:2])
+        assert_equal(subtrajectoriesA[1], trajectory[3:4])
+        assert_equal(subtrajectoriesA[2], trajectory[11:13])
+        assert_equal(subtrajectoriesB[0], trajectory[5:7])
+        assert_equal(subtrajectoriesB[1], trajectory[8:9])
+        assert_equal(subtrajectoriesABA[0], trajectory[3:12])
+        # now we run the subtrajectory_indices function and test it
+        indicesA = trajectory.subtrajectory_indices(subtrajectoriesA)
+        indicesB = trajectory.subtrajectory_indices(subtrajectoriesB)
+        indicesABA = trajectory.subtrajectory_indices(subtrajectoriesABA)
+        assert_equal(indicesA, [[0, 1], [3], [11, 12]])
+        assert_equal(indicesB, [[5, 6], [8]])
+        assert_equal(indicesABA, [[3, 4, 5, 6, 7, 8, 9, 10, 11]])

openpathsampling/trajectory.py

@@ -572,19 +572,9 @@ def is_correlated(self, other):
         bool
             returns True if at least one snapshot appears in both trajectories
         """
-
-        # if hasattr(self, 'idx') and hasattr(other, 'idx'):
-        #     shared_store = set(self.idx.keys()) & set(other.idx.keys())
-        #     # both are saved so use the snapshot idx as identifiers
-        #     if len(shared_store) > 0:
-        #         storage = list(shared_store)[0]
-        #         t1id = storage.trajectories.snapshot_indices(self.idx[storage])
-        #         t2id = storage.trajectories.snapshot_indices(other.idx[storage])
-        #         return bool(set(t1id) & set(t2id))
-
-        # Use some fallback
         return bool(self.shared_configurations(other))
 
+
     def shared_configurations(self, other):
         """
         Returns a set of shared snapshots
@@ -601,6 +591,7 @@ def shared_configurations(self, other):
         """
         return set([snap for snap in self]) & set([snap for snap in other])
 
+
     def shared_subtrajectory(self, other):
         """
         Returns a subtrajectory which only contains frames present in other
@@ -618,6 +609,25 @@ def shared_subtrajectory(self, other):
         shared = self.shared_configurations(other)
         return Trajectory([snap for snap in self if snap in shared])
 
+
+    def subtrajectory_indices(self, subtrajectories):
+        """
+        Returns a list of lists of indices for frames from subtrajectories.
+
+        Parameters
+        ----------
+        subtrajectories : list of :class:`.Trajectory`
+            input list of subtrajectories
+
+        Returns
+        -------
+        list of list of int
+            the indices within this trajectory of the frames in each
+            subtrajectory
+        """
+        return [[self.index(s) for s in subtrj] for subtrj in subtrajectories]
+
+
     # =============================================================================================
     # UTILITY FUNCTIONS
     # =============================================================================================
@@ -715,4 +725,4 @@ def topology(self):
             if self.atom_indices is not None:
                 topology = topology.subset(self.atom_indices)
 
-        return topology
+        return topology