1 //////////////////////////////////////////////////////////////////////////////////////
2 // This file is distributed under the University of Illinois/NCSA Open Source License.
3 // See LICENSE file in top directory for details.
4 //
5 // Copyright (c) 2020 QMCPACK developers.
6 //
7 // File developed by: Bryan Clark, bclark@Princeton.edu, Princeton University
8 //                    Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
9 //                    Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
10 //                    Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
11 //                    Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory
12 //                    Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
13 //
14 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
15 //////////////////////////////////////////////////////////////////////////////////////
16 
17 
18 /** @file EstimatorManagerBase.h
19  * @brief Manager class of scalar estimators
20  */
21 #ifndef QMCPLUSPLUS_ESTIMATORMANAGERBASE_H
22 #define QMCPLUSPLUS_ESTIMATORMANAGERBASE_H
23 
24 #include "Configuration.h"
25 #include "Utilities/Timer.h"
26 #include "Utilities/PooledData.h"
27 #include "Message/Communicate.h"
28 #include "Estimators/ScalarEstimatorBase.h"
29 #include "Particle/Walker.h"
30 #include "OhmmsPETE/OhmmsVector.h"
31 #include "OhmmsData/HDFAttribIO.h"
32 #include <bitset>
33 
34 namespace qmcplusplus
35 {
36 class MCWalkerConfiguration;
37 class QMCHamiltonian;
38 class CollectablesEstimator;
39 
40 namespace testing
41 {
42 class EstimatorManagerBaseTest;
43 } // namespace testing
44 
45 
46 /** Class to manage a set of ScalarEstimators */
47 class EstimatorManagerBase
48 {
49 public:
50   typedef QMCTraits::FullPrecRealType RealType;
51   using FullPrecRealType = QMCTraits::FullPrecRealType;
52 
53   typedef ScalarEstimatorBase EstimatorType;
54   typedef std::vector<RealType> BufferType;
55   using MCPWalker = Walker<QMCTraits, PtclOnLatticeTraits>;
56 
57   //enum { WEIGHT_INDEX=0, BLOCK_CPU_INDEX, ACCEPT_RATIO_INDEX, TOTAL_INDEX};
58 
59   ///name of the primary estimator name
60   std::string MainEstimatorName;
61   ///the root file name
62   std::string RootName;
63   ///energy
64   TinyVector<RealType, 4> RefEnergy;
65   // //Cummulative energy, weight and variance
66   // TinyVector<RealType,4>  EPSum;
67   ///default constructor
68   EstimatorManagerBase(Communicate* c = 0);
69   ///copy constructor
70   EstimatorManagerBase(EstimatorManagerBase& em);
71   ///destructor
72   virtual ~EstimatorManagerBase();
73 
74   /** set the communicator */
75   void setCommunicator(Communicate* c);
76 
77   /** return the communicator
78    */
getCommunicator()79   Communicate* getCommunicator() { return myComm; }
80 
81   /** return true if the rank == 0
82    */
is_manager()83   inline bool is_manager() const { return !myComm->rank(); }
84 
85   ///return the number of ScalarEstimators
size()86   inline int size() const { return Estimators.size(); }
87 
88   /** add a property with a name
89    * @param aname name of the column
90    * @return the property index so that its value can be set by setProperty(i)
91    *
92    * Append a named column. BlockProperties do not contain any meaning data
93    * but manages the name to index map for PropertyCache.
94    */
addProperty(const char * aname)95   inline int addProperty(const char* aname) { return BlockProperties.add(aname); }
96 
97   /** set the value of the i-th column with a value v
98    * @param i column index
99    * @param v value
100    */
setProperty(int i,RealType v)101   inline void setProperty(int i, RealType v) { PropertyCache[i] = v; }
102 
getProperty(int i)103   inline RealType getProperty(int i) const { return PropertyCache[i]; }
104 
105   int addObservable(const char* aname);
106 
getObservable(int i)107   inline RealType getObservable(int i) const { return TotalAverages[i]; }
108 
109   void getData(int i, std::vector<RealType>& values);
110 
111   /** add an Estimator
112    * @param newestimator New Estimator
113    * @param aname name of the estimator
114    * @return locator of newestimator
115    */
116   int add(EstimatorType* newestimator, const std::string& aname);
117   //int add(CompositeEstimatorBase* newestimator, const std::string& aname);
118 
119   /** add a main estimator
120    * @param newestimator New Estimator
121    * @return locator of newestimator
122    */
add(EstimatorType * newestimator)123   int add(EstimatorType* newestimator) { return add(newestimator, MainEstimatorName); }
124 
125   ///return a pointer to the estimator aname
126   EstimatorType* getEstimator(const std::string& a);
127 
128   ///return a pointer to the estimator
129   EstimatorType* getMainEstimator();
130 
131   void setCollectionMode(bool collect);
132   //void setAccumulateMode (bool setAccum) {AccumulateBlocks = setAccum;};
133 
134   ///process xml tag associated with estimators
135   //bool put(xmlNodePtr cur);
136   bool put(QMCHamiltonian& H, xmlNodePtr cur);
137 
138   void resetTargetParticleSet(ParticleSet& p);
139 
140   /** reset the estimator
141    */
142   void reset();
143 
144   /** start a run
145    * @param blocks number of blocks
146    * @param record if true, will write to a file
147    *
148    * Replace reportHeader and reset functon.
149    */
150   void start(int blocks, bool record = true);
151   /** stop a qmc run
152    *
153    * Replace finalize();
154    */
155   void stop();
156   /** stop a qmc run
157    */
158   void stop(const std::vector<EstimatorManagerBase*> m);
159 
160 
161   /** start  a block
162    * @param steps number of steps in a block
163    */
164   void startBlock(int steps);
165 
166   /** stop a block
167    * @param accept acceptance rate of this block
168    */
169   void stopBlock(RealType accept, bool collectall = true);
170 
171   /** stop a block
172    * @param m list of estimator which has been collecting data independently
173    */
174   void stopBlock(const std::vector<EstimatorManagerBase*>& m);
175 
176   /** accumulate the measurements
177    * @param W walkers
178    */
179   void accumulate(MCWalkerConfiguration& W);
180 
181   /** accumulate the measurements for a subset of walkers [it,it_end)
182    * @param W walkers
183    * @param it first walker
184    * @param it_end last walker
185    */
186   void accumulate(MCWalkerConfiguration& W, MCWalkerConfiguration::iterator it, MCWalkerConfiguration::iterator it_end);
187 
188   /** get the average of per-block energy and variance of all the blocks
189    * Note: this is not weighted average. It can be the same as weighted average only when block weights are identical.
190    */
191   void getApproximateEnergyVariance(RealType& e, RealType& var);
192 
193   template<class CT>
write(CT & anything,bool doappend)194   void write(CT& anything, bool doappend)
195   {
196     anything.write(h_file, doappend);
197   }
198 
get_AverageCache()199   auto& get_AverageCache() { return AverageCache; }
get_SquaredAverageCache()200   auto& get_SquaredAverageCache() { return SquaredAverageCache; }
201 
202 protected:
203   //  TODO: fix needless use of bitset instead of clearer more visible booleans
204   std::bitset<8> Options;
205   ///size of the message buffer
206   int BufferSize;
207   ///number of records in a block
208   int RecordCount;
209   ///index for the block weight PropertyCache(weightInd)
210   int weightInd;
211   ///index for the block cpu PropertyCache(cpuInd)
212   int cpuInd;
213   ///index for the acceptance rate PropertyCache(acceptInd)
214   int acceptInd;
215   ///hdf5 handler
216   hid_t h_file;
217   ///total weight accumulated in a block
218   RealType BlockWeight;
219   ///file handler to write data
220   std::ofstream* Archive;
221   ///file handler to write data for debugging
222   std::ofstream* DebugArchive;
223   ///communicator to handle communication
224   Communicate* myComm;
225   /** pointer to the primary ScalarEstimatorBase
226    */
227   ScalarEstimatorBase* MainEstimator;
228   /** pointer to the CollectablesEstimator
229    *
230    * Do not need to clone: owned by the master thread
231    */
232   CollectablesEstimator* Collectables;
233   /** accumulator for the energy
234    *
235    * @todo expand it for all the scalar observables to report the final results
236    */
237   ScalarEstimatorBase::accumulator_type energyAccumulator;
238   /** accumulator for the variance **/
239   ScalarEstimatorBase::accumulator_type varAccumulator;
240   ///cached block averages of the values
241   Vector<RealType> AverageCache;
242   ///cached block averages of the squared values
243   Vector<RealType> SquaredAverageCache;
244   ///cached block averages of properties, e.g. BlockCPU
245   Vector<RealType> PropertyCache;
246   ///manager of scalar data
247   RecordNamedProperty<RealType> BlockAverages;
248   ///manager of property data
249   RecordNamedProperty<RealType> BlockProperties;
250   ///block averages: name to value
251   RecordNamedProperty<RealType> TotalAverages;
252   ///data accumulated over the blocks
253   Matrix<RealType> TotalAveragesData;
254   ///index mapping between BlockAverages and TotalAverages
255   std::vector<int> Block2Total;
256   ///column map
257   std::map<std::string, int> EstimatorMap;
258   ///estimators of simple scalars
259   std::vector<EstimatorType*> Estimators;
260   ///convenient descriptors for hdf5
261   std::vector<observable_helper*> h5desc;
262   /////estimators of composite data
263   //CompositeEstimatorSet* CompEstimators;
264   ///Timer
265   Timer MyTimer;
266 
267 private:
268   ///number of maximum data for a scalar.dat
269   int max4ascii;
270   //Data for communication
271   std::vector<BufferType*> RemoteData;
272   ///collect data and write
273   void collectBlockAverages();
274   ///add header to an std::ostream
275   void addHeader(std::ostream& o);
276   size_t FieldWidth;
277 
278   friend class qmcplusplus::testing::EstimatorManagerBaseTest;
279 };
280 } // namespace qmcplusplus
281 #endif
282