Implemented perambulator read/write ... but in binary format. Will switch to Hdf5 when I have Antonins feedback

Hadrons/Modules/MDistil/Distil.hpp

@@ -201,28 +201,139 @@ inline GridCartesian * MakeLowerDimGrid( GridCartesian * gridHD )
  ******************************************************************************/
 
 template<typename Scalar_, int NumIndices_>
-class Perambulator : public Eigen::Tensor<Scalar_, NumIndices_, Eigen::RowMajor | Eigen::DontAlign>
+class NamedTensor : public Eigen::Tensor<Scalar_, NumIndices_, Eigen::RowMajor | Eigen::DontAlign>
 {
-protected:
 public:
+  typedef Eigen::Tensor<Scalar_, NumIndices_, Eigen::RowMajor | Eigen::DontAlign> ET;
   std::array<std::string,NumIndices_> IndexNames;
 public:
   template<typename... IndexTypes>
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Perambulator(std::array<std::string,NumIndices_> &IndexNames_, Eigen::Index firstDimension, IndexTypes... otherDimensions)
+  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NamedTensor(std::array<std::string,NumIndices_> &IndexNames_, Eigen::Index firstDimension, IndexTypes... otherDimensions)
   : IndexNames{IndexNames_}, Eigen::Tensor<Scalar_, NumIndices_, Eigen::RowMajor | Eigen::DontAlign>(firstDimension, otherDimensions...)
   {
     // The number of dimensions used to construct a tensor must be equal to the rank of the tensor.
     EIGEN_STATIC_ASSERT(sizeof...(otherDimensions) + 1 == NumIndices_, YOU_MADE_A_PROGRAMMING_MISTAKE)
   }
 
-  inline void WriteTemporary(const std::string &FileName){}
-
   // Share data for timeslices we calculated with other nodes
   inline void SliceShare( GridCartesian * gridLowDim, GridCartesian * gridHighDim ) {
     Grid::SliceShare( gridLowDim, gridHighDim, this->data(), (int) (this->size() * sizeof(Scalar_)));
   }
+
+  // load and save - not virtual - probably all changes
+  inline void load(const std::string filename);
+  inline void save(const std::string filename) const;
+  inline void ReadTemporary(const std::string filename);
+  inline void WriteTemporary(const std::string filename) const;
 };
 
+/******************************************************************************
+ Save NamedTensor binary format
+ ******************************************************************************/
+
+template<typename Scalar_, int NumIndices_>
+void NamedTensor<Scalar_, NumIndices_>::WriteTemporary(const std::string filename) const {
+  std::cout << GridLogMessage << "Writing NamedTensor to \"" << filename << "\"" << std::endl;
+  std::ofstream w(filename, std::ios::binary);
+  // total number of elements
+  uint32_t ul = htonl( static_cast<uint32_t>( this->size() ) );
+  w.write(reinterpret_cast<const char *>(&ul), sizeof(ul));
+  // number of dimensions
+  uint16_t us = htons( static_cast<uint16_t>( this->NumIndices ) );
+  w.write(reinterpret_cast<const char *>(&us), sizeof(us));
+  // dimensions together with names
+  int d = 0;
+  for( auto dim : this->dimensions() ) {
+    // size of this dimension
+    us = htons( static_cast<uint16_t>( dim ) );
+    w.write(reinterpret_cast<const char *>(&us), sizeof(us));
+    // length of this dimension name
+    us = htons( static_cast<uint16_t>( IndexNames[d].size() ) );
+    w.write(reinterpret_cast<const char *>(&us), sizeof(us));
+    // dimension name
+    w.write(IndexNames[d].c_str(), IndexNames[d].size());
+    d++;
+  }
+  // Actual data
+  w.write(reinterpret_cast<const char *>(this->data()),(int) (this->size() * sizeof(Scalar_)));
+}
+
+/******************************************************************************
+ Load NamedTensor binary format
+ ******************************************************************************/
+
+template<typename Scalar_, int NumIndices_>
+void NamedTensor<Scalar_, NumIndices_>::ReadTemporary(const std::string filename) {
+  std::cout << GridLogMessage << "Reading NamedTensor from \"" << filename << "\"" << std::endl;
+  std::ifstream r(filename, std::ios::binary);
+  // total number of elements
+  uint32_t ul;
+  r.read(reinterpret_cast<char *>(&ul), sizeof(ul));
+  assert( this->size() == ntohl( ul ) && "Error: total number of elements" );
+  // number of dimensions
+  uint16_t us;
+  r.read(reinterpret_cast<char *>(&us), sizeof(us));
+  assert( this->NumIndices == ntohs( us ) && "Error: number of dimensions" );
+  // dimensions together with names
+  int d = 0;
+  for( auto dim : this->dimensions() ) {
+    // size of dimension
+    r.read(reinterpret_cast<char *>(&us), sizeof(us));
+    assert( dim == ntohs( us ) && "size of dimension" );
+    // length of dimension name
+    r.read(reinterpret_cast<char *>(&us), sizeof(us));
+    size_t l = ntohs( us );
+    assert( l == IndexNames[d].size() && "length of dimension name" );
+    // dimension name
+    std::string s( l, '?' );
+    r.read(&s[0], l);
+    assert( s == IndexNames[d] && "dimension name" );
+    d++;
+  }
+  // Actual data
+  r.read(reinterpret_cast<char *>(this->data()),(int) (this->size() * sizeof(Scalar_)));
+}
+
+/******************************************************************************
+ Save NamedTensor Hdf5 format
+ ******************************************************************************/
+
+template<typename Scalar_, int NumIndices_>
+void NamedTensor<Scalar_, NumIndices_>::save(const std::string filename) const {
+  std::cout << GridLogMessage << "Writing NamedTensor to \"" << filename << "\"" << std::endl;
+#ifndef HAVE_HDF5
+  std::cout << GridErrorMessage << "Error: I/O for NamedTensor requires HDF5" << std::endl;
+#else
+  Hdf5Writer w(filename);
+  //w << this->NumIndices << this->dimensions() << this->IndexNames;
+#endif
+}
+
+/******************************************************************************
+ Load NamedTensor Hdf5 format
+ ******************************************************************************/
+
+template<typename Scalar_, int NumIndices_>
+void NamedTensor<Scalar_, NumIndices_>::load(const std::string filename) {
+  std::cout << GridLogMessage << "Reading NamedTensor from \"" << filename << "\"" << std::endl;
+#ifndef HAVE_HDF5
+  std::cout << GridErrorMessage << "Error: I/O for NamedTensor requires HDF5" << std::endl;
+#else
+  Hdf5Reader r(filename);
+  typename ET::Dimensions d;
+  std::array<std::string,NumIndices_> n;
+  //r >> this->NumIndices >> d >> n;
+  //this->IndexNames = n;
+#endif
+}
+
+/******************************************************************************
+ Perambulator object
+ ******************************************************************************/
+
+template<typename Scalar_, int NumIndices_>
+using Perambulator = NamedTensor<Scalar_, NumIndices_>;
+
 /*************************************************************************************
  
  Rotate eigenvectors into our phase convention

tests/hadrons/Test_hadrons_distil.cc

@@ -252,16 +252,16 @@ bool bNumber( int &ri, const char * & pstr, bool bGobbleWhiteSpace = true )
 
 #ifdef DEBUG
 
-typedef Eigen::Tensor<Complex,3,Eigen::RowMajor | Eigen::DontAlign> MyTensor;
+typedef Grid::Hadrons::MDistil::NamedTensor<Complex,3> MyTensor;
 
-void DebugShowTensor(MyTensor &x)
+void DebugShowTensor(MyTensor &x, const char * n)
 {
   const MyTensor::Index s{x.size()};
-  std::cout << "x.size() = " << s << std::endl;
-  std::cout << "x.NumDimensions = " << x.NumDimensions << " (TensorBase)" << std::endl;
-  std::cout << "x.NumIndices = " << x.NumIndices << std::endl;
+  std::cout << n << ".size() = " << s << std::endl;
+  std::cout << n << ".NumDimensions = " << x.NumDimensions << " (TensorBase)" << std::endl;
+  std::cout << n << ".NumIndices = " << x.NumIndices << std::endl;
   const MyTensor::Dimensions & d{x.dimensions()};
-  std::cout << "d.size() = " << d.size() << std::endl;
+  std::cout << n << ".dimensions().size() = " << d.size() << std::endl;
   std::cout << "Dimensions are ";
   for(auto i : d ) std::cout << "[" << i << "]";
   std::cout << std::endl;
@@ -287,24 +287,30 @@ void DebugShowTensor(MyTensor &x)
   Complex * p = x.data();
   for( auto i = 0 ; i < s ; i++ ) {
     if( i ) std::cout << ", ";
-    std::cout << "x.data()[" << i << "]=" << * p++;
+    std::cout << n << ".data()[" << i << "]=" << * p++;
   }
   std::cout << std::endl;
 }
 
 bool DebugEigenTest()
 {
-  MyTensor x(2,3,4);
-  DebugShowTensor(x);
-  // Test initialisation of an array of strings
+  const char pszTestFileName[] = "test_tensor.bin";
   std::array<std::string,3> as={"Alpha", "Beta", "Gamma"};
+  MyTensor x(as, 2,3,4);
+  DebugShowTensor(x, "x");
+  x.WriteTemporary(pszTestFileName);
+  // Test initialisation of an array of strings
   for( auto a : as )
     std::cout << a << std::endl;
   Grid::Hadrons::MDistil::Perambulator<Complex,3> p{as,2,7,2};
-  DebugShowTensor(p);
+  DebugShowTensor(p, "p");
   std::cout << "p.IndexNames follow" << std::endl;
   for( auto a : p.IndexNames )
     std::cout << a << std::endl;
+  // Now see whether we can read a tensor back
+  MyTensor y(as, 2,3,4);
+  y.ReadTemporary(pszTestFileName);
+  DebugShowTensor(y, "y");
   return true;
 }
 #endif