diff --git a/EigenRand/Dists/Basic.h b/EigenRand/Dists/Basic.h index 1021823..b45af31 100644 --- a/EigenRand/Dists/Basic.h +++ b/EigenRand/Dists/Basic.h @@ -333,6 +333,52 @@ namespace Eigen }; + /** + * @brief Generator of Bernoulli distribution + * + * @tparam _Scalar + */ + template + class BernoulliGen : public GenBase, _Scalar> + { + uint32_t p; + public: + using Scalar = _Scalar; + + BernoulliGen(double _p = 0.5) + { + eigen_assert(0 <= _p && _p <= 1 ); + p = (uint32_t)(_p * 0x80000000); + } + + BernoulliGen(const BernoulliGen&) = default; + BernoulliGen(BernoulliGen&&) = default; + + BernoulliGen& operator=(const BernoulliGen&) = default; + BernoulliGen& operator=(BernoulliGen&&) = default; + + template + EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) + { + using namespace Eigen::internal; + return (((uint32_t)pfirst(std::forward(rng)()) & 0x7FFFFFFF) < p) ? 1 : 0; + } + + template + EIGEN_STRONG_INLINE const Packet packetOp(Rng&& rng) + { + using namespace Eigen::internal; + using IPacket = decltype(reinterpret_to_int(std::declval())); + using RUtils = RawbitsMaker; + auto one = pset1(1); + auto zero = pset1(0); + auto r = RUtils{}.rawbits(std::forward(rng)); + r = pand(r, pset1(0x7FFFFFFF)); + return pblendv(pcmplt(r, pset1(p)), one, zero); + } + }; + + template using RandBitsType = CwiseNullaryOp, typename Derived::Scalar, Urng, true>, const Derived>; @@ -464,6 +510,48 @@ namespace Eigen o.rows(), o.cols(), { std::forward(urng) } }; } + + template + using BernoulliType = CwiseNullaryOp, typename Derived::Scalar, Urng, true>, const Derived>; + + /** + * @brief generates 1 with probability `p` and 0 with probability `1 - p` + * + * @tparam Derived + * @tparam Urng + * @param rows the number of rows being generated + * @param cols the number of columns being generated + * @param urng c++11-style random number generator + * @param p a probability of generating 1 + * @return a random matrix expression with a shape (`rows`, `cols`) + */ + template + inline const BernoulliType + bernoulli(Index rows, Index cols, Urng&& urng, double p = 0.5) + { + return { + rows, cols, { std::forward(urng), BernoulliGen{ p } } + }; + } + + /** + * @brief generates 1 with probability `p` and 0 with probability `1 - p` + * + * @tparam Derived + * @tparam Urng + * @param o an instance of any type of Eigen::DenseBase + * @param urng c++11-style random number generator + * @param p a probability of generating 1 + * @return a random matrix expression of the same shape as `o` + */ + template + inline const BernoulliType + bernoulli(Derived& o, Urng&& urng, double p = 0.5) + { + return { + o.rows(), o.cols(), { std::forward(urng), BernoulliGen{ p } } + }; + } } } diff --git a/EigenRand/Dists/Discrete.h b/EigenRand/Dists/Discrete.h index 2e242d6..667054f 100644 --- a/EigenRand/Dists/Discrete.h +++ b/EigenRand/Dists/Discrete.h @@ -228,7 +228,7 @@ namespace Eigen * * @param _min, _max the range of integers being generated */ - UniformIntGen(_Scalar _min, _Scalar _max) + UniformIntGen(_Scalar _min = 0, _Scalar _max = 0) : pmin{ _min }, pdiff{ (size_t)(_max - _min) } { if ((pdiff + 1) > pdiff) @@ -245,6 +245,9 @@ namespace Eigen UniformIntGen(const UniformIntGen&) = default; UniformIntGen(UniformIntGen&&) = default; + UniformIntGen& operator=(const UniformIntGen&) = default; + UniformIntGen& operator=(UniformIntGen&&) = default; + template EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) { @@ -389,6 +392,17 @@ namespace Eigen { } + DiscreteGen() + : DiscreteGen({ 1 }) + { + } + + DiscreteGen(const DiscreteGen&) = default; + DiscreteGen(DiscreteGen&&) = default; + + DiscreteGen& operator=(const DiscreteGen&) = default; + DiscreteGen& operator=(DiscreteGen&&) = default; + template EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) { @@ -513,6 +527,17 @@ namespace Eigen { } + DiscreteGen() + : DiscreteGen({ 1 }) + { + } + + DiscreteGen(const DiscreteGen&) = default; + DiscreteGen(DiscreteGen&&) = default; + + DiscreteGen& operator=(const DiscreteGen&) = default; + DiscreteGen& operator=(DiscreteGen&&) = default; + template EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) { @@ -619,6 +644,17 @@ namespace Eigen { } + DiscreteGen() + : DiscreteGen({ 1 }) + { + } + + DiscreteGen(const DiscreteGen&) = default; + DiscreteGen(DiscreteGen&&) = default; + + DiscreteGen& operator=(const DiscreteGen&) = default; + DiscreteGen& operator=(DiscreteGen&&) = default; + template EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) { @@ -708,7 +744,7 @@ namespace Eigen * * @param _mean mean of the distribution */ - PoissonGen(double _mean) + PoissonGen(double _mean = 1) : mean{ _mean }, ne_mean{ std::exp(-_mean) } { sqrt_tmean = std::sqrt(2 * mean); @@ -716,6 +752,12 @@ namespace Eigen g1 = mean * log_mean - std::lgamma(mean + 1); } + PoissonGen(const PoissonGen&) = default; + PoissonGen(PoissonGen&&) = default; + + PoissonGen& operator=(const PoissonGen&) = default; + PoissonGen& operator=(PoissonGen&&) = default; + template EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) { @@ -835,6 +877,12 @@ namespace Eigen } } + BinomialGen(const BinomialGen&) = default; + BinomialGen(BinomialGen&&) = default; + + BinomialGen& operator=(const BinomialGen&) = default; + BinomialGen& operator=(BinomialGen&&) = default; + template EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) { @@ -959,11 +1007,17 @@ namespace Eigen * * @param _p probability of a trial generating true */ - GeometricGen(double _p) + GeometricGen(double _p = 0.5) : p{ _p }, rlog_q{ 1 / std::log(1 - p) } { } + GeometricGen(const GeometricGen&) = default; + GeometricGen(GeometricGen&&) = default; + + GeometricGen& operator=(const GeometricGen&) = default; + GeometricGen& operator=(GeometricGen&&) = default; + template EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) { diff --git a/EigenRand/Dists/NormalExp.h b/EigenRand/Dists/NormalExp.h index 6663c51..2bc8218 100644 --- a/EigenRand/Dists/NormalExp.h +++ b/EigenRand/Dists/NormalExp.h @@ -110,6 +110,9 @@ namespace Eigen NormalGen(const NormalGen&) = default; NormalGen(NormalGen&&) = default; + NormalGen& operator=(const NormalGen&) = default; + NormalGen& operator=(NormalGen&&) = default; + template EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) { @@ -156,6 +159,9 @@ namespace Eigen LognormalGen(const LognormalGen&) = default; LognormalGen(LognormalGen&&) = default; + LognormalGen& operator=(const LognormalGen&) = default; + LognormalGen& operator=(LognormalGen&&) = default; + template EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) { @@ -199,6 +205,9 @@ namespace Eigen StudentTGen(const StudentTGen&) = default; StudentTGen(StudentTGen&&) = default; + StudentTGen& operator=(const StudentTGen&) = default; + StudentTGen& operator=(StudentTGen&&) = default; + template EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) { @@ -267,6 +276,9 @@ namespace Eigen ExponentialGen(const ExponentialGen&) = default; ExponentialGen(ExponentialGen&&) = default; + ExponentialGen& operator=(const ExponentialGen&) = default; + ExponentialGen& operator=(ExponentialGen&&) = default; + template EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) { @@ -320,6 +332,9 @@ namespace Eigen GammaGen(const GammaGen&) = default; GammaGen(GammaGen&&) = default; + GammaGen& operator=(const GammaGen&) = default; + GammaGen& operator=(GammaGen&&) = default; + template EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) { @@ -487,6 +502,9 @@ namespace Eigen WeibullGen(const WeibullGen&) = default; WeibullGen(WeibullGen&&) = default; + WeibullGen& operator=(const WeibullGen&) = default; + WeibullGen& operator=(WeibullGen&&) = default; + template EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) { @@ -533,6 +551,9 @@ namespace Eigen ExtremeValueGen(const ExtremeValueGen&) = default; ExtremeValueGen(ExtremeValueGen&&) = default; + ExtremeValueGen& operator=(const ExtremeValueGen&) = default; + ExtremeValueGen& operator=(ExtremeValueGen&&) = default; + template EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) { @@ -577,6 +598,9 @@ namespace Eigen ChiSquaredGen(const ChiSquaredGen&) = default; ChiSquaredGen(ChiSquaredGen&&) = default; + ChiSquaredGen& operator=(const ChiSquaredGen&) = default; + ChiSquaredGen& operator=(ChiSquaredGen&&) = default; + template EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) { @@ -619,6 +643,9 @@ namespace Eigen CauchyGen(const CauchyGen&) = default; CauchyGen(CauchyGen&&) = default; + CauchyGen& operator=(const CauchyGen&) = default; + CauchyGen& operator=(CauchyGen&&) = default; + template EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) { @@ -674,6 +701,9 @@ namespace Eigen BetaGen(const BetaGen&) = default; BetaGen(BetaGen&&) = default; + BetaGen& operator=(const BetaGen&) = default; + BetaGen& operator=(BetaGen&&) = default; + template EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) { @@ -752,6 +782,9 @@ namespace Eigen FisherFGen(const FisherFGen&) = default; FisherFGen(FisherFGen&&) = default; + FisherFGen& operator=(const FisherFGen&) = default; + FisherFGen& operator=(FisherFGen&&) = default; + template EIGEN_STRONG_INLINE const _Scalar operator() (Rng&& rng) { diff --git a/EigenRand/Macro.h b/EigenRand/Macro.h index f1f6864..356a8e2 100644 --- a/EigenRand/Macro.h +++ b/EigenRand/Macro.h @@ -14,7 +14,7 @@ #define EIGENRAND_WORLD_VERSION 0 #define EIGENRAND_MAJOR_VERSION 3 -#define EIGENRAND_MINOR_VERSION 0 +#define EIGENRAND_MINOR_VERSION 2 #if EIGEN_VERSION_AT_LEAST(3,3,7) #else diff --git a/EigenRand/MorePacketMath.h b/EigenRand/MorePacketMath.h index b43d511..49403fc 100644 --- a/EigenRand/MorePacketMath.h +++ b/EigenRand/MorePacketMath.h @@ -195,8 +195,8 @@ namespace Eigen template EIGEN_STRONG_INLINE Packet pcmple(const Packet& a, const Packet& b); - template - EIGEN_STRONG_INLINE Packet pblendv(const Packet& ifPacket, const Packet& thenPacket, const Packet& elsePacket); + template + EIGEN_STRONG_INLINE Packet pblendv(const PacketIf& ifPacket, const Packet& thenPacket, const Packet& elsePacket); template EIGEN_STRONG_INLINE Packet pgather(const int* addr, const Packet& index); @@ -527,7 +527,14 @@ namespace Eigen template<> EIGEN_STRONG_INLINE Packet8i pcmplt(const Packet8i& a, const Packet8i& b) { +#ifdef EIGEN_VECTORIZE_AVX2 return _mm256_cmpgt_epi32(b, a); +#else + Packet4i a1, a2, b1, b2; + split_two(a, a1, a2); + split_two(b, b1, b2); + return combine_two((Packet4i)_mm_cmpgt_epi32(b1, a1), (Packet4i)_mm_cmpgt_epi32(b2, a2)); +#endif } template<> @@ -560,6 +567,12 @@ namespace Eigen return _mm256_blendv_ps(elsePacket, thenPacket, ifPacket); } + template<> + EIGEN_STRONG_INLINE Packet8f pblendv(const Packet8i& ifPacket, const Packet8f& thenPacket, const Packet8f& elsePacket) + { + return pblendv(_mm256_castsi256_ps(ifPacket), thenPacket, elsePacket); + } + template<> EIGEN_STRONG_INLINE Packet8i pblendv(const Packet8i& ifPacket, const Packet8i& thenPacket, const Packet8i& elsePacket) { @@ -576,6 +589,12 @@ namespace Eigen return _mm256_blendv_pd(elsePacket, thenPacket, ifPacket); } + template<> + EIGEN_STRONG_INLINE Packet4d pblendv(const Packet8i& ifPacket, const Packet4d& thenPacket, const Packet4d& elsePacket) + { + return pblendv(_mm256_castsi256_pd(ifPacket), thenPacket, elsePacket); + } + template<> EIGEN_STRONG_INLINE Packet8i pgather(const int* addr, const Packet8i& index) { @@ -842,6 +861,12 @@ namespace Eigen #endif } + template<> + EIGEN_STRONG_INLINE Packet4f pblendv(const Packet4i& ifPacket, const Packet4f& thenPacket, const Packet4f& elsePacket) + { + return pblendv(_mm_castsi128_ps(ifPacket), thenPacket, elsePacket); + } + template<> EIGEN_STRONG_INLINE Packet4i pblendv(const Packet4i& ifPacket, const Packet4i& thenPacket, const Packet4i& elsePacket) { @@ -862,6 +887,13 @@ namespace Eigen #endif } + + template<> + EIGEN_STRONG_INLINE Packet2d pblendv(const Packet4i& ifPacket, const Packet2d& thenPacket, const Packet2d& elsePacket) + { + return pblendv(_mm_castsi128_pd(ifPacket), thenPacket, elsePacket); + } + template<> EIGEN_STRONG_INLINE Packet4i pgather(const int* addr, const Packet4i& index) { diff --git a/EigenRand/MvDists/Multinomial.h b/EigenRand/MvDists/Multinomial.h index e31a089..17de6ec 100644 --- a/EigenRand/MvDists/Multinomial.h +++ b/EigenRand/MvDists/Multinomial.h @@ -52,6 +52,9 @@ namespace Eigen MultinomialGen(const MultinomialGen&) = default; MultinomialGen(MultinomialGen&&) = default; + MultinomialGen& operator=(const MultinomialGen&) = default; + MultinomialGen& operator=(MultinomialGen&&) = default; + Index dims() const { return probs.rows(); } template diff --git a/EigenRand/MvDists/MvNormal.h b/EigenRand/MvDists/MvNormal.h index adf958c..b26986d 100644 --- a/EigenRand/MvDists/MvNormal.h +++ b/EigenRand/MvDists/MvNormal.h @@ -93,6 +93,9 @@ namespace Eigen MvNormalGen(const MvNormalGen&) = default; MvNormalGen(MvNormalGen&&) = default; + MvNormalGen& operator=(const MvNormalGen&) = default; + MvNormalGen& operator=(MvNormalGen&&) = default; + Index dims() const { return mean.rows(); } template @@ -119,7 +122,7 @@ namespace Eigen * @param cov covariance matrix (should be positive semi-definite) */ template - inline auto makeMvNormGen(const MatrixBase& mean, const MatrixBase& cov) + inline auto makeMvNormalGen(const MatrixBase& mean, const MatrixBase& cov) -> MvNormalGen::Scalar, MatrixBase::RowsAtCompileTime> { static_assert( @@ -143,7 +146,7 @@ namespace Eigen * @param lt lower triangular matrix of decomposed covariance */ template - inline auto makeMvNormGenFromLt(const MatrixBase& mean, const MatrixBase& lt) + inline auto makeMvNormalGenFromLt(const MatrixBase& mean, const MatrixBase& lt) -> MvNormalGen::Scalar, MatrixBase::RowsAtCompileTime> { static_assert( @@ -209,6 +212,9 @@ namespace Eigen WishartGen(const WishartGen&) = default; WishartGen(WishartGen&&) = default; + WishartGen& operator=(const WishartGen&) = default; + WishartGen& operator=(WishartGen&&) = default; + Index dims() const { return chol.rows(); } template @@ -366,6 +372,9 @@ namespace Eigen InvWishartGen(const InvWishartGen&) = default; InvWishartGen(InvWishartGen&&) = default; + InvWishartGen& operator=(const InvWishartGen&) = default; + InvWishartGen& operator=(InvWishartGen&&) = default; + Index dims() const { return chol.rows(); } template diff --git a/README.md b/README.md index 1e55a15..e584f9a 100644 --- a/README.md +++ b/README.md @@ -347,6 +347,10 @@ The results of EigenRand and C++ std appear to be equivalent within the margin o MIT License ## History + +### 0.3.2 (2021-03-26) +* A default constructor for `DiscreteGen` was added. + ### 0.3.1 (2020-11-15) * Compiling errors in the environment `EIGEN_COMP_MINGW && __GXX_ABI_VERSION < 1004` was fixed. diff --git a/test/benchmark.cpp b/test/benchmark.cpp index 3bd3bfc..3e56ac1 100644 --- a/test/benchmark.cpp +++ b/test/benchmark.cpp @@ -190,6 +190,39 @@ std::map test_eigenrand(size_t size, const std::string& suf auto scope = bh.measure("discrete/int(s=250)" + suffix, xi); xi = Eigen::Rand::discreteLike(xi, urng, ws.begin(), ws.end()); } + + { + auto scope = bh.measure("bernoulli/int" + suffix, xi); + xi = Eigen::Rand::bernoulli(xi, urng, 1. / 3); + } + + { + auto scope = bh.measure("bernoulli/int/gen" + suffix, xi); + Eigen::Rand::BernoulliGen gen{ 1. / 3}; + xi = gen.generateLike(xi, urng); + } + + { + auto scope = bh.measure("bernoulli" + suffix, x); + x = Eigen::Rand::bernoulli(x, urng, 1. / 3); + } + + { + auto scope = bh.measure("bernoulli/gen" + suffix, x); + Eigen::Rand::BernoulliGen gen{ 1. / 3 }; + x = gen.generateLike(x, urng); + } + + { + auto scope = bh.measure("bernoulli/double" + suffix, xd); + xd = Eigen::Rand::bernoulli(xd, urng, 1. / 3); + } + + { + auto scope = bh.measure("bernoulli/double/gen" + suffix, x); + Eigen::Rand::BernoulliGen gen{ 1. / 3 }; + xd = gen.generateLike(xd, urng); + } { auto scope = bh.measure("balanced" + suffix, x); diff --git a/test/benchmark_mv.cpp b/test/benchmark_mv.cpp index fa78ad3..c10f9b0 100644 --- a/test/benchmark_mv.cpp +++ b/test/benchmark_mv.cpp @@ -133,7 +133,7 @@ std::map test_eigenrand(size_t size, const std::string& suf 1, 2, 0, 0, 0, 0, 3, 1, 0, 0, 1, 2; - test_single_and_batch(size, bh, "MvNormal/4/float" + suffix, Eigen::Rand::makeMvNormGen(mean, cov), urng); + test_single_and_batch(size, bh, "MvNormal/4/float" + suffix, Eigen::Rand::makeMvNormalGen(mean, cov), urng); } { @@ -142,7 +142,7 @@ std::map test_eigenrand(size_t size, const std::string& suf cov = Eigen::Rand::normalLike(cov, urng) * 0.1f; cov = cov * cov.transpose(); cov += Eigen::MatrixXf::Identity(100, 100) * 4; - test_single_and_batch(size, bh, "MvNormal/100/float" + suffix, Eigen::Rand::makeMvNormGen(mean, cov), urng); + test_single_and_batch(size, bh, "MvNormal/100/float" + suffix, Eigen::Rand::makeMvNormalGen(mean, cov), urng); } {