diff --git a/_codeql_detected_source_root b/_codeql_detected_source_root new file mode 120000 index 000000000..945c9b46d --- /dev/null +++ b/_codeql_detected_source_root @@ -0,0 +1 @@ +. \ No newline at end of file diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index f40875985..552cd0a15 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -21,6 +21,7 @@ #include "BiorthBess.H" #include "BasisFactory.H" #include "BiorthCube.H" +#include "Covariance.H" #include "sltableMP2.H" #include "SLGridMP2.H" #include "YamlCheck.H" @@ -106,22 +107,14 @@ namespace BasisClasses //! Coefficient variance computation enabled bool pcavar = false; - //! Data type for covariance - bool floatType = false; - //@{ //! Sample counts and masses for covariance computation Eigen::VectorXi sampleCounts; Eigen::VectorXd sampleMasses; //@} - //@{ - //! HDF5 compression settings - unsigned H5compress = 5; - unsigned H5chunk = 1024*1024; // (1 MB) - bool H5shuffle = true; - bool H5szip = false; - //@} + //! Covariance storage instance + std::shared_ptr covarStore; //! Round time key to emulated fixed-point arithmetic double roundTime(double time) @@ -131,19 +124,6 @@ namespace BasisClasses return std::floor(time * multiplier + 0.5) / multiplier; } - //! Write H5 covariance; returns updated group count - virtual unsigned writeCovarH5(HighFive::Group& group, unsigned count, double time); - - //! Write H5 parameters for coefficient covariance - virtual void writeCovarH5Params(HighFive::File& file) {} - - //! Add coefficient covariance data to an HDF5 file - virtual void extendCoefCovariance(const std::string& filename, double time); - - //! Variance file versioning (Version 1.0 has more granularity and - //! poor compression) - inline static const std::string CovarianceFileVersion = "1.1"; - //! Store full coavariance matrix? bool covar = true; @@ -290,51 +270,36 @@ namespace BasisClasses return varray; } - //! Return a vector of tuples of basis functions and the - //! covariance matrix for subsamples of particles - using CoefCovarType = std::tuple; - virtual std::vector> - getCoefCovariance() - { - // Must be overriden; base implementation throws error - throw std::runtime_error("BiorthBasis::getCoefCovariance: Not implemented for this basis"); - } - - //! Get sample counts for the covariance computation - virtual std::tuple - getCovarSamples() - { - return std::make_tuple(sampleCounts, sampleMasses); - } - //! Write coefficient covariance data to an HDF5 file - void writeCoefCovariance(const std::string& compname, - const std::string& runtag, - double time=0.0); - - //! Choose between float and double storage for covariance - void setCovarFloatType(bool flt) { floatType = flt; } - - //! HDF5 compression settings - void setCovarH5Compress(unsigned level, unsigned chunksize, bool shuffle, bool szip=false) + virtual void writeCoefCovariance(const std::string& compname, const std::string& runtag, + double time=0.0) { - H5compress = level; - H5chunk = chunksize; - H5shuffle = shuffle; - H5szip = szip; + // Must be overriden; base implementation throws error + throw std::runtime_error("BiorthBasis::writeCoefCovariance: " + "Not implemented for this basis"); } //! Make covariance after accumulation virtual void makeCoefCovariance(void) {} //! Enable covariance computation with optional sample time - virtual void enableCoefCovariance(bool pcavar, int sampT_in, bool covar_in) + virtual void enableCoefCovariance(bool pcavar, int sampT_in, bool ftype, bool covar_in) { // Must be overriden; base implementation throws error throw std::runtime_error("BiorthBasis::enableCoefCovariance: " "Not implemented for this basis"); } + //! HDF5 compression settings + void setCovarH5Compress(unsigned level, unsigned chunksize, bool shuffle, bool szip=false) + { + if (covarStore) { + covarStore->setCovarH5Compress(level, chunksize, shuffle, szip); + } else { + throw std::runtime_error("BiorthBasis::setCovarH5Compress: covariance storage not initialized"); + } + } + //! Evaluate acceleration in Cartesian coordinates in centered //! coordinate system for a collections of points virtual RowMatrixXd& getAccel(Eigen::Ref pos) @@ -518,19 +483,39 @@ namespace BasisClasses return ret; } - //! Return a vector of tuples of basis functions and the - //! covariance matrix for subsamples of particles - std::vector> - getCoefCovariance(); + //! Write coefficient covariance data to an HDF5 file + virtual void writeCoefCovariance(const std::string& compname, const std::string& runtag, double time=0.0) + { + if (covarStore) { + std::vector> covarData(meanV.size()); + for (int T=0; TwriteCoefCovariance(compname, runtag, elements, time); + } else { + throw std::runtime_error("Spherical::writeCoefCovariance: covariance storage not initialized"); + } + } + + //! Request and initialize covariance computation virtual void enableCoefCovariance(bool pcavar_in, int sampT_in=100, - bool covar_in=true) + bool ftype=false, bool covar_in=true) { pcavar = pcavar_in; sampT = sampT_in; covar = covar_in; - if (pcavar) init_covariance(); + if (pcavar) { + init_covariance(); + covarStore = std::make_shared + ([this](HighFive::File& file){this->writeCovarH5Params(file);}, BasisID, ftype, covar); + } } + //! Create and the coefficients from a function callback with the //! provided time value. Set potential=true if function is a //! potential field. Density is assumed if false (default). @@ -990,7 +975,6 @@ namespace BasisClasses //! For coefficient writing typedef Eigen::Matrix EigenColMajor; - //@{ //! Basis construction parameters @@ -1103,24 +1087,21 @@ namespace BasisClasses return ret; } - //! Return a vector of tuples of basis functions and the - //! covariance matrix for subsamples of particles - std::vector> - getCoefCovariance() - { - return sl->getCoefCovariance(); - } - - std::tuple - getCovarSamples() + //! Write coefficient covariance data to an HDF5 file + virtual void writeCoefCovariance(const std::string& compname, const std::string& runtag, double time=0.0) { - // Copy to base class members - std::tie(sampleCounts, sampleMasses) = sl->getCovarSamples(); - return {sampleCounts, sampleMasses}; + if (covarStore) { + auto covarData = sl->getCoefCovariance(); + std::tie(sampleCounts, sampleMasses) = sl->getCovarSamples(); + SubsampleCovariance::CovarData elements = std::make_tuple(sampleCounts, sampleMasses, covarData); + covarStore->writeCoefCovariance(compname, runtag, elements, time); + } else { + throw std::runtime_error("Spherical::writeCoefCovariance: covariance storage not initialized"); + } } virtual void enableCoefCovariance(bool pcavar_in, int sampT_in=100, - bool covar_in=true) + bool ftype=false, bool covar_in=true) { pcavar = pcavar_in; sampT = sampT_in; @@ -1128,6 +1109,8 @@ namespace BasisClasses if (pcavar) { sl->setSampT(std::max(1, sampT)); sl->set_covar(true); + covarStore = std::make_shared + ([this](HighFive::File& file){this->writeCovarH5Params(file);}, BasisID, ftype, covar); } } //! Create and the coefficients from a function callback with the @@ -1443,25 +1426,35 @@ namespace BasisClasses return true; } - //! Return a vector of tuples of basis functions and the - //! covariance matrix for subsamples of particles - std::vector> getCoefCovariance(); - - //! Get sample count and mass arrays - std::tuple - getCovarSamples() + //! Write coefficient covariance data to an HDF5 file + virtual void writeCoefCovariance(const std::string& compname, const std::string& runtag, double time=0.0) { - return {sampleCounts, sampleMasses}; + if (covarStore) { + std::vector> covarData(meanV.size()); + for (int T=0; TwriteCoefCovariance(compname, runtag, elements, time); + } else { + throw std::runtime_error("Spherical::writeCoefCovariance: covariance storage not initialized"); + } } + //! Enable coefficient covariance computation virtual void enableCoefCovariance(bool pcavar_in, int sampT_in=100, - bool covar_in=false) + bool ftype=false, bool covar_in=false) { pcavar = pcavar_in; sampT = sampT_in; covar = covar_in; - if (pcavar) { init_covariance(); } + if (pcavar) { + init_covariance(); + covarStore = std::make_shared + ([this](HighFive::File& file){this->writeCovarH5Params(file);}, BasisID, ftype, covar); + } } //! Return wave-number indices from flattened index @@ -1471,69 +1464,6 @@ namespace BasisClasses unsigned index1D(int kx, int ky, int kz); }; - class CovarianceReader - { - public: - - using CoefCovarType = std::tuple; - - protected: - - std::vector times; - std::vector sampleCounts; - std::vector sampleMasses; - std::vector>> covarData; - - const double fixedPointPrecision = 1.0e08; - std::map timeMap; - unsigned timeIndex(double time) - { - int itime = static_cast(time * fixedPointPrecision + 0.5); - auto it = timeMap.find(itime); - if (it == timeMap.end()) - throw std::runtime_error("CovarianceReader: time not found"); - return it->second; - } - - std::string basisID; - - public: - - //! Constructor from HDF5 file - CovarianceReader(const std::string& filename, int stride=1); - - //! Get time list - std::vector Times() { return times; } - - //@{ - //! Return a vector of tuples of basis functions and the - //! covariance matrix for subsamples of particles - std::vector> - getCoefCovariance(unsigned index) { return covarData[index]; } - - std::vector> - getCoefCovariance(double time) { return covarData[timeIndex(time)]; } - //@} - - - //@{ - //! Get sample counts for the covariance computation - std::tuple - getCovarSamples(unsigned index) - { - return {sampleCounts[index], sampleMasses[index]}; - } - - std::tuple - getCovarSamples(double time) - { - unsigned index = timeIndex(time); - return {sampleCounts[index], sampleMasses[index]}; - } - //@} - - std::string basisIDname() { return basisID; } - }; //! Time-dependent potential-density model using BasisCoef = std::tuple, std::shared_ptr>; diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index d8574afcb..38ada210d 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1024,29 +1024,6 @@ namespace BasisClasses return ret; } - - /** Return a vector of tuples of basis functions and the covariance - matrix for subsamples of particles */ - std::vector> - Spherical::getCoefCovariance() - { - std::vector> ret; - if (pcavar) { - ret.resize(sampT); - for (int T=0; T(ret[T][l]) = meanV[T][l]; - std::get<1>(ret[T][l]) = covrV[T][l]; - } - } - } - - return ret; - } - - #define MINEPS 1.0e-10 void BiorthBasis::legendre_R(int lmax, double x, Eigen::MatrixXd& p) @@ -3893,7 +3870,8 @@ namespace BasisClasses "check", "method", "pcavar," - "subsamp" + "subsamp", + "nint" }; Cube::Cube(const YAML::Node& CONF) : BiorthBasis(CONF, "cube") @@ -4200,6 +4178,7 @@ namespace BasisClasses /** Return a vector of tuples of basis functions and the covariance matrix for subsamples of particles for Cube type */ + /* std::vector> Cube::getCoefCovariance() { @@ -4220,6 +4199,7 @@ namespace BasisClasses return ret; } + */ std::vector Cube::crt_eval(double x, double y, double z) @@ -5136,564 +5116,6 @@ namespace BasisClasses file.createAttribute("nmaxy", HighFive::DataSpace::From(nmaxy)).write(nmaxy); file.createAttribute("nmaxz", HighFive::DataSpace::From(nmaxz)).write(nmaxz); } - - unsigned BiorthBasis::writeCovarH5(HighFive::Group& snaps, unsigned count, double time) - { - std::ostringstream stim; - stim << std::setw(8) << std::setfill('0') << std::right << count++; - - // Make a new group for this time - // - HighFive::Group stanza = snaps.createGroup(stim.str()); - - // Add a time attribute - // - time = roundTime(time); - stanza.createAttribute("Time", HighFive::DataSpace::From(time)).write(time); - - // Enable compression - // - auto dcpl1 = HighFive::DataSetCreateProps{}; // sample stats - auto dcpl2 = HighFive::DataSetCreateProps{}; // coefficients - auto dcpl3 = HighFive::DataSetCreateProps{}; // covariance - - // Properties for sample stats - if (H5compress) { - unsigned int csz = sampleCounts.size(); - dcpl1.add(HighFive::Chunking({csz, 1})); - if (H5shuffle) dcpl1.add(HighFive::Shuffle()); - dcpl1.add(HighFive::Deflate(H5compress)); - } - - // Add the sample statistics - // - HighFive::DataSet s1data = stanza.createDataSet("sampleCounts", sampleCounts, dcpl1); - HighFive::DataSet s2data = stanza.createDataSet("sampleMasses", sampleMasses, dcpl1); - - // Covariance data - // - auto covdata = getCoefCovariance(); - - // Number of samples - // - unsigned sampleSize = covdata.size(); - unsigned ltot = covdata[0].size(); - unsigned nmax = std::get<0>(covdata[0][0]).rows(); - unsigned diagonalSize = nmax; - - // Save variance or full covariance - if (covar) diagonalSize = nmax*(nmax + 1)/2; - - // Add data dimensions - // - stanza.createAttribute - ("sampleSize", HighFive::DataSpace::From(sampleSize)).write(sampleSize); - - stanza.createAttribute - ("angularSize", HighFive::DataSpace::From(ltot)).write(ltot); - - stanza.createAttribute - ("rankSize", HighFive::DataSpace::From(nmax)).write(nmax); - - int icov = covar ? 1 : 0; - stanza.createAttribute - ("fullCovar", HighFive::DataSpace::From(icov)).write(icov); - - if (H5compress) { - // Szip parameters - const int options_mask = H5_SZIP_NN_OPTION_MASK; - const int pixels_per_block = 8; - - // Properties for coefficients - // - unsigned int csz2 = nmax * ltot * sampleSize; - HighFive::Chunking data_dims2{std::min(csz2, H5chunk), 1}; - - dcpl2.add(data_dims2); - if (H5shuffle) dcpl2.add(HighFive::Shuffle()); - if (H5szip) { - dcpl2.add(HighFive::Szip(options_mask, pixels_per_block)); - } else { - dcpl2.add(HighFive::Deflate(H5compress)); - } - - // Properties for covariance - // - unsigned int csz3 = ltot * diagonalSize * sampleSize; - HighFive::Chunking data_dims3{std::min(csz3, H5chunk), 1}; - - dcpl3.add(data_dims3); - if (H5shuffle) dcpl3.add(HighFive::Shuffle()); - if (H5szip) { - dcpl3.add(HighFive::Szip(options_mask, pixels_per_block)); - } else { - dcpl3.add(HighFive::Deflate(H5compress)); - } - } - - // Pack the coefficient data - // - if (floatType) { - // Create a vector of doubles for the real and imaginary parts - Eigen::VectorXf real_part(nmax*ltot*sampleSize); - Eigen::VectorXf imag_part(nmax*ltot*sampleSize); - - for (size_t T=0, c=0; T(covdata[T][l])(n)); - imag_part(c) = std::imag(std::get<0>(covdata[T][l])(n)); - } - } - } - - // Create two separate, compressed datasets - stanza.createDataSet("coefficients_real", real_part, dcpl2); - stanza.createDataSet("coefficients_imag", imag_part, dcpl2); - - if (std::get<1>(covdata[0][0]).size()) { - - real_part.resize(ltot*diagonalSize*sampleSize); - imag_part.resize(ltot*diagonalSize*sampleSize); - - for (size_t T=0, c=0; T(covdata[T][l])(n1, n2)); - imag_part(c) = std::imag(std::get<1>(covdata[T][l])(n1, n2)); - } - } - // Pack the diagonal only - // - else { - real_part(c ) = std::real(std::get<1>(covdata[T][l])(n1, n1)); - imag_part(c++) = std::imag(std::get<1>(covdata[T][l])(n1, n1)); - } - } - } - } - // Create two separate, compressed datasets - stanza.createDataSet("covariance_real", real_part, dcpl3); - stanza.createDataSet("covariance_imag", imag_part, dcpl3); - } - - } else { - Eigen::VectorXd real_part(ltot*nmax*sampleSize); - Eigen::VectorXd imag_part(ltot*nmax*sampleSize); - - for (size_t T=0, c=0; T(covdata[T][l])(n)); - imag_part(c) = std::imag(std::get<0>(covdata[T][l])(n)); - } - } - } - - // Create two separate, compressed datasets - // - stanza.createDataSet("coefficients_real", real_part, dcpl2); - stanza.createDataSet("coefficients_imag", imag_part, dcpl2); - - if (std::get<1>(covdata[0][0]).size()) { - - real_part.resize(ltot*diagonalSize*sampleSize); - imag_part.resize(ltot*diagonalSize*sampleSize); - - for (size_t T=0, c=0; T(covdata[T][l])(n1, n2)); - imag_part(c) = std::imag(std::get<1>(covdata[T][l])(n1, n2)); - } - } - // Pack the diagonal only - // - else { - real_part(c ) = std::real(std::get<1>(covdata[T][l])(n1, n1)); - imag_part(c++) = std::imag(std::get<1>(covdata[T][l])(n1, n1)); - } - } - } - } - - // Create two separate, compressed datasets - // - stanza.createDataSet("covariance_real", real_part, dcpl3); - stanza.createDataSet("covariance_imag", imag_part, dcpl3); - } - } - // END: sample loop - - return count; - } - - void BiorthBasis::writeCoefCovariance(const std::string& compname, const std::string& runtag, double time) - { - // Check that variance computation is on - // - if (not pcavar) { - std::cout << "BiorthBasis::writeCoefCovariance: covariance computation is disabled. " - << "Set 'pcavar: true' to enable." << std::endl; - return; - } - - // Only root process writes - // - if (myid) return; - - // Check that there is something to write - // - int totalCount = 0; - std::tie(sampleCounts, sampleMasses) = getCovarSamples(); - totalCount += sampleCounts.sum(); - - if (totalCount==0) { - std::cout << "BiorthBasis::writeCoefCovariance: no data" << std::endl; - return; - } - - // Round time - // - time = roundTime(time); - - // The H5 filename - // - std::string fname = "coefcovar." + compname + "." + runtag + ".h5"; - - // Check if file exists? - // - try { - // Open the HDF5 file in read-write mode, creating if it doesn't - // exist - HighFive::File file(fname, - HighFive::File::ReadWrite | - HighFive::File::Create); - - // Check for version string - std::string path = "CovarianceFileVersion"; - - // Check for valid HDF file by attribute - if (file.hasAttribute(path)) { - extendCoefCovariance(fname, time); - return; - } - - // Write the Version string - // - file.createAttribute("CovarianceFileVersion", HighFive::DataSpace::From(CovarianceFileVersion)).write(CovarianceFileVersion); - - // Write the basis identifier string - // - file.createAttribute("BasisID", HighFive::DataSpace::From(BasisID)).write(BasisID); - - // Write the data type size - // - int sz = 8; if (floatType) sz = 4; - file.createAttribute("FloatSize", HighFive::DataSpace::From(sz)).write(sz); - - // Write the specific parameters - // - writeCovarH5Params(file); - - // Group count variable - // - unsigned count = 0; - HighFive::DataSet dataset = file.createDataSet("count", count); - - // Create a new group for coefficient snapshots - // - HighFive::Group group = file.createGroup("snapshots"); - - // Write the coefficients - // - count = writeCovarH5(group, count, time); - - // Update the count - // - dataset.write(count); - - } catch (const HighFive::Exception& err) { - // Handle HighFive specific errors (e.g., file not found) - throw std::runtime_error - (std::string("BiorthBasis::writeCoefCovariance HighFive Error: ") + err.what()); - } catch (const std::exception& err) { - // Handle other general exceptions - throw std::runtime_error - (std::string("BiorthBasis::writeCoefCovariance Error: ") + err.what()); - } - } - - void BiorthBasis::extendCoefCovariance(const std::string& fname, double time) - { - try { - // Open an hdf5 file - // - HighFive::File file(fname, HighFive::File::ReadWrite); - - // Get the dataset - HighFive::DataSet dataset = file.getDataSet("count"); - - unsigned count; - dataset.read(count); - - HighFive::Group group = file.getGroup("snapshots"); - - // Write the coefficients - // - count = writeCovarH5(group, count, time); - - // Update the count - // - dataset.write(count); - - } catch (HighFive::Exception& err) { - throw std::runtime_error - (std::string("BiorthBasis::extendCoefCovariance: HighFive error: ") + err.what()); - } - } - - // Read covariance data - CovarianceReader::CovarianceReader(const std::string& filename, int stride) - { - try { - // Open an existing hdf5 file for reading - // - HighFive::File file(filename, HighFive::File::ReadOnly); - - // Write the Version string - // - std::string version; - file.getAttribute("CovarianceFileVersion").read(version); - // Check for alpha version - if (version == std::string("1.0")) { - throw std::runtime_error("CovarianceReader: this is an early alpha test version. Please remake your files"); - } - // Test for current version - if (version != std::string("1.1")) { - throw std::runtime_error(std::string("CovarianceReader: unsupported file version, ") + version); - } - - // Read the basis identifier string - // - file.getAttribute("BasisID").read(basisID); - - // Get the float size - int sz = 8; - file.getAttribute("FloatSize").read(sz); - if (sz != 4 and sz != 8) { - std::ostringstream sout; - sout << "CovarianceReader: unsupported float size, " << sz; - throw std::runtime_error(sout.str()); - } - - int lmax, nmax, ltot; - - // Current implemented spherical types - const std::set sphereType = {"Spherical", "SphereSL", "Bessel"}; - - // Currently implemented cylindrical types - const std::set cylinderType = {"Cylindrical"}; - - if (sphereType.find(basisID) != sphereType.end()) { - file.getAttribute("lmax").read(lmax); - file.getAttribute("nmax").read(nmax); - ltot = (lmax+1)*(lmax+2)/2; - } else if (cylinderType.find(basisID) != cylinderType.end()) { - file.getAttribute("mmax").read(lmax); - file.getAttribute("nmax").read(nmax); - ltot = lmax + 1; - } else if (basisID == "Cube") { - int nmaxx, nmaxy, nmaxz; - file.getAttribute("nmaxx").read(nmaxx); - file.getAttribute("nmaxy").read(nmaxy); - file.getAttribute("nmaxz").read(nmaxz); - ltot = (2*nmaxx + 1) * (2*nmaxy + 1) * (2*nmaxz + 1); - } else { - throw std::runtime_error(std::string("CovarianceReader: unknown or unimplemented covariance for basis type, ") + basisID); - } - - // Group count variable - // - unsigned count = 0; - file.getDataSet("count").read(count); - - // Open the snapshot group - // - auto snaps = file.getGroup("snapshots"); - - for (unsigned n=0; n(Time * fixedPointPrecision + 0.5); - timeMap[itime] = times.size(); - times.push_back(Time); - - // Get sample properties - // - sampleCounts.push_back(Eigen::VectorXi()); - stanza.getDataSet("sampleCounts").read(sampleCounts.back()); - - sampleMasses.push_back(Eigen::VectorXd()); - stanza.getDataSet("sampleMasses").read(sampleMasses.back()); - - // Get data attributes - // - int nT, lSize, rank, icov=1; - stanza.getAttribute("sampleSize") .read(nT); - stanza.getAttribute("angularSize").read(lSize); - stanza.getAttribute("rankSize") .read(rank); - - if (stanza.hasAttribute("fullCovar")) - stanza.getAttribute("fullCovar").read(icov); - - // Full covariance or variance? - // - bool covar = icov==1 ? true : false; - - // Allocate sample vector for current time - // - covarData.push_back(std::vector>(nT)); - - // Storage - // - Eigen::VectorXcd data0, data1; - - // Get the flattened coefficient array - // - if (sz==4) { - // Get the real and imaginary parts - // - Eigen::VectorXf data_real = - stanza.getDataSet("coefficients_real").read(); - - Eigen::VectorXf data_imag = - stanza.getDataSet("coefficients_imag").read(); - - // Resize the complex array and assign - // - data0.resize(data_real.size()); - data0.real() = data_real.cast(); - data0.imag() = data_imag.cast(); - - // Check for existence of covariance - // - if (stanza.exist("covariance_real")) { - - data_real = - stanza.getDataSet("covariance_real").read(); - - data_imag = - stanza.getDataSet("covariance_imag").read(); - - // Resize the complex array and assign - data1.resize(data_real.size()); - data1.real() = data_real.cast(); - data1.imag() = data_imag.cast(); - } - } else { - // Get the real and imaginary parts - Eigen::VectorXd data_real = - stanza.getDataSet("coefficients_real").read(); - - Eigen::VectorXd data_imag = - stanza.getDataSet("coefficients_imag").read(); - - // Resize the complex array and assign - data0.resize(data_real.size()); - data0.real() = data_real; - data0.imag() = data_imag; - - // Check for existence of covariance - // - if (stanza.exist("covariance_real")) { - - // Get the real and imaginary parts - data_real = - stanza.getDataSet("covariance_real").read(); - - data_imag = - stanza.getDataSet("covariance_imag").read(); - - // Resize the complex array and assign - data1.resize(data_real.size()); - data1.real() = data_real; - data1.imag() = data_imag; - } - } - - // Positions in data stanzas - int sCof = 0, sCov = 0; - - // Loop through all indices and repack - for (int T=0; T elem(lSize); - for (auto & e : elem) { - // Coefficients - std::get<0>(e).resize(rank); - // Covariance matrix - if (data1.size()) std::get<1>(e).resize(rank, rank); - } - - // Pack the coefficient data - int c = 0; - for (size_t l=0; l(elem[l])(n) = data0(sCof + c++); - } - } - sCof += c; - - // Pack the covariance data - c = 0; - for (size_t l=0; l(elem[l])(n1, n2) = data1(sCov + c++); - if (n1 != n2) - std::get<1>(elem[l])(n2, n1) = std::get<1>(elem[l])(n1, n2); - } - else { - if (n1==n2) - std::get<1>(elem[l])(n1, n2) = data1(sCov + c++); - else - std::get<1>(elem[l])(n1, n2) = 0.0; - } - } - } - } - } - sCov += c; - - // Add the data - covarData.back()[T] = std::move(elem); - } - // END: sample loop - } - // END: snapshot loop - } catch (HighFive::Exception& err) { - std::cerr << err.what() << std::endl; - } - } CoefClasses::CoefStrPtr Spherical::makeFromFunction (std::function func, diff --git a/expui/CMakeLists.txt b/expui/CMakeLists.txt index ba44a03eb..e23b72b7d 100644 --- a/expui/CMakeLists.txt +++ b/expui/CMakeLists.txt @@ -36,10 +36,11 @@ endif() # Make the expui shared library # set(expui_SOURCES BasisFactory.cc BiorthBasis.cc FieldBasis.cc - CoefContainer.cc CoefStruct.cc FieldGenerator.cc expMSSA.cc - Coefficients.cc KMeans.cc Centering.cc ParticleIterator.cc - Koopman.cc BiorthBess.cc SvdSignChoice.cc KDdensity.cc - UnitValidator.cc) + Covariance.cc CoefContainer.cc CoefStruct.cc FieldGenerator.cc + expMSSA.cc Coefficients.cc KMeans.cc Centering.cc + ParticleIterator.cc Koopman.cc BiorthBess.cc SvdSignChoice.cc + KDdensity.cc UnitValidator.cc) + add_library(expui ${expui_SOURCES}) set_target_properties(expui PROPERTIES OUTPUT_NAME expui) target_include_directories(expui PUBLIC ${common_INCLUDE}) diff --git a/expui/Coefficients.H b/expui/Coefficients.H index 9effd7c5b..b6f21c6e1 100644 --- a/expui/Coefficients.H +++ b/expui/Coefficients.H @@ -97,6 +97,9 @@ namespace CoefClasses //! Write parameter attributes (needed for derived classes) virtual void WriteH5Params(HighFive::File& file) = 0; + //! Checking file for parameter consistency + virtual bool CheckH5Params(HighFive::File& file) = 0; + //! Write coefficient data in H5 virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count) = 0; @@ -311,6 +314,9 @@ namespace CoefClasses //! Write parameter attributes virtual void WriteH5Params(HighFive::File& file); + //! Check parameter attributes + virtual bool CheckH5Params(HighFive::File& file); + //! Write coefficient data in H5 virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count); @@ -447,6 +453,9 @@ namespace CoefClasses //! Write parameter attributes virtual void WriteH5Params(HighFive::File& file); + //! Check parameter attributes; returns true for success + virtual bool CheckH5Params(HighFive::File& file); + //! Write coefficient data in H5 virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count); @@ -593,6 +602,9 @@ namespace CoefClasses //! Write parameter attributes virtual void WriteH5Params(HighFive::File& file); + //! Check parameter attributes + virtual bool CheckH5Params(HighFive::File& file); + //! Write coefficient data in H5 virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count); @@ -724,6 +736,9 @@ namespace CoefClasses //! Write parameter attributes virtual void WriteH5Params(HighFive::File& file); + //! Check parameter attributes + virtual bool CheckH5Params(HighFive::File& file); + //! Write coefficient data in H5 virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count); @@ -849,6 +864,9 @@ namespace CoefClasses //! Write parameter attributes virtual void WriteH5Params(HighFive::File& file); + //! Check parameter attributes + virtual bool CheckH5Params(HighFive::File& file); + //! Write coefficient data in H5 virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count); @@ -958,6 +976,9 @@ namespace CoefClasses //! Write parameter attributes virtual void WriteH5Params(HighFive::File& file); + //! Check parameter attributes + virtual bool CheckH5Params(HighFive::File& file); + //! Write coefficient data in H5 virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count); @@ -1072,6 +1093,9 @@ namespace CoefClasses //! Write parameter attributes virtual void WriteH5Params(HighFive::File& file); + //! Check parameter attributes + virtual bool CheckH5Params(HighFive::File& file); + //! Write coefficient data in H5 virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count); @@ -1202,6 +1226,9 @@ namespace CoefClasses //! Write parameter attributes virtual void WriteH5Params(HighFive::File& file); + //! Check parameter attributes + virtual bool CheckH5Params(HighFive::File& file); + //! Write coefficient data in H5 virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count); diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc index dc7444dca..45c7d57a1 100644 --- a/expui/Coefficients.cc +++ b/expui/Coefficients.cc @@ -283,6 +283,18 @@ namespace CoefClasses auto in = stanza.getDataSet("coefficients").read(); + // Sanity check on shape + // + if (in.rows() != (Lmax+1)*(Lmax+2)/2 or in.cols() != Nmax) { + if (myid==0) { + std::cout << "---- SphCoefs WARNING: skipping time=" << Time + << " with shape (" << in.rows() << ", " << in.cols() << ")" + << " expected (" << (Lmax+1)*(Lmax+2)/2 + << ", " << Nmax << ")" << std::endl; + } + continue; + } + // If we have a legacy set of coefficients, re-order the // coefficients to match the new HighFive/Eigen ordering // @@ -840,6 +852,50 @@ namespace CoefClasses file.createAttribute("forceID", HighFive::DataSpace::From(forceID)).write(forceID); } + bool SphCoefs::CheckH5Params(HighFive::File& file) + { + double scale0 = coefs.begin()->second->scale, scale1; + std::string forceID0(coefs.begin()->second->id), forceID1; + int Lmax1, Nmax1; + + file.getAttribute("lmax" ).read(Lmax1 ); + file.getAttribute("nmax" ).read(Nmax1 ); + file.getAttribute("scale" ).read(scale1 ); + file.getAttribute("forceID").read(forceID1); + + bool ret = true; + + if (Lmax != Lmax1) { + if (myid==0) + std::cout << "---- SphCoefs::CheckH5Params: Lmax mismatch " << Lmax + << " != " << Lmax1 << std::endl; + ret = false; + } + + if (Nmax != Nmax1) { + if (myid==0) + std::cout << "---- SphCoefs::CheckH5Params: Nmax mismatch " << Nmax + << " != " << Nmax1 << std::endl; + ret = false; + } + + if (std::abs(scale0 - scale1)/scale0 > 1.e-8) { + if (myid==0) + std::cout << "---- SphCoefs::CheckH5Params: scale mismatch " << scale0 + << " != " << scale1 << std::endl; + ret = false; + } + + if (forceID0 != forceID1) { + if (myid==0) + std::cout << "---- SphCoefs::CheckH5Params: forceID mismatch <" << forceID0 + << "> != <" << forceID1 << ">" << std::endl; + ret = false; + } + + return ret; + } + unsigned SphCoefs::WriteH5Times(HighFive::Group& snaps, unsigned count) { for (auto c : coefs) { @@ -1064,6 +1120,17 @@ namespace CoefClasses auto in = stanza.getDataSet("coefficients").read(); + // Sanity check on shape + // + if (in.rows() != Mmax+1 or in.cols() != Nmax) { + if (myid==0) + std::cout << "---- CylCoefs WARNING: skipping time=" << Time + << " with shape (" << in.rows() << ", " << in.cols() << ")" + << " expected (" << Mmax+1 << ", " << Nmax << ")" + << std::endl; + continue; + } + // If we have a legacy set of coefficients, re-order the // coefficients to match the new HighFive/Eigen ordering // @@ -1256,6 +1323,41 @@ namespace CoefClasses file.createAttribute("forceID", HighFive::DataSpace::From(forceID)).write(forceID); } + bool CylCoefs::CheckH5Params(HighFive::File& file) + { + std::string forceID0(coefs.begin()->second->id), forceID1; + int Mmax1, Nmax1; + + file.getAttribute("mmax" ).read(Mmax1 ); + file.getAttribute("nmax" ).read(Nmax1 ); + file.getAttribute("forceID").read(forceID1); + + bool ret = true; + + if (Mmax != Mmax1) { + if (myid==0) + std::cout << "---- CylCoefs::CheckH5Params: Mmax mismatch " << Mmax + << " != " << Mmax1 << std::endl; + ret = false; + } + + if (Nmax != Nmax1) { + if (myid==0) + std::cout << "---- CylCoefs::CheckH5Params: Nmax mismatch " << Nmax + << " != " << Nmax1 << std::endl; + ret = false; + } + + if (forceID0 != forceID1) { + if (myid==0) + std::cout << "---- CylCoefs::CheckH5Params: forceID mismatch <" << forceID0 + << "> != <" << forceID1 << ">" << std::endl; + ret = false; + } + + return ret; + } + unsigned CylCoefs::WriteH5Times(HighFive::Group& snaps, unsigned count) { for (auto c : coefs) { @@ -1591,6 +1693,48 @@ namespace CoefClasses file.createAttribute("forceID", HighFive::DataSpace::From(forceID)).write(forceID); } + bool SlabCoefs::CheckH5Params(HighFive::File& file) + { + std::string forceID0(coefs.begin()->second->id), forceID1; + int NmaxX1, NmaxY1, NmaxZ1; + + file.getAttribute("nmaxx" ).read(NmaxX1 ); + file.getAttribute("nmaxy" ).read(NmaxY1 ); + file.getAttribute("nmaxz" ).read(NmaxZ1 ); + file.getAttribute("forceID").read(forceID1); + + bool ret = true; + + if (NmaxX != NmaxX1) { + if (myid==0) + std::cout << "---- SlabCoefs::CheckH5Params: NmaxX mismatch " << NmaxX + << " != " << NmaxX1 << std::endl; + ret = false; + } + + if (NmaxY != NmaxY1) { + if (myid==0) + std::cout << "---- SlabCoefs::CheckH5Params: NmaxY mismatch " << NmaxY + << " != " << NmaxY1 << std::endl; + ret = false; + } + + if (NmaxZ != NmaxZ1) { + if (myid==0) + std::cout << "---- SlabCoefs::CheckH5Params: NmaxZ mismatch " << NmaxZ + << " != " << NmaxZ1 << std::endl; + ret = false; + } + + if (forceID0 != forceID1) { + if (myid==0) + std::cout << "---- SlabCoefs::CheckH5Params: forceID mismatch <" << forceID0 + << "> != <" << forceID1 << ">" << std::endl; + ret = false; + } + return ret; + } + unsigned SlabCoefs::WriteH5Times(HighFive::Group& snaps, unsigned count) { for (auto c : coefs) { @@ -1954,6 +2098,48 @@ namespace CoefClasses file.createAttribute("forceID", HighFive::DataSpace::From(forceID)).write(forceID); } + bool CubeCoefs::CheckH5Params(HighFive::File& file) + { + std::string forceID0(coefs.begin()->second->id), forceID1; + int NmaxX1, NmaxY1, NmaxZ1; + + file.getAttribute("nmaxx" ).read(NmaxX1 ); + file.getAttribute("nmaxy" ).read(NmaxY1 ); + file.getAttribute("nmaxz" ).read(NmaxZ1 ); + file.getAttribute("forceID").read(forceID1); + + bool ret = true; + + if (NmaxX != NmaxX1) { + if (myid==0) + std::cout << "---- CubeCoefs::CheckH5Params: NmaxX mismatch " << NmaxX + << " != " << NmaxX1 << std::endl; + ret = false; + } + + if (NmaxY != NmaxY1) { + if (myid==0) + std::cout << "---- CubeCoefs::CheckH5Params: NmaxY mismatch " << NmaxY + << " != " << NmaxY1 << std::endl; + ret = false; + } + + if (NmaxZ != NmaxZ1) { + if (myid==0) + std::cout << "---- CubeCoefs::CheckH5Params: NmaxZ mismatch " << NmaxZ + << " != " << NmaxZ1 << std::endl; + ret = false; + } + if (forceID0 != forceID1) { + if (myid==0) + std::cout << "---- CubeCoefs::CheckH5Params: forceID mismatch <" << forceID0 + << "> != <" << forceID1 << ">" << std::endl; + ret = false; + } + + return ret; + } + unsigned CubeCoefs::WriteH5Times(HighFive::Group& snaps, unsigned count) { for (auto c : coefs) { @@ -2308,6 +2494,35 @@ namespace CoefClasses file.createAttribute("rank", HighFive::DataSpace::From(rank)).write(rank); } + bool TrajectoryData::CheckH5Params(HighFive::File& file) + { + int traj1, rank1; + + file.getAttribute("traj").read(traj1); + file.getAttribute("rank").read(rank1); + + int traj0 = coefs.begin()->second->traj; + int rank0 = coefs.begin()->second->rank; + + bool ret = true; + + if (traj0 != traj1) { + if (myid==0) + std::cout << "---- TrajectoryData::CheckH5Params: traj mismatch " << traj0 + << " != " << traj1 << std::endl; + ret = false; + } + + if (rank0 != rank1) { + if (myid==0) + std::cout << "---- TrajectoryData::CheckH5Params: rank mismatch " << rank0 + << " != " << rank1 << std::endl; + ret = false; + } + + return ret; + } + unsigned TrajectoryData::WriteH5Times(HighFive::Group& snaps, unsigned count) { snaps.createDataSet("times", times); @@ -2552,6 +2767,26 @@ namespace CoefClasses file.createAttribute("cols", HighFive::DataSpace::From(cols)).write(cols); } + bool TableData::CheckH5Params(HighFive::File& file) + { + int cols1; + + file.getAttribute("cols").read(cols1 ); + + int cols0 = coefs.begin()->second->cols; + + bool ret = true; + + if (cols0 != cols1) { + if (myid==0) + std::cout << "---- TableData::CheckH5Params: cols mismatch " << cols0 + << " != " << cols1 << std::endl; + ret = false; + } + + return ret; + } + unsigned TableData::WriteH5Times(HighFive::Group& snaps, unsigned count) { snaps.createDataSet("times", times); @@ -2899,7 +3134,15 @@ namespace CoefClasses // ReadH5Units(file); + // Check the specific parameters + // + if (not CheckH5Params(file)) { + throw std::runtime_error("Coefs::ExtendH5Coefs: " + "H5 parameter check failed, aborting extension"); + } + // Get the dataset + // HighFive::DataSet dataset = file.getDataSet("count"); unsigned count; @@ -3181,6 +3424,61 @@ namespace CoefClasses file.createAttribute ("dof", HighFive::DataSpace::From(dof) ).write(dof); } + bool SphFldCoefs::CheckH5Params(HighFive::File& file) + { + double scale0 = coefs.begin()->second->scale, scale1; + int nfld1, lmax1, nmax1, dof1; + + file.getAttribute("nfld" ).read(nfld1 ); + file.getAttribute("lmax" ).read(lmax1 ); + file.getAttribute("nmax" ).read(nmax1 ); + file.getAttribute("scale").read(scale1); + file.getAttribute("dof" ).read(dof1 ); + + int nfld0 = Nfld; + int lmax0 = Lmax; + int nmax0 = Nmax; + int dof0 = dof; + + bool ret = true; + + if (nfld0 != nfld1) { + if (myid==0) + std::cout << "---- SphFldCoefs::CheckH5Params: nfld mismatch " << nfld0 + << " != " << nfld1 << std::endl; + ret = false; + } + + if (lmax0 != lmax1) { + if (myid==0) + std::cout << "---- SphFldCoefs::CheckH5Params: lmax mismatch " << lmax0 + << " != " << lmax1 << std::endl; + ret = false; + } + if (nmax0 != nmax1) { + if (myid==0) + std::cout << "---- SphFldCoefs::CheckH5Params: nmax mismatch " << nmax0 + << " != " << nmax1 << std::endl; + ret = false; + } + + if (scale0 != scale1) { + if (myid==0) + std::cout << "---- SphFldCoefs::CheckH5Params: scale mismatch " << scale0 + << " != " << scale1 << std::endl; + ret = false; + } + + if (dof0 != dof1) { + if (myid==0) + std::cout << "---- SphFldCoefs::CheckH5Params: dof mismatch " << dof0 + << " != " << dof1 << std::endl; + ret = false; + } + + return ret; + } + unsigned SphFldCoefs::WriteH5Times(HighFive::Group& snaps, unsigned count) { for (auto c : coefs) { @@ -3306,6 +3604,62 @@ namespace CoefClasses file.createAttribute ("dof", HighFive::DataSpace::From(dof) ).write(dof); } + bool CylFldCoefs::CheckH5Params(HighFive::File& file) + { + double scale0 = coefs.begin()->second->scale, scale1; + int nfld1, mmax1, nmax1, dof1; + + file.getAttribute("nfld" ).read(nfld1 ); + file.getAttribute("mmax" ).read(mmax1 ); + file.getAttribute("nmax" ).read(nmax1 ); + file.getAttribute("scale").read(scale1); + file.getAttribute("dof" ).read(dof1 ); + + int nfld0 = Nfld; + int mmax0 = Mmax; + int nmax0 = Nmax; + int dof0 = dof; + + bool ret = true; + + if (nfld0 != nfld1) { + if (myid==0) + std::cout << "---- CylFldCoefs::CheckH5Params: nfld mismatch " << nfld0 + << " != " << nfld1 << std::endl; + ret = false; + } + + if (mmax0 != mmax1) { + if (myid==0) + std::cout << "---- CylFldCoefs::CheckH5Params: mmax mismatch " << mmax0 + << " != " << mmax1 << std::endl; + ret = false; + } + + if (nmax0 != nmax1) { + if (myid==0) + std::cout << "---- CylFldCoefs::CheckH5Params: nmax mismatch " << nmax0 + << " != " << nmax1 << std::endl; + ret = false; + } + + if (scale0 != scale1) { + if (myid==0) + std::cout << "---- CylFldCoefs::CheckH5Params: scale mismatch " << scale0 + << " != " << scale1 << std::endl; + ret = false; + } + + if (dof0 != dof1) { + if (myid==0) + std::cout << "---- CylFldCoefs::CheckH5Params: dof mismatch " << dof0 + << " != " << dof1 << std::endl; + ret = false; + } + + return ret; + } + unsigned CylFldCoefs::WriteH5Times(HighFive::Group& snaps, unsigned count) { for (auto c : coefs) { diff --git a/expui/Covariance.cc b/expui/Covariance.cc new file mode 100644 index 000000000..b7910fe24 --- /dev/null +++ b/expui/Covariance.cc @@ -0,0 +1,569 @@ +#include + +#include "Covariance.H" + +namespace BasisClasses +{ + + unsigned SubsampleCovariance::writeCovarH5 + (HighFive::Group& snaps, CovarData& elem, unsigned count, double time) + { + std::ostringstream stim; + stim << std::setw(8) << std::setfill('0') << std::right << count++; + + // Make a new group for this time + // + HighFive::Group stanza = snaps.createGroup(stim.str()); + + // Add a time attribute + // + time = roundTime(time); + stanza.createAttribute("Time", HighFive::DataSpace::From(time)).write(time); + + // Enable compression + // + auto dcpl1 = HighFive::DataSetCreateProps{}; // sample stats + auto dcpl2 = HighFive::DataSetCreateProps{}; // coefficients + auto dcpl3 = HighFive::DataSetCreateProps{}; // covariance + + // Number of samples + // + unsigned sampleSize = std::get<2>(elem).size(); + unsigned ltot = std::get<2>(elem)[0].size(); + unsigned nmax = std::get<0>(std::get<2>(elem)[0][0]).size(); + unsigned diagonalSize = nmax; + + // Properties for sample stats + if (H5compress) { + unsigned int csz = sampleSize; + dcpl1.add(HighFive::Chunking({csz, 1})); + if (H5shuffle) dcpl1.add(HighFive::Shuffle()); + dcpl1.add(HighFive::Deflate(H5compress)); + } + + // Add the sample statistics + // + HighFive::DataSet s1data = stanza.createDataSet("sampleCounts", std::get<0>(elem), dcpl1); + HighFive::DataSet s2data = stanza.createDataSet("sampleMasses", std::get<1>(elem), dcpl1); + + // Save variance or full covariance + if (covar) diagonalSize = nmax*(nmax + 1)/2; + + // Add data dimensions + // + stanza.createAttribute + ("sampleSize", HighFive::DataSpace::From(sampleSize)).write(sampleSize); + + stanza.createAttribute + ("angularSize", HighFive::DataSpace::From(ltot)).write(ltot); + + stanza.createAttribute + ("rankSize", HighFive::DataSpace::From(nmax)).write(nmax); + + int icov = covar ? 1 : 0; + stanza.createAttribute + ("fullCovar", HighFive::DataSpace::From(icov)).write(icov); + + if (H5compress) { + // Szip parameters + const int options_mask = H5_SZIP_NN_OPTION_MASK; + const int pixels_per_block = 8; + + // Properties for coefficients + // + unsigned int csz2 = nmax * ltot * sampleSize; + HighFive::Chunking data_dims2{std::min(csz2, H5chunk), 1}; + + dcpl2.add(data_dims2); + if (H5shuffle) dcpl2.add(HighFive::Shuffle()); + if (H5szip) { + dcpl2.add(HighFive::Szip(options_mask, pixels_per_block)); + } else { + dcpl2.add(HighFive::Deflate(H5compress)); + } + + // Properties for covariance + // + unsigned int csz3 = ltot * diagonalSize * sampleSize; + HighFive::Chunking data_dims3{std::min(csz3, H5chunk), 1}; + + dcpl3.add(data_dims3); + if (H5shuffle) dcpl3.add(HighFive::Shuffle()); + if (H5szip) { + dcpl3.add(HighFive::Szip(options_mask, pixels_per_block)); + } else { + dcpl3.add(HighFive::Deflate(H5compress)); + } + } + + // Pack the coefficient data + // + if (floatType) { + // Create a vector of doubles for the real and imaginary parts + Eigen::VectorXf real_part(nmax*ltot*sampleSize); + Eigen::VectorXf imag_part(nmax*ltot*sampleSize); + + for (size_t T=0, c=0; T(std::get<2>(elem)[T][l])(n)); + imag_part(c) = std::imag(std::get<0>(std::get<2>(elem)[T][l])(n)); + } + } + } + + // Create two separate, compressed datasets + stanza.createDataSet("coefficients_real", real_part, dcpl2); + stanza.createDataSet("coefficients_imag", imag_part, dcpl2); + + if (std::get<1>(std::get<2>(elem)[0][0]).size()) { + + real_part.resize(ltot*diagonalSize*sampleSize); + imag_part.resize(ltot*diagonalSize*sampleSize); + + for (size_t T=0, c=0; T(std::get<2>(elem)[T][l])(n1, n2)); + imag_part(c) = std::imag(std::get<1>(std::get<2>(elem)[T][l])(n1, n2)); + } + } + // Pack the diagonal only + // + else { + real_part(c ) = std::real(std::get<1>(std::get<2>(elem)[T][l])(n1, n1)); + imag_part(c++) = std::imag(std::get<1>(std::get<2>(elem)[T][l])(n1, n1)); + } + } + } + } + // Create two separate, compressed datasets + stanza.createDataSet("covariance_real", real_part, dcpl3); + stanza.createDataSet("covariance_imag", imag_part, dcpl3); + } + + } else { + Eigen::VectorXd real_part(ltot*nmax*sampleSize); + Eigen::VectorXd imag_part(ltot*nmax*sampleSize); + + for (size_t T=0, c=0; T(std::get<2>(elem)[T][l])(n)); + imag_part(c) = std::imag(std::get<0>(std::get<2>(elem)[T][l])(n)); + } + } + } + + // Create two separate, compressed datasets + // + stanza.createDataSet("coefficients_real", real_part, dcpl2); + stanza.createDataSet("coefficients_imag", imag_part, dcpl2); + + if (std::get<1>(std::get<2>(elem)[0][0]).size()) { + + real_part.resize(ltot*diagonalSize*sampleSize); + imag_part.resize(ltot*diagonalSize*sampleSize); + + for (size_t T=0, c=0; T(std::get<2>(elem)[T][l])(n1, n2)); + imag_part(c) = std::imag(std::get<1>(std::get<2>(elem)[T][l])(n1, n2)); + } + } + // Pack the diagonal only + // + else { + real_part(c ) = std::real(std::get<1>(std::get<2>(elem)[T][l])(n1, n1)); + imag_part(c++) = std::imag(std::get<1>(std::get<2>(elem)[T][l])(n1, n1)); + } + } + } + } + + // Create two separate, compressed datasets + // + stanza.createDataSet("covariance_real", real_part, dcpl3); + stanza.createDataSet("covariance_imag", imag_part, dcpl3); + } + } + // END: sample loop + + return count; + } + + void SubsampleCovariance::writeCoefCovariance + (const std::string& compname, const std::string& runtag, CovarData& elem, double time) + { + // Only root process writes + // + if (myid) return; + + // The H5 filename + // + std::string fname = "coefcovar." + compname + "." + runtag + ".h5"; + + writeCoefCovariance(fname, elem, time); + } + + void SubsampleCovariance::writeCoefCovariance + (const std::string& fname, CovarData& elem, double time) + { + // Only root process writes + // + if (myid) return; + + // Check that there is something to write + // + int totalCount = 0; + totalCount += std::get<0>(elem).sum(); + + if (totalCount==0) { + std::cout << "SubsampleCovariance::writeCoefCovariance: no data" << std::endl; + return; + } + + // Round time + // + time = roundTime(time); + + // Check if file exists? + // + try { + // Open the HDF5 file in read-write mode, creating if it doesn't + // exist + HighFive::File file(fname, + HighFive::File::ReadWrite | + HighFive::File::Create); + + // Check for version string + std::string path = "CovarianceFileVersion"; + + // Check for valid HDF file by attribute + if (file.hasAttribute(path)) { + extendCoefCovariance(fname, elem, time); + return; + } + + // Write the Version string + // + file.createAttribute("CovarianceFileVersion", HighFive::DataSpace::From(CovarianceFileVersion)).write(CovarianceFileVersion); + + // Write the basis identifier string + // + file.createAttribute("BasisID", HighFive::DataSpace::From(BasisID)).write(BasisID); + + // Write the data type size + // + int sz = 8; if (floatType) sz = 4; + file.createAttribute("FloatSize", HighFive::DataSpace::From(sz)).write(sz); + + // Write the specific parameters + // + writeCovarH5Params(file); + + // Group count variable + // + unsigned count = 0; + HighFive::DataSet dataset = file.createDataSet("count", count); + + // Create a new group for coefficient snapshots + // + HighFive::Group group = file.createGroup("snapshots"); + + // Write the coefficients + // + count = writeCovarH5(group, elem, count, time); + + // Update the count + // + dataset.write(count); + + } catch (const HighFive::Exception& err) { + // Handle HighFive specific errors (e.g., file not found) + throw std::runtime_error + (std::string("SubsampleCovariance::writeCoefCovariance HighFive Error: ") + err.what()); + } catch (const std::exception& err) { + // Handle other general exceptions + throw std::runtime_error + (std::string("SubsampleCovariance::writeCoefCovariance Error: ") + err.what()); + } + } + + void SubsampleCovariance::extendCoefCovariance + (const std::string& fname, CovarData& elem, double time) + { + try { + // Open an hdf5 file + // + HighFive::File file(fname, HighFive::File::ReadWrite); + + // Get the dataset + HighFive::DataSet dataset = file.getDataSet("count"); + + unsigned count; + dataset.read(count); + + HighFive::Group group = file.getGroup("snapshots"); + + // Write the coefficients + // + count = writeCovarH5(group, elem, count, time); + + // Update the count + // + dataset.write(count); + + } catch (HighFive::Exception& err) { + throw std::runtime_error + (std::string("SubsampleCovariance::extendCoefCovariance: HighFive error: ") + err.what()); + } + } + + // Read covariance data + SubsampleCovariance::SubsampleCovariance + (const std::string& filename, int stride) + { + try { + // Open an existing hdf5 file for reading + // + HighFive::File file(filename, HighFive::File::ReadOnly); + + // Write the Version string + // + std::string version; + file.getAttribute("CovarianceFileVersion").read(version); + // Check for alpha version + if (version == std::string("1.0")) { + throw std::runtime_error("SubsampleCovariance: this is an early alpha test version. Please remake your files"); + } + // Test for current version + if (version != std::string("1.1")) { + throw std::runtime_error(std::string("SubsampleCovariance: unsupported file version, ") + version); + } + + // Read the basis identifier string + // + file.getAttribute("BasisID").read(BasisID); + + // Get the float size + int sz = 8; + file.getAttribute("FloatSize").read(sz); + if (sz != 4 and sz != 8) { + std::ostringstream sout; + sout << "SubsampleCovariance: unsupported float size, " << sz; + throw std::runtime_error(sout.str()); + } + + int lmax, nmax, ltot; + + // Current implemented spherical types + const std::set sphereType = {"Spherical", "SphereSL", "Bessel"}; + + // Currently implemented cylindrical types + const std::set cylinderType = {"Cylindrical"}; + + std::cout << "SubsampleCovariance: reading basis type " << BasisID << std::endl; + + if (sphereType.find(BasisID) != sphereType.end()) { + file.getAttribute("lmax").read(lmax); + file.getAttribute("nmax").read(nmax); + ltot = (lmax+1)*(lmax+2)/2; + } else if (cylinderType.find(BasisID) != cylinderType.end()) { + file.getAttribute("mmax").read(lmax); + file.getAttribute("nmax").read(nmax); + ltot = lmax + 1; + } else if (BasisID == "Cube") { + int nmaxx, nmaxy, nmaxz; + file.getAttribute("nmaxx").read(nmaxx); + file.getAttribute("nmaxy").read(nmaxy); + file.getAttribute("nmaxz").read(nmaxz); + ltot = (2*nmaxx + 1) * (2*nmaxy + 1) * (2*nmaxz + 1); + } else { + throw std::runtime_error(std::string("SubsampleCovariance: unknown or unimplemented covariance for basis type, ") + BasisID); + } + + // Group count variable + // + unsigned count = 0; + file.getDataSet("count").read(count); + + // Open the snapshot group + // + auto snaps = file.getGroup("snapshots"); + + for (unsigned n=0; n>(nT) }; + + // Storage + // + Eigen::VectorXcd data0, data1; + + // Get the flattened coefficient array + // + if (sz==4) { + // Get the real and imaginary parts + // + Eigen::VectorXf data_real = + stanza.getDataSet("coefficients_real").read(); + + Eigen::VectorXf data_imag = + stanza.getDataSet("coefficients_imag").read(); + + // Resize the complex array and assign + // + data0.resize(data_real.size()); + data0.real() = data_real.cast(); + data0.imag() = data_imag.cast(); + + // Check for existence of covariance + // + if (stanza.exist("covariance_real")) { + + data_real = + stanza.getDataSet("covariance_real").read(); + + data_imag = + stanza.getDataSet("covariance_imag").read(); + + // Resize the complex array and assign + data1.resize(data_real.size()); + data1.real() = data_real.cast(); + data1.imag() = data_imag.cast(); + } + } else { + // Get the real and imaginary parts + Eigen::VectorXd data_real = + stanza.getDataSet("coefficients_real").read(); + + Eigen::VectorXd data_imag = + stanza.getDataSet("coefficients_imag").read(); + + // Resize the complex array and assign + data0.resize(data_real.size()); + data0.real() = data_real; + data0.imag() = data_imag; + + // Check for existence of covariance + // + if (stanza.exist("covariance_real")) { + + // Get the real and imaginary parts + data_real = + stanza.getDataSet("covariance_real").read(); + + data_imag = + stanza.getDataSet("covariance_imag").read(); + + // Resize the complex array and assign + data1.resize(data_real.size()); + data1.real() = data_real; + data1.imag() = data_imag; + } + } + + // Positions in data stanzas + int sCof = 0, sCov = 0; + + // Loop through all indices and repack + for (int T=0; T elem(lSize); + for (auto & e : elem) { + // Coefficients + std::get<0>(e).resize(rank); + // Covariance matrix + if (data1.size()) std::get<1>(e).resize(rank, rank); + } + + // Pack the coefficient data + int c = 0; + for (size_t l=0; l(elem[l])(n) = data0(sCof + c++); + } + } + sCof += c; + + // Pack the covariance data + c = 0; + for (size_t l=0; l(elem[l])(n1, n2) = data1(sCov + c++); + if (n1 != n2) + std::get<1>(elem[l])(n2, n1) = std::get<1>(elem[l])(n1, n2); + } + else { + if (n1==n2) + std::get<1>(elem[l])(n1, n2) = data1(sCov + c++); + else + std::get<1>(elem[l])(n1, n2) = 0.0; + } + } + } + } + } + sCov += c; + + // Add the data + std::get<2>(covarData[Time])[T] = std::move(elem); + } + // END: sample loop + } + // END: snapshot loop + } catch (HighFive::Exception& err) { + std::cerr << err.what() << std::endl; + } + } + +} +// END: namespace BasisClasses diff --git a/include/Covariance.H b/include/Covariance.H new file mode 100644 index 000000000..d894288ae --- /dev/null +++ b/include/Covariance.H @@ -0,0 +1,171 @@ +#ifndef _Covariance_H +#define _Covariance_H + +#include +#include +#include +#include + +#include +#include // For 3d rectangular grids +#include + +#include +#include +#include +#include +#include + +#include "localmpi.H" +#include "exputils.H" + +namespace BasisClasses +{ + /** + Class to maintain subsample covariance of basis functions + */ + class SubsampleCovariance + { + public: + + using ParamCallback = std::function; + + using CoefCovarType = std::tuple; + + using CovarData = std::tuple>>; + + protected: + + ParamCallback paramCallback; + + //! Variance file versioning (Version 1.0 has more granularity and + //! poor compression) + inline static const std::string CovarianceFileVersion = "1.1"; + + //! Store full covariance matrix? + bool covar = true; + + //! Use 32-bit float storage for covariance? + bool floatType = false; + + //! Fixed-point multiplier for time mapping. Eight decimal places + //! should be enough here... + const double multiplier = 1.0e+08; + + //! Round time key to emulated fixed-point arithmetic + double roundTime(double time) + { + return std::floor(time * multiplier + 0.5) / multiplier; + } + + //! Round time key to emulated fixed-point arithmetic + int roundTimeInt(double time) + { + return std::floor(time * multiplier + 0.5); + } + + //@{ + //! HDF5 compression settings + unsigned H5compress = 5; + unsigned H5chunk = 1024*1024; // (1 MB) + bool H5shuffle = true; + bool H5szip = false; + //@} + + //! Basis identifier string + std::string BasisID; + + //! Time list + std::vector times; + + //! Imported data members + std::map covarData; + + //! Write H5 covariance; returns updated group count + virtual unsigned writeCovarH5 + (HighFive::Group& group, CovarData& elem, unsigned count, double time); + + //! Write H5 parameters for coefficient covariance + virtual void writeCovarH5Params(HighFive::File& file) + { + if (paramCallback) { + paramCallback(file); + } else { + throw std::runtime_error("SubsampleCovariance::writeCovarH5Params: no paramCallback defined"); + } + } + + //! Add coefficient covariance data to an HDF5 file + virtual void extendCoefCovariance + (const std::string& filename, CovarData& elem, double time); + + public: + + //! Constructor from scratch + SubsampleCovariance(ParamCallback func, const std::string& BasisID, + bool flt, bool cov=false) + : paramCallback(func), BasisID(BasisID), floatType(flt), covar(cov) {} + + //! Constructor from HDF5 file + SubsampleCovariance(const std::string& filename, int stride=1); + + //! Get time list + std::vector Times() { return times; } + + //! Return a vector of tuples of basis functions and the + //! covariance matrix for subsamples of particles + virtual std::vector> + getCoefCovariance(double time) + { + auto it = covarData.find(roundTime(time)); + if (it == covarData.end()) + throw std::runtime_error("SubsampleCovariance::getCoefCovariance: " + "time not found"); + return std::get<2>(it->second); + } + + //! Get sample counts for the covariance computation + virtual std::tuple + getCovarSamples(double time) + { + auto it = covarData.find(roundTime(time)); + if (it == covarData.end()) + throw std::runtime_error("SubsampleCovariance::getCovarSamples: " + "time not found"); + return std::make_tuple(std::get<0>(it->second), std::get<1>(it->second)); + } + + //! Write coefficient covariance data to an HDF5 file, creating a + //! filename based on component name and run tag + void writeCoefCovariance(const std::string& compname, + const std::string& runtag, + CovarData& covarElem, + double time=0.0); + + //! Write coefficient covariance data to an HDF5 file with a + //! supplied path + filename + void writeCoefCovariance(const std::string& filename, + CovarData& covarElem, + double time=0.0); + + //! Choose between float and double storage for covariance + void setCovarFloatType(bool flt) { floatType = flt; } + + //! HDF5 compression settings + void setCovarH5Compress(unsigned level, unsigned chunksize, bool shuffle, bool szip=false) + { + H5compress = level; + H5chunk = chunksize; + H5shuffle = shuffle; + H5szip = szip; + } + + //! Return the basis ID string that generated this covariance sample + const std::string& basisIDname() const { return BasisID; } + }; +} +// END: namespace BasisClasses + +#endif // _Covariance_H + diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H index ebb2af947..386e0b975 100644 --- a/include/EmpCylSL.H +++ b/include/EmpCylSL.H @@ -205,7 +205,7 @@ protected: std::vector numbT; std::vector massT; unsigned sampT, defSampT; - + bool make_covar = false; std::vector vc, vs; diff --git a/include/plummer.H b/include/plummer.H index e1520ef92..4711b2500 100644 --- a/include/plummer.H +++ b/include/plummer.H @@ -6,6 +6,10 @@ class PlummerSphere : public AxiSymModel { + using AxiSymModel::get_mass; + using AxiSymModel::get_density; + using AxiSymModel::get_pot; + private: double rmin; diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index f16390a64..016b63026 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -439,23 +439,6 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE_PURE(void, BiorthBasis, set_coefs, coefs); } - using CCelement = std::tuple; - using CCreturn = std::vector>; - - CCreturn getCoefCovariance(void) override { - PYBIND11_OVERRIDE(CCreturn, BiorthBasis, getCoefCovariance,); - } - - using Selement = std::tuple; - - Selement getCovarSamples() override { - PYBIND11_OVERRIDE(Selement, BiorthBasis, getCovarSamples,); - } - - void enableCoefCovariance(bool pcavar, int nsamples, bool covar) override { - PYBIND11_OVERRIDE(void, BiorthBasis, enableCoefCovariance, pcavar, nsamples, covar); - } - }; class PySpherical : public Spherical @@ -521,8 +504,8 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE_PURE(std::vector, Spherical, orthoCheck, knots); } - void enableCoefCovariance(bool pcavar, int nsamples, bool covar) override { - PYBIND11_OVERRIDE(void, Spherical, enableCoefCovariance, pcavar, nsamples, covar); + void enableCoefCovariance(bool pcavar, int nsamples, bool ftype, bool covar) override { + PYBIND11_OVERRIDE(void, Spherical, enableCoefCovariance, pcavar, nsamples, ftype, covar); } }; @@ -576,8 +559,8 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(void, Cylindrical, make_coefs,); } - void enableCoefCovariance(bool pcavar, int nsamples, bool covar) override { - PYBIND11_OVERRIDE(void, Cylindrical, enableCoefCovariance, pcavar, nsamples, covar); + void enableCoefCovariance(bool pcavar, int nsamples, bool ftype, bool covar) override { + PYBIND11_OVERRIDE(void, Cylindrical, enableCoefCovariance, pcavar, nsamples, ftype, covar); } }; @@ -871,21 +854,8 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(void, Cube, make_coefs,); } - using CCelement = std::tuple; - using CCreturn = std::vector>; - - CCreturn getCoefCovariance(void) override { - PYBIND11_OVERRIDE(CCreturn, Cube, getCoefCovariance,); - } - - using Selement = std::tuple; - - Selement getCovarSamples() override { - PYBIND11_OVERRIDE(Selement, Cube, getCovarSamples,); - } - - void enableCoefCovariance(bool pcavar, int nsamples, bool covar) override { - PYBIND11_OVERRIDE(void, Cube, enableCoefCovariance, pcavar, nsamples, covar); + void enableCoefCovariance(bool pcavar, int nsamples, bool ftype, bool covar) override { + PYBIND11_OVERRIDE(void, Cube, enableCoefCovariance, pcavar, nsamples, ftype, covar); } }; @@ -1362,62 +1332,6 @@ void BasisFactoryClasses(py::module &m) makeFromArray : create coefficients contributions )", py::arg("mass"), py::arg("pos")) - .def("getCoefCovariance", &BasisClasses::BiorthBasis::getCoefCovariance, - R"( - Return the coefficient vectors and covariance matrices for each partition of the - accumulated particles. The number of partitions is set by the configuration - parameter 'sampT' (default: 100). - - The first dimension are the time samples. The second dimension is the angular - index. Each element is a tuple of the coefficient vector and covariance. - The values are complex128 and represents the full amplitude and phase information. - - Returns - ------- - arrays: list[list[tuple[np.ndarray, np.ndarray]]] - Each element is a tuple (coef_vector, coef_covariance_matrix), - where coef_vector and coef_covariance_matrix are numpy.ndarray. - Each coef_covariance_matrix is of shape (nmax, nmax). All values are - complex128. - - Shape and Indexing - ------------------ - - The first list index is the number of time samples. - - The second list index is the angular elements. For spherical bases, all (l, m) pairs - are in triangular index order with l in [0,...,lmax] and m in [0,...,l] for a total of - (lmax+1)*(lmax+2)/2 entries. For cylindrical bases, there are (mmax+1) harmonic - entries for each value m in [0,...,mmax]. - - Each covariance matrix is of shape (nmax, nmax), where nmax is the number of basis - functions - - See also - -------- - getCovarSamples : get counts and mass in each partition - writeCoefCovariance : write the coefficient vectors, covariance - matrices, and metadata to an HDF5 file. - )") - .def("getCovarSamples", &BasisClasses::BiorthBasis::getCovarSamples, - R"( - Return a vector of counts and mass used to compute the partitioned - vectors and covariance. The length of both returned vectors equals 'sampT' - - Parameters - ---------- - None - - Returns - ------- - arrays: tuple(ndarray, ndarray) - a tuple of arrays containing the counts and mass in each - partitioned sample - - See also - -------- - getCoefCovariance : get the coefficient vectors and covariance matrices - for the partitioned phase space. - writeCoefCovariance : write the coefficient vectors, covariance - matrices, and metadata to an HDF5 file. - )") .def("writeCoefCovariance", &BasisClasses::BiorthBasis::writeCoefCovariance, R"( Write the partitioned coefficient vectors and covariance matrices @@ -1449,21 +1363,6 @@ void BasisFactoryClasses(py::module &m) for the partitioned phase space. getCovarSamples : get counts and mass in each partition )", py::arg("compname"), py::arg("runtag"), py::arg("time")=0.0) - .def("setCovarFloatType", &BasisClasses::BiorthBasis::setCovarFloatType, - R"( - Set the floating point type used for covariance storage in HDF5. Set to - true for float4. Type float8 is the default (false). This does not affect - the return type in getCoefCovariance(). - - Parameters - ---------- - floatType : bool - Use 'float32' rather than 'float64' if true - - Returns - ------- - None - )") .def("enableCoefCovariance", &BasisClasses::BiorthBasis::enableCoefCovariance, R"( Enable or disable the coefficient covariance computation and set the @@ -1475,6 +1374,9 @@ void BasisFactoryClasses(py::module &m) enable (true) or disable (false) the covariance computation nsamples : int number of time partitions to use for covariance computation + ftype: bool + if true, use float32 for covariance storage; if false, + use float64 (default: false) covar: bool if true, compute and save covariance to the HDF5 file; if false, save mean and variance vectors only (default: true) @@ -1482,7 +1384,8 @@ void BasisFactoryClasses(py::module &m) Returns ------- None - )", py::arg("pcavar"), py::arg("nsamples")=100, py::arg("covar")=true) + )", py::arg("pcavar"), py::arg("nsamples")=100, py::arg("ftype")=false, + py::arg("covar")=true) .def("setCovarH5Compress", &BasisClasses::BiorthBasis::setCovarH5Compress, R"( Set the HDF5 compression level for covariance storage in HDF5. The Szip @@ -2498,6 +2401,9 @@ void BasisFactoryClasses(py::module &m) enable (true) or disable (false) the covariance computation nsamples : int number of time partitions to use for covariance computation + ftype: bool + if true, use float32 for covariance storage; if false, + use float64 (default: false) covar: bool if true, compute and save covariance to the HDF5 file; if false, save mean and variance vectors only (default: false) @@ -2512,7 +2418,8 @@ void BasisFactoryClasses(py::module &m) because the number of basis functions can be large. To save disk space, covariance computation is disabled by default. The user can enable it by calling this member function with covar=True. - )", py::arg("pcavar"), py::arg("nsamples")=100, py::arg("covar")=false) + )", py::arg("pcavar"), py::arg("nsamples")=100, py::arg("ftype")=false, + py::arg("covar")=false) .def("index1D", &BasisClasses::Cube::index1D, R"( Returns a flattened 1-d index into the arrays and matrices returned by the @@ -2962,7 +2869,7 @@ void BasisFactoryClasses(py::module &m) py::arg("ps"), py::arg("basiscoef"), py::arg("func"), py::arg("nout")=0); - py::class_> + py::class_> (m, "CovarianceReader") .def(py::init(), R"( @@ -2981,7 +2888,7 @@ void BasisFactoryClasses(py::module &m) the new instance )", py::arg("filename"), py::arg("stride")=1) - .def("Times", &BasisClasses::CovarianceReader::Times, + .def("Times", &BasisClasses::SubsampleCovariance::Times, R"( Get the list of evaluation times @@ -2995,23 +2902,7 @@ void BasisFactoryClasses(py::module &m) a list of evaluation times )") .def("getCoefCovariance", static_cast>> - (BasisClasses::CovarianceReader::*)(unsigned)>(&BasisClasses::CovarianceReader::getCoefCovariance), - R"( - Get the covariance matrices for the basis coefficients - - Parameters - ---------- - index : int - the time index - Returns - ------- - list(list(tuple(numpy.ndarray, numpy.ndarray))) - list of partitioned coefficients and their covariance matrices for - each subsample - )", - py::arg("index")) - .def("getCoefCovariance", static_cast>> - (BasisClasses::CovarianceReader::*)(double)>(&BasisClasses::CovarianceReader::getCoefCovariance), + (BasisClasses::SubsampleCovariance::*)(double)>(&BasisClasses::SubsampleCovariance::getCoefCovariance), R"( Get the covariance matrices for the basis coefficients @@ -3027,23 +2918,7 @@ void BasisFactoryClasses(py::module &m) each subsample. The returns are complex-valued. )", py::arg("time")) - .def("getCovarSamples", - py::overload_cast(&BasisClasses::CovarianceReader::getCovarSamples), - R"( - Get sample counts for the covariance computation - - Parameters - ---------- - index : unsigned int - the time index - - Returns - ------- - tuple(numpy.ndarray, numpy.ndarray) - sample counts and masses for the covariance computation - )", - py::arg("index")) - .def("getCovarSamples", py::overload_cast(&BasisClasses::CovarianceReader::getCovarSamples), + .def("getCovarSamples", &BasisClasses::SubsampleCovariance::getCovarSamples, R"( Get sample counts for the covariance computation @@ -3058,7 +2933,7 @@ void BasisFactoryClasses(py::module &m) sample counts and masses for the covariance computation )", py::arg("time")) - .def("basisIDname", &BasisClasses::CovarianceReader::basisIDname, + .def("basisIDname", &BasisClasses::SubsampleCovariance::basisIDname, R"( Get the basis ID name diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc index d81054bb1..1164ebf6b 100644 --- a/pyEXP/CoefWrappers.cc +++ b/pyEXP/CoefWrappers.cc @@ -37,7 +37,8 @@ void CoefficientClasses(py::module &m) { coefficient structure using Python. To do this, use the constructor to make a blank instance, assign the dimensions and use assign() to create a data matrix with the supplied matrix or - array. The dimen- sions are: + array. The dimensions are: + 1. (lmax, nmax) for SphStruct 2. (mmax, nmax) for a CylStruct 3. (nmaxx, nmaxy, nmaxz) for a SlabStruct @@ -143,6 +144,10 @@ void CoefficientClasses(py::module &m) { PYBIND11_OVERRIDE_PURE(void, Coefs, WriteH5Params, file); } + bool CheckH5Params(HighFive::File& file) override { + PYBIND11_OVERRIDE_PURE(bool, Coefs, CheckH5Params, file); + } + unsigned WriteH5Times(HighFive::Group& group, unsigned count) override { PYBIND11_OVERRIDE_PURE(unsigned, Coefs, WriteH5Times, group, count); } @@ -232,6 +237,10 @@ void CoefficientClasses(py::module &m) { PYBIND11_OVERRIDE(void, SphCoefs, WriteH5Params, file); } + bool CheckH5Params(HighFive::File& file) override { + PYBIND11_OVERRIDE(bool, SphCoefs, CheckH5Params, file); + } + unsigned WriteH5Times(HighFive::Group& group, unsigned count) override { PYBIND11_OVERRIDE(unsigned, SphCoefs, WriteH5Times, group, count); } @@ -311,6 +320,10 @@ void CoefficientClasses(py::module &m) { PYBIND11_OVERRIDE(void, CylCoefs, WriteH5Params, file); } + bool CheckH5Params(HighFive::File& file) override { + PYBIND11_OVERRIDE(bool, CylCoefs, CheckH5Params, file); + } + unsigned WriteH5Times(HighFive::Group& group, unsigned count) override { PYBIND11_OVERRIDE(unsigned, CylCoefs, WriteH5Times, group, count); } @@ -389,6 +402,10 @@ void CoefficientClasses(py::module &m) { PYBIND11_OVERRIDE(void, SlabCoefs, WriteH5Params, file); } + bool CheckH5Params(HighFive::File& file) override { + PYBIND11_OVERRIDE(bool, SlabCoefs, CheckH5Params, file); + } + unsigned WriteH5Times(HighFive::Group& group, unsigned count) override { PYBIND11_OVERRIDE(unsigned, SlabCoefs, WriteH5Times, group, count); } @@ -468,6 +485,10 @@ void CoefficientClasses(py::module &m) { PYBIND11_OVERRIDE(void, CubeCoefs, WriteH5Params, file); } + bool CheckH5Params(HighFive::File& file) override { + PYBIND11_OVERRIDE(bool, CubeCoefs, CheckH5Params, file); + } + unsigned WriteH5Times(HighFive::Group& group, unsigned count) override { PYBIND11_OVERRIDE(unsigned, CubeCoefs, WriteH5Times, group, count); } @@ -547,6 +568,10 @@ void CoefficientClasses(py::module &m) { PYBIND11_OVERRIDE(void, TableData, WriteH5Params, file); } + bool CheckH5Params(HighFive::File& file) override { + PYBIND11_OVERRIDE(bool, TableData, CheckH5Params, file); + } + unsigned WriteH5Times(HighFive::Group& group, unsigned count) override { PYBIND11_OVERRIDE(unsigned, TableData, WriteH5Times, group, count); } @@ -625,6 +650,10 @@ void CoefficientClasses(py::module &m) { PYBIND11_OVERRIDE(void, TrajectoryData, WriteH5Params, file); } + bool CheckH5Params(HighFive::File& file) override { + PYBIND11_OVERRIDE(bool, TrajectoryData, CheckH5Params, file); + } + unsigned WriteH5Times(HighFive::Group& group, unsigned count) override { PYBIND11_OVERRIDE(unsigned, TrajectoryData, WriteH5Times, group, count); } diff --git a/src/AxisymmetricBasis.H b/src/AxisymmetricBasis.H index fdf2cbb3d..c4b473453 100644 --- a/src/AxisymmetricBasis.H +++ b/src/AxisymmetricBasis.H @@ -43,6 +43,9 @@ @param tk_type is the smoothing type, one of: Hall, VarianceCut, CumulativeCut, VarianceWeighted @param subsamp true sets partition variance computation (default: false) + + @param nint is the step cadence for subsample variance computation + (default: 0, no subsampling) */ class AxisymmetricBasis : public Basis { @@ -77,6 +80,9 @@ protected: //! First step for PCA computation int npca0; + // Frequency of subsample computation + int nint = 0; + //! Hall smoothing exponent (default: 1.0) double hexp; @@ -125,6 +131,7 @@ protected: //@{ //! Mass and counts for subsample covariance + std::vector countT, countT1; std::vector massT, massT1; unsigned sampT, defSampT; //@} diff --git a/src/AxisymmetricBasis.cc b/src/AxisymmetricBasis.cc index 6409501b0..18fb43000 100644 --- a/src/AxisymmetricBasis.cc +++ b/src/AxisymmetricBasis.cc @@ -17,10 +17,12 @@ AxisymmetricBasis::valid_keys = "dof", "npca", "npca0", + "nint", "pcavar", "pcaeof", "pcadiag", "pcavtk", + "covar", "subsamp", "hexp", "snr", @@ -44,6 +46,7 @@ AxisymmetricBasis:: AxisymmetricBasis(Component* c0, const YAML::Node& conf) : dof = 3; npca = 500; npca0 = 0; + nint = 0; pcavar = false; pcaeof = false; pcadiag = false; @@ -82,6 +85,15 @@ AxisymmetricBasis:: AxisymmetricBasis(Component* c0, const YAML::Node& conf) : if (conf["tk_type"]) tk_type = setTK(conf["tk_type"].as()); if (conf["Mmax"] and not conf["Lmax"]) Lmax = Mmax; + + if (conf["nint"]) { // For OutSample support + nint = conf["nint"].as(); + if (nint>0) { + pcavar = true; + subsamp = true; + } + } + } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameters in AxisymmetricBasis: " @@ -1206,10 +1218,12 @@ void AxisymmetricBasis::parallel_gather_coef2(void) // Report mass used [with storage sanity checks] // - if (massT1.size() == massT.size() and massT.size()==sampT) + if (massT1.size() == massT.size() and massT.size()==sampT) { + MPI_Allreduce(&countT1[0], &countT[0], sampT, + MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD); MPI_Allreduce(&massT1[0], &massT[0], sampT, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - else { + } else { std::cout << "[" << myid << "] AxisymmetricBasis: " << "coef2 out of bounds in mass" << std::endl; } diff --git a/src/Basis.H b/src/Basis.H index 119467da6..51229aa38 100644 --- a/src/Basis.H +++ b/src/Basis.H @@ -3,6 +3,7 @@ #include "PotAccel.H" #include +#include //! Defines a basis-based potential and acceleration class class Basis : public PotAccel @@ -39,7 +40,6 @@ public: double *tdens0, double *dpotl0, double *tdens, double *tpotl, double *tpotr, double *tpotz, double *tpotp) = 0; - /** @name Utility functions */ // @{ diff --git a/src/Basis.cc b/src/Basis.cc index 8eb381a1f..5d94de5a9 100644 --- a/src/Basis.cc +++ b/src/Basis.cc @@ -110,5 +110,3 @@ void Basis::sinecosine_R(int mmax, double phi, } } } - - diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 42682157f..fa7ea2057 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -8,8 +8,8 @@ set(exp_SOURCES Basis.cc Bessel.cc Component.cc Cube.cc Cylinder.cc generateRelaxation.cc HaloBulge.cc incpos.cc incvel.cc ComponentContainer.cc OutAscii.cc OutHDF5.cc OutMulti.cc OutRelaxation.cc OrbTrace.cc OutDiag.cc OutLog.cc OutVel.cc - OutCoef.cc multistep.cc parse.cc SlabSL.cc step.cc tidalField.cc - ultra.cc ultrasphere.cc MPL.cc OutFrac.cc OutCalbr.cc + OutCoef.cc OutSample.cc multistep.cc parse.cc SlabSL.cc step.cc + tidalField.cc ultra.cc ultrasphere.cc MPL.cc OutFrac.cc OutCalbr.cc ParticleFerry.cc chkSlurm.c chkTimer.cc GravKernel.cc CenterFile.cc PolarBasis.cc FlatDisk.cc signals.cc) diff --git a/src/Cube.H b/src/Cube.H index 38c27a4ec..80b3ad5a8 100644 --- a/src/Cube.H +++ b/src/Cube.H @@ -47,6 +47,24 @@ dimensions at once and planes, axes, and 1d which compute wave vectors by plane. This may be extended to specific algorithms in the future. This does nothing for the CPU-only computation. + + @param wrap boolean to enforce periodic wrapping (default: true) + + @param nint integer number of steps between subsample computations + (default: 0, no subsampling) + + @param samplesz integer number of particles in subsample (default: + 100) + + @param subsampleFloat boolean to use single precision in HDF5 + files (default: false) + + @note The subsampling for the cpu implemtation is round-robin + while the subsampling for the gpu implementation is a partition. + This implies that the coefficient values will agree but the + subsamples contain different particles. Since the goal of + subsampling is statistical assessment, this should not be a + limitation. */ class Cube : public PotAccel { @@ -65,15 +83,35 @@ private: //@{ //! Variables int imx, imy, imz, osize; - int use1, use0; + int use1, use0, sampT=0, nint=0; std::complex kfac; - double dfac; + double dfac, last=-std::numeric_limits::max(); bool wrap = true; // Enforce periodic wrapping //@} //! Valid keys for YAML configurations static const std::set valid_keys; + //@{ + //! Covariance structures. First index is T, second is the + //! flattened 3-d k vector + std::vector meanV; + std::vector covrV; + + //@{ + //! Per thread covariance structures + std::vector> meanV1; + std::vector> covrV1; + std::vector workV1; + std::vector countV1; + std::vector massV1; + //@} + + void init_covariance(); + void zero_covariance(); + //@} + + #if HAVE_LIBCUDA==1 virtual void determine_coefficients_cuda(); virtual void determine_acceleration_cuda(); @@ -104,7 +142,12 @@ private: thrust::device_vector> dw_coef; thrust::device_vector> df_coef; - void resize_coefs(int N, int osize, int gridSize, int stride); + thrust::device_vector u_d; + + std::vector>> T_samp; + + void resize_coefs(int N, int osize, int psize, int gridSize, + int stride, int sampT); }; //! A storage instance @@ -249,18 +292,20 @@ private: to = tmp; } + //! Write basis-specific parameters to HDF5 covariance file + void writeCovarH5Params(HighFive::File& file); public: - //! Id string - std::string id; - //! Constructor Cube(Component* c0, const YAML::Node& conf); //! Destructor virtual ~Cube(); + //! Bring in base-class overloads + using PotAccel::determine_coefficients; + //! Compute the coefficients void determine_coefficients(void); @@ -269,6 +314,15 @@ public: //! Coefficient output void dump_coefs_h5(const std::string& file); + + /** Return a subsample of the basis coefficients and covariance + elements for matrices for analysis. The default implementation + returns empty vectors. This parallels the implementation in + expui. */ + CovarData getSubsample(); + + //! Indicate whether or not subsampling is supported + virtual bool hasSubsample() { return true; } }; diff --git a/src/Cube.cc b/src/Cube.cc index 7d55ddd14..994c23112 100644 --- a/src/Cube.cc +++ b/src/Cube.cc @@ -15,7 +15,10 @@ Cube::valid_keys = { "nmaxy", "nmaxz", "method", - "wrap" + "wrap", + "nint", + "samplesz", + "subsampleFloat" }; //@{ @@ -35,6 +38,7 @@ Cube::Cube(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) coef_dump = true; byPlanes = true; cuMethod = "planes"; + sampT = 100; // Default parameter values // @@ -93,6 +97,15 @@ Cube::Cube(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) // dfac = 2.0*M_PI; kfac = std::complex(0.0, dfac); + + // Initialize covariance + // + init_covariance(); + +#if HAVE_LIBCUDA==1 + cuda_initialize(); +#endif + } Cube::~Cube(void) @@ -119,6 +132,16 @@ void Cube::initialize(void) if (conf["nmaxz" ]) nmaxz = conf["nmaxz" ].as(); if (conf["method"]) cuMethod = conf["method"].as(); if (conf["wrap" ]) wrap = conf["wrap" ].as(); + + if (conf["nint"]) { + nint = conf["nint"].as(); + if (nint>0) computeSubsample = true; + } + + if (conf["samplesz"]) { + sampT = conf["samplesz"].as(); + } + } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameters in Cube: " @@ -134,12 +157,81 @@ void Cube::initialize(void) if (myid==0) std::cout << "---- Cube::initialize: wrap=" << std::boolalpha << wrap << std::endl; +} -#if HAVE_LIBCUDA==1 - cuda_initialize(); -#endif +void Cube::init_covariance() +{ + if (computeSubsample) { + + meanV.resize(sampT); + for (auto& v : meanV) { + v.resize(osize); + } + + workV1.resize(nthrds); + for (auto& v : workV1) v.resize(osize); + + if (fullCovar) { + covrV.resize(sampT); + for (auto& v : covrV) { + v.resize(osize, osize); + } + } else { + covrV.clear(); + } + + sampleCounts.resize(sampT); + sampleMasses.resize(sampT); + + meanV1.resize(nthrds); + covrV1.resize(nthrds); + countV1.resize(nthrds); + massV1.resize(nthrds); + + for (int n=0; n facx, facy, facz; int ix, iy, iz; + if (requestSubsample) workV1[id].setZero(); + for (facx=startx, ix=0; ixPart(i); if (part->indx < 10) { @@ -244,9 +343,31 @@ void * Cube::determine_coefficients_thread(void * arg) << std::endl; } } + // END: deepDebug } + // END: iz loop } + // END: iy loop } + // END: ix loop + + if (requestSubsample) { + // Which subsample bin? + // + int T = q % sampT; + + // Accumulate counts and masses + // + countV1[id](T) += 1; + massV1[id](T) += mass; + + // Accumulate subsample contributions + // + meanV1[id][T] += workV1[id] * mass; + if (fullCovar) + covrV1[id][T] += workV1[id] * workV1[id].adjoint() * mass; + } + } // END: particle loop } @@ -287,6 +408,19 @@ void Cube::determine_coefficients(void) for (int n=0; n::max()) { + if (nint>0 && this_step % nint == 0) { + if (tnow > last) { + requestSubsample = true; + last = tnow; + zero_covariance(); + } + } + } else { + subsampleComputed = false; + } + #if HAVE_LIBCUDA==1 (*barrier)("Cube::entering cuda coefficients", __FILE__, __LINE__); if (component->cudaDevice>=0 and use_cuda) { @@ -328,9 +462,51 @@ void Cube::determine_coefficients(void) // Last level? // if (multistep and mlevel==multistep) { - compute_multistep_coefficients(); + compute_multistep_coefficients(); } + // Accumulate mean and covariance subsample contributions + // + if (requestSubsample) { + + // Only finalize at the last multistep level + // + if ( (multistep and mlevel==multistep) or multistep==0 ) { + + // Sum over threads + // + for (int n=1; n("nminx", HighFive::DataSpace::From(nminx)).write(nminx); + file.createAttribute("nminy", HighFive::DataSpace::From(nminy)).write(nminy); + file.createAttribute("nminz", HighFive::DataSpace::From(nminz)).write(nminz); + file.createAttribute("nmaxx", HighFive::DataSpace::From(nmaxx)).write(nmaxx); + file.createAttribute("nmaxy", HighFive::DataSpace::From(nmaxy)).write(nmaxy); + file.createAttribute("nmaxz", HighFive::DataSpace::From(nmaxz)).write(nmaxz); +} + + +PotAccel::CovarData Cube::getSubsample() +{ + using covTuple = std::tuple; + // Prepare the covariance structure + std::vector< std::vector > covar(sampT); + + // Resize the inner vector to 1 (one component for Cube) + for (auto & v : covar) v.resize(1); + + // Fill the covariance structure with subsamples + for (int T=0; T valid_keys; + //! Write basis-specific parameters to HDF5 covariance file + void writeCovarH5Params(HighFive::File& file); + public: //! Mutexes for multithreading @@ -480,6 +487,14 @@ public: //! Sanity check on grid: dumps SM-style images of initial field void dump_mzero(const string& name, int step); + /** Return a subsample of the basis coefficients and covariance + elements for matrices for analysis. The default implementation + returns empty vectors. This parallels the implementation in + expui. */ + CovarData getSubsample(); + + //! Indicate whether or not subsampling is supported + virtual bool hasSubsample() { return true; } }; diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 448a5f769..3c07a73f5 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -62,6 +62,7 @@ Cylinder::valid_keys = { "pcavtk", "pcadiag", "subsamp", + "nint", "try_cache", "density", "EVEN_M", @@ -137,6 +138,7 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : pcadiag = false; pcaeof = false; subsamp = false; + nint = 0; nvtk = 1; pcainit = true; coef_dump = true; @@ -183,7 +185,11 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : if (mlim>=0) ortho->set_mlim(mlim); if (EVEN_M) ortho->setEven(EVEN_M); ortho->setSampT(defSampT); - + if (nint>0) { + EmpCylSL::PCAVAR = true; + ortho->init_pca(); + } + try { if (conf["tk_type"]) ortho->setTK(conf["tk_type"].as()); } @@ -451,6 +457,7 @@ void Cylinder::initialize() if (conf["pcavtk" ]) pcavtk = conf["pcavtk" ].as(); if (conf["pcadiag" ]) pcadiag = conf["pcadiag" ].as(); if (conf["subsamp" ]) subsamp = conf["subsamp" ].as(); + if (conf["nint" ]) nint = conf["nint" ].as(); if (conf["try_cache" ]) try_cache = conf["try_cache" ].as(); if (conf["EVEN_M" ]) EVEN_M = conf["EVEN_M" ].as(); if (conf["cmap" ]) cmapR = conf["cmap" ].as(); @@ -943,6 +950,24 @@ void Cylinder::determine_coefficients_particles(void) pcainit = false; } + // No covarance by default + // + compute = false; + + // Subsample computation flag + // + if (nint) { + // Computed only for mstep==0 and every nint steps + if (mstep==0 and this_step % nint == 0) { + compute = true; + ortho->set_covar(true); + requestSubsample = true; + } else { + ortho->set_covar(false); + subsampleComputed = false; + } + } + if (pcavar or pcaeof) { if (this_step >= npca0) compute = (mstep == 0) && !( (this_step-npca0) % npca); @@ -1041,6 +1066,12 @@ void Cylinder::determine_coefficients_particles(void) if ((pcavar or pcaeof) and mlevel==0) ortho->pca_hall(compute, subsamp); + // If subsample requested and computed, turn off for next time + if (nint and compute and mlevel==multistep) { + requestSubsample = false; + subsampleComputed = true; + } + //========================= // Apply Hall smoothing //========================= @@ -1816,3 +1847,22 @@ void Cylinder::occt_output() } } } + + +PotAccel::CovarData Cylinder::getSubsample() +{ + std::tie(sampleCounts, sampleMasses) = ortho->getCovarSamples(); + auto covarData = ortho->getCoefCovariance(); + return {sampleCounts, sampleMasses, covarData}; +} + +void Cylinder::writeCovarH5Params(HighFive::File& file) +{ + file.createAttribute("mmax", HighFive::DataSpace::From(mmax)).write(mmax); + file.createAttribute("nmax", HighFive::DataSpace::From(nmax)).write(nmax); + file.createAttribute("rcylmin", HighFive::DataSpace::From(rcylmin)).write(rcylmin); + file.createAttribute("rcylmax", HighFive::DataSpace::From(rcylmax)).write(rcylmax); + file.createAttribute("acyl", HighFive::DataSpace::From(acyl)).write(acyl); + file.createAttribute("hcyl", HighFive::DataSpace::From(hcyl)).write(hcyl); +} + diff --git a/src/Direct.H b/src/Direct.H index c4d705147..ec15423f1 100644 --- a/src/Direct.H +++ b/src/Direct.H @@ -48,6 +48,9 @@ class Direct : public PotAccel { + //! Pull in base-class overrides + using PotAccel::determine_coefficients; + private: int soft_indx; diff --git a/src/ExternalForce.H b/src/ExternalForce.H index c76c3d4a8..2c8744853 100644 --- a/src/ExternalForce.H +++ b/src/ExternalForce.H @@ -11,6 +11,9 @@ */ class ExternalForce : public PotAccel { + //! Bring in base-class overloads + using PotAccel::determine_coefficients; + private: void determine_coefficients(void); diff --git a/src/NoForce.H b/src/NoForce.H index 97f4f6682..5226e1003 100644 --- a/src/NoForce.H +++ b/src/NoForce.H @@ -14,6 +14,9 @@ */ class NoForce : public PotAccel { + //! Pull in base-class overrides + using PotAccel::determine_coefficients; + private: void initialize() {} diff --git a/src/Orient.cc b/src/Orient.cc index a2f4b48a0..9f0eb37c6 100644 --- a/src/Orient.cc +++ b/src/Orient.cc @@ -18,21 +18,20 @@ void EL3::debug() const { if (myid==0) { - cerr << left << setfill('-'); - ostringstream ostr; - ostr << "--- EL3 [" << myid << "] "; - cerr << setw(60) << ostr.str() << endl << setfill(' '); - cerr.precision(4); - cerr << setw(12) << T - << setw(12) << M - << setw(12) << E; - cerr << " ["; - for (int i=1; i<=3; i++) cerr << setw(12) << L[i]; - cerr << "] ["; - for (int i=1; i<=3; i++) cerr << setw(12) << R[i]; - cerr << "]" << endl; - cerr << left << setfill('-') - << setw(60) << '-' << endl << setfill(' '); + std::cerr << left << setfill('-'); + std::ostringstream ostr; ostr << "--- EL3 [" << myid << "] "; + std::cerr << setw(60) << ostr.str() << std::endl << std::setfill(' '); + std::cerr.precision(4); + std::cerr << std::setw(12) << T + << std::setw(12) << M + << std::setw(12) << E; + std::cerr << " ["; + for (int i=0; i<3; i++) std::cerr << std::setw(12) << L[i]; + std::cerr << "] ["; + for (int i=0; i<3; i++) std::cerr << std::setw(12) << R[i]; + std::cerr << "]" << std::endl; + std::cerr << std::left << std::setfill('-') + << std::setw(60) << '-' << std::endl << std::setfill(' '); } } @@ -50,9 +49,9 @@ Orient::Orient(int n, int nwant, int naccel, unsigned Oflg, unsigned Cflg, damp = damping; linear = false; - pos = vector(3); - psa = vector(3); - vel = vector(3); + pos = std::vector(3); + psa = std::vector(3); + vel = std::vector(3); center .setZero(); center0.setZero(); diff --git a/src/OutSample.H b/src/OutSample.H new file mode 100644 index 000000000..eb7407641 --- /dev/null +++ b/src/OutSample.H @@ -0,0 +1,61 @@ +#ifndef _OutSample_H +#define _OutSample_H + +#include "Covariance.H" + +/** @file OutSample.H + + @brief Dump coefficients at each interval + + @param filename of the coefficient file + + @param name of the component to dump + + @note This output class requires a different strategy than most of + the others. Rather than maintaining its own step counter, it + checks the component's subsampling status to determine whether + subsample data is available. Then it updates the HDF5 file. The + `mstep` and `last` parameters in the `Run`. The component's + `PotAccel` instance must support subsampling by overriding and + implementing `subsampleReady()` as well as the underlying + `getSubsample()` member to retrieve the data. +*/ +class OutSample : public Output +{ + +private: + + std::string filename; + double prev = -std::numeric_limits::max(); + Component *tcomp; + unsigned level=5; + unsigned chunksize=1024*1024; // 1 MB + bool shuffle=true; + bool szip=false; + bool floatType = false; + + void initialize(void); + + //! Valid keys for YAML configurations + static const std::set valid_keys; + + //! Coefficient covariance storage instance + std::shared_ptr covarStore; + +public: + + //! Constructor + OutSample(const YAML::Node& conf); + + //! Generate the output + /*! + \param nstep is the current time step used to decide whether or not + to dump + \param last should be true on final step to force phase space dump + indepentently of whether or not the frequency criterion is met + */ + void Run(int nstep, int mstep, bool last); + +}; + +#endif diff --git a/src/OutSample.cc b/src/OutSample.cc new file mode 100644 index 000000000..b0b0ff857 --- /dev/null +++ b/src/OutSample.cc @@ -0,0 +1,122 @@ +#include +#include +#include + +#include "expand.H" +#include "Basis.H" +#include "OutSample.H" + +const std::set +OutSample::valid_keys = { + "floatType", + "filename", + "name", + "level", + "chunksize", + "compress", + "szip" +}; + +OutSample::OutSample(const YAML::Node& conf) : Output(conf) +{ + nintsub = std::numeric_limits::max(); + tcomp = NULL; + + initialize(); + + if (!(tcomp->force->hasSubsample())) { + throw std::runtime_error("OutSample: no subsampling for this force"); + } + + covarStore = std::make_shared + ([this](HighFive::File& file){this->tcomp->force->writeCovarH5Params(file);}, + this->tcomp->force->id, floatType, this->tcomp->force->FullCovar()); + + covarStore->setCovarH5Compress(level, chunksize, shuffle, szip); + + if (myid==0) + std::cout << "---- OutSample: writing subsample coefficients for component '" + << tcomp->name << "' to file " << filename << " for force '" + << tcomp->force->id << "'" << std::endl; +} + +void OutSample::initialize() +{ + // Remove matched keys + // + for (auto v : valid_keys) current_keys.erase(v); + + // Assign values from YAML + // + try { + // Search for desired component by name + if (conf["name"]) { + std::string tmp = conf["name"].as(); + for (auto c : comp->components) { + if (!(c->name.compare(tmp))) tcomp = c; + } + } + + if (!tcomp) { + std::string message = "OutSample: no component to trace. Please specify " + "the component name using the 'name' parameter."; + throw std::runtime_error(message); + } + + if (conf["filename"]) { + filename = outdir + conf["filename"].as(); + } else { + filename = outdir + "coefcovar." + tcomp->name + "." + runtag; + } + + if (conf["floatType"]) { + floatType = conf["floatType"].as(); + } + + if (conf["level"]) { + level = conf["level"].as(); + } + + if (conf["chunksize"]) { + chunksize = conf["chunksize"].as(); + } + + if (conf["shuffle"]) { + shuffle = conf["shuffle"].as(); + } + + if (conf["szip"]) { + szip = conf["szip"].as(); + } + + } + catch (YAML::Exception & error) { + if (myid==0) std::cout << "Error parsing parameters in OutSample: " + << error.what() << std::endl + << std::string(60, '-') << std::endl + << "Config node" << std::endl + << std::string(60, '-') << std::endl + << conf << std::endl + << std::string(60, '-') << std::endl; + throw std::runtime_error("OutSample::initialize: error parsing YAML"); + } +} + +void OutSample::Run(int n, int mstep, bool last) +{ + // Subsampling is not available for output + // + if (not tcomp->force->subsampleReady()) return; + + // Check for repeat time + // + if (tnow <= prev) return; + + // Cache the current simulation time + prev = tnow; + + // Write the subsample data + auto elem = tcomp->force->getSubsample(); + covarStore->writeCoefCovariance(filename, elem, tnow); + +} diff --git a/src/OutputContainer.cc b/src/OutputContainer.cc index c26235e94..64673ea3c 100644 --- a/src/OutputContainer.cc +++ b/src/OutputContainer.cc @@ -20,6 +20,7 @@ #include "OutFrac.H" #include "OutCalbr.H" #include "OutMulti.H" +#include "OutSample.H" OutputContainer::OutputContainer() { @@ -112,6 +113,10 @@ void OutputContainer::initialize(void) out.push_back(new OutCalbr(node)); } + else if ( !name.compare("outsamp") ) { + out.push_back(new OutSample(node)); + } + else { string msg("I don't know about the output type: "); msg += name; diff --git a/src/PotAccel.H b/src/PotAccel.H index 39f10fd3e..6f1f56460 100644 --- a/src/PotAccel.H +++ b/src/PotAccel.H @@ -9,10 +9,18 @@ #include #include +#include + +#include +#include +#include +#include +#include #include "Particle.H" #include "StringTok.H" #include "YamlCheck.H" +#include "localmpi.H" #include "config_exp.h" @@ -34,6 +42,16 @@ public: //! For timing typedef std::vector TList; + /** Define the type returned by getSubsample + + The outermost vectors contain subsamples of mass per sample, + counts per sample, and arrays of coefficients and covariance + elements per sample indexed by harmonic + */ + using CovarElement = + std::tuple< Eigen::VectorXd, Eigen::VectorXi, + std::vector>> >; + private: // Threading stuff @@ -93,6 +111,26 @@ protected: //! Current YAML keys to check configuration std::set current_keys; + //! @name Subsampling + //@{ + + //! Reset subsample data + virtual void subsampleReset() { subsampleComputed = false; } + + //! Save covariance subsamples to file + bool fullCovar = false; + + //! Request to compute subsample data + bool requestSubsample = false; + + //! Enable computation of subsample data + bool computeSubsample = false; + + //! Is subsample computed and ready to write? + bool subsampleComputed = false; + + //@} + public: //! For timing data @@ -103,7 +141,7 @@ public: }; //! Id string - string id; + std::string id; //! Possible geometries enum Geometry {sphere, cylinder, cube, slab, table, other}; @@ -181,6 +219,49 @@ public: virtual void dump_coefs(ostream &out) {}; virtual void dump_coefs_h5(const std::string &file) {}; + //@name Subsampling methods + //@{ + //! Set flag to compute subsample data + virtual void subsampleRequest() { requestSubsample = true; } + + //! Indicate whether or not subsamples are available + virtual bool subsampleReady() { return subsampleComputed; } + + //! Indicate whether or not subsampling is supported + virtual bool hasSubsample() { return false; } + + //! Write basis-specific parameters to HDF5 covariance file + virtual void writeCovarH5Params(HighFive::File& file) + { + if (myid==0) + std::cout << "Basis::writeCovarH5Params: " + << "not implemented for basis type <" << id << ">" + << std::endl; + } + //@} + + //@{ + //! Sample counts and masses for covariance computation + Eigen::VectorXi sampleCounts; + Eigen::VectorXd sampleMasses; + //@} + + //@{ + //! Type for coefficient and covariance matrix tuple + + using CoefCovarType = std::tuple; + + using CovarData = std::tuple>>; + + virtual CovarData getSubsample() + { + // Must be overriden; base implementation throws error + throw std::runtime_error("PotAccel::getSubsample: Not implemented for this basis"); + } + //@} + + /** Update the multi time step force algorithm when moving particle i from level cur to level next @@ -218,6 +299,9 @@ public: //! Check whether coeffcients are not available bool NoCoefs() { return play_back and !play_cnew; } + //! Get fullCovar flag + bool FullCovar() { return fullCovar; } + //! For timing the fraction spent in threaded force methods //@{ void print_timings(const string&, TList& tlist); @@ -234,7 +318,6 @@ public: //! Get unaccounted keys std::set unmatched() { return current_keys; } - }; //! Helper class to pass info to threaded member diff --git a/src/PotAccel.cc b/src/PotAccel.cc index 8449b2ef3..c1d6f9d9f 100644 --- a/src/PotAccel.cc +++ b/src/PotAccel.cc @@ -2,6 +2,11 @@ #include #include +#include + +#include +#include + #include "expand.H" #include "PotAccel.H" diff --git a/src/Shells.H b/src/Shells.H index af15e0ad7..f97cab9bb 100644 --- a/src/Shells.H +++ b/src/Shells.H @@ -27,6 +27,9 @@ class Shells : public PotAccel { + //! Pull in base-class overrides + using PotAccel::determine_coefficients; + private: bool self_consistent, firstime_accel; diff --git a/src/SlabSL.H b/src/SlabSL.H index 573f04859..693ded977 100644 --- a/src/SlabSL.H +++ b/src/SlabSL.H @@ -23,17 +23,20 @@ class SlabSL : public PotAccel { -//! Header structure for Sturm-Liouville slab expansion -//! Used for deprecated native coefficient files -struct SlabSLCoefHeader { - double time; - double zmax; - double h; - int type; - int nmaxx, nmaxy, nmaxz; - int jmax; -}; - + //! Pull in base-class overrides + using PotAccel::determine_coefficients; + + //! Header structure for Sturm-Liouville slab expansion + //! Used for deprecated native coefficient files + struct SlabSLCoefHeader { + double time; + double zmax; + double h; + int type; + int nmaxx, nmaxy, nmaxz; + int jmax; + }; + private: std::shared_ptr grid; diff --git a/src/Sphere.H b/src/Sphere.H index b53aaa575..4c369bc78 100644 --- a/src/Sphere.H +++ b/src/Sphere.H @@ -69,6 +69,9 @@ typedef std::shared_ptr SLGridSphPtr; class Sphere : public SphericalBasis { + //! Pull in base-class overrides + using PotAccel::determine_coefficients; + private: SLGridSphPtr ortho; diff --git a/src/SphericalBasis.H b/src/SphericalBasis.H index c4dc2730f..24bcdb63f 100644 --- a/src/SphericalBasis.H +++ b/src/SphericalBasis.H @@ -407,6 +407,9 @@ protected: void biorthogonality_check(); //@} + //! Write basis-specific parameters to HDF5 covariance file + void writeCovarH5Params(HighFive::File& file); + public: //! Use YAML header in coefficient file @@ -544,6 +547,15 @@ public: //! Dump current coefficients into named HDF5 file void dump_coefs_h5(const std::string& file); + + /** Return a subsample of the basis coefficients and covariance + elements for matrices for analysis. The default implementation + returns empty vectors. This parallels the implementation in + expui. */ + CovarData getSubsample(); + + //! Indicate whether or not subsampling is supported + virtual bool hasSubsample() { return true; } }; #endif diff --git a/src/SphericalBasis.cc b/src/SphericalBasis.cc index 406e748ca..5f2ce4319 100644 --- a/src/SphericalBasis.cc +++ b/src/SphericalBasis.cc @@ -45,7 +45,8 @@ SphericalBasis::valid_keys = { "playback", "coefCompute", "coefMaster", - "orthocheck" + "orthocheck", + "subsampleFloat" }; SphericalBasis::SphericalBasis(Component* c0, const YAML::Node& conf, MixtureBasis *m) : @@ -491,7 +492,8 @@ void * SphericalBasis::determine_coefficients_thread(void * arg) if (pcavar) { whch = indx % sampT; pthread_mutex_lock(&cc_lock); - massT1[whch] += mass; + countT1[whch] += 1; + massT1 [whch] += mass; pthread_mutex_unlock(&cc_lock); } } @@ -679,23 +681,43 @@ void SphericalBasis::determine_coefficients_particles(void) // if (!self_consistent && !firstime_coef && !initializing) return; + // No covarance by default + // + compute = false; + + // Subsample computation flag + // + if (nint) { + // Computed only for mstep==0 and every nint steps + if (mstep==0 and this_step % nint == 0) { + compute = true; + requestSubsample = true; + } + else { + requestSubsample = false; + subsampleComputed = false; + + } + } + if (pcavar or pcaeof) { if (this_step >= npca0) compute = (mstep == 0) && !( (this_step-npca0) % npca); - else - compute = false; } - int loffset, moffset, use1; if (compute) { + requestSubsample = true; + if (massT.size() == 0) { // Allocate storage for subsampling if (defSampT) sampT = defSampT; else sampT = floor(sqrt(component->CurTotal())); massT .resize(sampT, 0); massT1 .resize(sampT, 0); + countT .resize(sampT, 0); + countT1 .resize(sampT, 0); expcoefT .resize(sampT); for (auto & t : expcoefT ) { @@ -734,6 +756,7 @@ void SphericalBasis::determine_coefficients_particles(void) if (pcavar) { for (auto & t : expcoefT1) { for (auto & v : t) v->setZero(); } for (auto & t : expcoefM1) { for (auto & v : t) v->setZero(); } + for (auto & v : countT1) v = 0; for (auto & v : massT1) v = 0; } @@ -889,6 +912,9 @@ void SphericalBasis::determine_coefficients_particles(void) for (int i=0; i; + // Prepare the covariance structure + std::vector< std::vector > covar(sampT); + + // Number of real angular terms l in [0, Lmax], m in [0, l] + int totL = (Lmax+1)*(Lmax+2)/2; + + // Sanity check + if (expcoefT[0].size() != totL) { + std::ostringstream ss; + ss << "SphericalBasis::getSubsample() " + << "internal error: unexpected number of absolute angular terms, " + << expcoefT[0].size() << " != " << totL; + throw std::runtime_error(ss.str()); + } + + // Resize the covariance structure + for (auto & v : covar) v.resize(totL); + + // Resize the sample counts and masses + sampleCounts.resize(sampT); + sampleMasses.resize(sampT); + + // Fill the covariance structure with subsamples + for (int T=0; T((*expcoefT[T][k])(n)); + + if (fullCovar) { + // inner n loop + for (int nn=0; nn((*expcoefM[T][k])(n, nn)); + } + } + } + // Assign data to return structure + covar[T][k++] = std::make_tuple(meanV, covrV); + } + // END: m loop + } + // END: l loop + + sampleCounts[T] = countT[T]; + sampleMasses[T] = massT[T]; + } + // END: T loop + + if (myid==0) { + std::cout << "SphericalBasis::getSubsample(): " + << "returning " << sampT << " subsamples, " + << (fullCovar ? "full" : "diagonal") << " covariance" + << ", with counts=" << sampleCounts.size()*sizeof(int) + << ", masses=" << sampleMasses.size()*sizeof(double) + << ", data=" + << covar.size()*covar[0].size()* + (std::get<0>(covar[0][0]).size() + std::get<1>(covar[0][0]).size())*sizeof(std::complex) + << std::endl; + } + + return {sampleCounts, sampleMasses, covar}; +} + + +void SphericalBasis::writeCovarH5Params(HighFive::File& file) +{ + file.createAttribute("lmax", HighFive::DataSpace::From(Lmax)).write(Lmax); + file.createAttribute("nmax", HighFive::DataSpace::From(nmax)).write(nmax); + file.createAttribute("scale", HighFive::DataSpace::From(scale)).write(scale); + file.createAttribute("rmin", HighFive::DataSpace::From(rmin)).write(rmin); + file.createAttribute("rmax", HighFive::DataSpace::From(rmax)).write(rmax); +} + diff --git a/src/cudaCube.cu b/src/cudaCube.cu index e97aadba5..9f3098066 100644 --- a/src/cudaCube.cu +++ b/src/cudaCube.cu @@ -38,29 +38,16 @@ int cubeIndex(int i, int j, int k) i += cubeNumX; j += cubeNumY; k += cubeNumZ; - return k*cubeNX*cubeNY + j*cubeNX + i; -} -/* -// Index function for modulus coefficients -// -__device__ -thrust::tuple TensorIndices(int indx) -{ - int k = indx/(cubeNX*cubeNY); - int j = (indx - k*cubeNX*cubeNY)/cubeNX; - int i = indx - (j + k*cubeNY)*cubeNX; - - return {i, j, k}; + return (i*cubeNY + j)*cubeNZ + k; } -*/ __device__ thrust::tuple cubeWaveNumbers(int indx) { - int k = indx/(cubeNX*cubeNY); - int j = (indx - k*cubeNX*cubeNY)/cubeNX; - int i = indx - (j + k*cubeNY)*cubeNX; + int i = indx/(cubeNY*cubeNZ); + int j = (indx - i*cubeNY*cubeNZ)/cubeNY; + int k = indx - (j + i*cubeNY)*cubeNX; return {i-cubeNumX, j-cubeNumY, k-cubeNumZ}; } @@ -90,10 +77,10 @@ void testConstantsCube() // void Cube::cuda_initialize() { - // Default + // Turn off full covariance // - byPlanes = true; - + fullCovar = false; + // Make method string lower case // std::transform(cuMethod.begin(), cuMethod.end(), cuMethod.begin(), @@ -114,7 +101,8 @@ void Cube::cuda_initialize() if (cuMethod.find("1d") != std::string::npos) byPlanes = true; if (myid==0) - std::cout << "---- Cube::cuda_initialize: byPlanes=" + std::cout << "---- Cube::cuda_initialize: cuMethod=" << cuMethod << std::endl + << "---- Cube::cuda_initialize: byPlanes=" << std::boolalpha << byPlanes << std::endl; } @@ -163,7 +151,7 @@ void Cube::initialize_constants() __global__ void coefKernelCube (dArray P, dArray I, dArray coef, - int stride, PII lohi) + dArray used, int stride, PII lohi) { // Thread ID // @@ -189,6 +177,8 @@ __global__ void coefKernelCube cuFP_t pos[3] = {p.pos[0], p.pos[1], p.pos[2]}; cuFP_t mm = p.mass; + used._v[i] = mm; // Mark particle as used + // Enforce periodicity in the unit cube // if (cubeWrap) { @@ -269,7 +259,7 @@ __global__ void coefKernelCube __global__ void coefKernelCubeX (dArray P, dArray I, dArray coef, - int jj, int kk, int stride, PII lohi) + dArray used, int ii, int jj, int stride, PII lohi) { // Thread ID // @@ -295,6 +285,8 @@ __global__ void coefKernelCubeX cuFP_t pos[3] = {p.pos[0], p.pos[1], p.pos[2]}; cuFP_t mm = p.mass; + used._v[i] = mm; // Mark particle as used + // Enforce periodicity in the unit cube // if (cubeWrap) { @@ -318,10 +310,10 @@ __global__ void coefKernelCubeX #ifdef NORECURSION // Index loop // - for (int s=0; s P, dArray I, CmplxT acc[3] = {0.0, 0.0, 0.0}, pot = 0.0; cuFP_t pos[3] = {p.pos[0], p.pos[1], p.pos[2]}; - cuFP_t mm = p.mass; - int ind[3]; #ifdef NORECURSION // Index loop @@ -519,21 +509,30 @@ public: } }; -void Cube::cudaStorage::resize_coefs(int N, int osize, int gridSize, int stride) +void Cube::cudaStorage::resize_coefs +(int N, int osize, int psize, int gridSize, int stride, int sampT) { // Reserve space for coefficient reduction // - if (dN_coef.capacity() < osize*N) - dN_coef.reserve(osize*N); + if (dN_coef.capacity() < psize*N) + dN_coef.reserve(psize*N); - if (dc_coef.capacity() < osize*gridSize) - dc_coef.reserve(osize*gridSize); + if (dc_coef.capacity() < psize*gridSize) + dc_coef.reserve(psize*gridSize); // Set space for current step // - dN_coef.resize(osize*N); - dc_coef.resize(osize*gridSize); - dw_coef.resize(osize); // This will stay fixed + dN_coef.resize(psize*N); + dc_coef.resize(psize*gridSize); + dw_coef.resize(psize); // This will stay fixed + + u_d.resize(N); + + // Resize covariance storage, if sampT>0 + if (sampT) { + T_samp.resize(sampT); + for (int T=0; Tstream), cuS.df_coef.begin(), cuS.df_coef.end(), 0.0); + + if (requestSubsample) { + + cuS.T_samp.resize(sampT); + + for (int T=0; Tstream), + cuS.T_samp[T].begin(), cuS.T_samp[T].end(), 0.0); + } + } } void Cube::determine_coefficients_cuda() @@ -597,6 +609,10 @@ void Cube::determine_coefficients_cuda() // cuda_zero_coefs(); + // Zero mass accumulation array + // + thrust::fill(cuS.u_d.begin(), cuS.u_d.end(), 0.0); + // Get sorted particle range for mlevel // PII lohi = component->CudaGetLevelRange(mlevel, mlevel), cur; @@ -660,23 +676,30 @@ void Cube::determine_coefficients_cuda() // int sMemSize = BLOCK_SIZE * sizeof(CmplxT); + std::vector::iterator> bm; + if (requestSubsample) { + for (int T=0; Tstream>>> (toKernel(cs->cuda_particles), toKernel(cs->indx1), - toKernel(cuS.dN_coef), jj, kk, stride, cur); + toKernel(cuS.dN_coef), toKernel(cuS.u_d), ii, jj, stride, cur); // Begin the reduction by blocks [perhaps this should use a // stride?] @@ -686,12 +709,12 @@ void Cube::determine_coefficients_cuda() reduceSum <<stream>>> - (toKernel(cuS.dc_coef), toKernel(cuS.dN_coef), imx, N); + (toKernel(cuS.dc_coef), toKernel(cuS.dN_coef), imz, N); // Finish the reduction for this order in parallel // thrust::counting_iterator index_begin(0); - thrust::counting_iterator index_end(gridSize1*imx); + thrust::counting_iterator index_end(gridSize1*imz); // The key_functor indexes the sum reduced series by array index // @@ -707,17 +730,91 @@ void Cube::determine_coefficients_cuda() cuS.dw_coef.begin(), cuS.dw_coef.end(), beg, beg, thrust::plus()); - thrust::advance(beg, imx); + thrust::advance(beg, imz); + + if (requestSubsample) { + + int sN = N/sampT; + int nT = sampT; + + if (sN==0) { // Fail-safe underrun + sN = 1; + nT = N; + } + + for (int T=0; T(mend - mbeg); + massV1 [0][T] += static_cast(thrust::reduce(mbeg, mend)); + } + + // Begin the reduction per grid block + // + /* A reminder to consider implementing strides in reduceSum */ + /* + unsigned int stride1 = s/BLOCK_SIZE/deviceProp.maxGridSize[0] + 1; + unsigned int gridSize1 = s/BLOCK_SIZE/stride1; + + if (s > gridSize1*BLOCK_SIZE*stride1) gridSize1++; + */ + + unsigned int gridSize1 = s/BLOCK_SIZE; + if (s > gridSize1*BLOCK_SIZE) gridSize1++; + + // Sum reduce into gridsize1 blocks taking advantage of + // GPU warp structure + // + + // Mean computation only + // + reduceSumS + <<stream>>> + (toKernel(cuS.dc_coef), toKernel(cuS.dN_coef), imz, N, k, k+s); + + // Finish the reduction for this order in parallel + // + thrust::counting_iterator index_begin(0); + thrust::counting_iterator index_end(gridSize1*imz); + + thrust::reduce_by_key + ( + thrust::cuda::par.on(cs->stream), + thrust::make_transform_iterator(index_begin, key_functor(gridSize1)), + thrust::make_transform_iterator(index_end, key_functor(gridSize1)), + cuS.dc_coef.begin(), thrust::make_discard_iterator(), cuS.dw_coef.begin() + ); + + + thrust::transform(thrust::cuda::par.on(cs->stream), + cuS.dw_coef.begin(), cuS.dw_coef.end(), + bm[T], bm[T], thrust::plus()); + + thrust::advance(bm[T], imz); + } + } } // END: y wave number loop } // END: z wave number loop } else { - + // Adjust cached storage, if necessary // - cuS.resize_coefs(N, osize, gridSize, stride); + cuS.resize_coefs(N, osize, osize, gridSize, stride, sampT); // Compute the coefficient contribution for each order // @@ -725,7 +822,7 @@ void Cube::determine_coefficients_cuda() coefKernelCube<<stream>>> (toKernel(cs->cuda_particles), toKernel(cs->indx1), - toKernel(cuS.dN_coef), stride, cur); + toKernel(cuS.dN_coef), toKernel(cuS.u_d), stride, cur); // Begin the reduction by blocks [perhaps this should use a // stride?] @@ -756,7 +853,74 @@ void Cube::determine_coefficients_cuda() cuS.dw_coef.begin(), cuS.dw_coef.end(), beg, beg, thrust::plus()); - thrust::advance(beg, osize); + if (requestSubsample) { + + int sN = N/sampT; + int nT = sampT; + + if (sN==0) { // Fail-safe underrun + sN = 1; + nT = N; + } + + for (int T=0; T(mend - mbeg); + massV1 [0][T] += static_cast(thrust::reduce(mbeg, mend)); + + // Begin the reduction per grid block + // + /* A reminder to consider implementing strides in reduceSum */ + /* + unsigned int stride1 = s/BLOCK_SIZE/deviceProp.maxGridSize[0] + 1; + unsigned int gridSize1 = s/BLOCK_SIZE/stride1; + + if (s > gridSize1*BLOCK_SIZE*stride1) gridSize1++; + */ + + unsigned int gridSize1 = s/BLOCK_SIZE; + if (s > gridSize1*BLOCK_SIZE) gridSize1++; + + // Sum reduce into gridsize1 blocks taking advantage of + // GPU warp structure + // + + // Mean computation only + // + reduceSumS + <<stream>>> + (toKernel(cuS.dc_coef), toKernel(cuS.dN_coef), osize, N, k, k+s); + + // Finish the reduction for this order in parallel + // + thrust::counting_iterator index_begin(0); + thrust::counting_iterator index_end(gridSize1*osize); + + thrust::reduce_by_key + ( + thrust::cuda::par.on(cs->stream), + thrust::make_transform_iterator(index_begin, key_functor(gridSize1)), + thrust::make_transform_iterator(index_end, key_functor(gridSize1)), + cuS.dc_coef.begin(), thrust::make_discard_iterator(), cuS.dw_coef.begin() + ); + + + thrust::transform(thrust::cuda::par.on(cs->stream), + cuS.dw_coef.begin(), cuS.dw_coef.end(), + bm[T], bm[T], thrust::plus()); + } + } } use1 += N; // Increment particle count @@ -766,6 +930,18 @@ void Cube::determine_coefficients_cuda() // host_coefs = cuS.df_coef; + if (requestSubsample) { + // T loop + // + for (int T=0; T retV = cuS.T_samp[T]; + + for (int n=0; n(retV[n]); + } + } + } + // Create a wavenumber tuple from a flattened index // auto indices = [&](int indx) @@ -1162,13 +1338,13 @@ void Cube::multistep_update_cuda() if (byPlanes) { // Adjust cached storage, if necessary // - cuS.resize_coefs(N, imx, gridSize, stride); + cuS.resize_coefs(N, osize, imz, gridSize, stride, sampT); // Compute the coefficient contribution for each order // auto beg = cuS.df_coef.begin(); - for (int kk=-nmaxz; kk<=nmaxz; kk++) { + for (int ii=-nmaxx; ii<=nmaxx; ii++) { for (int jj=-nmaxy; jj<=nmaxy; jj++) { @@ -1176,19 +1352,19 @@ void Cube::multistep_update_cuda() // coefKernelCubeX<<stream>>> (toKernel(cs->cuda_particles), toKernel(cs->indx1), - toKernel(cuS.dN_coef), jj, kk, stride, cur); + toKernel(cuS.dN_coef), toKernel(cuS.u_d), ii, jj, stride, cur); unsigned int gridSize1 = N/BLOCK_SIZE; if (N > gridSize1*BLOCK_SIZE) gridSize1++; reduceSum <<stream>>> - (toKernel(cuS.dc_coef), toKernel(cuS.dN_coef), imx, N); + (toKernel(cuS.dc_coef), toKernel(cuS.dN_coef), imz, N); // Finish the reduction for this order in parallel // thrust::counting_iterator index_begin(0); - thrust::counting_iterator index_end(gridSize1*imx); + thrust::counting_iterator index_end(gridSize1*imz); // The key_functor indexes the sum reduced series by array index // @@ -1204,17 +1380,17 @@ void Cube::multistep_update_cuda() cuS.dw_coef.begin(), cuS.dw_coef.end(), beg, beg, thrust::plus()); - thrust::advance(beg, imx); + thrust::advance(beg, imz); } - // END: z wave numbers + // END: y wave numbers } - // END: y wave numbers + // END: x wave numbers } // END: bunch by planes else { // Adjust cached storage, if necessary // - cuS.resize_coefs(N, osize, gridSize, stride); + cuS.resize_coefs(N, osize, osize, gridSize, stride, sampT); // Shared memory size for the reduction // @@ -1228,7 +1404,7 @@ void Cube::multistep_update_cuda() // coefKernelCube<<stream>>> (toKernel(cs->cuda_particles), toKernel(cs->indx1), - toKernel(cuS.dN_coef), stride, cur); + toKernel(cuS.dN_coef), toKernel(cuS.u_d), stride, cur); unsigned int gridSize1 = N/BLOCK_SIZE; if (N > gridSize1*BLOCK_SIZE) gridSize1++;