diff --git a/Hadrons/Distil.hpp b/Hadrons/Distil.hpp index ef2b3cbf..a5f87c98 100644 --- a/Hadrons/Distil.hpp +++ b/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) - They can be persisted to disk in tag Name_, and IndexNames are validated on load. - TODO: WHAT TO SAVE / VALIDATE ON LOAD (Override to warn instead of assert on load) - Ensemble string - Configuration number - Noise unique string - Distillation parameters + Eigen::Tensor of type Scalar_ and rank NumIndices_ (row-major order), together with a name for each index. + Index names are mutable, but tensor dimensionality is not (size of each dimension is mutable). + They can be persisted to / restored from disk, by default using tag Name. + During restore from disk, these validations are performed: + 1) Tensor dimensionality must match + 2) IndexNames are validated against current values + 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 &IndexNames_> +template class NamedTensor : Serializable { public: @@ -59,40 +59,40 @@ public: GRID_SERIALIZABLE_CLASS_MEMBERS(NamedTensor, ET, tensor, std::vector, IndexNames ); - - // Get the default index names as std::vector - std::vector DefaultIndexNames() - { - std::vector names{NumIndices_}; - for (std::size_t i = 0; i < NumIndices_; i++) - names[i] = IndexNames_[i]; - return names; - } - + + // Name of the object and Index names as set in the constructor + const std::string &Name; + const std::vector &DefaultIndexNames; + + virtual ~NamedTensor(){}; // 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 &IndexNames_) + : IndexNames{IndexNames_}, Name{Name_}, DefaultIndexNames{IndexNames_} {} // Construct a named tensor explicitly specifying size of each dimension template - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NamedTensor(Eigen::Index firstDimension, IndexTypes... otherDimensions) - : tensor(firstDimension, otherDimensions...), IndexNames{DefaultIndexNames()} + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NamedTensor(const std::string &Name_, const std::vector &IndexNames_, + 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 &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() - && std::equal( IndexNames[i].begin(), IndexNames[i].end(), IndexNames_[i].begin(), + bSame = IndexNames[i].size() == CheckNames[i].size() + && 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 void write(Writer &w, const std::string &Tag = Name_) const - { write(w, Tag, *this); } - - inline void write(const std::string &filename, const std::string &Tag = Name_) const + void write(const std::string &FileName, const std::string &Tag) const { - std::string sFileName{filename}; - sFileName.append( NamedTensorFileExtension ); - LOG(Message) << "Writing " << Name_ << " to file " << sFileName << " tag " << Tag << std::endl; - Default_Writer w( sFileName ); - write( w, Tag ); + std::string FileName_{FileName}; + FileName_.append( NamedTensorFileExtension ); + LOG(Message) << "Writing " << Name << " to file " << FileName_ << " tag " << Tag << std::endl; + Default_Writer w( FileName_ ); + write( w, Tag, *this ); } - - // Read and validate index names - template void read(Reader &r, bool bValidate = true, const std::string &Tag = Name_) + void write(const std::string &FileName) const { return write(FileName, Name); } + + // Read tensor. + // Validate: + // 1) index names (if requested) + // 2) index dimensions (if they are non-zero when called) + template void read(Reader &r, bool bValidate, const std::string &Tag) { // Grab index names and dimensions std::vector 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"); } - - inline void read (const std::string &filename, bool bValidate = true, const std::string &Tag = Name_) + template void read(Reader &r, bool bValidate = true) { read(r, bValidate, Name); } + + inline void read (const std::string &FileName, bool bValidate, const std::string &Tag) { - std::string sFileName{filename}; - sFileName.append( NamedTensorFileExtension ); - LOG(Message) << "Reading " << Name_ << " from file " << sFileName << " tag " << Tag << std::endl; - Default_Reader r(sFileName); + Default_Reader r(FileName + NamedTensorFileExtension); 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; // Noise vector (index order: nnoise, nt, nvec, ns) -using NoiseTensor = Eigen::Tensor; -extern const std::string PerambTensorName; -extern const std::array PerambIndexNames; -using PerambTensor = NamedTensor; +class NoiseTensor : public NamedTensor +{ + static const std::string Name_; + static const std::vector 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 + 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 +{ + static const std::string Name_; + static const std::vector 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 + 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 diff --git a/Hadrons/Modules/MDistil/DistilVectors.hpp b/Hadrons/Modules/MDistil/DistilVectors.hpp index 0b549913..643e7d63 100644 --- a/Hadrons/Modules/MDistil/DistilVectors.hpp +++ b/Hadrons/Modules/MDistil/DistilVectors.hpp @@ -169,7 +169,7 @@ void TDistilVectors::setup(void) const int SI{ static_cast( perambulator.tensor.dimension(5) ) }; const int Nt_inv{ static_cast( perambulator.tensor.dimension(4) ) }; const int nnoise{ static_cast( perambulator.tensor.dimension(3) ) }; - assert( nnoise >= static_cast( noise.dimension(0) ) && "Not enough noise vectors for perambulator" ); + assert( nnoise >= static_cast( 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( perambulator.tensor.dimension(1) ) }; const int nvec{Hadrons::MDistil::DistilParameters::ParameterDefault(par().nvec, nvec_per, true) }; @@ -254,7 +254,7 @@ void TDistilVectors::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; diff --git a/Hadrons/Modules/MDistil/Noises.hpp b/Hadrons/Modules/MDistil/Noises.hpp index 8b7d1f45..c17a29fa 100644 --- a/Hadrons/Modules/MDistil/Noises.hpp +++ b/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 @@ -113,10 +113,6 @@ void TNoises::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::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 diff --git a/Hadrons/Modules/MDistil/Perambulator.cc b/Hadrons/Modules/MDistil/Perambulator.cc index b141eb95..7b1d68d9 100644 --- a/Hadrons/Modules/MDistil/Perambulator.cc +++ b/Hadrons/Modules/MDistil/Perambulator.cc @@ -35,13 +35,23 @@ using namespace MDistil; template class Grid::Hadrons::MDistil::TPerambulator; +BEGIN_HADRONS_NAMESPACE + // Global constants for distillation -const std::string Grid::Hadrons::MDistil::PerambTensorName{ "Perambulator" }; -const std::array 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 NoiseTensor::DefaultIndexNames_{"nNoise", "nT", "nVec", "nS"}; + +const std::string PerambTensor::Name_{"Perambulator"}; +const std::vector PerambTensor::DefaultIndexNames_{"nT", "nVec", "LI", "nNoise", "nT_inv", "SI"}; + +END_MODULE_NAMESPACE +END_HADRONS_NAMESPACE diff --git a/Hadrons/Modules/MDistil/Perambulator.hpp b/Hadrons/Modules/MDistil/Perambulator.hpp index e2ab7e0d..17ea41a2 100644 --- a/Hadrons/Modules/MDistil/Perambulator.hpp +++ b/Hadrons/Modules/MDistil/Perambulator.hpp @@ -208,7 +208,7 @@ void TPerambulator::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::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());