Re: Add uniform_inside_sphere_distribution

[Ed Smith-Rowland <3dw4rd at verizon dot net>]

On 11/06/2014 02:23 AM, Marc Glisse wrote:
> On Wed, 5 Nov 2014, Ed Smith-Rowland wrote:
>
> > On 11/05/2014 04:25 PM, Marc Glisse wrote:
> > > On Wed, 5 Nov 2014, Ed Smith-Rowland wrote:
> > >
> > > > Like the uniform_on_sphere_distribution which is used inside, the 
> > > > 2-dimensional case uses rejection
> > >
> > > Could you point out where in the code you are special-casing 
> > > dimension 2? Somehow I can't see it in the patch.
> > My miscommunication.  Sorry.
> >
> > The special casing - and all the hard work - is done by the contained 
> > uniform_on_sphere_distribution.
>
> In dimension 2, uniform_on_sphere_distribution generates a uniform 
> random point in the unit disk and projects it to the circle. Using 
> that + one more random number to regenerate a uniform random point in 
> the disk seems wasteful to me.
OK.

That was a sobering exercise in the difference between asymptotic 
analysis and reality.  I measured transform vs. rejection from Dim=1 to 
9 on a couple boxen.
Time for rejection is a power law like t = A.2^D and transform was 
linear t = B + C.D.  On the machine with results below A ~ 12, B ~ 20, C 
~ 880.  Rather different numbers but same results on an older box.  
Basically, the linear coefficient is huge.  I can't blame pow or sqrt 
though by the numbers.  My timer is crude though.  i'll try more tools.

ed@bad-horse:~/tr2$ LD_LIBRARY_PATH=/home/ed/bin/lib64:$LD_LIBRARY_PATH 
./deathmatch 100000

  New (Rejection) Distribution - Dim=1
  Total time per iter = 38 ms
  Mean time per iter  = 0.00038 ms

  Old (Transform) Distribution - Dim=1
  Total time per iter = 905 ms
  Mean time per iter  = 0.00905 ms

  New (Rejection) Distribution - Dim=2
  Total time per iter = 48 ms
  Mean time per iter  = 0.00048 ms

  Old (Transform) Distribution - Dim=2
  Total time per iter = 1776 ms
  Mean time per iter  = 0.01776 ms

  New (Rejection) Distribution - Dim=3
  Total time per iter = 104 ms
  Mean time per iter  = 0.00104 ms

  Old (Transform) Distribution - Dim=3
  Total time per iter = 2659 ms
  Mean time per iter  = 0.02659 ms

  New (Rejection) Distribution - Dim=4
  Total time per iter = 227 ms
  Mean time per iter  = 0.00227 ms

  Old (Transform) Distribution - Dim=4
  Total time per iter = 3532 ms
  Mean time per iter  = 0.03532 ms

  New (Rejection) Distribution - Dim=5
  Total time per iter = 525 ms
  Mean time per iter  = 0.00525 ms

  Old (Transform) Distribution - Dim=5
  Total time per iter = 4402 ms
  Mean time per iter  = 0.04402 ms

  New (Rejection) Distribution - Dim=6
  Total time per iter = 1275 ms
  Mean time per iter  = 0.01275 ms

  Old (Transform) Distribution - Dim=6
  Total time per iter = 5299 ms
  Mean time per iter  = 0.05299 ms

  New (Rejection) Distribution - Dim=7
  Total time per iter = 3269 ms
  Mean time per iter  = 0.03269 ms

  Old (Transform) Distribution - Dim=7
  Total time per iter = 6193 ms
  Mean time per iter  = 0.06193 ms

<---  BREAK EVEN  --->

  New (Rejection) Distribution - Dim=8
  Total time per iter = 8722 ms
  Mean time per iter  = 0.08722 ms

  Old (Transform) Distribution - Dim=8
  Total time per iter = 7047 ms
  Mean time per iter  = 0.07047 ms

  New (Rejection) Distribution - Dim=9
  Total time per iter = 23814 ms
  Mean time per iter  = 0.23814 ms

  Old (Transform) Distribution - Dim=9
  Total time per iter = 7893 ms
  Mean time per iter  = 0.07893 ms

  Bunch of calls to pow and sqrt
  Total time per iter = 6 ms
  Mean time per iter  = 6e-05 ms

  uniform_on_sphere_distribution - Dim=3
  Total time per iter = 87 ms
  Mean time per iter  = 0.00087 ms

  default_random_engine - Dim=3
  Total time per iter = 4 ms
  Mean time per iter  = 4e-05 ms

The new patch dispatches to rejection for Dim<8, transform otherwise.

I include deathmatch.  You can hack bits from the patch to make the old 
and new headers.

I have to admit I'm puzzled as to why the transform method is so slow 
relative to rejection.

Ed

Attachment: CL_uisd
Description: Text document

Attachment: patch_uisd_2
Description: Text document

// /home/ed/bin/bin/g++ -std=c++1z -o deathmatch deathmatch.cpp

// LD_LIBRARY_PATH=/home/ed/bin/lib64:$LD_LIBRARY_PATH ./deathmatch 100000

#include "uniform_inside_sphere_distribution.h"
#include "uniform_inside_sphere_distribution_old.h"
#include <ext/random>
#include "timer.h"

#include <iostream>

template <std::size_t _Dim>
  void
  test(double radius = 1.0, int num = 1000000)
  {
    Timer timer;

    std::default_random_engine re; // the default engine

    __gnu_cxx::uniform_inside_sphere_distribution<_Dim> insphnew(radius);

    timer.start();
    for (int i = 0; i < num; ++i)
      insphnew(re);
    timer.stop();

    std::cout << std::endl << "  New (Rejection) Distribution - Dim=" << _Dim;
    std::cout << std::endl << "  Total time per iter = "
              << 1.0 * timer.time_elapsed() << " ms";
    std::cout << std::endl << "  Mean time per iter  = "
              << 1.0 * timer.time_elapsed() / num << " ms" << std::endl;

    __old::uniform_inside_sphere_distribution<_Dim> insphold(radius);

    timer.start();
    for (int i = 0; i < num; ++i)
      insphold(re);
    timer.stop();

    std::cout << std::endl << "  Old (Transform) Distribution - Dim=" << _Dim;
    std::cout << std::endl << "  Total time per iter = "
              << 1.0 * timer.time_elapsed() << " ms";
    std::cout << std::endl << "  Mean time per iter  = "
              << 1.0 * timer.time_elapsed() / num << " ms" << std::endl;
  }

int
main(int n_app_args, char ** app_arg)
{
  int num = 10000;
  if (n_app_args > 1)
    num = atoi(app_arg[1]);

  //  Segment

  test<1>(1.0, num);

  //  Circle

  test<2>(2.0, num);

  //  Sphere

  test<3>(3.0, num);

  //  Hypersphere

  test<4>(4.0, num);

  //  Etc...

  test<5>(5.0, num);

  test<6>(6.0, num);

  test<7>(7.0, num);

  test<8>(8.0, num);

  test<9>(9.0, num);

  Timer timer;

  //

  timer.start();
  double sum = 0.0;
  for (int i = 0; i < num; ++i)
  {
    sum += std::sqrt(i*1.0);
    sum -= std::pow(i*1.0, 0.333333);
  }
  timer.stop();
  std::cout << std::endl << "  Bunch of calls to pow";
  std::cout << std::endl << "  Total time per iter = "
            << 1.0 * timer.time_elapsed() << " ms";
  std::cout << std::endl << "  Mean time per iter  = "
            << 1.0 * timer.time_elapsed() / num << " ms" << std::endl;

  //

  std::default_random_engine re; // the default engine

  __gnu_cxx::uniform_on_sphere_distribution<3> onsphnew;

  timer.start();
  for (int i = 0; i < num; ++i)
    onsphnew(re);
  timer.stop();

  std::cout << std::endl << "  uniform_on_sphere_distribution - Dim=" << 3;
  std::cout << std::endl << "  Total time per iter = "
            << 1.0 * timer.time_elapsed() << " ms";
  std::cout << std::endl << "  Mean time per iter  = "
            << 1.0 * timer.time_elapsed() / num << " ms" << std::endl;

  //

  timer.start();
  for (int i = 0; i < num * 3; ++i)
    re();
  timer.stop();

  std::cout << std::endl << "  default_random_engine - Dim=" << 3;
  std::cout << std::endl << "  Total time per iter = "
            << 1.0 * timer.time_elapsed() << " ms";
  std::cout << std::endl << "  Mean time per iter  = "
            << 1.0 * timer.time_elapsed() / num << " ms" << std::endl;
}