diff --git a/src/Estimators/EstimatorManagerNew.cpp b/src/Estimators/EstimatorManagerNew.cpp index d80bd112b2..32a87a5869 100644 --- a/src/Estimators/EstimatorManagerNew.cpp +++ b/src/Estimators/EstimatorManagerNew.cpp @@ -55,7 +55,7 @@ bool EstimatorManagerNew::areThereListeners() const return std::any_of(operator_ests_.begin(), operator_ests_.end(), [](auto& oper_est) { return oper_est->isListenerRequired(); }); } - + template bool EstimatorManagerNew::createEstimator(EstimatorInput& input, Args&&... args) { @@ -373,7 +373,7 @@ void EstimatorManagerNew::reduceOperatorEstimators() { operator_data_sizes[iop] = operator_ests_[iop]->get_data().size(); } - // 1 larger because we put the weight in to avoid dependence of the Scalar estimators being reduced firt. + // 1 larger because we put the weight in to avoid dependence of the Scalar estimators being reduced first. size_t nops = *(std::max_element(operator_data_sizes.begin(), operator_data_sizes.end())) + 1; std::vector operator_send_buffer; std::vector operator_recv_buffer; diff --git a/src/Estimators/InputSection.h b/src/Estimators/InputSection.h index 8e6c29e4c7..15827594a5 100644 --- a/src/Estimators/InputSection.h +++ b/src/Estimators/InputSection.h @@ -18,7 +18,6 @@ #include #include "Configuration.h" -#include "OhmmsData/ParameterSet.h" #include "Containers/OhmmsPETE/TinyVector.h" #include "Message/UniformCommunicateError.h" @@ -44,12 +43,12 @@ class InputSection // Becuase it hurts to read all the trailing _ in the constructors of input section subtypes // NOLINTBEGIN(readability-indentifier-naming) - std::string section_name; // name of the input section + std::string section_name; // name of the input section - std::unordered_set attributes; // list of attribute variables - std::unordered_set parameters; // list of parameter variables - std::unordered_set delegates; // input nodes delegate to next level of input parsing. - std::unordered_set required; // list of required variables + std::unordered_set attributes; // list of attribute variables + std::unordered_set parameters; // list of parameter variables + std::unordered_set delegates; // input nodes delegate to next level of input parsing. + std::unordered_set required; // list of required variables std::unordered_set strings; // list of string variables that can have one value std::unordered_set multi_strings; // list of string variables that can one or more values @@ -146,8 +145,7 @@ class InputSection * \param[in] tag parmater name or node ename delgation is controlled by * \param[in] delegate_handler factory function for delegated input function. */ - void registerDelegate(const std::string& tag, - DelegateHandler delegate_handler); + void registerDelegate(const std::string& tag, DelegateHandler delegate_handler); /** Do validation for a particular subtype of InputSection * Called by check_valid. @@ -176,7 +174,9 @@ class InputSection /** Derived class can overrides this to do custom parsing of the element values for Custom elements * These can have a name attribute only. */ - [[noreturn]] virtual void setFromStreamCustom(const std::string& ename, const std::string& name, std::istringstream& svalue) + [[noreturn]] virtual void setFromStreamCustom(const std::string& ename, + const std::string& name, + std::istringstream& svalue) { throw std::runtime_error("derived class must provide handleCustom method if custom parameters are used"); } diff --git a/src/Estimators/OperatorEstBase.h b/src/Estimators/OperatorEstBase.h index f59ec16926..8aa99f2f14 100644 --- a/src/Estimators/OperatorEstBase.h +++ b/src/Estimators/OperatorEstBase.h @@ -102,11 +102,13 @@ class OperatorEstBase virtual std::unique_ptr spawnCrowdClone() const = 0; /** Write to previously registered observable_helper hdf5 wrapper. + * + * You need to override this if your estimator does anything more than just write its data to a single h5 descriptor. * * if you haven't registered Operator Estimator - * this will do nothing. + * this will/should do nothing. */ - void write(hdf_archive& file); + virtual void write(hdf_archive& file); /** zero data appropriately for the DataLocality */ diff --git a/src/Estimators/SpinDensityNew.cpp b/src/Estimators/SpinDensityNew.cpp index 11d1cc18c9..1741b2096f 100644 --- a/src/Estimators/SpinDensityNew.cpp +++ b/src/Estimators/SpinDensityNew.cpp @@ -130,11 +130,25 @@ void SpinDensityNew::accumulate(const RefVector& walkers, walkers_weight_ += weight; int p = 0; size_t offset = 0; + QMCT::PosType u; for (int s = 0; s < species_.size(); ++s, offset += dp_.npoints) for (int ps = 0; ps < species_size_[s]; ++ps, ++p) { - QMCT::PosType u = lattice_.toUnit(pset.R[p] - dp_.corner); - size_t point = offset; + switch (pset.R.getUnit()) + { + case PosUnit::Cartesian: + u = lattice_.toUnit(pset.R[p] - dp_.corner); + break; + case PosUnit::Lattice: +#ifndef NDEBUG + if (!(p < pset.getTotalNum())) + throw std::runtime_error("p < pset.getTotalNum(): " + std::to_string(p) + " " + + std::to_string(pset.getTotalNum())); +#endif + u = pset.R[p] - lattice_.toUnit(dp_.corner); + break; + } + size_t point = offset; for (int d = 0; d < QMCT::DIM; ++d) point += dp_.gdims[d] * ((int)(dp_.grid[d] * (u[d] - std::floor(u[d])))); //periodic only accumulateToData(point, weight); @@ -195,41 +209,64 @@ void SpinDensityNew::collect(const RefVector& type_erased_opera } } -void SpinDensityNew::report(const std::string& pad) +void SpinDensityNew::report(const std::string& pad) const { report(pad, app_log()); } + +void SpinDensityNew::report(const std::string& pad, std::ostream& out) const { auto& dp_ = derived_parameters_; - app_log() << pad << "SpinDensity report" << std::endl; - app_log() << pad << " dim = " << QMCT::DIM << std::endl; - app_log() << pad << " npoints = " << dp_.npoints << std::endl; - app_log() << pad << " grid = " << dp_.grid << std::endl; - app_log() << pad << " gdims = " << dp_.gdims << std::endl; - app_log() << pad << " corner = " << dp_.corner << std::endl; - app_log() << pad << " center = " << dp_.corner + lattice_.Center << std::endl; - app_log() << pad << " cell " << std::endl; + out << pad << "SpinDensity report" << std::endl; + out << pad << " dim = " << QMCT::DIM << std::endl; + out << pad << " npoints = " << dp_.npoints << std::endl; + out << pad << " grid = " << dp_.grid << std::endl; + out << pad << " gdims = " << dp_.gdims << std::endl; + out << pad << " corner = " << dp_.corner << std::endl; + out << pad << " center = " << dp_.corner + lattice_.Center << std::endl; + out << pad << " cell " << std::endl; for (int d = 0; d < QMCT::DIM; ++d) - app_log() << pad << " " << d << " " << lattice_.Rv[d] << std::endl; - app_log() << pad << " end cell " << std::endl; - app_log() << pad << " nspecies = " << species_.size() << std::endl; + out << pad << " " << d << " " << lattice_.Rv[d] << std::endl; + out << pad << " end cell " << std::endl; + out << pad << " nspecies = " << species_.size() << std::endl; for (int s = 0; s < species_.size(); ++s) - app_log() << pad << " species[" << s << "]" - << " = " << species_.speciesName[s] << " " << species_size_[s] << std::endl; - app_log() << pad << "end SpinDensity report" << std::endl; + out << pad << " species[" << s << "]" + << " = " << species_.speciesName[s] << " " << species_size_[s] << std::endl; + out << pad << "end SpinDensity report" << std::endl; } -void SpinDensityNew::registerOperatorEstimator(hdf_archive& file) + +void SpinDensityNew::registerOperatorEstimator(hdf_archive& file) {} + +void SpinDensityNew::write(hdf_archive& file) { - std::vector my_indexes; + auto writeFullPrecData = [&](auto& file, auto& fp_data) { + std::vector my_indexes; + std::array ng{1, derived_parameters_.npoints}; + hdf_path hdf_name{my_name_}; + file.push(hdf_name); + Vector data; + std::size_t offset = 0; + for (int s = 0; s < species_.size(); ++s) + { + data.attachReference(fp_data.data() + offset, derived_parameters_.npoints); + hid_t g_id = file.push(species_.speciesName[s], true); + file.append(data, "value", append_index_); + offset += derived_parameters_.npoints; + file.pop(); + } + ++append_index_; + file.pop(); + }; - std::vector ng(1, derived_parameters_.npoints); +#ifdef MIXED_PRECISION + std::vector expanded_data(data_.size(), 0.0); + std::copy_n(data_.begin(), data_.size(), expanded_data.begin()); + assert(!data_.empty()); + writeFullPrecData(file, expanded_data); + // auto total = std::accumulate(data_->begin(), data_->end(), 0.0); + // std::cout << "data size: " << data_->size() << " : " << total << '\n'; - hdf_path hdf_name{my_name_}; - for (int s = 0; s < species_.size(); ++s) - { - h5desc_.emplace_back(hdf_name / species_.speciesName[s]); - auto& oh = h5desc_.back(); - oh.set_dimensions(ng, 0); - } +#else + writeFullPrecData(file, data_); +#endif } - } // namespace qmcplusplus diff --git a/src/Estimators/SpinDensityNew.h b/src/Estimators/SpinDensityNew.h index a74d3189bb..15d16ba24f 100644 --- a/src/Estimators/SpinDensityNew.h +++ b/src/Estimators/SpinDensityNew.h @@ -109,6 +109,10 @@ class SpinDensityNew : public OperatorEstBase */ void registerOperatorEstimator(hdf_archive& file) override; + /** Custom write method since spin density new has its own hdf5 format + */ + void write(hdf_archive& file) override; + private: SpinDensityNew(const SpinDensityNew& sdn) = default; @@ -118,7 +122,8 @@ class SpinDensityNew : public OperatorEstBase size_t getFullDataSize(); void accumulateToData(size_t point, QMCT::RealType weight); void reset(); - void report(const std::string& pad); + void report(const std::string& pad) const; + void report(const std::string& pad, std::ostream& out) const; //data members const SpinDensityInput input_; @@ -138,6 +143,9 @@ class SpinDensityNew : public OperatorEstBase SpinDensityInput::DerivedParameters derived_parameters_; /**}@*/ + /// current index for h5d_append + hsize_t append_index_ = 0; + friend class testing::SpinDensityNewTests; }; diff --git a/src/Estimators/tests/test_SpinDensityNew.cpp b/src/Estimators/tests/test_SpinDensityNew.cpp index d357104ab4..2a0122c594 100644 --- a/src/Estimators/tests/test_SpinDensityNew.cpp +++ b/src/Estimators/tests/test_SpinDensityNew.cpp @@ -2,7 +2,7 @@ // This file is distributed under the University of Illinois/NCSA Open Source License. // See LICENSE file in top directory for details. // -// Copyright (c) 2021 QMCPACK developers. +// Copyright (c) 2023 QMCPACK developers. // // File developed by: Peter Doak, doakpw@ornl.gov, Oak Ridge National Lab // @@ -15,10 +15,13 @@ #include "SpinDensityInput.h" #include "ValidSpinDensityInput.h" #include "SpinDensityNew.h" +#include "hdf/hdf_archive.h" +#include "hdf/hdf_path.h" #include "RandomForTest.h" #include "ParticleSet.h" #include "TrialWaveFunction.h" #include "EstimatorTesting.h" +#include "NativeInitializerPrint.hpp" #include "OhmmsData/Libxml2Doc.h" @@ -28,11 +31,13 @@ namespace qmcplusplus { -using QMCT = QMCTraits; +using QMCT = QMCTraits; +using Real = QMCTraits::RealType; +using FullPrecReal = QMCTraits::FullPrecRealType; namespace testing { -/** class to preserve access control in MomentumDistribution +/** class to preserve access control in SpinDensityNew */ class SpinDensityNewTests { @@ -44,32 +49,67 @@ class SpinDensityNewTests CHECK(sdn.species_size_ == sdn2.species_size_); CHECK(sdn.data_ != sdn2.data_); } + + const std::string& getName(const SpinDensityNew& sdn) { return sdn.my_name_; } + const SpinDensityInput::DerivedParameters& getDerivedParameters(const SpinDensityNew& sdn) + { + return sdn.derived_parameters_; + } + static void report(const std::string& padding, const SpinDensityNew& sdn) { sdn.report(padding, std::cout); } }; } // namespace testing - -void accumulateFromPsets(int ncrowds, SpinDensityNew& sdn, UPtrVector& crowd_sdns) +void accumulateFromPsets(int n_crowds, + int n_walkers, + SpinDensityNew& sdn, + UPtrVector& crowd_sdns, + int n_electrons, + PosUnit pos_unit = PosUnit::Lattice) { - const SimulationCell simulation_cell; - for (int iops = 0; iops < ncrowds; ++iops) + // aka PtclOnLatticeTraits::ParticleLayout CrystalLattice + CrystalLattice lattice; + lattice.set(Tensor(3.37316115, 3.37316115, 0.00000000, 0.00000000, 3.37316115, 3.37316115, + 3.37316115, 0.00000000, 3.37316115)); + const SimulationCell simulation_cell(lattice); + for (int iops = 0; iops < n_crowds; ++iops) { std::vector walkers; - int nwalkers = 4; - for (int iw = 0; iw < nwalkers; ++iw) - walkers.emplace_back(2); + for (int iw = 0; iw < n_walkers; ++iw) + walkers.emplace_back(n_electrons); std::vector psets; crowd_sdns.emplace_back(sdn.spawnCrowdClone()); SpinDensityNew& crowd_sdn = dynamic_cast(*(crowd_sdns.back())); - for (int iw = 0; iw < nwalkers; ++iw) + for (int iw = 0; iw < n_walkers; ++iw) { psets.emplace_back(simulation_cell); ParticleSet& pset = psets.back(); - pset.create({2}); - pset.R[0] = ParticleSet::PosType(0.00000000, 0.00000000, 0.00000000); - pset.R[1] = ParticleSet::PosType(0.68658058, 0.68658058, 0.68658058); + pset.create({n_electrons}); + pset.R.setUnit(pos_unit); + Real inv_ne = 1 / static_cast(n_electrons); + const auto& lattice = simulation_cell.getLattice(); + std::cout << "inv_ne: " << inv_ne << '\n'; + // There is a strong assumption that the electron coordinates will only be from -0.5 to 0.5 in lattice coords + for (int ie = 0; ie < n_electrons; ++ie) + { + ParticleSet::SingleParticlePos pos(ie * inv_ne, ie * inv_ne, ie * inv_ne); + std::cout << "position" << pos[0] << " " << pos[1] << " " << pos[2] << '\n'; + ParticleSet::SingleParticlePos corner{-0.5, -0.5, -0.5}; + ParticleSet::SingleParticlePos adjusted_pos = pos - corner; + switch (pos_unit) + { + case PosUnit::Lattice: + pset.R[ie] = adjusted_pos; + break; + case PosUnit::Cartesian: + std::cout << "Cartesian R: " << lattice.R << '\n'; + pset.R[ie] = lattice.toCart(adjusted_pos); + std::cout << "cart position:" << pset.R[ie][0] << " " << pset.R[ie][1] << " " << pset.R[ie][2] << '\n'; + break; + } + } } std::vector wfns; @@ -84,30 +124,37 @@ void accumulateFromPsets(int ncrowds, SpinDensityNew& sdn, UPtrVector& rft, UPtrVector& crowd_sdns) +void randomUpdateAccumulate(int n_walkers, + testing::RandomForTest& rft, + UPtrVector& crowd_sdns, + int n_elecs) { - const SimulationCell simulation_cell; + // aka PtclOnLatticeTraits::ParticleLayout CrystalLattice + CrystalLattice lattice; + lattice.set(Tensor(3.37316115, 3.37316115, 0.00000000, 0.00000000, 3.37316115, 3.37316115, + 3.37316115, 0.00000000, 3.37316115)); + + const SimulationCell simulation_cell(lattice); for (auto& uptr_crowd_sdn : crowd_sdns) { std::vector walkers; - int nwalkers = 4; - for (int iw = 0; iw < nwalkers; ++iw) - walkers.emplace_back(2); + for (int iw = 0; iw < n_walkers; ++iw) + walkers.emplace_back(n_elecs); std::vector psets; SpinDensityNew& crowd_sdn = dynamic_cast(*(uptr_crowd_sdn)); - std::vector rng_reals(nwalkers * QMCT::DIM * 2); + std::vector rng_reals(n_walkers * QMCT::DIM * n_elecs); rft.fillVecRng(rng_reals); auto it_rng_reals = rng_reals.begin(); - for (int iw = 0; iw < nwalkers; ++iw) + for (int iw = 0; iw < n_walkers; ++iw) { psets.emplace_back(simulation_cell); ParticleSet& pset = psets.back(); - pset.create({2}); - pset.R[0] = ParticleSet::PosType(*it_rng_reals++, *it_rng_reals++, *it_rng_reals++); - pset.R[1] = ParticleSet::PosType(*it_rng_reals++, *it_rng_reals++, *it_rng_reals++); + pset.create({n_elecs}); + for (int i = 0; i < n_elecs; ++i) + pset.R[i] = ParticleSet::PosType(*it_rng_reals++, *it_rng_reals++, *it_rng_reals++); } std::vector wfns; @@ -227,13 +274,11 @@ TEST_CASE("SpinDensityNew::accumulate", "[estimators]") std::vector& data_ref = sdn.get_data(); // There should be a check that the discretization of particle locations expressed in lattice coords // is correct. This just checks it hasn't changed from how it was in SpinDensity which lacked testing. - CHECK(data_ref[555] == 4); - CHECK(data_ref[1777] == 4); } TEST_CASE("SpinDensityNew::collect(DataLocality::crowd)", "[estimators]") { - { + auto checkCollectForNParticles = [](int n_elec) { using MCPWalker = OperatorEstBase::MCPWalker; using QMCT = QMCTraits; @@ -246,16 +291,17 @@ TEST_CASE("SpinDensityNew::collect(DataLocality::crowd)", "[estimators]") int ispecies = species_set.addSpecies("u"); ispecies = species_set.addSpecies("d"); CHECK(ispecies == 1); - int iattribute = species_set.addAttribute("membersize"); - species_set(iattribute, 0) = 1; - species_set(iattribute, 1) = 1; + int i_msize = species_set.addAttribute("membersize"); + species_set(i_msize, 0) = n_elec / 2 + n_elec % 2; + species_set(i_msize, 1) = n_elec / 2; SpinDensityNew sdn(std::move(sdi), species_set); UPtrVector crowd_sdns; - int ncrowds = 2; + int n_crowds = 2; + int n_walkers = 4; - accumulateFromPsets(ncrowds, sdn, crowd_sdns); + accumulateFromPsets(n_crowds, n_walkers, sdn, crowd_sdns, n_elec); RefVector crowd_oeb_refs = convertUPtrToRefVector(crowd_sdns); sdn.collect(crowd_oeb_refs); @@ -263,17 +309,28 @@ TEST_CASE("SpinDensityNew::collect(DataLocality::crowd)", "[estimators]") std::vector& data_ref = sdn.get_data(); // There should be a check that the discretization of particle locations expressed in lattice coords // is correct. This just checks it hasn't changed from how it was in SpinDensity which lacked testing. - CHECK(data_ref[555] == 4 * ncrowds); - CHECK(data_ref[1666] == 4 * ncrowds); - } + Real up_sum{0}; + Real dn_sum{0}; + for (int i = 0; i < 1000; ++i) + up_sum += data_ref[i]; + for (int i = 0; i < 1000; ++i) + dn_sum += data_ref[i + 1000]; + // \todo why is it 8? + CHECK(up_sum == n_crowds * n_walkers * species_set(i_msize, 0)); + CHECK(dn_sum == n_crowds * n_walkers * species_set(i_msize, 1)); + }; + checkCollectForNParticles(1); + checkCollectForNParticles(2); + checkCollectForNParticles(3); + checkCollectForNParticles(11); } TEST_CASE("SpinDensityNew::collect(DataLocality::rank)", "[estimators]") { - { - using MCPWalker = OperatorEstBase::MCPWalker; - using QMCT = QMCTraits; + using MCPWalker = OperatorEstBase::MCPWalker; + using QMCT = QMCTraits; + auto checkCollectAtPoint = [](int n_elec, PosUnit pos_unit) { Libxml2Document doc; bool okay = doc.parseFromString(testing::valid_spin_density_input_sections[0]); REQUIRE(okay); @@ -283,18 +340,26 @@ TEST_CASE("SpinDensityNew::collect(DataLocality::rank)", "[estimators]") int ispecies = species_set.addSpecies("u"); ispecies = species_set.addSpecies("d"); CHECK(ispecies == 1); - int iattribute = species_set.addAttribute("membersize"); - species_set(iattribute, 0) = 1; - species_set(iattribute, 1) = 1; + int i_msize = species_set.addAttribute("membersize"); + std::vector species_size(species_set.size()); + species_size[0] = n_elec / 2 + n_elec % 2; + species_size[1] = n_elec / 2; + species_set(i_msize, 0) = species_size[0]; + species_set(i_msize, 1) = species_size[1]; + SpinDensityNew sdn(std::move(sdi), species_set, DataLocality::rank); + testing::SpinDensityNewTests sdnt; + const auto& derived_parameters = sdnt.getDerivedParameters(sdn); + auto lattice = testing::makeTestLattice(); UPtrVector crowd_sdns; - int ncrowds = 2; + int n_crowds = 2; + int n_walkers = 4; - accumulateFromPsets(ncrowds, sdn, crowd_sdns); + accumulateFromPsets(n_crowds, n_walkers, sdn, crowd_sdns, n_elec, pos_unit); RefVector crowd_oeb_refs = convertUPtrToRefVector(crowd_sdns); sdn.collect(crowd_oeb_refs); @@ -302,9 +367,38 @@ TEST_CASE("SpinDensityNew::collect(DataLocality::rank)", "[estimators]") std::vector& data_ref = sdn.get_data(); // There should be a check that the discretization of particle locations expressed in lattice coords // is correct. This just checks it hasn't changed from how it was in SpinDensity which lacked testing. - CHECK(data_ref[555] == 4 * ncrowds); - CHECK(data_ref[1666] == 4 * ncrowds); - } + // + INFO("n_elec: " << n_elec); + CHECK(data_ref.size() == 2000); + sdnt.report("", sdn); + std::cout << NativePrint(data_ref) << '\n'; + + + size_t offset = 0; + int p = 0; + Real inv_ne = 1 / static_cast(n_elec); + for (int s = 0; s < species_set.size(); ++s, offset += derived_parameters.npoints) + for (int ps = 0; ps < species_size[s]; ++ps, ++p) + { + //QMCT::PosType u = lattice_.toUnit(pset.R[p] - dp_.corner); + ParticleSet::SingleParticlePos pos(p * inv_ne, p * inv_ne, p * inv_ne); + std::cout << "position" << pos[0] << " " << pos[1] << " " << pos[2] << '\n'; + //pos = lattice.toCart(pos); + //std::cout << "position" << pos[0] << " " << pos[1] << " " << pos[2] <<'\n'; + + //pset.R[ie] = pos; + size_t point = offset; + for (int d = 0; d < QMCT::DIM; ++d) + point += derived_parameters.gdims[d] * + ((int)(derived_parameters.grid[d] * (pos[d] - std::floor(pos[d])))); //periodic only + INFO("species:" << s << " ps: " << ps << " point:" << point); + CHECK(data_ref[point] == n_walkers * n_crowds); + } + }; + checkCollectAtPoint(2, PosUnit::Lattice); + checkCollectAtPoint(3, PosUnit::Lattice); + checkCollectAtPoint(2, PosUnit::Cartesian); + checkCollectAtPoint(3, PosUnit::Cartesian); } TEST_CASE("SpinDensityNew algorithm comparison", "[estimators]") @@ -312,51 +406,145 @@ TEST_CASE("SpinDensityNew algorithm comparison", "[estimators]") using MCPWalker = OperatorEstBase::MCPWalker; using QMCT = QMCTraits; - Libxml2Document doc; - bool okay = doc.parseFromString(testing::valid_spin_density_input_sections[0]); - REQUIRE(okay); - xmlNodePtr node = doc.getRoot(); - SpinDensityInput sdi(node); - SpeciesSet species_set; - int ispecies = species_set.addSpecies("u"); - ispecies = species_set.addSpecies("d"); - CHECK(ispecies == 1); - int iattribute = species_set.addAttribute("membersize"); - species_set(iattribute, 0) = 1; - species_set(iattribute, 1) = 1; + auto checkAlgorithmsForNParticles = [](int n_elec) { + Libxml2Document doc; + bool okay = doc.parseFromString(testing::valid_spin_density_input_sections[0]); + REQUIRE(okay); + xmlNodePtr node = doc.getRoot(); + SpinDensityInput sdi(node); + SpeciesSet species_set; + int ispecies = species_set.addSpecies("u"); + ispecies = species_set.addSpecies("d"); + CHECK(ispecies == 1); + int i_msize = species_set.addAttribute("membersize"); + species_set(i_msize, 0) = n_elec / 2 + n_elec % 2; + species_set(i_msize, 1) = n_elec / 2; + + SpinDensityInput sdi_copy = sdi; + + int n_crowds = 3; + int n_walkers = 4; + int nsteps = 4; + + SpinDensityNew sdn_rank(std::move(sdi), species_set, DataLocality::rank); + UPtrVector crowd_sdns_rank; + accumulateFromPsets(n_crowds, n_walkers, sdn_rank, crowd_sdns_rank, n_elec); + testing::RandomForTest rng_for_test_rank; + for (int i = 0; i < nsteps; ++i) + randomUpdateAccumulate(n_walkers, rng_for_test_rank, crowd_sdns_rank, n_elec); + RefVector crowd_oeb_refs_rank = convertUPtrToRefVector(crowd_sdns_rank); + sdn_rank.collect(crowd_oeb_refs_rank); + std::vector& data_ref_rank = sdn_rank.get_data(); + SpinDensityNew sdn_crowd(std::move(sdi), species_set, DataLocality::crowd); + UPtrVector crowd_sdns_crowd; + accumulateFromPsets(n_crowds, n_walkers, sdn_crowd, crowd_sdns_crowd, n_elec); + testing::RandomForTest rng_for_test_crowd; + for (int i = 0; i < nsteps; ++i) + randomUpdateAccumulate(n_walkers, rng_for_test_crowd, crowd_sdns_crowd, n_elec); + RefVector crowd_oeb_refs_crowd = convertUPtrToRefVector(crowd_sdns_crowd); + sdn_crowd.collect(crowd_oeb_refs_crowd); + std::vector& data_ref_crowd = sdn_crowd.get_data(); + + for (size_t i = 0; i < data_ref_rank.size(); ++i) + { + if (data_ref_crowd[i] != data_ref_rank[i]) + FAIL_CHECK("crowd local " << data_ref_crowd[i] << " != rank local " << data_ref_rank[i] << " at index " << i); + break; + } - SpinDensityInput sdi_copy = sdi; - - int ncrowds = 3; - int nsteps = 4; - - SpinDensityNew sdn_rank(std::move(sdi), species_set, DataLocality::rank); - UPtrVector crowd_sdns_rank; - accumulateFromPsets(ncrowds, sdn_rank, crowd_sdns_rank); - testing::RandomForTest rng_for_test_rank; - for (int i = 0; i < nsteps; ++i) - randomUpdateAccumulate(rng_for_test_rank, crowd_sdns_rank); - RefVector crowd_oeb_refs_rank = convertUPtrToRefVector(crowd_sdns_rank); - sdn_rank.collect(crowd_oeb_refs_rank); - std::vector& data_ref_rank = sdn_rank.get_data(); - - SpinDensityNew sdn_crowd(std::move(sdi), species_set, DataLocality::crowd); - UPtrVector crowd_sdns_crowd; - accumulateFromPsets(ncrowds, sdn_crowd, crowd_sdns_crowd); - testing::RandomForTest rng_for_test_crowd; - for (int i = 0; i < nsteps; ++i) - randomUpdateAccumulate(rng_for_test_crowd, crowd_sdns_crowd); - RefVector crowd_oeb_refs_crowd = convertUPtrToRefVector(crowd_sdns_crowd); - sdn_crowd.collect(crowd_oeb_refs_crowd); - std::vector& data_ref_crowd = sdn_crowd.get_data(); - - for (size_t i = 0; i < data_ref_rank.size(); ++i) - { - if (data_ref_crowd[i] != data_ref_rank[i]) - FAIL_CHECK("crowd local " << data_ref_crowd[i] << " != rank local " << data_ref_rank[i] << " at index " << i); - break; - } + auto sumAndCheck = [i_msize, nsteps, n_crowds, n_walkers, &species_set](const auto& data_ref) { + Real up_sum{0}; + Real dn_sum{0}; + for (int i = 0; i < 1000; ++i) + up_sum += data_ref[i]; + for (int i = 0; i < 1000; ++i) + dn_sum += data_ref[i + 1000]; + // +1 is for the accumulateFromPsets + CHECK(up_sum == (nsteps + 1) * n_crowds * n_walkers * species_set(i_msize, 0)); + CHECK(dn_sum == (nsteps + 1) * n_crowds * n_walkers * species_set(i_msize, 1)); + }; + sumAndCheck(data_ref_rank); + sumAndCheck(data_ref_crowd); + }; + + checkAlgorithmsForNParticles(2); + checkAlgorithmsForNParticles(3); } +TEST_CASE("SpinDensityNew::registerAndWrite", "[estimators]") +{ + auto checkWriteForNParticles = [](int n_elec) { + using MCPWalker = OperatorEstBase::MCPWalker; + using QMCT = QMCTraits; + + Libxml2Document doc; + bool okay = doc.parseFromString(testing::valid_spin_density_input_sections[0]); + REQUIRE(okay); + xmlNodePtr node = doc.getRoot(); + SpinDensityInput sdi(node); + SpeciesSet species_set; + int ispecies = species_set.addSpecies("u"); + ispecies = species_set.addSpecies("d"); + CHECK(ispecies == 1); + int i_msize = species_set.addAttribute("membersize"); + species_set(i_msize, 0) = n_elec / 2 + n_elec % 2; + species_set(i_msize, 1) = n_elec / 2; + + SpinDensityNew sdn(std::move(sdi), species_set); + + UPtrVector crowd_sdns; + int n_crowds = 2; + int n_walkers = 4; + + accumulateFromPsets(n_crowds, n_walkers, sdn, crowd_sdns, n_elec); + + RefVector crowd_oeb_refs = convertUPtrToRefVector(crowd_sdns); + sdn.collect(crowd_oeb_refs); + + hdf_archive hd; + std::string test_file{"spin_density_new_test.hdf"}; + okay = hd.create(test_file); + REQUIRE(okay); + + sdn.registerOperatorEstimator(hd); + sdn.write(hd); + hd.close(); + + hdf_archive read_hd; + okay = read_hd.open(test_file); + REQUIRE(okay); + std::vector up_data; + std::vector dn_data; + + testing::SpinDensityNewTests sdnt; + hdf_path sdn_name_up{sdnt.getName(sdn)}; + + sdn_name_up /= "u/value"; + std::cout << sdn_name_up.string() << '\n'; + std::array select_one_grid{0, -1}; + read_hd.readSlabSelection(up_data, select_one_grid, sdn_name_up.string()); //sdn_name_up.string()); + CHECK(up_data.size() == 1000); + hdf_path sdn_name_dn{sdnt.getName(sdn)}; + sdn_name_dn /= "d/value"; + std::cout << sdn_name_dn.string() << '\n'; + read_hd.readSlabSelection(dn_data, select_one_grid, sdn_name_dn.string()); + CHECK(up_data.size() == 1000); + + auto sumAndCheckData = [n_crowds, n_walkers](const auto& data, int num_particles) { + FullPrecReal sum = std::accumulate(data.begin(), data.end(), 0.0); + // +1 is for the accumulateFromPsets + INFO("num_particles: " << num_particles); + CHECK(sum == n_crowds * n_walkers * num_particles); + }; + + std::cout << NativePrint(up_data); + std::cout << NativePrint(dn_data); + + sumAndCheckData(up_data, species_set(i_msize, 0)); + sumAndCheckData(dn_data, species_set(i_msize, 1)); + }; + checkWriteForNParticles(2); + checkWriteForNParticles(3); +} } // namespace qmcplusplus diff --git a/src/Particle/ParticleSet.h b/src/Particle/ParticleSet.h index a769a371c3..1970b38a65 100644 --- a/src/Particle/ParticleSet.h +++ b/src/Particle/ParticleSet.h @@ -72,7 +72,13 @@ class ParticleSet : public QMCTraits, public OhmmsElementBase, public PtclOnLatt quantum_domains quantum_domain; //@{ public data members - ///Species ID + /** Species ID + * This forms a redundant and unchecked for consistency pair with the species_set_ membersize attribute. + * This can result in out of bounds access depending on which the client of ParticleSet trusts.! + * Group id is more widely used but implementation specific and a per particle value. + * the membersize attribute is a per group value and allows code to process by group. + * \todo Fix this mess. + */ ParticleIndex GroupID; ///Position ParticlePos R; diff --git a/src/io/hdf/hdf_archive.h b/src/io/hdf/hdf_archive.h index 651e9deadd..f0b407a576 100644 --- a/src/io/hdf/hdf_archive.h +++ b/src/io/hdf/hdf_archive.h @@ -43,6 +43,8 @@ class Communicate; namespace qmcplusplus { /** class to handle hdf file + * + * Wrapper functions necessary because hdf5 c interface h5d_xxx semantics. */ class hdf_archive { @@ -247,6 +249,45 @@ class hdf_archive } } + /** Append data to a dynamically sized dataspace. + * + * \param[in] data reference to actual data object, hdf5_proxy, hdf_stl, hdf_pete provide proxies for most. + * \param[in] aname hdf5 path from the current file group on. + * \param[in] append_index In a dynamic dataseries which entry will this be, updated to next entry. + * + * This has no legacy usage and exceptions are thrown in case of error. + * + * This is pointless unless you do something with the current_append_index + * which is the necessary state to determine where in the data series you + * are writing. + */ + template + unsigned long long append(T& data, const std::string& aname, const unsigned long long current_append_index) + { + auto local_append_index = current_append_index; + if (!appendEntry(data, aname, local_append_index)) + { + throw std::runtime_error("HDF5 write failure in hdf_archive::write " + aname); + } + return local_append_index; + } + + /** write the data to the group aname and return status + * use write() for inbuilt error checking + * @return true if successful + */ + template + bool appendEntry(T& data, const std::string& aname, unsigned long long& current_append_index) + { + if (Mode[NOIO]) + return true; + if (!(Mode[IS_PARALLEL] || Mode[IS_MASTER])) + throw std::runtime_error("Only write data in parallel or by master but not every rank!"); + hid_t p = group_id.empty() ? file_id : group_id.top(); + h5data_proxy::type> e(data); + return e.append(data, p, aname, current_append_index, xfer_plist); + } + /** write the container data with a specific shape and check status * @param data container, linear storage required. * @param shape shape on the hdf file diff --git a/src/io/hdf/hdf_dataproxy.h b/src/io/hdf/hdf_dataproxy.h index fce4ac36bd..abb9e5b491 100644 --- a/src/io/hdf/hdf_dataproxy.h +++ b/src/io/hdf/hdf_dataproxy.h @@ -48,6 +48,16 @@ struct h5data_proxy : public h5_space_type { return h5d_write(grp, aname.c_str(), FileSpace::rank, dims, get_address(&ref), xfer_plist); } + + inline bool append(const data_type& ref, hid_t grp, const std::string& aname, hsize_t& current_leading_index, hid_t xfer_plist = H5P_DEFAULT) const + { + std::vector my_dims; + hsize_t rank = dims.size() + 1; + my_dims.resize(rank, 1); + copy(dims.begin(), dims.end(), my_dims.begin() + 1); + return h5d_append(grp, aname.c_str(), current_leading_index, FileSpace::rank, my_dims.data(), get_address(&ref), xfer_plist); + } + }; /** specialization for bool, convert to int diff --git a/src/io/hdf/hdf_pete.h b/src/io/hdf/hdf_pete.h index 55a2eb3001..92e53acc87 100644 --- a/src/io/hdf/hdf_pete.h +++ b/src/io/hdf/hdf_pete.h @@ -45,6 +45,12 @@ struct h5data_proxy> : public h5_space_type { return h5d_write(grp, aname.c_str(), FileSpace::rank, dims, get_address(ref.data()), xfer_plist); } + + inline bool append(const data_type& ref, hid_t grp, const std::string& aname, hsize_t& current_append_index, hid_t xfer_plist = H5P_DEFAULT) + { + std::array my_dims{1,dims[0]}; + return h5d_append(grp, aname.c_str(), current_append_index, 2, my_dims.data(), get_address(ref.data()), 1, xfer_plist); + } }; diff --git a/src/io/hdf/tests/test_hdf_archive.cpp b/src/io/hdf/tests/test_hdf_archive.cpp index f9410d91dd..8e2832a53c 100644 --- a/src/io/hdf/tests/test_hdf_archive.cpp +++ b/src/io/hdf/tests/test_hdf_archive.cpp @@ -199,6 +199,41 @@ TEST_CASE("hdf_archive_vector", "[hdf]") CHECK(v_const[i] == v2_for_const[i]); } +TEST_CASE("hdf_archive_append_pete", "[hdf]") +{ + using qmcplusplus::Vector; + hdf_archive hd; + hd.create("test_append_pete.hdf"); + + Vector v1{2.3, 3.4, 5.6}; + hsize_t append_index = 0; + + std::string name = "pete_vector_series"; + hd.append(v1, name, append_index); + + Vector v2{4.0, 5.0, 6.0}; + + hd.append(v2, name, append_index); + hd.close(); + + hdf_archive read_hd; + bool okay = read_hd.open("test_append_pete.hdf"); + REQUIRE(okay); + + Vector v_read; + std::array select_one_vector{0, -1}; + read_hd.readSlabSelection(v_read, select_one_vector, name); + + CHECK(v_read.size() == 3); + CHECK(v1 == v_read); + + select_one_vector = {1, -1}; + read_hd.readSlabSelection(v_read, select_one_vector, name); + + CHECK(v_read.size() == 3); + CHECK(v2 == v_read); +} + TEST_CASE("hdf_archive_group", "[hdf]") { hdf_archive hd;