MetaData.h

Go to the documentation of this file.

/*************************************************************************************
 
    Grid physics library, www.github.com/paboyle/Grid 
 
    Source file: ./lib/parallelIO/Metadata.h
 
    Copyright (C) 2015, 2026
 
    Author: Peter Boyle <paboyle@ph.ed.ac.uk>
    Author: Jamie Hudspith <renwick.james.hudspth@gmail.com>
    Author: Gaurav Ray <gaurav.sinharay@swansea.ac.uk>
 
    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.
 
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
 
    You should have received a copy of the GNU General Public License along
    with this program; if not, write to the Free Software Foundation, Inc.,
    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
    See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/*  END LEGAL */
 
#include <algorithm>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <map>
#include <unistd.h>
#include <sys/utsname.h>
#include <pwd.h>
 
NAMESPACE_BEGIN(Grid);

// Precision mapping
template<class vobj> static std::string getFormatString (void)
{
  std::string format;
  typedef typename getPrecision<vobj>::real_scalar_type stype;
  if ( sizeof(stype) == sizeof(float) ) {
    format = std::string("IEEE32BIG");
  }
  if ( sizeof(stype) == sizeof(double) ) {
    format = std::string("IEEE64BIG");
  }
  return format;
};

  // header specification/interpretation
    class FieldNormMetaData : Serializable {
    public:
      GRID_SERIALIZABLE_CLASS_MEMBERS(FieldNormMetaData, double, norm2);
    };
    class FieldMetaData : Serializable {
    public:
 
      GRID_SERIALIZABLE_CLASS_MEMBERS(FieldMetaData,
                      int, nd,
                      std::vector<int>, dimension,
                      std::vector<std::string>, boundary,
                      int, data_start,
                      std::string, hdr_version,
                      std::string, storage_format,
                      double, link_trace,
                      double, plaquette,
                      uint32_t, checksum,
                      uint32_t, scidac_checksuma,
                      uint32_t, scidac_checksumb,
                      unsigned int, sequence_number,
                      std::string, data_type,
                      std::string, ensemble_id,
                      std::string, ensemble_label,
                      std::string, ildg_lfn,
                      std::string, creator,
                      std::string, creator_hardware,
                      std::string, creation_date,
                      std::string, archive_date,
                      std::string, floating_point);
      // WARNING: non-initialised values might lead to twisted parallel IO
      // issues, std::string are fine because they initliase to size 0
      // as per C++ standard.
      FieldMetaData(void) 
      : nd(4), dimension(4,0), boundary(4, ""), data_start(0),
      link_trace(0.), plaquette(0.), checksum(0),
      scidac_checksuma(0), scidac_checksumb(0), sequence_number(0)
      {}
  };
 
// PB disable using namespace - this is a header and forces namesapce visibility for all 
// including files
//using namespace Grid;

// Bit and Physical Checksumming and QA of data
inline void GridMetaData(GridBase *grid,FieldMetaData &header)
{
  int nd = grid->_ndimension;
  header.nd = nd;
  header.dimension.resize(nd);
  header.boundary.resize(nd);
  header.data_start = 0;
  for(int d=0;d<nd;d++) {
    header.dimension[d] = grid->_fdimensions[d];
  }
  for(int d=0;d<nd;d++) {
    header.boundary[d] = std::string("PERIODIC");
  }
}
 
inline void MachineCharacteristics(FieldMetaData &header)
{
  // Who
  struct passwd *pw = getpwuid (getuid());
  if (pw) header.creator = std::string(pw->pw_name); 
 
  // When
  std::time_t t = std::time(nullptr);
  std::tm tm_ = *std::localtime(&t);
  std::ostringstream oss; 
  oss << std::put_time(&tm_, "%c %Z");
  header.creation_date = oss.str();
  header.archive_date  = header.creation_date;
 
  // What
  struct utsname name;  uname(&name);
  header.creator_hardware = std::string(name.nodename)+"-";
  header.creator_hardware+= std::string(name.machine)+"-";
  header.creator_hardware+= std::string(name.sysname)+"-";
  header.creator_hardware+= std::string(name.release);
}
 
#define dump_meta_data(field, s)                    \
  s << "BEGIN_HEADER"      << std::endl;                \
  s << "HDR_VERSION = "    << field.hdr_version    << std::endl;    \
  s << "DATATYPE = "       << field.data_type      << std::endl;    \
  s << "STORAGE_FORMAT = " << field.storage_format << std::endl;    \
  for(int i=0;i<4;i++){                         \
    s << "DIMENSION_" << i+1 << " = " << field.dimension[i] << std::endl ; \
  }                                 \
  s << "LINK_TRACE = " << std::setprecision(10) << field.link_trace << std::endl; \
  s << "PLAQUETTE  = " << std::setprecision(10) << field.plaquette  << std::endl; \
  for(int i=0;i<4;i++){                         \
    s << "BOUNDARY_"<<i+1<<" = " << field.boundary[i] << std::endl; \
  }                                 \
                                    \
  s << "CHECKSUM = "<< std::hex << std::setw(10) << field.checksum << std::dec<<std::endl; \
  s << "SCIDAC_CHECKSUMA = "<< std::hex << std::setw(10) << field.scidac_checksuma << std::dec<<std::endl; \
  s << "SCIDAC_CHECKSUMB = "<< std::hex << std::setw(10) << field.scidac_checksumb << std::dec<<std::endl; \
  s << "ENSEMBLE_ID = "     << field.ensemble_id      << std::endl; \
  s << "ENSEMBLE_LABEL = "  << field.ensemble_label   << std::endl; \
  s << "SEQUENCE_NUMBER = " << field.sequence_number  << std::endl; \
  s << "CREATOR = "         << field.creator          << std::endl; \
  s << "CREATOR_HARDWARE = "<< field.creator_hardware << std::endl; \
  s << "CREATION_DATE = "   << field.creation_date    << std::endl; \
  s << "ARCHIVE_DATE = "    << field.archive_date     << std::endl; \
  s << "FLOATING_POINT = "  << field.floating_point   << std::endl; \
  s << "END_HEADER"         << std::endl;
 
template<class vobj> inline void PrepareMetaData(Lattice<vobj> & field, FieldMetaData &header)
{
  GridBase *grid = field.Grid();
  std::string format = getFormatString<vobj>();
  header.floating_point = format;
  header.checksum = 0x0; // Nersc checksum unused in ILDG, Scidac
  GridMetaData(grid,header); 
  MachineCharacteristics(header);
}
template<class Impl>
class GaugeStatistics
{
public:
  void operator()(Lattice<vLorentzColourMatrixD> & data,FieldMetaData &header)
  {
    header.link_trace = WilsonLoops<Impl>::linkTrace(data);
    header.plaquette  = WilsonLoops<Impl>::avgPlaquette(data);
  }
};
typedef GaugeStatistics<PeriodicGimplD> PeriodicGaugeStatistics;
typedef GaugeStatistics<ConjugateGimplD> ConjugateGaugeStatistics;
template<> inline void PrepareMetaData<vLorentzColourMatrixD>(Lattice<vLorentzColourMatrixD> & field, FieldMetaData &header)
{
  GridBase *grid = field.Grid();
  std::string format = getFormatString<vLorentzColourMatrixD>();
  header.floating_point = format;
  header.checksum = 0x0; // Nersc checksum unused in ILDG, Scidac
  GridMetaData(grid,header); 
  MachineCharacteristics(header);
}

// Utilities ; these are QCD aware
inline void reconstruct3(LorentzColourMatrix & cm)
{
  assert( Nc < 4 && Nc > 1 ) ;
  for(int mu=0;mu<Nd;mu++){
    #if Nc == 2
      cm(mu)()(1,0) = -adj(cm(mu)()(0,y)) ;
      cm(mu)()(1,1) =  adj(cm(mu)()(0,x)) ;
    #else
      const int x=0 , y=1 , z=2 ; // a little disingenuous labelling
      cm(mu)()(2,x) = adj(cm(mu)()(0,y)*cm(mu)()(1,z)-cm(mu)()(0,z)*cm(mu)()(1,y)); //x= yz-zy
      cm(mu)()(2,y) = adj(cm(mu)()(0,z)*cm(mu)()(1,x)-cm(mu)()(0,x)*cm(mu)()(1,z)); //y= zx-xz
      cm(mu)()(2,z) = adj(cm(mu)()(0,x)*cm(mu)()(1,y)-cm(mu)()(0,y)*cm(mu)()(1,x)); //z= xy-yx
    #endif
  }
}
 
// the elements of the final row of an SU(N) matrix
// can be obtained from the determinants of the N N-1 by N-1 matrices
// formed by deleting each column and the Nth row of the SU(N) matrix.
// 1) create a ColourSubMatrix object to store the N-1 dim matrix.
// 2) peek the colour matrix in a particular direction. 
// 3) fill the ColourSubMatrix object from the peeked matrix.
// 4) take the Determinant and fill the corresponding element
// 5) repeat for each Lorentz index 
inline void reconstructSU(LorentzColourMatrix &cm)
{
  using ColourSubMatrix = iScalar<iScalar<iMatrix<Complex, Nc-1> > > ;
 
  ColourSubMatrix Su; // for the Nc-1 by Nc-1 matrix
  ColourMatrix SU;   // for the Nc-1 by N matrix read from disk
 
  for(int mu=0; mu<Nd; mu++) {
    Su = Zero();
    SU = peekIndex<LorentzIndex>(cm,mu); 
    for(int k=0; k<Nc; k++) {
      for(int i=0; i<Nc-1; i++) {
        for(int j=0; j<Nc-1; j++) {
          int J = (j<k) ? j : j+1; // for correct indexing of columns in SU
          Su()()(i,j) = SU()()(i,J);
        }
      }
      SU()()(Nc-1,k) = std::pow(-1,k+Nc-1) * adj(Determinant(Su));
    } 
    pokeIndex<LorentzIndex>(cm,SU,mu);
  }
}
 
// this function determines if a sequence of integers {i0...ik}
// is an even or odd parity permutation. even/odd if the number of 
// inversions needed to get back to the lexicographic 1st sequence is
// even or odd. ex: {1,2,0} --> {1,0,2} --> {0,1,2} implies {1,2,0} is even. 
inline bool is_perm_even(std::vector<int> &v) {
 
    int n = v.size();
    std::vector<int> a(n,0);
    int c = 0;
 
    for(int j=0; j<n; j++) {
        if(a[j]==0) {
            c += 1;
            a[j] = 1;
            int i = j;
            while(v[i]!=j) {
                i = v[i];
                a[i] = 1;
            }
        }
    }
 
    if ((n-c)%2==0) {
      return true;
    } else {
      return false;
    }
}

//
//  this function follows the prescription laid out by the
//  ildg spec, forming a sum of products in lexicographic order.
//
// SU[N-1][k] = sum ( SU[0][i0].SU[1][i1]...SU[N-2][iN-2] )*
//         {i0...iN-2!=k}
//
// see appendix A.2 of
// https://www-zeuthen.desy.de/apewww/ILDG/specifications/ildg-file-format-1.2.pdf
inline void unique_reconstructSU(LorentzColourMatrix &cm)
{
 
  std::vector<int> indices(Nc);
  std::iota(indices.begin(), indices.end(), 0); // {0,1,...,Nc-1}
 
  for(int mu=0; mu<Nd; mu++) {
    for(int k=0; k<Nc; k++) {
      std::vector<int> cols = indices;
      cols.erase(cols.begin()+k);     // {0,...,k-1,k+1,...Nc-1}
      cm(mu)()(Nc-1,k) = 0;
            // construct sum of products
      // as we cycle through permutations of cols
      do
      {
        std::vector<int> perm = cols;
        perm.emplace_back(k);
        ComplexD prod(1,0);
        for(int i=0;i<Nc-1;i++) {
          prod *= cm(mu)()(i,cols[i]) ;
        }
        if( is_perm_even(perm) ) { cm(mu)()(Nc-1,k) += conjugate(prod); }
        else                                     { cm(mu)()(Nc-1,k) -= conjugate(prod); }
      }
      while (std::next_permutation(cols.begin(), cols.end())); 
    }
  }
}

//  this function takes a reduced format
//  Sp(2N) field with N rows and 2N columns 
//  and reconstructs the full 2Nx2N matrix.
//  the full matrix has block structure: 
//      | A  B |
//      |-B* A*|
//  where A and B are NxN matrices.
//  see appendix A.2 of
// https://www-zeuthen.desy.de/apewww/ILDG/specifications/ildg-file-format-1.2.pdf
inline void reconstructSp(LorentzColourMatrix & cm) 
{
  assert( Nc%2==0 ); 
  int N = Nc/2;
  for(int mu=0;mu<Nd;mu++){
    for(int i=0;i<N;i++){
      for(int j=0;j<N;j++){
        cm(mu)()(i+N,j) = -adj( cm(mu)()(i,j+N) ); //B to -B*
        cm(mu)()(i+N,j+N) = adj( cm(mu)()(i,j) );   //A to A*
      }
    }
  }
}

// Some data types for intermediate storage
template<typename vtype> using iLorentzColour2x3 = iVector<iVector<iVector<vtype, Nc>, Nc-1>, Nd >;
typedef iLorentzColour2x3<Complex>  LorentzColour2x3;
typedef iLorentzColour2x3<ComplexF> LorentzColour2x3F;
typedef iLorentzColour2x3<ComplexD> LorentzColour2x3D;
 
// intermediate types for Sp(2N) fields
template<typename vtype> using iLorentzColourNx2N = iVector<iVector<iVector<vtype, Nc>, Nc/2>, Nd >;
typedef iLorentzColourNx2N<ComplexF>  LorentzColourNx2NF;
typedef iLorentzColourNx2N<ComplexD>  LorentzColourNx2ND;

// Simple classes for precision conversion
template <class fobj, class sobj>
struct BinarySimpleUnmunger {
  typedef typename getPrecision<fobj>::real_scalar_type fobj_stype;
  typedef typename getPrecision<sobj>::real_scalar_type sobj_stype;
  
  void operator()(sobj &in, fobj &out) {
    // take word by word and transform accoding to the status
    fobj_stype *out_buffer = (fobj_stype *)&out;
    sobj_stype *in_buffer = (sobj_stype *)&in;
    size_t fobj_words = sizeof(out) / sizeof(fobj_stype);
    size_t sobj_words = sizeof(in) / sizeof(sobj_stype);
    assert(fobj_words == sobj_words);
    
    for (unsigned int word = 0; word < sobj_words; word++)
      out_buffer[word] = in_buffer[word];  // type conversion on the fly
    
  }
};
 
template <class fobj, class sobj>
struct BinarySimpleMunger {
  typedef typename getPrecision<fobj>::real_scalar_type fobj_stype;
  typedef typename getPrecision<sobj>::real_scalar_type sobj_stype;
 
  void operator()(fobj &in, sobj &out) {
    // take word by word and transform accoding to the status
    fobj_stype *in_buffer = (fobj_stype *)&in;
    sobj_stype *out_buffer = (sobj_stype *)&out;
    size_t fobj_words = sizeof(in) / sizeof(fobj_stype);
    size_t sobj_words = sizeof(out) / sizeof(sobj_stype);
    assert(fobj_words == sobj_words);
    
    for (unsigned int word = 0; word < sobj_words; word++)
      out_buffer[word] = in_buffer[word];  // type conversion on the fly
    
  }
};
 
 
template<class fobj,class sobj>
struct GaugeSimpleMunger{
  void operator()(fobj &in, sobj &out) {
    for (int mu = 0; mu < Nd; mu++) {
      for (int i = 0; i < Nc; i++) {
    for (int j = 0; j < Nc; j++) {
    out(mu)()(i, j) = in(mu)()(i, j);
    }}
    }
  };
};
 
template <class fobj, class sobj>
struct GaugeSimpleUnmunger {
  void operator()(sobj &in, fobj &out) {
    for (int mu = 0; mu < Nd; mu++) {
      for (int i = 0; i < Nc; i++) {
    for (int j = 0; j < Nc; j++) {
      out(mu)()(i, j) = in(mu)()(i, j);
    }}
    }
  };
};
 
template<class fobj,class sobj>
struct GaugeDoubleStoredMunger{
  void operator()(fobj &in, sobj &out) {
    for (int mu = 0; mu < Nds; mu++) {
      for (int i = 0; i < Nc; i++) {
        for (int j = 0; j < Nc; j++) {
          out(mu)()(i, j) = in(mu)()(i, j);
        }}
    }
  };
};
 
template <class fobj, class sobj>
struct GaugeDoubleStoredUnmunger {
  void operator()(sobj &in, fobj &out) {
    for (int mu = 0; mu < Nds; mu++) {
      for (int i = 0; i < Nc; i++) {
        for (int j = 0; j < Nc; j++) {
          out(mu)()(i, j) = in(mu)()(i, j);
        }}
    }
  };
};
 
template<class fobj,class sobj>
struct Gauge3x2munger{
  void operator() (fobj &in,sobj &out){
    for(int mu=0;mu<Nd;mu++){
      for(int i=0;i<Nc-1;i++){
    for(int j=0;j<Nc;j++){
      out(mu)()(i,j) = in(mu)(i)(j);
    }}
    }
    reconstruct3(out);
  }
};
 
template<class fobj,class sobj>
struct Gauge3x2unmunger{
  void operator() (sobj &in,fobj &out){
    for(int mu=0;mu<Nd;mu++){
      for(int i=0;i<Nc-1;i++){
    for(int j=0;j<Nc;j++){
      out(mu)(i)(j) = in(mu)()(i,j);
    }}
    }
  }
};
 
template<class fobj,class sobj, bool unique_su>
struct GaugeSUmunger;
// two specialisations for two ways to reconstruct
// NcxNc matrix from (Nc-1)xNc matrix
 
template<class fobj,class sobj>
struct GaugeSUmunger<fobj,sobj,false>{
  void operator() (fobj &in,sobj &out){
    for(int mu=0;mu<Nd;mu++){
      for(int i=0;i<Nc-1;i++){
        for(int j=0;j<Nc;j++){
          out(mu)()(i,j) = in(mu)(i)(j);
        }
      }
    }
  reconstructSU(out);
  }
};
 
template<class fobj,class sobj>
struct GaugeSUmunger<fobj,sobj,true>{
  void operator() (fobj &in,sobj &out){
    for(int mu=0;mu<Nd;mu++){
      for(int i=0;i<Nc-1;i++){
        for(int j=0;j<Nc;j++){
          out(mu)()(i,j) = in(mu)(i)(j);
        }
      }
    }
    unique_reconstructSU(out);
  }
};
 
 
template<class fobj,class sobj>
struct GaugeSUunmunger{
  void operator() (sobj &in,fobj &out){
    for(int mu=0;mu<Nd;mu++){
      for(int i=0;i<Nc-1;i++){
          for(int j=0;j<Nc;j++){
            out(mu)(i)(j) = in(mu)()(i,j);
          }
      }
    }
  }
};
 
template<class fobj,class sobj>
struct GaugeSpmunger{
  void operator() (fobj &in,sobj &out){
    for(int mu=0;mu<Nd;mu++){
      for(int i=0;i<Nc/2;i++){
            for(int j=0;j<Nc;j++){
            out(mu)()(i,j) = in(mu)(i)(j);
            }
        }
    }
    reconstructSp(out);
  }
};
 
// transform Sp(2N) fields into reduced Nx2N format.
template<class fobj,class sobj>
struct GaugeSpunmunger{
  void operator() (sobj &in,fobj &out){
    for(int mu=0;mu<Nd;mu++){
      for(int i=0;i<Nc/2;i++){
          for(int j=0;j<Nc;j++){
            out(mu)(i)(j) = in(mu)()(i,j);
            }
        }
    }
  }
};
 
// these are used as non-type template parameters when
// writing with GaugeUnMunger and as regular parameters
// when reading with IldgReader.readConfiguration 
enum struct FloatingPointFormat { IEEE64BIG, IEEE32BIG };
enum struct MatrixFormat { FULL, REDUCED };
// this struct is used to choose the appropriate
// unmunger (for writing) at compile time.
// there are 3 partial template specialisations,
// > no group reduction - treat SU and Sp the same
// > group reduction for SU fields
// > group reduction for Sp fields
// It also exposes the intermediate out_type for use in 
// IldgWriter.writeConfiguration
template<class vobj, class group_name, MatrixFormat m_fmt, FloatingPointFormat fp_fmt>
struct GaugeUnMunger;
 
// no group reduction
template<class vobj, class group_name, FloatingPointFormat fp_fmt>
struct GaugeUnMunger<vobj, group_name, MatrixFormat::FULL, fp_fmt>
{   
  using in_type  = typename vobj::scalar_object; 
  using out_type = typename std::tuple_element_t<static_cast<int>(fp_fmt), std::tuple<LorentzColourMatrixD,LorentzColourMatrixF>>;
  BinarySimpleUnmunger<out_type, in_type> unmunger;
 
  void operator() (in_type &in, out_type &out){
    unmunger(in,out);
  }
};
 
//template specialisation for SU
template<class vobj, FloatingPointFormat fp_fmt>
struct GaugeUnMunger<vobj, GroupName::SU, MatrixFormat::REDUCED, fp_fmt>
{
  using in_type  = typename vobj::scalar_object; 
  using tmp_type = typename std::tuple_element_t<static_cast<int>(fp_fmt), std::tuple<LorentzColourMatrixD,LorentzColourMatrixF>>;
  using out_type = typename std::tuple_element_t<static_cast<int>(fp_fmt), std::tuple<LorentzColour2x3D,LorentzColour2x3F>>;
 
  BinarySimpleUnmunger<tmp_type, in_type> binary_unmunger;
  GaugeSUunmunger<out_type, tmp_type> gauge_unmunger;
 
  void operator() (in_type &in, out_type &out){
    tmp_type tmp;
    binary_unmunger(in, tmp);
    gauge_unmunger(tmp,out);
  }
};
 
//template specialisation for Sp
template<class vobj, FloatingPointFormat fp_fmt>
struct GaugeUnMunger<vobj, GroupName::Sp, MatrixFormat::REDUCED, fp_fmt>
{
  using in_type  = typename vobj::scalar_object;
  using tmp_type = typename std::tuple_element_t<static_cast<int>(fp_fmt), std::tuple<LorentzColourMatrixD,LorentzColourMatrixF>>;
  using out_type = typename std::tuple_element_t<static_cast<int>(fp_fmt), std::tuple<LorentzColourNx2ND,LorentzColourNx2NF>>;
 
  BinarySimpleUnmunger<tmp_type, in_type> binary_unmunger;
  GaugeSpunmunger<out_type, tmp_type> gauge_unmunger;
 
  void operator() (in_type &in, out_type &out){
    tmp_type tmp;
    binary_unmunger(in, tmp);
    gauge_unmunger(tmp, out);
  }
};
 
 
NAMESPACE_END(Grid);