Skip to content

Commit 9014cad

Browse files
IthanilIthanil
committed
Added possibility to use fixed blocking and fixed number of findMRT2step/decorrelation steps in MCI
1 parent 6daa792 commit 9014cad

File tree

2 files changed

+111
-59
lines changed

2 files changed

+111
-59
lines changed

‎src/MCIntegrator.cpp

Lines changed: 99 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -6,51 +6,73 @@
66
#include <algorithm>
77
#include <vector>
88
#include <random>
9+
#include <stdexcept>
910

1011

1112
// --- Integrate
1213

13-
void MCI::integrate(const long &Nmc, double * average, double * error, bool findMRT2step, bool initialdecorrelation)
14+
void MCI::integrate(const long &Nmc, double * average, double * error, bool findMRT2step, bool initialdecorrelation, size_t nblocks)
1415
{
15-
long i;
16+
MCI::integrate(Nmc, average, error, findMRT2step ? -1 : 0, initialdecorrelation ? -1 : 0, nblocks);
17+
}
18+
19+
void MCI::integrate(const long &Nmc, double * average, double * error, int NfindMRT2stepIterations, int NdecorrelationSteps, size_t nblocks)
20+
{
21+
long i,j;
22+
const bool fixedBlocks = (nblocks>0);
23+
const long stepsPerBlock = fixedBlocks ? Nmc/nblocks : 1;
24+
const long trueNmc = fixedBlocks ? stepsPerBlock*nblocks : Nmc;
25+
const long ndatax = fixedBlocks ? nblocks : Nmc;
1626

1727
if ( _flagpdf )
1828
{
1929
//find the optimal mrt2 step
20-
if (findMRT2step) this->findMRT2Step();
30+
this->findMRT2Step(NfindMRT2stepIterations);
2131
// take care to do the initial decorrelation of the walker
22-
if (initialdecorrelation) this->initialDecorrelation();
32+
this->initialDecorrelation(NdecorrelationSteps);
2333
}
2434

2535
//allocation of the array where the data will be stored
26-
_datax = new double*[Nmc];
27-
for (i=0; i<Nmc; ++i){ *(_datax+i) = new double[_nobsdim]; }
36+
_datax = new double*[ndatax];
37+
for (i=0; i<ndatax; ++i) { *(_datax+i) = new double[_nobsdim]; }
2838

2939
//sample the observables
3040
if (_flagobsfile) _obsfile.open(_pathobsfile);
3141
if (_flagwlkfile) _wlkfile.open(_pathwlkfile);
3242
_flagMC = true;
33-
this->sample(Nmc, true);
43+
this->sample(trueNmc, true, stepsPerBlock);
3444
_flagMC = false;
3545
if (_flagobsfile) _obsfile.close();
3646
if (_flagwlkfile) _wlkfile.close();
3747

48+
//reduce block averages
49+
if (fixedBlocks) {
50+
const double fac = 1./stepsPerBlock;
51+
for (i=0; i<ndatax; ++i) {
52+
for (j=0; j<_nobsdim; ++j) {
53+
*(*(_datax+i)+j) *= fac;
54+
}
55+
}
56+
}
57+
3858
//estimate average and standard deviation
39-
if ( _flagpdf )
59+
if ( _flagpdf && !fixedBlocks)
4060
{
41-
mci::MultiDimCorrelatedEstimator(Nmc, _nobsdim, _datax, average, error);
61+
mci::MultiDimCorrelatedEstimator(ndatax, _nobsdim, _datax, average, error);
4262
}
4363
else
4464
{
45-
mci::MultiDimUncorrelatedEstimator(Nmc, _nobsdim, _datax, average, error);
46-
for (i=0; i<_nobsdim; ++i)
47-
{
48-
*(average+i) *=_vol; *(error+i) *=_vol;
65+
mci::MultiDimUncorrelatedEstimator(ndatax, _nobsdim, _datax, average, error);
66+
if (!_flagpdf) {
67+
for (i=0; i<_nobsdim; ++i) {
68+
*(average+i) *=_vol;
69+
*(error+i) *=_vol;
4970
}
71+
}
5072
}
5173

5274
//deallocation of the data array
53-
for (int i=0; i<Nmc; ++i){ delete [] *(_datax+i); }
75+
for (int i=0; i<ndatax; ++i){ delete [] *(_datax+i); }
5476
delete [] _datax;
5577
}
5678

@@ -64,10 +86,12 @@ void MCI::storeObservables()
6486
if ( _ridx%_freqobsfile == 0 )
6587
{
6688
_obsfile << _ridx;
67-
for (int j=0; j<_nobsdim; ++j)
89+
for (size_t i=0; i<_obs.size(); ++i) {
90+
for (int j=0; j<_obs[i]->getNObs(); ++j)
6891
{
69-
_obsfile << " " << *(*(_datax+_ridx)+j) ;
92+
_obsfile << " " << _obs[i]->getObservable(j);
7093
}
94+
}
7195
_obsfile << std::endl;
7296
}
7397
}
@@ -87,52 +111,72 @@ void MCI::storeWalkerPositions()
87111
}
88112

89113

90-
void MCI::initialDecorrelation()
114+
void MCI::initialDecorrelation(const int &NdecorrelationSteps)
91115
{
92-
long i;
93-
//constants
94-
const long MIN_NMC=100;
95-
//allocate the data array that will be used
96-
_datax = new double*[MIN_NMC];
97-
for (i=0; i<MIN_NMC; ++i) {*(_datax+i) = new double[_nobsdim]; }
98-
//do a first estimate of the observables
99-
this->sample(MIN_NMC, true);
100-
double * oldestimate = new double[_nobsdim];
101-
double * olderrestim = new double[_nobsdim];
102-
mci::MultiDimCorrelatedEstimator(MIN_NMC, _nobsdim, _datax, oldestimate, olderrestim);
103-
//start a loop which will stop when the observables are stabilized
104-
bool flag_loop=true;
105-
double * newestimate = new double[_nobsdim];
106-
double * newerrestim = new double[_nobsdim];
107-
double * foo;
108-
while ( flag_loop )
109-
{
116+
if (NdecorrelationSteps < 0) {
117+
long i;
118+
//constants
119+
const long MIN_NMC=100;
120+
//allocate the data array that will be used
121+
_datax = new double*[MIN_NMC];
122+
for (i=0; i<MIN_NMC; ++i) { *(_datax+i) = new double[_nobsdim]; }
123+
//do a first estimate of the observables
124+
this->sample(MIN_NMC, true);
125+
double * oldestimate = new double[_nobsdim];
126+
double * olderrestim = new double[_nobsdim];
127+
mci::MultiDimCorrelatedEstimator(MIN_NMC, _nobsdim, _datax, oldestimate, olderrestim);
128+
//start a loop which will stop when the observables are stabilized
129+
bool flag_loop=true;
130+
double * newestimate = new double[_nobsdim];
131+
double * newerrestim = new double[_nobsdim];
132+
double * foo;
133+
while ( flag_loop ) {
134+
flag_loop = false;
110135
this->sample(MIN_NMC, true);
111136
mci::MultiDimCorrelatedEstimator(MIN_NMC, _nobsdim, _datax, newestimate, newerrestim);
112-
for (i=0; i<_nobsdim; ++i)
113-
{
114-
if ( abs( *(oldestimate+i) - *(newestimate+i) ) <= *(olderrestim+i) + *(newerrestim+i) ) flag_loop=false;
115-
if (flag_loop) *(oldestimate+i) = *(newestimate+i);
137+
for (i=0; i<_nobsdim; ++i) {
138+
if ( std::abs( *(oldestimate+i) - *(newestimate+i) ) > *(olderrestim+i) + *(newerrestim+i) ) {
139+
flag_loop=true;
116140
}
141+
}
117142
foo=oldestimate;
118143
oldestimate=newestimate;
119144
newestimate=foo;
120145
}
121-
//memory deallocation
122-
delete [] newestimate;
123-
delete [] newerrestim;
124-
delete [] oldestimate;
125-
delete [] olderrestim;
126-
for (i=0; i<MIN_NMC; ++i){ delete[] *(_datax+i); }
127-
delete [] _datax;
146+
//memory deallocation
147+
delete [] newestimate;
148+
delete [] newerrestim;
149+
delete [] oldestimate;
150+
delete [] olderrestim;
151+
for (i=0; i<MIN_NMC; ++i){ delete[] *(_datax+i); }
152+
delete [] _datax;
153+
}
154+
else {
155+
this->sample(NdecorrelationSteps, false);
156+
}
128157
}
129158

130159

131-
void MCI::sample(const long &npoints, const bool &flagobs)
160+
void MCI::sample(const long &npoints, const bool &flagobs, const long &stepsPerBlock)
132161
{
162+
if (flagobs && stepsPerBlock>0 && npoints%stepsPerBlock!=0) {
163+
throw std::invalid_argument("If fixed blocking is used, npoints must be a multiple of stepsPerBlock.");
164+
}
165+
133166
int i;
134-
//initialize the running index
167+
//initialize the running indices
135168
_ridx=0;
169+
_bidx=0;
170+
171+
if (flagobs) {
172+
// set the data to 0
173+
for (i=0; i<npoints/stepsPerBlock; ++i) {
174+
for (long j=0; j<_nobsdim; ++j) {
175+
*(*(_datax+i)+j) = 0.;
176+
}
177+
}
178+
}
179+
136180
//initialize the pdf at x
137181
computeOldSamplingFunction();
138182
//initialize the observables values at x
@@ -170,13 +214,15 @@ void MCI::sample(const long &npoints, const bool &flagobs)
170214
if (_flagwlkfile) this->storeWalkerPositions();
171215
}
172216
_ridx++;
217+
_bidx = _ridx / stepsPerBlock;
173218
}
174219
} else
175220
{
176221
for (i=0; i<npoints; ++i)
177222
{
178223
this->doStepMRT2(&flagacc);
179224
_ridx++;
225+
_bidx = _ridx / stepsPerBlock;
180226
}
181227
}
182228
}
@@ -194,20 +240,22 @@ void MCI::sample(const long &npoints, const bool &flagobs)
194240
if (_flagwlkfile) this->storeWalkerPositions();
195241
}
196242
_ridx++;
243+
_bidx = _ridx / stepsPerBlock;
197244
}
198245
} else
199246
{
200247
for (i=0; i<npoints; ++i)
201248
{
202249
this->newRandomX();
203250
_ridx++;
251+
_bidx = _ridx / stepsPerBlock;
204252
}
205253
}
206254
}
207255
}
208256

209257

210-
void MCI::findMRT2Step()
258+
void MCI::findMRT2Step(const int &NfindMRT2stepIterations)
211259
{
212260
int j;
213261
//constants
@@ -221,7 +269,7 @@ void MCI::findMRT2Step()
221269
int cons_count = 0; //number of consecutive loops without need of changing mrt2step
222270
int counter = 0; //counter of loops
223271
double fact;
224-
while ( cons_count < MIN_CONS )
272+
while ( (cons_count < MIN_CONS && NfindMRT2stepIterations < 0) || counter < NfindMRT2stepIterations)
225273
{
226274
counter++;
227275
//do MIN_STAT M(RT)^2 steps
@@ -258,7 +306,6 @@ void MCI::findMRT2Step()
258306
{
259307
if ( *(_mrt2step+j) - ( *(*(_irange+j)+1) - *(*(_irange+j)) ) > 0. )
260308
{
261-
cons_count = MIN_CONS; //make the main loop terminate
262309
*(_mrt2step+j) = ( *(*(_irange+j)+1) - *(*(_irange+j)) );
263310
}
264311
}
@@ -267,7 +314,6 @@ void MCI::findMRT2Step()
267314
{
268315
if ( *(_mrt2step+j) < SMALLEST_ACCEPTABLE_DOUBLE )
269316
{
270-
cons_count = MIN_CONS; //make the main loop terminate
271317
*(_mrt2step+j) = SMALLEST_ACCEPTABLE_DOUBLE;
272318
}
273319
}
@@ -419,7 +465,7 @@ void MCI::saveObservables()
419465
{
420466
for (int j=0; j<_obs[i]->getNObs(); ++j)
421467
{
422-
_datax[_ridx][idx]=_obs[i]->getObservable(j);
468+
_datax[_bidx][idx]+=_obs[i]->getObservable(j);
423469
idx++;
424470
}
425471
}

‎src/MCIntegrator.hpp

Lines changed: 12 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -39,8 +39,9 @@ class MCI
3939

4040
std::vector<MCICallBackOnAcceptanceInterface *> _cback; // Vector of observable functions
4141

42-
int _ridx; // running index, which keeps track of the index of datax
43-
double ** _datax; // array that will contain all the measured observable
42+
int _ridx; // running index, which keeps track of the number of MC steps
43+
int _bidx; // index of the current block/datax element
44+
double ** _datax; // array that will contain all the measured observable (or block averages if used)
4445
bool _flagMC; //flag that is true only when MCI is accumulating data for the integral
4546

4647
std::ofstream _obsfile; //ofstream for storing obs values while sampling
@@ -73,11 +74,11 @@ class MCI
7374
void updateX();
7475
void doStepMRT2(bool * flagacc); //use this if there is a pdf
7576

76-
void findMRT2Step();
77+
void findMRT2Step(const int &NfindMRT2stepIterations = -1);
7778

78-
void sample(const long &npoints, const bool &flagobs);
79+
void sample(const long &npoints, const bool &flagobs, const long &stepsPerBlock = 1);
7980

80-
void initialDecorrelation();
81+
void initialDecorrelation(const int &NdecorrelationSteps = -1);
8182

8283
void storeObservables();
8384
void storeWalkerPositions();
@@ -129,7 +130,12 @@ class MCI
129130
double getAcceptanceRate(){return (double(_acc)/(double(_acc+_rej)));}
130131

131132
// --- Integrate
132-
void integrate(const long &Nmc, double * average, double * error, bool findMRT2step=true, bool initialdecorrelation=true);
133+
134+
// Wrapper for easiness/compatibility. Using automatic methods for findMRT2step and initial decorrelation (if enabled).
135+
void integrate(const long &Nmc, double * average, double * error, bool findMRT2step=true, bool initialdecorrelation = true, size_t nblocks = 0); // nblocks == 0 means using auto-blocking (no RAM benefit)
136+
137+
// Actual integrate implemention. With integer controls: -1 -> auto, 0 -> disabled, >0 -> use fixed step/iteration counts, e.g. useful for parallel computation
138+
void integrate(const long &Nmc, double * average, double * error, int NfindMRT2stepIterations, int NdecorrelationSteps, size_t nblocks = 0);
133139

134140

135141
};

0 commit comments

Comments
 (0)