@@ -74,9 +74,8 @@ namespace grb {
7474 /*
7575 * Do a Parallel Tempering pass.
7676 * This means exchanging states at low temperature with states at higher temperature.
77- * To make the code simpler, this will be done by exchanging the temperatures instead.
7877 *
79- * TODO: Fix this documentation.
78+ * TODO: Complete this documentation.
8079 *
8180 * @param[in,out] states On input: initial states.
8281 * @param[in,out] energies The initial energy of each state.
@@ -101,18 +100,22 @@ namespace grb {
101100 pt (
102101 std::vector< grb::Vector< StateType, backend > > &states,
103102 grb::Vector< EnergyType, backend > &energies,
104- const grb::Vector< TempType, backend > &betas
103+ const grb::Vector< TempType, backend > &betas,
104+ const int seed = 42
105105 ){
106106
107107 const size_t n_replicas = states.size ();
108108 // const size_t s = spmd<>::pid();
109109 // const size_t nprocs = spmd<>::nprocs();
110+ std::srand ( seed );
110111 grb::RC rc = grb::SUCCESS;
112+ std::minstd_rand rng ( seed );
113+ std::exponential_distribution< EnergyType > rand ( 1.0 );
111114
112115 for ( size_t i = n_replicas - 1 ; i > 0 ; --i ){
113116 const EnergyType de = ( energies[ i ] - energies[ i-1 ]) * (betas[ i ] - betas[ i-1 ]);
114117
115- if ( de >= 0 || std:: rand () < RAND_MAX * internal::exp ( de ) ){
118+ if ( - rand ( rng ) < de ){
116119 std::swap ( states[i], states[i-1 ] );
117120 std::swap ( energies[i], energies[i-1 ] );
118121 }
@@ -136,7 +139,8 @@ namespace grb {
136139 pt (
137140 std::vector< grb::Vector< StateType, backend > > &states,
138141 grb::Vector< EnergyType, backend > &energies,
139- const grb::Vector< TempType, backend > &betas
142+ const grb::Vector< TempType, backend > &betas,
143+ const int seed = 42
140144 ){
141145 static_assert ( backend != grb::BSP1D );
142146 // static_assert( grb::_GRB_BACKEND == grb::BSP1D );
@@ -151,29 +155,30 @@ namespace grb {
151155 assert ( grb::size (energies) == n_replicas );
152156 assert ( grb::size (betas) == n_replicas );
153157#endif
158+ std::minstd_rand rng ( seed + s );
159+ std::exponential_distribution< EnergyType > rand ( 1.0 );
154160 struct data {
155161 EnergyType e;
156162 TempType b;
157- int r;
163+ EnergyType r;
158164 };
159165 grb::Vector< StateType, backend > s0 ( n );
160166 grb::Vector< StateType, backend > s1 ( n );
161167 grb::set ( s0, static_cast < StateType >( 0 ) );
162168 grb::set ( s1, static_cast < StateType >( 0 ) );
163169
164-
165170 struct data msg[ 2 ];
166171 rc = rc ? rc : grb::resize ( s0, n );
167172 rc = rc ? rc : grb::resize ( s1, n );
168173 if ( rc != grb::SUCCESS ) return rc;
169- int rand = std:: rand ();
174+ const auto myrand = - rand ( rng );
170175
171176 for ( size_t si = nprocs ; rc == grb::SUCCESS && si > 0 ; --si ){
172177 if ( si-1 == s ){
173178 for ( size_t i = n_replicas - 1 ; i > 0 ; --i ){
174179 const EnergyType de = ( energies[ i ] - energies[ i-1 ]) * (betas[ i ] - betas[ i-1 ]);
175180
176- if ( de >= 0 || std:: rand () < RAND_MAX * internal::exp ( de ) ){
181+ if ( - rand ( rng ) < de ){
177182 std::swap ( states[i], states[i-1 ] );
178183 std::swap ( energies[i], energies[i-1 ] );
179184 }
@@ -186,7 +191,7 @@ namespace grb {
186191 grb::set ( s0, states[ n_replicas - 1 ] );
187192 msg[ 0 ].e = energies[ n_replicas - 1 ];
188193 msg[ 0 ].b = betas[ n_replicas - 1 ];
189- msg[ 0 ].r = rand ;
194+ msg[ 0 ].r = myrand ;
190195 }
191196 if ( si == 1 ) continue ;
192197
@@ -216,9 +221,8 @@ namespace grb {
216221
217222 const EnergyType de = ( msg[ 1 ].e - msg[ 0 ].e ) * ( msg[ 1 ].b - msg[ 0 ].b );
218223
219- if ( rc == grb::SUCCESS && ( de >= 0 || msg[ 0 ].r < RAND_MAX * internal::exp ( de ) ) ){
224+ if ( rc == grb::SUCCESS && ( msg[ 0 ].r < de ) ){
220225 if ( si == s+1 ){
221-
222226 rc = rc ? rc : grb::set ( states[ n_replicas - 1 ], s0 );
223227 rc = rc ? rc : grb::setElement (energies, msg[ 0 ].e , n_replicas - 1 );
224228 }else if ( si == s+2 ){
@@ -294,6 +298,7 @@ namespace grb {
294298 ){
295299
296300 const size_t s = spmd<>::pid ();
301+ const size_t n_procs = spmd<>::nprocs ();
297302 const size_t n_replicas = states.size ();
298303 const size_t n = grb::size (states[0 ]);
299304 (void ) n;
@@ -337,7 +342,7 @@ namespace grb {
337342 } // n_replicas
338343 if ( rc == SUCCESS && use_pt ){
339344 // do a Parallel Tempering move
340- rc = pt< backend >( states, energies, betas );
345+ rc = pt< backend >( states, energies, betas, i_sweep*n_procs + s );
341346 }
342347#ifndef NDEBUG
343348 if ( s == 0 ) {
@@ -355,6 +360,7 @@ namespace grb {
355360 if ( rc == SUCCESS ){
356361 rc = rc ? rc : grb::collectives<>::allreduce (
357362 best_energy, grb::operators::min< EnergyType >() );
363+ // TODO: update best state to match best energy
358364 }
359365
360366 return rc;
@@ -481,8 +487,6 @@ namespace grb {
481487 *
482488 * states should be a vector of already initialized and filled dense grb::Vector.
483489 *
484- * TODO: expand and complete documentation
485- *
486490 * Warning: This function allocates $O(n)$ memory for temporary vectors.
487491 *
488492 * @param[in,out] states On input: initial (dense) states.
@@ -536,7 +540,6 @@ namespace grb {
536540 (void ) s;
537541 grb::RC rc = grb::SUCCESS;
538542
539- assert ( grb::size (states[0 ]) == n );
540543 assert ( grb::nnz (states[0 ]) == n ); // state is dense
541544 assert ( states.size () == n_replicas );
542545
@@ -576,21 +579,20 @@ namespace grb {
576579
577580 grb::Vector< QType, backend > h ( n );
578581 grb::Vector< QType, backend > rand ( n );
579- grb::Vector< StateType , backend > delta ( n );
580- grb::Vector< EnergyType , backend > dn ( n );
582+ grb::Vector< QType , backend > delta ( n );
583+ grb::Vector< QType , backend > dn ( n );
581584 grb::Vector< bool , backend > accept ( n );
582- std::srand ( static_cast <unsigned >( seed ) );
583585 std::minstd_rand rng ( seed ); // minstd_rand or std::mt19937
584586
585- grb::resize ( h, n );
586- grb::resize ( rand, n );
587- grb::resize ( delta, n );
588- grb::resize ( dn, n );
589- grb::resize ( accept, n );
587+ rc = rc ? rc : grb::resize ( h, n );
588+ rc = rc ? rc : grb::resize ( rand, n );
589+ rc = rc ? rc : grb::resize ( delta, n );
590+ rc = rc ? rc : grb::resize ( dn, n );
591+ rc = rc ? rc : grb::resize ( accept, n );
590592
591593 std::vector< grb::Vector< bool , backend > > masks ;
592594 rc = rc ? rc : matrix_partition< descr >( masks, couplings, h, rand, seed );
593- grb::clear (h);
595+ rc = rc ? rc : grb::clear (h);
594596 constexpr auto dense_descr = descr | grb::descriptors::dense;
595597
596598 auto sweep_data = std::tie (
@@ -632,7 +634,9 @@ namespace grb {
632634 const size_t n = grb::size ( state );
633635 EnergyType delta_energy = static_cast < EnergyType >(0.0 );
634636 grb::RC rc = grb::SUCCESS;
635- (void ) n;
637+
638+ assert ( grb::nnz (state) == n ); // state has to be dense!
639+ assert ( grb::nnz (local_fields) == n );
636640
637641 if ( !empty_local_fields) {
638642 rc = rc ? rc : grb::set< descr >( h, local_fields );
@@ -644,7 +648,7 @@ namespace grb {
644648 std::exponential_distribution< EnergyType > rand_gen ( beta );
645649 for ( size_t i = 0 ; i < n; ++i ){
646650 const auto rnd = -rand_gen ( rng );
647- grb::setElement ( rand, rnd, i );
651+ rc = rc ? rc : grb::setElement ( rand, rnd, i );
648652 }
649653
650654 const grb::operators::leq< EnergyType > leq_operator;
0 commit comments