CHROMA
inline_prop_3pt_w.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Inline measurement 3pt_prop
3  *
4  */
5 
6 #include "handle.h"
15 #include "meas/sources/zN_src.h"
23 #include "meas/glue/mesplq.h"
24 #include "util/ft/sftmom.h"
25 #include "util/info/proginfo.h"
27 #include "util/info/unique_id.h"
28 #include "util/ferm/transf.h"
30 
31 namespace Chroma{
32  namespace InlineProp3ptEnv{
33  namespace{
34  AbsInlineMeasurement* createMeasurement(XMLReader& xml_in,
35  const std::string& path)
36  {
37  return new InlineMeas(Params(xml_in, path));
38  }
39 
40  //! Local registration flag
41  bool registered = false;
42  }
43 
44  const std::string name = "PROPAGATOR_3PT";
45 
46  //! Register all the factories
47  bool registerAll()
48  {
49  bool success = true;
50  if (! registered)
51  {
52  //success &= BaryonOperatorEnv::registerAll();
53  success &= TheInlineMeasurementFactory::Instance().registerObject(name, createMeasurement);
54  registered = true;
55  }
56  return success;
57  }
58 
59  void read(XMLReader& xml, const std::string& path, Params::Operator_t& op){
60  XMLReader paramtop(xml, path);
61 
62  read(paramtop, "gamma", op.gamma);
63  read(paramtop, "p", op.p);
64  read(paramtop, "t", op.t);
65  if(paramtop.count("factor")>0)
66  read(paramtop, "factor", op.f);
67  else
68  op.f = 1.0 ;
69  }
70 
71  // Writter for input parameters
72  void write(XMLWriter& xml, const std::string& path, const Params::Operator_t& op){
73  push(xml, path);
74 
75  write(xml, "p", op.p);
76  write(xml, "t", op.t);
77  write(xml, "gamma", op.gamma);
78  write(xml, "factor", op.f);
79 
80  pop(xml);
81  }
82 
83  // Reader for input parameters
84  void read(XMLReader& xml, const std::string& path, Params::Param_t& param){
85  XMLReader paramtop(xml, path);
86 
87  int version;
88  read(paramtop, "version", version);
89 
90  switch (version)
91  {
92  case 1:
93  /************************************************************/
94  read(paramtop,"op",param.op);
95  param.chi = readXMLArrayGroup(paramtop, "Quarks", "DilutionType");
96 
97  break;
98 
99  default :
100  /**************************************************************/
101 
102  QDPIO::cerr << "Input parameter version " << version << " unsupported." << std::endl;
103  QDP_abort(1);
104  }
105 
106 
107  }
108 
109 
110  // Writter for input parameters
111  void write(XMLWriter& xml, const std::string& path, const Params::Param_t& param){
112  push(xml, path);
113 
114  int version = 1;
115 
116  write(xml, "version", version);
117  write(xml, "op", param.op);
118 
119  push(xml,"Quarks");
120  for( int t(0);t<param.chi.size();t++){
121  push(xml,"elem");
122  xml<<param.chi[t].xml;
123  pop(xml);
124  }
125  pop(xml);
126 
127  pop(xml); // final pop
128  }
129 
130 
131  //! Gauge field parameters
132  void read(XMLReader& xml, const std::string& path, Params::NamedObject_t& input)
133  {
134  XMLReader inputtop(xml, path);
135 
136  read(inputtop, "gauge_id", input.gauge_id);
137  read(inputtop, "prop_id", input.prop_id);
138  read(inputtop, "prop3pt_id", input.prop3pt_id);
139  }
140 
141  //! Gauge field parameters
142  void write(XMLWriter& xml, const std::string& path, const Params::NamedObject_t& input){
143  push(xml, path);
144 
145  write(xml, "gauge_id", input.gauge_id);
146  write(xml, "prop_id", input.prop_id);
147  write(xml, "prop3pt_id", input.prop3pt_id);
148  pop(xml);
149  }
150 
151 
152  // Param stuff
154  frequency = 0;
155  param.op.gamma = 0 ;
156  param.op.f = 1.0 ;
157  }
158 
159  Params::Params(XMLReader& xml_in, const std::string& path)
160  {
161  try
162  {
163  XMLReader paramtop(xml_in, path);
164 
165  if (paramtop.count("Frequency") == 1)
166  read(paramtop, "Frequency", frequency);
167  else
168  frequency = 1;
169 
170  // Read program parameters
171  read(paramtop, "Param", param);
172 
173  // Read in the output propagator/source configuration info
174  read(paramtop, "NamedObject", named_obj);
175 
176  // Possible alternate XML file pattern
177  if (paramtop.count("xml_file") != 0)
178  {
179  read(paramtop, "xml_file", xml_file);
180  }
181  }
182  catch(const std::string& e)
183  {
184  QDPIO::cerr << __func__ << ": Caught Exception reading XML: " << e << std::endl;
185  QDP_abort(1);
186  }
187  }
188 
189 
190  void Params::write(XMLWriter& xml_out, const std::string& path)
191  {
192  push(xml_out, path);
193 
194  // Parameters for source construction
195  InlineProp3ptEnv::write(xml_out, "Param", param);
196 
197  // Write out the output propagator/source configuration info
198  InlineProp3ptEnv::write(xml_out, "NamedObject", named_obj);
199 
200  pop(xml_out);
201  }
202 
203 
204 
205  //--------------------------------------------------------------
206  // Function call
207  // void
208  //InlineProp3pt::operator()(unsigned long update_no,
209  // XMLWriter& xml_out)
210  void InlineMeas::operator()(unsigned long update_no,
211  XMLWriter& xml_out)
212  {
213  // If xml file not empty, then use alternate
214  if (params.xml_file != ""){
215  std::string xml_file = makeXMLFileName(params.xml_file, update_no);
216 
217  push(xml_out, "propagator_3pt");
218  write(xml_out, "update_no", update_no);
219  write(xml_out, "xml_file", xml_file);
220  pop(xml_out);
221 
222  XMLFileWriter xml(xml_file);
223  func(update_no, xml);
224  }
225  else{
226  func(update_no, xml_out);
227  }
228  }
229 
230 
231  // Function call
232  //void
233  //InlineProp3pt::func(unsigned long update_no,
234  // XMLWriter& xml_out)
235  void InlineMeas::func(unsigned long update_no,
236  XMLWriter& xml_out)
237  {
238  START_CODE();
239 
240  StopWatch snoop;
241  snoop.reset();
242  snoop.start();
243 
244  // Test and grab a reference to the gauge field
245  XMLBufferWriter gauge_xml;
246  try
247  {
248  TheNamedObjMap::Instance().getData< multi1d<LatticeColorMatrix> >(params.named_obj.gauge_id);
249  TheNamedObjMap::Instance().get(params.named_obj.gauge_id).getRecordXML(gauge_xml);
250  }
251  catch( std::bad_cast )
252  {
253  QDPIO::cerr << InlineProp3ptEnv::name << ": caught dynamic cast error"
254  << std::endl;
255  QDP_abort(1);
256  }
257  catch (const std::string& e)
258  {
259  QDPIO::cerr << name << ": std::map call failed: " << e
260  << std::endl;
261  QDP_abort(1);
262  }
263  const multi1d<LatticeColorMatrix>& u =
264  TheNamedObjMap::Instance().getData< multi1d<LatticeColorMatrix> >(params.named_obj.gauge_id);
265 
266  push(xml_out, "prop_3pt");
267  write(xml_out, "update_no", update_no);
268 
269  QDPIO::cout << name << ": Propagator for 3pt functions" << std::endl;
270 
271  proginfo(xml_out); // Print out basic program info
272 
273  // Write out the input
274  params.write(xml_out, "Input");
275 
276  // Write out the config info
277  write(xml_out, "Config_info", gauge_xml);
278 
279  push(xml_out, "Output_version");
280  write(xml_out, "out_version", 1);
281  pop(xml_out);
282 
283  // First calculate some gauge invariant observables just for info.
284  // This is really cheap.
285  MesPlq(xml_out, "Observables", u);
286 
287  // Save current seed
288  Seed ran_seed;
289  QDP::RNG::savern(ran_seed);
290 
291  //
292  // Read the source and solutions
293  //
294  StopWatch swatch;
295  swatch.reset();
296  swatch.start();
297 
298  int N_quarks = params.param.chi.size() ;
299 
300  multi1d< Handle< DilutionScheme<LatticeFermion> > > quarks(N_quarks);
301 
302  try{
303  // Loop over quark dilutions
304  for(int n(0); n < params.param.chi.size(); ++n){
305  const GroupXML_t& dil_xml = params.param.chi[n];
306  std::istringstream xml_d(dil_xml.xml);
307  XMLReader diltop(xml_d);
308  QDPIO::cout << "Dilution type = " << dil_xml.id << std::endl;
309  quarks[n] =
310  TheFermDilutionSchemeFactory::Instance().createObject(dil_xml.id,
311  diltop,
312  dil_xml.path);
313  }
314  }
315  catch(const std::string& e){
316  QDPIO::cerr << name << ": Caught Exception constructing dilution scheme: " << e << std::endl;
317  QDP_abort(1);
318  }
319 
320 
321  //-------------------------------------------------------------------
322  //Sanity checks
323 
324  //Another Sanity check, the three quarks must all be
325  //inverted on the same cfg
326  for (int n = 1 ; n < N_quarks ; ++n){
327  if (quarks[0]->getCfgInfo() != quarks[n]->getCfgInfo()){
328  QDPIO::cerr << name << " : Quarks do not contain the same cfg info";
329  QDPIO::cerr << ", quark "<< n << std::endl;
330  QDP_abort(1);
331  }
332  }
333 
334 
335  //Also ensure that the cfg on which the inversions were performed
336  //is the same as the cfg that we are using
337  {
338  std::string cfgInfo;
339 
340  //Really ugly way of doing this(Robert Heeeelp!!)
341  XMLBufferWriter top;
342  write(top, "Config_info", gauge_xml);
343  XMLReader from(top);
344  XMLReader from2(from, "/Config_info");
345  std::ostringstream os;
346  from2.print(os);
347 
348  cfgInfo = os.str();
349 
350  if (cfgInfo != quarks[0]->getCfgInfo()){
351  QDPIO::cerr << name << " : Quarks do not contain the same";
352  QDPIO::cerr << " cfg info as the gauge field." ;
353  QDPIO::cerr << "gauge: XX"<<cfgInfo<<"XX quarks: XX" ;
354  QDPIO::cerr << quarks[0]->getCfgInfo()<<"XX"<< std::endl;
355  QDP_abort(1);
356  }
357  }
358 
359 
360  int decay_dir = quarks[0]->getDecayDir();
361  //
362  // Initialize the slow Fourier transform phases
363  //
364  int mom2(0);
365  for(int i(0);i<Nd-1;i++)
366  mom2 += params.param.op.p[i]*params.param.op.p[i] ;
367 
368  //SftMom phases(params.param.mom2_max, false, decay_dir);
369  SftMom phases(mom2, false, decay_dir);
370 
371  // Another sanity check. The seeds of all the quarks must be different
372  // and thier decay directions must be the same
373  for(int n = 1 ; n < quarks.size(); ++n){
374  if(toBool(quarks[n]->getSeed()==quarks[0]->getSeed())){
375  QDPIO::cerr << name << ": error, quark seeds are the same" << std::endl;
376  QDP_abort(1);
377  }
378 
379  if(toBool(quarks[n]->getDecayDir()!=quarks[0]->getDecayDir())){
380  QDPIO::cerr<<name<< ": error, quark decay dirs do not match" <<std::endl;
381  QDP_abort(1);
382  }
383  }
384 
385  // Read the quark propagator and extract headers
387  PropSourceConst_t source_header;
388  QDPIO::cout << "Attempt to read propagator info" << std::endl;
389  try
390  {
391  // Try the cast to see if this is a valid source
392  LatticePropagator& tt_prop =
393  TheNamedObjMap::Instance().getData<LatticePropagator>(params.named_obj.prop_id);
394 
395  //Snarf the source info.
396  //This is will throw if the source_id is not there
397  XMLReader prop_file_xml, prop_record_xml;
398  TheNamedObjMap::Instance().get(params.named_obj.prop_id).getFileXML(prop_file_xml);
399  TheNamedObjMap::Instance().get(params.named_obj.prop_id).getRecordXML(prop_record_xml);
400  // Try to invert this record XML into a ChromaProp struct
401  // Also pull out the id of this source
402  {
403  read(prop_record_xml, "/Propagator/ForwardProp", prop_header);
404  read(prop_record_xml, "/Propagator/PropSource", source_header);
405  }
406  }
407  catch (std::bad_cast){
408  QDPIO::cerr << name << ": caught dynamic cast error" ;
409  QDPIO::cerr << std::endl;
410  QDP_abort(1);
411  }
412  catch (const std::string& e){
413  QDPIO::cerr << name << ": error extracting prop_header: " << e << std::endl;
414  QDP_abort(1);
415  }
416  const LatticePropagator& prop =
417  TheNamedObjMap::Instance().getData<LatticePropagator>(params.named_obj.prop_id);
418 
419  QDPIO::cout << "Propagator successfully read and parsed" << std::endl;
420 
421  //create the 3pt prop
422  try{
423  TheNamedObjMap::Instance().create<LatticePropagator>(params.named_obj.prop3pt_id);
424  }
425  catch (std::bad_cast){
426  QDPIO::cerr << name << ": caught dynamic cast error" << std::endl;
427  QDP_abort(1);
428  }
429  catch (const std::string& e){
430  QDPIO::cerr << name << ": error creating prop: " << e << std::endl;
431  QDP_abort(1);
432  }
433 
434  // Cast should be valid now
435  LatticePropagator& prop_3pt =
436  TheNamedObjMap::Instance().getData<LatticePropagator>(params.named_obj.prop3pt_id);
437 
438  // Need to check if the propagator gauge field is the same as the rest...
439  // NEED TO IMPLEMENT THIS
440 
441  //Print out here the operator details
442  write(xml_out, "op", params.param.op);
443  //write(xml_out, "t0", t0);
444  {
445  XMLBufferWriter top;
446  push(top, "tt");
447  write(top, "Operator", params.param.op);
448  //write(top, "t0", t0);
449  pop(top);
450  XMLReader from(top);
451  XMLReader from2(from, "/tt/Operator");
452  std::ostringstream os;
453  QDPIO::cout<<name<<" Operator: "<<std::endl ;
454  from2.print(os);
455  QDPIO::cout<<os.str()<<std::endl ;
456  //QDPIO::cout<< " Time slice: "<<t0<<std::endl ;
457  QDPIO::cout<<name<<" End Operator "<<std::endl ;
458  }
459 
460  LatticePropagator tmp_prop = Gamma(params.param.op.gamma)*prop ;
461  tmp_prop = params.param.op.f*tmp_prop;
462 
463  LatticeFermion ferm_3pt = zero ;
464  LatticeFermion ferm ;
465  LatticeComplex phase = phases[phases.momToNum(params.param.op.p)];
466  int t0 = params.param.op.t;
467  for(int s(0);s<Ns;s++)
468  for(int c(0);c<Nc;c++){
469  QDPIO::cout<<" Doing quark color and spin: "<<c<<" "<<s <<std::endl ;
470  PropToFerm(tmp_prop,ferm,c,s) ;
471  ferm_3pt = zero;
472  int count(0);
473  for(int n(0);n<quarks.size();n++){
474  int i_t0(-100) ;
475  for (int tt(0) ; tt < quarks[n]->getNumTimeSlices() ; ++tt)
476  if(quarks[n]->getT0(tt) == t0)
477  i_t0 = tt;
478  if(i_t0!=-100){ //this quark contains the appropriate t
479  count++;
480 
481  QDPIO::cout<<" Doing quark: "<<n <<std::endl ;
482  QDPIO::cout<<" quark: "<<n <<" has "<<quarks[n]->getDilSize(i_t0);
483  QDPIO::cout<<" dilutions on time slice "<<t0<<std::endl ;
484  for(int i = 0 ; i < quarks[n]->getDilSize(i_t0) ; ++i){
485  QDPIO::cout<<" Doing dilution : "<<i<<std::endl ;
486  LatticeComplex cc =
487  phase*localInnerProduct(quarks[n]->dilutedSource(i_t0,i),ferm);
488  ferm_3pt += sum(cc,phases.getSet()[t0])*quarks[n]->dilutedSolution(i_t0,i);
489  }
490  QDPIO::cout<<" Done with dilutions for quark: "<<n <<std::endl ;
491  }
492  }
493  if(count==0){
494  QDPIO::cerr<<name<< ": error, no appropriate time slice found" <<std::endl;
495  QDP_abort(1);
496 
497  }
498  ferm_3pt = ferm_3pt/Double(count) ; // compute the mean over noises
499  FermToProp(ferm_3pt,prop_3pt,c,s);
500  QDPIO::cout<<" Done with quark color and spin: "<<c<<" "<<s <<std::endl ;
501  }
502 
503  // Sanity check -
504  // write out the propagator (pion) correlator in the Nd-1 direction
505  {
506  // Initialize the slow Fourier transform phases
507  SftMom phases(0, true, Nd-1);
508 
509  multi1d<Double> corr = sumMulti(localNorm2(prop_3pt),phases.getSet());
510 
511  push(xml_out, "Prop_correlator");
512  write(xml_out, "prop_corr", corr);
513  pop(xml_out);
514  }
515 
516 
517  // Save the propagator info
518  try
519  {
520  QDPIO::cout << "Start writing propagator info" << std::endl;
521 
522  XMLBufferWriter file_xml;
523  push(file_xml, "propagator");
524  write(file_xml, "id", uniqueId()); // NOTE: new ID form
525  pop(file_xml);
526 
527  XMLBufferWriter record_xml;
528  push(record_xml , "Propagator");
529  write(record_xml, "ForwardProp", prop_header);
530  write(record_xml, "PropSource", source_header);
531  write(record_xml, "Config_info", gauge_xml);
532  //write(record_xml, "t0",t0);
533  write(record_xml, "Operator",params.param.op);
534  pop(record_xml);
535  // Write the propagator xml info
536  TheNamedObjMap::Instance().get(params.named_obj.prop3pt_id).setFileXML(file_xml);
537  TheNamedObjMap::Instance().get(params.named_obj.prop3pt_id).setRecordXML(record_xml);
538 
539  QDPIO::cout << "Propagator successfully updated" << std::endl;
540  }
541  catch (std::bad_cast){
542  QDPIO::cerr << name << ": caught dynamic cast error" << std::endl;
543  QDP_abort(1);
544  }
545  catch (const std::string& e){
546  QDPIO::cerr << name << ": error extracting prop_header: " << e << std::endl;
547  QDP_abort(1);
548  }
549 
550 
551  // Close the namelist output file XMLDAT
552  pop(xml_out); // Prop3pt
553 
554  snoop.stop();
555  QDPIO::cout << name << ": total time = "
556  << snoop.getTimeInSeconds()
557  << " secs" << std::endl;
558 
559  QDPIO::cout << name << ": ran successfully" << std::endl;
560 
561  END_CODE();
562  }
563  } // namespace InlineProp3ptEnv
564 }// namespace chroma
Inline measurement factory.
Baryon spin and projector matrices.
All baryon operators.
Factory for producing baryon operators.
Inline measurement of stochastic baryon operators.
void operator()(const unsigned long update_no, XMLWriter &xml_out)
Do the measurement.
void func(const unsigned long update_no, XMLWriter &xml_out)
Do the measurement.
Fourier transform phase factor support.
Definition: sftmom.h:35
const Set & getSet() const
The set to be used in sumMulti.
Definition: sftmom.h:57
static T & Instance()
Definition: singleton.h:432
Random Z(N) source construction using dilution.
All dilution scheme factories.
Factory for dilution schemes.
Parallel transport a lattice field.
void PropToFerm(const LatticePropagatorF &b, LatticeFermionF &a, int color_index, int spin_index)
Extract a LatticeFermion from a LatticePropagator.
Definition: transf.cc:226
void FermToProp(const LatticeFermionF &a, LatticePropagatorF &b, int color_index, int spin_index)
Insert a LatticeFermion into a LatticePropagator.
Definition: transf.cc:98
std::string uniqueId()
Return a unique id.
Definition: unique_id.cc:18
void proginfo(XMLWriter &xml)
Print out basic information about this program.
Definition: proginfo.cc:24
std::string makeXMLFileName(std::string xml_file, unsigned long update_no)
Return a xml file name for inline measurements.
multi1d< GroupXML_t > readXMLArrayGroup(XMLReader &xml_in, const std::string &path, const std::string &type_name)
Read group and return as a std::string.
Class for counted reference semantics.
std::map< std::string, SinkPropContainer_t > prop
ForwardProp_t prop_header
Inline measurement of stochastic 3pt functions.
unsigned n
Definition: ldumul_w.cc:36
void savern(int iseed[4])
Definition: make_seeds.cc:46
Make xml file writer.
int t
Definition: meslate.cc:37
Nd
Definition: meslate.cc:74
Named object function std::map.
static bool registered
Local registration flag.
void read(XMLReader &xml, const std::string &path, Params::Operator_t &op)
void write(XMLWriter &xml, const std::string &path, const Params::Operator_t &op)
bool registerAll()
Register all the factories.
QDPSubTypeTrait< typename BinaryReturn< C1, C2, FnLocalInnerProduct >::Type_t >::Type_t localInnerProduct(const QDPSubType< T1, C1 > &l, const QDPType< T2, C2 > &r)
const int N_quarks
Number of quarks to be used in this construction.
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
static multi1d< LatticeColorMatrix > u
push(xml_out,"Condensates")
Double c
Definition: invbicg.cc:108
int i
Definition: pbg5p_w.cc:55
pop(xml_out)
void MesPlq(const multi1d< LatticeColorMatrixF3 > &u, multi2d< Double > &plane_plaq, Double &link)
Definition: mesplq.cc:83
START_CODE()
Double zero
Definition: invbicg.cc:106
multi1d< LatticeFermion > s(Ncb)
int count
Definition: octave.h:14
FloatingPoint< double > Double
Definition: gtest.h:7351
::std::string string
Definition: gtest.h:1979
Print out basic info about this program.
Double sum
Definition: qtopcor.cc:37
All quark smearing constructors.
Factory for producing quark smearing objects.
Fourier transform phase factor support.
All make sink constructors.
Factory for producing quark prop sinks.
All source smearing.
Factory for producing quark smearing objects.
Propagator parameters.
Definition: qprop_io.h:75
Hold group xml and type id.
struct Chroma::InlineProp3ptEnv::Params::Param_t param
void write(XMLWriter &xml_out, const std::string &path)
struct Chroma::InlineProp3ptEnv::Params::NamedObject_t named_obj
Propagator source construction parameters.
Definition: qprop_io.h:27
Generate a unique id.
Volume source of Z(N) noise.