text:
/* boost random/xor_combine.hpp header file
 *
 * Copyright Jens Maurer 2002
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * See http://www.boost.org for most recent version including documentation.
 *
 * $Id: xor_combine.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $
 *
 */

#ifndef BOOST_RANDOM_XOR_COMBINE_HPP
#define BOOST_RANDOM_XOR_COMBINE_HPP

#include <istream>
#include <iosfwd>
#include <cassert>
#include <algorithm> // for std::min and std::max
#include <boost/config.hpp>
#include <boost/limits.hpp>
#include <boost/cstdint.hpp>     // uint32_t
#include <boost/random/detail/config.hpp>
#include <boost/random/detail/seed.hpp>
#include <boost/random/detail/seed_impl.hpp>

namespace boost {
namespace random {

/**
 * Instantiations of @c xor_combine_engine model a
 * \pseudo_random_number_generator.  To produce its output it
 * invokes each of the base generators, shifts their results
 * and xors them together.
 */
template<class URNG1, int s1, class URNG2, int s2>
class xor_combine_engine
{
public:
    typedef URNG1 base1_type;
    typedef URNG2 base2_type;
    typedef typename base1_type::result_type result_type;

    BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
    BOOST_STATIC_CONSTANT(int, shift1 = s1);
    BOOST_STATIC_CONSTANT(int, shift2 = s2);

    /**
     * Constructors a @c xor_combine_engine by default constructing
     * both base generators.
     */
    xor_combine_engine() : _rng1(), _rng2() { }

    /** Constructs a @c xor_combine by copying two base generators. */
    xor_combine_engine(const base1_type & rng1, const base2_type & rng2)
      : _rng1(rng1), _rng2(rng2) { }

    /**
     * Constructs a @c xor_combine_engine, seeding both base generators
     * with @c v.
     *
     * @xmlwarning
     * The exact algorithm used by this function may change in the future.
     * @endxmlwarning
     */
    BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(xor_combine_engine,
        result_type, v)
    { seed(v); }

    /**
     * Constructs a @c xor_combine_engine, seeding both base generators
     * with values produced by @c seq.
     */
    BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(xor_combine_engine,
        SeedSeq, seq)
    { seed(seq); }

    /**
     * Constructs a @c xor_combine_engine, seeding both base generators
     * with values from the iterator range [first, last) and changes
     * first to point to the element after the last one used.  If there
     * are not enough elements in the range to seed both generators,
     * throws @c std::invalid_argument.
     */
    template<class It> xor_combine_engine(It& first, It last)
      : _rng1(first, last), _rng2( /* advanced by other call */ first, last) { }

    /** Calls @c seed() for both base generators. */
    void seed() { _rng1.seed(); _rng2.seed(); }

    /** @c seeds both base generators with @c v. */
    BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(xor_combine_engine, result_type, v)
    { _rng1.seed(v); _rng2.seed(v); }

    /** @c seeds both base generators with values produced by @c seq. */
    BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(xor_combine_engine, SeedSeq, seq)
    { _rng1.seed(seq); _rng2.seed(seq); }

    /**
     * seeds both base generators with values from the iterator
     * range [first, last) and changes first to point to the element
     * after the last one used.  If there are not enough elements in
     * the range to seed both generators, throws @c std::invalid_argument.
     */
    template<class It> void seed(It& first, It last)
    {
        _rng1.seed(first, last);
        _rng2.seed(first, last);
    }

    /** Returns the first base generator. */
    const base1_type& base1() const { return _rng1; }

    /** Returns the second base generator. */
    const base2_type& base2() const { return _rng2; }

    /** Returns the next value of the generator. */
    result_type operator()()
    {
        return (_rng1() << s1) ^ (_rng2() << s2);
    }

    /** Fills a range with random values */
    template<class Iter>
    void generate(Iter first, Iter last)
    { detail::generate_from_int(*this, first, last); }

    /** Advances the state of the generator by @c z. */
    void discard(boost::uintmax_t z)
    {
        _rng1.discard(z);
        _rng2.discard(z);
    }

    /** Returns the smallest value that the generator can produce. */
    static result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () { return (std::min)((URNG1::min)(), (URNG2::min)()); }
    /** Returns the largest value that the generator can produce. */
    static result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () { return (std::max)((URNG1::min)(), (URNG2::max)()); }

    /**
     * Writes the textual representation of the generator to a @c std::ostream.
     */
    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, xor_combine_engine, s)
    {
        os << s._rng1 << ' ' << s._rng2;
        return os;
    }

    /**
     * Reads the textual representation of the generator from a @c std::istream.
     */
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, xor_combine_engine, s)
    {
        is >> s._rng1 >> std::ws >> s._rng2;
        return is;
    }

    /** Returns true if the two generators will produce identical sequences. */
    BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(xor_combine_engine, x, y)
    { return x._rng1 == y._rng1 && x._rng2 == y._rng2; }

    /** Returns true if the two generators will produce different sequences. */
    BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(xor_combine_engine)

private:
    base1_type _rng1;
    base2_type _rng2;
};

#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
//  A definition is required even for integral static constants
template<class URNG1, int s1, class URNG2, int s2>
const bool xor_combine_engine<URNG1, s1, URNG2, s2>::has_fixed_range;
template<class URNG1, int s1, class URNG2, int s2>
const int xor_combine_engine<URNG1, s1, URNG2, s2>::shift1;
template<class URNG1, int s1, class URNG2, int s2>
const int xor_combine_engine<URNG1, s1, URNG2, s2>::shift2;
#endif

/// \cond show_private

/** Provided for backwards compatibility. */
template<class URNG1, int s1, class URNG2, int s2,
    typename URNG1::result_type v = 0>
class xor_combine : public xor_combine_engine<URNG1, s1, URNG2, s2>
{
    typedef xor_combine_engine<URNG1, s1, URNG2, s2> base_type;
public:
    typedef typename base_type::result_type result_type;
    xor_combine() {}
    xor_combine(result_type val) : base_type(val) {}
    template<class It>
    xor_combine(It& first, It last) : base_type(first, last) {}
    xor_combine(const URNG1 & rng1, const URNG2 & rng2)
      : base_type(rng1, rng2) { }

    result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return (std::min)((this->base1().min)(), (this->base2().min)()); }
    result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return (std::max)((this->base1().min)(), (this->base2().max)()); }
};

/// \endcond

} // namespace random
} // namespace boost

#endif // BOOST_RANDOM_XOR_COMBINE_HPP

max_stars_repo_name: PXLVision/opengv

max_issues_repo_name: isuhao/klib

max_forks_repo_name: isuhao/klib


text:
#include "wregoolar/ret.hpp"

#define BOOST_TEST_MODULE basic template operations
#include <boost/test/unit_test.hpp>

BOOST_AUTO_TEST_CASE( cset_cons_0 )
{
    using namespace wregoolar::ret;

    using R = cset<'a'>;
    using S = cset<>;
    using T = cons_cset<'a', S>::type;
    using U = std::is_same<R, T>;

    BOOST_CHECK_EQUAL(U::value, true);
}

BOOST_AUTO_TEST_CASE( cset_cons_1 )
{
    using namespace wregoolar::ret;

    using R = cset<'a', 'b', 'c'>;
    using S = cset<'b', 'c'>;
    using T = cons_cset<'a', S>::type;
    using U = std::is_same<R, T>;

    BOOST_CHECK_EQUAL(U::value, true);
}

BOOST_AUTO_TEST_CASE( cset_contains_0 )
{
    using namespace wregoolar::ret;

    using S = cset<>;
    BOOST_CHECK_EQUAL(S::contains<'a'>::value, false);
}

BOOST_AUTO_TEST_CASE( cset_contains_1 )
{
    using namespace wregoolar::ret;

    using S = cset<'a'>;
    BOOST_CHECK_EQUAL(S::contains<'a'>::value, true);
}

BOOST_AUTO_TEST_CASE( cset_contains_2 )
{
    using namespace wregoolar::ret;

    using S = cset<'b'>;
    BOOST_CHECK_EQUAL(S::contains<'a'>::value, false);
}

BOOST_AUTO_TEST_CASE( cset_contains_3 )
{
    using namespace wregoolar::ret;

    using S = cset<'a', 'b', 'c', 'd'>;
    BOOST_CHECK_EQUAL(S::contains<'a'>::value, true);
}

BOOST_AUTO_TEST_CASE( cset_contains_4 )
{
    using namespace wregoolar::ret;

    using S = cset<'a', 'b', 'c', 'd'>;
    BOOST_CHECK_EQUAL(S::contains<'c'>::value, true);
}

BOOST_AUTO_TEST_CASE( cset_is_contained_in_0 )
{
    using namespace wregoolar::ret;

    using S = cset<'a', 'b', 'c', 'd'>;
    using T = cset<'b', 'd'>;
    BOOST_CHECK_EQUAL(T::is_contained_in<S>::value, true);
    BOOST_CHECK_EQUAL(S::is_contained_in<T>::value, false);
}

BOOST_AUTO_TEST_CASE( equal_csets_0 )
{
    using namespace wregoolar::ret;

    using S = cset<'a', 'b', 'c', 'd'>;
    using T = cset<'b', 'd'>;
    using X = equal_csets<S,T>;
    BOOST_CHECK_EQUAL(X::value, false);
    using Y = equal_csets<T,S>;
    BOOST_CHECK_EQUAL(Y::value, false);
}

BOOST_AUTO_TEST_CASE( equal_csets_1 )
{
    using namespace wregoolar::ret;

    using S = cset<'a', 'b', 'c', 'd'>;
    using T = cset<'b', 'd', 'a', 'c'>;
    using X = equal_csets<S,T>;
    BOOST_CHECK_EQUAL(X::value, true);
    using Y = equal_csets<T,S>;
    BOOST_CHECK_EQUAL(Y::value, true);
}

BOOST_AUTO_TEST_CASE( intersection_0 )
{
    using namespace wregoolar::ret;

    using R = cset<'a', 'c'>;
    using S = cset<'a', 'b', 'c', 'd'>;
    using T = cset<'g', 'a', 'f', 'c'>;
    using X = intersection<S,T>::type;
    using Y = equal_csets<X,R>;
    static_assert(Y::value);
}

BOOST_AUTO_TEST_CASE( intersection_1 )
{
    using namespace wregoolar::ret;

    using R = cset<'a', 'c'>;
    using S = cset<'a', 'b', 'c', 'd'>;
    using T = cset<'g', 'a', 'f', 'c'>;
    using X = intersection<S,T>::type;
    using Y = equal_csets<X,R>;
    static_assert(Y::value);
}

BOOST_AUTO_TEST_CASE( union_0 )
{
    using namespace wregoolar::ret;

    using R = cset<'a', 'b', 'c', 'd', 'f', 'g'>;
    using S = cset<'a', 'b', 'c', 'd'>;
    using T = cset<'g', 'a', 'f', 'c'>;
    using X1 = union_<S,T>::type;
    using Y1 = equal_csets<X1,R>;
    static_assert(Y1::value);
    using X2 = union_<T,S>::type;
    using Y2 = equal_csets<X2,R>;
    static_assert(Y2::value);
}

BOOST_AUTO_TEST_CASE( has_epsilon_0 )
{
    using namespace wregoolar::ret;

    using R0 = empty;
    static_assert(R0::has_epsilon::value == false);

    using R1 = epsilon;
    static_assert(R1::has_epsilon::value == true);

    using R2 = atom<'a'>;
    static_assert(R2::has_epsilon::value == false);

    using R3 = star<atom<'a'>>;
    static_assert(R3::has_epsilon::value == true);

    using R4 = conj<atom<'a'>, atom<'b'>>;
    static_assert(R4::has_epsilon::value == false);

    using R5 = conj<star<atom<'a'>>, atom<'b'>>;
    static_assert(R5::has_epsilon::value == false);

    using R6 = conj<star<atom<'a'>>, star<atom<'b'>>>;
    static_assert(R6::has_epsilon::value == true);

    using R7 = disj<atom<'a'>, atom<'b'>>;
    static_assert(R7::has_epsilon::value == false);

    using R8 = disj<star<atom<'a'>>, atom<'b'>>;
    static_assert(R8::has_epsilon::value == true);

    using R9 = disj<star<atom<'a'>>, star<atom<'b'>>>;
    static_assert(R9::has_epsilon::value == true);
}

BOOST_AUTO_TEST_CASE( derivative_0 )
{
    using namespace wregoolar::ret;

    using R0 = atom<'a'>;

    using Sa = R0::derivative<'a'>;
    static_assert(std::is_same<Sa, epsilon>::value);

    using Sb = R0::derivative<'b'>;
    static_assert(std::is_same<Sb, empty>::value);
}

BOOST_AUTO_TEST_CASE( derivative_1 )
{
    using namespace wregoolar::ret;

    using R0 = conj<atom<'a'>, atom<'b'>>;

    using Sa = R0::derivative<'a'>;
    static_assert(std::is_same<Sa, atom<'b'>>::value);

    using Sb = R0::derivative<'b'>;
    static_assert(std::is_same<Sb, empty>::value);
}

BOOST_AUTO_TEST_CASE( derivative_2 )
{
    using namespace wregoolar::ret;

    using R0 = disj<atom<'a'>, atom<'b'>>;

    using Sa = R0::derivative<'a'>;
    static_assert(std::is_same<Sa, epsilon>::value);

    using Sb = R0::derivative<'b'>;
    static_assert(std::is_same<Sb, epsilon>::value);
}

ext: cpp

max_stars_repo_name: drtconway/wregoolar

max_issues_repo_name: drtconway/wregoolar

max_forks_repo_name: drtconway/wregoolar


text:
#define BOOST_TEST_MODULE pcraster model_engine objectlink
#include <boost/test/unit_test.hpp>
#include "geo_rasterspace.h"
#include "com_exception.h"
#include "calc_findsymbol.h"
#include "calc_operations.h"
#include "calc_operator.h"
#include "calc_runtimeenv.h"
#include "calc_spatial.h"
#include "calc_nonspatial.h"

// NOTE use string failureExpected in files expected to fail, see style guide

// getMeta()
#include "calcLibWrapper.cc"
//! to trace execution of ctor and methods
static size_t execCallMarker=0;

void calc::CalcLibDemoObjectLink::setDem(
    const REAL4         *dem)
{
  execCallMarker=100*(size_t)dem[0];
}
void calc::CalcLibDemoObjectLink::getDem(
    REAL4 *             dem)
{
  dem[2]=10;
}

void calc::CalcLibDemoObjectLink::noArguments()
{
}

void calc::CalcLibDemoObjectLink::testOrder(
    REAL4*                           result,
    const std::vector<const REAL4 *> dems,
    UINT1                            singleLdd)
{
  result[0]=dems[0][0];
  result[1]=dems[1][0];
  result[2]=singleLdd;
}

void calc::CalcLibDemoObjectLink::testOrder2(
    std::vector<REAL4*>&             result,
    const std::vector<const REAL4 *> dems,
    UINT1                            singleLdd) const
{
  PRECOND(result.size()==2);
  result[1][0]=dems[0][0];
  result[1][1]=dems[1][0];
  result[1][2]=singleLdd;
  for(size_t i=0; i < 3; ++i)
    result[0][i]=result[1][i]*10;
}

calc::CalcLibDemoObjectLink::CalcLibDemoObjectLink(
    const geo::RasterSpace& /* rs */)
{
  execCallMarker=99;
}

calc::CalcLibDemoObjectLink::~CalcLibDemoObjectLink()
{
  execCallMarker=77;
}

BOOST_AUTO_TEST_CASE(testLoadLink)
{
  using namespace calc;

  ObjectLinkMeta olm=getMeta();

  BOOST_CHECK( olm.className()=="CalcLibDemoObjectLink");
  ObjectLinkMeta::MethodMap::const_iterator pos;
  // ctor
  pos = olm.methodMap().find("");
  BOOST_CHECK(pos != olm.methodMap().end());

  pos = olm.methodMap().find("getDem");
  BOOST_CHECK(pos != olm.methodMap().end());

  pos = olm.methodMap().find("failureExpected");
  BOOST_CHECK(pos == olm.methodMap().end());
}

BOOST_AUTO_TEST_CASE(testLoadLink2)
{
  using namespace calc;

  // multiple loads harmless?
  globalOperations.load(getMeta);

  bool catched=false;
  try {
  globalOperations.load(getMeta);
  } catch(...) {
    catched=true;
  }
  BOOST_CHECK(catched);

  // ctor has className
  const Operator *o;
  o=globalOperations["CalcLibDemoObjectLink"];
  BOOST_CHECK(o);
  BOOST_CHECK(o->nrResults()       == 1);
  BOOST_CHECK(o->resultType(0).vs()== VS_OBJECT);

  // method
  o=globalOperations["CalcLibDemoObjectLink::getDem"];
  BOOST_CHECK(o);
  BOOST_CHECK(o->nrResults()       == 1);
  BOOST_CHECK(o->resultType(0).vs()== VS_S);

  //! unknown method
  BOOST_CHECK(globalOperations["CalcLibDemoObjectLink::failureExpected"] == 0);
}

BOOST_AUTO_TEST_CASE(testExec)
{
  using namespace calc;

 RunTimeEnv      rte(geo::RasterSpace(1,3));
 // test only the calls discarding RunTimeEnv info
 BOOST_CHECK(execCallMarker==0);

 ObjectLinkProxy<CalcLibDemoObjectLink>* t=
       new ObjectLinkProxy<CalcLibDemoObjectLink>("",rte.rasterSpace(),0);

 {
   std::vector<Field *> f;
   REAL4 demIn[3]= { 4,8,16};
   f.push_back(new Spatial(VS_S,demIn,3));
   t->exec1("setDem",f);
   BOOST_CHECK(execCallMarker==400);
   delete f[0];
 }
 {
   std::vector<Field *> f;
   f.push_back(new Spatial(VS_S,CRI_f,3));
   f[0]->dest_f()[2] = -1;
   t->exec1("getDem",f);
   BOOST_CHECK(f[0]->src_f()[2] == 10);
   delete f[0];
 }
 {
   // test order of arguments
   std::vector<Field *> f;
   f.push_back(new Spatial(VS_S,CRI_f,3)); // result

   REAL4 demIn[3]= { 10,16,16};
   f.push_back(new Spatial(VS_S,demIn,3)); // dems 0
   demIn[0]= 20;
   f.push_back(new Spatial(VS_S,demIn,3)); // dems 1

   UINT1 l(5);
   f.push_back(new NonSpatial(VS_L,l));    // bool arg
   BOOST_CHECK(f.back()->cri() == CRI_1);
   t->exec1("testOrder",f);

   BOOST_CHECK(f[0]->src_f()[0] == 10);
   BOOST_CHECK(f[0]->src_f()[1] == 20);
   BOOST_CHECK(f[0]->src_f()[2] == 5);
   for (size_t i=0; i < f.size(); ++i)
     delete f[i];
 }

 bool catched=false;
 // unknown method
 try {
  std::vector<Field *> f;
  t->exec1("failureExpected",f);
 } catch(const ObjectLink::UnknownMethod& ) {
   catched=true;
 }
 BOOST_CHECK(catched);

 catched=false;
 // not enough arguments
 try {
  std::vector<Field *> f;
  t->exec1("testOrder",f);
 } catch(const std::out_of_range&) {
   catched=true;
 }
 BOOST_CHECK(catched);

 delete t;
 BOOST_CHECK(execCallMarker==77);
}

BOOST_AUTO_TEST_CASE(testExec2)
{
  using namespace calc;

  {
   RunTimeEnv      rte(geo::RasterSpace(2,2));

   // ctor
   const Operator *o;
   o=opName2op("CalcLibDemoObjectLink");
   BOOST_CHECK(o);
   o->exec(&rte,0);
   BOOST_CHECK(execCallMarker==99);

   // pop created ObjectLink
   DataValue *dv = rte.popDataValue();

   // now no ObjectLink present
   bool catched(false);
   try {
    o=opName2op("CalcLibDemoObjectLink::getDem");
    BOOST_CHECK(o);
    o->exec(&rte,0);
   } catch(const com::Exception& e) {
BOOST_CHECK(e.messages().find( "no ObjectLink present") != std::string::npos);
     catched=true;
   }
   BOOST_CHECK(catched);

   // delete ObjectLink
   delete dv;
   BOOST_CHECK(execCallMarker==77);
  }
}

BOOST_AUTO_TEST_CASE(testNoArguments)
{
  using namespace calc;

  RunTimeEnv rte(geo::RasterSpace(2,2));
  Operator const* op = opName2op("CalcLibDemoObjectLink");
  BOOST_CHECK(op);
  op->exec(&rte, 0);
  op = opName2op("CalcLibDemoObjectLink::noArguments");
  BOOST_CHECK(op);
  DataValue *dv = rte.popDataValue();
  BOOST_CHECK(dv);
  delete dv;
}

BOOST_AUTO_TEST_CASE(testCheckAndExec)
{
  using namespace calc;

 RunTimeEnv      rte(geo::RasterSpace(1,3));

 // ctor
 const Operator *o;
 o=globalOperations["CalcLibDemoObjectLink"];
 BOOST_CHECK(o);

 // wrong ctor call
 bool catched(false);
 try {
   // use createResultField to create the (wrong) input argument
  rte.pushField(rte.createResultField(DataType(VS_S,ST_NON)));
  o->checkAndExec(&rte,1);
 } catch (const com::Exception& e) {
BOOST_CHECK(e.messages().find(
         "'CalcLibDemoObjectLink' too many arguments specified") != std::string::npos);
   catched=true;
 }
 BOOST_CHECK(catched);

 // good ctor
 o->checkAndExec(&rte,0);

 // pop created ObjectLink
 DataValue *dv = rte.popDataValue();
 BOOST_CHECK(dynamic_cast<ObjectLink*>(dv));

 {
   REAL4 demIn[3]= { 10,16,16};
   rte.pushDataValue(dv);
   rte.pushField(new Spatial(VS_S,demIn,3));
   o=opName2op("CalcLibDemoObjectLink::setDem");
   BOOST_CHECK(o);
   BOOST_CHECK(o->nrResults()==0);

   o->checkAndExec(&rte,1);
 }

 // exec something
 {
   rte.pushDataValue(dv);

   REAL4 demIn[3]= { 10,16,16};
   rte.pushField(new Spatial(VS_S,demIn,3)); // dems 0
   demIn[0]= 20;
   rte.pushField(new Spatial(VS_S,demIn,3)); // dems 1

   UINT1 l(5);
   rte.pushField(new NonSpatial(VS_L,l));    // bool arg

   o=globalOperations["CalcLibDemoObjectLink::testOrder"];
   BOOST_CHECK(o);
   o->checkAndExec(&rte,3);

   Field *res= rte.popField();

   BOOST_CHECK(res->isSpatial());
   BOOST_CHECK(res->cri()==CRI_f);
   BOOST_CHECK(res->src_f()[0] == 10);
   BOOST_CHECK(res->src_f()[1] == 20);
   BOOST_CHECK(res->src_f()[2] == 5);
   delete res;
 }

 {
   rte.pushDataValue(dv);

   REAL4 demIn[3]= { 10,16,16};
   rte.pushField(new Spatial(VS_S,demIn,3)); // dems 0
   demIn[0]= 20;
   rte.pushField(new Spatial(VS_S,demIn,3)); // dems 1

   UINT1 l(5);
   rte.pushField(new NonSpatial(VS_L,l));    // bool arg

   o=globalOperations["CalcLibDemoObjectLink::testOrder2"];
   BOOST_CHECK(o);
   o->checkAndExec(&rte,3);

   // reversed stack order
   Field *res1= rte.popField();
   Field *res0= rte.popField();

   BOOST_CHECK(res1->isSpatial());
   BOOST_CHECK(res1->cri()==CRI_f);
   BOOST_CHECK(res1->src_f()[0] == 10);
   BOOST_CHECK(res1->src_f()[1] == 20);
   BOOST_CHECK(res1->src_f()[2] == 5);
   BOOST_CHECK(res0->isSpatial());
   BOOST_CHECK(res0->cri()==CRI_f);
   BOOST_CHECK(res0->src_f()[0] == 100);
   BOOST_CHECK(res0->src_f()[1] == 200);
   BOOST_CHECK(res0->src_f()[2] == 50);
   delete res0;
   delete res1;
 }

 // now no ObjectLink present
 catched=false;
 try {
  o=globalOperations["CalcLibDemoObjectLink::getDem"];
  BOOST_CHECK(o);
  o->exec(&rte,0);
 } catch(const com::Exception& e) {
BOOST_CHECK(e.messages().find( "no ObjectLink present") != std::string::npos);
   catched=true;
 }
 BOOST_CHECK(catched);

 // delete ObjectLink
 delete dv;
 BOOST_CHECK(execCallMarker==77);
}

ext: cc

max_stars_repo_name: quanpands/wflow

max_issues_repo_name: quanpands/wflow

max_forks_repo_name: quanpands/wflow


text:
// This file is part of SWGANH which is released under the MIT license.
// See file LICENSE or go to http://swganh.com/LICENSE

#include "command_factory.h"

#include <boost/thread/lock_guard.hpp>

#include "swganh/app/swganh_kernel.h"
#include "swganh/service/service_manager.h"

#include "command_interface.h"
#include "command_properties.h"
#include "command_service_interface.h"

using namespace swganh::command;
using swganh::app::SwganhKernel;

CommandFactory::CommandFactory() {}

CommandFactory::~CommandFactory() {}

void CommandFactory::Initialize(SwganhKernel* kernel)
{
    kernel_ = kernel;
    command_service_ = kernel_->GetServiceManager()->GetService<CommandServiceInterface>("CommandService");
}

bool CommandFactory::IsRegistered(swganh::HashString command)
{
    return command_creators_.find(command) != command_creators_.end();
}

void CommandFactory::AddCommandCreator(swganh::HashString command, swganh::command::CommandCreator creator)
{
    auto properties = command_service_->FindPropertiesForCommand(command);

    if (!properties)
    {
        CommandProperties new_properties;
        new_properties.command_name = command.ident_string();
        new_properties.allow_in_locomotion = 0xffffffffffffffff;
        new_properties.allow_in_state = 0xffffffffffffffff;

        properties.reset(new_properties);
    }

    {
        boost::lock_guard<boost::mutex> lg(creators_mutex_);
        command_creators_.insert(std::make_pair(command, std::make_shared<CreatorData>(std::move(creator), *properties)));
    }
}

void CommandFactory::RemoveCommandCreator(swganh::HashString command)
{
    boost::lock_guard<boost::mutex> lg(creators_mutex_);
    command_creators_.erase(command);
}

std::shared_ptr<CommandInterface> CommandFactory::CreateCommand(swganh::HashString command)
{
    std::shared_ptr<CommandInterface> new_command = nullptr;

    boost::lock_guard<boost::mutex> lg(creators_mutex_);

    auto creators_iter = command_creators_.find(command);
    if (creators_iter != command_creators_.end())
    {
        auto& creator = creators_iter->second;
        new_command = creator->creator_func();
        new_command->Initialize(kernel_, creator->properties);
    }

    return new_command;
}

max_stars_repo_name: apathyboy/swganh

max_issues_repo_name: apathyboy/swganh

max_forks_repo_name: apathyboy/swganh


text:
//
// FILE NAME:       $HeadURL: svn+ssh://svn.cm.aol.com/advertising/adlearn/gen1/trunk/lib/cpp/DataProxy/ShellStreamTransformer.hpp $
//
// REVISION:        $Revision: 281532 $
//
// COPYRIGHT:       (c) 2008 Advertising.com All Rights Reserved.
//
// LAST UPDATED:    $Date: 2013-06-19 21:05:52 -0400 (Wed, 19 Jun 2013) $
// UPDATED BY:      $Author: esaxe $

#ifndef _SHELL_STREAM_TRANSFORMER_HPP_
#define _SHELL_STREAM_TRANSFORMER_HPP_

#include "MVException.hpp"
#include "ITransformFunction.hpp"
#include <boost/shared_ptr.hpp>
#include <istream>
#include <map>
#include <string>

MV_MAKEEXCEPTIONCLASS( ShellStreamTransformerException, MVException );

class ShellStreamTransformer : public ITransformFunction
{
public:
	ShellStreamTransformer();
	virtual ~ShellStreamTransformer();

	virtual boost::shared_ptr< std::istream > TransformInput( boost::shared_ptr< std::istream > i_pInputStream, const std::map< std::string, std::string >& i_rParameters );
};

#endif //_SHELL_STREAM_TRANSFORMER_HPP_

ext: hpp

max_stars_repo_name: aol/DataProxyLibrary

max_issues_repo_name: aol/DataProxyLibrary

max_forks_repo_name: aol/DataProxyLibrary


text:
//==============================================================================
//         Copyright 2003 - 2011 LASMEA UMR 6602 CNRS/Univ. Clermont II
//         Copyright 2009 - 2011 LRI    UMR 8623 CNRS/Univ Paris Sud XI
//
//          Distributed under the Boost Software License, Version 1.0.
//                 See accompanying file LICENSE.txt or copy at
//                     http://www.boost.org/LICENSE_1_0.txt
//==============================================================================
#ifndef NT2_TRIGONOMETRIC_FUNCTIONS_GENERIC_REM_PIO2_MEDIUM_HPP_INCLUDED
#define NT2_TRIGONOMETRIC_FUNCTIONS_GENERIC_REM_PIO2_MEDIUM_HPP_INCLUDED
#include <nt2/trigonometric/functions/rem_pio2_medium.hpp>
#include <nt2/include/functions/simd/round.hpp>
#include <nt2/include/functions/simd/fast_toint.hpp>
#include <nt2/sdk/meta/as_integer.hpp>
#include <boost/fusion/include/std_pair.hpp>
#include <nt2/include/constants/three.hpp>
#include <nt2/include/constants/pio2_1.hpp>
#include <nt2/include/constants/pio2_1t.hpp>
#include <nt2/include/constants/pio2_2.hpp>
#include <nt2/include/constants/pio2_2t.hpp>
#include <nt2/include/constants/pio2_3.hpp>
#include <nt2/include/constants/pio2_3t.hpp>
#include <nt2/include/constants/invpio_2.hpp>
#include <nt2/include/functions/simd/minus.hpp>
#include <nt2/include/functions/simd/multiplies.hpp>
#include <nt2/include/functions/simd/bitwise_and.hpp>
#include <nt2/include/functions/simd/minus.hpp>

namespace nt2 { namespace ext
{
  NT2_FUNCTOR_IMPLEMENTATION(nt2::tag::rem_pio2_medium_, tag::cpu_,
                             (A0)(A1),
                             (generic_ < floating_<A0> > )
                             (generic_ < integer_<A1>  > )
                             (generic_ < floating_<A0> > )
                            )
  {
    typedef void result_type;
    inline result_type operator()(A0 const& x, A1 & n, A0 & xr)
    {
      n = nt2::rem_pio2_medium(x, xr);
    }
  };

  NT2_FUNCTOR_IMPLEMENTATION(nt2::tag::rem_pio2_medium_, tag::cpu_,
                             (A0),
                             (generic_ < floating_<A0> > )
                             )
  {
    typedef typename nt2::meta::as_integer<A0>::type int_t;
    typedef std::pair<int_t, A0>               result_type;
    NT2_FUNCTOR_CALL(1)
    {
      result_type res;
      res.first = nt2::rem_pio2_medium(a0,res.second);
      return res;
    }
  };

  NT2_FUNCTOR_IMPLEMENTATION(nt2::tag::rem_pio2_medium_, tag::cpu_,
                             (A0),
                             (generic_ < floating_<A0> > )
                             (generic_ < floating_<A0> > )
                             )
  {
    typedef typename meta::as_integer<A0>::type result_type;
    inline result_type operator()(A0 const& t, A0 & xr) const
    {
      const A0 fn = nt2::round(t*Invpio_2<A0>());
      A0 r  = t-fn*nt2::Pio2_1<A0>();
      A0 w  = fn*nt2::Pio2_1t<A0>();
      A0 t2 = r;
      w  = fn*nt2::Pio2_2<A0>();
      r  = t2-w;
      w  = fn*nt2::Pio2_2t<A0>()-((t2-r)-w);
      t2 = r;
      w  = fn*nt2::Pio2_3<A0>();
      r  = t2-w;
      w  = fn*nt2::Pio2_3t<A0>()-((t2-r)-w);
      xr = r-w;
      return  bitwise_and(fast_toint(fn), Three<result_type>());
    }
  };
} }
#endif


text:
/**
 *-----------------------------------------------------------------------------
 * Title      : Memory Slave
 * ----------------------------------------------------------------------------
 * File       : Slave.cpp
 * Created    : 2016-09-20
 * ----------------------------------------------------------------------------
 * Description:
 * Memory slave interface.
 * ----------------------------------------------------------------------------
 * This file is part of the rogue software platform. It is subject to
 * the license terms in the LICENSE.txt file found in the top-level directory
 * of this distribution and at:
 *    https://confluence.slac.stanford.edu/display/ppareg/LICENSE.html.
 * No part of the rogue software platform, including this file, may be
 * copied, modified, propagated, or distributed except according to the terms
 * contained in the LICENSE.txt file.
 * ----------------------------------------------------------------------------
**/
#include <rogue/interfaces/memory/Slave.h>
#include <rogue/interfaces/memory/Master.h>
#include <rogue/interfaces/memory/Constants.h>
#include <rogue/interfaces/memory/Transaction.h>
#include <rogue/GeneralError.h>
#include <memory>
#include <rogue/GilRelease.h>
#include <rogue/ScopedGil.h>

#ifndef NO_PYTHON
#define BOOST_BIND_GLOBAL_PLACEHOLDERS
#include <boost/python.hpp>
namespace bp  = boost::python;
#endif

namespace rim = rogue::interfaces::memory;

// Init class counter
uint32_t rim::Slave::classIdx_ = 0;

//! Class instance lock
std::mutex rim::Slave::classMtx_;

//! Create a slave container
rim::SlavePtr rim::Slave::create (uint32_t min, uint32_t max) {
   rim::SlavePtr s = std::make_shared<rim::Slave>(min,max);
   return(s);
}

//! Create object
rim::Slave::Slave(uint32_t min, uint32_t max) {
   min_ = min;
   max_ = max;

   classMtx_.lock();
   if ( classIdx_ == 0 ) classIdx_ = 1;
   id_ = classIdx_;
   classIdx_++;
   classMtx_.unlock();

   name_ = std::string("Unnamed_") + std::to_string(id_);
}

//! Destroy object
rim::Slave::~Slave() { }

//! Stop the interface
void rim::Slave::stop() {}

//! Register a master.
void rim::Slave::addTransaction(rim::TransactionPtr tran) {
   rogue::GilRelease noGil;
   std::lock_guard<std::mutex> lock(slaveMtx_);
   tranMap_[tran->id()] = tran;
}

//! Get transaction with index, called by sub classes
rim::TransactionPtr rim::Slave::getTransaction(uint32_t index) {
   rim::TransactionPtr ret;
   TransactionMap::iterator it;
   TransactionMap::iterator exp;

   rogue::GilRelease noGil;
   std::lock_guard<std::mutex> lock(slaveMtx_);

   if ( (it = tranMap_.find(index)) != tranMap_.end() ) {
      ret = it->second;

      // Remove from list
      tranMap_.erase(it);

      // Update any timers for transactions started after received transaction
      exp = tranMap_.end();
      for ( it = tranMap_.begin(); it != tranMap_.end(); ++it ) {
         if ( it->second->expired() ) exp = it;
         else it->second->refreshTimer(ret);
      }

      // Clean up if we found an expired transaction, overtime this will clean up
      // the list, even if it deletes one expired transaction per call
      if ( exp != tranMap_.end() ) {
          tranMap_.erase(exp);
      }

   }
   return ret;
}

//! Get min size from slave
uint32_t rim::Slave::min() {
   return min_;
}

//! Get min size from slave
uint32_t rim::Slave::max() {
   return max_;
}

//! Get id from slave
uint32_t rim::Slave::id() {
   return id_;
}

//! Get name from slave
std::string rim::Slave::name() {
   return name_;
}

//! Set name
void rim::Slave::setName(std::string name) {
   name_ = name;
}

//! Return id to requesting master
uint32_t rim::Slave::doSlaveId() {
   return(id());
}

//! Return name to requesting master
std::string rim::Slave::doSlaveName() {
   return(name());
}

//! Return min access size to requesting master
uint32_t rim::Slave::doMinAccess() {
   return(min());
}

//! Return max access size to requesting master
uint32_t rim::Slave::doMaxAccess() {
   return(max());
}

//! Return offset
uint64_t rim::Slave::doAddress() {
   return(0);
}

//! Post a transaction
void rim::Slave::doTransaction(rim::TransactionPtr transaction) {
   transaction->error("Unsupported transaction using unconnected memory bus");
}

void rim::Slave::setup_python() {
#ifndef NO_PYTHON
   bp::class_<rim::SlaveWrap, rim::SlaveWrapPtr, boost::noncopyable>("Slave",bp::init<uint32_t,uint32_t>())
      .def("setName",         &rim::Slave::setName)
      .def("_addTransaction", &rim::Slave::addTransaction)
      .def("_getTransaction", &rim::Slave::getTransaction)
      .def("_doMinAccess",    &rim::Slave::doMinAccess,   &rim::SlaveWrap::defDoMinAccess)
      .def("_doMaxAccess",    &rim::Slave::doMaxAccess,   &rim::SlaveWrap::defDoMaxAccess)
      .def("_doAddress",      &rim::Slave::doAddress,     &rim::SlaveWrap::defDoAddress)
      .def("_doTransaction",  &rim::Slave::doTransaction, &rim::SlaveWrap::defDoTransaction)
      .def("__lshift__",      &rim::Slave::lshiftPy)
      .def("_stop",           &rim::Slave::stop)
   ;
#endif
}

#ifndef NO_PYTHON

//! Constructor
rim::SlaveWrap::SlaveWrap(uint32_t min, uint32_t max) : rim::Slave(min,max) {}

//! Return min access size to requesting master
uint32_t rim::SlaveWrap::doMinAccess() {
   {
      rogue::ScopedGil gil;

      if (boost::python::override pb = this->get_override("_doMinAccess")) {
         try {
            return(pb());
         } catch (...) {
            PyErr_Print();
         }
      }
   }
   return(rim::Slave::doMinAccess());
}

//! Return min access size to requesting master
uint32_t rim::SlaveWrap::defDoMinAccess() {
   return(rim::Slave::doMinAccess());
}

//! Return max access size to requesting master
uint32_t rim::SlaveWrap::doMaxAccess() {
   {
      rogue::ScopedGil gil;

      if (boost::python::override pb = this->get_override("_doMaxAccess")) {
         try {
            return(pb());
         } catch (...) {
            PyErr_Print();
         }
      }
   }
   return(rim::Slave::doMaxAccess());
}

//! Return max access size to requesting master
uint32_t rim::SlaveWrap::defDoMaxAccess() {
   return(rim::Slave::doMinAccess());
}

//! Return offset
uint64_t rim::SlaveWrap::doAddress() {
   {
      rogue::ScopedGil gil;

      if (boost::python::override pb = this->get_override("_doAddress")) {
         try {
            return(pb());
         } catch (...) {
            PyErr_Print();
         }
      }
   }
   return(rim::Slave::doAddress());
}

//! Return offset
uint64_t rim::SlaveWrap::defDoAddress() {
   return(rim::Slave::doAddress());
}

//! Post a transaction. Master will call this method with the access attributes.
void rim::SlaveWrap::doTransaction(rim::TransactionPtr transaction) {
   {
      rogue::ScopedGil gil;

      if (boost::python::override pb = this->get_override("_doTransaction")) {
         try {
            pb(transaction);
            return;
         } catch (...) {
            PyErr_Print();
         }
      }
   }
   rim::Slave::doTransaction(transaction);
}

//! Post a transaction. Master will call this method with the access attributes.
void rim::SlaveWrap::defDoTransaction(rim::TransactionPtr transaction) {
   rim::Slave::doTransaction(transaction);
}

void rim::Slave::lshiftPy ( boost::python::object p ) {
   rim::MasterPtr mst;

   // First Attempt to access object as a memory master
   boost::python::extract<rim::MasterPtr> get_master(p);

   // Test extraction
   if ( get_master.check() ) mst = get_master();

   // Otherwise look for indirect call
   else if ( PyObject_HasAttrString(p.ptr(), "_getMemoryMaster" ) ) {

      // Attempt to convert returned object to master pointer
      boost::python::extract<rim::MasterPtr> get_master(p.attr("_getMemoryMaster")());

      // Test extraction
      if ( get_master.check() ) mst = get_master();
   }

   // Success
   if ( mst != NULL ) mst->setSlave(shared_from_this());
   else throw(rogue::GeneralError::create("memory::Slave::lshiftPy","Attempt to use << operator with incompatible memory master"));
}

#endif

//! Support << operator in C++
rim::MasterPtr & rim::Slave::operator <<(rim::MasterPtr & other) {
   other->setSlave(shared_from_this());
   return other;
}

max_stars_repo_name: slaclab/rouge

max_issues_repo_name: slaclab/rouge

max_forks_repo_name: slaclab/rouge


text:
//
// Copyright (C) 2019 by Yuri Victorovich. All rights reserved.
//

#include <stdlib.h>

#include <string>
#include <array>
#include <list>
#include <vector>
#include <set>
#include <iostream>
#include <iomanip>
#include <sstream>
#include <fstream>
#include <cmath>
#include <random>
#include <chrono>
#include <ctime>
#include <limits>
#include <csignal>
#include <exception>
#include <functional>
#include <stdexcept> // runtime_error

#include <cxxopts.hpp>
#include <nlohmann/json.hpp>
#include <boost/format.hpp>

#if USE_PARALLELISM
#include <thread>
#include "ThreadPool.h"
#endif

#include "Vec3.h"

//
// TODO
// * gravity
// * thermal flow
//

//
// types
//

typedef double Float;
typedef uint64_t Count;
typedef uint64_t CpuCycles;

//
// physics constants and laws
//

namespace PhysicsConsts {

static constexpr Float Na = 6.022e+23;  // Avogadro constant
static constexpr Float R = 8.3144598;   // 8.3144598(48) in kg m² s⁻² K−1 mol−1 , see https://en.wikipedia.org/wiki/Gas_constant
static constexpr Float k = R/Na;        // Boltzmann constant, see https://en.wikipedia.org/wiki/Boltzmann_constant
static constexpr Float Troom = 293.15;  // room temperature in Kelvins: 20°C=293.15K
static constexpr Float Patm = 101325;   // 101.325 kPa

}; // PhysicsConsts

namespace PhysicsFormulas {

using namespace PhysicsConsts;

static Float constexpr energyToTemperature(Float E) {return E/k/(3./2.);}
static Float constexpr temperatureToEnergy(Float T) {return k*T*(3./2.);}

}; // PhysicsFormulas

// quick gas computations:
// Vmol(0°C)  = 1mol*R*T=273.15K/101325 ≈ 0.0224m³ = 22.4dm³
// Vmol(20°C) = 1mol*R*T=293.15K/101325 ≈ 0.0241m³ = 24.1dm³
// V(1mil, 20°C) = (n=1mil/Na)RT/Patm ≈ 0.342e-6 (~0.3μm)

#if USE_PARALLELISM
static ThreadPool *threadPool = nullptr; // static thread pool user by all code here
#endif

#if DBG_TRACK_PARTICLES
#define DBG_TRACK_PARTICLE_MSG(msg) std::cout << "...[track]: " << msg << std::endl;
#define DBG_TRACK_PARTICLE_CODE(cmds ...) cmds
#else
#define DBG_TRACK_PARTICLE_MSG(msg) {/*do nothing*/}
#define DBG_TRACK_PARTICLE_CODE(cmds ...) /*do nothing*/
#endif

#if USE_PARALLELISM
#define TID "[" << ThreadPool::tid() << "]"
#else
#define TID ""
#endif

template<typename Tsrc, typename Tdst> Tdst convertType(const Tsrc &v);
template<> std::string convertType<std::string,std::string>(const std::string &v) {return v;}
template<> uint64_t convertType<std::string,uint64_t>(const std::string &v) {return std::stoull(v);}
template<> Float convertType<std::string,Float>(const std::string &v) {return std::stod(v);}

template<typename T>
static std::vector<T> splitString(const std::string &strToSplit, char delimeter) {
  std::stringstream ss(strToSplit);
  std::string item;
  std::vector<T> split;
  while (std::getline(ss, item, delimeter))
     split.push_back(convertType<std::string,T>(item));
  return split;
}

//
// helper classes
//

class exception : public std::exception {
  std::runtime_error m;
public:
  exception(const std::string &newMsg) : m(newMsg) { } // we use exception directly so constructor is public
  const char* what() const noexcept override {
    return m.what();
  }
}; // exception

class AOut : public std::ostringstream { // atomic message printer
public:
  ~AOut() {
    std::cout << str();
  }
};

//
// helper functions
//

namespace constexpr_funcs {

// borrowed from the answer in https://stackoverflow.com/questions/23781506/compile-time-computing-of-number-of-bits-needed-to-encode-n-different-states/23795363
constexpr unsigned floorlog2(unsigned x) {
  return x == 1 ? 0 : 1+floorlog2(x >> 1);
}

constexpr unsigned ceillog2(unsigned x) {
  return x == 1 ? 0 : floorlog2(x - 1) + 1;
}

namespace Detail {
  static double constexpr sqrtNewtonRaphson(double x, double curr, double prev) {
    return curr == prev
      ? curr
      : sqrtNewtonRaphson(x, 0.5 * (curr + x / curr), curr);
  }

  static double constexpr cbrtNewtonRaphson(double x, double curr, double prev) {
    return curr == prev
      ? curr
      : cbrtNewtonRaphson(x, curr - (curr*curr*curr - x)/(3*curr*curr), curr);
  }
}

static double constexpr sqrt(double x) {
  return x >= 0 && x < std::numeric_limits<double>::infinity()
    ? Detail::sqrtNewtonRaphson(x, x, 0)
    : std::numeric_limits<double>::quiet_NaN();
}

static double constexpr cbrt(double x) {
  return x >= 0 && x < std::numeric_limits<double>::infinity()
    ? Detail::cbrtNewtonRaphson(x, x, 0)
    : std::numeric_limits<double>::quiet_NaN();
}

unsigned constexpr lim1(unsigned i) {
  return i >= 1 ? i : 1;
}

}; // constexpr_funcs

namespace xasm {

  static inline CpuCycles getCpuCycles() {
    unsigned int lo,hi;
    __asm__ __volatile__ ("rdtsc" : "=a" (lo), "=d" (hi));
    return ((CpuCycles)hi << 32) | lo;
  }

}; // xasm

static void formatUInt64(uint64_t x, std::ostream &os) {
  if (x >= 1000) {
    formatUInt64(x/uint64_t(1000), os);
    os << ',';
    os << std::setw(3) << std::setfill('0') << x%1000;
  } else {
    os << x;
  }
}

static std::string formatUInt64(uint64_t x) {
  std::ostringstream ss;
  formatUInt64(x, ss);
  return ss.str();
}

//
// relevant physics laws
//

// Ideal gas law: https://en.wikipedia.org/wiki/Ideal_gas_law
// PV = nRT    // n is number of moles
// E = ½mv² = 3/2 kT      // relation between the average kinetic energy and temperature

//
// params
//
static constexpr Count N = 1000000;
static constexpr Count Ncreate = N; // (DEBUG) should be =N, but allows to override the number of particles for debugging
static constexpr Float particlesPerBucket = 1.; // average number of particles per bucket. Performance drops when the number is >1 (like 2)
static constexpr Float m = 6.646476410e-27; // He atomic mass in kg
static constexpr unsigned Pairwise_RandomizeRounds = 1500*N;
static constexpr Float SZa[3] = {4e-6,1e-7,1e-7};    // size of the box, when the particles are outside this box they reflect to the other side
static constexpr Float SZ(int i) {return SZa[i-1];}
static constexpr Float particleRadius = 140e-12; // He atomic radius
static constexpr Float particleRadius2 = (2*particleRadius)*(2*particleRadius);
static constexpr Float Teff = PhysicsConsts::Troom;  // effective temperature (same as default temperature)
static constexpr Float Vthermal = constexpr_funcs::sqrt(PhysicsFormulas::temperatureToEnergy(Teff)*2/m);  // thermal velocity
static constexpr Float Vcutoff = Vthermal*4.; // velocity over which there is a neglible number of particles XXX TODO 5. should be percentage estimate based
static constexpr Float penetrationCoefficient = 0.3; // fraction of the particle radius that we allow to be penetrated at worst, considering Vcutoff
static constexpr Float dt = particleRadius*penetrationCoefficient/Vcutoff;
static constexpr unsigned NumOutputBuckets = 200; // how many buckets we build
static constexpr Float InteractionPct = 0.001; // how much energy is transferred
#if DBG_SAVE_IMAGES
static constexpr std::array<Float,2> imageAreaX = {{0.45,0.55}};
static constexpr std::array<Float,2> imageAreaY = {{0.45,0.55}};
static constexpr std::array<Float,2> imageAreaZ = {{0.45,0.55}}; //{{0.5-particleRadius,0.5+particleRadius}};// {{0.45,0.55}};
#endif

// XXX XXX Where does this belong?
static Float maxwellSigma(Float T) {return std::sqrt(PhysicsConsts::k*T/m);}

//
// Misc variables
//

static bool signalOccurred = false;

//
// Stats of the run
//
#if USE_PARALLELISM
static std::mutex statsLock;  // protects all but statsNumBucketMoves that is protected by ParticlesIndex::instanceLock
#endif
static Count statsNumBucketMoves = 0;
static Count statsNumCollisionsPP = 0;    // particle-particle collisions
static Count statsNumCollisionsPW = 0;
static Count statsNumCollisionsPWtemp = 0; // PW collisions with temperature adjustment
static Float statsWallDeltaEnergy = 0;     // how much energy is injected through the wall (see the 'temperature' option, negative value means that the energy is taken)

//
// Distributions
//
/*
namespace Distributions {
namespace Maxwell {

using namespace PhysicsConsts;

static Float distributionPerVelocity(Float T, Float V) {
  auto V2 = V*V;
  return std::pow(m/(2*M_PI*k*T), 3./2.)*(4*M_PI*V2)*std::exp(-m*V2/(2*k*T));
}

}; // Maxwell
}; // Distributions
*/

//
// Wall functors
//

static Float wallFunctorElastic(Float v) {return -v;}
static std::function<Float(Float)> wallFunctors[3][2];

//
// Classes
//

class ParticlesIndex;

class Particle {
public: // data
  Vec3  pos;
  Vec3  v;
#if DBG_TRACK_PARTICLES
  unsigned pno;   // 1-based index
  bool     track;
#endif
public: // methods
  Particle() { }
  Particle(const Vec3 &newPos, const Vec3 &newV)
    : pos(newPos), v(newV) DBG_TRACK_PARTICLE_CODE(, track(false))
  {
  }
  ~Particle() { }
  Float energy() const {
    return m*v.len2()/2;
  }
  Float velocity() const {
    return v.len();
  }
  void move(Float dt) {
    auto collideWall = [](unsigned coord, Float &partCoord, Float &partVelocity, Float wall1, Float wall2) {
      bool collided = false;
      if (partCoord < wall1) {
        partCoord = wall1 + (wall1-partCoord);
        partVelocity = wallFunctors[coord-1][0](partVelocity);
        statsNumCollisionsPW++;
        collided = true;
      }
      if (partCoord > wall2) {
        partCoord = wall2 - (partCoord-wall2);
        partVelocity = wallFunctors[coord-1][1](partVelocity);
        statsNumCollisionsPW++;
        collided = true;
      }
      return collided;
    };
    // move the point
    Vec3 pt(pos(X) + dt*v(X), pos(Y) + dt*v(Y), pos(Z) + dt*v(Z));
    // collide with walls
    bool collided = false;
    collided |= collideWall(X, pt(X), v(X), 0.+particleRadius, SZ(X)-particleRadius);
    collided |= collideWall(Y, pt(Y), v(Y), 0.+particleRadius, SZ(Y)-particleRadius);
    collided |= collideWall(Z, pt(Z), v(Z), 0.+particleRadius, SZ(Z)-particleRadius);
#if DBG_TRACK_PARTICLES
    if (track && collided)
      DBG_TRACK_PARTICLE_MSG(str(boost::format("particle#%1% collided with a wall") % pno)) // TODO report which wall
#endif
    // move to point
    move(pt);
  }
  void move(const Vec3 &newPos);
  static Float energyToVelocity(Float E) { // a reverse of energy()
    return std::sqrt(E*2/m);
  }
  void addEnergy(Float deltaE) {
    auto E1 = energy();
    auto E2 = E1 + deltaE;
    Float r = std::sqrt(E2/E1);
    v *= r;
  }
  static Float adjustCountPerUnit(Float V, unsigned cnt) {
    // count adjusted for unit size in the phase space
    return Float(cnt)/(4/3*V*V);
  }
  bool collisionCourse(const Particle &other) const {
    auto va = (v+other.v)/2;
    return (v-va)*(other.v-va) < 0; // when velocities are in opposite directions the scalar product is negative
  }
  Float distance2(const Particle &other) const {
    return (pos-other.pos).len2();
  }
  void collide(Particle &other) {
    //AOut() << "collide > this=" << *this << " other=" << other << "Etotal=" << (energy()+other.energy()) << std::endl;
    Vec3 dv = (other.v - v).project((other.pos - pos).normalizeZ());
    v += dv;
    other.v -= dv;
    //AOut() << "collide < this=" << *this << " other=" << other << "Etotal=" << (energy()+other.energy()) << std::endl;
#if DBG_TRACK_PARTICLES
    if (track)
      DBG_TRACK_PARTICLE_MSG(str(boost::format("particle#%1% collided with a particle #%2% (particle#%3% is the current particle in collision)") % pno % other.pno % pno))
    if (other.track)
      DBG_TRACK_PARTICLE_MSG(str(boost::format("particle#%1% collided with a particle #%2% (particle#%3% is the other particle in collision)") % other.pno % pno % other.pno))
#endif
  }
  friend std::ostream& operator<<(std::ostream &os, const Particle &p) {
    os << "{pos=" << p.pos << " v=" << p.v << "}";
    return os;
  }
}; // Particle

class ParticlesIndex { // index of particles in XYZ slot space
public:
  class Slot : public std::list<Particle*> {
    friend std::ostream& operator<<(std::ostream &os, const Slot &s) {
      unsigned ix = (&s - &instance.index[0][0][0])/(NSpaceSlots[1]*NSpaceSlots[2]);
      unsigned iy = (&s - &instance.index[ix][0][0])/NSpaceSlots[2];
      unsigned iz = (&s - &instance.index[ix][iy][0]);
      os << "{slot[" << ix << "][" << iy << "][" << iz << "] sz=" << s.size() << "}";
      return os;
    }
  }; // Slot
  static constexpr unsigned NSpaceSlots[3] = {constexpr_funcs::lim1((unsigned)(SZ(X)/constexpr_funcs::cbrt(SZ(X)*SZ(Y)*SZ(Z)/N*particlesPerBucket)+0.5)),
                                              constexpr_funcs::lim1((unsigned)(SZ(Y)/constexpr_funcs::cbrt(SZ(X)*SZ(Y)*SZ(Z)/N*particlesPerBucket)+0.5)),
                                              constexpr_funcs::lim1((unsigned)(SZ(Z)/constexpr_funcs::cbrt(SZ(X)*SZ(Y)*SZ(Z)/N*particlesPerBucket)+0.5))};
  static constexpr Float SZdivNSpaceSlots[3] = {SZa[0]/NSpaceSlots[0], SZa[1]/NSpaceSlots[1], SZa[2]/NSpaceSlots[2]};
  // leads to a fractional average occupancy of the bucket, seems to be the choice leading to the fastest computation
  static ParticlesIndex instance;
#if USE_PARALLELISM
  static std::mutex instanceLock;
#endif
private:
  std::array<std::array<std::array<Slot,NSpaceSlots[2]>,NSpaceSlots[1]>,NSpaceSlots[0]> index;
public:
  ParticlesIndex() {
    std::cout << "ParticlesIndex: NSpaceSlots={" << NSpaceSlots[0] << "," << NSpaceSlots[1] << "," << NSpaceSlots[2] << "}"
              << " (total=" << formatUInt64(NSpaceSlots[0]*NSpaceSlots[1]*NSpaceSlots[2]) << ")" << std::endl;
  }
  auto& get() {return index;}
  void add(const Vec3 &pos, Particle *p) {
    add(findSlot(pos), p);
  }
  void add(Slot &slot, Particle *p) {
    slot.push_front(p);
  }
  void add(Particle *p) {
    add(findSlot(p), p);
  }
  void del(const Vec3 &pos, Particle *p) {
    del(findSlot(pos), p);
  }
  void del(Slot &slot, Particle *p) {
    for (auto i = slot.begin(); i != slot.end(); i++)
      if (*i == p) {
        slot.erase(i);
        return;
      }
    assert(0); // unreachable
  }
  void del(Particle *p) {
    del(findSlot(p), p);
  }
  Slot& findSlot(const Particle *p) {return findSlot(p->pos);}
  Slot& findSlot(const Vec3 &pos) {return index[coordToSlot(X, pos(X))][coordToSlot(Y, pos(Y))][coordToSlot(Z, pos(Z))];}
private:
  static unsigned coordToSlot(unsigned ci, Float c) {
    return c/SZdivNSpaceSlots[ci-1];
  }
}; // ParticlesIndex

ParticlesIndex ParticlesIndex::instance;
#if USE_PARALLELISM
std::mutex ParticlesIndex::instanceLock;
#endif

//
// History
//

class History {
  typedef std::chrono::time_point<std::chrono::system_clock> Tm;
  using json = nlohmann::json;
  class Record {
    std::string action;
    Tm tm1; // time it started
    Tm tm2; // time it ended
    std::map<std::string,std::string> props; // arbitrary properties of the run
    auto serialize() const {
      auto tmToStr = [](Tm tm) {
        std::time_t tm_c = std::chrono::system_clock::to_time_t(tm);
        std::ostringstream ss;
        ss << std::put_time(std::localtime(&tm_c), "%F %T");
        return ss.str();
      };
      json j = json::object();
      j["action"] = action;
      j["tm1"] = std::chrono::duration_cast<std::chrono::milliseconds>(tm1.time_since_epoch()).count();
      j["tm1-str"] = tmToStr(tm1);
      j["tm2"] = std::chrono::duration_cast<std::chrono::milliseconds>(tm2.time_since_epoch()).count();
      j["tm2-str"] = tmToStr(tm2);
      if (!props.empty()) {
        auto jprops = json::object();
        for (auto const &kv : props)
          jprops[kv.first] = kv.second;
        j["props"] = jprops;
      }
      return j;
    }
    static Record unserialize(const json &j) {
      Record rec;
      rec.action = j["action"];
      rec.tm1 = Tm(std::chrono::milliseconds(j["tm1"]));
      rec.tm2 = Tm(std::chrono::milliseconds(j["tm2"]));
      if (j.find("props") != j.end()) {
        auto jprops = j["props"];
        for (auto it = jprops.begin(), ite = jprops.end(); it != ite; it++)
          rec.props[it.key()] = it.value();
      }
      return rec;
    }
    friend class History;
  }; // Record
  std::vector<Record>   lst;
public:
  History() { }
  void begin(const std::string &action) {
    Record rec;
    rec.action = action;
    rec.tm1 = std::chrono::system_clock::now();
    lst.push_back(rec);
  }
  void addProperty(const std::string &key, const std::string &val) {
    lst.rbegin()->props[key] = val;
  }
  void end() {
    lst.rbegin()->tm2 = std::chrono::system_clock::now();
  }
  auto serialize() const {
    json j = json::array();
    for (auto rec : lst)
      j.insert(j.end(), rec.serialize());
    return j;
  }
  static History unserialize(const json &j) {
    History h;
    for (auto jr : j)
      h.lst.push_back(Record::unserialize(jr));
    return h;
  }
}; // History

static History history;

//
// Helper classes: iterators
//

namespace Iterate {

template<typename Fn>
static void CellNeighborsForward(int ix, int iy, int iz, Fn &&fn);

namespace Detail {

template<typename Fn>
static void ThroughOverlapsYZ(unsigned ix, Fn &&fn) {
  const auto& sx = ParticlesIndex::instance.get()[ix];
  for (unsigned iy = 0; iy < ParticlesIndex::NSpaceSlots[1]; iy++) {
    const auto& sy = sx[iy];
    for (unsigned iz = 0; iz < ParticlesIndex::NSpaceSlots[2]; iz++) {
      const auto& sz = sy[iz];
      // process all pairs: same bucket
      for (auto it1 = sz.begin(), ite = sz.end(); it1 != ite; it1++) {
        auto p1 = *it1;
        for (auto it2 = it1; ++it2 != ite;) {
          auto p2 = *it2;
#if DBG_TRACK_PARTICLES
          //if (p1->track || p2->track)
          //  DBG_TRACK_PARTICLE_MSG(str(boost::format("checking same-bucket (%1%) for a collision between particle#%2% (%3%) and particle#%4% (%5%) ...")
          //                                           % sz % p1->pno % *p1 % p2->pno % *p2))
#endif
          if (p1->distance2(*p2) <= particleRadius2)
            fn(p1, p2);
        }
      } // same bucket
      // process all pairs: cross-bucket
      CellNeighborsForward(ix,iy,iz, [&sz,fn](int ix, int iy, int iz) {
        auto &slot2 = ParticlesIndex::instance.get()[ix][iy][iz];
        for (auto it1 = sz.begin(), it1e = sz.end(); it1 != it1e; it1++) {
          auto p1 = *it1;
          for (auto it2 = slot2.begin(), it2e = slot2.end(); it2 != it2e; it2++) {
            auto p2 = *it2;
#if DBG_TRACK_PARTICLES
          //if (p1->track || p2->track)
          //  DBG_TRACK_PARTICLE_MSG(str(boost::format("checking cross-bucket (%1%->%2%) for a collision between particle#%3% (%4%) and particle#%5% (%6%) ...")
          //                                           % sz % slot2 % p1->pno % *p1 % p2->pno % *p2))
#endif
            if (p1->distance2(*p2) <= particleRadius2)
              fn(p1, p2);
          }
        }
      });
    }
  }
}

}; // Detail

template<typename Fn>
static void CellNeighborsForward(int ix, int iy, int iz, Fn &&fn) {
  for (auto [dix,diy,diz] : std::array<std::array<int,3>,13>( // all forward-oriented neighboring buckets: (3^3-1)/2 = 13
        {{{{ 0, 0,+1}}, {{ 0,+1,-1}}, {{ 0,+1, 0}}, {{ 0,+1,+1}},
          {{+1, 0,-1}}, {{+1, 0, 0}}, {{+1, 0,+1}}, {{+1,-1,-1}}, {{+1,-1, 0}}, {{+1,-1,+1}}, {{+1,+1,-1}}, {{+1,+1, 0}}, {{+1,+1,+1}}}}))
    if (/*0 <= ix+dix &&*/ ix+dix < (int)ParticlesIndex::NSpaceSlots[0] &&
        0 <= iy+diy &&     iy+diy < (int)ParticlesIndex::NSpaceSlots[1] &&
        0 <= iz+diz &&     iz+diz < (int)ParticlesIndex::NSpaceSlots[2])
      fn(ix+dix, iy+diy, iz+diz);
}

template<typename Fn>
static void ThroughOverlapsSer(Fn &&fn) {
  for (unsigned ix = 0; ix < ParticlesIndex::NSpaceSlots[0]; ix++)
    Detail::ThroughOverlapsYZ(ix, fn);
}

#if USE_PARALLELISM
template<typename Fn>
static void ThroughOverlapsPar(Fn &&fn) {
  unsigned dix = (ParticlesIndex::NSpaceSlots[0]+1)/NCPU;
  std::mutex lock;
  bool done[1+NCPU] = {false};
  auto fnStitchJob = [fn](unsigned xStitch) {
    //AOut() << TID << ">> stitch @ xStitch=" << xStitch << std::endl;
    Detail::ThroughOverlapsYZ(xStitch, fn);
    //AOut() << TID << "<< stitch @ xStitch=" << xStitch << std::endl;
  };
  for (unsigned ix = 0, cpu = 1; ix < ParticlesIndex::NSpaceSlots[0]; cpu++) {
    unsigned ix1 = ix + dix, ixe;
    if (cpu == NCPU && ix1 != ParticlesIndex::NSpaceSlots[0])
      ix1 = ParticlesIndex::NSpaceSlots[0];
    ixe = cpu < NCPU ? ix1-1 : ix1;
    threadPool->doJob([cpu, ix, ixe, fn, fnStitchJob, &done, &lock]() {
      //AOut() << TID << ">> {cpu#" << cpu << "} main @ ix=" << ix << ".." << ixe-1 << std::endl;
      for (unsigned x = ix; x < ixe; x++)
        Detail::ThroughOverlapsYZ(x, fn);
      // done
      std::unique_lock<std::mutex> l(lock);
      done[cpu] = true;
      if (cpu > 1 && done[cpu-1])
        // stitch back
        threadPool->doJob([ix,fnStitchJob]() {
          fnStitchJob(ix-1);
        });
      if (cpu < NCPU && done[cpu+1])
        // stitch front
        threadPool->doJob([ixe,fnStitchJob]() {
          fnStitchJob(ixe);
        });
      //AOut() << TID << "<< main @ ix=" << ix << ".." << ixe-1 << std::endl;
    });
    ix = ix1;
  }
  // wait for all tasks to finish
  //AOut() << ">> waitForAll" << std::endl;
  threadPool->waitForAll();
  //AOut() << "<< waitForAll" << std::endl;
}
#endif

template<typename Fn>
static void ThroughOverlaps(Fn &&fn) {
#if USE_PARALLELISM
  if (2*NCPU-1 <= ParticlesIndex::NSpaceSlots[0]) // need to have at least 2 rows in each, except for the last one
    ThroughOverlapsPar(fn);
  else
    ThroughOverlapsSer(fn);
#else
  ThroughOverlapsSer(fn);
#endif
}

}; // Iterate

//
// Separate methods
//

void Particle::move(const Vec3 &newPos) {
  auto &slot1 = ParticlesIndex::instance.findSlot(this);
  pos = newPos;
  auto &slot2 = ParticlesIndex::instance.findSlot(this);
  if (&slot1 != &slot2) {
#if USE_PARALLELISM
    std::unique_lock<std::mutex> l(ParticlesIndex::instanceLock);
#endif
    ParticlesIndex::instance.del(slot1, this);
    ParticlesIndex::instance.add(slot2, this);
#if DBG_TRACK_PARTICLES
    if (track)
      DBG_TRACK_PARTICLE_MSG(str(boost::format("particle#%1% move between slots %2% -> %3%") % pno % slot1 % slot2)) // TODO find/keep slot numbers and print them
#endif
    //AOut() << "MOVE p=" << this << " x=" << x << " y=" << y << " z=" << z << std::endl;
    statsNumBucketMoves++;
  }
}

class DataBucket {
public:
  Float    V;    // velocity
  Count    cnt;  // how many particles fall into this bucket
  DataBucket(Float newV) : V(newV), cnt(0) { }
}; // DataBucket

//
// Data
//

static std::array<Particle, Ncreate> particles;

//
// Random value generator
//

class Random {
  // types
  typedef std::uniform_real_distribution<Float> urdFloat;
  typedef std::uniform_int_distribution<unsigned> uidUnsigned;
  typedef std::normal_distribution<Float> normFloat;
  // fields
  unsigned                                 seed;
  std::mt19937                             generator;
  urdFloat     uniformFloat01;
  uidUnsigned  uniform1N;
  uidUnsigned  uniformUInt01;
  urdFloat     uniformCoord[3]; // per dimension
  urdFloat     uniformEnergy;
  normFloat    normal;
public:
  Random()
  : seed(std::chrono::system_clock::now().time_since_epoch().count()),
    generator(seed),
    uniformFloat01(0.0, 1.0),
    uniform1N(1,Ncreate),
    uniformUInt01(0,1),
    uniformCoord{urdFloat(0.0+particleRadius, SZ(X)-particleRadius),
                 urdFloat(0.0+particleRadius, SZ(Y)-particleRadius),
                 urdFloat(0.0+particleRadius, SZ(Z)-particleRadius)},
    uniformEnergy(PhysicsFormulas::temperatureToEnergy(Teff)*0.95, PhysicsFormulas::temperatureToEnergy(Teff)*1.05)
  { }
  auto rand01() {
    return uniformFloat01(generator);
  }
  auto boolean() {
    return uniformUInt01(generator) == 1 ? true : false;
  }
  auto sphereAngles() {
    // from http://corysimon.github.io/articles/uniformdistn-on-sphere/
    Float theta = 2 * M_PI * rand01();
    Float phi = std::acos(1 - 2 * rand01());
    return std::array<Float,2>({{theta,phi}});
  }
  auto sphereCoords() {
    // from http://corysimon.github.io/articles/uniformdistn-on-sphere/
    auto [theta,phi] = sphereAngles();
    // convert to x/y/z
    double x = std::sin(phi)*std::cos(theta);
    double y = std::sin(phi)*std::sin(theta);
    double z = cos(phi);
    //
    return std::array<Float,3>({{x,y,z}});
  }
  auto coord(unsigned idx) {
    return uniformCoord[idx](generator);
  }
  auto energyUniformArea() {
    // one area
    return uniformEnergy(generator);
  }
  auto energySpike50_150() {
    // more dissipated pattern
    return PhysicsFormulas::temperatureToEnergy(Teff)*(boolean() ? 0.5 : 1.5);
  }
  auto particlePair() {
    unsigned i1, i2;
    while ((i1 = uniform1N(generator)) == (i2 = uniform1N(generator))) { }
    return std::array<unsigned,2>({{i1,i2}});
  }
  Vec3 maxwell3(Float sigma) {
    return Vec3(sigma*normal(generator), sigma*normal(generator), sigma*normal(generator));
  }
  Float maxwell1Plus(Float sigma) {
    auto v = sigma*normal(generator);
    return v >= 0 ? v : -v;
  }
}; // Random

static Random rg;

template<typename VGen>
static void generateParticles(VGen &&vgen) {
  // record in history
  history.begin(str(boost::format("generate %1% particles") % Ncreate));

  // TODO eliminate initial overlaps
  auto genParticleEnergyRange = [vgen]() {
    return Particle(Vec3(rg.coord(0), rg.coord(1), rg.coord(2)), vgen());
  };
  for (Count i = 0; i < Ncreate; i++)
    particles[i] = genParticleEnergyRange();

#if DBG_TRACK_PARTICLES
  {
    //Float v = 0.001/*per-frame evolution percentage*/*0.08/*imageAreaX span*/ / 0.00001/*dt*/;
    //particles[0] = Particle(Vec3(0.5,0.5,0.5), Vec3(v,v,0));
    //particles[0].track = true;
    //particles[1] = Particle(Vec3(0.5+0.053+particleRadius,0.5+0.053,0.5), Vec3(-v,-v,0));
    //particles[2] = Particle(Vec3(-0.05+0.5,-0.05+0.5,0.5), Vec3(2*v,2*v,0));
  }
#endif

  // add to index
  for (auto &p : particles)
    ParticlesIndex::instance.add(p.pos, &p);

  // remove overlaps
  bool hadOverlaps;
  do {
    hadOverlaps = false;
    // find overlaps
    std::set<Particle*> overlaps;
    Iterate::ThroughOverlapsSer([&overlaps](Particle *p1, Particle *p2) {
      overlaps.insert(p2);
    });
    // regenerate collided particles
    hadOverlaps = !overlaps.empty();
    for (auto p : overlaps) {
      ParticlesIndex::instance.del(p);
      *p = genParticleEnergyRange();
      ParticlesIndex::instance.add(p);
    }
    AOut() << "generate.overlap.iteration: overlapsCount=" << overlaps.size() << std::endl;
  } while (hadOverlaps);

  // record in history
  history.end();
}

static std::vector<DataBucket> particlesToBuckets(const std::array<Particle,Ncreate> &particles) {
  std::array<Float,2> energyLimits = {{0., 6.*PhysicsFormulas::temperatureToEnergy(Teff)}}; // TODO replace the ad-hoc 6. with the analytical expression based on the percentage of ignored particles
  Float minV = Particle::energyToVelocity(energyLimits[0]);
  Float maxV = Particle::energyToVelocity(energyLimits[1]);
  Float deltaV = (maxV - minV)/NumOutputBuckets;

  // generate initial buckets
  std::vector<DataBucket> buckets;
  for (unsigned i = 0; i < NumOutputBuckets; i++)
    buckets.push_back(DataBucket(minV + deltaV*i + deltaV/2));

  // count particles
  for (auto const &p : particles) {
    auto V = p.velocity();
    if (V < minV)
      continue;
    unsigned bucketNo = (V - minV)/deltaV;
    if (bucketNo >= NumOutputBuckets)
      continue; // discard out-of-range particles
    buckets[bucketNo].cnt++;
  }

  return buckets;
}

template<typename Rec, typename FnPrn>
static void writeFile(const std::string &fname, const std::vector<Rec> &recs, FnPrn&& prn) {
  std::ofstream file;
  file.open(fname);
  for (auto const &r : recs)
    file << prn(r) << std::endl;
  file.close();
}

static void transferPercentageOfEnergy(Particle &p1, Particle &p2, Float frac) {
  // p1 transfers pct of enerty to p2
  //AOut() << "transferPercentageOfEnergy > p1.e=" << p1.energy() << " p2.e=" << p2.energy() << std::endl;
  auto deltaE = p1.energy()*frac;
  p1.addEnergy(-deltaE);
  p2.addEnergy(+deltaE);
  //AOut() << "transferPercentageOfEnergy < p1.e=" << p1.energy() << " p2.e=" << p2.energy() << std::endl;
}

namespace ParticlesTotal {

static Float energy() {
  Float e = 0.;
  for (auto &p : particles)
    e += p.energy();
  return e;
}

static Float temperature() {
  return PhysicsFormulas::energyToTemperature(energy()/Ncreate);
}

}; // ParticlesTotal

//
// Images (save particles as an image for visual inspection)
//

#if DBG_SAVE_IMAGES

// tried CImg, but it is incompatible, because images having more than 1 bits per pixel are handled by calling external apps, like ImageMagick
#include <FreeImage.h>

class ImageSaver {
  // consts
  enum {SzX = 512, SzY = 512};
  // sub-area bounds that we image (we can't image all pixels because there's a million of them)
public:
  ImageSaver() { // have a constructor primarily for initialize/deinitialize functions
    FreeImage_Initialise();
    // delete previous images
    system("rm -f image-t=*");
  }
  ~ImageSaver() {
    FreeImage_DeInitialise();
  }
  void save(Float t, unsigned digitsInTime) const {
    // fname
    char fname[32];
    sprintf(fname, "image-t=%.*lf.png", digitsInTime, t);
    // write images
    FIBITMAP *bitmap = FreeImage_Allocate(SzX,SzY,3*8/*BitsPerPixel*/);
    RGBQUAD pixel;
    pixel = {255,255,255};
    FreeImage_FillBackground(bitmap, &pixel);
    for (auto &p : particles) {
      if (inBounds(p)) {
        auto [x,y,zc/*0..199*/] = mapToImage(p);
        //pixel = {255,BYTE(255-zc-56),BYTE(255-zc-56)}; // BGR
        pixel = {0,0,0}; // BGR
#if DBG_TRACK_PARTICLES
        if (p.track) {
          pixel = {0,0,255};
          //AOut() << "Display tracked particle @t=" << t << ": " << p << std::endl;
        }
#endif
        FreeImage_SetPixelColor(bitmap,x,y,&pixel);
      }
    }
    FreeImage_Save(FIF_PNG, bitmap, fname, 0);
    FreeImage_Unload(bitmap);
  }
private: // internals
  static bool inBounds(const Particle &p) {
    return inBounds(p.pos(X), imageAreaX) &&
           inBounds(p.pos(Y), imageAreaY) &&
           inBounds(p.pos(Z), imageAreaZ);
  }
  static bool inBounds(Float particleCoord, const std::array<Float,2> &bounds) {
    return bounds[0] <= particleCoord && particleCoord < bounds[1];
  }
  static std::tuple<unsigned/*X: 0..SzX-1*/,unsigned/*Y: 0..SzY-1*/,unsigned/*Z: 0..199*/> mapToImage(const Particle &p) {
    return {mapToImage(p.pos(X), imageAreaX, SzX),
            mapToImage(p.pos(Y), imageAreaY, SzY),
            mapToImage(p.pos(Z), imageAreaZ, 200)};
  }
  static unsigned mapToImage(Float particleCoord, const std::array<Float,2> &bounds, unsigned Sz) {
    return (particleCoord - bounds[0])/(bounds[1] - bounds[0])*Sz;
  }
}; // ImageSaver

static ImageSaver imageSaver;
#endif

//
// Serializer/unserializer
//

class Serializer {
  using json = nlohmann::json;
public:
  static auto serialize() {
    json j = json::object();
    j["particles"] = serializeParticles();
    j["history"] = history.serialize();
    return j;
  }
  static void unserialize(const json &j) {
    unserializeParticles(j["particles"]);
    history = History::unserialize(j["history"]);
  }
private:
  static json serializeParticles() {
    auto serializeVec = [](const Vec3 &v){
      json j = json::array();
      j.insert(j.end(), v(X));
      j.insert(j.end(), v(Y));
      j.insert(j.end(), v(Z));
      return j;
    };

    json arrParticles = json::array();
    DBG_TRACK_PARTICLE_CODE(unsigned n = 1;)
    for (auto &p : particles) {
      json particle = json::object();
      DBG_TRACK_PARTICLE_CODE(particle["n"] = n++;)
      particle["x"] = serializeVec(p.pos);
      particle["v"] = serializeVec(p.v);
      arrParticles.insert(arrParticles.end(), particle);
    }

    return arrParticles;
  }
  static void unserializeParticles(const json &j) {
    auto unserializeVec = [](json &j){
      Vec3 v;
      v(X) = j[0];
      v(Y) = j[1];
      v(Z) = j[2];
      return v;
    };
    unsigned i = 0;
    for (auto &p : particles) {
      auto particle = j[i++];
      p = Particle(unserializeVec(particle["x"]), unserializeVec(particle["v"]));
#if DBG_TRACK_PARTICLES
      p.track = false; // not serialized
#endif
    }
  }
}; // Serializer

//
// Evolve
//

class Evolver {
public:
  static void evolvePairwiseVelocity() {
    for (unsigned i = 0; i < Pairwise_RandomizeRounds; i++) {
      auto pp = rg.particlePair();
      transferPercentageOfEnergy(particles[pp[0]-1], particles[pp[1]-1], InteractionPct);
    }
  }
  template<typename FnBefore, typename FnAfter>
  static void evolvePhysically(Count numCycles, FnBefore &&fnBeforeCycle, FnAfter &&fnAfterCycle) {
    // record in history
    history.begin(str(boost::format("evolve %1% particles, numCycles=%2%") % Ncreate % numCycles));
    history.addProperty("temperatureBefore",  str(boost::format("%1%") % ParticlesTotal::temperature()));
    // evolve
#if DBG_SAVE_IMAGES
    imageSaver.save(0./*t*/, 5/*digits in time*/);
#endif
    // evolve
    Float t = 0.;
    for (Count cycle = 1; cycle <= numCycles; cycle++) {
      fnBeforeCycle(cycle);

      //
      // move
      //

      moveIterateByParticle();
      //moveIterateByBucket();
      t += dt;

      //
      // iter-particle collisions
      //
      Iterate::ThroughOverlaps([](Particle *p1, Particle *p2) {
        assert(p1->collisionCourse(*p2)); // overshoot the center due to too high speed?
        p1->collide(*p2);
        statsNumCollisionsPP++;
      });
#if DBG_SAVE_IMAGES
      imageSaver.save(t, 5/*digits in time*/);
#endif
      // cycle is done: allow the caller to perform custom actions
      fnAfterCycle(cycle);
    }
    // record in history
    history.addProperty("numCycles",         str(boost::format("%1%") % numCycles));
    history.addProperty("dt",                str(boost::format("%1%") % dt));
    history.addProperty("temperatureAfter",  str(boost::format("%1%") % ParticlesTotal::temperature()));
    history.end();
  }
private:
  static void moveIterateByParticle() { // faster when occupancy percentage is low
#if USE_PARALLELISM
    if (particles.size() >= NCPU)
      moveIterateByParticlePar();
    else
      moveIterateByParticleSeq();
#else
    moveIterateByParticleSeq();
#endif
  }
#if USE_PARALLELISM
  static void moveIterateByParticlePar() {
    Count ie = particles.size(), di = ie / NCPU;
    for (Count idx1 = 0, cpu = 1; idx1 < ie; cpu++) {
      Count idx2 = idx1 + di;
      if (idx2 > ie)
        idx2 = ie;
      if (cpu == NCPU && idx2 < ie)
        idx2 = ie;
      threadPool->doJob([idx1,idx2]() {
        for (Count i = idx1; i < idx2; i++)
          particles[i].move(dt);
      });
      idx1 = idx2;
    }
    // wait for all tasks to finish
    threadPool->waitForAll();
  }
#endif
  static void moveIterateByParticleSeq() {
    for (auto &p : particles)
      p.move(dt);
  }
  //static void moveIterateByBucketHelper(ParticlesIndex::Slot::iterator it, const ParticlesIndex::Slot::iterator &ite) {
  //  auto p = *it;
  //  it++;
  //  if (it != ite)
  //    moveIterateByBucketHelper(it, ite);
  //  p->move(dt, fnWallX1, fnWallX2, fnWallY1, fnWallY2, fnWallZ1, fnWallZ2);
  //}
  //static void moveIterateByBucket() { // slower when occupancy percentage is low
  //  for (auto slot = &ParticlesIndex::instance.get()[0][0][0], slote = slot + ParticlesIndex::NSpaceSlots[0]*ParticlesIndex::NSpaceSlots[1]*ParticlesIndex::NSpaceSlots[2]; slot < slote; slot++)
  //    if (!slot->empty()) {
  //      moveIterateByBucketHelper(slot->begin(), slot->end());
  //      //sz += slot->size();
  //    }
  //}
}; // Evolver

//
// MAIN
//

int mainGuarded(int argc, char *argv[]) {

  if (DBG_TRACK_PARTICLES) {
    std::cout << "!!!DEBUG CODE!!!" << std::endl;
    std::cout << "!!!DEBUG CODE!!! (DBG_TRACK_PARTICLES is enabled)" << std::endl;
    std::cout << "!!!DEBUG CODE!!!" << std::endl;
  }

#if USE_PARALLELISM
  threadPool = new ThreadPool(NCPU);
#endif

  //
  // local functions
  //

  auto fnSaveOutput = [](const std::string &fname) {
    AOut() << "serializing to " << fname << " ..." << std::endl;
    std::ofstream file;
    file.open(fname, std::ofstream::out);
    file << Serializer::serialize();
    file.close();
    AOut() << "done serializing to " << fname << " ..." << std::endl;
  };
  auto fnSaveBuckets = [](Count cycle) {
    writeFile<DataBucket>(str(boost::format("buckets-cycle%1%.txt") % cycle), particlesToBuckets(particles), [](const DataBucket &b) {
      std::ostringstream os;
      os << b.V << " " << b.cnt;
      return os.str();
    });
  };

  //
  // cycles
  //

  auto cpuCycles0 = xasm::getCpuCycles();

  //
  // parse arguments
  //

  cxxopts::Options options("Particle Simulator", "Simulator of the gas particles motion");
  options.add_options()
    ("d,distribution", "Velocity distribution used to generate particles, and to inject energy (maxwell,spike-50-150,uniform-area)", cxxopts::value<std::string>()->default_value("maxwell"))
    ("r,restart",      "Restart using the snapshot of particle positions/velocities", cxxopts::value<std::string>())
    ("o,output",       "Write the final particle state into this file (default would be set to particles-{time}.json)", cxxopts::value<std::string>()->default_value(
      str(boost::format("particles-%1%.json") % std::chrono::duration_cast<std::chrono::nanoseconds>(std::chrono::system_clock::now().time_since_epoch()).count())))
    ("c,cycles",       "How many cycles to perform (default is 4000)", cxxopts::value<Count>()->default_value("4000"))
    ("p,print",        "How frequently to print cycles (default is 1, which means print every cycle)", cxxopts::value<unsigned>()->default_value("1"))
    ("b,buckets",      "How frequently to print buckets (default is 0, which means print only in the end)", cxxopts::value<unsigned>()->default_value("0"))
    ("T,temperature",  "Set temperature of the walls, ex. -X:271:+X:272", cxxopts::value<std::string>()->default_value(""))
#if DBG_TRACK_PARTICLES
    ("t,track",        "Track particles: 1-based, colon-separated values", cxxopts::value<std::string>())
      // could use std::vector<unsigned> here to do '-t pno1 -t pno 2 ...', but maybe it's better (more compact) this way with colon-separated values
#endif
    ;

  auto result = options.parse(argc, argv);
  std::string optDistribution = result["distribution"].as<std::string>();
  if (optDistribution!="maxwell" && optDistribution!="spike-50-150" && optDistribution!="uniform-area")
    throw exception(str(boost::format("unknown velocity distribution '%1%' supplied to --generate argument") % optDistribution));
  bool optRestart = result.count("restart") > 0;
  std::string optRestartFile = optRestart ? result["restart"].as<std::string>() : "";
  std::string optOutputFile = result["output"].as<std::string>();
  Count numCycles = result["cycles"].as<Count>();
  unsigned cyclePrintPeriod = result["print"].as<unsigned>();
  unsigned cycleWriteBuckets = result["buckets"].as<unsigned>();
  auto optWallTemperatures = splitString<std::string>(result["temperature"].as<std::string>(), ':');
  if (optWallTemperatures.size()%2 != 0)
    throw exception(str(boost::format("wall temperature specifier is expected to consist of pairs, supplied '%1%' parameters") % optWallTemperatures.size()));

  Float optWallTempSigma[3][2] = {{0., 0.}, {0., 0.}, {0., 0.}};
  for (unsigned i = 0, ie = optWallTemperatures.size(); i!=ie; i++) {
    auto val = optWallTemperatures[i];
    if (val.size() != 2 || (val[0]!='-' && val[0]!='+') || (val[1]!='X' && val[1]!='Y' && val[1]!='Z'))
      throw exception(str(boost::format("wall temperature specifier is expected to have {sign}{Coord}, supplied '%1%'") % val));
    switch (val[1]) {
    case 'X': *(val[0]=='-' ? &optWallTempSigma[X-1][0] : &optWallTempSigma[X-1][1]) = maxwellSigma(convertType<std::string,Float>(optWallTemperatures[++i])); break;
    case 'Y': *(val[0]=='-' ? &optWallTempSigma[Y-1][0] : &optWallTempSigma[Y-1][1]) = maxwellSigma(convertType<std::string,Float>(optWallTemperatures[++i])); break;
    case 'Z': *(val[0]=='-' ? &optWallTempSigma[Z-1][0] : &optWallTempSigma[Z-1][1]) = maxwellSigma(convertType<std::string,Float>(optWallTemperatures[++i])); break;
    }
  }
#if DBG_TRACK_PARTICLES
  std::vector<Count> optTrack;
  if (result.count("track") > 0)
    optTrack = splitString<Count>(result["track"].as<std::string>(), ':');
  { // assign pno in particles
    unsigned pno = 1;
    for (auto &p : particles)
      p.pno = pno++;
  }
#endif

  //
  // checks
  //
  for (auto c : {X,Y,Z})
    assert(2*particleRadius < Float(SZ(c))/ParticlesIndex::NSpaceSlots[c-1]); // particle diameter should be < index slot size, because otherwise collision detection
                                                                              // needs to look out more than in 1 slot away which makes it impractical
  //
  // params & stats
  //

  AOut() << "params(init): dt=" << dt
                      << " Vthermal=" << Vthermal
                      << " numCycles=" << numCycles
                      << std::endl;
  AOut() << "stats(init): spacePercentageOccupiedByParticles=" << Ncreate*(4./3.*M_PI*std::pow(particleRadius,3))/(SZ(X)*SZ(Y)*SZ(Z))*100. << "%"
                     << " avgParticlePerBucket=" << Float(Ncreate)/(ParticlesIndex::NSpaceSlots[0]*ParticlesIndex::NSpaceSlots[1]*ParticlesIndex::NSpaceSlots[2])
                     << " P/Patm (at Teff)=" << ((Ncreate/PhysicsConsts::Na)*PhysicsConsts::R*Teff/(SZ(X)*SZ(Y)*SZ(Z))/PhysicsConsts::Patm)
                     << std::endl;

  //
  // generate or restart
  //
  if (optRestart) {
    std::cout << "unserializing from " << optRestartFile << " ..." << std::endl;
    std::ifstream file;
    file.open(optRestartFile, std::ifstream::in);
    Serializer::unserialize(nlohmann::json::parse(file));
    file.close();
    for (auto &p : particles)
      ParticlesIndex::instance.add(&p);
    std::cout << "done unserializing from " << optRestartFile << " ..." << std::endl;
  } else {
    if (optDistribution == "maxwell") {
      Float sigma = maxwellSigma(Teff);
      generateParticles([sigma]() {
        return rg.maxwell3(sigma);
      });
    } else if (optDistribution == "spike-50-150") {
      generateParticles([]() {
        auto [sphX,sphY,sphZ] = rg.sphereCoords();
        auto E = rg.energySpike50_150();
        auto V = Particle::energyToVelocity(E);
        return Vec3(sphX*V, sphY*V, sphZ*V);
      });
    } else if (optDistribution == "area") {
      generateParticles([]() {
        auto [sphX,sphY,sphZ] = rg.sphereCoords();
        auto E = rg.energyUniformArea();
        auto V = Particle::energyToVelocity(E);
        return Vec3(sphX*V, sphY*V, sphZ*V);
      });
    }
  }
#if DBG_TRACK_PARTICLES
  for (auto track : optTrack)
    if (track != 0)
      particles[track-1].track = true;
#endif

  //
  // initial log
  //
  {
    auto energy = ParticlesTotal::energy();
    std::cout << "log(init): energy(before)=" << energy
                        << " temperature(before)=" << PhysicsFormulas::energyToTemperature(energy/Ncreate)
                        << " numCycles=" << numCycles
                        << std::endl;
  }

  //
  // signal processing
  //

  std::signal(SIGHUP, [](int signum) { signalOccurred = true; });

  //
  // evolve
  //
  if (0)
    Evolver::evolvePairwiseVelocity();
  if (1) {
    unsigned prevStatsNumCollisionsPP = 0;
    auto cpuCyclesBefore = xasm::getCpuCycles();
    // fill wall collision functors
    for (auto c : std::array<unsigned,3>({{0,1,2}})) {
      if (optWallTempSigma[c][0]==0.)
        wallFunctors[c][0] = &wallFunctorElastic;
      else
        wallFunctors[c][0] = [=](Float v) {
          auto newV = rg.maxwell1Plus(optWallTempSigma[c][0]);
          Float deltaEnergy = m*(newV*newV - v*v)/2;
          {
#if USE_PARALLELISM
            std::unique_lock<std::mutex> l(statsLock);
#endif
            statsNumCollisionsPWtemp++;
            statsWallDeltaEnergy += deltaEnergy;
          }
          return +newV;
        };
      if (optWallTempSigma[c][1]==0.)
        wallFunctors[c][1] = &wallFunctorElastic;
      else
        wallFunctors[c][1] = [=](Float v) {
          auto newV = rg.maxwell1Plus(optWallTempSigma[c][1]);
          Float deltaEnergy = m*(newV*newV - v*v)/2;
          {
#if USE_PARALLELISM
            std::unique_lock<std::mutex> l(statsLock);
#endif
            statsNumCollisionsPWtemp++;
            statsWallDeltaEnergy += deltaEnergy;
          }
          return -newV;
        };
    }
    // evolve
    Evolver::evolvePhysically(numCycles, [&prevStatsNumCollisionsPP](unsigned cycle) {
        prevStatsNumCollisionsPP = statsNumCollisionsPP;
      }, [numCycles,cyclePrintPeriod,cycleWriteBuckets,&optOutputFile,&prevStatsNumCollisionsPP,&cpuCyclesBefore,fnSaveOutput,fnSaveBuckets](unsigned cycle) {
        // print tick and stats
        if (cycle%cyclePrintPeriod == 0) {
          Count cpuAtCycle = xasm::getCpuCycles();
          AOut() << "tick#" << cycle << ":evolvePhysically:"
                 << " avgCpuCyclesPerTick=" << formatUInt64((cpuAtCycle - cpuCyclesBefore)/Count(cycle))
                 << " statsNumCollisionsPP=" << statsNumCollisionsPP << " (+" << (statsNumCollisionsPP-prevStatsNumCollisionsPP) << ")"
                 << " statsNumCollisionsPW=" << statsNumCollisionsPW
                 << (statsNumCollisionsPWtemp>0 ? str(boost::format(" statsNumCollisionsPWtemp=%1%") % statsNumCollisionsPWtemp) : "")
                 << (statsWallDeltaEnergy!=0. ? str(boost::format(" statsWallDeltaEnergy=%1% (ΔT=%2%)")
                                                                  % statsWallDeltaEnergy
                                                                  % PhysicsFormulas::energyToTemperature(statsWallDeltaEnergy/Ncreate))
                                              : "")
                 << " statsNumBucketMoves=" << statsNumBucketMoves
                 << std::endl;
        }

        // save output if requested
        bool savedOnSignal = false;
        if (signalOccurred && cycle != numCycles) {
          fnSaveOutput(str(boost::format("%1%.cycle%2%.json") % optOutputFile % cycle));
          fnSaveBuckets(cycle);
          signalOccurred = false;
          savedOnSignal = true;
        }

        // write buckets
        if (!savedOnSignal && ((cycleWriteBuckets!=0 && cycle%cycleWriteBuckets==0) || cycle == numCycles))
          fnSaveBuckets(cycle);
      }
    );
  }

  //
  // final log & stats
  //
  {
    auto energy = ParticlesTotal::energy();
    std::cout << "log(fini): energy(after)=" << energy << std::endl;
    std::cout << "stats(fini): collisionsPerParticle(PP)=" << Float(statsNumCollisionsPP)/N
                          << " temperature(after)=" << PhysicsFormulas::energyToTemperature(energy/Ncreate)
                          << " collisions(PP)=" << formatUInt64(statsNumCollisionsPP)
                          << " collisions(PW)=" << formatUInt64(statsNumCollisionsPW)
                          << std::endl;
  }

  //
  // serialize for later reuse
  //

  fnSaveOutput(optOutputFile);
  if (numCycles == 0)
    fnSaveBuckets(0);

  //
  // cycles
  //
  auto cpuCyclesEnd = xasm::getCpuCycles();
  std::cout << "consumed cycles: " << formatUInt64(cpuCyclesEnd-cpuCycles0) << " {now at " << formatUInt64(cpuCyclesEnd) << "}" << std::endl;

#if USE_PARALLELISM
  delete threadPool;
#endif

  return 0;
}

int main(int argc, char *argv[]) {
  try {
    return mainGuarded(argc, argv);
  } catch (cxxopts::OptionException const& e) {
    std::cerr << "option parsing error: " << e.what() << std::endl;
    return 1;
  } catch (std::ifstream::failure const& e) {
    std::cerr << "file error: " << e.what() << std::endl;
    return 1;
  } catch (nlohmann::json::exception const& e) {
    std::cerr << "json error: " << e.what() << std::endl;
    return 1;
  } catch (exception const& e) {
    std::cerr << "application error: " << e.what() << std::endl;
    return 1;
  } catch (std::exception const& e) {
    std::cerr << "unknown exception of type '" << typeid(e).name() << "' caught: " << e.what() << std::endl;
    return 1;
  }
}

max_stars_repo_name: yurivict/gas-simulator

max_issues_repo_name: yurivict/gas-simulator

max_forks_repo_name: yurivict/gas-simulator