Development

avtraj/__init__.py

@@ -3,6 +3,7 @@
 
 """
 from __future__ import annotations
+__version__ = "0.0.10"
 
 import os
 

avtraj/av.py

@@ -243,7 +243,6 @@ def dRDA(
         25.816984253459424
         >>> av_1.dRDA(av_2, distance_sampling_method="sobol_sequence", distance_samples=100)
         """
-        #return av_functions.average_distance(self, av, **kwargs)
         return av_functions.average_distance_labellib(self, av, **kwargs)
 
     def mean_fret_efficiency(
@@ -271,7 +270,6 @@ def mean_fret_efficiency(
         >>> av_1.mean_fret_efficiency(av_2)
         0.9421520823306685
         """
-        #return av_functions.mean_fret_efficiency(self, av, forster_radius, **kwargs)
         return av_functions.mean_fret_efficiency_label_lib(self, av, forster_radius, **kwargs)
 
     def dRDAE(
@@ -626,13 +624,6 @@ def __init__(
         d2 = sd[:, 0] + sd[:, 1] + sd[:, 2]
         xyzr.T[np.where(d2 < asr ** 2)[0]] *= np.array([1.0, 1.0, 1.0, 0.0])
 
-        # calculate an interaction array if no array was provided
-        #if self._xyzrq is None:
-        #    xyzc = np.copy(self.xyzr)
-        #    xyzc[3] += self.contact_volume_thickness
-        #    w = np.ones(xyzc.shape[1], dtype=np.float64)
-        #    self._xyzrq = np.vstack([xyzc, w])
-
         self.attachment_coordinate = attachment_coordinate
         self._min_linker_length = None
         self._grid_density = None

avtraj/av_functions.py

@@ -46,8 +46,6 @@
     'f8', 'f8', 'f8', 'f8'
 ]
 
-sobol_sequence = []#sobol_lib.i4_sobol_generate(6, DISTANCE_SAMPLES)
-
 
 def calculate_av(
         xyzr: np.ndarray,

avtraj/trajectory.py

@@ -3,11 +3,14 @@
 
 import mdtraj
 import numpy as np
+import functools
 
 from .base import PythonBase
 from . import av_functions
 from .av import AccessibleVolume
 
+MAX_CACHE_AV = 100
+
 
 class AVTrajectory(PythonBase):
     """Calculates for a an MD trajectory a corresponding trajectory of
@@ -16,9 +19,7 @@ class AVTrajectory(PythonBase):
 
     Attributes
     ----------
-
     trajectory : mdtraj trajectory
-
     position_name : position name
 
 
@@ -33,7 +34,6 @@ class AVTrajectory(PythonBase):
     >>> import mdtraj as md
     >>> import avtraj as avt
     >>> traj = md.load('./doc/examples/traj.h5')
-
     >>> av_parameters = dict()
     >>> av_parameters['chain_identifier'] = 'A'
     >>> av_parameters['residue_seq_number'] = 344
@@ -71,7 +71,6 @@ def __init__(
         Parameters
         ----------
         name : str
-
         traj : mdtraj Trajectory object or str
             Either a mdtraj Trajectory object or a string containing a the path
             to a trajectory file.
@@ -90,24 +89,16 @@ def __init__(
                     'radius1': 3.5,
                     'simulation_grid_resolution': 0.5,
                 }
-        cache_avs : bool (optional)
-            If cache_avs is True the AVs are stored in an dictionary where the
-            keys of the dictionary correspond to the frame numbers to prevent
-            excessive recalculations. If this parameter is False the AVs for a
-            frame are recalculated and not stored.
         attachment_atom_selection: str (optional)
             This is a MDTraj selection string that is used to define the
             attachment atom. If this selection string is not provided,
             the attachment atom is determined using the parameters provided
             in the av_parameters dictionary.
-
-
         """
         kwargs['name'] = name
         kwargs['verbose'] = kwargs.pop('verbose', False)
         av_parameters['labeling_site_name'] = name
         kwargs['av_parameters'] = av_parameters
-        kwargs['cache_avs'] = kwargs.get('cache_avs', True)
         super().__init__(self, **kwargs)
 
         # Load the trajectory or use the passed mdtraj Trajectory object
@@ -163,6 +154,53 @@ def __init__(
     def __len__(self):
         return len(self.trajectory)
 
+    @functools.lru_cache(maxsize=MAX_CACHE_AV)
+    def get_av(self, frame_i: int):
+        frame = self.trajectory[frame_i]
+        parameters = self.av_parameters
+        parameters['labeling_site_name'] += "_" + str(frame_i)
+        xyz = (frame.xyz[0] * 10.0).astype(np.float64)
+        vdw = self.vdw
+        xyzr = np.vstack([xyz.T, vdw])
+        topology = self.trajectory.topology
+        r1 = parameters.get('radius1', 0.0)
+        r2 = parameters.get('radius2', 0.0)
+        r3 = parameters.get('radius3', 0.0)
+        inter_action_radius = max(r1, r2, r3)
+
+        # compile array xyzrq of acv atoms with radii and weights
+        try:
+            label_interaction_sites = parameters['label_interaction_sites']
+            xyzrq = list()
+            for interaction_site in label_interaction_sites:
+                selection_type, selection = interaction_site["selection"].split(": ")
+                weight = interaction_site["weight"]
+                radius = interaction_site["radius"]
+                if selection_type == "MDTraj":
+                    ai = topology.select(selection)
+                else:
+                    raise AttributeError(
+                        'Only MDTraj selections are allowed as strip mask'
+                    )
+                if len(ai) > 0:
+                    xyzi = xyz[ai]
+                    ri = vdw[ai] + radius + inter_action_radius
+                    wi = np.zeros_like(ri) + float(weight)
+                    xyzrq.append(np.vstack([xyzi.T, ri, wi]).T)
+                    xyzrq_array = np.vstack(xyzrq).T
+                else:
+                    xyzrq_array = []
+        except KeyError:
+            xyzrq_array = []
+        parameters['interaction_sites_xyzrq'] = xyzrq_array
+        attachment_coordinate = xyz[self.attachment_atom_index]
+        av = AccessibleVolume(
+            xyzr,
+            attachment_coordinate,
+            **parameters
+        )
+        return av
+
     def __getitem__(self, key):
         if isinstance(key, int):
             frame_idx = [key]
@@ -173,71 +211,17 @@ def __getitem__(self, key):
             frame_idx = range(start, min(stop, len(self)), step)
         else:
             raise ValueError("")
-
         re = []
         for frame_i in frame_idx:
-            frame = self.trajectory[frame_i]
-
-            # If an AV was already calculated use pre-calculated av
-            if frame_i in self._avs.keys() and self.cache_avs:
-                av = self._avs[frame_i]
-            else:
-                parameters = self.av_parameters
-                parameters['labeling_site_name'] += "_" + str(frame_i)
-                xyz = (frame.xyz[0] * 10.0).astype(np.float64)
-                vdw = self.vdw
-                xyzr = np.vstack([xyz.T, vdw])
-                topology = self.trajectory.topology
-                r1 = parameters.get('radius1', 0.0)
-                r2 = parameters.get('radius2', 0.0)
-                r3 = parameters.get('radius3', 0.0)
-                inter_action_radius = max(r1, r2, r3)
-
-                # compile array xyzrq of acv atoms with radii and weights
-                try:
-                    label_interaction_sites = parameters['label_interaction_sites']
-                    xyzrq = list()
-                    for interaction_site in label_interaction_sites:
-                        selection_type, selection = interaction_site["selection"].split(": ")
-                        weight = interaction_site["weight"]
-                        radius = interaction_site["radius"]
-                        if selection_type == "MDTraj":
-                            ai = topology.select(selection)
-                        else:
-                            raise AttributeError(
-                                'Only MDTraj selections are allowed as strip mask'
-                            )
-                        if len(ai) > 0:
-                            xyzi = xyz[ai]
-                            ri = vdw[ai] + radius + inter_action_radius
-                            wi = np.zeros_like(ri) + float(weight)
-                            xyzrq.append(np.vstack([xyzi.T, ri, wi]).T)
-                            xyzrq_array = np.vstack(xyzrq).T
-                        else:
-                            xyzrq_array = []
-                except KeyError:
-                    xyzrq_array = []
-                parameters['interaction_sites_xyzrq'] = xyzrq_array
-
-                attachment_coordinate = xyz[self.attachment_atom_index]
-                av = AccessibleVolume(
-                    xyzr,
-                    attachment_coordinate,
-                    **parameters
-                )
-                self._avs[frame_i] = av
-
+            av = self.get_av(frame_i)
             re.append(av)
-
         if len(re) == 1:
             return re[0]
         else:
             return re
 
 
-class AvDistanceTrajectory(
-    PythonBase
-):
+class AvDistanceTrajectory(PythonBase):
     """
     The AvPotential class provides the possibility to calculate the reduced
     or unreduced chi2 given a set of labeling positions and experimental
@@ -266,7 +250,6 @@ def __init__(self, trajectory, labeling, **kwargs):
         self.distances = labeling['Distances']
         self.positions = labeling['Positions']
         self.trajectory = trajectory
-
         # Initialize the AV trajectories with the parameters provided by the
         # labeling dictionary
         arguments = [
@@ -290,66 +273,68 @@ def __init__(self, trajectory, labeling, **kwargs):
     def __len__(self):
         return len(self.trajectory)
 
-    def __getitem__(self, key):
+    @functools.lru_cache(maxsize=MAX_CACHE_AV)
+    def get_frame_results(
+            self, frame_idx: int
+    ) -> (float, float, float, float):
+        """ Compute distances and chi2 for a frame
 
+        :param frame_idx: index of the frame
+        :return: mean DA distance rDA, FRET averaged distance rDAE, distance
+        between mean positions rMP, chi2 of the frame
+        """
+        avs = self.avs
+        rDA, rDAE, rMP, chi2 = None, None, None, None
+        for distance_key in self.distances:
+            distance = self.distances[distance_key]
+            av1 = avs[distance['position1_name']][frame_idx]
+            av2 = avs[distance['position2_name']][frame_idx]
+            R0 = distance['Forster_radius']
+
+            rMP = av_functions.dRmp(av1, av2)
+            rDAE = av1.dRDAE(av2, forster_radius=R0)
+            rDA = av1.dRDA(av2)
+
+            if self.verbose:
+                print("RDA: %s" % rDA)
+                print("RDA_E: %s" % rDAE)
+                print("RDA_mp: %s" % rMP)
+            if self.distances[distance_key]['distance_type'] == 'RDAMean':
+                self.distances[distance_key]['model_distance'] = rDA
+            elif self.distances[distance_key]['distance_type'] == 'RDAMeanE':
+                self.distances[distance_key]['model_distance'] = rDAE
+            elif self.distances[distance_key]['distance_type'] == 'Rmp':
+                self.distances[distance_key]['model_distance'] = rMP
+        # compare to experiment: calculate chi2
+        chi2 = 0.0
+        for distance in list(self.distances.values()):
+            dm = distance['model_distance']
+            de = distance['distance']
+            error_neg = distance['error_neg']
+            error_pos = distance['error_pos']
+            d = dm - de
+            chi2 += (d / error_neg) ** 2 if d < 0 else (d / error_pos) ** 2
+        return rDA, rDAE, rMP, chi2
+
+    def __getitem__(self, key):
         if isinstance(key, int):
             frame_idx = [key]
         else:
             start = 0 if key.start is None else key.start
             stop = None if key.stop is None else key.stop
             step = 1 if key.step is None else key.step
             frame_idx = range(start, min(stop, len(self)), step)
-
         re = dict(
             (
                 key,
                 {'rMP': [], 'rDA': [], 'rDAE': [], 'chi2': []}
             ) for key in self.distances.keys()
         )
         for frame_i in frame_idx:
-            # Don't repeat calculations
-            if frame_i in self._d.keys():
-                rDA, rDAE, rMP, chi2 = self._d[frame_i]
-            else:
-                # calculate the AVs of the frame
-                avs = self.avs
-                for distance_key in self.distances:
-                    distance = self.distances[distance_key]
-                    av1 = avs[distance['position1_name']][frame_i]
-                    av2 = avs[distance['position2_name']][frame_i]
-                    R0 = distance['Forster_radius']
-
-                    rMP = av_functions.dRmp(av1, av2)
-                    rDAE = av1.dRDAE(av2, forster_radius=R0)
-                    rDA = av1.dRDA(av2)
-
-                    if self.verbose:
-                        print("RDA: %s" % rDA)
-                        print("RDA_E: %s" % rDAE)
-                        print("RDA_mp: %s" % rMP)
-
-                    if self.distances[distance_key]['distance_type'] == 'RDAMean':
-                        self.distances[distance_key]['model_distance'] = rDA
-                    elif self.distances[distance_key]['distance_type'] == 'RDAMeanE':
-                        self.distances[distance_key]['model_distance'] = rDAE
-                    elif self.distances[distance_key]['distance_type'] == 'Rmp':
-                        self.distances[distance_key]['model_distance'] = rMP
-
-                # compare to experiment: calculate chi2
-                chi2 = 0.0
-                for distance in list(self.distances.values()):
-                    dm = distance['model_distance']
-                    de = distance['distance']
-                    error_neg = distance['error_neg']
-                    error_pos = distance['error_pos']
-                    d = dm - de
-                    chi2 += (d / error_neg) ** 2 if d < 0 else (d / error_pos) ** 2
-                self._d[frame_i] = (rDA, rDAE, rMP, chi2)
-
+            rMP, rDA, rDAE, chi2 = self.get_frame_results(frame_i)
             for distance_key in self.distances:
                 re[distance_key]['rMP'].append(rMP)
                 re[distance_key]['rDA'].append(rDA)
                 re[distance_key]['rDAE'].append(rDAE)
                 re[distance_key]['chi2'].append(chi2)
-
-        return re
+        return re

conda-recipe/bld.bat

@@ -0,0 +1 @@
+"%PYTHON%" setup.py install --single-version-externally-managed --record=record.txt

conda-recipe/build.sh

@@ -0,0 +1,10 @@
+#!/usr/bin/env bash
+# build documentation
+#sphinx-apidoc -o doc avtraj
+#cd doc
+#make html
+#cd ..
+
+$PYTHON setup.py install --single-version-externally-managed --record=record.txt  # Python command to install the script.
+#python setup.py sdist
+#twine upload dist/*

conda-recipe/conda_build_config.yaml

@@ -0,0 +1,3 @@
+python: 
+  - 3.6
+  - 3.7