CHROMA
inline_disco_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 #include "util/ferm/key_val_db.h"
32 #include <vector>
33 #include <map>
34 
35 namespace Chroma{
36  namespace InlineDiscoEnv{
37  namespace{
38  AbsInlineMeasurement* createMeasurement(XMLReader& xml_in,
39  const std::string& path)
40  {
41  return new InlineMeas(Params(xml_in, path));
42  }
43 
44  //! Local registration flag
45  bool registered = false;
46  }
47 
48  const std::string name = "DISCO";
49 
50  //! Register all the factories
51  bool registerAll()
52  {
53  bool success = true;
54  if (! registered)
55  {
56  //success &= BaryonOperatorEnv::registerAll();
57  success &= TheInlineMeasurementFactory::Instance().registerObject(name, createMeasurement);
58  registered = true;
59  }
60  return success;
61  }
62 
63  // Reader for input parameters
64  void read(XMLReader& xml, const std::string& path, Params::Param_t& param){
65  XMLReader paramtop(xml, path);
66 
67  int version;
68  read(paramtop, "version", version);
69 
70  switch (version)
71  {
72  case 1:
73  /************************************************************/
74  read(paramtop,"max_path_length",param.max_path_length);
75  read(paramtop,"p2_max",param.p2_max);
76  read(paramtop,"mass_label",param.mass_label);
77  param.chi = readXMLArrayGroup(paramtop, "Quarks", "DilutionType");
78 
79  break;
80 
81  default :
82  /**************************************************************/
83 
84  QDPIO::cerr << "Input parameter version " << version << " unsupported." << std::endl;
85  QDP_abort(1);
86  }
87 
88 
89  }
90 
91 
92  // Writter for input parameters
93  void write(XMLWriter& xml, const std::string& path, const Params::Param_t& param){
94  push(xml, path);
95 
96  int version = 1;
97 
98  write(xml, "version", version);
99 
100  write(xml,"max_path_length",param.max_path_length);
101  write(xml,"p2_max",param.p2_max);
102  write(xml,"mass_label",param.mass_label);
103 
104  push(xml,"Quarks");
105  for( int t(0);t<param.chi.size();t++){
106  push(xml,"elem");
107  xml<<param.chi[t].xml;
108  pop(xml);
109  }
110  pop(xml);
111 
112  pop(xml); // final pop
113  }
114 
115 
116  //! Gauge field parameters
117  void read(XMLReader& xml, const std::string& path, Params::NamedObject_t& input)
118  {
119  XMLReader inputtop(xml, path);
120 
121  read(inputtop, "gauge_id", input.gauge_id);
122  read(inputtop, "op_db_file", input.op_db_file);
123  }
124 
125  //! Gauge field parameters
126  void write(XMLWriter& xml, const std::string& path, const Params::NamedObject_t& input){
127  push(xml, path);
128 
129  write(xml, "gauge_id", input.gauge_id);
130  write(xml, "op_db_file", input.op_db_file);
131  pop(xml);
132  }
133 
134 
135  // Param stuff
137  frequency = 0;
138  }
139 
140  Params::Params(XMLReader& xml_in, const std::string& path)
141  {
142  try
143  {
144  XMLReader paramtop(xml_in, path);
145 
146  if (paramtop.count("Frequency") == 1)
147  read(paramtop, "Frequency", frequency);
148  else
149  frequency = 1;
150 
151  // Read program parameters
152  read(paramtop, "Param", param);
153 
154  // Read in the output propagator/source configuration info
155  read(paramtop, "NamedObject", named_obj);
156 
157  // Possible alternate XML file pattern
158  if (paramtop.count("xml_file") != 0)
159  {
160  read(paramtop, "xml_file", xml_file);
161  }
162  }
163  catch(const std::string& e)
164  {
165  QDPIO::cerr << __func__ << ": Caught Exception reading XML: " << e << std::endl;
166  QDP_abort(1);
167  }
168  }
169 
170 
171  void Params::write(XMLWriter& xml_out, const std::string& path)
172  {
173  push(xml_out, path);
174 
175  // Parameters for source construction
176  InlineDiscoEnv::write(xml_out, "Param", param);
177 
178  // Write out the output propagator/source configuration info
179  InlineDiscoEnv::write(xml_out, "NamedObject", named_obj);
180 
181  pop(xml_out);
182  }
183 
184 
185  //! Meson operator
187  {
188  unsigned short int t_slice ; /*!< Meson operator time slice */
189  multi1d<short int> disp ; /*!< Displacement dirs of quark (right)*/
190  multi1d<short int> mom ; /*!< D-1 momentum of this operator */
191 
193  mom.resize(Nd-1);
194  }
195  };
196 
197  bool operator<(const KeyOperator_t& a, const KeyOperator_t& b){
198  return ((a.t_slice<b.t_slice)||(a.mom<b.mom)||(a.disp<b.disp));
199  }
200 
201  std::ostream& operator<<(std::ostream& os, const KeyOperator_t& d)
202  {
203  os << "KeyOperator_t:"
204  << " t_slice = " << d.t_slice
205  << ", disp = ";
206  for (int i=0; i<d.disp.size();i++){
207  os << d.disp[i] << " " ;
208  }
209  os << ", mom = ";
210  for (int i=0; i<d.mom.size();i++){
211  os << d.mom[i] << " " ;
212  }
213  os << std::endl;
214 
215  return os;
216  }
218  public:
219  multi1d<ComplexD> op ;
220  ValOperator_t(){op.resize(Ns*Ns);} // Here go the 16 gamma matrices
222  } ;
223 
224  //-------------------------------------------------------------------------
225  //! stream IO
226  std::ostream& operator<<(std::ostream& os, const ValOperator_t& d)
227  {
228  os << "ValOperator_t:\n";
229  for (int i=0; i<d.op.size();i++){
230  os <<" gamma["<<i<<"] = "<< d.op[i] << std::endl ;
231  }
232 
233  return os;
234  }
235 
236  struct KeyVal_t{
239  };
240 
241  //! KeyOperator reader
242  void read(BinaryReader& bin, KeyOperator_t& d){
243  read(bin,d.t_slice);
244  unsigned short int n ;
245  read(bin,n);
246  d.disp.resize(n);
247  read(bin,d.disp);
248  d.mom.resize(Nd-1) ;
249  read(bin,d.mom);
250  }
251  //! KeyOperator writer
252  void write(BinaryWriter& bin, const KeyOperator_t& d){
253  write(bin,d.t_slice);
254  unsigned short int n ;
255  n = d.disp.size();
256  write(bin,n);
257  write(bin,d.disp);
258  write(bin,d.mom);
259  }
260 
261  //! ValOperator reader
262  void read(BinaryReader& bin, ValOperator_t& d){
263  d.op.resize(Ns*Ns);
264  read(bin,d.op);
265  }
266  //! ValOperator writer
267  void write(BinaryWriter& bin, const ValOperator_t& d){
268  write(bin,d.op);
269  }
270 
271  namespace{
272  StandardOutputStream& operator<<(StandardOutputStream& os, const multi1d<short int>& d){
273  if (d.size() > 0){
274  os << d[0];
275  for(int i=1; i < d.size(); ++i)
276  os << " " << d[i];
277  }
278  return os;
279  }
280  }
281 
282  void do_disco(std::map< KeyOperator_t, ValOperator_t >& db,
283  const LatticeFermion& qbar,
284  const LatticeFermion& q,
285  const SftMom& p,
286  const int& t,
287  const multi1d<short int>& path,
288  const int& max_path_length ){
289  QDPIO::cout<<" Computing Operator with path length "<<path.size()
290  <<" on timeslice "<<t<<". Path: "<<path <<std::endl;
291 
292  ValOperator_t val ;
293  KeyOperator_t key ;
294  std::pair<KeyOperator_t, ValOperator_t> kv ;
295  kv.first.t_slice = t ;
296  if(path.size()==0){
297  kv.first.disp.resize(1);
298  kv.first.disp[0] = 0 ;
299  }
300  else
301  kv.first.disp = path ;
302 
303  multi1d< multi1d<ComplexD> > foo(p.numMom()) ;
304  for (int m(0); m < p.numMom(); m++)
305  foo[m].resize(Ns*Ns);
306  for(int g(0);g<Ns*Ns;g++){
307  LatticeComplex cc = localInnerProduct(qbar,Gamma(g)*q);
308  for (int m(0); m < p.numMom(); m++){
309  foo[m][g] = sum(p[m]*cc,p.getSet()[t]) ;
310  }
311  }
312  for (int m(0); m < p.numMom(); m++){
313  for(int i(0);i<(Nd-1);i++)
314  kv.first.mom[i] = p.numToMom(m)[i] ;
315 
316  kv.second.op = foo[m];
317  std::pair<std::map< KeyOperator_t, ValOperator_t >::iterator, bool> itbo;
318 
319  itbo = db.insert(kv);
320  if( itbo.second ){
321  QDPIO::cout<<"Inserting new entry in std::map\n";
322  }
323  else{ // if insert fails, key already exists, so add result
324  std::cout<<"Key = "<<kv.first<<std::endl;
325  QDPIO::cout<<"Adding result to value already there"<<std::endl;
326  for(int i(0);i<kv.second.op.size();i++){
327  itbo.first->second.op[i] += kv.second.op[i] ;
328  }
329  }
330  }
331 
332  if(path.size()<max_path_length){
333  QDPIO::cout<<" attempt to add new path. "
334  <<" current path length is : "<<path.size();
335  multi1d<short int> new_path(path.size()+1);
336  QDPIO::cout<<" new path length is : "<<new_path.size()<<std::endl;
337  for(int i(0);i<path.size();i++)
338  new_path[i] = path[i] ;
339  for(int sign(-1);sign<2;sign+=2)
340  for(int mu(0);mu<Nd;mu++){
341  new_path[path.size()]= sign*(mu+1) ;
342  //skip back tracking
343  bool back_track=false ;
344  if(path.size()>0)
345  if(path[path.size()-1] == -new_path[path.size()])
346  back_track=true;
347  if(!back_track){
348  QDPIO::cout<<" Added path: "<<new_path<<std::endl;
349  LatticeFermion q_mu ;
350  if(sign>0)
351  q_mu = shift(q, FORWARD, mu);
352  else
353  q_mu = shift(q, BACKWARD, mu);
354 
355  do_disco(db, qbar, q_mu, p, t, new_path, max_path_length);
356  } // skip backtracking
357  } // mu
358  }
359 
360  }// do_disco
361 
362 
363  //--------------------------------------------------------------
364  // Function call
365  // void
366  //InlineDisco::operator()(unsigned long update_no,
367  // XMLWriter& xml_out)
368  void InlineMeas::operator()(unsigned long update_no,
369  XMLWriter& xml_out)
370  {
371  // If xml file not empty, then use alternate
372  if (params.xml_file != ""){
373  std::string xml_file = makeXMLFileName(params.xml_file, update_no);
374 
375  push(xml_out, "propagator_3pt");
376  write(xml_out, "update_no", update_no);
377  write(xml_out, "xml_file", xml_file);
378  pop(xml_out);
379 
380  XMLFileWriter xml(xml_file);
381  func(update_no, xml);
382  }
383  else{
384  func(update_no, xml_out);
385  }
386  }
387 
388 
389  // Function call
390  //void
391  //InlineDisco::func(unsigned long update_no,
392  // XMLWriter& xml_out)
393  void InlineMeas::func(unsigned long update_no,
394  XMLWriter& xml_out)
395  {
396  START_CODE();
397 
398  StopWatch snoop;
399  snoop.reset();
400  snoop.start();
401 
402  // Test and grab a reference to the gauge field
403  XMLBufferWriter gauge_xml;
404  try
405  {
406  TheNamedObjMap::Instance().getData< multi1d<LatticeColorMatrix> >(params.named_obj.gauge_id);
407  TheNamedObjMap::Instance().get(params.named_obj.gauge_id).getRecordXML(gauge_xml);
408  }
409  catch( std::bad_cast )
410  {
411  QDPIO::cerr << InlineDiscoEnv::name << ": caught dynamic cast error"
412  << std::endl;
413  QDP_abort(1);
414  }
415  catch (const std::string& e)
416  {
417  QDPIO::cerr << name << ": std::map call failed: " << e
418  << std::endl;
419  QDP_abort(1);
420  }
421  const multi1d<LatticeColorMatrix>& u =
422  TheNamedObjMap::Instance().getData< multi1d<LatticeColorMatrix> >(params.named_obj.gauge_id);
423 
424  push(xml_out, "disco");
425  write(xml_out, "update_no", update_no);
426 
427  QDPIO::cout << name << ": Disconnected diagrams" << std::endl;
428 
429  proginfo(xml_out); // Print out basic program info
430 
431  // Write out the input
432  params.write(xml_out, "Input");
433 
434  // Write out the config info
435  write(xml_out, "Config_info", gauge_xml);
436 
437  push(xml_out, "Output_version");
438  write(xml_out, "out_version", 1);
439  pop(xml_out);
440 
441  // First calculate some gauge invariant observables just for info.
442  // This is really cheap.
443  MesPlq(xml_out, "Observables", u);
444 
445  // Save current seed
446  Seed ran_seed;
447  QDP::RNG::savern(ran_seed);
448 
449  //
450  // Read the source and solutions
451  //
452  StopWatch swatch;
453  swatch.reset();
454  swatch.start();
455 
456  int N_quarks = params.param.chi.size() ;
457 
458  multi1d< Handle< DilutionScheme<LatticeFermion> > > quarks(N_quarks);
459 
460  try{
461  // Loop over quark dilutions
462  for(int n(0); n < params.param.chi.size(); ++n){
463  const GroupXML_t& dil_xml = params.param.chi[n];
464  std::istringstream xml_d(dil_xml.xml);
465  XMLReader diltop(xml_d);
466  QDPIO::cout << "Dilution type = " << dil_xml.id << std::endl;
467  quarks[n] =
468  TheFermDilutionSchemeFactory::Instance().createObject(dil_xml.id,
469  diltop,
470  dil_xml.path);
471  }
472  }
473  catch(const std::string& e){
474  QDPIO::cerr << name << ": Caught Exception constructing dilution scheme: " << e << std::endl;
475  QDP_abort(1);
476  }
477 
478 
479  //-------------------------------------------------------------------
480  //Sanity checks
481 
482  //Another Sanity check, the three quarks must all be
483  //inverted on the same cfg
484  for (int n = 1 ; n < N_quarks ; ++n){
485  if (quarks[0]->getCfgInfo() != quarks[n]->getCfgInfo()){
486  QDPIO::cerr << name << " : Quarks do not contain the same cfg info";
487  QDPIO::cerr << ", quark "<< n << std::endl;
488  QDP_abort(1);
489  }
490  }
491 
492 
493  //Also ensure that the cfg on which the inversions were performed
494  //is the same as the cfg that we are using
495  {
496  std::string cfgInfo;
497 
498  //Really ugly way of doing this(Robert Heeeelp!!)
499  XMLBufferWriter top;
500  write(top, "Config_info", gauge_xml);
501  XMLReader from(top);
502  XMLReader from2(from, "/Config_info");
503  std::ostringstream os;
504  from2.print(os);
505 
506  cfgInfo = os.str();
507 
508  if (cfgInfo != quarks[0]->getCfgInfo()){
509  QDPIO::cerr << name << " : Quarks do not contain the same";
510  QDPIO::cerr << " cfg info as the gauge field." ;
511  QDPIO::cerr << "gauge: XX"<<cfgInfo<<"XX quarks: XX" ;
512  QDPIO::cerr << quarks[0]->getCfgInfo()<<"XX"<< std::endl;
513  QDP_abort(1);
514  }
515  }
516 
517 
518  int decay_dir = quarks[0]->getDecayDir();
519  //
520  // Initialize the slow Fourier transform phases
521  //
522  //SftMom phases(params.param.mom2_max, false, decay_dir);
523  SftMom phases(params.param.p2_max, false, decay_dir);
524 
525  // Another sanity check. The seeds of all the quarks must be different
526  // and thier decay directions must be the same
527  for(int n = 1 ; n < quarks.size(); ++n){
528  if(toBool(quarks[n]->getSeed()==quarks[0]->getSeed())){
529  QDPIO::cerr << name << ": error, quark seeds are the same" << std::endl;
530  QDP_abort(1);
531  }
532 
533  if(toBool(quarks[n]->getDecayDir()!=quarks[0]->getDecayDir())){
534  QDPIO::cerr<<name<< ": error, quark decay dirs do not match" <<std::endl;
535  QDP_abort(1);
536  }
537  }
538 
539  std::map< KeyOperator_t, ValOperator_t > data ;
540 
541  for(int n(0);n<quarks.size();n++){
542  for (int it(0) ; it < quarks[n]->getNumTimeSlices() ; ++it){
543  int t = quarks[n]->getT0(it) ;
544  QDPIO::cout<<" Doing quark: "<<n <<std::endl ;
545  QDPIO::cout<<" quark: "<<n <<" has "<<quarks[n]->getDilSize(it);
546  QDPIO::cout<<" dilutions on time slice "<<t<<std::endl ;
547  for(int i = 0 ; i < quarks[n]->getDilSize(it) ; ++i){
548  QDPIO::cout<<" Doing dilution : "<<i<<std::endl ;
549  multi1d<short int> d ;
550  LatticeFermion qbar = quarks[n]->dilutedSource(it,i);
551  LatticeFermion q = quarks[n]->dilutedSolution(it,i);
552  QDPIO::cout<<" Starting recursion "<<std::endl ;
553  do_disco(data, qbar, q, phases, t, d, params.param.max_path_length);
554  QDPIO::cout<<" done with recursion! "
555  <<" The length of the path is: "<<d.size()<<std::endl ;
556  }
557  QDPIO::cout<<" Done with dilutions for quark: "<<n <<std::endl ;
558  }
559  }
560 
561  // DB storage
562  BinaryStoreDB<SerialDBKey<KeyOperator_t>,SerialDBData<ValOperator_t> > qdp_db;
563 
564  // Open the file, and write the meta-data and the binary for this operator
565  {
566  XMLBufferWriter file_xml;
567 
568  push(file_xml, "DBMetaData");
569  write(file_xml, "id", std::string("eigElemOp"));
570  write(file_xml, "lattSize", QDP::Layout::lattSize());
571  write(file_xml, "decay_dir", decay_dir);
572  write(file_xml, "Params", params.param);
573  write(file_xml, "Config_info", gauge_xml);
574  pop(file_xml);
575 
576  std::string file_str(file_xml.str());
577  qdp_db.setMaxUserInfoLen(file_str.size());
578 
579  qdp_db.open(params.named_obj.op_db_file, O_RDWR | O_CREAT, 0664);
580 
581  qdp_db.insertUserdata(file_str);
582  }
583 
584  // Write the data
587  std::map< KeyOperator_t, ValOperator_t >::iterator it;
588  for(it=data.begin();it!=data.end();it++){
589  key.key() = it->first ;
590  val.data().op.resize(it->second.op.size()) ;
591  // normalize to number of quarks
592  for(int i(0);i<it->second.op.size();i++)
593  val.data().op[i] = it->second.op[i]/toDouble(quarks.size());
594  qdp_db.insert(key,val) ;
595  }
596 
597  // Close the namelist output file XMLDAT
598  pop(xml_out); // Disco
599 
600  snoop.stop();
601  QDPIO::cout << name << ": total time = "
602  << snoop.getTimeInSeconds()
603  << " secs" << std::endl;
604 
605  QDPIO::cout << name << ": ran successfully" << std::endl;
606 
607  END_CODE();
608  }
609  } // namespace InlineDiscoEnv
610 }// 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.
Serializable value harness.
Definition: key_val_db.h:69
D & data()
Setter.
Definition: key_val_db.h:78
Serializable key harness.
Definition: key_val_db.h:21
K & key()
Setter.
Definition: key_val_db.h:30
Fourier transform phase factor support.
Definition: sftmom.h:35
static T & Instance()
Definition: singleton.h:432
int mu
Definition: cool.cc:24
Random Z(N) source construction using dilution.
All dilution scheme factories.
Factory for dilution schemes.
Parallel transport a lattice field.
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.
Inline measurement of stochastic 3pt functions.
Key and values for DB.
unsigned n
Definition: ldumul_w.cc:36
static int m[4]
Definition: make_seeds.cc:16
void savern(int iseed[4])
Definition: make_seeds.cc:46
Make xml file writer.
int it
Definition: meslate.cc:33
int t
Definition: meslate.cc:37
Nd
Definition: meslate.cc:74
Double q
Definition: mesq.cc:17
Named object function std::map.
static bool registered
Local registration flag.
void read(XMLReader &xml, const std::string &path, Params::Param_t &param)
void write(XMLWriter &xml, const std::string &path, const Params::Param_t &param)
const std::string name
std::ostream & operator<<(std::ostream &os, const KeyOperator_t &d)
bool registerAll()
Register all the factories.
bool operator<(const KeyOperator_t &a, const KeyOperator_t &b)
void do_disco(std::map< KeyOperator_t, ValOperator_t > &db, const LatticeFermion &qbar, const LatticeFermion &q, const SftMom &p, const int &t, const multi1d< short int > &path, const int &max_path_length)
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")
int i
Definition: pbg5p_w.cc:55
Complex a
Definition: invbicg.cc:95
pop(xml_out)
void MesPlq(const multi1d< LatticeColorMatrixF3 > &u, multi2d< Double > &plane_plaq, Double &link)
Definition: mesplq.cc:83
DComplex d
Definition: invbicg.cc:99
START_CODE()
Complex b
Definition: invbicg.cc:96
::std::string string
Definition: gtest.h:1979
#define FORWARD
Definition: primitives.h:82
#define BACKWARD
Definition: primitives.h:83
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.
Hold group xml and type id.
SerialDBData< ValOperator_t > v
SerialDBKey< KeyOperator_t > k
struct Chroma::InlineDiscoEnv::Params::NamedObject_t named_obj
struct Chroma::InlineDiscoEnv::Params::Param_t param
void write(XMLWriter &xml_out, const std::string &path)
Generate a unique id.
Volume source of Z(N) noise.