"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "include/armadillo_bits/arma_rng_cxx11.hpp" between
armadillo-10.8.2.tar.xz and armadillo-11.0.0.tar.xz

About: Armadillo is a C++ linear algebra library (matrix maths) aiming towards a good balance between speed and ease of use.

arma_rng_cxx11.hpp  (armadillo-10.8.2.tar.xz):arma_rng_cxx11.hpp  (armadillo-11.0.0.tar.xz)
skipping to change at line 18 skipping to change at line 18
// You may obtain a copy of the License at // You may obtain a copy of the License at
// http://www.apache.org/licenses/LICENSE-2.0 // http://www.apache.org/licenses/LICENSE-2.0
// //
// Unless required by applicable law or agreed to in writing, software // Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS, // distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and // See the License for the specific language governing permissions and
// limitations under the License. // limitations under the License.
// ------------------------------------------------------------------------ // ------------------------------------------------------------------------
//! \addtogroup arma_rng_cxx11 //! \addtogroup arma_rng_cxx03
//! @{ //! @{
class arma_rng_cxx11 class arma_rng_cxx03
{ {
public: public:
typedef std::mt19937_64::result_type seed_type; typedef unsigned int seed_type;
inline void set_seed(const seed_type val); inline static void set_seed(const seed_type val);
arma_inline int randi_val(); arma_inline static int randi_val();
arma_inline double randu_val(); arma_inline static double randu_val();
arma_inline double randn_val(); inline static double randn_val();
template<typename eT> template<typename eT>
arma_inline void randn_dual_val(eT& out1, eT& out2); inline static void randn_dual_val(eT& out1, eT& out2);
template<typename eT> template<typename eT>
inline void randi_fill(eT* mem, const uword N, const int a, const int b); inline static void randi_fill(eT* mem, const uword N, const int a, const int b );
inline static int randi_max_val(); inline static int randi_max_val();
template<typename eT>
inline void randg_fill_simple(eT* mem, const uword N, const double a, const do
uble b);
template<typename eT>
inline void randg_fill(eT* mem, const uword N, const double a, const double b)
;
private:
arma_aligned std::mt19937_64 engine; // typedef for
std::mersenne_twister_engine with preset parameters
arma_aligned std::uniform_int_distribution<int> i_distr; // by default u
ses a=0, b=std::numeric_limits<int>::max()
arma_aligned std::uniform_real_distribution<double> u_distr; // by default u
ses [0,1) interval
arma_aligned std::normal_distribution<double> n_distr; // by default u
ses mean=0.0 and stddev=1.0
}; };
inline inline
void void
arma_rng_cxx11::set_seed(const arma_rng_cxx11::seed_type val) arma_rng_cxx03::set_seed(const arma_rng_cxx03::seed_type val)
{ {
engine.seed(val); std::srand(val);
i_distr.reset();
u_distr.reset();
n_distr.reset();
} }
arma_inline arma_inline
int int
arma_rng_cxx11::randi_val() arma_rng_cxx03::randi_val()
{ {
return i_distr(engine); #if (RAND_MAX == 32767)
} {
// NOTE: this is a better-than-nothing solution
// NOTE: see also arma_rng_cxx03::randi_max_val()
arma_inline u32 val1 = u32(std::rand());
double u32 val2 = u32(std::rand());
arma_rng_cxx11::randu_val()
{
return u_distr(engine);
}
arma_inline val1 <<= 15;
double
arma_rng_cxx11::randn_val() return (val1 | val2);
{ }
return n_distr(engine); #else
{
return std::rand();
}
#endif
} }
template<typename eT>
arma_inline arma_inline
void double
arma_rng_cxx11::randn_dual_val(eT& out1, eT& out2) arma_rng_cxx03::randu_val()
{ {
out1 = eT( n_distr(engine) ); return double( double(randi_val()) * ( double(1) / double(randi_max_val()) ) )
out2 = eT( n_distr(engine) ); ;
} }
template<typename eT>
inline inline
void double
arma_rng_cxx11::randi_fill(eT* mem, const uword N, const int a, const int b) arma_rng_cxx03::randn_val()
{ {
std::uniform_int_distribution<int> local_i_distr(a, b); // polar form of the Box-Muller transformation:
// http://en.wikipedia.org/wiki/Box-Muller_transformation
// http://en.wikipedia.org/wiki/Marsaglia_polar_method
for(uword i=0; i<N; ++i) double tmp1 = double(0);
double tmp2 = double(0);
double w = double(0);
do
{ {
mem[i] = eT(local_i_distr(engine)); tmp1 = double(2) * double(randi_val()) * (double(1) / double(randi_max_val()
)) - double(1);
tmp2 = double(2) * double(randi_val()) * (double(1) / double(randi_max_val()
)) - double(1);
w = tmp1*tmp1 + tmp2*tmp2;
} }
} while( w >= double(1) );
inline return double( tmp1 * std::sqrt( (double(-2) * std::log(w)) / w) );
int
arma_rng_cxx11::randi_max_val()
{
return std::numeric_limits<int>::max();
} }
template<typename eT> template<typename eT>
inline inline
void void
arma_rng_cxx11::randg_fill_simple(eT* mem, const uword N, const double a, const double b) arma_rng_cxx03::randn_dual_val(eT& out1, eT& out2)
{ {
std::gamma_distribution<double> g_distr(a,b); // make sure we are internally using at least floats
typedef typename promote_type<eT,float>::result eTp;
for(uword i=0; i<N; ++i) eTp tmp1 = eTp(0);
eTp tmp2 = eTp(0);
eTp w = eTp(0);
do
{ {
mem[i] = eT(g_distr(engine)); tmp1 = eTp(2) * eTp(randi_val()) * (eTp(1) / eTp(randi_max_val())) - eTp(1);
tmp2 = eTp(2) * eTp(randi_val()) * (eTp(1) / eTp(randi_max_val())) - eTp(1);
w = tmp1*tmp1 + tmp2*tmp2;
} }
while( w >= eTp(1) );
const eTp k = std::sqrt( (eTp(-2) * std::log(w)) / w);
out1 = eT(tmp1*k);
out2 = eT(tmp2*k);
} }
template<typename eT> template<typename eT>
inline inline
void void
arma_rng_cxx11::randg_fill(eT* mem, const uword N, const double a, const double b) arma_rng_cxx03::randi_fill(eT* mem, const uword N, const int a, const int b)
{ {
#if defined(ARMA_USE_OPENMP) if( (a == 0) && (b == RAND_MAX) )
{ {
if((N < 512) || omp_in_parallel()) { (*this).randg_fill_simple(mem, N, a, b for(uword i=0; i<N; ++i)
); return; }
typedef std::mt19937_64 motor_type;
typedef std::mt19937_64::result_type ovum_type;
typedef std::gamma_distribution<double> distr_type;
const uword n_threads = uword( mp_thread_limit::get() );
std::vector<motor_type> g_motor(n_threads);
std::vector<distr_type> g_distr(n_threads);
const distr_type g_distr_base(a,b);
for(uword t=0; t < n_threads; ++t)
{ {
motor_type& g_motor_t = g_motor[t]; mem[i] = eT(std::rand());
distr_type& g_distr_t = g_distr[t];
g_motor_t.seed( ovum_type(t) + ovum_type((*this).randi_val()) );
g_distr_t.param( g_distr_base.param() );
} }
}
else
{
const uword length = uword(b - a + 1);
const uword chunk_size = N / n_threads; const double scale = double(length) / double(randi_max_val());
#pragma omp parallel for schedule(static) num_threads(int(n_threads)) for(uword i=0; i<N; ++i)
for(uword t=0; t < n_threads; ++t)
{ {
const uword start = (t+0) * chunk_size; mem[i] = eT((std::min)( b, (int( double(randi_val()) * scale ) + a) ));
const uword endp1 = (t+1) * chunk_size;
motor_type& g_motor_t = g_motor[t];
distr_type& g_distr_t = g_distr[t];
for(uword i=start; i < endp1; ++i) { mem[i] = eT( g_distr_t(g_motor_t));
}
} }
motor_type& g_motor_0 = g_motor[0];
distr_type& g_distr_0 = g_distr[0];
for(uword i=(n_threads*chunk_size); i < N; ++i) { mem[i] = eT( g_distr_0(g_
motor_0)); }
} }
}
inline
int
arma_rng_cxx03::randi_max_val()
{
#if (RAND_MAX == 32767)
return ( (32767 << 15) + 32767);
#else #else
{ return RAND_MAX;
(*this).randg_fill_simple(mem, N, a, b);
}
#endif #endif
} }
//! @} //! @}
 End of changes. 40 change blocks. 
115 lines changed or deleted 87 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)