C++ std::mt19937类代码示例

void SeedFeatureFactory::train( const std::vector< std::shared_ptr<ImageOverSegmentation> > &ios, const std::vector<VectorXs> & lbl ) {
	printf("  * Training SeedFeature\n");
	static std::mt19937 rand;
	const int N_SAMPLES = 5000;
	int n_pos=0, n_neg=0;
	for( VectorXs l: lbl ) {
		n_pos += (l.array()>=0).cast<int>().sum();
		n_neg += (l.array()==-1).cast<int>().sum();
	}
	
	// Collect training examples
	float sampling_freq[] = {0.5f*N_SAMPLES / n_neg, 0.5f*N_SAMPLES / n_pos};
	std::vector<RowVectorXf> f;
	std::vector<float> l;
#pragma omp parallel for
	for( int i=0; i<ios.size(); i++ ) {
		RMatrixXf ftr = SeedFeature::computeObjFeatures( *ios[i] );
		for( int j=0; j<ios[i]->Ns(); j++ )
			if( lbl[i][j] >= -1 && rand() < rand.max()*sampling_freq[ lbl[i][j]>=0 ] ) {
#pragma omp critical
				{
					l.push_back( lbl[i][j]>=0 );
					f.push_back( ftr.row(j) );
				}
			}
	}
	
	printf("    - Computing parameters\n");
	// Fit the ranking functions
	RMatrixXf A( f.size(), f[0].size() );
	VectorXf b( l.size() );
	for( int i=0; i<f.size(); i++ ) {
		A.row(i) = f[i];
		b[i] = l[i];
	}
	
	// Solve A*x = b
	param_ = A.colPivHouseholderQr().solve(b);
	printf("    - done %f\n",(A*param_-b).array().abs().mean());
}

namespace rng
{
	std::mt19937 r;

	void init() {
		r.seed((unsigned)std::chrono::system_clock::now().time_since_epoch().count());
	}
	
	int d6(std::mt19937& urng) {
		return d6Gen(urng);
	}
	
	int degrees(std::mt19937& urng) {
		return degreesGen(urng);
	}
	
	double radians(std::mt19937& urng) {
		return radiansGen(urng);
	}
	
	bool boolean(std::mt19937& urng) {
		return booleanGen(urng);
	}

	int getInt(int lo, int hi, std::mt19937& urng) {
		return std::uniform_int_distribution<int>(lo, hi)(urng);
	}

	int getInt(int hi, std::mt19937& urng) {
		return std::uniform_int_distribution<int>(0, hi)(urng);
	}

	bool oneInX(double x, std::mt19937& urng) {
		return std::bernoulli_distribution(1.0 / x)(urng);
	}

	bool xInY(double x, double y, std::mt19937& urng) {
		return std::bernoulli_distribution(x / y)(urng);
	}
}

namespace math /** Miscellaneous math routines. */ {

// Global random object.
extern std::mt19937 randGen;
// Global uniform distribution.
extern std::uniform_real_distribution<> randUniformDist;
// Global normal distribution.
extern std::normal_distribution<> randNormalDist;

/**
 * Set the random seed used by the random functions (Random() and RandInt()).
 * The seed is casted to a 32-bit integer before being given to the random
 * number generator, but a size_t is taken as a parameter for API consistency.
 *
 * @param seed Seed for the random number generator.
 */
inline void RandomSeed(const size_t seed)
{
  randGen.seed((uint32_t) seed);
  srand((unsigned int) seed);
#if ARMA_VERSION_MAJOR > 3 || \
    (ARMA_VERSION_MAJOR == 3 && ARMA_VERSION_MINOR >= 930)
  // Armadillo >= 3.930 has its own random number generator internally that we
  // need to set the seed for also.
  arma::arma_rng::set_seed(seed);
#endif
}

/**
 * Generates a uniform random number between 0 and 1.
 */
inline double Random()
{
  return randUniformDist(randGen);
}

/**
 * Generates a uniform random number in the specified range.
 */
inline double Random(const double lo, const double hi)
{
  return lo + (hi - lo) * randUniformDist(randGen);
}

/**
 * Generates a uniform random integer.
 */
inline int RandInt(const int hiExclusive)
{
  return (int) std::floor((double) hiExclusive * randUniformDist(randGen));
}

/**
 * Generates a uniform random integer.
 */
inline int RandInt(const int lo, const int hiExclusive)
{
  return lo + (int) std::floor((double) (hiExclusive - lo)
                               * randUniformDist(randGen));
}

/**
 * Generates a normally distributed random number with mean 0 and variance 1.
 */
inline double RandNormal()
{
  return randNormalDist(randGen);
}

/**
 * Generates a normally distributed random number with specified mean and
 * variance.
 *
 * @param mean Mean of distribution.
 * @param variance Variance of distribution.
 */
inline double RandNormal(const double mean, const double variance)
{
  return variance * randNormalDist(randGen) + mean;
}

} // namespace math

int GameMain()
{
	s_rng.seed(1729);

	TestConstructors();
	TestMagnAndNorm();
	TestMembers();
	TestNonMembers();
	TestOperators();

	Sleep(1000);
	return 0;
}

void init(dynamo::Simulation& Sim, const double density)
{
  RNG.seed(std::random_device()());
  Sim.ranGenerator.seed(std::random_device()());

  const double elasticity = 1.0;
  const size_t cells = 7;
  const double wallkT = 1.0;

  std::unique_ptr<dynamo::UCell> packptr(new dynamo::CUFCC(std::array<long, 3>{{cells, cells, cells}}, dynamo::Vector(1,1,1), new dynamo::UParticle()));
  packptr->initialise();
  std::vector<dynamo::Vector> latticeSites(packptr->placeObjects(dynamo::Vector(0,0,0)));
  const size_t N = latticeSites.size();
  const double boxL = std::cbrt(N / density);
  Sim.primaryCellSize = dynamo::Vector(boxL, boxL, boxL);

  dynamo::shared_ptr<dynamo::ParticleProperty> D(new dynamo::ParticleProperty(N, dynamo::Property::Units::Length(), "D", 1.0));
  dynamo::shared_ptr<dynamo::ParticleProperty> M(new dynamo::ParticleProperty(N, dynamo::Property::Units::Mass(), "M", 1.0));
  Sim._properties.push(D);
  Sim._properties.push(M);

  Sim.dynamics = dynamo::shared_ptr<dynamo::Dynamics>(new dynamo::DynGravity(&Sim, dynamo::Vector(0, -1, 0)));
  Sim.BCs = dynamo::shared_ptr<dynamo::BoundaryCondition>(new dynamo::BCNone(&Sim));
  Sim.ptrScheduler = dynamo::shared_ptr<dynamo::SNeighbourList>(new dynamo::SNeighbourList(&Sim, new dynamo::FELBoundedPQ<dynamo::PELMinMax<3> >()));

  Sim.interactions.push_back(dynamo::shared_ptr<dynamo::Interaction>(new dynamo::IHardSphere(&Sim, "D", elasticity, new dynamo::IDPairRangeAll(), "Bulk")));
  Sim.addSpecies(dynamo::shared_ptr<dynamo::Species>(new dynamo::SpPoint(&Sim, new dynamo::IDRangeAll(&Sim), "M", "Bulk", 0)));

  Sim.locals.push_back(dynamo::shared_ptr<dynamo::Local>(new dynamo::LWall(&Sim, elasticity, "D", dynamo::Vector(1, 0, 0), dynamo::Vector(-0.5 * Sim.primaryCellSize[0] - 1, 0, 0), "XwallLow", new dynamo::IDRangeAll(&Sim))));
  Sim.locals.push_back(dynamo::shared_ptr<dynamo::Local>(new dynamo::LWall(&Sim, elasticity, "D", dynamo::Vector(-1, 0, 0), dynamo::Vector(0.5 * Sim.primaryCellSize[0] + 1, 0, 0), "XwallHigh", new dynamo::IDRangeAll(&Sim))));

  Sim.locals.push_back(dynamo::shared_ptr<dynamo::Local>(new dynamo::LWall(&Sim, elasticity, "D", dynamo::Vector(0, 0, 1), dynamo::Vector(0, 0, -0.5 * Sim.primaryCellSize[2] - 1), "ZwallLow", new dynamo::IDRangeAll(&Sim))));
  Sim.locals.push_back(dynamo::shared_ptr<dynamo::Local>(new dynamo::LWall(&Sim, elasticity, "D", dynamo::Vector(0, 0, -1), dynamo::Vector(0, 0, 0.5 * Sim.primaryCellSize[2] + 1), "ZwallHigh", new dynamo::IDRangeAll(&Sim))));

  Sim.locals.push_back(dynamo::shared_ptr<dynamo::Local>(new dynamo::LWall(&Sim, elasticity, "D", dynamo::Vector(0, 1, 0), dynamo::Vector(0, - 0.5 * Sim.primaryCellSize[1] - 1, 0), "GroundPlate", new dynamo::IDRangeAll(&Sim), wallkT)));
  
  for (size_t i(0); i < latticeSites.size(); ++i)
    {
      Sim.particles.push_back(dynamo::Particle(boxL * latticeSites[i], getRandVelVec(), Sim.particles.size()));
      D->getProperty(i) = (i < 100) ? 1 : 0.5;
      M->getProperty(i) = (i < 100) ? 1 : 0.5 * 0.5 * 0.5;
    }

  Sim.ensemble = dynamo::Ensemble::loadEnsemble(Sim);

  dynamo::InputPlugin(&Sim, "Rescaler").zeroMomentum();
  dynamo::InputPlugin(&Sim, "Rescaler").rescaleVels(1.0);

  BOOST_CHECK_EQUAL(Sim.N(), 1372);
  BOOST_CHECK_CLOSE(Sim.getNumberDensity() * Sim.units.unitVolume(), density, 0.000000001);
}

void init(dynamo::Simulation& Sim, const double density)
{
  RNG.seed(std::random_device()());
  Sim.ranGenerator.seed(std::random_device()());

  Sim.dynamics = dynamo::shared_ptr<dynamo::Dynamics>(new dynamo::DynNewtonian(&Sim));
  Sim.BCs = dynamo::shared_ptr<dynamo::BoundaryCondition>(new dynamo::BCPeriodic(&Sim));
  Sim.ptrScheduler = dynamo::shared_ptr<dynamo::SNeighbourList>(new dynamo::SNeighbourList(&Sim, new DefaultSorter()));

  std::unique_ptr<dynamo::UCell> packptr(new dynamo::CUSC(std::array<long, 3>{{128, 128, 1}}, dynamo::Vector{1,1,1}, new dynamo::UParticle()));
  packptr->initialise();
  std::vector<dynamo::Vector> latticeSites(packptr->placeObjects(dynamo::Vector{0,0,0}));
  
  //As we're in 2D we need to take the square root
  double particleDiam = std::sqrt(density / latticeSites.size());
  Sim.units.setUnitLength(particleDiam);
  Sim.units.setUnitTime(particleDiam); 

  Sim.primaryCellSize = dynamo::Vector{1, 1, 4 * particleDiam};
  
  typedef std::pair<double,double> Step;
  std::vector<Step> steps;
  steps.push_back(Step{1.0,0.1});
  steps.push_back(Step{0.9,0.2});
  steps.push_back(Step{0.8,0.3});
  steps.push_back(Step{0.7,0.4});
  steps.push_back(Step{0.6,0.5});
  steps.push_back(Step{0.5,0.6});
  steps.push_back(Step{0.4,0.7});
  steps.push_back(Step{0.3,0.8});
  steps.push_back(Step{0.2,0.9});
  steps.push_back(Step{0.1,1.0});
  Sim.interactions.push_back(dynamo::shared_ptr<dynamo::Interaction>(new dynamo::IStepped(&Sim, dynamo::shared_ptr<dynamo::Potential>(new dynamo::PotentialStepped(steps, false)), new dynamo::IDPairRangeAll(), "Bulk", particleDiam, 1.0)));
    
  Sim.addSpecies(dynamo::shared_ptr<dynamo::Species>(new dynamo::SpPoint(&Sim, new dynamo::IDRangeAll(&Sim), 1.0, "Bulk", 0)));

  unsigned long nParticles = 0;
  Sim.particles.reserve(latticeSites.size());
  for (const dynamo::Vector & position : latticeSites)
    Sim.particles.push_back(dynamo::Particle(position, getRandVelVec() * Sim.units.unitVelocity(), nParticles++));

  for (auto& particle : Sim.particles)
    particle.getVelocity()[2] = 0;

  dynamo::InputPlugin(&Sim, "Rescaler").zeroMomentum();
  dynamo::InputPlugin(&Sim, "Rescaler").rescaleVels(2.0 / 3.0);
    
  Sim.ensemble = dynamo::Ensemble::loadEnsemble(Sim);

  BOOST_CHECK_EQUAL(Sim.N(), 128 * 128);
}

void init(dynamo::Simulation& Sim, const double density)
{
  RNG.seed(std::random_device()());
  Sim.ranGenerator.seed(std::random_device()());

  double massFrac = 0.001, sizeRatio = 0.5;
  size_t Na=100;
  
  Sim.dynamics = dynamo::shared_ptr<dynamo::Dynamics>(new dynamo::DynNewtonian(&Sim));
  Sim.BCs = dynamo::shared_ptr<dynamo::BoundaryCondition>(new dynamo::BCPeriodic(&Sim));
  Sim.ptrScheduler = dynamo::shared_ptr<dynamo::SNeighbourList>(new dynamo::SNeighbourList(&Sim, new DefaultSorter()));

  std::unique_ptr<dynamo::UCell> packptr(new dynamo::CUFCC(std::array<long, 3>{{10, 10, 10}}, dynamo::Vector(1, 1, 1), new dynamo::UParticle()));
  packptr->initialise();
  std::vector<dynamo::Vector> latticeSites(packptr->placeObjects(dynamo::Vector(0,0,0)));
  Sim.primaryCellSize = dynamo::Vector(1,1,1);

  double simVol = 1.0;
  for (size_t iDim = 0; iDim < NDIM; ++iDim)
    simVol *= Sim.primaryCellSize[iDim];

  double particleDiam = std::cbrt(simVol * density / latticeSites.size());
  
  Sim.interactions.push_back(dynamo::shared_ptr<dynamo::Interaction>(new dynamo::IHardSphere(&Sim, particleDiam, new dynamo::IDPairRangeSingle(new dynamo::IDRangeRange(0, Na - 1)), "AAInt")));
  
  Sim.interactions.push_back(dynamo::shared_ptr<dynamo::Interaction>(new dynamo::IHardSphere(&Sim, ((1.0 + sizeRatio) / 2.0) * particleDiam, new dynamo::IDPairRangePair(new dynamo::IDRangeRange(0, Na - 1), new dynamo::IDRangeRange(Na, latticeSites.size() - 1)), "ABInt")));

  Sim.interactions.push_back(dynamo::shared_ptr<dynamo::Interaction>(new dynamo::IHardSphere(&Sim, sizeRatio * particleDiam, new dynamo::IDPairRangeAll(), "BBInt")));

  Sim.addSpecies(dynamo::shared_ptr<dynamo::Species>(new dynamo::SpPoint(&Sim, new dynamo::IDRangeRange(0, Na - 1), 1.0, "A", 0)));

  Sim.addSpecies(dynamo::shared_ptr<dynamo::Species>(new dynamo::SpPoint(&Sim, new dynamo::IDRangeRange(Na, latticeSites.size() - 1), massFrac, "B", 0)));

  Sim.units.setUnitLength(particleDiam);

  unsigned long nParticles = 0;
  Sim.particles.reserve(latticeSites.size());

  for (const dynamo::Vector & position : latticeSites)
    Sim.particles.push_back(dynamo::Particle(position, getRandVelVec() * Sim.units.unitVelocity(), nParticles++));

  Sim.ensemble = dynamo::Ensemble::loadEnsemble(Sim);

  dynamo::InputPlugin(&Sim, "Rescaler").zeroMomentum();
  dynamo::InputPlugin(&Sim, "Rescaler").rescaleVels(1.0);

  BOOST_CHECK_EQUAL(Sim.N(), 4000);
  //BOOST_CHECK_CLOSE(Sim.getNumberDensity() * Sim.units.unitVolume(), density, 0.000000001);
  //BOOST_CHECK_CLOSE(Sim.getPackingFraction(), Sim.getNumberDensity() * Sim.units.unitVolume() * M_PI / 6.0, 0.000000001);
}

int main(int argc, char** argv) {
  rng.seed();
  CALL(TEST_HashFunctions);
  CALL(overview::TEST_Construction);
  CALL(overview::TEST_Retreival);
  CALL(overview::TEST_Resizing);
  CALL(separate_chaining::TEST_Construction);
  CALL(separate_chaining::TEST_Retreival);
  CALL(open_addressing::TEST_Construction);
  CALL(open_addressing::TEST_Retreival);
  CALL(cuckoo::TEST_CuckooHashing);
  CALL(TEST_MostCommon);
  CALL(TEST_Cache);
  return 0;
}

// -------------------------------------------------------------------------
int hpx_main(boost::program_options::variables_map& vm)
{
    std::size_t device     = vm["device"].as<std::size_t>();
    std::size_t sizeMult   = vm["sizemult"].as<std::size_t>();
    std::size_t iterations = vm["iterations"].as<std::size_t>();
    //
    unsigned int seed = std::random_device{}();
     if (vm.count("seed"))
        seed = vm["seed"].as<unsigned int>();

    gen.seed(seed);
    std::cout << "using seed: " << seed << std::endl;

    //
    sizeMult = (std::min)(sizeMult, std::size_t(100));
    sizeMult = (std::max)(sizeMult, std::size_t(1));
    //
    // use a larger block size for Fermi and above, query default cuda target properties
    hpx::compute::cuda::target target(device);

    std::cout
        << "GPU Device " << target.native_handle().get_device()
        << ": \"" << target.native_handle().processor_name() << "\" "
        << "with compute capability " << target.native_handle().processor_family()
        << "\n";

    int block_size = (target.native_handle().processor_family() < 2) ? 16 : 32;

    sMatrixSize matrix_size;
    matrix_size.uiWA = 2 * block_size * sizeMult;
    matrix_size.uiHA = 4 * block_size * sizeMult;
    matrix_size.uiWB = 2 * block_size * sizeMult;
    matrix_size.uiHB = 4 * block_size * sizeMult;
    matrix_size.uiWC = 2 * block_size * sizeMult;
    matrix_size.uiHC = 4 * block_size * sizeMult;

    printf("MatrixA(%u,%u), MatrixB(%u,%u), MatrixC(%u,%u)\n\n",
           matrix_size.uiWA, matrix_size.uiHA,
           matrix_size.uiWB, matrix_size.uiHB,
           matrix_size.uiWC, matrix_size.uiHC);

    matrixMultiply<float>(matrix_size, device, iterations);
    return hpx::finalize();
}

void init(dynamo::Simulation& Sim, double density = 0.5)
{
  RNG.seed(std::random_device()());
  Sim.ranGenerator.seed(std::random_device()());

  const double elasticity = 1.0;
  const double lambda = 1.5;
  const double welldepth = 1.0;

  Sim.dynamics = dynamo::shared_ptr<dynamo::Dynamics>(new dynamo::DynNewtonian(&Sim));
  Sim.BCs = dynamo::shared_ptr<dynamo::BoundaryCondition>(new dynamo::BCPeriodic(&Sim));
  Sim.ptrScheduler = dynamo::shared_ptr<dynamo::SNeighbourList>(new dynamo::SNeighbourList(&Sim, new DefaultSorter()));

  std::unique_ptr<dynamo::UCell> packptr(new dynamo::CUFCC(std::array<long, 3>{{7,7,7}}, dynamo::Vector(1,1,1), new dynamo::UParticle()));
  packptr->initialise();
  std::vector<dynamo::Vector> latticeSites(packptr->placeObjects(dynamo::Vector(0,0,0)));
  Sim.primaryCellSize = dynamo::Vector(1,1,1);

  double simVol = 1.0;
  for (size_t iDim = 0; iDim < NDIM; ++iDim)
    simVol *= Sim.primaryCellSize[iDim];

  double particleDiam = std::cbrt(simVol * density / latticeSites.size());

  Sim.interactions.push_back(dynamo::shared_ptr<dynamo::Interaction>(new dynamo::ISquareWell(&Sim, particleDiam, lambda, welldepth, elasticity, new dynamo::IDPairRangeAll(), "Bulk")));
  Sim.addSpecies(dynamo::shared_ptr<dynamo::Species>(new dynamo::SpPoint(&Sim, new dynamo::IDRangeAll(&Sim), 1.0, "Bulk", 0)));
  Sim.units.setUnitLength(particleDiam);
  Sim.units.setUnitTime(particleDiam); 

  unsigned long nParticles = 0;
  Sim.particles.reserve(latticeSites.size());
  for (const dynamo::Vector & position : latticeSites)
    Sim.particles.push_back(dynamo::Particle(position, getRandVelVec() * Sim.units.unitVelocity(), nParticles++));

  Sim.ensemble = dynamo::Ensemble::loadEnsemble(Sim);

  dynamo::InputPlugin(&Sim, "Rescaler").zeroMomentum();
  dynamo::InputPlugin(&Sim, "Rescaler").rescaleVels(1.0);

  BOOST_CHECK_EQUAL(Sim.N(), 1372);
  BOOST_CHECK_CLOSE(Sim.getNumberDensity() * Sim.units.unitVolume(), density, 0.000000001);
  BOOST_CHECK_CLOSE(Sim.getPackingFraction(), Sim.getNumberDensity() * Sim.units.unitVolume() * M_PI / 6.0, 0.000000001);
}

Random::Random(UInt64 seed) {
  if (seed == 0) {
    if( !static_gen_seeded ) {
      #if NDEBUG
        unsigned int static_seed = (unsigned int)std::chrono::system_clock::now().time_since_epoch().count();
      #else
        unsigned int static_seed = DEBUG_RANDOM_SEED;
      #endif
      static_gen.seed( static_seed );
      static_gen_seeded = true;
      NTA_INFO << "Random seed: " << static_seed;
    }
    seed_ = static_gen(); //generate random value from HW RNG
  } else {
    seed_ = seed;
  }
  // if seed is zero at this point, there is a logic error.
  NTA_CHECK(seed_ != 0);
  gen.seed((unsigned int)seed_); //seed the generator
  steps_ = 0;
}

    void use_graph_generator(const vle::devs::InitEventList& events)
    {
        model_number = events.getInt("model-number");
        std::string generator_name = events.getString("graph-generator");

        if (events.exist("graph-seed"))
            generator.seed(events.getInt("graph-seed"));

        if (generator_name == "defined")
            use_defined_graph_generator(events);
        else if (generator_name == "small-world")
            use_smallworld_graph_generator(events);
        else if (generator_name == "scale-free")
            use_scalefree_graph_generator(events);
        else if (generator_name == "sorted-erdos-renyi")
            use_sortederdesrenyi_graph_generator(events);
        else if (generator_name == "erdos_renyi")
            use_erdosrenyi_graph_generator(events);
        else
            throw vle::utils::ModellingError("Unknown graph gererator: %s",
                                             generator_name.c_str());
    }

void init(dynamo::Simulation& Sim, const double density)
{
  RNG.seed(std::random_device()());
  Sim.ranGenerator.seed(std::random_device()());

  const double elasticity = 1.0;
  const size_t N = 1000;
  Sim.dynamics = dynamo::shared_ptr<dynamo::Dynamics>(new dynamo::DynNewtonian(&Sim));
  Sim.BCs = dynamo::shared_ptr<dynamo::BoundaryCondition>(new dynamo::BCPeriodic(&Sim));
  Sim.ptrScheduler = dynamo::shared_ptr<dynamo::SNeighbourList>(new dynamo::SNeighbourList(&Sim, new DefaultSorter()));
  
  dynamo::CURandom packroutine(N, dynamo::Vector{1,1,1}, new dynamo::UParticle());
  packroutine.initialise();
  std::vector<dynamo::Vector> latticeSites(packroutine.placeObjects(dynamo::Vector{0,0,0}));
  Sim.BCs = dynamo::shared_ptr<dynamo::BoundaryCondition>(new dynamo::BCPeriodic(&Sim));
  double particleDiam = std::cbrt(density / N);
  Sim.interactions.push_back(dynamo::shared_ptr<dynamo::Interaction>(new dynamo::ILines(&Sim, particleDiam, elasticity, new dynamo::IDPairRangeAll(), "Bulk")));
  Sim.addSpecies(dynamo::shared_ptr<dynamo::Species>(new dynamo::SpSphericalTop(&Sim, new dynamo::IDRangeAll(&Sim), 1.0, "Bulk", 0, particleDiam * particleDiam / 12.0)));
  Sim.units.setUnitLength(particleDiam);

  unsigned long nParticles = 0;
  Sim.particles.reserve(latticeSites.size());
  for (const dynamo::Vector & position : latticeSites)
    Sim.particles.push_back(dynamo::Particle(position, getRandVelVec() * Sim.units.unitVelocity(), nParticles++));

  Sim.dynamics->initOrientations();

  Sim.ensemble = dynamo::Ensemble::loadEnsemble(Sim);

  dynamo::InputPlugin(&Sim, "Rescaler").zeroMomentum();
  dynamo::InputPlugin(&Sim, "Rescaler").rescaleVels(1.0);

  BOOST_CHECK_EQUAL(Sim.N(), N);
  BOOST_CHECK_CLOSE(Sim.getNumberDensity() * Sim.units.unitVolume(), density, 0.000000001);
  BOOST_CHECK_CLOSE(Sim.getPackingFraction(), 0, 0.000000001);
}

/*******************************************************************************
Function to initialize ODE.
*******************************************************************************/
void initODE()
{
	///////////////// Initializing the ODE general features ////////////////////

	dInitODE();								//Initialize library.
	World = dWorldCreate();					//Crate a new dynamics, empty world.
	Space = dSimpleSpaceCreate(0);			//Create a new space for collision (independent).
	ContactGroup = dJointGroupCreate(0);	//Create a joints container, without specifying size.

	dWorldSetGravity( World, 0.0, -9.81, 0 );	//Add gravity to this World.

	//Define error conrrection constants.
	dWorldSetERP( World, 0.2 );
	dWorldSetCFM( World, 1e-5 );

	//Set the velocity that interpenetrating objects will separate at.
	dWorldSetContactMaxCorrectingVel( World, 0.9 );

	//Set the safety area for contacts to be considered at rest.
	dWorldSetContactSurfaceLayer( World, 0.001 );

	//Automatically disable objects that have come to a rest.
	dWorldSetAutoDisableFlag( World, false );

	/////////////// Initializing the rigid bodies in the world /////////////////

	//Create a collision plane and add it to space. The next parameters are the
	//literal components of ax + by + cz = d.
	dCreatePlane( Space, 0.0, 1.0, 0.0, 0.0 );

	const dReal xPos = 0;
	const dReal yPos = 5;
	const dReal zPos = 0;
	const dReal xRot = 0;
	const dReal yRot = 0;
	const dReal zRot = 0;
	const dReal radius = .75;
	const dReal length = 4;
	const dReal sides[3] = { 2, 2, 2 };
	

	//Create body
	body.Body = dBodyCreate(World);
	dBodySetPosition(body.Body, xPos, yPos, zPos);
	dMatrix3 Orient3;
	//dRFromAxisAndAngle(Orient3, 0, 0, 1, 3.14/4);
	dRFromEulerAngles(Orient3, xRot, yRot, zRot);
	//dRFromEulerAngles(Orient3, 0, 0, 0);
	dBodySetRotation(body.Body, Orient3);
	dBodySetLinearVel(body.Body, 0, 0, 0);
	dBodySetData(body.Body, (void *)0);
	dMass bodyMass;
	dMassSetCapsule(&bodyMass, 1.0, 3, radius, length);
	//dMassSetBox(&bodyMass, 10, sides[0], sides[1], sides[2]);
	dBodySetMass(body.Body, &bodyMass);
	body.Geom = dCreateCapsule(Space, radius, length);
	//body.Geom = dCreateBox(Space, sides[0], sides[1], sides[2]);
	dGeomSetBody(body.Geom, body.Body);

	float head_pos[ 3 ] = { 0.0f, 5.0f, -1.0f };
	createInvisibleHead( head_pos );
	
	createFrontLegs();
	createMiddleLegs();
	createBackLegs();

	rng_engine.seed( std::chrono::duration_cast< std::chrono::microseconds >( std::chrono::system_clock::now().time_since_epoch() ).count() );

	const int foodPrize = 10;

	dReal position[3] = { 10.0f, 0.0f, -20.0f };
	addFoodParticle( position, foodPrize, &World, &Space );
}

void init_random_seed(int seed)
{
  generator.seed(seed);
}