Merge pull request #358 from mmphys/feature/serialisation-test

Add a ragged std::vector to the serialisation test

.gitignore

@@ -88,6 +88,7 @@ Thumbs.db
 # build directory #
 ###################
 build*/*
+Documentation/_build
 
 # IDE related files #
 #####################

Grid/qcd/smearing/WilsonFlow.h

@@ -7,7 +7,6 @@ Source file: ./lib/qcd/modules/plaquette.h
 Copyright (C) 2017
 
 Author: Guido Cossu <guido.cossu@ed.ac.uk>
-Author: Chulwoo Jung <chulwoo@bnl.gov>
 
 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
@@ -36,7 +35,7 @@ template <class Gimpl>
 class WilsonFlow: public Smear<Gimpl>{
   unsigned int Nstep;
   unsigned int measure_interval;
-  mutable RealD epsilon, taus,tolerance;
+  mutable RealD epsilon, taus;
 
 
   mutable WilsonGaugeAction<Gimpl> SG;
@@ -48,15 +47,13 @@ class WilsonFlow: public Smear<Gimpl>{
 public:
   INHERIT_GIMPL_TYPES(Gimpl)
 
-  explicit WilsonFlow(unsigned int Nstep, RealD epsilon, unsigned int interval = 1, RealD tol = 1e-3):
+  explicit WilsonFlow(unsigned int Nstep, RealD epsilon, unsigned int interval = 1):
   Nstep(Nstep),
     epsilon(epsilon),
-    tolerance(tol),
     measure_interval(interval),
     SG(WilsonGaugeAction<Gimpl>(3.0)) {
     // WilsonGaugeAction with beta 3.0
     assert(epsilon > 0.0);
-    assert(tolerance > 0.0);
     LogMessage();
   }
 
@@ -67,8 +64,6 @@ public:
 	      << "[WilsonFlow] epsilon : " << epsilon << std::endl;
     std::cout << GridLogMessage
 	      << "[WilsonFlow] full trajectory : " << Nstep * epsilon << std::endl;
-    std::cout << GridLogMessage
-	      << "[WilsonFlow] tolerance : " << tolerance << std::endl;
   }
 
   virtual void smear(GaugeField&, const GaugeField&) const;
@@ -111,14 +106,11 @@ void WilsonFlow<Gimpl>::evolve_step_adaptive(typename Gimpl::GaugeField &U, Real
   if (maxTau - taus < epsilon){
     epsilon = maxTau-taus;
   }
-  std::cout << GridLogMessage << "Integration epsilon : " << epsilon << std::endl;
+  //std::cout << GridLogMessage << "Integration epsilon : " << epsilon << std::endl;
   GaugeField Z(U.Grid());
   GaugeField Zprime(U.Grid());
-  GaugeField tmp(U.Grid()), Uprime(U.Grid()),Usave(U.Grid());
-
+  GaugeField tmp(U.Grid()), Uprime(U.Grid());
   Uprime = U;
-  Usave = U;
-
   SG.deriv(U, Z);
   Zprime = -Z;
   Z *= 0.25;                                  // Z0 = 1/4 * F(U)
@@ -136,33 +128,18 @@ void WilsonFlow<Gimpl>::evolve_step_adaptive(typename Gimpl::GaugeField &U, Real
   Z *= 3.0/4.0;                               // Z = 17/36*Z0 -8/9*Z1 +3/4*Z2
   Gimpl::update_field(Z, U, -2.0*epsilon);    // V(t+e) = exp(ep*Z)*W2
 
-  // Ramos arXiv:1301.4388
+  // Ramos 
   Gimpl::update_field(Zprime, Uprime, -2.0*epsilon); // V'(t+e) = exp(ep*Z')*W0
   // Compute distance as norm^2 of the difference
   GaugeField diffU = U - Uprime;
-// Wrong
-//  RealD diff = norm2(diffU);
-//  std::cout << GridLogMessage << "norm2: " << diff << std::endl;
+  RealD diff = norm2(diffU);
+  // adjust integration step
     
-//  RealD tol=1e-3;
-  
-  RealD diff = real(rankInnerMax(diffU,diffU));
-  diff = sqrt(diff)/18.; // distance defined in Ramos 
-
-  GridBase *grid = diffU.Grid();
-  std::cout << GridLogMessage << "max: " << diff << std::endl;
-  grid->GlobalMax(diff);
-  std::cout << GridLogMessage << "max: " << diff << std::endl;
-    
-  if(diff < tolerance) {
   taus += epsilon;
-//  std::cout << GridLogMessage << "Adjusting integration step with distance: " << diff << std::endl;
-  } else {
-    U = Usave;
-  }
+  //std::cout << GridLogMessage << "Adjusting integration step with distance: " << diff << std::endl;
     
-  epsilon = epsilon*0.95*std::pow(tolerance/diff,1./3.);
-  std::cout << GridLogMessage << "Distance : "<<diff<<"New epsilon : " << epsilon << std::endl;
+  epsilon = epsilon*0.95*std::pow(1e-4/diff,1./3.);
+  //std::cout << GridLogMessage << "New epsilon : " << epsilon << std::endl;
 
 }
 
@@ -207,11 +184,8 @@ void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const {
 template <class Gimpl>
 void WilsonFlow<Gimpl>::smear_adaptive(GaugeField& out, const GaugeField& in, RealD maxTau){
   out = in;
-//  taus = epsilon;
-  taus = 0.;
+  taus = epsilon;
   unsigned int step = 0;
-  double measTau = epsilon*measure_interval;
-  std::cout << GridLogMessage << "measTau :"<< measTau << std::endl;
   do{
     step++;
     //std::cout << GridLogMessage << "Evolution time :"<< taus << std::endl;
@@ -219,12 +193,10 @@ void WilsonFlow<Gimpl>::smear_adaptive(GaugeField& out, const GaugeField& in, Re
     std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : "
 		  << step << "  " << taus << "  "
 	      << energyDensityPlaquette(out) << std::endl;
-//    if( step % measure_interval == 0){
-    if( taus >  measTau ) {
+    if( step % measure_interval == 0){
       std::cout << GridLogMessage << "[WilsonFlow] Top. charge           : "
 		<< step << "  " 
 		<< WilsonLoops<PeriodicGimplR>::TopologicalCharge(out) << std::endl;
-      measTau += epsilon*measure_interval;
     }
   } while (taus < maxTau);
 

Grid/serialisation/BaseIO.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/serialisation/BaseIO.h
+
+Copyright (C) 2015
+
+Author: Michael Marshall <michael.marshall@ed.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 <Grid/GridCore.h>
+
+NAMESPACE_BEGIN(Grid)
+
+std::uint64_t EigenIO::EigenResizeCounter(0);
+
+NAMESPACE_END(Grid)

Grid/serialisation/BaseIO.h

@@ -9,6 +9,7 @@
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: Guido Cossu <guido.cossu@ed.ac.uk>
+Author: Michael Marshall <michael.marshall@ed.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
@@ -30,6 +31,7 @@ Author: Guido Cossu <guido.cossu@ed.ac.uk>
 #ifndef GRID_SERIALISATION_ABSTRACT_READER_H
 #define GRID_SERIALISATION_ABSTRACT_READER_H
 
+#include <atomic>
 #include <type_traits>
 #include <Grid/tensors/Tensors.h>
 #include <Grid/serialisation/VectorUtils.h>
@@ -110,6 +112,10 @@ namespace Grid {
     template <typename ET>
     inline typename std::enable_if<is_tensor_of_container<ET>::value, typename Traits<ET>::scalar_type *>::type
     getFirstScalar(ET &eigenTensor) { return eigenTensor.data()->begin(); }
+
+    // Counter for resized EigenTensors (poor man's substitute for allocator)
+    // Defined in BinaryIO.cc
+    extern std::uint64_t EigenResizeCounter;
   }
 
   // Abstract writer/reader classes ////////////////////////////////////////////
@@ -497,8 +503,14 @@ namespace Grid {
   typename std::enable_if<EigenIO::is_tensor_variable<ETensor>::value, void>::type
   Reader<T>::Reshape(ETensor &t, const std::array<typename ETensor::Index, ETensor::NumDimensions> &dims )
   {
+#ifdef GRID_OMP
+    // The memory counter is the reason this must be done from the primary thread
+    assert(omp_in_parallel()==0 && "Deserialisation which resizes Eigen tensor must happen from primary thread");
+#endif
+    EigenIO::EigenResizeCounter -= static_cast<uint64_t>(t.size()) * sizeof(typename ETensor::Scalar);
     //t.reshape( dims );
     t.resize( dims );
+    EigenIO::EigenResizeCounter += static_cast<uint64_t>(t.size()) * sizeof(typename ETensor::Scalar);
   }
 
   template <typename T>

Grid/serialisation/Hdf5IO.cc

@@ -1,3 +1,34 @@
+/*************************************************************************************
+ 
+ Grid physics library, www.github.com/paboyle/Grid
+ 
+ Source file: ./Grid/serialisation/VectorUtils.h
+ 
+ Copyright (C) 2015
+ 
+ Author: Antonin Portelli <antonin.portelli@me.com>
+ Author: Peter Boyle <paboyle@ed.ac.uk>
+ Author: Guido Cossu <guido.cossu@ed.ac.uk>
+ Author: Michael Marshall <michael.marshall@ed.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 <Grid/Grid.h>
 
 using namespace Grid;

Grid/serialisation/Hdf5IO.h

@@ -1,3 +1,34 @@
+/*************************************************************************************
+ 
+ Grid physics library, www.github.com/paboyle/Grid
+ 
+ Source file: ./Grid/serialisation/VectorUtils.h
+ 
+ Copyright (C) 2015
+ 
+ Author: Peter Boyle <paboyle@ed.ac.uk>
+ Author: Antonin Portelli <antonin.portelli@me.com>
+ Author: Guido Cossu <guido.cossu@ed.ac.uk>
+ Author: Michael Marshall <michael.marshall@ed.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 */
+
 #ifndef GRID_SERIALISATION_HDF5_H
 #define GRID_SERIALISATION_HDF5_H
 
@@ -34,11 +65,13 @@ namespace Grid
     template <typename U>
     void writeDefault(const std::string &s, const U &x);
     template <typename U>
-    typename std::enable_if<element<std::vector<U>>::is_number, void>::type
+    void writeRagged(const std::string &s, const std::vector<U> &x);
+    template <typename U>
+    typename std::enable_if<is_flattenable<std::vector<U>>::value>::type
     writeDefault(const std::string &s, const std::vector<U> &x);
     template <typename U>
-    typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
-    writeDefault(const std::string &s, const std::vector<U> &x);
+    typename std::enable_if<!is_flattenable<std::vector<U>>::value>::type
+    writeDefault(const std::string &s, const std::vector<U> &x) { writeRagged(s, x); }
     template <typename U>
     void writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const U * pDataRowMajor, size_t NumElements);
     H5NS::Group & getGroup(void);
@@ -64,11 +97,13 @@ namespace Grid
     template <typename U>
     void readDefault(const std::string &s, U &output);
     template <typename U>
-    typename std::enable_if<element<std::vector<U>>::is_number, void>::type
+    void readRagged(const std::string &s, std::vector<U> &x);
+    template <typename U>
+    typename std::enable_if<is_flattenable<std::vector<U>>::value>::type
     readDefault(const std::string &s, std::vector<U> &x);
     template <typename U>
-    typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
-    readDefault(const std::string &s, std::vector<U> &x);
+    typename std::enable_if<!is_flattenable<std::vector<U>>::value>::type
+    readDefault(const std::string &s, std::vector<U> &x) { readRagged(s, x); }
     template <typename U>
     void readMultiDim(const std::string &s, std::vector<U> &buf, std::vector<size_t> &dim);
     H5NS::Group & getGroup(void);
@@ -176,11 +211,13 @@ namespace Grid
   }
 
   template <typename U>
-  typename std::enable_if<element<std::vector<U>>::is_number, void>::type
+  typename std::enable_if<is_flattenable<std::vector<U>>::value>::type
   Hdf5Writer::writeDefault(const std::string &s, const std::vector<U> &x)
+  {
+    if (isRegularShape(x))
     {
       // alias to element type
-    typedef typename element<std::vector<U>>::type Element;
+      using Scalar = typename is_flattenable<std::vector<U>>::type;
       
       // flatten the vector and getting dimensions
       Flatten<std::vector<U>> flat(x);
@@ -188,12 +225,16 @@ namespace Grid
       const auto           &flatx = flat.getFlatVector();
       for (auto &d: flat.getDim())
         dim.push_back(d);
-    writeMultiDim<Element>(s, dim, &flatx[0], flatx.size());
+      writeMultiDim<Scalar>(s, dim, &flatx[0], flatx.size());
+    }
+    else
+    {
+      writeRagged(s, x);
+    }
   }
   
   template <typename U>
-  typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
-  Hdf5Writer::writeDefault(const std::string &s, const std::vector<U> &x)
+  void Hdf5Writer::writeRagged(const std::string &s, const std::vector<U> &x)
   {
     push(s);
     writeSingleAttribute(x.size(), HDF5_GRID_GUARD "vector_size",
@@ -229,7 +270,7 @@ namespace Grid
   void Hdf5Reader::readMultiDim(const std::string &s, std::vector<U> &buf, std::vector<size_t> &dim)
   {
     // alias to element type
-    typedef typename element<std::vector<U>>::type Element;
+    using Scalar = typename is_flattenable<std::vector<U>>::type;
     
     // read the dimensions
     H5NS::DataSpace       dataSpace;
@@ -260,26 +301,33 @@ namespace Grid
       H5NS::DataSet dataSet;
       
       dataSet = group_.openDataSet(s);
-      dataSet.read(buf.data(), Hdf5Type<Element>::type());
+      dataSet.read(buf.data(), Hdf5Type<Scalar>::type());
     }
     else
     {
       H5NS::Attribute attribute;
       
       attribute = group_.openAttribute(s);
-      attribute.read(Hdf5Type<Element>::type(), buf.data());
+      attribute.read(Hdf5Type<Scalar>::type(), buf.data());
     }
   }
 
   template <typename U>
-  typename std::enable_if<element<std::vector<U>>::is_number, void>::type
+  typename std::enable_if<is_flattenable<std::vector<U>>::value>::type
   Hdf5Reader::readDefault(const std::string &s, std::vector<U> &x)
+  {
+    if (H5Lexists        (group_.getId(), s.c_str(), H5P_DEFAULT) > 0
+     && H5Aexists_by_name(group_.getId(), s.c_str(), HDF5_GRID_GUARD "vector_size", H5P_DEFAULT ) > 0)
+    {
+      readRagged(s, x);
+    }
+    else
     {
       // alias to element type
-    typedef typename element<std::vector<U>>::type Element;
+      using Scalar = typename is_flattenable<std::vector<U>>::type;
 
       std::vector<size_t>   dim;
-    std::vector<Element>  buf;
+      std::vector<Scalar>   buf;
       readMultiDim( s, buf, dim );
 
       // reconstruct the multidimensional vector
@@ -287,10 +335,10 @@ namespace Grid
 
       x = r.getVector();
     }
+  }
   
   template <typename U>
-  typename std::enable_if<!element<std::vector<U>>::is_number, void>::type
-  Hdf5Reader::readDefault(const std::string &s, std::vector<U> &x)
+  void Hdf5Reader::readRagged(const std::string &s, std::vector<U> &x)
   {
     uint64_t size;
     

Grid/serialisation/MacroMagic.h

@@ -118,13 +118,13 @@ static inline std::string SerialisableClassName(void) {return std::string(#cname
 static constexpr bool isEnum = false; \
 GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_MEMBER,__VA_ARGS__))\
 template <typename T>\
-static inline void write(Writer<T> &WR,const std::string &s, const cname &obj){ \
+static inline void write(::Grid::Writer<T> &WR,const std::string &s, const cname &obj){ \
   push(WR,s);\
   GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_WRITE_MEMBER,__VA_ARGS__))	\
   pop(WR);\
 }\
 template <typename T>\
-static inline void read(Reader<T> &RD,const std::string &s, cname &obj){	\
+static inline void read(::Grid::Reader<T> &RD,const std::string &s, cname &obj){	\
   if (!push(RD,s))\
   {\
     std::cout << ::Grid::GridLogWarning << "IO: Cannot open node '" << s << "'" << std::endl; \

Grid/serialisation/VectorUtils.h

@@ -9,6 +9,7 @@
  Author: Antonin Portelli <antonin.portelli@me.com>
  Author: Peter Boyle <paboyle@ph.ed.ac.uk>
  Author: paboyle <paboyle@ph.ed.ac.uk>
+ Author: Michael Marshall <michael.marshall@ed.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
@@ -236,21 +237,36 @@ namespace Grid {
     }
   }
 
-  // Vector element trait //////////////////////////////////////////////////////  
-  template <typename T>
-  struct element
+  // is_flattenable<T>::value is true if T is a std::vector<> which can be flattened //////////////////////
+  template <typename T, typename V = void>
+  struct is_flattenable : std::false_type
   {
-    typedef T type;
-    static constexpr bool is_number = false;
+    using type      = T;
+    using grid_type = T;
+    static constexpr int vecRank = 0;
+    static constexpr bool isGridTensor = false;
+    static constexpr bool children_flattenable = std::is_arithmetic<T>::value or is_complex<T>::value;
   };
 
   template <typename T>
-  struct element<std::vector<T>>
+  struct is_flattenable<T, typename std::enable_if<isGridTensor<T>::value>::type> : std::false_type
   {
-    typedef typename element<T>::type type;
-    static constexpr bool is_number = std::is_arithmetic<T>::value
-                                      or is_complex<T>::value
-                                      or element<T>::is_number;
+    using type      = typename GridTypeMapper<T>::scalar_type;
+    using grid_type = T;
+    static constexpr int vecRank = 0;
+    static constexpr bool isGridTensor = true;
+    static constexpr bool children_flattenable = true;
+  };
+
+  template <typename T>
+  struct is_flattenable<std::vector<T>, typename std::enable_if<is_flattenable<T>::children_flattenable>::type>
+  : std::true_type
+  {
+    using type      = typename is_flattenable<T>::type;
+    using grid_type = typename is_flattenable<T>::grid_type;
+    static constexpr bool isGridTensor = is_flattenable<T>::isGridTensor;
+    static constexpr int vecRank = is_flattenable<T>::vecRank + 1;
+    static constexpr bool children_flattenable = true;
   };
   
   // Vector flattening utility class ////////////////////////////////////////////
@@ -259,22 +275,29 @@ namespace Grid {
   class Flatten
   {
   public:
-    typedef typename element<V>::type Element;
+    using Scalar  = typename is_flattenable<V>::type;
+    static constexpr bool isGridTensor = is_flattenable<V>::isGridTensor;
   public:
     explicit                    Flatten(const V &vector);
-    const V &                    getVector(void);
-    const std::vector<Element> & getFlatVector(void);
-    const std::vector<size_t>  & getDim(void);
+    const V &                   getVector(void)     const { return vector_; }
+    const std::vector<Scalar> & getFlatVector(void) const { return flatVector_; }
+    const std::vector<size_t> & getDim(void)        const { return dim_; }
   private:
-    void accumulate(const Element &e);
-    template <typename W>
-    void accumulate(const W &v);
-    void accumulateDim(const Element &e);
-    template <typename W>
-    void accumulateDim(const W &v);
+    template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
+    accumulate(const W &e);
+    template <typename W> typename std::enable_if<!is_flattenable<W>::value &&  is_flattenable<W>::isGridTensor>::type
+    accumulate(const W &e);
+    template <typename W> typename std::enable_if< is_flattenable<W>::value>::type
+    accumulate(const W &v);
+    template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
+    accumulateDim(const W &e) {} // Innermost is a scalar - do nothing
+    template <typename W> typename std::enable_if<!is_flattenable<W>::value &&  is_flattenable<W>::isGridTensor>::type
+    accumulateDim(const W &e);
+    template <typename W> typename std::enable_if< is_flattenable<W>::value>::type
+    accumulateDim(const W &v);
   private:
     const V             &vector_;
-    std::vector<Element> flatVector_;
+    std::vector<Scalar> flatVector_;
     std::vector<size_t> dim_;
   };
   
@@ -283,23 +306,32 @@ namespace Grid {
   class Reconstruct
   {
   public:
-    typedef typename element<V>::type Element;
+    using Scalar  = typename is_flattenable<V>::type;
+    static constexpr bool isGridTensor = is_flattenable<V>::isGridTensor;
   public:
-    Reconstruct(const std::vector<Element> &flatVector,
+    Reconstruct(const std::vector<Scalar> &flatVector,
                 const std::vector<size_t> &dim);
-    const V &                    getVector(void);
-    const std::vector<Element> & getFlatVector(void);
-    const std::vector<size_t>  & getDim(void);
+    const V &                   getVector(void)     const { return vector_; }
+    const std::vector<Scalar> & getFlatVector(void) const { return flatVector_; }
+    const std::vector<size_t> & getDim(void)        const { return dim_; }
   private:
-    void fill(std::vector<Element> &v);
-    template <typename W>
-    void fill(W &v);
-    void resize(std::vector<Element> &v, const unsigned int dim);
-    template <typename W>
-    void resize(W &v, const unsigned int dim);
+    template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
+    fill(W &v);
+    template <typename W> typename std::enable_if<!is_flattenable<W>::value &&  is_flattenable<W>::isGridTensor>::type
+    fill(W &v);
+    template <typename W> typename std::enable_if< is_flattenable<W>::value>::type
+    fill(W &v);
+    template <typename W> typename std::enable_if< is_flattenable<W>::value &&  is_flattenable<W>::vecRank==1>::type
+    resize(W &v, const unsigned int dim);
+    template <typename W> typename std::enable_if< is_flattenable<W>::value && (is_flattenable<W>::vecRank>1)>::type
+    resize(W &v, const unsigned int dim);
+    template <typename W> typename std::enable_if<!is_flattenable<W>::isGridTensor>::type
+    checkInnermost(const W &e) {} // Innermost is a scalar - do nothing
+    template <typename W> typename std::enable_if< is_flattenable<W>::isGridTensor>::type
+    checkInnermost(const W &e);
   private:
     V                         vector_;
-    const std::vector<Element> &flatVector_;
+    const std::vector<Scalar> &flatVector_;
     std::vector<size_t>       dim_;
     size_t                    ind_{0};
     unsigned int              dimInd_{0};
@@ -307,14 +339,24 @@ namespace Grid {
 
   // Flatten class template implementation
   template <typename V>
-  void Flatten<V>::accumulate(const Element &e)
+  template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
+  Flatten<V>::accumulate(const W &e)
   {
     flatVector_.push_back(e);
   }
   
   template <typename V>
-  template <typename W>
-  void Flatten<V>::accumulate(const W &v)
+  template <typename W> typename std::enable_if<!is_flattenable<W>::value && is_flattenable<W>::isGridTensor>::type
+  Flatten<V>::accumulate(const W &e)
+  {
+    for (const Scalar &x: e) {
+      flatVector_.push_back(x);
+    }
+  }
+
+  template <typename V>
+  template <typename W> typename std::enable_if<is_flattenable<W>::value>::type
+  Flatten<V>::accumulate(const W &v)
   {
     for (auto &e: v)
     {
@@ -323,11 +365,17 @@ namespace Grid {
   }
   
   template <typename V>
-  void Flatten<V>::accumulateDim(const Element &e) {};
+  template <typename W> typename std::enable_if<!is_flattenable<W>::value && is_flattenable<W>::isGridTensor>::type
+  Flatten<V>::accumulateDim(const W &e)
+  {
+    using Traits = GridTypeMapper<typename is_flattenable<W>::grid_type>;
+    for (int rank=0; rank < Traits::Rank; ++rank)
+      dim_.push_back(Traits::Dimension(rank));
+  }
   
   template <typename V>
-  template <typename W>
-  void Flatten<V>::accumulateDim(const W &v)
+  template <typename W> typename std::enable_if<is_flattenable<W>::value>::type
+  Flatten<V>::accumulateDim(const W &v)
   {
     dim_.push_back(v.size());
     accumulateDim(v[0]);
@@ -337,32 +385,26 @@ namespace Grid {
   Flatten<V>::Flatten(const V &vector)
   : vector_(vector)
   {
-    accumulate(vector_);
     accumulateDim(vector_);
+    std::size_t TotalSize{ dim_[0] };
+    for (int i = 1; i < dim_.size(); ++i) {
+      TotalSize *= dim_[i];
     }
-  
-  template <typename V>
-  const V & Flatten<V>::getVector(void)
-  {
-    return vector_;
-  }
-  
-  template <typename V>
-  const std::vector<typename Flatten<V>::Element> &
-  Flatten<V>::getFlatVector(void)
-  {
-    return flatVector_;
-  }
-  
-  template <typename V>
-  const std::vector<size_t> & Flatten<V>::getDim(void)
-  {
-    return dim_;
+    flatVector_.reserve(TotalSize);
+    accumulate(vector_);
   }
   
   // Reconstruct class template implementation
   template <typename V>
-  void Reconstruct<V>::fill(std::vector<Element> &v)
+  template <typename W> typename std::enable_if<!is_flattenable<W>::value && !is_flattenable<W>::isGridTensor>::type
+  Reconstruct<V>::fill(W &v)
+  {
+    v = flatVector_[ind_++];
+  }
+  
+  template <typename V>
+  template <typename W> typename std::enable_if<!is_flattenable<W>::value &&  is_flattenable<W>::isGridTensor>::type
+  Reconstruct<V>::fill(W &v)
   {
     for (auto &e: v)
     {
@@ -371,8 +413,8 @@ namespace Grid {
   }
 
   template <typename V>
-  template <typename W>
-  void Reconstruct<V>::fill(W &v)
+  template <typename W> typename std::enable_if<is_flattenable<W>::value>::type
+  Reconstruct<V>::fill(W &v)
   {
     for (auto &e: v)
     {
@@ -381,14 +423,15 @@ namespace Grid {
   }
   
   template <typename V>
-  void Reconstruct<V>::resize(std::vector<Element> &v, const unsigned int dim)
+  template <typename W> typename std::enable_if<is_flattenable<W>::value && is_flattenable<W>::vecRank==1>::type
+  Reconstruct<V>::resize(W &v, const unsigned int dim)
   {
     v.resize(dim_[dim]);
   }
   
   template <typename V>
-  template <typename W>
-  void Reconstruct<V>::resize(W &v, const unsigned int dim)
+  template <typename W> typename std::enable_if<is_flattenable<W>::value && (is_flattenable<W>::vecRank>1)>::type
+  Reconstruct<V>::resize(W &v, const unsigned int dim)
   {
     v.resize(dim_[dim]);
     for (auto &e: v)
@@ -398,34 +441,31 @@ namespace Grid {
   }
   
   template <typename V>
-  Reconstruct<V>::Reconstruct(const std::vector<Element> &flatVector,
+  template <typename W> typename std::enable_if<is_flattenable<W>::isGridTensor>::type
+  Reconstruct<V>::checkInnermost(const W &)
+  {
+    using Traits = GridTypeMapper<typename is_flattenable<W>::grid_type>;
+    const int gridRank{Traits::Rank};
+    const int dimRank{static_cast<int>(dim_.size())};
+    assert(dimRank >= gridRank && "Tensor rank too low for Grid tensor");
+    for (int i=0; i<gridRank; ++i) {
+      assert(dim_[dimRank - gridRank + i] == Traits::Dimension(i) && "Tensor dimension doesn't match Grid tensor");
+    }
+    dim_.resize(dimRank - gridRank);
+  }
+
+  template <typename V>
+  Reconstruct<V>::Reconstruct(const std::vector<Scalar> &flatVector,
                               const std::vector<size_t> &dim)
   : flatVector_(flatVector)
   , dim_(dim)
   {
+    checkInnermost(vector_);
+    assert(dim_.size() == is_flattenable<V>::vecRank && "Tensor rank doesn't match nested std::vector rank");
     resize(vector_, 0);
     fill(vector_);
   }
   
-  template <typename V>
-  const V & Reconstruct<V>::getVector(void)
-  {
-    return vector_;
-  }
-  
-  template <typename V>
-  const std::vector<typename Reconstruct<V>::Element> &
-  Reconstruct<V>::getFlatVector(void)
-  {
-    return flatVector_;
-  }
-  
-  template <typename V>
-  const std::vector<size_t> & Reconstruct<V>::getDim(void)
-  {
-    return dim_;
-  }
-
   // Vector IO utilities ///////////////////////////////////////////////////////
   // helper function to read space-separated values
   template <typename T>
@@ -459,6 +499,64 @@ namespace Grid {
     
     return os;
   }
+
+  // In general, scalar types are considered "flattenable" (regularly shaped)
+  template <typename T>
+  bool isRegularShapeHelper(const std::vector<T> &, std::vector<std::size_t> &, int, bool)
+  {
+    return true;
+  }
+
+  template <typename T>
+  bool isRegularShapeHelper(const std::vector<std::vector<T>> &v, std::vector<std::size_t> &Dims, int Depth, bool bFirst)
+  {
+    if( bFirst)
+    {
+      assert( Dims.size() == Depth     && "Bug: Delete this message after testing" );
+      Dims.push_back(v[0].size());
+      if (!Dims[Depth])
+        return false;
+    }
+    else
+    {
+      assert( Dims.size() >= Depth + 1 && "Bug: Delete this message after testing" );
+    }
+    for (std::size_t i = 0; i < v.size(); ++i)
+    {
+      if (v[i].size() != Dims[Depth] || !isRegularShapeHelper(v[i], Dims, Depth + 1, bFirst && i==0))
+      {
+        return false;
+      }
+    }
+    return true;
+  }
+
+  template <typename T>
+  bool isRegularShape(const T &t) { return true; }
+
+  template <typename T>
+  bool isRegularShape(const std::vector<T> &v) { return !v.empty(); }
+
+  // Return non-zero if all dimensions of this std::vector<std::vector<T>> are regularly shaped
+  template <typename T>
+  bool isRegularShape(const std::vector<std::vector<T>> &v)
+  {
+    if (v.empty() || v[0].empty())
+      return false;
+    // Make sure all of my rows are the same size
+    std::vector<std::size_t> Dims;
+    Dims.reserve(is_flattenable<T>::vecRank);
+    Dims.push_back(v.size());
+    Dims.push_back(v[0].size());
+    for (std::size_t i = 0; i < Dims[0]; ++i)
+    {
+      if (v[i].size() != Dims[1] || !isRegularShapeHelper(v[i], Dims, 2, i==0))
+      {
+        return false;
+      }
+    }
+    return true;
+  }
 }
 
 // helper function to read space-separated values

Grid/tensors/Tensor_class.h

@@ -417,7 +417,7 @@ public:
       stream << "{";
       for (int j = 0; j < N; j++) {
 	stream << o._internal[i][j];
-	if (i < N - 1) stream << ",";
+	if (j < N - 1) stream << ",";
       }
       stream << "}";
       if (i != N - 1) stream << "\n\t\t";

documentation/manual.rst

@@ -1787,7 +1787,7 @@ Hdf5Writer     Hdf5Reader       HDF5
 
 Write interfaces, similar to the XML facilities in QDP++ are presented. However,
 the serialisation routines are automatically generated by the macro, and a virtual
-reader adn writer interface enables writing to any of a number of formats.
+reader and writer interface enables writing to any of a number of formats.
 
 **Example**::
 
@@ -1814,6 +1814,91 @@ reader adn writer interface enables writing to any of a number of formats.
    }
 
 
+Eigen tensor support -- added 2019H1
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The Serialisation library was expanded in 2019 to support de/serialisation of
+Eigen tensors. De/serialisation of existing types was not changed. Data files
+without Eigen tensors remain compatible with earlier versions of Grid and other readers.
+Conversely, data files containing serialised Eigen tensors is a breaking change.
+
+Eigen tensor serialisation support was added to BaseIO, which was modified to provide a Traits class
+to recognise Eigen tensors with elements that are either: primitive scalars (arithmetic and complex types);
+or Grid tensors.
+
+**Traits determining de/serialisable scalars**::
+
+  // Is this an Eigen tensor
+  template<typename T> struct is_tensor : std::integral_constant<bool,
+    std::is_base_of<Eigen::TensorBase<T, Eigen::ReadOnlyAccessors>, T>::value> {};
+  // Is this an Eigen tensor of a supported scalar
+  template<typename T, typename V = void> struct is_tensor_of_scalar : public std::false_type {};
+  template<typename T> struct is_tensor_of_scalar<T, typename std::enable_if<is_tensor<T>::value && is_scalar<typename T::Scalar>::value>::type> : public std::true_type {};
+  // Is this an Eigen tensor of a supported container
+  template<typename T, typename V = void> struct is_tensor_of_container : public std::false_type {};
+  template<typename T> struct is_tensor_of_container<T, typename std::enable_if<is_tensor<T>::value && isGridTensor<typename T::Scalar>::value>::type> : public std::true_type {};
+
+
+Eigen tensors are regular, multidimensional objects, and each Reader/Writer
+was extended to support this new datatype. Where the Eigen tensor contains
+a Grid tensor, the dimensions of the data written are the dimensions of the
+Eigen tensor plus the dimensions of the underlying Grid scalar. Dimensions
+of size 1 are preserved.
+
+**New Reader/Writer methods for multi-dimensional data**::
+
+  template <typename U>
+  void readMultiDim(const std::string &s, std::vector<U> &buf, std::vector<size_t> &dim);
+  template <typename U>
+  void writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const U * pDataRowMajor, size_t NumElements);
+
+
+On readback, the Eigen tensor rank must match the data being read, but the tensor
+dimensions will be resized if necessary. Resizing is not possible for Eigen::TensorMap<T>
+because these tensors use a buffer provided at construction, and this buffer cannot be changed.
+Deserialisation failures cause Grid to assert.
+
+
+HDF5 Optimisations -- added June 2021
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Grid serialisation is intended to be light, deterministic and provide a layer of abstraction over
+multiple file formats. HDF5 excels at handling multi-dimensional data, and the Grid HDF5Reader/HDF5Writer exploits this.
+When serialising nested ``std::vector<T>``, where ``T`` is an arithmetic or complex type,
+the Hdf5Writer writes the data as an Hdf5 DataSet object.
+
+However, nested ``std::vector<std::vector<...T>>`` might be "ragged", i.e. not necessarily regular. E.g. a 3d nested
+``std::vector`` might contain 2 rows, the first being a 2x2 block and the second row being a 1 x 2 block.
+A bug existed whereby this was not checked on write, so nested, ragged vectors
+were written as a regular dataset, with a buffer under/overrun and jumbled contents.
+
+Clearly this was not used in production, as the bug went undetected until now. Fixing this bug
+is an opportunity to further optimise the HDF5 file format.
+
+The goals of this change are to:
+
+* Make changes to the Hdf5 file format only -- i.e. do not impact other file formats
+
+* Implement file format changes in such a way that they are transparent to the Grid reader
+
+* Correct the bug for ragged vectors of numeric / complex types
+
+* Extend the support of nested std::vector<T> to arbitrarily nested Grid tensors
+
+
+The trait class ``element`` has been redefined to ``is_flattenable``, which is a trait class for
+potentially "flattenable" objects. These are (possibly nested) ``std::vector<T>`` where ``T`` is
+an arithmetic, complex or Grid tensor type. Flattenable objects are tested on write
+(with the function ``isRegularShape``) to see whether they actually are regular.
+
+Flattenable, regular objects are written to a multidimensional HDF5 DataSet.
+Otherwise, an Hdf5 sub group is created with the object "name", and each element of the outer dimension is
+recursively written to as object "name_n", where n is a 0-indexed number.
+
+On readback (by Grid)), the presence of a subgroup containing the attribute ``Grid_vector_size`` triggers a
+"ragged read", otherwise a read from a DataSet is attempted.
+
+
 Data parallel field IO
 -----------------------
 

tests/IO/Test_serialisation.cc

@@ -48,7 +48,9 @@ public:
                           std::vector<double>, array,
                           std::vector<std::vector<double> >, twodimarray,
 			  std::vector<std::vector<std::vector<std::complex<double>> > >, cmplx3darray,
-			  SpinColourMatrix, scm
+			  SpinColourMatrix, scm,
+                          std::vector<std::vector<std::vector<int> > >, ragged,
+                          std::vector<std::vector<SpinColourMatrix> >, vscm
                           );
   myclass() {}
   myclass(int i)
@@ -56,6 +58,10 @@ public:
   , twodimarray(3,std::vector<double>(5, 1.23456))
   , cmplx3darray(3,std::vector<std::vector<std::complex<double>>>(5, std::vector<std::complex<double>>(7, std::complex<double>(1.2, 3.4))))
   , ve(2, myenum::blue)
+  , ragged( {{{i+1},{i+2,i+3}}, // ragged
+            {{i+4,i+5,i+6,i+7},{i+8,i+9,i+10,i+11},{i+12,i+13,i+14,i+15}}, // block
+            {{i+16,i+17},{i+18,i+19,i+20}}} ) //ragged
+  , vscm(3, std::vector<SpinColourMatrix>(5))
   {
     e=myenum::red;
     x=i;
@@ -68,6 +74,13 @@ public:
     scm()(0, 2)(1, 1) = 6.336;
     scm()(2, 1)(2, 2) = 7.344;
     scm()(1, 1)(2, 0) = 8.3534;
+    int Counter = i;
+    for( auto & v : vscm ) {
+      for( auto & j : v ) {
+        j = std::complex<double>(Counter, -Counter);
+        Counter++;
+      }
+    }
   }
 };
 

tests/smearing/Test_WilsonFlow.cc

@@ -33,7 +33,6 @@ namespace Grid{
     GRID_SERIALIZABLE_CLASS_MEMBERS(WFParameters,
             int, steps,
             double, step_size,
-            double, tol,
             int, meas_interval,
             double, maxTau); // for the adaptive algorithm
        
@@ -83,27 +82,13 @@ int main(int argc, char **argv) {
   SU<Nc>::HotConfiguration(pRNG, Umu);
   
   typedef Grid::XmlReader       Serialiser;
-//  Serialiser Reader("input.xml");
-//  WFParameters WFPar(Reader);
-//  ConfParameters CPar(Reader);
-//  WFParameters WFPar;
-  int steps = 800;
-  double step_size=0.02;
-  double tol=1e-4;
-  int meas_interval=50;
-  double maxTau = 16;
-//  ConfParameters CPar;
-//  CPar. conf_prefix="configurations/ckpoint_lat";
-//  CPar. rng_prefix="rngs/ckpoint_rng";
-//  CPar. StartConfiguration=100,
-//  CPar. EndConfiguration=110,
-//  CPar. Skip=1;
-//  CheckpointerParameters CPPar(CPar.conf_prefix, CPar.rng_prefix);
-  CheckpointerParameters CPPar("configurations/ckpoint_lat","rngs/ckpoint_rng");
+  Serialiser Reader("input.xml");
+  WFParameters WFPar(Reader);
+  ConfParameters CPar(Reader);
+  CheckpointerParameters CPPar(CPar.conf_prefix, CPar.rng_prefix);
   BinaryHmcCheckpointer<PeriodicGimplR> CPBin(CPPar);
 
-//  for (int conf = CPar.StartConfiguration; conf <= CPar.EndConfiguration; conf+= CPar.Skip){
-  for (int conf = 100; conf <= 110; conf+= 1){
+  for (int conf = CPar.StartConfiguration; conf <= CPar.EndConfiguration; conf+= CPar.Skip){
 
   CPBin.CheckpointRestore(conf, Umu, sRNG, pRNG);
 
@@ -111,10 +96,9 @@ int main(int argc, char **argv) {
   std::cout << GridLogMessage << "Initial plaquette: "
     << WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu) << std::endl;
 
-  WilsonFlow<PeriodicGimplR> WF(steps, step_size, meas_interval);
+  WilsonFlow<PeriodicGimplR> WF(WFPar.steps, WFPar.step_size, WFPar.meas_interval);
 
-//  WF.smear_adaptive(Uflow, Umu, maxTau);
-  WF.smear(Uflow, Umu);
+  WF.smear_adaptive(Uflow, Umu, WFPar.maxTau);
 
   RealD WFlow_plaq = WilsonLoops<PeriodicGimplR>::avgPlaquette(Uflow);
   RealD WFlow_TC   = WilsonLoops<PeriodicGimplR>::TopologicalCharge(Uflow);