Merge pull request #73 from pele-python/store_trajectory

Store trajectory: This has been sitting here for a while, and I have been using it without any issues, so I am merging

mcpele/monte_carlo/__init__.py

@@ -1,11 +1,13 @@
 from _accept_test_cpp import MetropolisTest
 from _conf_test_cpp import CheckSphericalContainer
-from _conf_test_cpp import CheckSphericalContainerConfig
+from _conf_test_cpp import CheckSphericalContainerConfig, ConfTestOR
 from _takestep_cpp import RandomCoordsDisplacement, SampleGaussian
 from _takestep_cpp import GaussianCoordsDisplacement, ParticlePairSwap
 from _takestep_cpp import TakeStepPattern, TakeStepProbabilities
 from _takestep_cpp import UniformSphericalSampling, UniformCubicSampling
 from _monte_carlo_cpp import _BaseMCRunner
 from _action_cpp import RecordEnergyHistogram, RecordEnergyTimeseries, RecordPairDistHistogram, RecordLowestEValueTimeseries, RecordDisplacementPerParticleTimeseries
+from _action_cpp import RecordCoordsTimeseries
 from _nullpotential_cpp import NullPotential
 from mcrunner import Metropolis_MCrunner
+

mcpele/monte_carlo/_action_cpp.pxd

@@ -2,6 +2,7 @@ cimport pele.potentials._pele as _pele
 #from pele.potentials._pele cimport array_wrap_np
 from _pele_mc cimport cppAction,_Cdef_Action, shared_ptr
 from libcpp cimport bool as cbool
+from libcpp.deque cimport deque
 
 # cython has no support for integer template argument.  This is a hack to get around it
 # https://groups.google.com/forum/#!topic/cython-users/xAZxdCFw6Xs
@@ -35,7 +36,7 @@ cdef extern from "mcpele/record_scalar_timeseries.h" namespace "mcpele":
         _pele.Array[double] get_time_series() except +
         void clear() except +
         cbool moving_average_is_stable(size_t, double) except +
-        
+
 cdef extern from "mcpele/record_energy_timeseries.h" namespace "mcpele":
     cdef cppclass cppRecordEnergyTimeseries "mcpele::RecordEnergyTimeseries":
         cppRecordEnergyTimeseries(size_t, size_t) except +
@@ -50,4 +51,17 @@ cdef extern from "mcpele/record_displacement_per_particle_timeseries.h" namespac
     cdef cppclass cppRecordDisplacementPerParticleTimeseries "mcpele::RecordDisplacementPerParticleTimeseries":
         cppRecordDisplacementPerParticleTimeseries(size_t, size_t,
             _pele.Array[double], size_t) except +
-
+
+cdef extern from "mcpele/record_vector_timeseries.h" namespace "mcpele":
+    cdef cppclass cppRecordVectorTimeseries "mcpele::RecordVectorTimeseries":
+        deque[_pele.Array[double]] get_time_series() except +
+        void clear() except +
+        size_t get_record_every() except +
+
+cdef extern from "mcpele/record_coords_timeseries.h" namespace "mcpele":
+    cdef cppclass cppRecordCoordsTimeseries "mcpele::RecordCoordsTimeseries":
+        cppRecordCoordsTimeseries(size_t, size_t, size_t) except +
+        _pele.Array[double] get_mean_coordinate_vector() except +
+        _pele.Array[double] get_mean2_coordinate_vector() except +
+        _pele.Array[double] get_variance_coordinate_vector() except +
+        size_t get_count() except +

mcpele/monte_carlo/_action_cpp.pyx

@@ -380,3 +380,89 @@ class RecordDisplacementPerParticleTimeseries(_Cdef_RecordDisplacementPerParticl
     boxdimension: int
         dimensionality of the space (dimensionality of box)
     """
+
+cdef class _Cdef_RecordCoordsTimeseries(_Cdef_Action):
+    cdef cppRecordVectorTimeseries* newptr
+    cdef cppRecordCoordsTimeseries* newptr2
+    def __cinit__(self, ndof, record_every=1, eqsteps=0):
+        self.thisptr = shared_ptr[cppAction](<cppAction*> new cppRecordCoordsTimeseries(ndof, record_every, eqsteps))
+        self.newptr = <cppRecordVectorTimeseries*> self.thisptr.get()
+        self.newptr2 = <cppRecordCoordsTimeseries*> self.thisptr.get()
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    def get_time_series(self):
+        """get a trajectory
+        
+        Returns
+        -------
+        numpy.array
+            deque of arrays
+        """
+        cdef deque[_pele.Array[double]] dq = self.newptr.get_time_series()
+        cdef double *seriesdata
+        cdef np.ndarray[double, ndim=2, mode="c"] series = np.zeros((dq.size(), dq[0].size()))
+        cdef size_t i, j
+        for i in xrange(dq.size()):
+            seriesdata = dq[i].data()
+            for j in xrange(dq[i].size()):
+                series[i][j] = seriesdata[j]
+        return series
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    def get_mean_variance_time_series(self):
+        """returns the mean coordinate vector
+        
+        Returns
+        -------
+        numpy.array
+            mean coordinate vector
+        numpy.array
+            element-wise variance coordinate vector
+        """
+        cdef _pele.Array[double] coordi = self.newptr2.get_mean_coordinate_vector()
+        cdef _pele.Array[double] vari = self.newptr2.get_variance_coordinate_vector()
+        cdef double *coorddata = coordi.data()
+        cdef double *vardata = vari.data()
+        cdef np.ndarray[double, ndim=1, mode="c"] coord = np.zeros(coordi.size())
+        cdef np.ndarray[double, ndim=1, mode="c"] var = np.zeros(vari.size())
+        cdef size_t i
+        for i in xrange(coordi.size()):
+            coord[i] = coorddata[i]
+            var[i] = vardata[i]
+        return coord, var
+
+    def get_avg_count(self):
+        """get number of steps over which the average and variance were computed
+        
+        Returns
+        -------
+        count : integer
+            sample size
+        """
+        count = self.newptr2.get_count()
+        return count
+
+    def clear(self):
+        """clear time series container
+        
+        deletes the entries in the c++ container
+        """
+        self.newptr.clear()
+
+class RecordCoordsTimeseries(_Cdef_RecordCoordsTimeseries):
+    """Record a trajectory of the system coordinates
+       
+    This class is the Python interface for the c++ bv::RecordCoordsTimeseries 
+    :class:`Action` class implementation.
+    
+    Parameters
+    ----------
+    ndof : int
+        dimensionality of coordinate array
+    record_every : int
+        interval every which the coordinates are recorded
+    eqsteps : int
+        number of equilibration steps to skip when computing averages
+    """

mcpele/monte_carlo/_conf_test_cpp.pxd

@@ -1,5 +1,5 @@
 from libcpp cimport bool as cbool
-from _pele_mc cimport cppConfTest,_Cdef_ConfTest, shared_ptr
+from _pele_mc cimport cppConfTest, _Cdef_ConfTest, shared_ptr
 
 cdef extern from "mcpele/check_spherical_container.h" namespace "mcpele":
     cdef cppclass cppCheckSphericalContainer "mcpele::CheckSphericalContainer":
@@ -8,3 +8,8 @@ cdef extern from "mcpele/check_spherical_container.h" namespace "mcpele":
 cdef extern from "mcpele/check_spherical_container_config.h" namespace "mcpele":
     cdef cppclass cppCheckSphericalContainerConfig "mcpele::CheckSphericalContainerConfig":
         cppCheckSphericalContainerConfig(double) except +
+
+cdef extern from "mcpele/conf_test_OR.h" namespace "mcpele":
+    cdef cppclass cppConfTestOR "mcpele::ConfTestOR":
+        cppConfTestOR() except +
+        void add_test(shared_ptr[cppConfTest]) except +

mcpele/monte_carlo/_conf_test_cpp.pyx

@@ -46,3 +46,33 @@ class CheckSphericalContainerConfig(_Cdef_CheckSphericalContainerConfig):
     radius : double
         radius of the spherical container, centered at **0**
     """
+
+#===============================================================================
+# Union of configurational tests
+#===============================================================================
+
+cdef class _Cdef_ConfTestOR(_Cdef_ConfTest):
+    cdef cppConfTestOR* newptr
+    def __cinit__(self):
+        self.thisptr = shared_ptr[cppConfTest](<cppConfTest*> new cppConfTestOR())
+        self.newptr = <cppConfTestOR*> self.thisptr.get()
+
+    def add_test(self, _Cdef_ConfTest test):
+        """add a conf test to the union
+        
+        Parameters
+        ----------
+        step : :class:`ConfTest`
+            object of class :class:`ConfTest` constructed beforehand
+        """
+        self.newptr.add_test(test.thisptr)
+
+class ConfTestOR(_Cdef_ConfTestOR):
+    """Creates a union of multiple configurational tests
+    
+    Python interface for c++ ConfTestOR. This configurational
+    test takes multiple conf tests and generates an union,
+    therefore it is sufficient to satisfy one of these
+    subtests to pass their union.
+    
+    """

mcpele/utils/__init__.py

@@ -0,0 +1 @@
+from _utils import *

mcpele/utils/_utils.py

@@ -0,0 +1,23 @@
+from __future__ import division
+import numpy as np
+import pandas as pd
+
+def write_2d_array_to_hf5(array, key, path):
+    """
+    this function can be used to dump any 2d array to a hf5 database
+    use this to dump a trajectory to a database from python
+    """
+    assert array.ndim == 2
+    nind , ncol = array.shape
+    ind = [i for i in xrange(nind)]
+    col = [i for i in xrange(ncol)]
+    df = pd.DataFrame(np.array(array), index=ind, columns=col)
+    df.to_hdf(path, key)
+
+def read_hf5_to_2d_array(path, key):
+    """
+    read a 2d array from a hf5 database
+    """
+    df = pd.read_hdf(path, key)
+    array = np.array(df.values)
+    return array

source/conf_test_OR.cpp

@@ -0,0 +1,27 @@
+#include "mcpele/conf_test_OR.h"
+
+namespace mcpele {
+
+ConfTestOR::ConfTestOR(){}
+
+void ConfTestOR::add_test(std::shared_ptr<ConfTest> test_input)
+{
+    m_tests.push_back(test_input);
+    m_tests.swap(m_tests);
+}
+
+bool ConfTestOR::conf_test(pele::Array<double> &trial_coords, MC * mc)
+{
+    if (m_tests.size() == 0) {
+        throw std::runtime_error("ConfTestOR::conf_test: no conf test specified");
+    }
+    for(auto & test : m_tests){
+        bool result = test->conf_test(trial_coords, mc);
+        if (result){
+            return true;
+        }
+    }
+    return false;
+}
+
+} // namespace mcpele

source/mcpele/conf_test_OR.h

@@ -0,0 +1,28 @@
+#ifndef _MCPELE_CONF_TEST_OR_H__
+#define _MCPELE_CONF_TEST_OR_H__
+
+#include <vector>
+#include <random>
+
+#include "mc.h"
+
+namespace mcpele {
+
+/**
+ * Create union of two configurational tests,
+ * it is sufficient for one of them to be true in order to pass the overall test.
+ *
+ */
+class ConfTestOR : public ConfTest {
+private:
+    std::vector<std::shared_ptr<ConfTest> > m_tests;
+public:
+    virtual ~ConfTestOR(){}
+    ConfTestOR();
+    void add_test(std::shared_ptr<ConfTest> test_input);
+    bool conf_test(pele::Array<double> &trial_coords, MC * mc);
+};
+
+} // namespace mcpele
+
+#endif // #ifndef _MCPELE_CONF_TEST_OR_H__

source/mcpele/record_coords_timeseries.h

@@ -0,0 +1,36 @@
+#ifndef _MCPELE_RECORD_COORDS_TIMESERIES_H__
+#define _MCPELE_RECORD_COORDS_TIMESERIES_H__
+
+#include "record_vector_timeseries.h"
+
+namespace mcpele {
+
+class RecordCoordsTimeseries : public RecordVectorTimeseries {
+private:
+    pele::Array<double> m_mcv, m_mcv2;
+    const size_t m_ndof;
+    size_t m_count;
+    double m_update_average(double avg, double x);
+    void m_update_mean_coord_vector(pele::Array<double> &new_coords);
+public:
+    RecordCoordsTimeseries(const size_t ndof, const size_t record_every, const size_t eqsteps);
+    virtual ~RecordCoordsTimeseries() {}
+    virtual pele::Array<double> get_recorded_vector(pele::Array<double> &coords,
+          const double energy, const bool accepted, MC* mc) { return coords; }
+    virtual void action(pele::Array<double> &coords, double energy, bool accepted, MC* mc);
+    pele::Array<double> get_mean_coordinate_vector(){return m_mcv.copy();}
+    pele::Array<double> get_mean2_coordinate_vector(){return m_mcv2.copy();}
+    pele::Array<double> get_variance_coordinate_vector()
+    {
+        pele::Array<double> var = m_mcv2.copy();
+        for(size_t i=0; i<m_ndof; ++i){
+            var[i] -= m_mcv[i]*m_mcv[i];
+        }
+        return var.copy();
+    }
+    size_t get_count(){return m_count;}
+};
+
+} // namespace mcpele
+
+#endif // #ifndef _MCPELE_RECORD_COORDS_TIMESERIES_H__

source/mcpele/record_vector_timeseries.h

@@ -0,0 +1,43 @@
+#ifndef _MCPELE_RECORD_VECTOR_TIMESERIES_H__
+#define _MCPELE_RECORD_VECTOR_TIMESERIES_H__
+
+#include "mc.h"
+#include <deque>
+#include <cstdlib>
+
+namespace mcpele {
+
+/**
+ * Record vector time series, every record_every-th step.
+ */
+class RecordVectorTimeseries : public Action {
+protected:
+    const size_t m_record_every, m_eqsteps;
+    std::deque<pele::Array<double>> m_time_series;
+    void m_record_vector_value(pele::Array<double> input)
+    {
+        try{
+            m_time_series.push_back(input.copy());
+        }
+        catch(std::bad_alloc &ba){
+            std::cerr<< "mcpele::RecordVectorTimeseries: bad_alloc caught: " << ba.what() << std::endl;
+            std::exit(EXIT_FAILURE);
+        }
+    }
+public:
+    RecordVectorTimeseries(const size_t record_every, const size_t eqsteps);
+    virtual ~RecordVectorTimeseries(){}
+    virtual void action(pele::Array<double> &coords, double energy, bool accepted, MC* mc);
+    virtual pele::Array<double> get_recorded_vector(pele::Array<double> &coords, const double energy, const bool accepted, MC* mc)=0;
+    std::deque<pele::Array<double>> get_time_series()
+    {
+        m_time_series.shrink_to_fit();
+        return m_time_series;
+    }
+    void clear() { m_time_series.clear(); }
+    size_t get_record_every(){return m_record_every;}
+};
+
+} // namespace mcpele
+
+#endif // #ifndef _MCPELE_RECORD_VECTOR_TIMESERIES_H__