diff --git a/sprout/adapt/sscrisk/cel/array.hpp b/sprout/adapt/sscrisk/cel/array.hpp index 199dc515..74d3c34c 100644 --- a/sprout/adapt/sscrisk/cel/array.hpp +++ b/sprout/adapt/sscrisk/cel/array.hpp @@ -2,9 +2,15 @@ #define SPROUT_ADAPT_SSCRISK_CEL_ARRAY_HPP #include +#include #include #include #include +#if SPROUT_USE_INDEX_ITERATOR_IMPLEMENTATION +# include +# include +# include +#endif namespace sprout { // @@ -19,6 +25,85 @@ namespace sprout { typedef sscrisk::cel::array type; }; }; + +#if SPROUT_USE_INDEX_ITERATOR_IMPLEMENTATION + // + // fixed_container_traits + // + template + struct fixed_container_traits > + : public sprout::detail::fixed_container_traits_base > + { + public: + typedef sscrisk::cel::array fixed_container_type; + typedef fixed_container_type internal_type; + typedef internal_type clone_type; + typedef sprout::index_iterator iterator; + typedef sprout::index_iterator const_iterator; + public: + SPROUT_STATIC_CONSTEXPR typename sprout::detail::fixed_container_traits_base::size_type fixed_size + = std::tuple_size::type>::value + ; + }; + // + // fixed_container_traits + // + template + struct fixed_container_traits const> + : public sprout::detail::fixed_container_traits_base const> + { + public: + typedef sscrisk::cel::array const fixed_container_type; + typedef typename std::remove_const::type internal_type; + typedef internal_type clone_type; + typedef sprout::index_iterator iterator; + typedef sprout::index_iterator const_iterator; + public: + SPROUT_STATIC_CONSTEXPR typename sprout::detail::fixed_container_traits_base::size_type fixed_size + = std::tuple_size::type>::value + ; + }; + + // + // begin + // + template + typename sprout::fixed_container_traits >::iterator begin(sscrisk::cel::array& range) { + return typename sprout::fixed_container_traits >::iterator(range, 0); + } + template + SPROUT_CONSTEXPR typename sprout::fixed_container_traits >::const_iterator begin(sscrisk::cel::array const& range) { + return typename sprout::fixed_container_traits >::const_iterator(range, 0); + } + + // + // cbegin + // + template + SPROUT_CONSTEXPR typename sprout::fixed_container_traits >::const_iterator cbegin(sscrisk::cel::array const& range) { + return sprout::fixed_container_traits >::const_iterator(range, 0); + } + + // + // end + // + template + typename sprout::fixed_container_traits >::iterator end(sscrisk::cel::array& range) { + return typename sprout::fixed_container_traits >::iterator(range, range.size()); + } + template + SPROUT_CONSTEXPR typename sprout::fixed_container_traits >::const_iterator end(sscrisk::cel::array const& range) { + return typename sprout::fixed_container_traits >::const_iterator(range, range.size()); + } + + // + // cend + // + template + SPROUT_CONSTEXPR typename sprout::fixed_container_traits >::const_iterator cend(sscrisk::cel::array const& range) { + return typename sprout::fixed_container_traits >::const_iterator(range, range.size()); + } +#endif } // namespace sprout #endif // #ifndef SPROUT_ADAPT_SSCRISK_CEL_ARRAY_HPP diff --git a/sprout/random.hpp b/sprout/random.hpp index 8e18490f..4e2a1375 100644 --- a/sprout/random.hpp +++ b/sprout/random.hpp @@ -9,10 +9,9 @@ #include #include #include +#include #include #include #include #endif // #ifndef SPROUT_RANDOM_HPP - - diff --git a/sprout/random/binomial_distribution.hpp b/sprout/random/binomial_distribution.hpp new file mode 100644 index 00000000..009d4d28 --- /dev/null +++ b/sprout/random/binomial_distribution.hpp @@ -0,0 +1,425 @@ +#ifndef SPROUT_RANDOM_BINOMIAL_DISTRIBUTION_HPP +#define SPROUT_RANDOM_BINOMIAL_DISTRIBUTION_HPP + +#include +#include +#include +#include +#include +#include + +namespace sprout { + namespace random { + namespace detail { + template + struct binomial_table { + public: + SPROUT_STATIC_CONSTEXPR sprout::array table = sprout::array{{ + 0.08106146679532726, + 0.04134069595540929, + 0.02767792568499834, + 0.02079067210376509, + 0.01664469118982119, + 0.01387612882307075, + 0.01189670994589177, + 0.01041126526197209, + 0.009255462182712733, + 0.008330563433362871 + }}; + }; + template + SPROUT_CONSTEXPR sprout::array sprout::random::detail::binomial_table::table; + } // namespace detail + // + // binomial_distribution + // + template + class binomial_distribution { + public: + typedef RealType input_type; + typedef IntType result_type; + private: + struct btrd_type { + public: + RealType r; + RealType nr; + RealType npq; + RealType b; + RealType a; + RealType c; + RealType alpha; + RealType v_r; + RealType u_rv_r; + }; + private: + static SPROUT_CONSTEXPR IntType arg_check(IntType t_arg, RealType p_arg) { + return t_arg >= IntType(0) && RealType(0) <= p_arg && p_arg <= RealType(1) + ? t_arg + : throw "assert(t_arg >= IntType(0) && RealType(0) <= p_arg && p_arg <= RealType(1))" + ; + } + public: + // + // param_type + // + class param_type { + public: + typedef binomial_distribution distribution_type; + private: + IntType t_; + RealType p_; + public: + SPROUT_CONSTEXPR param_type() + : t_(1) + , p_(0.5) + {} + SPROUT_CONSTEXPR explicit param_type(IntType t_arg, RealType p_arg = RealType(0.5)) + : t_(arg_check(t_arg, p_arg)) + , p_(p_arg) + {} + SPROUT_CONSTEXPR IntType t() const { + return t_; + } + SPROUT_CONSTEXPR RealType p() const { + return p_; + } + template + friend std::basic_ostream& operator>>( + std::basic_istream& lhs, + param_type const& rhs + ) + { + return lhs >> rhs.t_ >> std::ws >> rhs.p_; + } + template + friend std::basic_ostream& operator<<( + std::basic_ostream& lhs, + param_type const& rhs + ) + { + return lhs << rhs.t_ << " " << rhs.p_; + } + SPROUT_CONSTEXPR friend bool operator==(param_type const& lhs, param_type const& rhs) { + return lhs.t_ == rhs.t_ && lhs.p_ == rhs.p_; + } + SPROUT_CONSTEXPR friend bool operator!=(param_type const& lhs, param_type const& rhs) { + return !(lhs == rhs); + } + }; + private: + static SPROUT_CONSTEXPR IntType init_t(IntType t) { + return t; + } + static SPROUT_CONSTEXPR RealType init_p(RealType p) { + return RealType(0.5) < p ? 1 - p : p; + } + static SPROUT_CONSTEXPR IntType init_m(IntType t, RealType p) { + return static_cast((init_t(t) + 1) * init_p(p)); + } + static SPROUT_CONSTEXPR btrd_type init_btrd_6(IntType t, RealType p, RealType r, RealType nr, RealType npq, RealType b, RealType a, RealType c, RealType alpha, RealType v_r) { + return btrd_type { + r, + nr, + npq, + b, + a, + c, + alpha, + v_r, + RealType(0.86) * v_r + }; + } + static SPROUT_CONSTEXPR btrd_type init_btrd_5(IntType t, RealType p, RealType r, RealType nr, RealType npq, RealType sqrt_npq, RealType b) { + return init_btrd_6(t, p, r, nr, npq, b, RealType(-0.0873) + RealType(0.0248) * b + RealType(0.01) * p, t * p + RealType(0.5), (RealType(2.83) + RealType(5.1) / b) * sqrt_npq, RealType(0.92) - RealType(4.2) / b); + } + static SPROUT_CONSTEXPR btrd_type init_btrd_4(IntType t, RealType p, RealType r, RealType nr, RealType npq, RealType sqrt_npq) { + return init_btrd_5(t, p, r, nr, npq, sqrt_npq, RealType(1.15) + RealType(2.53) * sqrt_npq); + } + static SPROUT_CONSTEXPR btrd_type init_btrd_3(IntType t, RealType p, RealType r, RealType nr, RealType npq) { + using std::sqrt; + return init_btrd_4(t, p, r, nr, npq, sqrt(npq)); + } + static SPROUT_CONSTEXPR btrd_type init_btrd_2(IntType t, RealType p, RealType r) { + return init_btrd_3(t, p, r, (t + 1) * r, t * p * (1 - p)); + } + static SPROUT_CONSTEXPR btrd_type init_btrd_1(IntType t, RealType p) { + return init_btrd_2(t, p, p / (1 - p)); + } + static SPROUT_CONSTEXPR btrd_type init_btrd(IntType t, RealType p) { + return init_btrd_1(init_t(t), init_p(p)); + } + static SPROUT_CONSTEXPR RealType init_q_n(IntType t, RealType p) { + using std::pow; + return pow(1 - init_p(p), static_cast(init_t(t))); + } + static SPROUT_CONSTEXPR bool init_use_inversion(IntType t, RealType p) { + return init_m(t, p) < 11; + } + static SPROUT_CONSTEXPR RealType fc_1(RealType ikp1) { + return (RealType(1) / 12 - (RealType(1) / 360 - (RealType(1) / 1260) * (ikp1 * ikp1)) * (ikp1 * ikp1)) * ikp1; + } + static SPROUT_CONSTEXPR RealType fc(IntType k) { + return k < 10 + ? sprout::random::detail::binomial_table::table[k] + : fc_1(RealType(1) / (k + 1)) + ; + } + private: + IntType t_; + RealType p_; + IntType m_; + btrd_type btrd_; + RealType q_n_; + private: + SPROUT_CONSTEXPR bool use_inversion() const { + return m_ < 11; + } + template + SPROUT_CONSTEXPR sprout::random::random_result invert_4(Engine const& eng, RealType u, RealType q, RealType s, RealType a, RealType r, IntType x = 0) const { + return u > r + ? invert_4(eng, u - r, q, s, a, ((a / (x + 1)) - s) * r, x + 1) + : sprout::random::random_result(x, eng, *this) + ; + } + template + SPROUT_CONSTEXPR sprout::random::random_result invert_3(IntType t, Engine const& eng, RealType u, RealType q, RealType s) const { + return invert_4(eng, u, q, s, (t + 1) * s, q_n_); + } + template + SPROUT_CONSTEXPR sprout::random::random_result invert_2(IntType t, RealType p, Engine const& eng, RealType u, RealType q) const { + return invert_3(t, eng, u, q, p / q); + } + template + SPROUT_CONSTEXPR sprout::random::random_result invert_1(IntType t, RealType p, Engine const& eng, RealType u) const { + return invert_2(t, p, eng, u, 1 - p); + } + template + SPROUT_CONSTEXPR sprout::random::random_result invert_0(IntType t, RealType p, Random const& rnd) const { + return invert_1(t, p, rnd.engine(), rnd.result()); + } + template + SPROUT_CONSTEXPR sprout::random::random_result invert(IntType t, RealType p, Engine const& eng) const { + return invert_0(t, p, sprout::random::uniform_01()(eng)); + } + template + SPROUT_CONSTEXPR sprout::random::random_result invert2_0(IntType t, sprout::random::random_result const& rnd) const { + return sprout::random::random_result(t - rnd.result(), rnd.engine(), rnd.distribution()); + } + template + SPROUT_CONSTEXPR sprout::random::random_result invert2(IntType t, RealType p, Engine const& eng) const { + return invert2_0(t, invert(t, p, eng)); + } + template + SPROUT_CONSTEXPR sprout::random::random_result generate_10(Engine const& eng, RealType v, IntType k, IntType nm, RealType h, IntType nk) const { + using std::log; + return v <= h + (t_ + 1) * log(static_cast(nm) / nk) + (k + RealType(0.5)) * log(nk * btrd_.r / (k + 1))- fc(k)- fc(t_ - k) + ? sprout::random::random_result(k, eng, *this) + : generate(eng) + ; + } + template + SPROUT_CONSTEXPR sprout::random::random_result generate_9(Engine const& eng, RealType v, IntType k, IntType nm) const { + using std::log; + return generate_10(eng, v, k, nm, (m_ + RealType(0.5)) * log((m_ + 1) / (btrd_.r * nm)) + fc(m_) + fc(t_ - m_), t_ - k + 1); + } + template + SPROUT_CONSTEXPR sprout::random::random_result generate_8(Engine const& eng, RealType v, IntType k, RealType rho, RealType t) const { + return v < t - rho ? sprout::random::random_result(k, eng, *this) + : v > t + rho ? generate(eng) + : generate_9(eng, v, k, t_ - m_ + 1) + ; + } + template + SPROUT_CONSTEXPR sprout::random::random_result generate_7_3(Engine const& eng, RealType v, IntType k, RealType f) const { + return v <= f + ? sprout::random::random_result(k, eng, *this) + : generate(eng) + ; + } + template + SPROUT_CONSTEXPR sprout::random::random_result generate_7_2(Engine const& eng, RealType v, IntType k, RealType f, IntType i) const { + return i != k + ? generate_7_2(eng, v * (btrd_.nr / (i + 1) - btrd_.r), k, f, i + 1) + : generate_7_3(eng, v, k, f) + ; + } + template + SPROUT_CONSTEXPR sprout::random::random_result generate_7_1(Engine const& eng, RealType v, IntType k, RealType f, IntType i) const { + return i != k + ? generate_7_1(eng, v, k, f * (btrd_.nr / (i + 1) - btrd_.r), i + 1) + : generate_7_3(eng, v, k, f) + ; + } + template + SPROUT_CONSTEXPR sprout::random::random_result generate_7(Engine const& eng, RealType v, IntType k, RealType f = RealType(1)) const { + return m_ < k ? generate_7_1(eng, v, k, f * (btrd_.nr / (m_ + 1) - btrd_.r), m_ + 1) + : m_ > k ? generate_7_2(eng, v * (btrd_.nr / (k + 1) - btrd_.r), k, f, k + 1) + : generate_7_3(eng, v, k, f) + ; + } + template + SPROUT_CONSTEXPR sprout::random::random_result generate_6(Engine const& eng, RealType v, IntType k, RealType km) const { + using std::log; + return km <= 15 + ? generate_7(eng, v, k) + : generate_8(eng, log(v), k, (km / btrd_.npq) * (((km / RealType(3.0) + RealType(0.625)) * km + RealType(1.0) / 6) / btrd_.npq + RealType(0.5)), -km * km / (2 * btrd_.npq)) + ; + } + template + SPROUT_CONSTEXPR sprout::random::random_result generate_5(Engine const& eng, RealType v, RealType u, RealType us, IntType k) const { + using std::abs; + return k < 0 || k > t_ + ? generate(eng) + : generate_6(eng, v * btrd_.alpha / (btrd_.a / (us * us) + btrd_.b), k, abs(k - m_)) + ; + } + template + SPROUT_CONSTEXPR sprout::random::random_result generate_4(Engine const& eng, RealType v, RealType u, RealType us) const { + using std::floor; + return generate_5(eng, v, u, us, static_cast(floor((2 * btrd_.a / us + btrd_.b) * u + btrd_.c))); + } + template + SPROUT_CONSTEXPR sprout::random::random_result generate_3(Engine const& eng, RealType v, RealType u) const { + using std::abs; + return generate_4(eng, v, u, 0.5 - abs(u)); + } + template + SPROUT_CONSTEXPR sprout::random::random_result generate_2(Random const& rnd, RealType v) const { + return v >= btrd_.v_r + ? generate_3( + rnd.engine(), + v, + rnd.result() - RealType(0.5) + ) + : generate_3( + rnd.engine(), + rnd.result() * btrd_.v_r, + ((v / btrd_.v_r - RealType(0.93)) < 0 ? RealType(-0.5) : RealType(0.5)) - (v / btrd_.v_r - RealType(0.93)) + ) + ; + } + template + SPROUT_CONSTEXPR sprout::random::random_result generate_1_1(Engine const& eng, RealType u) const { + using std::floor; + using std::abs; + return sprout::random::random_result( + static_cast(floor((2 * btrd_.a / (RealType(0.5) - abs(u)) + btrd_.b) * u + btrd_.c)), + eng, + *this + ); + } + template + SPROUT_CONSTEXPR sprout::random::random_result generate_1(Engine const& eng, RealType v) const { + return v <= btrd_.u_rv_r + ? generate_1_1(eng, v / btrd_.v_r - RealType(0.43)) + : generate_2(sprout::random::uniform_01()(eng), v) + ; + } + template + SPROUT_CONSTEXPR sprout::random::random_result generate_0(Random const& rnd) const { + return generate_1(rnd.engine(), rnd.result()); + } + template + SPROUT_CONSTEXPR sprout::random::random_result generate(Engine const& eng) const { + return generate_0(sprout::random::uniform_01()(eng)); + } + template + SPROUT_CONSTEXPR sprout::random::random_result generate2_0(IntType t, sprout::random::random_result const& rnd) const { + return sprout::random::random_result(t - rnd.result(), rnd.engine(), rnd.distribution()); + } + template + SPROUT_CONSTEXPR sprout::random::random_result generate2(IntType t, Engine const& eng) const { + return generate2_0(t, generate(eng)); + } + void init() { + m_ = init_m(t_, p_); + if (use_inversion()) { + q_n_ = init_q_n(t_, p_); + } else { + btrd_ = init_btrd(t_, p_); + } + } + public: + SPROUT_CONSTEXPR binomial_distribution() + : t_(1) + , p_(RealType(0.5)) + , m_(init_m(1, RealType(0.5))) + , btrd_(!init_use_inversion(1, RealType(0.5)) ? init_btrd(1, RealType(0.5)) : btrd_type()) + , q_n_(init_use_inversion(1, RealType(0.5)) ? init_q_n(1, RealType(0.5)) : RealType()) + {} + SPROUT_CONSTEXPR explicit binomial_distribution(IntType t_arg, RealType p_arg = RealType(0.5)) + : t_(arg_check(t_arg, p_arg)) + , p_(p_arg) + , m_(init_m(t_arg, p_arg)) + , btrd_(!init_use_inversion(t_arg, p_arg) ? init_btrd(t_arg, p_arg) : btrd_type()) + , q_n_(init_use_inversion(t_arg, p_arg) ? init_q_n(t_arg, p_arg) : RealType()) + {} + SPROUT_CONSTEXPR explicit binomial_distribution(param_type const& parm) + : t_(parm.t()) + , p_(parm.p()) + , m_(init_m(parm.t(), parm.p())) + , btrd_(!init_use_inversion(parm.t(), parm.p()) ? init_btrd(parm.t(), parm.p()) : btrd_type()) + , q_n_(init_use_inversion(parm.t(), parm.p()) ? init_q_n(parm.t(), parm.p()) : RealType()) + {} + SPROUT_CONSTEXPR result_type t() const { + return t_; + } + SPROUT_CONSTEXPR result_type p() const { + return p_; + } + SPROUT_CONSTEXPR result_type min() const { + return 0; + } + SPROUT_CONSTEXPR result_type max() const { + return t_; + } + SPROUT_CONSTEXPR param_type param() const { + return param_type(t_, p_); + } + void param(param_type const& parm) { + t_ = parm.a(); + p_ = parm.b(); + init(); + } + template + SPROUT_CONSTEXPR sprout::random::random_result operator()(Engine const& eng) const { + return use_inversion() ? RealType(0.5) < p_ + ? invert2(t_, 1 - p_, eng) + : invert(t_, p_, eng) + : RealType(0.5) < p_ ? generate2(t_, eng) + : generate(eng) + ; + } + template + friend std::basic_ostream& operator>>( + std::basic_istream& lhs, + binomial_distribution const& rhs + ) + { + param_type parm; + return lhs >> parm; + param(parm); + return lhs; + } + template + friend std::basic_ostream& operator<<( + std::basic_ostream& lhs, + binomial_distribution const& rhs + ) + { + return lhs << param(); + } + SPROUT_CONSTEXPR friend bool operator==(binomial_distribution const& lhs, binomial_distribution const& rhs) { + return lhs.param() == rhs.param(); + } + SPROUT_CONSTEXPR friend bool operator!=(binomial_distribution const& lhs, binomial_distribution const& rhs) { + return !(lhs == rhs); + } + }; + } // namespace random + + using sprout::random::binomial_distribution; +} // namespace sprout + +#endif // #ifndef SPROUT_RANDOM_BINOMIAL_DISTRIBUTION_HPP +