Implemented specialisations of NamedTensor as derived classes, however this suffers a number of problems:

1) virtual functions not available in base class constructor where I'd like to use them - e.g. IndexNames
2) Must define new constructors in derived classes
... so the specialisations are fatter than I'd like. Would prefer to revert to specifying tensor name and index name defaults in template

Hadrons/Distil.hpp

@@ -37,18 +37,18 @@ BEGIN_HADRONS_NAMESPACE
 
 /******************************************************************************
  NamedTensor
- This is an Eigen::Tensor of type Scalar_ and rank NumIndices_ (row-major order)
+ Eigen::Tensor of type Scalar_ and rank NumIndices_ (row-major order), together with a name for each index.
- They can be persisted to disk in tag Name_, and IndexNames are validated on load.
+ Index names are mutable, but tensor dimensionality is not (size of each dimension is mutable).
- TODO: WHAT TO SAVE / VALIDATE ON LOAD (Override to warn instead of assert on load)
+ They can be persisted to / restored from disk, by default using tag Name.
- Ensemble string
+ During restore from disk, these validations are performed:
- Configuration number
+   1) Tensor dimensionality must match
- Noise unique string
+   2) IndexNames are validated against current values
- Distillation parameters
+   3) If the tensor has non-zero size, the tensor being loaded must have same extent in each dimension
  ******************************************************************************/
 
 extern const std::string NamedTensorFileExtension;
 
-template<typename Scalar_, int NumIndices_, const std::string &Name_, const std::array<std::string,NumIndices_> &IndexNames_>
+template<typename Scalar_, int NumIndices_>
 class NamedTensor : Serializable
 {
 public:
@@ -59,40 +59,40 @@ public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(NamedTensor,
                                     ET,                       tensor,
                                     std::vector<std::string>, IndexNames );
-    
+
-    // Get the default index names as std::vector
+    // Name of the object and Index names as set in the constructor
-    std::vector<std::string> DefaultIndexNames()
+    const std::string              &Name;
-    {
+    const std::vector<std::string> &DefaultIndexNames;
-        std::vector<std::string> names{NumIndices_};
+
-        for (std::size_t i = 0; i < NumIndices_; i++)
+    virtual ~NamedTensor(){};
-            names[i] = IndexNames_[i];
-        return names;
-    }
-    
     // Default constructor (assumes tensor will be loaded from file)
-    EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NamedTensor() : IndexNames{DefaultIndexNames()} {}
+    EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NamedTensor(const std::string &Name_, const std::vector<std::string> &IndexNames_)
+    : IndexNames{IndexNames_}, Name{Name_}, DefaultIndexNames{IndexNames_} {}
     
     // Construct a named tensor explicitly specifying size of each dimension
     template<typename... IndexTypes>
-    EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NamedTensor(Eigen::Index firstDimension, IndexTypes... otherDimensions)
+    EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NamedTensor(const std::string &Name_, const std::vector<std::string> &IndexNames_,
-    : tensor(firstDimension, otherDimensions...), IndexNames{DefaultIndexNames()}
+                                                      Eigen::Index firstDimension, IndexTypes... otherDimensions)
+    : tensor(firstDimension, otherDimensions...), IndexNames{IndexNames_}, Name{Name_}, DefaultIndexNames{IndexNames_}
     {
         assert(sizeof...(otherDimensions) + 1 == NumIndices_ && "NamedTensor: dimensions != tensor rank");
     }
-    
+
     // Do my index names match the default for my type?
-    bool ValidateIndexNames() const
+    bool ValidateIndexNames( const std::vector<std::string> &CheckNames ) const
     {
+        assert( CheckNames.size() == NumIndices_ && "Bug: CheckNames don't match NumIndices_" );
         bool bSame{ IndexNames.size() == NumIndices_ };
         for( std::size_t i = 0; bSame && i < NumIndices_; i++ )
         {
-            bSame = IndexNames[i].size() == IndexNames_[i].size()
+            bSame = IndexNames[i].size() == CheckNames[i].size()
-            && std::equal( IndexNames[i].begin(), IndexNames[i].end(), IndexNames_[i].begin(),
+            && std::equal( IndexNames[i].begin(), IndexNames[i].end(), CheckNames[i].begin(),
                           [](const char & c1, const char & c2){ return c1 == c2 || std::toupper(c1) == std::toupper(c2); });
         }
         return bSame;
     }
-    
+    bool ValidateIndexNames() const { return ValidateIndexNames(DefaultIndexNames); }
+
 #ifdef HAVE_HDF5
     using Default_Reader = Grid::Hdf5Reader;
     using Default_Writer = Grid::Hdf5Writer;
@@ -101,40 +101,40 @@ public:
     using Default_Writer = Grid::BinaryWriter;
 #endif
     
-    template<typename Writer> void write(Writer &w, const std::string &Tag = Name_) const
+    void write(const std::string &FileName, const std::string &Tag) const
-    { write(w, Tag, *this); }
-    
-    inline void write(const std::string &filename, const std::string &Tag = Name_) const
     {
-        std::string sFileName{filename};
+        std::string FileName_{FileName};
-        sFileName.append( NamedTensorFileExtension );
+        FileName_.append( NamedTensorFileExtension );
-        LOG(Message) << "Writing " << Name_ << " to file " << sFileName << " tag " << Tag << std::endl;
+        LOG(Message) << "Writing " << Name << " to file " << FileName_ << " tag " << Tag << std::endl;
-        Default_Writer w( sFileName );
+        Default_Writer w( FileName_ );
-        write( w, Tag );
+        write( w, Tag, *this );
     }
-    
+    void write(const std::string &FileName) const { return write(FileName, Name); }
-    // Read and validate index names
+
-    template<typename Reader> void read(Reader &r, bool bValidate = true, const std::string &Tag = Name_)
+    // Read tensor.
+    // Validate:
+    //  1) index names (if requested)
+    //  2) index dimensions (if they are non-zero when called)
+    template<typename Reader> void read(Reader &r, bool bValidate, const std::string &Tag)
     {
         // Grab index names and dimensions
         std::vector<std::string> OldIndexNames{std::move(IndexNames)};
-        typename ET::Dimensions OldDimensions{tensor.dimensions()};
+        const typename ET::Dimensions OldDimensions{tensor.dimensions()};
         read(r, Tag, *this);
         const typename ET::Dimensions & NewDimensions{tensor.dimensions()};
         for (int i = 0; i < NumIndices_; i++)
             assert(OldDimensions[i] == 0 || OldDimensions[i] == NewDimensions[i] && "NamedTensor::read dimension size");
         if (bValidate)
-            assert(ValidateIndexNames() && "NamedTensor::read dimension name");
+            assert(ValidateIndexNames(OldIndexNames) && "NamedTensor::read dimension name");
     }
-    
+    template<typename Reader> void read(Reader &r, bool bValidate = true) { read(r, bValidate, Name); }
-    inline void read (const std::string &filename, bool bValidate = true, const std::string &Tag = Name_)
+
+    inline void read (const std::string &FileName, bool bValidate, const std::string &Tag)
     {
-        std::string sFileName{filename};
+        Default_Reader r(FileName + NamedTensorFileExtension);
-        sFileName.append( NamedTensorFileExtension );
-        LOG(Message) << "Reading " << Name_ << " from file " << sFileName << " tag " << Tag << std::endl;
-        Default_Reader r(sFileName);
         read(r, bValidate, Tag);
     }
+    inline void read (const std::string &FileName, bool bValidate= true) { return read(FileName, bValidate, Name); }
 };
 
 /******************************************************************************
@@ -147,11 +147,34 @@ BEGIN_MODULE_NAMESPACE(MDistil)
 using LapEvecs = Grid::Hadrons::EigenPack<LatticeColourVector>;
 
 // Noise vector (index order: nnoise, nt, nvec, ns)
-using NoiseTensor = Eigen::Tensor<Complex, 4, Eigen::RowMajor>;
 
-extern const std::string PerambTensorName;
+class NoiseTensor : public NamedTensor<Complex, 4>
-extern const std::array<std::string, 6> PerambIndexNames;
+{
-using PerambTensor = NamedTensor<SpinVector, 6, PerambTensorName, PerambIndexNames>;
+    static const std::string               Name_;
+    static const std::vector<std::string>  DefaultIndexNames_;
+    public:
+    // Default constructor (assumes tensor will be loaded from file)
+    EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NoiseTensor() : NamedTensor{Name_, DefaultIndexNames_} {}
+
+    // Construct a named tensor explicitly specifying size of each dimension
+    template<typename... IndexTypes>
+    EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NoiseTensor(Eigen::Index nNoise, Eigen::Index nT, Eigen::Index nVec, Eigen::Index nS)
+    : NamedTensor{Name_, DefaultIndexNames_, nNoise, nT, nVec, nS} {}
+};
+
+class PerambTensor : public NamedTensor<SpinVector, 6>
+{
+    static const std::string               Name_;
+    static const std::vector<std::string>  DefaultIndexNames_;
+    public:
+    // Default constructor (assumes tensor will be loaded from file)
+    EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE PerambTensor() : NamedTensor{Name_, DefaultIndexNames_} {}
+
+    // Construct a named tensor explicitly specifying size of each dimension
+    template<typename... IndexTypes>
+    EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE PerambTensor(Eigen::Index nT, Eigen::Index nVec, Eigen::Index LI, Eigen::Index nNoise, Eigen::Index nT_inv, Eigen::Index SI)
+    : NamedTensor{Name_, DefaultIndexNames_, nT, nVec, LI, nNoise, nT_inv, SI} {}
+};
 
 END_MODULE_NAMESPACE
 END_HADRONS_NAMESPACE

Hadrons/Modules/MDistil/DistilVectors.hpp

@@ -169,7 +169,7 @@ void TDistilVectors<FImpl>::setup(void)
     const int SI{ static_cast<int>( perambulator.tensor.dimension(5) ) };
     const int Nt_inv{ static_cast<int>( perambulator.tensor.dimension(4) ) };
     const int nnoise{ static_cast<int>( perambulator.tensor.dimension(3) ) };
-    assert( nnoise >= static_cast<int>( noise.dimension(0) ) && "Not enough noise vectors for perambulator" );
+    assert( nnoise >= static_cast<int>( noise.tensor.dimension(0) ) && "Not enough noise vectors for perambulator" );
     // Nvec defaults to what's in the perambulator unless overriden
     const int nvec_per{ static_cast<int>( perambulator.tensor.dimension(1) ) };
     const int nvec{Hadrons::MDistil::DistilParameters::ParameterDefault(par().nvec, nvec_per, true) };
@@ -254,7 +254,7 @@ void TDistilVectors<FImpl>::execute(void)
                                 for (int ik = dk; ik < nvec; ik += LI){
                                     for (int is = ds; is < Ns; is += SI){
                                         ExtractSliceLocal(evec3d,epack.evec[ik],0,t_inv-Ntfirst,Tdir);
-                                        tmp3d_nospin = evec3d * noise(inoise, t_inv, ik, is);
+                                        tmp3d_nospin = evec3d * noise.tensor(inoise, t_inv, ik, is);
                                         tmp3d=0;
                                         pokeSpin(tmp3d,tmp3d_nospin,is);
                                         tmp2=0;

Hadrons/Modules/MDistil/Noises.hpp

@@ -45,9 +45,9 @@ public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(NoisesPar,
                                     int, nnoise,
                                     int, nvec,
-                                    std::string, UniqueIdentifier,
                                     std::string, TI,
-                                    std::string, LI)
+                                    std::string, LI,
+                                    std::string, NoiseFileName)
 };
 
 template <typename FImpl>
@@ -113,10 +113,6 @@ void TNoises<FImpl>::execute(void)
     const int LI{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().LI, nvec, false) };
     const bool full_tdil{ TI == Nt }; \
     const bool exact_distillation{ full_tdil && LI == nvec }; \
-    std::string UniqueIdentifier{par().UniqueIdentifier};
-    if (UniqueIdentifier.empty())
-        UniqueIdentifier = getName();
-    UniqueIdentifier.append( std::to_string( vm().getTrajectory() ) );
     
     // We use our own seeds so we can specify different noises per quark
     Real rn;
@@ -126,17 +122,24 @@ void TNoises<FImpl>::execute(void)
             for (int ivec = 0; ivec < nvec; ivec++) {
                 for (int is = 0; is < Ns; is++) {
                     if (exact_distillation)
-                        noise(inoise, t, ivec, is) = 1.;
+                        noise.tensor(inoise, t, ivec, is) = 1.;
                     else{
                         random(rngSerial(),rn);
                         // We could use a greater number of complex roots of unity
                         // ... but this seems to work well
-                        noise(inoise, t, ivec, is) = (rn > 0.5) ? -1 : 1;
+                        noise.tensor(inoise, t, ivec, is) = (rn > 0.5) ? -1 : 1;
                     }
                 }
             }
         }
     }
+    if (env().getGrid()->IsBoss())
+    {
+        std::string sName {par().NoiseFileName};
+        sName.append(".");
+        sName.append(std::to_string(vm().getTrajectory()));
+        noise.write(sName.c_str());
+    }
 }
 
 END_MODULE_NAMESPACE

Hadrons/Modules/MDistil/Perambulator.cc

@@ -35,13 +35,23 @@ using namespace MDistil;
 
 template class Grid::Hadrons::MDistil::TPerambulator<FIMPL>;
 
+BEGIN_HADRONS_NAMESPACE
+
 // Global constants for distillation
 
-const std::string Grid::Hadrons::MDistil::PerambTensorName{ "Perambulator" };
-const std::array<std::string, 6> Grid::Hadrons::MDistil::PerambIndexNames{"nT", "nVec", "LI", "nNoise", "nT_inv", "SI"};
-
 #ifdef HAVE_HDF5
-extern const std::string Grid::Hadrons::NamedTensorFileExtension{".h5"};
+extern const std::string NamedTensorFileExtension{".h5"};
 #else
-extern const std::string Grid::Hadrons::NamedTensorFileExtension{".dat"};
+extern const std::string NamedTensorFileExtension{".dat"};
 #endif
+
+BEGIN_MODULE_NAMESPACE(MDistil)
+
+const std::string              NoiseTensor::Name_{"Noises"};
+const std::vector<std::string> NoiseTensor::DefaultIndexNames_{"nNoise", "nT", "nVec", "nS"};
+
+const std::string              PerambTensor::Name_{"Perambulator"};
+const std::vector<std::string> PerambTensor::DefaultIndexNames_{"nT", "nVec", "LI", "nNoise", "nT_inv", "SI"};
+
+END_MODULE_NAMESPACE
+END_HADRONS_NAMESPACE

Hadrons/Modules/MDistil/Perambulator.hpp

@@ -208,7 +208,7 @@ void TPerambulator<FImpl>::execute(void)
                                 for (int is = ds; is < Ns; is += SI)
                                 {
                                     ExtractSliceLocal(evec3d,epack.evec[ik],0,t_inv-Ntfirst,Tdir);
-                                    tmp3d_nospin = evec3d * noise(inoise, t_inv, ik, is);
+                                    tmp3d_nospin = evec3d * noise.tensor(inoise, t_inv, ik, is);
                                     tmp3d=0;
                                     pokeSpin(tmp3d,tmp3d_nospin,is);
                                     tmp2=0;
@@ -276,8 +276,6 @@ void TPerambulator<FImpl>::execute(void)
     if (grid4d->IsBoss())
     {
         std::string sPerambName {par().PerambFileName};
-        if (sPerambName.empty())
-            sPerambName = getName();
         sPerambName.append(".");
         sPerambName.append(std::to_string(vm().getTrajectory()));
         perambulator.write(sPerambName.c_str());