Merge pull request #7 from Countoscope/new_structure_and_names

New class structure and names, and example

.gitignore

@@ -5,4 +5,5 @@
 
 # pip install folders
 countoscope.egg-info/
-build/
+build/
+examples/example_xyz.png

README.rst

@@ -11,26 +11,19 @@
 .. image:: https://github.com/Countoscope/countoscope/actions/workflows/python-app.yml/badge.svg
    :target: https://github.com/Countoscope/countoscope/actions/workflows/python-app.yml
 
-Py-Countoscope
-##############
+Countoscope (Python)
+####################
 
-This is the repository containing the Python code of the countoscope.
+This is the repository containing the Python code of the Countoscope.
 
+Installation (users)
+--------------------
 Clone the repository using:
 
 ..  code:: bash
 
   git clone https://github.com/Countoscope/countoscope.git
 
-Optionally, clone the repository with its submodule using:
-
-..  code:: bash
-
-  git clone https://github.com/Countoscope/countoscope.git --recurse-submodules
-
-Installation
-------------
-
 Install all the required Python packages using:
 
 ..  code:: bash
@@ -43,66 +36,32 @@ Then, type:
 
     pip install .
 
-Use
----
-
-In a Python script or a Jupyter Notebook, define the path to the data trajectory
-file. For instance, using the `3D-closed` dataset provided in the `datasets`
-repository:
-
-..  code:: python
-
-  path_to_data = "/mpath/datasets/datasets/3D-closed/trajectory.xyz"
-
-Import the Countoscope as well as NumPy by typing:
-
-..  code:: python
-
-  from countoscope import Countoscope
-  import numpy as np
-
-The trajectory file `trajectory.xyz` corresponds to a system of 190 particles in
-a :math:`(30 Å)^3` box. Let us define the system size as a NumPy array:
+Installation (developers)
+-------------------------
+Clone the repository with its submodule (containing datasets for testing) using:
 
-..  code:: python
-
-    system_size = np.array([30, 30, 30])
-
-Finally, let us choose a grid size for the Countoscope measurement:
-
-..  code:: python
-
-    box_size=np.array([10, 10, 10])
-
-Then, launch the Countoscope calculation using *trajectory_file*, *system_size*,
-and *box_size* as input parameters:
+..  code:: bash
 
-..  code:: python
+  git clone https://github.com/Countoscope/countoscope.git --recurse-submodules
 
-    results = Countoscope(trajectory_file = path_to_data,
-                        system_size=system_size,
-                        box_size=box_size)
-    results.run()
+You will need Git LFS installed (https://git-lfs.com/)
 
-After the calculation is done, all the computed data can be obtained from the
-`results`object. For instance, for :math:`<N>`, type:
+Install all the required Python packages using:
 
-..  code:: python
+..  code:: bash
 
-    print(np.round(results.mean_of_N,2))
+    pip install -r requirements.txt
 
-which will return:
+Then, type:
 
 ..  code:: bash
 
-    0.84
+    pip install -e .
 
-To plot :math:`<\Delta N^2>`, let us import Pyplot first:
-
-..  code:: python
+Use
+---
 
-    import matplotlib.pyplot as plt
-    plt.loglog(results.delta_n2)
+See examples/example_xyz.py
 
 Run the tests
 -------------

countoscope/__init__.py

@@ -7,3 +7,4 @@
 # Released under MIT Licence
 
 from .countoscope import *
+from .utilities import convert_homemade_format

countoscope/box.py → countoscope/boxcountingtool.py

@@ -5,29 +5,27 @@
 # (see the AUTHORS.rst file for the full list of names)
 #
 # Released under MIT Licence
-"""Module for preparing and cutting the box."""
-from .trajectory import Trajectory
+"""Module for counting the particles into the boxes."""
+from .trajectorytool import TrajectoryTool
 import numpy as np
 
 
-class Box(Trajectory):
-    r"""Class providing options for dividing the system into small boxes.
-    
-    The box size can be provided as an integer or an array of same 
-    size as the system dimension. 
-    """
-
+class BoxCountingTool(TrajectoryTool):
+    """ Class for tools based on counting particles in boxes """
     def __init__(self,
                  box_size: None,
+                 start = None,
+                 stop = None,
                  *args,
                  **kwargs
                  ):
         super().__init__(*args, **kwargs)
         self.box_size = box_size
-        self.run_box()
+        self.start = start
+        self.stop = stop
+        self.prepare_boxes()
 
-    def run_box(self):
-        self.run_trajectory()
+    def prepare_boxes(self):
         self.specify_measurement_box_size()
         self.select_the_slices()
         self.allocate_memory()
@@ -81,4 +79,50 @@ def allocate_memory(self):
                                       len(self.box_bounds[0])-1,
                                       len(self.box_bounds[1])-1,
                                       len(self.box_bounds[2])-1))
-        self.number_matrix = number_matrix
+        self.number_matrix = number_matrix
+
+    def run_counting(self):
+        self.prepare_run()
+        self.store_particle_into_boxes()
+
+    def prepare_run(self):
+        """Detect starting and stopping frames for analysis."""
+        if self.start is None:
+            self.start = 0
+        if self.stop is None:
+            self.stop = self.nb_steps
+
+    def store_particle_into_boxes(self):
+        """Run over the trajectory and place the particles into number_matrix"""
+        for i in np.arange(self.start, self.stop):
+            frame = self.trajectory.read_step(i)
+            #positions = frame.positions
+            #if self.variable == "x":
+            #    weights = None
+            #elif self.variable == "v":
+            #    weights = np.linalg.norm(self.atom_group.velocities, axis=1)
+            #elif self.variable == "f":
+            #    weights = np.linalg.norm(self.atom_group.forces, axis=1)
+            # todo: allow atom selection
+            if self.dimension == 2:
+                H, _ = np.histogramdd(frame.positions[:,0:2],
+                                      bins = (len(self.box_bounds[0])-1,
+                                              len(self.box_bounds[1])-1,),
+                                      range = ((np.min(self.box_bounds[0]),
+                                                np.max(self.box_bounds[0])),
+                                                (np.min(self.box_bounds[1]),
+                                                 np.max(self.box_bounds[1]))))
+                self.number_matrix[i, :, :] = H
+            elif self.dimension == 3:
+                H, _ = np.histogramdd(frame.positions,
+                                    bins = (len(self.box_bounds[0])-1,
+                                            len(self.box_bounds[1])-1,
+                                            len(self.box_bounds[2])-1),
+                                    range = ((np.min(self.box_bounds[0]),
+                                            np.max(self.box_bounds[0])),
+                                            (np.min(self.box_bounds[1]),
+                                            np.max(self.box_bounds[1])),
+                                            (np.min(self.box_bounds[2]),
+                                            np.max(self.box_bounds[2]))))
+                                    #weights=weights)
+                self.number_matrix[i, :, :, :] = H

countoscope/countoscope.py

@@ -6,22 +6,20 @@
 #
 # Released under MIT Licence
 """Module for sorting the particles in boxes."""
-from .sorting import Sorting
+from .boxcountingtool import BoxCountingTool
 from .utilities import autocorrelation_function, calculate_remaining_axis
 import numpy as np
 
-
-class Countoscope(Sorting):
+class Countoscope(BoxCountingTool):
     def __init__(self,
                  *args,
                  **kwargs
                  ):
         super().__init__(*args, **kwargs)
 
-    def run(self):
-        self.evaluate_correlation()
+    def count(self):
+        self.run_counting() # from BoxCountingTool
         self.evaluate_mean_particle_number()
-        self.evaluate_deltan2()
 
     def evaluate_correlation(self):
         """Evaluate the correlation function."""
@@ -45,6 +43,7 @@ def evaluate_correlation(self):
                                             axis = remaining_axis)
         self.err_correlation_function = np.std(per_box_correlation_functions,
                                                axis = remaining_axis)/np.sqrt(nb_boxes)
+        return self.correlation_function, self.err_correlation_function
 
 
     def evaluate_mean_particle_number(self):
@@ -81,7 +80,12 @@ def evaluate_mean_particle_number(self):
         self.mean_of_N_squared = np.mean(self.per_box_mean_of_N_squared,
                                          axis = remaining_axis)
 
-    def evaluate_deltan2(self):
+    def evaluate_deltaN2(self):
         """Evaluate <deltan2>"""
+        if not hasattr(self, 'correlation_function'):
+            # in case we have not yet computed the correlation function, do it now
+            self.evaluate_correlation()
+
         self.delta_n2 = 2*(self.mean_of_square_of_N-self.mean_of_N_squared) \
             - 2*(self.correlation_function - self.mean_of_N_squared)
+        return self.delta_n2