[Boost-commit] svn:boost r64521 - in sandbox/SOC/2010/quasi_random: boost/random libs/random/example libs/random/example/detail

Subject: [Boost-commit] svn:boost r64521 - in sandbox/SOC/2010/quasi_random: boost/random libs/random/example libs/random/example/detail
From: jvd_at_[hidden]
Date: 2010-08-01 12:16:46

Author: qrng
Date: 2010-08-01 12:16:42 EDT (Sun, 01 Aug 2010)
New Revision: 64521
URL: http://svn.boost.org/trac/boost/changeset/64521

Log:
Added the Cauchy-Crofton example that computes the surface area of a box, also a convergence graph, which 
clearly displays that a quasi-Monte Carlo really converges faster.

Added:
   sandbox/SOC/2010/quasi_random/libs/random/example/Cauchy-Crofton convergence.png   (contents, props changed)
   sandbox/SOC/2010/quasi_random/libs/random/example/box_intersection.hpp   (contents, props changed)
   sandbox/SOC/2010/quasi_random/libs/random/example/cauchy-crofton.cpp   (contents, props changed)
   sandbox/SOC/2010/quasi_random/libs/random/example/detail/
   sandbox/SOC/2010/quasi_random/libs/random/example/detail/constant_size_array.hpp   (contents, props changed)
   sandbox/SOC/2010/quasi_random/libs/random/example/detail/list_assignment.hpp   (contents, props changed)
   sandbox/SOC/2010/quasi_random/libs/random/example/vector_n.hpp   (contents, props changed)
Text files modified: 
   sandbox/SOC/2010/quasi_random/boost/random/faure.hpp |     5 ++---                                   
   1 files changed, 2 insertions(+), 3 deletions(-)

Modified: sandbox/SOC/2010/quasi_random/boost/random/faure.hpp
==============================================================================
--- sandbox/SOC/2010/quasi_random/boost/random/faure.hpp	(original)
+++ sandbox/SOC/2010/quasi_random/boost/random/faure.hpp	2010-08-01 12:16:42 EDT (Sun, 01 Aug 2010)
@@ -305,7 +305,7 @@
   {
     this->curr_elem = 0;
     this->seq_count = init;
-    this->lattice.update(this->seq_count, this->quasi_state);
+    this->lattice.update(init, this->quasi_state);
   }
 
   //=========================Doxygen needs this!==============================
@@ -350,8 +350,7 @@
 /** @cond hide_private_members */
   void compute_next_vector()
   {
-    this->seq_count++;
-    this->lattice.update(this->seq_count, this->quasi_state);
+    this->lattice.update(this->seq_count++, this->quasi_state);
   }
 /** @endcond */
 };

Added: sandbox/SOC/2010/quasi_random/libs/random/example/Cauchy-Crofton convergence.png
==============================================================================
Binary file. No diff available.

Added: sandbox/SOC/2010/quasi_random/libs/random/example/box_intersection.hpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/quasi_random/libs/random/example/box_intersection.hpp	2010-08-01 12:16:42 EDT (Sun, 01 Aug 2010)
@@ -0,0 +1,91 @@
+// Copyright Justinas Vygintas Daugmaudis, 2010.
+// Use, modification and distribution is subject to the
+// Boost Software License, Version 1.0. (See accompanying
+// file LICENSE-1.0 or http://www.boost.org/LICENSE-1.0)
+
+#ifndef HEADER_BOX_INTERSECTION_HPP_INCLUDED
+#define HEADER_BOX_INTERSECTION_HPP_INCLUDED
+
+#include "vector_n.hpp"
+
+namespace example {
+
+// A simple ray structure
+struct ray
+{
+  typedef math::vector_n<float, 3> vector_t;
+
+  explicit ray(const vector_t& o, const vector_t& d)
+    : origin(o)
+    , direction(d)
+    , inv_direction(1.f / d.x(), 1.f / d.y(), 1.f / d.z())
+  {
+    for(int i = 0; i < 3; ++i)
+      sign[i] = (inv_direction(i) < 0.f);
+  }
+
+  vector_t origin;
+  vector_t direction;
+  vector_t inv_direction;
+  int sign[3];
+};
+
+
+// Axis aligned box
+class box
+{
+public:
+  typedef math::vector_n<float, 3> vector_t;
+
+  explicit box(const vector_t& minc, const vector_t& maxc)
+  {
+    bounds[0] = minc;
+    bounds[1] = maxc;
+  }
+
+  explicit box(float width, float height, float depth)
+  {
+    bounds[0] = 0.f, 0.f, 0.f;
+    bounds[1] = width, height, depth;
+  }
+
+  const vector_t& min_corner() const { return bounds[0]; }
+  const vector_t& max_corner() const { return bounds[1]; }
+
+  bool intersect(const ray& r, float t0, float t1) const
+  {
+    float tmin, tmax, tymin, tymax, tzmin, tzmax;
+    tmin = (bounds[r.sign[0]].x() - r.origin.x()) * r.inv_direction.x();
+    tmax = (bounds[1-r.sign[0]].x() - r.origin.x()) * r.inv_direction.x();
+    tymin = (bounds[r.sign[1]].y() - r.origin.y()) * r.inv_direction.y();
+    tymax = (bounds[1-r.sign[1]].y() - r.origin.y()) * r.inv_direction.y();
+    if ( (tmin > tymax) || (tymin > tmax) )
+      return false;
+    if (tymin > tmin)
+      tmin = tymin;
+    if (tymax < tmax)
+      tmax = tymax;
+    tzmin = (bounds[r.sign[2]].z() - r.origin.z()) * r.inv_direction.z();
+    tzmax = (bounds[1-r.sign[2]].z() - r.origin.z()) * r.inv_direction.z();
+    if ( (tmin > tzmax) || (tzmin > tmax) )
+      return false;
+    if (tzmin > tmin)
+      tmin = tzmin;
+    if (tzmax < tmax)
+      tmax = tzmax;
+    return ( (tmin < t1) && (tmax > t0) );
+  }
+
+  float surface_area() const
+  {
+    vector_t diagonal = bounds[1] - bounds[0];
+    return 2.f * (diagonal.x()*diagonal.y() + diagonal.y()*diagonal.z() + diagonal.z()*diagonal.x());
+  }
+
+private:
+  vector_t bounds[2];
+};
+
+} // namespace example
+
+#endif // HEADER_BOX_INTERSECTION_HPP_INCLUDED

Added: sandbox/SOC/2010/quasi_random/libs/random/example/cauchy-crofton.cpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/quasi_random/libs/random/example/cauchy-crofton.cpp	2010-08-01 12:16:42 EDT (Sun, 01 Aug 2010)
@@ -0,0 +1,112 @@
+// Copyright Justinas Vygintas Daugmaudis, 2010.
+// Use, modification and distribution is subject to the
+// Boost Software License, Version 1.0. (See accompanying
+// file LICENSE-1.0 or http://www.boost.org/LICENSE-1.0)
+
+#include <boost/random/niederreiter_base2.hpp>
+#include <boost/random/uniform_real.hpp>
+#include <boost/math/constants/constants.hpp>
+
+#include <iostream>
+
+//
+// DESCRIPTION:
+// ~~~~~~~~~~~~
+//
+// Samples the surface area of a box using the Cauchy-Crofton method.
+// The example is intentionally simple. It is important to note that while
+// the surface area of a box is easy to compute, it is not always the case,
+// especially if one needs to compute the surface area of an isosurface. This,
+// for example, is useful when one computes solvent excluded protein surface area;
+// in that case one would replace a 'box' with, for example, a 'protein_model'.
+
+#include "box_intersection.hpp"
+
+
+template<typename Engine, typename Vector, typename T>
+void uniform_chords(Engine& gen, const Vector& center, T radius, std::size_t n_chords,
+    std::vector<Vector>& points)
+{
+  const T angle = 2 * boost::math::constants::pi<T>();
+
+  const std::size_t n_points = n_chords * 2;
+
+  boost::uniform_real<T> rnd;
+
+  Vector point_on_sphere;
+
+  // Using formula from Rovira 05 to compute points uniformly
+  // distributed on the surface of a sphere
+  for(std::size_t i = 0; i != n_points; ++i)
+  {
+    T cos_theta = 1 - 2 * rnd(gen);
+    T sin_theta = std::sqrt(1 - (cos_theta * cos_theta));
+    T phi = angle * rnd(gen);
+    T sin_phi = std::sin(phi), cos_phi = std::cos(phi); // consider using sincos function
+
+    point_on_sphere = sin_theta*sin_phi, cos_theta, sin_theta*cos_phi;
+
+    point_on_sphere *= radius;
+    point_on_sphere += center;
+
+    points.push_back(point_on_sphere);
+  } // surface point generation
+}
+
+
+template<typename Engine>
+void cauchy_crofton(Engine& gen, const example::box& model, std::size_t n_chords)
+{
+  typedef typename example::box::vector_t vector_t;
+  typedef typename vector_t::value_type value_t;
+
+  std::cout << "Sample surface area, number of chords = " << n_chords << "\n";
+
+  vector_t bounding_sphere_center = (model.min_corner() + model.max_corner()) * 0.5f;
+  value_t bounding_sphere_radius = math::distance(bounding_sphere_center, model.max_corner());
+  value_t bounding_sphere_area = 4 * boost::math::constants::pi<value_t>() *
+      (bounding_sphere_radius * bounding_sphere_radius);
+
+  std::vector<vector_t> points;
+
+  // Generate uniformly distributed chord entry and exit points on a S^2 ball
+  uniform_chords(gen, bounding_sphere_center, bounding_sphere_radius, n_chords, points);
+
+  // Compute the number of chords that intersect a given Model
+  std::size_t points_on_model = 0;
+  for(std::size_t i = 0; i != points.size(); i += 2)
+  {
+    example::ray r(points[i], points[i+1] - points[i]);
+    if( model.intersect(r, 0, 1) ) {
+      points_on_model += 2; // this is an approximation,
+    }                       // but sufficient enough for this particular example
+  }
+
+  value_t sampled_area =
+      static_cast<value_t>(points_on_model) /
+        static_cast<value_t>(points.size())
+    * bounding_sphere_area;
+
+  std::cout
+    << "Bounding sphere area: " << bounding_sphere_area << "\n"
+    << "Points on model: " << points_on_model << "\n"
+    << "Points on sphere: " << points.size() << "\n"
+    << "Sampled area: " << sampled_area << "\n"
+    ;
+}
+
+
+int main()
+{
+  example::box aab(1.f, 0.3f, 0.2f); // a simple slightly oblong box
+
+  std::cout << "The analytical surface area of a box: " << aab.surface_area() << "\n";
+
+  boost::niederreiter_base2_4d gen_4d; // QRNG must be 4-dimensional!
+  cauchy_crofton(gen_4d, aab, 500);
+
+  gen_4d.seed(); // restart the generator in order to repeat the experiment above
+  cauchy_crofton(gen_4d, aab, 5000);
+
+  return 0;
+}

Added: sandbox/SOC/2010/quasi_random/libs/random/example/detail/constant_size_array.hpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/quasi_random/libs/random/example/detail/constant_size_array.hpp	2010-08-01 12:16:42 EDT (Sun, 01 Aug 2010)
@@ -0,0 +1,37 @@
+// Copyright Justinas Vygintas Daugmaudis, 2010.
+// Use, modification and distribution is subject to the
+// Boost Software License, Version 1.0. (See accompanying
+// file LICENSE-1.0 or http://www.boost.org/LICENSE-1.0)
+
+#ifndef HEADER_DETAIL_CONSTANT_SIZE_ARRAY_HPP_INCLUDED
+#define HEADER_DETAIL_CONSTANT_SIZE_ARRAY_HPP_INCLUDED
+
+#include <boost/array.hpp>
+
+namespace detail {
+
+// This is a minimal wrapper around boost::array so that it can be
+// used as an underlaying container for uBLAS data structures
+// such as vector and matrix
+template<class T, std::size_t N>
+struct constant_size_array : public boost::array<T, N>
+{
+    constant_size_array(std::size_t m)
+    {
+        // Assert that this array holds exactly the same
+        // size of memory block as requested
+        BOOST_ASSERT(m == N);
+    }
+
+    void resize(std::size_t) {
+        BOOST_ASSERT("Cannot resize a constant sized array." == 0);
+    }
+
+    void resize(std::size_t size, T init) {
+        BOOST_ASSERT("Cannot resize a constant sized array." == 0);
+    }
+};
+
+} // namespace detail
+
+#endif // HEADER_DETAIL_CONSTANT_SIZE_ARRAY_HPP_INCLUDED

Added: sandbox/SOC/2010/quasi_random/libs/random/example/detail/list_assignment.hpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/quasi_random/libs/random/example/detail/list_assignment.hpp	2010-08-01 12:16:42 EDT (Sun, 01 Aug 2010)
@@ -0,0 +1,112 @@
+// Copyright Justinas Vygintas Daugmaudis, 2010.
+// Use, modification and distribution is subject to the
+// Boost Software License, Version 1.0. (See accompanying
+// file LICENSE-1.0 or http://www.boost.org/LICENSE-1.0)
+
+#ifndef HEADER_DETAIL_LIST_ASSIGNMENT_HPP_INCLUDED
+#define HEADER_DETAIL_LIST_ASSIGNMENT_HPP_INCLUDED
+
+namespace detail {
+
+// Forward declaration
+template<class T, class ForwardIterator>
+class scalar_or_list_assignment;
+
+
+template<class T, typename ForwardIterator>
+class list_assignment
+{
+private:
+  typedef list_assignment<T, ForwardIterator>
+      self_type;
+  typedef scalar_or_list_assignment<T, ForwardIterator>
+      switch_type;
+
+  ForwardIterator m_iterator_;
+
+public:
+  explicit list_assignment(switch_type& rhs);
+  explicit list_assignment(ForwardIterator pos)
+    : m_iterator_(pos)
+  {}
+
+  inline self_type operator,(T value)
+  {
+      *m_iterator_ = value;
+      return self_type(++m_iterator_);
+  }
+};
+
+
+template<class T, class ForwardIterator>
+class scalar_or_list_assignment
+{
+private:
+  friend class list_assignment<T, ForwardIterator>;
+
+  typedef scalar_or_list_assignment<T, ForwardIterator>
+      self_t;
+  typedef list_assignment<T, ForwardIterator>
+      assign_from_list;
+
+public:
+  explicit scalar_or_list_assignment(ForwardIterator first, ForwardIterator last, T value)
+      : m_first(first)
+      , m_last(last)
+      , m_value(value)
+  {
+      // Immediately initialize the first element
+      // and move iterator position forward.
+      *m_first = value; ++m_first;
+  }
+
+  scalar_or_list_assignment(const self_t& rhs)
+      : m_first(rhs.m_first)
+      , m_last(rhs.m_last)
+      , m_value(rhs.m_value)
+  {
+      // This is, admitedly, not very pretty, but quite necessary
+      // for optimization purposes; initialization would be O(n^2)
+      // instead of O(n) otherwise.
+      const_cast<self_t *>(&rhs)->disable();
+  }
+
+  ~scalar_or_list_assignment()
+  {
+      std::fill(m_first, m_last, m_value);
+  }
+
+  // Continues initialization of [first + 1, last)
+  // by consecutive assignment of given scalar values
+  // to the underlying data structure.
+  assign_from_list operator,(T value)
+  {
+      *m_first = value; ++m_first;
+      return assign_from_list(*this);
+  }
+
+private:
+  void disable() {
+      m_first = m_last;
+  }
+
+private:
+  ForwardIterator m_first;
+  ForwardIterator m_last;
+  T m_value;
+};
+
+
+//=============================================================================
+// list_assignment c-tor implementation
+
+template<class T, typename Iterator>
+inline list_assignment<T, Iterator>::list_assignment(scalar_or_list_assignment<T, Iterator>& rhs)
+   : m_iterator_(rhs.m_first)
+{
+  rhs.disable();
+}
+
+} // namespace detail
+
+#endif // HEADER_DETAIL_LIST_ASSIGNMENT_HPP_INCLUDED

Added: sandbox/SOC/2010/quasi_random/libs/random/example/vector_n.hpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2010/quasi_random/libs/random/example/vector_n.hpp	2010-08-01 12:16:42 EDT (Sun, 01 Aug 2010)
@@ -0,0 +1,112 @@
+// Copyright Justinas Vygintas Daugmaudis, 2010.
+// Use, modification and distribution is subject to the
+// Boost Software License, Version 1.0. (See accompanying
+// file LICENSE-1.0 or http://www.boost.org/LICENSE-1.0)
+
+#ifndef HEADER_VECTOR_N_HPP_INCLUDED
+#define HEADER_VECTOR_N_HPP_INCLUDED
+
+#include <boost/call_traits.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+#include "detail/list_assignment.hpp"
+#include "detail/constant_size_array.hpp"
+
+namespace math {
+
+template<class T, std::size_t N>
+class vector_n : public boost::numeric::ublas::vector<
+    T
+  , detail::constant_size_array<T, N> >
+{
+private:
+  typedef detail::constant_size_array<T, N>
+      array_type;
+  typedef typename boost::call_traits<T>::param_type
+      param_type;
+
+  typedef boost::numeric::ublas::vector<
+      T
+    , array_type > base_vector;
+
+  typedef detail::scalar_or_list_assignment<
+      param_type
+    , typename array_type::iterator > assignment_expression;
+public:
+  //! Default construction.
+  vector_n(): base_vector(N) {}
+
+  //! Convenience constructors
+  vector_n(T v0) : base_vector(N)
+  {
+    BOOST_STATIC_ASSERT( N == 1 );
+    this->data()[0] = v0;
+  }
+
+  vector_n(T v0, T v1) : base_vector(N)
+  {
+    BOOST_STATIC_ASSERT( N == 2 );
+    this->data()[0] = v0;
+    this->data()[1] = v1;
+  }
+
+  vector_n(T v0, T v1, T v2) : base_vector(N)
+  {
+    BOOST_STATIC_ASSERT( N == 3 );
+    this->data()[0] = v0;
+    this->data()[1] = v1;
+    this->data()[2] = v2;
+  }
+
+  vector_n(T v0, T v1, T v2, T v3) : base_vector(N)
+  {
+    BOOST_STATIC_ASSERT( N == 4 );
+    this->data()[0] = v0;
+    this->data()[1] = v1;
+    this->data()[2] = v2;
+    this->data()[3] = v3;
+  }
+
+  //! Construction and assignment
+  //! from a uBLAS vector expression
+  //! or copy assignment.
+  template<class R>
+  vector_n(const boost::numeric::ublas::vector_expression<R>& r)
+    : base_vector(r)
+  {}
+
+  template<class R>
+  vector_n& operator=(const boost::numeric::ublas::vector_expression<R>& r)
+  {
+    base_vector::operator=(r);
+    return *this;
+  }
+
+  // Assignment from a scalar operand or a list of scalar operands.
+  assignment_expression operator=(param_type value)
+  {
+    return assignment_expression(this->data().begin(), this->data().end(), value);
+  }
+
+  // Convenient getters.
+  T x() const { return this->data()[0]; }
+  T y() const { BOOST_STATIC_ASSERT(N > 1); return this->data()[1]; }
+  T z() const { BOOST_STATIC_ASSERT(N > 2); return this->data()[2]; }
+  T w() const { BOOST_STATIC_ASSERT(N > 3); return this->data()[3]; }
+};
+
+template<class T, std::size_t N>
+inline T length_sq(const vector_n<T, N>& v)
+{
+  return boost::numeric::ublas::inner_prod(v, v);
+}
+
+template<typename T, std::size_t N>
+inline T distance(vector_n<T, N> lhs, const vector_n<T, N>& rhs)
+{
+  lhs -= rhs;
+  return std::sqrt(length_sq(lhs));
+}
+
+} // namespace math
+
+#endif // HEADER_VECTOR_N_HPP_INCLUDED