@@ -575,21 +575,21 @@ namespace grb {
575575 }
576576
577577 grb::Vector< QType, backend > h ( n );
578- grb::Vector< QType, backend > log_rand ( n );
578+ grb::Vector< QType, backend > rand ( n );
579579 grb::Vector< StateType, backend > delta ( n );
580580 grb::Vector< EnergyType, backend > dn ( n );
581581 grb::Vector< bool , backend > accept ( n );
582582 std::srand ( static_cast <unsigned >( seed ) );
583583 std::minstd_rand rng ( seed ); // minstd_rand or std::mt19937
584584
585585 grb::resize ( h, n );
586- grb::resize ( log_rand , n );
586+ grb::resize ( rand , n );
587587 grb::resize ( delta, n );
588588 grb::resize ( dn, n );
589589 grb::resize ( accept, n );
590590
591591 std::vector< grb::Vector< bool , backend > > masks ;
592- rc = rc ? rc : matrix_partition< descr >( masks, couplings, h, log_rand , seed );
592+ rc = rc ? rc : matrix_partition< descr >( masks, couplings, h, rand , seed );
593593 grb::clear (h);
594594 constexpr auto dense_descr = descr | grb::descriptors::dense;
595595
@@ -598,15 +598,19 @@ namespace grb {
598598 (const typeof (local_fields)&) local_fields,
599599 (const typeof (masks)&) masks,
600600 h,
601- log_rand ,
601+ rand ,
602602 delta,
603603 dn,
604604 accept,
605605 rng,
606606 (const typeof (ring)&) ring
607607 );
608608
609- const auto ising_sweep = [](
609+ #ifdef NDEBUG
610+ const auto ising_sweep = [](
611+ #else
612+ const auto ising_sweep = [&get_energy](
613+ #endif
610614 grb::Vector< StateType, backend > &state,
611615 const TempType &beta,
612616 typeof (sweep_data) &data
@@ -618,7 +622,7 @@ namespace grb {
618622 const auto &local_fields = std::get<1 >(data);
619623 const auto &masks = std::get<2 >(data);
620624 auto &h = std::get<3 >(data);
621- auto &log_rand = std::get<4 >(data);
625+ auto &rand = std::get<4 >(data);
622626 auto &delta = std::get<5 >(data);
623627 auto &dn = std::get<6 >(data);
624628 auto &accept = std::get<7 >(data);
@@ -628,49 +632,45 @@ namespace grb {
628632 const size_t n = grb::size ( state );
629633 EnergyType delta_energy = static_cast < EnergyType >(0.0 );
630634 grb::RC rc = grb::SUCCESS;
635+ (void ) n;
631636
632637 if ( !empty_local_fields) {
633638 rc = rc ? rc : grb::set< descr >( h, local_fields );
634639 }else {
635640 rc = rc ? rc : grb::set< descr >( h, static_cast < QType >( 0.0 ) );
636641 }
637642 rc = rc ? rc : grb::mxv< dense_descr >( h, couplings, state , ring );
638- std::uniform_real_distribution< QType > rand ( 0.0 , 1.0 );
639- for ( size_t j = 0 ; j < n ; ++j ){
640- const auto rnd = rand ( rng );
641- rc = rc ? rc : grb::setElement (log_rand, internal::log ( rnd ), j );
643+ std::uniform_real_distribution< QType > rand_gen ( 0.0 , 1.0 );
644+ for ( size_t i = 0 ; i < n; ++i ){
645+ grb::setElement ( rand, rand_gen ( rng ), i );
642646 }
647+
643648#ifndef NDEBUG
644649 const grb::Vector< StateType > old_state = state;
645650#endif
646651 rc = rc ? rc : grb::wait ();
647652 for (const auto &mask : masks ){
648- rc = rc ? rc : grb::clear ( accept );
649- rc = rc ? rc : grb::clear ( delta );
650- rc = rc ? rc : grb::clear ( dn );
651-
652653 // dn = (2*state_slice - 1) * h_slice
653654 rc = rc ? rc : grb::set< descr >( dn, mask, state );
654655 rc = rc ? rc : grb::foldl< descr >( dn, static_cast < EnergyType >( 2 ), ring.getMultiplicativeMonoid () );
655656 rc = rc ? rc : grb::foldl< descr >( dn, static_cast < EnergyType >( -1 ), ring.getAdditiveMonoid () );
656657 rc = rc ? rc : grb::foldl< descr >( dn, h, ring.getMultiplicativeMonoid () );
657658
658- // ( dn >= 0 ) | ( log_rand < beta * dn )
659+ // ( dn >= 0 ) | ( rand < beta * dn )
659660 rc = rc ? rc : grb::set< descr >( accept, mask );
660661 rc = rc ? rc : grb::wait (); // needed to avoid ERROR: Segmentation Fault with nonblocking backend
661662 rc = rc ? rc : grb::eWiseLambda< descr >(
662- [ &mask, &accept, &dn, &log_rand , beta ]( const size_t i ){
663+ [ &mask, &accept, &dn, &rand , beta ]( const size_t i ){
663664 if ( mask[i] ){
664- accept[i] = ( dn[i] >= 0 ) || ( log_rand [i] < beta * dn[i] );
665+ accept[i] = ( dn[i] >= 0 ) || ( internal::log ( rand [i] ) < beta * dn[i] );
665666 }
666- }, mask, log_rand , dn, accept );
667+ }, mask, rand , dn, accept );
667668
668669 // new_state = np.where(accept, 1 - old, old)
669670 rc = rc ? rc : grb::foldl< descr >( state, accept, static_cast < StateType >( -1 ), ring.getMultiplicativeMonoid () );
670671 rc = rc ? rc : grb::foldl< descr >( state, accept, static_cast < StateType >( 1 ), ring.getAdditiveMonoid () );
671672
672673 // delta = new - old ==> delta[accept] = 2*new_state[accept]-1
673- rc = rc ? rc : grb::clear ( delta );
674674 rc = rc ? rc : grb::set< descr >( delta, accept, state );
675675 rc = rc ? rc : grb::foldl< descr >( delta, accept, static_cast < StateType >( 2 ), ring.getMultiplicativeMonoid () );
676676 rc = rc ? rc : grb::foldl< descr >( delta, accept, static_cast < StateType >( -1 ), ring.getAdditiveMonoid () );
@@ -694,8 +694,8 @@ namespace grb {
694694
695695 EnergyType e1 = static_cast < EnergyType >( 0.0 ),
696696 e2 = static_cast < EnergyType >( 0.0 );
697- // get_energy(e1, old_state);
698- // get_energy(e2, new_state);
697+ get_energy (e1 , old_state);
698+ get_energy (e2 , new_state);
699699 const auto real_delta = e2 - e1 ;
700700 if ( s == 0 ){
701701 std::cerr << " \n\t Delta_energy: " << delta_energy;
@@ -704,7 +704,7 @@ namespace grb {
704704 std::cerr << std::endl;
705705 }
706706
707- // assert( ISCLOSE(real_delta, delta_energy ) );
707+ assert ( ISCLOSE (real_delta, delta_energy ) );
708708#endif
709709 return delta_energy;
710710 };
0 commit comments