Merge pull request #71 from kjs73/uniform_sampling

Uniform sampling

cpp_tests/source/test_takestep.cpp

@@ -6,6 +6,7 @@
 #include <gtest/gtest.h>
 
 #include "pele/harmonic.h"
+#include "pele/meta_pow.h"
 
 #include "mcpele/adaptive_takestep.h"
 #include "mcpele/histogram.h"
@@ -14,6 +15,8 @@
 #include "mcpele/random_coords_displacement.h"
 #include "mcpele/take_step_pattern.h"
 #include "mcpele/take_step_probabilities.h"
+#include "mcpele/uniform_cubic_sampling.h"
+#include "mcpele/uniform_spherical_sampling.h"
 
 using pele::Array;
 using mcpele::MC;
@@ -123,6 +126,57 @@ TEST_F(TakeStepTest, Single_BasicFunctionality_AllParticlesSampledUniformly){
     EXPECT_NEAR_RELATIVE(hist_uniform_single.get_variance(), (ntot * ntot - 1) / 12, (ntot * ntot - 1) / (12 * sqrt(ntot)));
 }
 
+TEST_F(TakeStepTest, UniformCubic_CorrectMoments){
+    // test that first two moments of uniform cubic are correct
+    // http://mathworld.wolfram.com/UniformDistribution.html
+    const size_t ndim = 20;
+    const size_t nsamples = 1e4;
+    const double delta = 42.42;
+    std::vector<std::shared_ptr<mcpele::Histogram> > hist(100, std::make_shared<mcpele::Histogram>(0, nparticles - 1, 1));
+    mcpele::UniformCubicSampling sampler(42, delta);
+    pele::Array<double> x(ndim);
+    for (size_t i = 0; i < nsamples; ++i) {
+        sampler.displace(x, NULL);
+        for (size_t k = 0; k < ndim; ++k) {
+            hist[k]->add_entry(x[k]);
+        }
+    }
+    for (size_t i = 0; i < ndim; ++i) {
+        EXPECT_NEAR_RELATIVE(hist[i]->get_mean(), 0, 2 / sqrt(nsamples));
+        EXPECT_NEAR_RELATIVE(hist[i]->get_variance(), delta * delta / 3, 2 / sqrt(nsamples));
+    }
+}
+
+TEST_F(TakeStepTest, UniformSpherical_CorrectMoments){
+    // test uniform spherical is OK in 2d
+    const size_t nsamples = 1e5;
+    const double radius = 42.42;
+    mcpele::Histogram hist(0, nparticles - 1, 1);
+    mcpele::Histogram hist3(0, nparticles - 1, 1);
+    mcpele::UniformSphericalSampling sampler(42, radius);
+    pele::Array<double> x2(2);
+    pele::Array<double> x3(3);
+    for (size_t i = 0; i < nsamples; ++i) {
+        sampler.displace(x2, NULL);
+        hist.add_entry(pele::dot(x2, x2));
+        sampler.displace(x3, NULL);
+        hist3.add_entry(pele::dot(x3, x3));
+    }
+    /**For a random walk in a disk of radius R, the mean of the squared
+     * displacement from the origin is R^/2 and the variance of the
+     * squared displacement is R^4/12.
+     */
+    EXPECT_NEAR_RELATIVE(radius * radius / 2, hist.get_mean(), 1e-3);
+    EXPECT_NEAR_RELATIVE(pele::pos_int_pow<4>(radius) / 12, hist.get_variance(), 1e-2);
+    /**For a random walk in a 3d sphere of radius R, the mean of the
+     * squared displacement from the origin is 3 * R^2 / 5 and the
+     * variance of the squared displacement from the origin is
+     * 12 * R^4 / 175.
+     */
+     EXPECT_NEAR_RELATIVE(3 * pele::pos_int_pow<2>(radius) / 5, hist3.get_mean(), 1e-3);
+     EXPECT_NEAR_RELATIVE(12 * pele::pos_int_pow<4>(radius) / 175, hist3.get_variance(), 1e-2);
+}
+
 TEST_F(TakeStepTest, Single_BasicFunctionalityAveragingErasing_NIterations){
     //n iterations give expected variation
     mcpele::RandomCoordsDisplacementSingle displ(seed, nparticles, ndim, stepsize);

cpp_tests/source/test_time_series.cpp

@@ -109,7 +109,7 @@ TEST(EVTimeseries, Works) {
     EXPECT_EQ(series.size(), niter / record_every);
     const double eigenvalue_reference = series[0];
     for (size_t i = 0; i < series.size(); ++i) {
-        EXPECT_NEAR_RELATIVE(series[i], eigenvalue_reference, 1e-10);
+        EXPECT_NEAR_RELATIVE(series[i], eigenvalue_reference, 1e-9);
     }
 }
 

examples/compute_pi/compute_pi.py

@@ -0,0 +1,57 @@
+"""
+This example computes the numerical value of Pi by Monte Carlo.
+Points are sampled uniformly at random from an n-dimensional cube of
+side length 2R. By measuring the fraction of points contained within
+a ball of radius R contained in the cube, we determine Pi.
+"""
+from __future__ import division
+import numpy as np
+import copy
+from scipy.special import gamma
+from mcpele.monte_carlo import _BaseMCRunner
+from mcpele.monte_carlo import NullPotential
+from mcpele.monte_carlo import UniformCubicSampling
+from mcpele.monte_carlo import CheckSphericalContainerConfig
+
+def volume_nball(radius, n):
+    return np.power(np.pi, n / 2) * np.power(radius, n) / gamma(n / 2 + 1)
+
+def get_pi(accepted_fraction, ndim):
+    return np.power(2 ** ndim * accepted_fraction * gamma(ndim / 2 + 1), 2 / ndim)
+
+class MC(_BaseMCRunner):
+    def set_control(self, temp):
+        self.set_temperature(temp)
+
+class ComputePi(object):
+    def __init__(self, ndim=2, nsamples=1e4):
+        #
+        self.ndim = ndim
+        self.nsamples = nsamples
+        #
+        self.radius = 44
+        self.potential = NullPotential()
+        self.mc = MC(self.potential, np.ones(self.ndim), 1, self.nsamples)
+        self.step = UniformCubicSampling(42, self.radius)
+        self.mc.set_takestep(self.step)
+        self.mc.set_report_steps(0)
+        self.conftest_check_spherical_container = CheckSphericalContainerConfig(self.radius)
+        self.mc.add_conf_test(self.conftest_check_spherical_container)
+        self.mc.set_print_progress()
+        self.mc.run()
+        self.p = self.mc.get_accepted_fraction()
+        self.pi = get_pi(self.p, self.ndim)
+
+if __name__ == "__main__":
+    nsamples = 1e5
+    ndim_ = []
+    res = []
+    for ndim in xrange(2, 16):
+        print("computing pi in {} dimensions".format(ndim))
+        c = ComputePi(ndim=ndim, nsamples=nsamples)
+        res.append(c.pi)
+        ndim_.append(ndim)
+    for (i, p) in zip(ndim_, res):
+        print("dimension", i)
+        print("pi", p)
+        print("pi / np.pi", p / np.pi)

mcpele/monte_carlo/__init__.py

@@ -1,7 +1,10 @@
 from _accept_test_cpp import MetropolisTest
 from _conf_test_cpp import CheckSphericalContainer
 from _conf_test_cpp import CheckSphericalContainerConfig
-from _takestep_cpp import RandomCoordsDisplacement, SampleGaussian, GaussianCoordsDisplacement, ParticlePairSwap, TakeStepPattern, TakeStepProbabilities
+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 _nullpotential_cpp import NullPotential

mcpele/monte_carlo/_takestep_cpp.pxd

@@ -18,6 +18,16 @@ cdef extern from "mcpele/random_coords_displacement.h" namespace "mcpele":
         size_t get_seed() except +
         void set_generator_seed(size_t) except +
         double get_stepsize() except +
+
+cdef extern from "mcpele/uniform_spherical_sampling.h" namespace "mcpele":
+    cdef cppclass cppUniformSphericalSampling "mcpele::UniformSphericalSampling":
+        cppUniformSphericalSampling(size_t, double) except +
+        void set_generator_seed(size_t) except +
+
+cdef extern from "mcpele/uniform_cubic_sampling.h" namespace "mcpele":
+    cdef cppclass cppUniformCubicSampling "mcpele::UniformCubicSampling":
+        cppUniformCubicSampling(size_t, double) except +
+        void set_generator_seed(size_t) except +
 
 cdef extern from "mcpele/gaussian_coords_displacement.h" namespace "mcpele":
     cdef cppclass cppGaussianTakeStep "mcpele::GaussianTakeStep":

mcpele/monte_carlo/_takestep_cpp.pyx

@@ -96,6 +96,73 @@ class RandomCoordsDisplacement(_Cdef_RandomCoordsDisplacement):
         dimensionality of the space (box dimensionality)
     """
 
+#===============================================================================
+# UniformSphericalSampling
+#===============================================================================
+
+cdef class _Cdef_UniformSphericalSampling(_Cdef_TakeStep):
+    cdef cppUniformSphericalSampling* newptr
+    def __cinit__(self, rseed=42, radius=1):
+        self.thisptr = shared_ptr[cppTakeStep](<cppTakeStep*> new cppUniformSphericalSampling(rseed, radius))
+        self.newptr = <cppUniformSphericalSampling*> self.thisptr.get()
+    def set_generator_seed(self, input):
+        """sets the random number generator seed
+        
+        Parameters
+        ----------
+        input : pos int
+            random number generator seed
+        """
+        cdef inp = input
+        self.newptr.set_generator_seed(inp)
+
+class UniformSphericalSampling(_Cdef_UniformSphericalSampling):
+    """Sample uniformly at random inside N-ball.
+    
+    Implements the method described here:
+    http://math.stackexchange.com/questions/87230/picking-random-points-in-the-volume-of-sphere-with-uniform-probability
+    Variates $X_1$ to $X_N$ are sampled independently from a standard
+    normal distribution and then rescaled as described in the reference.
+    See also, e.g. Numerical Recipes, 3rd ed, page 1129.
+    
+    Parameters
+    ----------
+    rseed : pos int
+        seed for the random number generator (std:library 64 bits Merseene Twister)
+    radius : double
+        radius of ball
+    """
+
+#===============================================================================
+# UniformCubicSampling
+#===============================================================================
+
+cdef class _Cdef_UniformCubicSampling(_Cdef_TakeStep):
+    cdef cppUniformCubicSampling* newptr
+    def __cinit__(self, rseed=42, delta=1):
+        self.thisptr = shared_ptr[cppTakeStep](<cppTakeStep*> new cppUniformCubicSampling(rseed, delta))
+        self.newptr = <cppUniformCubicSampling*> self.thisptr.get()
+    def set_generator_seed(self, input):
+        """sets the random number generator seed
+        
+        Parameters
+        ----------
+        input : pos int
+            random number generator seed
+        """
+        cdef inp = input
+        self.newptr.set_generator_seed(inp)
+
+class UniformCubicSampling(_Cdef_UniformCubicSampling):
+    """Sample uniformly at random inside cube.
+    Parameters
+    ----------
+    rseed : pos int
+        seed for the random number generator (std:library 64 bits Merseene Twister)
+    delta : double
+        half side length of cube
+    """
+
 #===============================================================================
 # GaussianCoordsDisplacement
 #===============================================================================

source/mcpele/uniform_cubic_sampling.h

@@ -0,0 +1,27 @@
+#ifndef _MCPELE_UNIFORM_CUBIC_SAMPLING_H__
+#define _MCPELE_UNIFORM_CUBIC_SAMPLING_H__
+
+namespace mcpele {
+
+class UniformCubicSampling : public TakeStep {
+protected:
+    std::mt19937_64 m_gen;
+    std::uniform_real_distribution<double> m_dist;
+public:
+    virtual ~UniformCubicSampling() {}
+    UniformCubicSampling(const size_t seed=42, const double delta=1)
+        : m_gen(seed),
+          m_dist(-delta, delta)
+    {}
+    void set_generator_seed(const size_t inp) { m_gen.seed(inp); }
+    virtual void displace(pele::Array<double>& coords, MC* mc)
+    {
+        for (size_t i = 0; i < coords.size(); ++i) {
+            coords[i] = m_dist(m_gen);
+        }
+    }
+};    
+
+} // namespace mcpele
+
+#endif // #ifndef _MCPELE_UNIFORM_CUBIC_SAMPLING_H__

source/mcpele/uniform_spherical_sampling.h

@@ -0,0 +1,51 @@
+#ifndef _MCPELE_UNIFORM_SPHERICAL_SAMPLING_H__
+#define _MCPELE_UNIFORM_SPHERICAL_SAMPLING_H__
+
+namespace mcpele {
+
+/**
+ * Sample points uniformly at random within an N-ball.
+ * See also: http://math.stackexchange.com/questions/87230/picking-random-points-in-the-volume-of-sphere-with-uniform-probability
+ */    
+
+class UniformSphericalSampling : public TakeStep {
+protected:
+    std::mt19937_64 m_gen;
+    const double m_radius;    
+    std::normal_distribution<double> m_dist_normal;
+    std::uniform_real_distribution<double> m_dist_uniform;
+public:
+    virtual ~UniformSphericalSampling() {}
+    UniformSphericalSampling(const size_t seed=42, const double radius=1)
+        : m_gen(seed),
+          m_radius(radius),
+          m_dist_normal(0, 1),
+          m_dist_uniform(0, 1)
+    {}
+    void set_generator_seed(const size_t inp) { m_gen.seed(inp); }
+    virtual void displace(pele::Array<double>& coords, MC* mc)
+    {
+        for (size_t i = 0; i < coords.size(); ++i) {
+            coords[i] = m_dist_normal(m_gen);
+        }
+        /**
+         * From Numerical Recipes:
+         * Picking a random point on a sphere:
+         * 1) generate n independent, identically distributed, normal random numbers, y_0, ..., y_{n-1}
+         * 2) get point {x} on unit sphere in n dimensions by {x} = {y} / norm({y})
+         * Picking a random point inside a sphere:
+         * 3) generate an additional uniform random number u in [0,1]
+         * 4) compute point x_i = y_i * u^{1/n} / norm({y})
+         */
+        // This computes 1 / norm({y}).
+        double tmp = 1.0 / norm(coords);
+        // This computes u^{1/n} / norm({y}) and rescales to a sphere of radius m_radius.
+        tmp *= m_radius * std::pow(m_dist_uniform(m_gen), 1.0 / coords.size());
+        // This computes the sampled random point in the sphere.
+        coords *= tmp;
+    }
+};
+
+} // namespace mcpele
+
+#endif // #ifndef _MCPELE_UNIFORM_SPHERICAL_SAMPLING_H__