From 4de87bc16240b186c2117f05941febaa7332bc69 Mon Sep 17 00:00:00 2001 From: joeseaer <2969111832@qq.com> Date: Sun, 5 Jul 2026 14:30:56 +0800 Subject: [PATCH 1/4] ENH:Add 3D PDF calculator support --- .../srreal/ObjCrystStructureAdapter.cpp | 2 +- src/diffpy/srreal/PDF3DCalculator.cpp | 1308 +++++++++++++++++ src/diffpy/srreal/PDF3DCalculator.hpp | 257 ++++ src/diffpy/srreal/PDFCalculator.cpp | 8 + src/diffpy/srreal/PDFCalculator.hpp | 2 + src/diffpy/srreal/R3linalg.cpp | 84 ++ src/diffpy/srreal/R3linalg.hpp | 2 + 7 files changed, 1662 insertions(+), 1 deletion(-) create mode 100644 src/diffpy/srreal/PDF3DCalculator.cpp create mode 100644 src/diffpy/srreal/PDF3DCalculator.hpp diff --git a/src/diffpy/srreal/ObjCrystStructureAdapter.cpp b/src/diffpy/srreal/ObjCrystStructureAdapter.cpp index 2b5e5660..969f5715 100644 --- a/src/diffpy/srreal/ObjCrystStructureAdapter.cpp +++ b/src/diffpy/srreal/ObjCrystStructureAdapter.cpp @@ -86,7 +86,7 @@ fetchSymmetryOperations(const ObjCryst::SpaceGroup& spacegroup) assert(nbtran * last <= nbsym); for (int nt = 0; nt < nbtran; ++nt) { - const REAL* pt = sgtrans[nt].tr; + const auto* pt = sgtrans[nt].tr; R3::Vector sgt(pt[0], pt[1], pt[2]); for (int i = 0; i < last; ++i) { diff --git a/src/diffpy/srreal/PDF3DCalculator.cpp b/src/diffpy/srreal/PDF3DCalculator.cpp new file mode 100644 index 00000000..e839ab35 --- /dev/null +++ b/src/diffpy/srreal/PDF3DCalculator.cpp @@ -0,0 +1,1308 @@ +/************************************************************************************* +* +* libdiffpy by DANSE Diffraction group +* Simon J. L. Billinge +* (c) 2009 The Trustees of Columbia University +* in the City of New York. All rights reserved. +* +* File coded by: Hai Qiao +* +* See AUTHORS.txt for a list of people who contributed. +* See LICENSE_DANSE.txt for license information. +* +/************************************************************************************* +* +* class PDF3DCalculator -- Implementation of PDF3DCalculator +* +*************************************************************************************/ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +using namespace std; + +namespace diffpy { +namespace srreal { + +// Static constants +constexpr double PDF3DCalculator::DEFAULT_RMAX_3D; +constexpr double PDF3DCalculator::DEFAULT_GRID_STEP; + +// Constructor --------------------------------------------------------------- + +PDF3DCalculator::PDF3DCalculator() : + mdr(DEFAULT_GRID_STEP), + mnbins(0), + maccumblocksize(32), + mapplyrho0background3d(true), + musecqwindow3d(true), + mcalculationmode3d(0), + mhistogramweightmode3d(0), + menablennDelta3d(false), + mnnDelta3d(0.0), + mnnDelta3dUpperBound(0.0), + mnnDeltaPositiveEta3d(0.9), + mdistanceDelta1_3d(0.0), + mdistanceDelta2_3d(0.0), + mdeltaPairA3d("*"), + mdeltaPairB3d("*"), + mdeltaShellIndex3d(1), + mdeltaShellTolerance3d(0.05), + mdeltaKeyTolerance3d(1.0e-6), + museadpscaleSensitivity3d(false), + madpScale3d(1.0), + mrho0backgroundscale3d(1.0), + mtotaloccupancy(0.0), + msfaverage(0.0) +{ + this->registerDoubleAttribute("enable_nn_delta3d", + this, + &PDF3DCalculator::getEnableNNDelta3DAttr, + &PDF3DCalculator::setEnableNNDelta3DAttr); + this->registerDoubleAttribute("nn_delta3d", + this, + &PDF3DCalculator::getNNDelta3D, + &PDF3DCalculator::setNNDelta3D); + this->registerDoubleAttribute("nn_delta3d_upper_bound", + this, + &PDF3DCalculator::getNNDelta3DUpperBound); + this->registerDoubleAttribute("nn_delta_positive_eta3d", + this, + &PDF3DCalculator::getNNDeltaPositiveEta3D, + &PDF3DCalculator::setNNDeltaPositiveEta3D); + this->registerDoubleAttribute("delta1_3d", + this, + &PDF3DCalculator::getDistanceDelta1_3D, + &PDF3DCalculator::setDistanceDelta1_3D); + this->registerDoubleAttribute("delta2_3d", + this, + &PDF3DCalculator::getDistanceDelta2_3D, + &PDF3DCalculator::setDistanceDelta2_3D); + this->registerDoubleAttribute("delta_shell_index3d", + this, + &PDF3DCalculator::getDeltaShellIndex3DAttr, + &PDF3DCalculator::setDeltaShellIndex3DAttr); + this->registerDoubleAttribute("delta_shell_tolerance3d", + this, + &PDF3DCalculator::getDeltaShellTolerance3D, + &PDF3DCalculator::setDeltaShellTolerance3D); + this->registerDoubleAttribute("delta_key_tolerance3d", + this, + &PDF3DCalculator::getDeltaKeyTolerance3D, + &PDF3DCalculator::setDeltaKeyTolerance3D); + this->registerDoubleAttribute("use_adp_scale_sensitivity3d", + this, + &PDF3DCalculator::getUseADPScaleSensitivity3DAttr, + &PDF3DCalculator::setUseADPScaleSensitivity3DAttr); + this->registerDoubleAttribute("adp_scale3d", + this, + &PDF3DCalculator::getADPScale3D, + &PDF3DCalculator::setADPScale3D); + this->registerDoubleAttribute("rho0_background_scale3d", + this, + &PDF3DCalculator::getRho0BackgroundScale3D, + &PDF3DCalculator::setRho0BackgroundScale3D); + this->registerDoubleAttribute("calculation_mode3d", + this, + &PDF3DCalculator::getCalculationMode3DAttr, + &PDF3DCalculator::setCalculationMode3DAttr); + this->registerDoubleAttribute("histogram_weight_mode3d", + this, + &PDF3DCalculator::getHistogramWeightMode3DAttr, + &PDF3DCalculator::setHistogramWeightMode3DAttr); + + this->setRmax(DEFAULT_RMAX_3D); + this->setRstep(mdr); + this->setQmax(12.0); +} + +// Public Methods ------------------------------------------------------------ + +void PDF3DCalculator::setGridStep(double dr) +{ + if (dr <= 0) throw std::invalid_argument("Grid step must be positive."); + if (dr != mdr) + { + mdr = dr; + this->setRstep(dr); + this->resetValue(); // triggers re-allocation + } +} + +double PDF3DCalculator::getGridStep() const +{ + return mdr; +} + +void PDF3DCalculator::setAccumBlockSize(int bs) +{ + if (bs <= 0) throw std::invalid_argument("Accumulation block size must be positive."); + maccumblocksize = bs; +} + +int PDF3DCalculator::getAccumBlockSize() const +{ + return maccumblocksize; +} + +void PDF3DCalculator::setApplyRho0Background3D(bool v) +{ + mapplyrho0background3d = v; +} + +bool PDF3DCalculator::getApplyRho0Background3D() const +{ + return mapplyrho0background3d; +} + +void PDF3DCalculator::setUseCQWindow3D(bool v) +{ + musecqwindow3d = v; +} + +bool PDF3DCalculator::getUseCQWindow3D() const +{ + return musecqwindow3d; +} + +void PDF3DCalculator::setCalculationMode3D(int mode) +{ + if (mode != 0 && mode != 1) + throw std::invalid_argument("3D calculation mode must be 0 (ADP PDF) or 1 (vector histogram)."); + if (mcalculationmode3d != mode) mticker.click(); + mcalculationmode3d = mode; +} + +int PDF3DCalculator::getCalculationMode3D() const +{ + return mcalculationmode3d; +} + +double PDF3DCalculator::getCalculationMode3DAttr() const +{ + return static_cast(mcalculationmode3d); +} + +void PDF3DCalculator::setCalculationMode3DAttr(double v) +{ + const int mode = static_cast(std::lround(v)); + if (std::fabs(v - static_cast(mode)) > 1.0e-8) + throw std::invalid_argument("3D calculation mode must be an integer value."); + this->setCalculationMode3D(mode); +} + +void PDF3DCalculator::setHistogramWeightMode3D(int mode) +{ + if (mode != 0 && mode != 1) + throw std::invalid_argument("3D histogram weight mode must be 0 (scattering) or 1 (count)."); + if (mhistogramweightmode3d != mode) mticker.click(); + mhistogramweightmode3d = mode; +} + +int PDF3DCalculator::getHistogramWeightMode3D() const +{ + return mhistogramweightmode3d; +} + +double PDF3DCalculator::getHistogramWeightMode3DAttr() const +{ + return static_cast(mhistogramweightmode3d); +} + +void PDF3DCalculator::setHistogramWeightMode3DAttr(double v) +{ + const int mode = static_cast(std::lround(v)); + if (std::fabs(v - static_cast(mode)) > 1.0e-8) + throw std::invalid_argument("3D histogram weight mode must be an integer value."); + this->setHistogramWeightMode3D(mode); +} + +void PDF3DCalculator::setEnableNNDelta3D(bool v) +{ + if (menablennDelta3d != v) mticker.click(); + menablennDelta3d = v; +} + +bool PDF3DCalculator::getEnableNNDelta3D() const +{ + return menablennDelta3d; +} + +double PDF3DCalculator::getEnableNNDelta3DAttr() const +{ + return menablennDelta3d ? 1.0 : 0.0; +} + +void PDF3DCalculator::setEnableNNDelta3DAttr(double v) +{ + this->setEnableNNDelta3D(v != 0.0); +} + +void PDF3DCalculator::setNNDelta3D(double v) +{ + if (v < 0.0) throw std::invalid_argument("3D NN delta must be non-negative."); + const double bound = this->currentNNDeltaUpperBound(); + if (!mdeltaShellRecords.empty()) + { + if (bound <= 0.0 && v > 0.0) + throw std::invalid_argument("3D NN delta cannot be positive without a positive-definite shell bound."); + if (bound > 0.0 && v >= bound) + throw std::invalid_argument("3D NN delta exceeds the positive-definite shell bound."); + } + if (mnnDelta3d != v) mticker.click(); + mnnDelta3d = v; +} + +const double& PDF3DCalculator::getNNDelta3D() const +{ + return mnnDelta3d; +} + +double PDF3DCalculator::getNNDelta3DUpperBound() const +{ + return this->currentNNDeltaUpperBound(); +} + +const double& PDF3DCalculator::getNNDeltaPositiveEta3D() const +{ + return mnnDeltaPositiveEta3d; +} + +void PDF3DCalculator::setNNDeltaPositiveEta3D(double v) +{ + if (v <= 0.0 || v >= 1.0) + throw std::invalid_argument("Positive-definite eta must be between 0 and 1."); + const double minproj = (mnnDeltaPositiveEta3d > 0.0) ? + (mnnDelta3dUpperBound / mnnDeltaPositiveEta3d) : 0.0; + if (mnnDeltaPositiveEta3d != v) mticker.click(); + mnnDeltaPositiveEta3d = v; + if (minproj > 0.0) mnnDelta3dUpperBound = mnnDeltaPositiveEta3d * minproj; + this->validateNNDelta3D(); + this->validateDistanceDecayDelta3D(); +} + +void PDF3DCalculator::setDistanceDelta1_3D(double v) +{ + if (v < 0.0) throw std::invalid_argument("3D distance-delta delta1 must be non-negative."); + if (mdistanceDelta1_3d != v) mticker.click(); + mdistanceDelta1_3d = v; + this->validateDistanceDecayDelta3D(); +} + +const double& PDF3DCalculator::getDistanceDelta1_3D() const +{ + return mdistanceDelta1_3d; +} + +void PDF3DCalculator::setDistanceDelta2_3D(double v) +{ + if (v < 0.0) throw std::invalid_argument("3D distance-delta delta2 must be non-negative."); + if (mdistanceDelta2_3d != v) mticker.click(); + mdistanceDelta2_3d = v; + this->validateDistanceDecayDelta3D(); +} + +const double& PDF3DCalculator::getDistanceDelta2_3D() const +{ + return mdistanceDelta2_3d; +} + +void PDF3DCalculator::setDeltaPairTypes3D(const std::string& a, const std::string& b) +{ + if (a.empty() || b.empty()) + throw std::invalid_argument("Delta pair atom types must be non-empty."); + if (mdeltaPairA3d != a || mdeltaPairB3d != b) mticker.click(); + mdeltaPairA3d = a; + mdeltaPairB3d = b; +} + +const std::string& PDF3DCalculator::getDeltaPairA3D() const +{ + return mdeltaPairA3d; +} + +const std::string& PDF3DCalculator::getDeltaPairB3D() const +{ + return mdeltaPairB3d; +} + +void PDF3DCalculator::setDeltaShellIndex3D(int v) +{ + if (v <= 0) throw std::invalid_argument("Delta shell index must be positive."); + if (mdeltaShellIndex3d != v) mticker.click(); + mdeltaShellIndex3d = v; +} + +int PDF3DCalculator::getDeltaShellIndex3D() const +{ + return mdeltaShellIndex3d; +} + +double PDF3DCalculator::getDeltaShellIndex3DAttr() const +{ + return static_cast(mdeltaShellIndex3d); +} + +void PDF3DCalculator::setDeltaShellIndex3DAttr(double v) +{ + this->setDeltaShellIndex3D(static_cast(std::floor(v + 0.5))); +} + +const double& PDF3DCalculator::getDeltaShellTolerance3D() const +{ + return mdeltaShellTolerance3d; +} + +void PDF3DCalculator::setDeltaShellTolerance3D(double v) +{ + if (v <= 0.0) throw std::invalid_argument("Delta shell tolerance must be positive."); + if (mdeltaShellTolerance3d != v) mticker.click(); + mdeltaShellTolerance3d = v; +} + +const double& PDF3DCalculator::getDeltaKeyTolerance3D() const +{ + return mdeltaKeyTolerance3d; +} + +void PDF3DCalculator::setDeltaKeyTolerance3D(double v) +{ + if (v <= 0.0) throw std::invalid_argument("Delta key tolerance must be positive."); + if (mdeltaKeyTolerance3d != v) mticker.click(); + mdeltaKeyTolerance3d = v; +} + +int PDF3DCalculator::getDeltaEligiblePairCount3D() const +{ + return static_cast(mdeltaEligibleBondKeys.size()); +} + +void PDF3DCalculator::setUseADPScaleSensitivity3D(bool v) +{ + if (museadpscaleSensitivity3d != v) mticker.click(); + museadpscaleSensitivity3d = v; + this->validateNNDelta3D(); +} + +bool PDF3DCalculator::getUseADPScaleSensitivity3D() const +{ + return museadpscaleSensitivity3d; +} + +double PDF3DCalculator::getUseADPScaleSensitivity3DAttr() const +{ + return museadpscaleSensitivity3d ? 1.0 : 0.0; +} + +void PDF3DCalculator::setUseADPScaleSensitivity3DAttr(double v) +{ + this->setUseADPScaleSensitivity3D(v != 0.0); +} + +void PDF3DCalculator::setADPScale3D(double v) +{ + if (v <= 0.0) throw std::invalid_argument("3D ADP scale must be positive."); + if (madpScale3d != v) mticker.click(); + madpScale3d = v; + this->validateNNDelta3D(); +} + +const double& PDF3DCalculator::getADPScale3D() const +{ + return madpScale3d; +} + +void PDF3DCalculator::setRho0BackgroundScale3D(double v) +{ + if (v < 0.0) throw std::invalid_argument("3D rho0 background scale must be non-negative."); + if (mrho0backgroundscale3d != v) mticker.click(); + mrho0backgroundscale3d = v; +} + +const double& PDF3DCalculator::getRho0BackgroundScale3D() const +{ + return mrho0backgroundscale3d; +} + +QuantityType PDF3DCalculator::get3DPDF() const +{ + QuantityType result; + std::vector grid = mgrid3d; + if (!this->isVectorHistogramMode3D()) + { + if (musecqwindow3d) applyQWindow3D(grid); + + const double rdf_scale = (mtotaloccupancy * msfaverage == 0.0) ? 0.0 : + 1.0 / (mtotaloccupancy * msfaverage * msfaverage); + if (rdf_scale != 1.0) + { + for (double& val : grid) val *= rdf_scale; + } + + if (mapplyrho0background3d) + { + const double rho0_bg = computeRho0Background(); + if (rho0_bg != 0.0) + { + for (double& val : grid) val -= rho0_bg; + } + } + + double qdamp = 0.0; + try + { + qdamp = this->getEnvelopeByType("qresolution")->getDoubleAttr("qdamp"); + } + catch (...) + { + qdamp = 0.0; + } + if (qdamp > 0.0) + { + for (size_t i = 0; i < grid.size(); ++i) + { + if (grid[i] == 0.0) continue; + double x, y, z; + indexToCoord(i, x, y, z); + const double r = sqrt(x * x + y * y + z * z); + grid[i] *= exp(-0.5 * (r * qdamp) * (r * qdamp)); + } + } + } + + size_t nonzero_count = 0; + for (double val : grid) { + if (val != 0.0) ++nonzero_count; + } + result.reserve(nonzero_count * 4); + + for (size_t i = 0; i < grid.size(); ++i) + { + if (grid[i] != 0.0) + { + double x, y, z; + indexToCoord(i, x, y, z); + + result.push_back(x); + result.push_back(y); + result.push_back(z); + result.push_back(grid[i]); + } + } + return result; +} + +QuantityType PDF3DCalculator::getRadialHistogram3D() const +{ + QuantityType result; + result.reserve(mradialhistogram3d.size() * 2); + for (size_t i = 0; i < mradialhistogram3d.size(); ++i) + { + result.push_back(static_cast(i) * mdr); + result.push_back(mradialhistogram3d[i]); + } + return result; +} + +void PDF3DCalculator::exportGrid3DBinary(const std::string& path, bool usefloat32, bool applypost) const +{ + std::vector grid = mgrid3d; + if (applypost) + { + if (musecqwindow3d) applyQWindow3D(grid); + + const double rdf_scale = (mtotaloccupancy * msfaverage == 0.0) ? 0.0 : + 1.0 / (mtotaloccupancy * msfaverage * msfaverage); + if (rdf_scale != 1.0) + for (double& val : grid) val *= rdf_scale; + + if (mapplyrho0background3d) + { + const double rho0_bg = computeRho0Background(); + if (rho0_bg != 0.0) + for (double& val : grid) val -= rho0_bg; + } + + double qdamp = 0.0; + try { qdamp = this->getEnvelopeByType("qresolution")->getDoubleAttr("qdamp"); } + catch (...) { qdamp = 0.0; } + if (qdamp > 0.0) + { + for (size_t i = 0; i < grid.size(); ++i) + { + if (grid[i] == 0.0) continue; + double x, y, z; + indexToCoord(i, x, y, z); + const double r = sqrt(x * x + y * y + z * z); + grid[i] *= exp(-0.5 * (r * qdamp) * (r * qdamp)); + } + } + } + + std::ofstream ofs(path.c_str(), std::ios::binary | std::ios::trunc); + if (!ofs) throw std::runtime_error("Failed to open output file: " + path); + + const size_t nxy = static_cast(mnbins) * mnbins; + if (usefloat32) + { + std::vector row(nxy); + for (int iz = 0; iz < mnbins; ++iz) + { + const size_t base = static_cast(iz) * nxy; + for (size_t i = 0; i < nxy; ++i) row[i] = static_cast(grid[base + i]); + ofs.write(reinterpret_cast(row.data()), static_cast(nxy * sizeof(float))); + } + } + else + { + for (int iz = 0; iz < mnbins; ++iz) + { + const size_t base = static_cast(iz) * nxy; + ofs.write(reinterpret_cast(&grid[base]), static_cast(nxy * sizeof(double))); + } + } + if (!ofs) throw std::runtime_error("Failed while writing output file: " + path); +} + +double PDF3DCalculator::computeRho0Background() const +{ + const StructureAdapterPtr& structure = this->getStructure(); + if (!structure) return 0.0; + const double partialpdfscale = this->getPartialPDFScale(); + return mrho0backgroundscale3d * partialpdfscale * structure->numberDensity(); +} + +void PDF3DCalculator::applyQWindow3D(std::vector& grid) const +{ + if (grid.empty()) return; + + const double qmin = this->getQmin(); + const double qmax = this->getQmax(); + if (qmin <= 0.0 && qmax <= 0.0) return; + + const int n = mnbins; + int npad = 1; + while (npad < n) npad <<= 1; + + const size_t n3 = static_cast(npad) * npad * npad; + std::vector data(2 * n3, 0.0); + + const int center = n / 2; + for (int iz = 0; iz < n; ++iz) + { + const int sz = (iz - center + npad) % npad; + for (int iy = 0; iy < n; ++iy) + { + const int sy = (iy - center + npad) % npad; + const size_t base_src = static_cast(iz * n + iy) * n; + const size_t base_dst = static_cast(sz * npad + sy) * npad; + for (int ix = 0; ix < n; ++ix) + { + const int sx = (ix - center + npad) % npad; + const size_t src = base_src + ix; + const size_t dst = base_dst + sx; + data[2 * dst] = grid[src]; + } + } + } + + for (int iz = 0; iz < npad; ++iz) + { + for (int iy = 0; iy < npad; ++iy) + { + double* row = &data[2 * ((iz * npad + iy) * npad)]; + gsl_fft_complex_radix2_forward(row, 1, npad); + } + } + + for (int iz = 0; iz < npad; ++iz) + { + for (int ix = 0; ix < npad; ++ix) + { + double* col = &data[2 * (iz * npad * npad + ix)]; + gsl_fft_complex_radix2_forward(col, npad, npad); + } + } + + for (int iy = 0; iy < npad; ++iy) + { + for (int ix = 0; ix < npad; ++ix) + { + double* line = &data[2 * (iy * npad + ix)]; + gsl_fft_complex_radix2_forward(line, npad * npad, npad); + } + } + + const double qstep = (npad > 0) ? (2.0 * M_PI / (npad * mdr)) : 0.0; + for (int iz = 0; iz < npad; ++iz) + { + const int kz = (iz <= npad / 2) ? iz : iz - npad; + const double qz = kz * qstep; + for (int iy = 0; iy < npad; ++iy) + { + const int ky = (iy <= npad / 2) ? iy : iy - npad; + const double qy = ky * qstep; + for (int ix = 0; ix < npad; ++ix) + { + const int kx = (ix <= npad / 2) ? ix : ix - npad; + const double qx = kx * qstep; + const double q = sqrt(qx * qx + qy * qy + qz * qz); + if ((qmax > 0.0 && q > qmax) || (qmin > 0.0 && q < qmin)) + { + const size_t idx = (static_cast(iz) * npad + iy) * npad + ix; + data[2 * idx] = 0.0; + data[2 * idx + 1] = 0.0; + } + } + } + } + + for (int iy = 0; iy < npad; ++iy) + { + for (int ix = 0; ix < npad; ++ix) + { + double* line = &data[2 * (iy * npad + ix)]; + gsl_fft_complex_radix2_inverse(line, npad * npad, npad); + } + } + + for (int iz = 0; iz < npad; ++iz) + { + for (int ix = 0; ix < npad; ++ix) + { + double* col = &data[2 * (iz * npad * npad + ix)]; + gsl_fft_complex_radix2_inverse(col, npad, npad); + } + } + + for (int iz = 0; iz < npad; ++iz) + { + for (int iy = 0; iy < npad; ++iy) + { + double* row = &data[2 * ((iz * npad + iy) * npad)]; + gsl_fft_complex_radix2_inverse(row, 1, npad); + } + } + + for (int iz = 0; iz < n; ++iz) + { + const int sz = (iz - center + npad) % npad; + for (int iy = 0; iy < n; ++iy) + { + const int sy = (iy - center + npad) % npad; + const size_t base_src = static_cast(sz * npad + sy) * npad; + const size_t base_dst = static_cast(iz * n + iy) * n; + for (int ix = 0; ix < n; ++ix) + { + const int sx = (ix - center + npad) % npad; + const size_t src = base_src + sx; + const size_t dst = base_dst + ix; + grid[dst] = data[2 * src]; + } + } + } +} + +// Protected Methods --------------------------------------------------------- + +void PDF3DCalculator::resetValue() +{ + const StructureAdapterPtr& structure = this->getStructure(); + if (!structure) + { + msfCache.clear(); + mtotaloccupancy = 0.0; + msfaverage = 0.0; + } + else + { + int nsite = this->countSites(); + msfCache.resize(nsite); + const auto& sftable = this->getScatteringFactorTable(); + + mtotaloccupancy = structure->totalOccupancy(); + double totsf = 0.0; + for (int i = 0; i < nsite; ++i) + { + std::string atom_type = structure->siteAtomType(i); + const double sf = sftable->lookup(atom_type); + const double occ = structure->siteOccupancy(i); + const double mult = structure->siteMultiplicity(i); + msfCache[i] = sf * occ; + totsf += msfCache[i] * mult; + } + msfaverage = (mtotaloccupancy == 0.0) ? 0.0 : (totsf / mtotaloccupancy); + } + + // Ensure odd number of bins so origin is centered. + mnbins = static_cast(2 * ceil(this->getRmax() / mdr)) + 1; + size_t total = static_cast(mnbins) * mnbins * mnbins; + mgrid3d.assign(total, 0.0); + const size_t nr = static_cast(ceil(this->getRmax() / mdr)) + 1; + mradialhistogram3d.assign(nr, 0.0); + + PDFCalculator::resetValue(); + if (!this->isVectorHistogramMode3D()) + { + this->buildDeltaEligibleShellTable(); + this->validateNNDelta3D(); + this->validateDistanceDecayDelta3D(); + } + if (mevaluator) mevaluator->setFlag(USEFULLSUM, true); +} + +void PDF3DCalculator::addPairContribution(const BaseBondGenerator& bnds, int summationscale) +{ + if (bnds.distance() == 0.0) return; + + int i0 = bnds.site0(); + int i1 = bnds.site1(); + + if (i0 >= static_cast(msfCache.size()) || i1 >= static_cast(msfCache.size())) + return; + + const R3::Vector& rvec = bnds.r01(); + const double pairscale = bnds.multiplicity() * static_cast(summationscale); + double sfprod = msfCache[i0] * msfCache[i1] * pairscale; + + if (this->isVectorHistogramMode3D()) + { + const double weight = (mhistogramweightmode3d == 1) ? pairscale : sfprod; + addVectorHistogramToGrid(rvec, bnds.distance(), weight); + return; + } + + const R3::Matrix& U_i = bnds.Ucartesian0(); + const R3::Matrix& U_j = bnds.Ucartesian1(); + R3::Matrix Sigma = this->effectivePairCovariance(bnds, U_i, U_j); + + addAnisotropicGaussianToGrid(rvec, Sigma, sfprod); +} + +void PDF3DCalculator::addVectorHistogramToGrid(const R3::Vector& r_ij, double distance, double weight) +{ + if (weight == 0.0) return; + + const size_t idx = this->coordToIndex(r_ij[0], r_ij[1], r_ij[2]); + if (idx < mgrid3d.size()) + { +#ifdef _OPENMP +#pragma omp atomic +#endif + mgrid3d[idx] += weight; + } + + if (distance >= 0.0 && !mradialhistogram3d.empty()) + { + const int ir = static_cast(std::lround(distance / mdr)); + if (ir >= 0 && ir < static_cast(mradialhistogram3d.size())) + { +#ifdef _OPENMP +#pragma omp atomic +#endif + mradialhistogram3d[static_cast(ir)] += weight; + } + } +} + +void PDF3DCalculator::addAnisotropicGaussianToGrid(const R3::Vector& r_ij, const R3::Matrix& Sigma, double sfprod) +{ + // Eigenvalue decomposition + R3::Vector eigenvalues; + R3::Matrix eigenvectors; + + R3::eigen_solve_3x3(Sigma, eigenvalues, eigenvectors); + + // Check for positive definiteness (eigenvalues are sorted ascending) + if (eigenvalues[0] <= 1e-8) { + return; + } + + // 3. Invert the covariance matrix + R3::Matrix Sigma_inv = R3::inverse(Sigma); + double det_Sigma = R3::determinant(Sigma); + + // 4. Normalization factor + const double two_pi = 2.0 * M_PI; + double norm_factor = sfprod / (pow(two_pi, 1.5) * sqrt(det_Sigma)); + + // 5. Determine sampling bounding box + double max_sigma = 0.0; + for(int k=0; k<3; ++k) { + max_sigma = std::max(max_sigma, sqrt(eigenvalues[k])); + } + + // 4.0 sigma cutoff + double cutoff_radius = 4.0 * max_sigma; + + // Safety clamps + if (cutoff_radius < mdr) cutoff_radius = mdr; + if (cutoff_radius > 10.0) cutoff_radius = 10.0; + + // 6. Iterate over the local grid indices + double halfspan = (mnbins / 2) * mdr; + + auto get_index_range = [&](double center_val) -> std::pair { + double min_val = center_val - cutoff_radius; + double max_val = center_val + cutoff_radius; + + // Node-centered grid: points at k * mdr + int start = static_cast(ceil((min_val + halfspan) / mdr)); + int end = static_cast(floor((max_val + halfspan) / mdr)); + + start = std::max(0, start); + end = std::min(mnbins - 1, end); + + return {start, end}; + }; + + std::pair xr = get_index_range(r_ij[0]); + std::pair yr = get_index_range(r_ij[1]); + std::pair zr = get_index_range(r_ij[2]); + + double cutoff_sq = 16.0; + + const int bs = std::max(1, maccumblocksize); + const int nbz = (zr.second - zr.first + bs) / bs; + const int nby = (yr.second - yr.first + bs) / bs; + const int nbx = (xr.second - xr.first + bs) / bs; + +#ifdef _OPENMP +#pragma omp parallel +#endif + { + std::vector local; + +#ifdef _OPENMP +#pragma omp for collapse(3) schedule(dynamic, 1) +#endif + for (int tbz = 0; tbz < nbz; ++tbz) + { + for (int tby = 0; tby < nby; ++tby) + { + for (int tbx = 0; tbx < nbx; ++tbx) + { + const int bz = zr.first + tbz * bs; + const int by = yr.first + tby * bs; + const int bx = xr.first + tbx * bs; + + const int izhi = std::min(zr.second, bz + bs - 1); + const int iyhi = std::min(yr.second, by + bs - 1); + const int ixhi = std::min(xr.second, bx + bs - 1); + + const int tz = izhi - bz + 1; + const int ty = iyhi - by + 1; + const int tx = ixhi - bx + 1; + const size_t nloc = static_cast(tz) * ty * tx; + local.assign(nloc, 0.0); + + bool any = false; + for (int iz = bz; iz <= izhi; ++iz) + { + const double z_grid = iz * mdr - halfspan; + const double dz = z_grid - r_ij[2]; + for (int iy = by; iy <= iyhi; ++iy) + { + const double y_grid = iy * mdr - halfspan; + const double dy = y_grid - r_ij[1]; + for (int ix = bx; ix <= ixhi; ++ix) + { + const double x_grid = ix * mdr - halfspan; + const double dx = x_grid - r_ij[0]; + R3::Vector delta(dx, dy, dz); + R3::Vector tmp = R3::mxvecproduct(Sigma_inv, delta); + const double mahalanobis_sq = R3::dot(delta, tmp); + if (mahalanobis_sq > cutoff_sq) continue; + + const double val = norm_factor * exp(-0.5 * mahalanobis_sq); + const int lz = iz - bz; + const int ly = iy - by; + const int lx = ix - bx; + const size_t lidx = (static_cast(lz) * ty + ly) * tx + lx; + local[lidx] += val; + any = true; + } + } + } + + if (!any) continue; + +#ifdef _OPENMP +#pragma omp critical(pdf3d_tile_reduce) +#endif + { + for (int lz = 0; lz < tz; ++lz) + { + const int iz = bz + lz; + for (int ly = 0; ly < ty; ++ly) + { + const int iy = by + ly; + for (int lx = 0; lx < tx; ++lx) + { + const size_t lidx = (static_cast(lz) * ty + ly) * tx + lx; + const double v = local[lidx]; + if (v == 0.0) continue; + const int ix = bx + lx; + const size_t gidx = (static_cast(iz) * mnbins + iy) * mnbins + ix; + mgrid3d[gidx] += v; + } + } + } + } + } + } + } + } +} + +// Private Helpers ----------------------------------------------------------- + +bool PDF3DCalculator::DeltaBondKey::operator==(const DeltaBondKey& other) const +{ + return site0 == other.site0 && site1 == other.site1 && + rx == other.rx && ry == other.ry && rz == other.rz; +} + +size_t PDF3DCalculator::DeltaBondKeyHash::operator()(const DeltaBondKey& key) const +{ + size_t seed = 0; + seed ^= std::hash()(key.site0) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + seed ^= std::hash()(key.site1) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + seed ^= std::hash()(key.rx) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + seed ^= std::hash()(key.ry) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + seed ^= std::hash()(key.rz) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + return seed; +} + +bool PDF3DCalculator::matchesDeltaPairTypes( + const std::string& atom0, const std::string& atom1) const +{ + const bool allpairs = (mdeltaPairA3d == "*" && mdeltaPairB3d == "*"); + if (allpairs) return true; + const bool forward = + (mdeltaPairA3d == "*" || atom0 == mdeltaPairA3d) && + (mdeltaPairB3d == "*" || atom1 == mdeltaPairB3d); + const bool reverse = + (mdeltaPairA3d == "*" || atom1 == mdeltaPairA3d) && + (mdeltaPairB3d == "*" || atom0 == mdeltaPairB3d); + if (forward || reverse) return true; + return (atom0 == mdeltaPairA3d && atom1 == mdeltaPairB3d) || + (atom0 == mdeltaPairB3d && atom1 == mdeltaPairA3d); +} + +PDF3DCalculator::DeltaBondKey +PDF3DCalculator::makeDeltaBondKey(const BaseBondGenerator& bnds) const +{ + const R3::Vector& r = bnds.r01(); + DeltaBondKey key; + key.site0 = bnds.site0(); + key.site1 = bnds.site1(); + key.rx = static_cast(std::lround(r[0] / mdeltaKeyTolerance3d)); + key.ry = static_cast(std::lround(r[1] / mdeltaKeyTolerance3d)); + key.rz = static_cast(std::lround(r[2] / mdeltaKeyTolerance3d)); + return key; +} + +bool PDF3DCalculator::isDeltaEligiblePair(const BaseBondGenerator& bnds) const +{ + if (!menablennDelta3d || mnnDelta3d == 0.0) return false; + return mdeltaEligibleBondKeys.find(this->makeDeltaBondKey(bnds)) != + mdeltaEligibleBondKeys.end(); +} + +double PDF3DCalculator::projectedSigmaAlongBond( + const R3::Matrix& Sigma, + const R3::Vector& r, + double distance) const +{ + if (distance <= 0.0) return 0.0; + R3::Vector e(r[0] / distance, r[1] / distance, r[2] / distance); + R3::Vector tmp = R3::mxvecproduct(Sigma, e); + return R3::dot(e, tmp); +} + +double PDF3DCalculator::positiveDeltaBoundAlongBond( + const R3::Matrix& Sigma, + const R3::Vector& r, + double distance) const +{ + if (distance <= 0.0) return 0.0; + R3::Vector e(r[0] / distance, r[1] / distance, r[2] / distance); + + R3::Vector eigenvalues; + R3::Matrix eigenvectors; + R3::eigen_solve_3x3(Sigma, eigenvalues, eigenvectors); + if (eigenvalues[0] <= 1.0e-10) + return 0.0; + + try + { + R3::Matrix Sigma_inv = R3::inverse(Sigma); + R3::Vector tmp = R3::mxvecproduct(Sigma_inv, e); + const double denom = R3::dot(e, tmp); + return (std::isfinite(denom) && denom > 0.0) ? (1.0 / denom) : 0.0; + } + catch (...) + { + return 0.0; + } +} + +bool PDF3DCalculator::isDistanceDecayDeltaActive() const +{ + return menablennDelta3d && + (mdistanceDelta1_3d > 0.0 || mdistanceDelta2_3d > 0.0); +} + +double PDF3DCalculator::distanceDecayDeltaFraction(double distance) const +{ + if (distance <= 0.0) return 0.0; + return mdistanceDelta1_3d / distance + + mdistanceDelta2_3d / (distance * distance); +} + +double PDF3DCalculator::currentNNDeltaUpperBound() const +{ + if (mnnDelta3dUpperBound <= 0.0) return 0.0; + return museadpscaleSensitivity3d ? + (mnnDelta3dUpperBound * madpScale3d) : mnnDelta3dUpperBound; +} + +bool PDF3DCalculator::isVectorHistogramMode3D() const +{ + return mcalculationmode3d == 1; +} + +void PDF3DCalculator::validateNNDelta3D() const +{ + if (this->isVectorHistogramMode3D()) return; + if (mnnDelta3d < 0.0) + throw std::invalid_argument("3D NN delta must be non-negative."); + if (mnnDelta3d == 0.0 || mdeltaShellRecords.empty()) return; + const double bound = this->currentNNDeltaUpperBound(); + if (bound <= 0.0) + throw std::invalid_argument("3D NN delta cannot be positive without a positive-definite shell bound."); + if (mnnDelta3d >= bound) + throw std::invalid_argument("3D NN delta exceeds the positive-definite shell bound."); +} + +void PDF3DCalculator::validateDistanceDecayDelta3D() const +{ + if (this->isVectorHistogramMode3D()) return; + if (mdistanceDelta1_3d < 0.0 || mdistanceDelta2_3d < 0.0) + throw std::invalid_argument("3D distance-decay delta parameters must be non-negative."); + if (!this->isDistanceDecayDeltaActive() || mdeltaShellRecords.empty()) return; + + for (const DeltaShellRecord& record : mdeltaShellRecords) + { + if (record.distance <= 0.0) continue; + const double f = this->distanceDecayDeltaFraction(record.distance); + if (f <= 0.0) continue; + if (record.projectedSigma <= 0.0 || record.positiveDeltaBound <= 0.0) + continue; + const double delta = record.projectedSigma * f; + const double bound = mnnDeltaPositiveEta3d * record.positiveDeltaBound; + if (bound <= 0.0 || delta >= bound) + throw std::invalid_argument("3D distance-decay delta violates the pair positive-definite covariance bound."); + } +} + +void PDF3DCalculator::buildDeltaEligibleShellTable() +{ + mdeltaEligibleBondKeys.clear(); + mdeltaShellRecords.clear(); + mnnDelta3dUpperBound = 0.0; + + const StructureAdapterPtr& structure = this->getStructure(); + if (!structure) return; + + BaseBondGeneratorPtr bnds = structure->createBondGenerator(); + this->configureBondGenerator(*bnds); + const int cntsites = structure->countSites(); + + for (int i0 = 0; i0 < cntsites; ++i0) + { + bnds->selectAnchorSite(i0); + bnds->selectSiteRange(0, cntsites); + for (bnds->rewind(); !bnds->finished(); bnds->next()) + { + if (bnds->distance() == 0.0) continue; + const std::string& atom0 = structure->siteAtomType(bnds->site0()); + const std::string& atom1 = structure->siteAtomType(bnds->site1()); + if (!this->matchesDeltaPairTypes(atom0, atom1)) continue; + + const R3::Vector& r = bnds->r01(); + R3::Matrix Sigma = bnds->Ucartesian0() + bnds->Ucartesian1(); + if (museadpscaleSensitivity3d) + { + for (int row = 0; row < 3; ++row) + for (int col = 0; col < 3; ++col) + Sigma(row, col) *= madpScale3d; + } + const double projected = this->projectedSigmaAlongBond( + Sigma, r, bnds->distance()); + const double positive_bound = this->positiveDeltaBoundAlongBond( + Sigma, r, bnds->distance()); + + DeltaShellRecord record; + record.site0 = bnds->site0(); + record.site1 = bnds->site1(); + record.r01x = r[0]; + record.r01y = r[1]; + record.r01z = r[2]; + record.distance = bnds->distance(); + record.projectedSigma = projected; + record.positiveDeltaBound = positive_bound; + mdeltaShellRecords.push_back(record); + } + } + + std::vector order(mdeltaShellRecords.size()); + for (size_t i = 0; i < order.size(); ++i) order[i] = i; + std::sort(order.begin(), order.end(), + [&](size_t a, size_t b) { + return mdeltaShellRecords[a].distance < mdeltaShellRecords[b].distance; + }); + + int shell = 0; + double last_distance = 0.0; + bool first = true; + double min_positive_bound = std::numeric_limits::infinity(); + for (size_t idx : order) + { + const double d = mdeltaShellRecords[idx].distance; + if (first || d - last_distance > mdeltaShellTolerance3d) + { + ++shell; + first = false; + } + last_distance = d; + if (shell == mdeltaShellIndex3d) + { + if (mdeltaShellRecords[idx].positiveDeltaBound < min_positive_bound) + min_positive_bound = mdeltaShellRecords[idx].positiveDeltaBound; + DeltaBondKey key; + key.site0 = mdeltaShellRecords[idx].site0; + key.site1 = mdeltaShellRecords[idx].site1; + key.rx = static_cast(std::lround( + mdeltaShellRecords[idx].r01x / mdeltaKeyTolerance3d)); + key.ry = static_cast(std::lround( + mdeltaShellRecords[idx].r01y / mdeltaKeyTolerance3d)); + key.rz = static_cast(std::lround( + mdeltaShellRecords[idx].r01z / mdeltaKeyTolerance3d)); + mdeltaEligibleBondKeys.insert(key); + } + } + + if (min_positive_bound < std::numeric_limits::infinity() && + min_positive_bound > 0.0) + { + mnnDelta3dUpperBound = mnnDeltaPositiveEta3d * min_positive_bound; + } +} + +R3::Matrix PDF3DCalculator::effectivePairCovariance( + const BaseBondGenerator& bnds, + const R3::Matrix& U_i, + const R3::Matrix& U_j) const +{ + R3::Matrix Sigma = U_i + U_j; + if (museadpscaleSensitivity3d) + { + for (int row = 0; row < 3; ++row) + for (int col = 0; col < 3; ++col) + Sigma(row, col) *= madpScale3d; + } + + const double d = bnds.distance(); + if (d <= 0.0) return Sigma; + + R3::Vector e(bnds.r01()[0] / d, bnds.r01()[1] / d, bnds.r01()[2] / d); + double delta_to_subtract = 0.0; + if (this->isDistanceDecayDeltaActive()) + { + const StructureAdapterPtr& structure = this->getStructure(); + if (!structure) return Sigma; + const std::string& atom0 = structure->siteAtomType(bnds.site0()); + const std::string& atom1 = structure->siteAtomType(bnds.site1()); + if (!this->matchesDeltaPairTypes(atom0, atom1)) return Sigma; + + const double fraction = this->distanceDecayDeltaFraction(d); + if (fraction <= 0.0) return Sigma; + const double projected = this->projectedSigmaAlongBond(Sigma, bnds.r01(), d); + delta_to_subtract = projected * fraction; + } + else + { + if (!this->isDeltaEligiblePair(bnds)) return Sigma; + delta_to_subtract = mnnDelta3d; + } + + if (delta_to_subtract <= 0.0) return Sigma; + + const double positive_bound = this->positiveDeltaBoundAlongBond(Sigma, bnds.r01(), d); + if (positive_bound <= 0.0) return Sigma; + + const double local_bound = mnnDeltaPositiveEta3d * positive_bound; + if (local_bound <= 0.0 || delta_to_subtract >= local_bound) + throw std::invalid_argument("3D delta violates the pair positive-definite covariance bound."); + + for (int row = 0; row < 3; ++row) + { + for (int col = 0; col < 3; ++col) + { + Sigma(row, col) -= delta_to_subtract * e[row] * e[col]; + } + } + return Sigma; +} + +size_t PDF3DCalculator::coordToIndex(double x, double y, double z) const +{ + double halfspan = (mnbins / 2) * mdr; + if (fabs(x) > halfspan || fabs(y) > halfspan || fabs(z) > halfspan) + return mgrid3d.size(); // out of bounds + + int ix = static_cast(std::lround((x + halfspan) / mdr)); + int iy = static_cast(std::lround((y + halfspan) / mdr)); + int iz = static_cast(std::lround((z + halfspan) / mdr)); + + ix = std::max(0, std::min(ix, mnbins - 1)); + iy = std::max(0, std::min(iy, mnbins - 1)); + iz = std::max(0, std::min(iz, mnbins - 1)); + + return static_cast((iz * mnbins + iy) * mnbins + ix); +} + +void PDF3DCalculator::indexToCoord(size_t idx, double& x, double& y, double& z) const +{ + size_t iz = idx / (static_cast(mnbins) * mnbins); + size_t rem = idx % (static_cast(mnbins) * mnbins); + size_t iy = rem / mnbins; + size_t ix = rem % mnbins; + + double halfspan = (mnbins / 2) * mdr; + x = ix * mdr - halfspan; + y = iy * mdr - halfspan; + z = iz * mdr - halfspan; +} + +} // namespace srreal +} // namespace diffpy + +#include +DIFFPY_INSTANTIATE_SERIALIZATION(diffpy::srreal::PDF3DCalculator) diff --git a/src/diffpy/srreal/PDF3DCalculator.hpp b/src/diffpy/srreal/PDF3DCalculator.hpp new file mode 100644 index 00000000..c9c22515 --- /dev/null +++ b/src/diffpy/srreal/PDF3DCalculator.hpp @@ -0,0 +1,257 @@ +/************************************************************************************* +* +* libdiffpy by DANSE Diffraction group +* Simon J. L. Billinge +* (c) 2009 The Trustees of Columbia University +* in the City of New York. All rights reserved. +* +* File coded by: Hai Qiao +* +* See AUTHORS.txt for a list of people who contributed. +* See LICENSE_DANSE.txt for license information. +* +/************************************************************************************* +* +* class PDF3DCalculator -- 3D real-space pair distribution function calculator +* +*************************************************************************************/ + +#ifndef PDF3DCALCULATOR_HPP_INCLUDED +#define PDF3DCALCULATOR_HPP_INCLUDED + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace diffpy { +namespace srreal { + +class PDF3DCalculator : public PDFCalculator +{ +public: + // constructor + PDF3DCalculator(); + + // Public interface to retrieve 3D PDF data + // Returns a flat vector of [x0, y0, z0, G0, x1, y1, z1, G1, ...] + QuantityType get3DPDF() const; + + // Export dense 3D grid (nz, ny, nx) directly to binary file. + // usefloat32=true writes float32, otherwise float64. + // applypost=true applies q-window/rdf_scale/rho0/qdamp before export. + void exportGrid3DBinary(const std::string& path, bool usefloat32 = true, bool applypost = true) const; + QuantityType getRadialHistogram3D() const; + + // Grid configuration + void setGridStep(double dr); + double getGridStep() const; + + // Blocked accumulation configuration + void setAccumBlockSize(int bs); + int getAccumBlockSize() const; + + void setApplyRho0Background3D(bool); + bool getApplyRho0Background3D() const; + + void setUseCQWindow3D(bool v); + bool getUseCQWindow3D() const; + + // Nearest-neighbor correlated-motion correction for anisotropic 3D PDF models. + // Boolean switches are also exposed as double attributes with 0/1 + // values because libdiffpy's generic attribute system is double-only. + void setEnableNNDelta3D(bool); + bool getEnableNNDelta3D() const; + double getEnableNNDelta3DAttr() const; + void setEnableNNDelta3DAttr(double); + + void setNNDelta3D(double); + const double& getNNDelta3D() const; + double getNNDelta3DUpperBound() const; + const double& getNNDeltaPositiveEta3D() const; + void setNNDeltaPositiveEta3D(double); + void setDistanceDelta1_3D(double); + const double& getDistanceDelta1_3D() const; + void setDistanceDelta2_3D(double); + const double& getDistanceDelta2_3D() const; + + void setDeltaPairTypes3D(const std::string&, const std::string&); + const std::string& getDeltaPairA3D() const; + const std::string& getDeltaPairB3D() const; + void setDeltaShellIndex3D(int); + int getDeltaShellIndex3D() const; + double getDeltaShellIndex3DAttr() const; + void setDeltaShellIndex3DAttr(double); + const double& getDeltaShellTolerance3D() const; + void setDeltaShellTolerance3D(double); + const double& getDeltaKeyTolerance3D() const; + void setDeltaKeyTolerance3D(double); + int getDeltaEligiblePairCount3D() const; + + void setUseADPScaleSensitivity3D(bool); + bool getUseADPScaleSensitivity3D() const; + double getUseADPScaleSensitivity3DAttr() const; + void setUseADPScaleSensitivity3DAttr(double); + void setADPScale3D(double); + const double& getADPScale3D() const; + + void setRho0BackgroundScale3D(double); + const double& getRho0BackgroundScale3D() const; + + void setCalculationMode3D(int); + int getCalculationMode3D() const; + double getCalculationMode3DAttr() const; + void setCalculationMode3DAttr(double); + + void setHistogramWeightMode3D(int); + int getHistogramWeightMode3D() const; + double getHistogramWeightMode3DAttr() const; + void setHistogramWeightMode3DAttr(double); + +protected: + // Override PairQuantity virtual methods + virtual void resetValue() override; + virtual void addPairContribution(const BaseBondGenerator& bnds, int) override; + +private: + struct DeltaBondKey + { + int site0; + int site1; + long rx; + long ry; + long rz; + bool operator==(const DeltaBondKey&) const; + }; + + struct DeltaBondKeyHash + { + size_t operator()(const DeltaBondKey&) const; + }; + + struct DeltaShellRecord + { + int site0; + int site1; + double r01x; + double r01y; + double r01z; + double distance; + double projectedSigma; + double positiveDeltaBound; + }; + + // Helper to add anisotropic Gaussian to the grid + void addAnisotropicGaussianToGrid(const R3::Vector& r_ij, const R3::Matrix& Sigma, double sfprod); + void addVectorHistogramToGrid(const R3::Vector& r_ij, double distance, double weight); + + // Helper: map 3D position to linear index + size_t coordToIndex(double x, double y, double z) const; + void indexToCoord(size_t idx, double& x, double& y, double& z) const; + + void applyQWindow3D(std::vector& grid) const; + double computeRho0Background() const; + void buildDeltaEligibleShellTable(); + bool isDeltaEligiblePair(const BaseBondGenerator&) const; + bool matchesDeltaPairTypes(const std::string&, const std::string&) const; + DeltaBondKey makeDeltaBondKey(const BaseBondGenerator&) const; + R3::Matrix effectivePairCovariance( + const BaseBondGenerator&, + const R3::Matrix&, + const R3::Matrix&) const; + double projectedSigmaAlongBond( + const R3::Matrix&, + const R3::Vector&, + double distance) const; + double positiveDeltaBoundAlongBond( + const R3::Matrix&, + const R3::Vector&, + double distance) const; + bool isDistanceDecayDeltaActive() const; + double distanceDecayDeltaFraction(double distance) const; + void validateDistanceDecayDelta3D() const; + void validateNNDelta3D() const; + double currentNNDeltaUpperBound() const; + bool isVectorHistogramMode3D() const; + + // Data members + double mdr; // grid spacing (assumes cubic grid centered at origin) + int mnbins; // number of bins per dimension (odd number, center at 0) + std::vector mgrid3d; // flattened 3D histogram (size = mnbins^3) + std::vector mradialhistogram3d; + int maccumblocksize; // block size for tiled accumulation loops + + bool mapplyrho0background3d; + bool musecqwindow3d; + int mcalculationmode3d; + int mhistogramweightmode3d; + + bool menablennDelta3d; + double mnnDelta3d; + double mnnDelta3dUpperBound; + double mnnDeltaPositiveEta3d; + double mdistanceDelta1_3d; + double mdistanceDelta2_3d; + std::string mdeltaPairA3d; + std::string mdeltaPairB3d; + int mdeltaShellIndex3d; + double mdeltaShellTolerance3d; + double mdeltaKeyTolerance3d; + bool museadpscaleSensitivity3d; + double madpScale3d; + double mrho0backgroundscale3d; + std::unordered_set mdeltaEligibleBondKeys; + std::vector mdeltaShellRecords; + + std::vector msfCache; + double mtotaloccupancy; + double msfaverage; + + // Internal constants + static constexpr double DEFAULT_RMAX_3D = 10.0; // Angstrom + static constexpr double DEFAULT_GRID_STEP = 0.1; // Angstrom + + // serialization + friend class boost::serialization::access; + template + void serialize(Archive& ar, const unsigned int version) + { + using boost::serialization::base_object; + ar & base_object(*this); + ar & mdr; + ar & mnbins; + ar & mgrid3d; + ar & mradialhistogram3d; + ar & maccumblocksize; + ar & mapplyrho0background3d; + ar & musecqwindow3d; + ar & mcalculationmode3d; + ar & mhistogramweightmode3d; + ar & menablennDelta3d; + ar & mnnDelta3d; + ar & mnnDelta3dUpperBound; + ar & mnnDeltaPositiveEta3d; + ar & mdistanceDelta1_3d; + ar & mdistanceDelta2_3d; + ar & mdeltaPairA3d; + ar & mdeltaPairB3d; + ar & mdeltaShellIndex3d; + ar & mdeltaShellTolerance3d; + ar & mdeltaKeyTolerance3d; + ar & museadpscaleSensitivity3d; + ar & madpScale3d; + ar & mrho0backgroundscale3d; + ar & msfCache; + ar & mtotaloccupancy; + ar & msfaverage; + } +}; + +} // namespace srreal +} // namespace diffpy + +#endif // PDF3DCALCULATOR_HPP_INCLUDED diff --git a/src/diffpy/srreal/PDFCalculator.cpp b/src/diffpy/srreal/PDFCalculator.cpp index 47c7189a..adb3cb4a 100644 --- a/src/diffpy/srreal/PDFCalculator.cpp +++ b/src/diffpy/srreal/PDFCalculator.cpp @@ -677,6 +677,14 @@ double PDFCalculator::sfAverage() const } +double PDFCalculator::getPartialPDFScale() const +{ + const double totocc = mstructure_cache.totaloccupancy; + return (totocc == 0.0) ? 0.0 : + (mstructure_cache.activeoccupancy / totocc); +} + + void PDFCalculator::cacheStructureData() { int cntsites = this->countSites(); diff --git a/src/diffpy/srreal/PDFCalculator.hpp b/src/diffpy/srreal/PDFCalculator.hpp index 79220e0a..bcce3850 100644 --- a/src/diffpy/srreal/PDFCalculator.hpp +++ b/src/diffpy/srreal/PDFCalculator.hpp @@ -114,6 +114,8 @@ class PDFCalculator : // support for PQEvaluatorOptimized virtual void stashPartialValue(); virtual void restorePartialValue(); + /// activeoccupancy / totaloccupancy used by baseline background term + double getPartialPDFScale() const; private: diff --git a/src/diffpy/srreal/R3linalg.cpp b/src/diffpy/srreal/R3linalg.cpp index 19a70cc7..337acd4e 100644 --- a/src/diffpy/srreal/R3linalg.cpp +++ b/src/diffpy/srreal/R3linalg.cpp @@ -16,6 +16,8 @@ * *****************************************************************************/ +#include +#include #include #include #include @@ -82,6 +84,88 @@ const Matrix& inverse(const Matrix& A) } +void eigen_solve_3x3(const Matrix& A, Vector& w, Matrix& V) +{ + V = identity(); + Matrix m = A; + + const int max_iter = 50; + const double eps = 1e-10; + + for (int iter = 0; iter < max_iter; ++iter) + { + double max_off_diag = 0.0; + int p = 0; + int q = 1; + + for (int i = 0; i < Ndim; ++i) + { + for (int j = i + 1; j < Ndim; ++j) + { + if (std::abs(m(i, j)) > max_off_diag) + { + max_off_diag = std::abs(m(i, j)); + p = i; + q = j; + } + } + } + + if (max_off_diag < eps) break; + + double phi = 0.5 * std::atan2( + 2.0 * m(p, q), m(q, q) - m(p, p)); + double c = std::cos(phi); + double s = std::sin(phi); + + double m_pp = m(p, p); + double m_qq = m(q, q); + double m_pq = m(p, q); + + m(p, p) = c * c * m_pp - 2.0 * s * c * m_pq + s * s * m_qq; + m(q, q) = s * s * m_pp + 2.0 * s * c * m_pq + c * c * m_qq; + m(p, q) = 0.0; + m(q, p) = 0.0; + + for (int i = 0; i < Ndim; ++i) + { + if (i == p || i == q) continue; + double m_ip = m(i, p); + double m_iq = m(i, q); + m(i, p) = c * m_ip - s * m_iq; + m(p, i) = m(i, p); + m(i, q) = s * m_ip + c * m_iq; + m(q, i) = m(i, q); + } + + for (int i = 0; i < Ndim; ++i) + { + double v_ip = V(i, p); + double v_iq = V(i, q); + V(i, p) = c * v_ip - s * v_iq; + V(i, q) = s * v_ip + c * v_iq; + } + } + + w[0] = m(0, 0); + w[1] = m(1, 1); + w[2] = m(2, 2); + + for (int i = 0; i < Ndim - 1; ++i) + { + for (int j = 0; j < Ndim - 1 - i; ++j) + { + if (w[j] <= w[j + 1]) continue; + std::swap(w[j], w[j + 1]); + for (int k = 0; k < Ndim; ++k) + { + std::swap(V(k, j), V(k, j + 1)); + } + } + } +} + + size_t hash_value(const Vector& v) { return boost::hash_range(v.begin(), v.end()); diff --git a/src/diffpy/srreal/R3linalg.hpp b/src/diffpy/srreal/R3linalg.hpp index bffd5c7c..91888346 100644 --- a/src/diffpy/srreal/R3linalg.hpp +++ b/src/diffpy/srreal/R3linalg.hpp @@ -150,6 +150,8 @@ const Matrix& identity(); const Matrix& zeromatrix(); double determinant(const Matrix& A); const Matrix& inverse(const Matrix& A); +void eigen_solve_3x3(const Matrix& A, Vector& eigenvalues, + Matrix& eigenvectors); const Vector& floor(const Vector&); template double norm(const V&); From 705b955049eb579ffca108a58f7ed474f7c4fd09 Mon Sep 17 00:00:00 2001 From: joeseaer <2969111832@qq.com> Date: Thu, 9 Jul 2026 01:29:43 +0800 Subject: [PATCH 2/4] MAINT: address PDF3DCalculator review feedback --- src/diffpy/srreal/PDF3DCalculator.cpp | 165 +++++++---------------- src/diffpy/srreal/PDF3DCalculator.hpp | 26 +--- src/diffpy/srreal/PDFCalculator.cpp | 19 ++- src/diffpy/srreal/PDFCalculator.hpp | 10 +- src/tests/TestPDF3DCalculator.hpp | 180 ++++++++++++++++++++++++++ src/tests/TestR3linalg.hpp | 47 +++++++ 6 files changed, 294 insertions(+), 153 deletions(-) create mode 100644 src/tests/TestPDF3DCalculator.hpp diff --git a/src/diffpy/srreal/PDF3DCalculator.cpp b/src/diffpy/srreal/PDF3DCalculator.cpp index e839ab35..e3bbca03 100644 --- a/src/diffpy/srreal/PDF3DCalculator.cpp +++ b/src/diffpy/srreal/PDF3DCalculator.cpp @@ -1,21 +1,3 @@ -/************************************************************************************* -* -* libdiffpy by DANSE Diffraction group -* Simon J. L. Billinge -* (c) 2009 The Trustees of Columbia University -* in the City of New York. All rights reserved. -* -* File coded by: Hai Qiao -* -* See AUTHORS.txt for a list of people who contributed. -* See LICENSE_DANSE.txt for license information. -* -/************************************************************************************* -* -* class PDF3DCalculator -- Implementation of PDF3DCalculator -* -*************************************************************************************/ - #include #include #include @@ -65,9 +47,7 @@ PDF3DCalculator::PDF3DCalculator() : mdeltaKeyTolerance3d(1.0e-6), museadpscaleSensitivity3d(false), madpScale3d(1.0), - mrho0backgroundscale3d(1.0), - mtotaloccupancy(0.0), - msfaverage(0.0) + mrho0backgroundscale3d(1.0) { this->registerDoubleAttribute("enable_nn_delta3d", this, @@ -444,44 +424,7 @@ QuantityType PDF3DCalculator::get3DPDF() const std::vector grid = mgrid3d; if (!this->isVectorHistogramMode3D()) { - if (musecqwindow3d) applyQWindow3D(grid); - - const double rdf_scale = (mtotaloccupancy * msfaverage == 0.0) ? 0.0 : - 1.0 / (mtotaloccupancy * msfaverage * msfaverage); - if (rdf_scale != 1.0) - { - for (double& val : grid) val *= rdf_scale; - } - - if (mapplyrho0background3d) - { - const double rho0_bg = computeRho0Background(); - if (rho0_bg != 0.0) - { - for (double& val : grid) val -= rho0_bg; - } - } - - double qdamp = 0.0; - try - { - qdamp = this->getEnvelopeByType("qresolution")->getDoubleAttr("qdamp"); - } - catch (...) - { - qdamp = 0.0; - } - if (qdamp > 0.0) - { - for (size_t i = 0; i < grid.size(); ++i) - { - if (grid[i] == 0.0) continue; - double x, y, z; - indexToCoord(i, x, y, z); - const double r = sqrt(x * x + y * y + z * z); - grid[i] *= exp(-0.5 * (r * qdamp) * (r * qdamp)); - } - } + this->applyPostProcessing3D(grid); } size_t nonzero_count = 0; @@ -523,34 +466,7 @@ void PDF3DCalculator::exportGrid3DBinary(const std::string& path, bool usefloat3 std::vector grid = mgrid3d; if (applypost) { - if (musecqwindow3d) applyQWindow3D(grid); - - const double rdf_scale = (mtotaloccupancy * msfaverage == 0.0) ? 0.0 : - 1.0 / (mtotaloccupancy * msfaverage * msfaverage); - if (rdf_scale != 1.0) - for (double& val : grid) val *= rdf_scale; - - if (mapplyrho0background3d) - { - const double rho0_bg = computeRho0Background(); - if (rho0_bg != 0.0) - for (double& val : grid) val -= rho0_bg; - } - - double qdamp = 0.0; - try { qdamp = this->getEnvelopeByType("qresolution")->getDoubleAttr("qdamp"); } - catch (...) { qdamp = 0.0; } - if (qdamp > 0.0) - { - for (size_t i = 0; i < grid.size(); ++i) - { - if (grid[i] == 0.0) continue; - double x, y, z; - indexToCoord(i, x, y, z); - const double r = sqrt(x * x + y * y + z * z); - grid[i] *= exp(-0.5 * (r * qdamp) * (r * qdamp)); - } - } + this->applyPostProcessing3D(grid); } std::ofstream ofs(path.c_str(), std::ios::binary | std::ios::trunc); @@ -578,6 +494,47 @@ void PDF3DCalculator::exportGrid3DBinary(const std::string& path, bool usefloat3 if (!ofs) throw std::runtime_error("Failed while writing output file: " + path); } +void PDF3DCalculator::applyPostProcessing3D(std::vector& grid) const +{ + if (musecqwindow3d) applyQWindow3D(grid); + + const double rdf_scale = this->getRDFScale(); + if (rdf_scale != 1.0) + { + for (double& val : grid) val *= rdf_scale; + } + + if (mapplyrho0background3d) + { + const double rho0_bg = computeRho0Background(); + if (rho0_bg != 0.0) + { + for (double& val : grid) val -= rho0_bg; + } + } + + double qdamp = 0.0; + try + { + qdamp = this->getEnvelopeByType("qresolution")->getDoubleAttr("qdamp"); + } + catch (...) + { + qdamp = 0.0; + } + if (qdamp > 0.0) + { + for (size_t i = 0; i < grid.size(); ++i) + { + if (grid[i] == 0.0) continue; + double x, y, z; + indexToCoord(i, x, y, z); + const double r = sqrt(x * x + y * y + z * z); + grid[i] *= exp(-0.5 * (r * qdamp) * (r * qdamp)); + } + } +} + double PDF3DCalculator::computeRho0Background() const { const StructureAdapterPtr& structure = this->getStructure(); @@ -721,33 +678,6 @@ void PDF3DCalculator::applyQWindow3D(std::vector& grid) const void PDF3DCalculator::resetValue() { - const StructureAdapterPtr& structure = this->getStructure(); - if (!structure) - { - msfCache.clear(); - mtotaloccupancy = 0.0; - msfaverage = 0.0; - } - else - { - int nsite = this->countSites(); - msfCache.resize(nsite); - const auto& sftable = this->getScatteringFactorTable(); - - mtotaloccupancy = structure->totalOccupancy(); - double totsf = 0.0; - for (int i = 0; i < nsite; ++i) - { - std::string atom_type = structure->siteAtomType(i); - const double sf = sftable->lookup(atom_type); - const double occ = structure->siteOccupancy(i); - const double mult = structure->siteMultiplicity(i); - msfCache[i] = sf * occ; - totsf += msfCache[i] * mult; - } - msfaverage = (mtotaloccupancy == 0.0) ? 0.0 : (totsf / mtotaloccupancy); - } - // Ensure odd number of bins so origin is centered. mnbins = static_cast(2 * ceil(this->getRmax() / mdr)) + 1; size_t total = static_cast(mnbins) * mnbins * mnbins; @@ -771,13 +701,14 @@ void PDF3DCalculator::addPairContribution(const BaseBondGenerator& bnds, int sum int i0 = bnds.site0(); int i1 = bnds.site1(); + int cntsites = this->countSites(); - if (i0 >= static_cast(msfCache.size()) || i1 >= static_cast(msfCache.size())) + if (i0 >= cntsites || i1 >= cntsites) return; const R3::Vector& rvec = bnds.r01(); const double pairscale = bnds.multiplicity() * static_cast(summationscale); - double sfprod = msfCache[i0] * msfCache[i1] * pairscale; + double sfprod = this->sfSite(i0) * this->sfSite(i1) * pairscale; if (this->isVectorHistogramMode3D()) { @@ -860,7 +791,7 @@ void PDF3DCalculator::addAnisotropicGaussianToGrid(const R3::Vector& r_ij, const double min_val = center_val - cutoff_radius; double max_val = center_val + cutoff_radius; - // Node-centered grid: points at k * mdr + // Node-centered grid: points at k * md int start = static_cast(ceil((min_val + halfspan) / mdr)); int end = static_cast(floor((max_val + halfspan) / mdr)); diff --git a/src/diffpy/srreal/PDF3DCalculator.hpp b/src/diffpy/srreal/PDF3DCalculator.hpp index c9c22515..b14abc97 100644 --- a/src/diffpy/srreal/PDF3DCalculator.hpp +++ b/src/diffpy/srreal/PDF3DCalculator.hpp @@ -1,21 +1,3 @@ -/************************************************************************************* -* -* libdiffpy by DANSE Diffraction group -* Simon J. L. Billinge -* (c) 2009 The Trustees of Columbia University -* in the City of New York. All rights reserved. -* -* File coded by: Hai Qiao -* -* See AUTHORS.txt for a list of people who contributed. -* See LICENSE_DANSE.txt for license information. -* -/************************************************************************************* -* -* class PDF3DCalculator -- 3D real-space pair distribution function calculator -* -*************************************************************************************/ - #ifndef PDF3DCALCULATOR_HPP_INCLUDED #define PDF3DCALCULATOR_HPP_INCLUDED @@ -153,6 +135,7 @@ class PDF3DCalculator : public PDFCalculator size_t coordToIndex(double x, double y, double z) const; void indexToCoord(size_t idx, double& x, double& y, double& z) const; + void applyPostProcessing3D(std::vector& grid) const; void applyQWindow3D(std::vector& grid) const; double computeRho0Background() const; void buildDeltaEligibleShellTable(); @@ -207,10 +190,6 @@ class PDF3DCalculator : public PDFCalculator std::unordered_set mdeltaEligibleBondKeys; std::vector mdeltaShellRecords; - std::vector msfCache; - double mtotaloccupancy; - double msfaverage; - // Internal constants static constexpr double DEFAULT_RMAX_3D = 10.0; // Angstrom static constexpr double DEFAULT_GRID_STEP = 0.1; // Angstrom @@ -245,9 +224,6 @@ class PDF3DCalculator : public PDFCalculator ar & museadpscaleSensitivity3d; ar & madpScale3d; ar & mrho0backgroundscale3d; - ar & msfCache; - ar & mtotaloccupancy; - ar & msfaverage; } }; diff --git a/src/diffpy/srreal/PDFCalculator.cpp b/src/diffpy/srreal/PDFCalculator.cpp index adb3cb4a..e2df550f 100644 --- a/src/diffpy/srreal/PDFCalculator.cpp +++ b/src/diffpy/srreal/PDFCalculator.cpp @@ -187,10 +187,7 @@ QuantityType PDFCalculator::getExtendedPDF() const QuantityType PDFCalculator::getExtendedRDF() const { QuantityType rdf(this->countExtendedPoints()); - const double& totocc = mstructure_cache.totaloccupancy; - double sfavg = this->sfAverage(); - double rdf_scale = (totocc * sfavg == 0.0) ? 0.0 : - 1.0 / (totocc * sfavg * sfavg); + double rdf_scale = this->getRDFScale(); QuantityType::iterator iirdf = rdf.begin(); QuantityType::const_iterator iival, iival_last; iival = this->value().begin() + @@ -491,9 +488,7 @@ void PDFCalculator::resetValue() // when applicable, configure linear baseline if (this->getBaseline()->type() == "linear") { - double partialpdfscale = - (0.0 == mstructure_cache.totaloccupancy) ? 0.0 : - mstructure_cache.activeoccupancy / mstructure_cache.totaloccupancy; + double partialpdfscale = this->getPartialPDFScale(); double pnumdensity = partialpdfscale * mstructure->numberDensity(); PDFBaseline& bl = *(this->getBaseline()); bl.setDoubleAttr("slope", -4 * M_PI * pnumdensity); @@ -685,6 +680,16 @@ double PDFCalculator::getPartialPDFScale() const } +double PDFCalculator::getRDFScale() const +{ + const double& totocc = mstructure_cache.totaloccupancy; + double sfavg = this->sfAverage(); + double rv = (totocc * sfavg == 0.0) ? 0.0 : + 1.0 / (totocc * sfavg * sfavg); + return rv; +} + + void PDFCalculator::cacheStructureData() { int cntsites = this->countSites(); diff --git a/src/diffpy/srreal/PDFCalculator.hpp b/src/diffpy/srreal/PDFCalculator.hpp index bcce3850..5d5eba4b 100644 --- a/src/diffpy/srreal/PDFCalculator.hpp +++ b/src/diffpy/srreal/PDFCalculator.hpp @@ -116,6 +116,12 @@ class PDFCalculator : virtual void restorePartialValue(); /// activeoccupancy / totaloccupancy used by baseline background term double getPartialPDFScale() const; + /// RDF scale from the cached total occupancy and average scattering factor + double getRDFScale() const; + /// effective scattering factor at a given site scaled by occupancy + const double& sfSite(int) const; + /// average scattering factor + double sfAverage() const; private: @@ -147,10 +153,6 @@ class PDFCalculator : void cutRipplePoints(QuantityType& y) const; // structure factors - fast lookup by site index - /// effective scattering factor at a given site scaled by occupancy - const double& sfSite(int) const; - /// average scattering factor - double sfAverage() const; void cacheStructureData(); void cacheRlimitsData(); diff --git a/src/tests/TestPDF3DCalculator.hpp b/src/tests/TestPDF3DCalculator.hpp new file mode 100644 index 00000000..dfc58d7e --- /dev/null +++ b/src/tests/TestPDF3DCalculator.hpp @@ -0,0 +1,180 @@ +/***************************************************************************** +* +* class TestPDF3DCalculator -- unit tests for PDF3DCalculator class +* +*****************************************************************************/ + +#include + +#include +#include +#include +#include + +#include +#include +#include + +using namespace std; +using namespace diffpy::srreal; + +class TestPDF3DCalculator : public CxxTest::TestSuite +{ + private: + + boost::shared_ptr mpdfc; + AtomicStructureAdapterPtr mstru2; + double meps; + + public: + + void setUp() + { + meps = diffpy::mathutils::SQRT_DOUBLE_EPS; + mpdfc.reset(new PDF3DCalculator); + + mstru2 = boost::make_shared(); + Atom ai; + ai.atomtype = "Ni"; + ai.xyz_cartn[0] = 0.0; + ai.xyz_cartn[1] = 0.0; + ai.xyz_cartn[2] = 0.0; + ai.uij_cartn = R3::identity(); + ai.uij_cartn(0, 0) = ai.uij_cartn(1, 1) = + ai.uij_cartn(2, 2) = 0.004; + mstru2->append(ai); + ai.xyz_cartn[0] = 1.0; + mstru2->append(ai); + } + + + void test_defaults_and_attributes() + { + TS_ASSERT_DELTA(0.1, mpdfc->getGridStep(), meps); + TS_ASSERT_EQUALS(32, mpdfc->getAccumBlockSize()); + TS_ASSERT_EQUALS(true, mpdfc->getApplyRho0Background3D()); + TS_ASSERT_EQUALS(true, mpdfc->getUseCQWindow3D()); + TS_ASSERT_EQUALS(false, mpdfc->getEnableNNDelta3D()); + TS_ASSERT_EQUALS(0, mpdfc->getCalculationMode3D()); + TS_ASSERT_EQUALS(0, mpdfc->getHistogramWeightMode3D()); + + TS_ASSERT_DELTA(0.0, + mpdfc->getDoubleAttr("enable_nn_delta3d"), meps); + TS_ASSERT_DELTA(0.0, mpdfc->getDoubleAttr("nn_delta3d"), meps); + TS_ASSERT_DELTA(0.9, + mpdfc->getDoubleAttr("nn_delta_positive_eta3d"), meps); + TS_ASSERT_DELTA(0.0, mpdfc->getDoubleAttr("delta1_3d"), meps); + TS_ASSERT_DELTA(0.0, mpdfc->getDoubleAttr("delta2_3d"), meps); + TS_ASSERT_DELTA(1.0, + mpdfc->getDoubleAttr("delta_shell_index3d"), meps); + TS_ASSERT_DELTA(1.0, mpdfc->getDoubleAttr("adp_scale3d"), meps); + + mpdfc->setDoubleAttr("calculation_mode3d", 1.0); + mpdfc->setDoubleAttr("histogram_weight_mode3d", 1.0); + mpdfc->setDoubleAttr("enable_nn_delta3d", 1.0); + TS_ASSERT_EQUALS(1, mpdfc->getCalculationMode3D()); + TS_ASSERT_EQUALS(1, mpdfc->getHistogramWeightMode3D()); + TS_ASSERT_EQUALS(true, mpdfc->getEnableNNDelta3D()); + } + + + void test_invalid_configuration() + { + TS_ASSERT_THROWS(mpdfc->setGridStep(0.0), invalid_argument); + TS_ASSERT_THROWS(mpdfc->setAccumBlockSize(0), invalid_argument); + TS_ASSERT_THROWS(mpdfc->setCalculationMode3D(2), invalid_argument); + TS_ASSERT_THROWS(mpdfc->setHistogramWeightMode3D(2), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setDoubleAttr("calculation_mode3d", 0.5), + invalid_argument); + TS_ASSERT_THROWS(mpdfc->setADPScale3D(0.0), invalid_argument); + TS_ASSERT_THROWS(mpdfc->setRho0BackgroundScale3D(-1.0), + invalid_argument); + } + + + void test_vector_histogram_eval() + { + mpdfc->setRmax(2.0); + mpdfc->setGridStep(0.5); + mpdfc->setCalculationMode3D(1); + mpdfc->setHistogramWeightMode3D(1); + mpdfc->eval(mstru2); + + QuantityType pdf3d = mpdfc->get3DPDF(); + TS_ASSERT(!pdf3d.empty()); + TS_ASSERT_EQUALS(0u, pdf3d.size() % 4); + + QuantityType rh = mpdfc->getRadialHistogram3D(); + size_t nr = static_cast( + std::ceil(mpdfc->getRmax() / mpdfc->getGridStep())) + 1; + TS_ASSERT_EQUALS(2 * nr, rh.size()); + + bool found_distance_one = false; + for (size_t i = 0; i < rh.size(); i += 2) + { + if (std::fabs(rh[i] - 1.0) < meps && rh[i + 1] != 0.0) + { + found_distance_one = true; + break; + } + } + TS_ASSERT(found_distance_one); + } + + + void test_serialization() + { + mpdfc->setRmax(1.0); + mpdfc->setGridStep(0.25); + mpdfc->setAccumBlockSize(7); + mpdfc->setApplyRho0Background3D(false); + mpdfc->setUseCQWindow3D(false); + mpdfc->setCalculationMode3D(1); + mpdfc->setHistogramWeightMode3D(1); + mpdfc->setEnableNNDelta3D(true); + mpdfc->setNNDelta3D(0.02); + mpdfc->setNNDeltaPositiveEta3D(0.8); + mpdfc->setDistanceDelta1_3D(0.03); + mpdfc->setDistanceDelta2_3D(0.04); + mpdfc->setDeltaPairTypes3D("Ni", "Ni"); + mpdfc->setDeltaShellIndex3D(2); + mpdfc->setDeltaShellTolerance3D(0.06); + mpdfc->setDeltaKeyTolerance3D(1.0e-5); + mpdfc->setUseADPScaleSensitivity3D(true); + mpdfc->setADPScale3D(1.2); + mpdfc->setRho0BackgroundScale3D(0.5); + + stringstream storage(ios::in | ios::out | ios::binary); + diffpy::serialization::oarchive oa(storage, ios::binary); + oa << mpdfc; + diffpy::serialization::iarchive ia(storage, ios::binary); + boost::shared_ptr pdfc1; + ia >> pdfc1; + + TS_ASSERT_DIFFERS(pdfc1.get(), mpdfc.get()); + TS_ASSERT_DELTA(0.25, pdfc1->getGridStep(), meps); + TS_ASSERT_EQUALS(7, pdfc1->getAccumBlockSize()); + TS_ASSERT_EQUALS(false, pdfc1->getApplyRho0Background3D()); + TS_ASSERT_EQUALS(false, pdfc1->getUseCQWindow3D()); + TS_ASSERT_EQUALS(1, pdfc1->getCalculationMode3D()); + TS_ASSERT_EQUALS(1, pdfc1->getHistogramWeightMode3D()); + TS_ASSERT_EQUALS(true, pdfc1->getEnableNNDelta3D()); + TS_ASSERT_DELTA(0.02, pdfc1->getNNDelta3D(), meps); + TS_ASSERT_DELTA(0.8, pdfc1->getNNDeltaPositiveEta3D(), meps); + TS_ASSERT_DELTA(0.03, pdfc1->getDistanceDelta1_3D(), meps); + TS_ASSERT_DELTA(0.04, pdfc1->getDistanceDelta2_3D(), meps); + TS_ASSERT_EQUALS(string("Ni"), pdfc1->getDeltaPairA3D()); + TS_ASSERT_EQUALS(string("Ni"), pdfc1->getDeltaPairB3D()); + TS_ASSERT_EQUALS(2, pdfc1->getDeltaShellIndex3D()); + TS_ASSERT_DELTA(0.06, pdfc1->getDeltaShellTolerance3D(), meps); + TS_ASSERT_DELTA(1.0e-5, pdfc1->getDeltaKeyTolerance3D(), meps); + TS_ASSERT_EQUALS(true, pdfc1->getUseADPScaleSensitivity3D()); + TS_ASSERT_DELTA(1.2, pdfc1->getADPScale3D(), meps); + TS_ASSERT_DELTA(0.5, pdfc1->getRho0BackgroundScale3D(), meps); + } + +}; // class TestPDF3DCalculator + +// End of file diff --git a/src/tests/TestR3linalg.hpp b/src/tests/TestR3linalg.hpp index 8f779dfb..e84ee175 100644 --- a/src/tests/TestR3linalg.hpp +++ b/src/tests/TestR3linalg.hpp @@ -134,6 +134,53 @@ class TestR3linalg : public CxxTest::TestSuite } + void test_eigen_solve_3x3() + { + R3::Matrix A( + 3.0, 0.2, 0.0, + 0.2, 2.0, 0.1, + 0.0, 0.1, 1.0); + R3::Vector eigenvalues; + R3::Matrix eigenvectors; + R3::eigen_solve_3x3(A, eigenvalues, eigenvectors); + + TS_ASSERT(eigenvalues[0] <= eigenvalues[1]); + TS_ASSERT(eigenvalues[1] <= eigenvalues[2]); + + EpsilonEqual eigclose(1.0e-8); + for (int i = 0; i < R3::Ndim; ++i) + { + R3::Vector v( + eigenvectors(0, i), + eigenvectors(1, i), + eigenvectors(2, i)); + R3::Vector Av = R3::mxvecproduct(A, v); + R3::Vector lv( + eigenvalues[i] * v[0], + eigenvalues[i] * v[1], + eigenvalues[i] * v[2]); + TS_ASSERT(eigclose(Av, lv)); + TS_ASSERT_DELTA(1.0, R3::norm(v), 1.0e-8); + } + + for (int i = 0; i < R3::Ndim; ++i) + { + for (int j = i + 1; j < R3::Ndim; ++j) + { + R3::Vector vi( + eigenvectors(0, i), + eigenvectors(1, i), + eigenvectors(2, i)); + R3::Vector vj( + eigenvectors(0, j), + eigenvectors(1, j), + eigenvectors(2, j)); + TS_ASSERT_DELTA(0.0, R3::dot(vi, vj), 1.0e-8); + } + } + } + + }; // class TestR3linalg // End of file From e504d159cf94c3890a6254f4339c73c56c87bbce Mon Sep 17 00:00:00 2001 From: joeseaer <2969111832@qq.com> Date: Sun, 2 Aug 2026 14:18:51 +0800 Subject: [PATCH 3/4] MAINT: rename ThreeDPDFCalculator and expand tests --- .../srreal/ObjCrystStructureAdapter.cpp | 2 + ...Calculator.cpp => ThreeDPDFCalculator.cpp} | 220 +++---- ...Calculator.hpp => ThreeDPDFCalculator.hpp} | 12 +- src/tests/TestObjCrystStructureAdapter.hpp | 15 + src/tests/TestPDF3DCalculator.hpp | 180 ------ src/tests/TestThreeDPDFCalculator.hpp | 541 ++++++++++++++++++ 6 files changed, 674 insertions(+), 296 deletions(-) rename src/diffpy/srreal/{PDF3DCalculator.cpp => ThreeDPDFCalculator.cpp} (83%) rename src/diffpy/srreal/{PDF3DCalculator.hpp => ThreeDPDFCalculator.hpp} (96%) delete mode 100644 src/tests/TestPDF3DCalculator.hpp create mode 100644 src/tests/TestThreeDPDFCalculator.hpp diff --git a/src/diffpy/srreal/ObjCrystStructureAdapter.cpp b/src/diffpy/srreal/ObjCrystStructureAdapter.cpp index 969f5715..8190e816 100644 --- a/src/diffpy/srreal/ObjCrystStructureAdapter.cpp +++ b/src/diffpy/srreal/ObjCrystStructureAdapter.cpp @@ -86,6 +86,8 @@ fetchSymmetryOperations(const ObjCryst::SpaceGroup& spacegroup) assert(nbtran * last <= nbsym); for (int nt = 0; nt < nbtran; ++nt) { + // Keep this compatible with ObjCryst builds where REAL is either + // float or double. const auto* pt = sgtrans[nt].tr; R3::Vector sgt(pt[0], pt[1], pt[2]); for (int i = 0; i < last; ++i) diff --git a/src/diffpy/srreal/PDF3DCalculator.cpp b/src/diffpy/srreal/ThreeDPDFCalculator.cpp similarity index 83% rename from src/diffpy/srreal/PDF3DCalculator.cpp rename to src/diffpy/srreal/ThreeDPDFCalculator.cpp index e3bbca03..28230c7c 100644 --- a/src/diffpy/srreal/PDF3DCalculator.cpp +++ b/src/diffpy/srreal/ThreeDPDFCalculator.cpp @@ -1,4 +1,4 @@ -#include +#include #include #include #include @@ -21,12 +21,12 @@ namespace diffpy { namespace srreal { // Static constants -constexpr double PDF3DCalculator::DEFAULT_RMAX_3D; -constexpr double PDF3DCalculator::DEFAULT_GRID_STEP; +constexpr double ThreeDPDFCalculator::DEFAULT_RMAX_3D; +constexpr double ThreeDPDFCalculator::DEFAULT_GRID_STEP; // Constructor --------------------------------------------------------------- -PDF3DCalculator::PDF3DCalculator() : +ThreeDPDFCalculator::ThreeDPDFCalculator() : mdr(DEFAULT_GRID_STEP), mnbins(0), maccumblocksize(32), @@ -51,59 +51,59 @@ PDF3DCalculator::PDF3DCalculator() : { this->registerDoubleAttribute("enable_nn_delta3d", this, - &PDF3DCalculator::getEnableNNDelta3DAttr, - &PDF3DCalculator::setEnableNNDelta3DAttr); + &ThreeDPDFCalculator::getEnableNNDelta3DAttr, + &ThreeDPDFCalculator::setEnableNNDelta3DAttr); this->registerDoubleAttribute("nn_delta3d", this, - &PDF3DCalculator::getNNDelta3D, - &PDF3DCalculator::setNNDelta3D); + &ThreeDPDFCalculator::getNNDelta3D, + &ThreeDPDFCalculator::setNNDelta3D); this->registerDoubleAttribute("nn_delta3d_upper_bound", this, - &PDF3DCalculator::getNNDelta3DUpperBound); + &ThreeDPDFCalculator::getNNDelta3DUpperBound); this->registerDoubleAttribute("nn_delta_positive_eta3d", this, - &PDF3DCalculator::getNNDeltaPositiveEta3D, - &PDF3DCalculator::setNNDeltaPositiveEta3D); + &ThreeDPDFCalculator::getNNDeltaPositiveEta3D, + &ThreeDPDFCalculator::setNNDeltaPositiveEta3D); this->registerDoubleAttribute("delta1_3d", this, - &PDF3DCalculator::getDistanceDelta1_3D, - &PDF3DCalculator::setDistanceDelta1_3D); + &ThreeDPDFCalculator::getDistanceDelta1_3D, + &ThreeDPDFCalculator::setDistanceDelta1_3D); this->registerDoubleAttribute("delta2_3d", this, - &PDF3DCalculator::getDistanceDelta2_3D, - &PDF3DCalculator::setDistanceDelta2_3D); + &ThreeDPDFCalculator::getDistanceDelta2_3D, + &ThreeDPDFCalculator::setDistanceDelta2_3D); this->registerDoubleAttribute("delta_shell_index3d", this, - &PDF3DCalculator::getDeltaShellIndex3DAttr, - &PDF3DCalculator::setDeltaShellIndex3DAttr); + &ThreeDPDFCalculator::getDeltaShellIndex3DAttr, + &ThreeDPDFCalculator::setDeltaShellIndex3DAttr); this->registerDoubleAttribute("delta_shell_tolerance3d", this, - &PDF3DCalculator::getDeltaShellTolerance3D, - &PDF3DCalculator::setDeltaShellTolerance3D); + &ThreeDPDFCalculator::getDeltaShellTolerance3D, + &ThreeDPDFCalculator::setDeltaShellTolerance3D); this->registerDoubleAttribute("delta_key_tolerance3d", this, - &PDF3DCalculator::getDeltaKeyTolerance3D, - &PDF3DCalculator::setDeltaKeyTolerance3D); + &ThreeDPDFCalculator::getDeltaKeyTolerance3D, + &ThreeDPDFCalculator::setDeltaKeyTolerance3D); this->registerDoubleAttribute("use_adp_scale_sensitivity3d", this, - &PDF3DCalculator::getUseADPScaleSensitivity3DAttr, - &PDF3DCalculator::setUseADPScaleSensitivity3DAttr); + &ThreeDPDFCalculator::getUseADPScaleSensitivity3DAttr, + &ThreeDPDFCalculator::setUseADPScaleSensitivity3DAttr); this->registerDoubleAttribute("adp_scale3d", this, - &PDF3DCalculator::getADPScale3D, - &PDF3DCalculator::setADPScale3D); + &ThreeDPDFCalculator::getADPScale3D, + &ThreeDPDFCalculator::setADPScale3D); this->registerDoubleAttribute("rho0_background_scale3d", this, - &PDF3DCalculator::getRho0BackgroundScale3D, - &PDF3DCalculator::setRho0BackgroundScale3D); + &ThreeDPDFCalculator::getRho0BackgroundScale3D, + &ThreeDPDFCalculator::setRho0BackgroundScale3D); this->registerDoubleAttribute("calculation_mode3d", this, - &PDF3DCalculator::getCalculationMode3DAttr, - &PDF3DCalculator::setCalculationMode3DAttr); + &ThreeDPDFCalculator::getCalculationMode3DAttr, + &ThreeDPDFCalculator::setCalculationMode3DAttr); this->registerDoubleAttribute("histogram_weight_mode3d", this, - &PDF3DCalculator::getHistogramWeightMode3DAttr, - &PDF3DCalculator::setHistogramWeightMode3DAttr); + &ThreeDPDFCalculator::getHistogramWeightMode3DAttr, + &ThreeDPDFCalculator::setHistogramWeightMode3DAttr); this->setRmax(DEFAULT_RMAX_3D); this->setRstep(mdr); @@ -112,7 +112,7 @@ PDF3DCalculator::PDF3DCalculator() : // Public Methods ------------------------------------------------------------ -void PDF3DCalculator::setGridStep(double dr) +void ThreeDPDFCalculator::setGridStep(double dr) { if (dr <= 0) throw std::invalid_argument("Grid step must be positive."); if (dr != mdr) @@ -123,43 +123,43 @@ void PDF3DCalculator::setGridStep(double dr) } } -double PDF3DCalculator::getGridStep() const +double ThreeDPDFCalculator::getGridStep() const { return mdr; } -void PDF3DCalculator::setAccumBlockSize(int bs) +void ThreeDPDFCalculator::setAccumBlockSize(int bs) { if (bs <= 0) throw std::invalid_argument("Accumulation block size must be positive."); maccumblocksize = bs; } -int PDF3DCalculator::getAccumBlockSize() const +int ThreeDPDFCalculator::getAccumBlockSize() const { return maccumblocksize; } -void PDF3DCalculator::setApplyRho0Background3D(bool v) +void ThreeDPDFCalculator::setApplyRho0Background3D(bool v) { mapplyrho0background3d = v; } -bool PDF3DCalculator::getApplyRho0Background3D() const +bool ThreeDPDFCalculator::getApplyRho0Background3D() const { return mapplyrho0background3d; } -void PDF3DCalculator::setUseCQWindow3D(bool v) +void ThreeDPDFCalculator::setUseCQWindow3D(bool v) { musecqwindow3d = v; } -bool PDF3DCalculator::getUseCQWindow3D() const +bool ThreeDPDFCalculator::getUseCQWindow3D() const { return musecqwindow3d; } -void PDF3DCalculator::setCalculationMode3D(int mode) +void ThreeDPDFCalculator::setCalculationMode3D(int mode) { if (mode != 0 && mode != 1) throw std::invalid_argument("3D calculation mode must be 0 (ADP PDF) or 1 (vector histogram)."); @@ -167,17 +167,17 @@ void PDF3DCalculator::setCalculationMode3D(int mode) mcalculationmode3d = mode; } -int PDF3DCalculator::getCalculationMode3D() const +int ThreeDPDFCalculator::getCalculationMode3D() const { return mcalculationmode3d; } -double PDF3DCalculator::getCalculationMode3DAttr() const +double ThreeDPDFCalculator::getCalculationMode3DAttr() const { return static_cast(mcalculationmode3d); } -void PDF3DCalculator::setCalculationMode3DAttr(double v) +void ThreeDPDFCalculator::setCalculationMode3DAttr(double v) { const int mode = static_cast(std::lround(v)); if (std::fabs(v - static_cast(mode)) > 1.0e-8) @@ -185,7 +185,7 @@ void PDF3DCalculator::setCalculationMode3DAttr(double v) this->setCalculationMode3D(mode); } -void PDF3DCalculator::setHistogramWeightMode3D(int mode) +void ThreeDPDFCalculator::setHistogramWeightMode3D(int mode) { if (mode != 0 && mode != 1) throw std::invalid_argument("3D histogram weight mode must be 0 (scattering) or 1 (count)."); @@ -193,17 +193,17 @@ void PDF3DCalculator::setHistogramWeightMode3D(int mode) mhistogramweightmode3d = mode; } -int PDF3DCalculator::getHistogramWeightMode3D() const +int ThreeDPDFCalculator::getHistogramWeightMode3D() const { return mhistogramweightmode3d; } -double PDF3DCalculator::getHistogramWeightMode3DAttr() const +double ThreeDPDFCalculator::getHistogramWeightMode3DAttr() const { return static_cast(mhistogramweightmode3d); } -void PDF3DCalculator::setHistogramWeightMode3DAttr(double v) +void ThreeDPDFCalculator::setHistogramWeightMode3DAttr(double v) { const int mode = static_cast(std::lround(v)); if (std::fabs(v - static_cast(mode)) > 1.0e-8) @@ -211,28 +211,28 @@ void PDF3DCalculator::setHistogramWeightMode3DAttr(double v) this->setHistogramWeightMode3D(mode); } -void PDF3DCalculator::setEnableNNDelta3D(bool v) +void ThreeDPDFCalculator::setEnableNNDelta3D(bool v) { if (menablennDelta3d != v) mticker.click(); menablennDelta3d = v; } -bool PDF3DCalculator::getEnableNNDelta3D() const +bool ThreeDPDFCalculator::getEnableNNDelta3D() const { return menablennDelta3d; } -double PDF3DCalculator::getEnableNNDelta3DAttr() const +double ThreeDPDFCalculator::getEnableNNDelta3DAttr() const { return menablennDelta3d ? 1.0 : 0.0; } -void PDF3DCalculator::setEnableNNDelta3DAttr(double v) +void ThreeDPDFCalculator::setEnableNNDelta3DAttr(double v) { this->setEnableNNDelta3D(v != 0.0); } -void PDF3DCalculator::setNNDelta3D(double v) +void ThreeDPDFCalculator::setNNDelta3D(double v) { if (v < 0.0) throw std::invalid_argument("3D NN delta must be non-negative."); const double bound = this->currentNNDeltaUpperBound(); @@ -247,22 +247,22 @@ void PDF3DCalculator::setNNDelta3D(double v) mnnDelta3d = v; } -const double& PDF3DCalculator::getNNDelta3D() const +const double& ThreeDPDFCalculator::getNNDelta3D() const { return mnnDelta3d; } -double PDF3DCalculator::getNNDelta3DUpperBound() const +double ThreeDPDFCalculator::getNNDelta3DUpperBound() const { return this->currentNNDeltaUpperBound(); } -const double& PDF3DCalculator::getNNDeltaPositiveEta3D() const +const double& ThreeDPDFCalculator::getNNDeltaPositiveEta3D() const { return mnnDeltaPositiveEta3d; } -void PDF3DCalculator::setNNDeltaPositiveEta3D(double v) +void ThreeDPDFCalculator::setNNDeltaPositiveEta3D(double v) { if (v <= 0.0 || v >= 1.0) throw std::invalid_argument("Positive-definite eta must be between 0 and 1."); @@ -275,7 +275,7 @@ void PDF3DCalculator::setNNDeltaPositiveEta3D(double v) this->validateDistanceDecayDelta3D(); } -void PDF3DCalculator::setDistanceDelta1_3D(double v) +void ThreeDPDFCalculator::setDistanceDelta1_3D(double v) { if (v < 0.0) throw std::invalid_argument("3D distance-delta delta1 must be non-negative."); if (mdistanceDelta1_3d != v) mticker.click(); @@ -283,12 +283,12 @@ void PDF3DCalculator::setDistanceDelta1_3D(double v) this->validateDistanceDecayDelta3D(); } -const double& PDF3DCalculator::getDistanceDelta1_3D() const +const double& ThreeDPDFCalculator::getDistanceDelta1_3D() const { return mdistanceDelta1_3d; } -void PDF3DCalculator::setDistanceDelta2_3D(double v) +void ThreeDPDFCalculator::setDistanceDelta2_3D(double v) { if (v < 0.0) throw std::invalid_argument("3D distance-delta delta2 must be non-negative."); if (mdistanceDelta2_3d != v) mticker.click(); @@ -296,12 +296,12 @@ void PDF3DCalculator::setDistanceDelta2_3D(double v) this->validateDistanceDecayDelta3D(); } -const double& PDF3DCalculator::getDistanceDelta2_3D() const +const double& ThreeDPDFCalculator::getDistanceDelta2_3D() const { return mdistanceDelta2_3d; } -void PDF3DCalculator::setDeltaPairTypes3D(const std::string& a, const std::string& b) +void ThreeDPDFCalculator::setDeltaPairTypes3D(const std::string& a, const std::string& b) { if (a.empty() || b.empty()) throw std::invalid_argument("Delta pair atom types must be non-empty."); @@ -310,90 +310,90 @@ void PDF3DCalculator::setDeltaPairTypes3D(const std::string& a, const std::strin mdeltaPairB3d = b; } -const std::string& PDF3DCalculator::getDeltaPairA3D() const +const std::string& ThreeDPDFCalculator::getDeltaPairA3D() const { return mdeltaPairA3d; } -const std::string& PDF3DCalculator::getDeltaPairB3D() const +const std::string& ThreeDPDFCalculator::getDeltaPairB3D() const { return mdeltaPairB3d; } -void PDF3DCalculator::setDeltaShellIndex3D(int v) +void ThreeDPDFCalculator::setDeltaShellIndex3D(int v) { if (v <= 0) throw std::invalid_argument("Delta shell index must be positive."); if (mdeltaShellIndex3d != v) mticker.click(); mdeltaShellIndex3d = v; } -int PDF3DCalculator::getDeltaShellIndex3D() const +int ThreeDPDFCalculator::getDeltaShellIndex3D() const { return mdeltaShellIndex3d; } -double PDF3DCalculator::getDeltaShellIndex3DAttr() const +double ThreeDPDFCalculator::getDeltaShellIndex3DAttr() const { return static_cast(mdeltaShellIndex3d); } -void PDF3DCalculator::setDeltaShellIndex3DAttr(double v) +void ThreeDPDFCalculator::setDeltaShellIndex3DAttr(double v) { this->setDeltaShellIndex3D(static_cast(std::floor(v + 0.5))); } -const double& PDF3DCalculator::getDeltaShellTolerance3D() const +const double& ThreeDPDFCalculator::getDeltaShellTolerance3D() const { return mdeltaShellTolerance3d; } -void PDF3DCalculator::setDeltaShellTolerance3D(double v) +void ThreeDPDFCalculator::setDeltaShellTolerance3D(double v) { if (v <= 0.0) throw std::invalid_argument("Delta shell tolerance must be positive."); if (mdeltaShellTolerance3d != v) mticker.click(); mdeltaShellTolerance3d = v; } -const double& PDF3DCalculator::getDeltaKeyTolerance3D() const +const double& ThreeDPDFCalculator::getDeltaKeyTolerance3D() const { return mdeltaKeyTolerance3d; } -void PDF3DCalculator::setDeltaKeyTolerance3D(double v) +void ThreeDPDFCalculator::setDeltaKeyTolerance3D(double v) { if (v <= 0.0) throw std::invalid_argument("Delta key tolerance must be positive."); if (mdeltaKeyTolerance3d != v) mticker.click(); mdeltaKeyTolerance3d = v; } -int PDF3DCalculator::getDeltaEligiblePairCount3D() const +int ThreeDPDFCalculator::getDeltaEligiblePairCount3D() const { return static_cast(mdeltaEligibleBondKeys.size()); } -void PDF3DCalculator::setUseADPScaleSensitivity3D(bool v) +void ThreeDPDFCalculator::setUseADPScaleSensitivity3D(bool v) { if (museadpscaleSensitivity3d != v) mticker.click(); museadpscaleSensitivity3d = v; this->validateNNDelta3D(); } -bool PDF3DCalculator::getUseADPScaleSensitivity3D() const +bool ThreeDPDFCalculator::getUseADPScaleSensitivity3D() const { return museadpscaleSensitivity3d; } -double PDF3DCalculator::getUseADPScaleSensitivity3DAttr() const +double ThreeDPDFCalculator::getUseADPScaleSensitivity3DAttr() const { return museadpscaleSensitivity3d ? 1.0 : 0.0; } -void PDF3DCalculator::setUseADPScaleSensitivity3DAttr(double v) +void ThreeDPDFCalculator::setUseADPScaleSensitivity3DAttr(double v) { this->setUseADPScaleSensitivity3D(v != 0.0); } -void PDF3DCalculator::setADPScale3D(double v) +void ThreeDPDFCalculator::setADPScale3D(double v) { if (v <= 0.0) throw std::invalid_argument("3D ADP scale must be positive."); if (madpScale3d != v) mticker.click(); @@ -401,24 +401,24 @@ void PDF3DCalculator::setADPScale3D(double v) this->validateNNDelta3D(); } -const double& PDF3DCalculator::getADPScale3D() const +const double& ThreeDPDFCalculator::getADPScale3D() const { return madpScale3d; } -void PDF3DCalculator::setRho0BackgroundScale3D(double v) +void ThreeDPDFCalculator::setRho0BackgroundScale3D(double v) { if (v < 0.0) throw std::invalid_argument("3D rho0 background scale must be non-negative."); if (mrho0backgroundscale3d != v) mticker.click(); mrho0backgroundscale3d = v; } -const double& PDF3DCalculator::getRho0BackgroundScale3D() const +const double& ThreeDPDFCalculator::getRho0BackgroundScale3D() const { return mrho0backgroundscale3d; } -QuantityType PDF3DCalculator::get3DPDF() const +QuantityType ThreeDPDFCalculator::getThreeDPDF() const { QuantityType result; std::vector grid = mgrid3d; @@ -449,7 +449,7 @@ QuantityType PDF3DCalculator::get3DPDF() const return result; } -QuantityType PDF3DCalculator::getRadialHistogram3D() const +QuantityType ThreeDPDFCalculator::getRadialHistogram3D() const { QuantityType result; result.reserve(mradialhistogram3d.size() * 2); @@ -461,7 +461,7 @@ QuantityType PDF3DCalculator::getRadialHistogram3D() const return result; } -void PDF3DCalculator::exportGrid3DBinary(const std::string& path, bool usefloat32, bool applypost) const +void ThreeDPDFCalculator::exportGrid3DBinary(const std::string& path, bool usefloat32, bool applypost) const { std::vector grid = mgrid3d; if (applypost) @@ -494,7 +494,7 @@ void PDF3DCalculator::exportGrid3DBinary(const std::string& path, bool usefloat3 if (!ofs) throw std::runtime_error("Failed while writing output file: " + path); } -void PDF3DCalculator::applyPostProcessing3D(std::vector& grid) const +void ThreeDPDFCalculator::applyPostProcessing3D(std::vector& grid) const { if (musecqwindow3d) applyQWindow3D(grid); @@ -535,7 +535,7 @@ void PDF3DCalculator::applyPostProcessing3D(std::vector& grid) const } } -double PDF3DCalculator::computeRho0Background() const +double ThreeDPDFCalculator::computeRho0Background() const { const StructureAdapterPtr& structure = this->getStructure(); if (!structure) return 0.0; @@ -543,7 +543,7 @@ double PDF3DCalculator::computeRho0Background() const return mrho0backgroundscale3d * partialpdfscale * structure->numberDensity(); } -void PDF3DCalculator::applyQWindow3D(std::vector& grid) const +void ThreeDPDFCalculator::applyQWindow3D(std::vector& grid) const { if (grid.empty()) return; @@ -676,7 +676,7 @@ void PDF3DCalculator::applyQWindow3D(std::vector& grid) const // Protected Methods --------------------------------------------------------- -void PDF3DCalculator::resetValue() +void ThreeDPDFCalculator::resetValue() { // Ensure odd number of bins so origin is centered. mnbins = static_cast(2 * ceil(this->getRmax() / mdr)) + 1; @@ -695,7 +695,7 @@ void PDF3DCalculator::resetValue() if (mevaluator) mevaluator->setFlag(USEFULLSUM, true); } -void PDF3DCalculator::addPairContribution(const BaseBondGenerator& bnds, int summationscale) +void ThreeDPDFCalculator::addPairContribution(const BaseBondGenerator& bnds, int summationscale) { if (bnds.distance() == 0.0) return; @@ -724,7 +724,7 @@ void PDF3DCalculator::addPairContribution(const BaseBondGenerator& bnds, int sum addAnisotropicGaussianToGrid(rvec, Sigma, sfprod); } -void PDF3DCalculator::addVectorHistogramToGrid(const R3::Vector& r_ij, double distance, double weight) +void ThreeDPDFCalculator::addVectorHistogramToGrid(const R3::Vector& r_ij, double distance, double weight) { if (weight == 0.0) return; @@ -750,7 +750,7 @@ void PDF3DCalculator::addVectorHistogramToGrid(const R3::Vector& r_ij, double di } } -void PDF3DCalculator::addAnisotropicGaussianToGrid(const R3::Vector& r_ij, const R3::Matrix& Sigma, double sfprod) +void ThreeDPDFCalculator::addAnisotropicGaussianToGrid(const R3::Vector& r_ij, const R3::Matrix& Sigma, double sfprod) { // Eigenvalue decomposition R3::Vector eigenvalues; @@ -873,7 +873,7 @@ void PDF3DCalculator::addAnisotropicGaussianToGrid(const R3::Vector& r_ij, const if (!any) continue; #ifdef _OPENMP -#pragma omp critical(pdf3d_tile_reduce) +#pragma omp critical(threedpdf_tile_reduce) #endif { for (int lz = 0; lz < tz; ++lz) @@ -902,13 +902,13 @@ void PDF3DCalculator::addAnisotropicGaussianToGrid(const R3::Vector& r_ij, const // Private Helpers ----------------------------------------------------------- -bool PDF3DCalculator::DeltaBondKey::operator==(const DeltaBondKey& other) const +bool ThreeDPDFCalculator::DeltaBondKey::operator==(const DeltaBondKey& other) const { return site0 == other.site0 && site1 == other.site1 && rx == other.rx && ry == other.ry && rz == other.rz; } -size_t PDF3DCalculator::DeltaBondKeyHash::operator()(const DeltaBondKey& key) const +size_t ThreeDPDFCalculator::DeltaBondKeyHash::operator()(const DeltaBondKey& key) const { size_t seed = 0; seed ^= std::hash()(key.site0) + 0x9e3779b9 + (seed << 6) + (seed >> 2); @@ -919,7 +919,7 @@ size_t PDF3DCalculator::DeltaBondKeyHash::operator()(const DeltaBondKey& key) co return seed; } -bool PDF3DCalculator::matchesDeltaPairTypes( +bool ThreeDPDFCalculator::matchesDeltaPairTypes( const std::string& atom0, const std::string& atom1) const { const bool allpairs = (mdeltaPairA3d == "*" && mdeltaPairB3d == "*"); @@ -935,8 +935,8 @@ bool PDF3DCalculator::matchesDeltaPairTypes( (atom0 == mdeltaPairB3d && atom1 == mdeltaPairA3d); } -PDF3DCalculator::DeltaBondKey -PDF3DCalculator::makeDeltaBondKey(const BaseBondGenerator& bnds) const +ThreeDPDFCalculator::DeltaBondKey +ThreeDPDFCalculator::makeDeltaBondKey(const BaseBondGenerator& bnds) const { const R3::Vector& r = bnds.r01(); DeltaBondKey key; @@ -948,14 +948,14 @@ PDF3DCalculator::makeDeltaBondKey(const BaseBondGenerator& bnds) const return key; } -bool PDF3DCalculator::isDeltaEligiblePair(const BaseBondGenerator& bnds) const +bool ThreeDPDFCalculator::isDeltaEligiblePair(const BaseBondGenerator& bnds) const { if (!menablennDelta3d || mnnDelta3d == 0.0) return false; return mdeltaEligibleBondKeys.find(this->makeDeltaBondKey(bnds)) != mdeltaEligibleBondKeys.end(); } -double PDF3DCalculator::projectedSigmaAlongBond( +double ThreeDPDFCalculator::projectedSigmaAlongBond( const R3::Matrix& Sigma, const R3::Vector& r, double distance) const @@ -966,7 +966,7 @@ double PDF3DCalculator::projectedSigmaAlongBond( return R3::dot(e, tmp); } -double PDF3DCalculator::positiveDeltaBoundAlongBond( +double ThreeDPDFCalculator::positiveDeltaBoundAlongBond( const R3::Matrix& Sigma, const R3::Vector& r, double distance) const @@ -993,32 +993,32 @@ double PDF3DCalculator::positiveDeltaBoundAlongBond( } } -bool PDF3DCalculator::isDistanceDecayDeltaActive() const +bool ThreeDPDFCalculator::isDistanceDecayDeltaActive() const { return menablennDelta3d && (mdistanceDelta1_3d > 0.0 || mdistanceDelta2_3d > 0.0); } -double PDF3DCalculator::distanceDecayDeltaFraction(double distance) const +double ThreeDPDFCalculator::distanceDecayDeltaFraction(double distance) const { if (distance <= 0.0) return 0.0; return mdistanceDelta1_3d / distance + mdistanceDelta2_3d / (distance * distance); } -double PDF3DCalculator::currentNNDeltaUpperBound() const +double ThreeDPDFCalculator::currentNNDeltaUpperBound() const { if (mnnDelta3dUpperBound <= 0.0) return 0.0; return museadpscaleSensitivity3d ? (mnnDelta3dUpperBound * madpScale3d) : mnnDelta3dUpperBound; } -bool PDF3DCalculator::isVectorHistogramMode3D() const +bool ThreeDPDFCalculator::isVectorHistogramMode3D() const { return mcalculationmode3d == 1; } -void PDF3DCalculator::validateNNDelta3D() const +void ThreeDPDFCalculator::validateNNDelta3D() const { if (this->isVectorHistogramMode3D()) return; if (mnnDelta3d < 0.0) @@ -1031,7 +1031,7 @@ void PDF3DCalculator::validateNNDelta3D() const throw std::invalid_argument("3D NN delta exceeds the positive-definite shell bound."); } -void PDF3DCalculator::validateDistanceDecayDelta3D() const +void ThreeDPDFCalculator::validateDistanceDecayDelta3D() const { if (this->isVectorHistogramMode3D()) return; if (mdistanceDelta1_3d < 0.0 || mdistanceDelta2_3d < 0.0) @@ -1052,7 +1052,7 @@ void PDF3DCalculator::validateDistanceDecayDelta3D() const } } -void PDF3DCalculator::buildDeltaEligibleShellTable() +void ThreeDPDFCalculator::buildDeltaEligibleShellTable() { mdeltaEligibleBondKeys.clear(); mdeltaShellRecords.clear(); @@ -1146,7 +1146,7 @@ void PDF3DCalculator::buildDeltaEligibleShellTable() } } -R3::Matrix PDF3DCalculator::effectivePairCovariance( +R3::Matrix ThreeDPDFCalculator::effectivePairCovariance( const BaseBondGenerator& bnds, const R3::Matrix& U_i, const R3::Matrix& U_j) const @@ -1202,7 +1202,7 @@ R3::Matrix PDF3DCalculator::effectivePairCovariance( return Sigma; } -size_t PDF3DCalculator::coordToIndex(double x, double y, double z) const +size_t ThreeDPDFCalculator::coordToIndex(double x, double y, double z) const { double halfspan = (mnbins / 2) * mdr; if (fabs(x) > halfspan || fabs(y) > halfspan || fabs(z) > halfspan) @@ -1219,7 +1219,7 @@ size_t PDF3DCalculator::coordToIndex(double x, double y, double z) const return static_cast((iz * mnbins + iy) * mnbins + ix); } -void PDF3DCalculator::indexToCoord(size_t idx, double& x, double& y, double& z) const +void ThreeDPDFCalculator::indexToCoord(size_t idx, double& x, double& y, double& z) const { size_t iz = idx / (static_cast(mnbins) * mnbins); size_t rem = idx % (static_cast(mnbins) * mnbins); @@ -1236,4 +1236,4 @@ void PDF3DCalculator::indexToCoord(size_t idx, double& x, double& y, double& z) } // namespace diffpy #include -DIFFPY_INSTANTIATE_SERIALIZATION(diffpy::srreal::PDF3DCalculator) +DIFFPY_INSTANTIATE_SERIALIZATION(diffpy::srreal::ThreeDPDFCalculator) diff --git a/src/diffpy/srreal/PDF3DCalculator.hpp b/src/diffpy/srreal/ThreeDPDFCalculator.hpp similarity index 96% rename from src/diffpy/srreal/PDF3DCalculator.hpp rename to src/diffpy/srreal/ThreeDPDFCalculator.hpp index b14abc97..425fcf02 100644 --- a/src/diffpy/srreal/PDF3DCalculator.hpp +++ b/src/diffpy/srreal/ThreeDPDFCalculator.hpp @@ -1,5 +1,5 @@ -#ifndef PDF3DCALCULATOR_HPP_INCLUDED -#define PDF3DCALCULATOR_HPP_INCLUDED +#ifndef THREEDPDFCALCULATOR_HPP_INCLUDED +#define THREEDPDFCALCULATOR_HPP_INCLUDED #include #include @@ -13,15 +13,15 @@ namespace diffpy { namespace srreal { -class PDF3DCalculator : public PDFCalculator +class ThreeDPDFCalculator : public PDFCalculator { public: // constructor - PDF3DCalculator(); + ThreeDPDFCalculator(); // Public interface to retrieve 3D PDF data // Returns a flat vector of [x0, y0, z0, G0, x1, y1, z1, G1, ...] - QuantityType get3DPDF() const; + QuantityType getThreeDPDF() const; // Export dense 3D grid (nz, ny, nx) directly to binary file. // usefloat32=true writes float32, otherwise float64. @@ -230,4 +230,4 @@ class PDF3DCalculator : public PDFCalculator } // namespace srreal } // namespace diffpy -#endif // PDF3DCALCULATOR_HPP_INCLUDED +#endif // THREEDPDFCALCULATOR_HPP_INCLUDED diff --git a/src/tests/TestObjCrystStructureAdapter.hpp b/src/tests/TestObjCrystStructureAdapter.hpp index cd4bf902..1cacfbc9 100644 --- a/src/tests/TestObjCrystStructureAdapter.hpp +++ b/src/tests/TestObjCrystStructureAdapter.hpp @@ -144,6 +144,21 @@ class TestObjCrystStructureAdapter : public CxxTest::TestSuite } + void test_real_translation_vectors() + { + // Ni.cif is F-centered. Its four equivalent atoms require + // fetchSymmetryOperations to read ObjCryst translation vectors, + // including builds where ObjCryst::REAL is double. + CrystalStructureAdapterPtr ni = + boost::dynamic_pointer_cast(m_ni); + TS_ASSERT(ni); + TS_ASSERT_EQUALS(4, ni->siteMultiplicity(0)); + CrystalStructureAdapter::AtomVector equivalent = + ni->getEquivalentAtoms(0); + TS_ASSERT_EQUALS(4u, equivalent.size()); + } + + void test_numberDensity() { const double eps = 1.0e-7; diff --git a/src/tests/TestPDF3DCalculator.hpp b/src/tests/TestPDF3DCalculator.hpp deleted file mode 100644 index dfc58d7e..00000000 --- a/src/tests/TestPDF3DCalculator.hpp +++ /dev/null @@ -1,180 +0,0 @@ -/***************************************************************************** -* -* class TestPDF3DCalculator -- unit tests for PDF3DCalculator class -* -*****************************************************************************/ - -#include - -#include -#include -#include -#include - -#include -#include -#include - -using namespace std; -using namespace diffpy::srreal; - -class TestPDF3DCalculator : public CxxTest::TestSuite -{ - private: - - boost::shared_ptr mpdfc; - AtomicStructureAdapterPtr mstru2; - double meps; - - public: - - void setUp() - { - meps = diffpy::mathutils::SQRT_DOUBLE_EPS; - mpdfc.reset(new PDF3DCalculator); - - mstru2 = boost::make_shared(); - Atom ai; - ai.atomtype = "Ni"; - ai.xyz_cartn[0] = 0.0; - ai.xyz_cartn[1] = 0.0; - ai.xyz_cartn[2] = 0.0; - ai.uij_cartn = R3::identity(); - ai.uij_cartn(0, 0) = ai.uij_cartn(1, 1) = - ai.uij_cartn(2, 2) = 0.004; - mstru2->append(ai); - ai.xyz_cartn[0] = 1.0; - mstru2->append(ai); - } - - - void test_defaults_and_attributes() - { - TS_ASSERT_DELTA(0.1, mpdfc->getGridStep(), meps); - TS_ASSERT_EQUALS(32, mpdfc->getAccumBlockSize()); - TS_ASSERT_EQUALS(true, mpdfc->getApplyRho0Background3D()); - TS_ASSERT_EQUALS(true, mpdfc->getUseCQWindow3D()); - TS_ASSERT_EQUALS(false, mpdfc->getEnableNNDelta3D()); - TS_ASSERT_EQUALS(0, mpdfc->getCalculationMode3D()); - TS_ASSERT_EQUALS(0, mpdfc->getHistogramWeightMode3D()); - - TS_ASSERT_DELTA(0.0, - mpdfc->getDoubleAttr("enable_nn_delta3d"), meps); - TS_ASSERT_DELTA(0.0, mpdfc->getDoubleAttr("nn_delta3d"), meps); - TS_ASSERT_DELTA(0.9, - mpdfc->getDoubleAttr("nn_delta_positive_eta3d"), meps); - TS_ASSERT_DELTA(0.0, mpdfc->getDoubleAttr("delta1_3d"), meps); - TS_ASSERT_DELTA(0.0, mpdfc->getDoubleAttr("delta2_3d"), meps); - TS_ASSERT_DELTA(1.0, - mpdfc->getDoubleAttr("delta_shell_index3d"), meps); - TS_ASSERT_DELTA(1.0, mpdfc->getDoubleAttr("adp_scale3d"), meps); - - mpdfc->setDoubleAttr("calculation_mode3d", 1.0); - mpdfc->setDoubleAttr("histogram_weight_mode3d", 1.0); - mpdfc->setDoubleAttr("enable_nn_delta3d", 1.0); - TS_ASSERT_EQUALS(1, mpdfc->getCalculationMode3D()); - TS_ASSERT_EQUALS(1, mpdfc->getHistogramWeightMode3D()); - TS_ASSERT_EQUALS(true, mpdfc->getEnableNNDelta3D()); - } - - - void test_invalid_configuration() - { - TS_ASSERT_THROWS(mpdfc->setGridStep(0.0), invalid_argument); - TS_ASSERT_THROWS(mpdfc->setAccumBlockSize(0), invalid_argument); - TS_ASSERT_THROWS(mpdfc->setCalculationMode3D(2), invalid_argument); - TS_ASSERT_THROWS(mpdfc->setHistogramWeightMode3D(2), - invalid_argument); - TS_ASSERT_THROWS( - mpdfc->setDoubleAttr("calculation_mode3d", 0.5), - invalid_argument); - TS_ASSERT_THROWS(mpdfc->setADPScale3D(0.0), invalid_argument); - TS_ASSERT_THROWS(mpdfc->setRho0BackgroundScale3D(-1.0), - invalid_argument); - } - - - void test_vector_histogram_eval() - { - mpdfc->setRmax(2.0); - mpdfc->setGridStep(0.5); - mpdfc->setCalculationMode3D(1); - mpdfc->setHistogramWeightMode3D(1); - mpdfc->eval(mstru2); - - QuantityType pdf3d = mpdfc->get3DPDF(); - TS_ASSERT(!pdf3d.empty()); - TS_ASSERT_EQUALS(0u, pdf3d.size() % 4); - - QuantityType rh = mpdfc->getRadialHistogram3D(); - size_t nr = static_cast( - std::ceil(mpdfc->getRmax() / mpdfc->getGridStep())) + 1; - TS_ASSERT_EQUALS(2 * nr, rh.size()); - - bool found_distance_one = false; - for (size_t i = 0; i < rh.size(); i += 2) - { - if (std::fabs(rh[i] - 1.0) < meps && rh[i + 1] != 0.0) - { - found_distance_one = true; - break; - } - } - TS_ASSERT(found_distance_one); - } - - - void test_serialization() - { - mpdfc->setRmax(1.0); - mpdfc->setGridStep(0.25); - mpdfc->setAccumBlockSize(7); - mpdfc->setApplyRho0Background3D(false); - mpdfc->setUseCQWindow3D(false); - mpdfc->setCalculationMode3D(1); - mpdfc->setHistogramWeightMode3D(1); - mpdfc->setEnableNNDelta3D(true); - mpdfc->setNNDelta3D(0.02); - mpdfc->setNNDeltaPositiveEta3D(0.8); - mpdfc->setDistanceDelta1_3D(0.03); - mpdfc->setDistanceDelta2_3D(0.04); - mpdfc->setDeltaPairTypes3D("Ni", "Ni"); - mpdfc->setDeltaShellIndex3D(2); - mpdfc->setDeltaShellTolerance3D(0.06); - mpdfc->setDeltaKeyTolerance3D(1.0e-5); - mpdfc->setUseADPScaleSensitivity3D(true); - mpdfc->setADPScale3D(1.2); - mpdfc->setRho0BackgroundScale3D(0.5); - - stringstream storage(ios::in | ios::out | ios::binary); - diffpy::serialization::oarchive oa(storage, ios::binary); - oa << mpdfc; - diffpy::serialization::iarchive ia(storage, ios::binary); - boost::shared_ptr pdfc1; - ia >> pdfc1; - - TS_ASSERT_DIFFERS(pdfc1.get(), mpdfc.get()); - TS_ASSERT_DELTA(0.25, pdfc1->getGridStep(), meps); - TS_ASSERT_EQUALS(7, pdfc1->getAccumBlockSize()); - TS_ASSERT_EQUALS(false, pdfc1->getApplyRho0Background3D()); - TS_ASSERT_EQUALS(false, pdfc1->getUseCQWindow3D()); - TS_ASSERT_EQUALS(1, pdfc1->getCalculationMode3D()); - TS_ASSERT_EQUALS(1, pdfc1->getHistogramWeightMode3D()); - TS_ASSERT_EQUALS(true, pdfc1->getEnableNNDelta3D()); - TS_ASSERT_DELTA(0.02, pdfc1->getNNDelta3D(), meps); - TS_ASSERT_DELTA(0.8, pdfc1->getNNDeltaPositiveEta3D(), meps); - TS_ASSERT_DELTA(0.03, pdfc1->getDistanceDelta1_3D(), meps); - TS_ASSERT_DELTA(0.04, pdfc1->getDistanceDelta2_3D(), meps); - TS_ASSERT_EQUALS(string("Ni"), pdfc1->getDeltaPairA3D()); - TS_ASSERT_EQUALS(string("Ni"), pdfc1->getDeltaPairB3D()); - TS_ASSERT_EQUALS(2, pdfc1->getDeltaShellIndex3D()); - TS_ASSERT_DELTA(0.06, pdfc1->getDeltaShellTolerance3D(), meps); - TS_ASSERT_DELTA(1.0e-5, pdfc1->getDeltaKeyTolerance3D(), meps); - TS_ASSERT_EQUALS(true, pdfc1->getUseADPScaleSensitivity3D()); - TS_ASSERT_DELTA(1.2, pdfc1->getADPScale3D(), meps); - TS_ASSERT_DELTA(0.5, pdfc1->getRho0BackgroundScale3D(), meps); - } - -}; // class TestPDF3DCalculator - -// End of file diff --git a/src/tests/TestThreeDPDFCalculator.hpp b/src/tests/TestThreeDPDFCalculator.hpp new file mode 100644 index 00000000..d1d8e3d1 --- /dev/null +++ b/src/tests/TestThreeDPDFCalculator.hpp @@ -0,0 +1,541 @@ +/***************************************************************************** +* +* class TestThreeDPDFCalculator -- unit tests for ThreeDPDFCalculator +* +*****************************************************************************/ + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +using namespace std; +using namespace diffpy::srreal; + +// Local Helpers ------------------------------------------------------------- + +namespace { + +AtomicStructureAdapterPtr +makeDimer(const string& atomtype, double uiso) +{ + AtomicStructureAdapterPtr stru = + boost::make_shared(); + Atom atom; + atom.atomtype = atomtype; + atom.xyz_cartn = R3::Vector(0.0, 0.0, 0.0); + atom.uij_cartn = R3::zeromatrix(); + atom.uij_cartn(0, 0) = atom.uij_cartn(1, 1) = + atom.uij_cartn(2, 2) = uiso; + stru->append(atom); + atom.xyz_cartn[0] = 1.0; + stru->append(atom); + return stru; +} + + +PeriodicStructureAdapterPtr +makePeriodicDimer() +{ + PeriodicStructureAdapterPtr stru = + boost::make_shared(); + stru->setLatPar(10.0, 10.0, 10.0, 90.0, 90.0, 90.0); + AtomicStructureAdapterPtr atoms = makeDimer("Ni", 0.004); + stru->append(atoms->at(0)); + stru->append(atoms->at(1)); + return stru; +} + + +void +configureSmallStandardCalculation(ThreeDPDFCalculator& calc) +{ + calc.setRmax(2.0); + calc.setGridStep(0.25); + calc.setUseCQWindow3D(false); + calc.setApplyRho0Background3D(false); +} + + +bool +allSparseValuesFinite(const QuantityType& data) +{ + if (data.size() % 4) return false; + for (size_t i = 0; i < data.size(); ++i) + { + if (!std::isfinite(data[i])) return false; + } + return true; +} + + +bool +findSparseValue( + const QuantityType& data, + double x, + double y, + double z, + double eps, + double& value) +{ + for (size_t i = 0; i + 3 < data.size(); i += 4) + { + if (std::fabs(data[i] - x) <= eps && + std::fabs(data[i + 1] - y) <= eps && + std::fabs(data[i + 2] - z) <= eps) + { + value = data[i + 3]; + return true; + } + } + return false; +} + + +bool +quantitiesClose( + const QuantityType& a, + const QuantityType& b, + double eps) +{ + if (a.size() != b.size()) return false; + for (size_t i = 0; i < a.size(); ++i) + { + if (std::fabs(a[i] - b[i]) > eps) return false; + } + return true; +} + + +size_t +binaryFileSize(const string& path) +{ + ifstream stream(path.c_str(), ios::binary | ios::ate); + if (!stream) return 0; + return static_cast(stream.tellg()); +} + +} // namespace + +////////////////////////////////////////////////////////////////////////////// +// class TestThreeDPDFCalculator +////////////////////////////////////////////////////////////////////////////// + +class TestThreeDPDFCalculator : public CxxTest::TestSuite +{ + private: + + boost::shared_ptr mpdfc; + AtomicStructureAdapterPtr mstru2; + double meps; + const string mfloat32path = "threedpdf_test_float32.bin"; + const string mfloat64path = "threedpdf_test_float64.bin"; + + public: + + void setUp() + { + meps = diffpy::mathutils::SQRT_DOUBLE_EPS; + mpdfc.reset(new ThreeDPDFCalculator); + mstru2 = makeDimer("Ni", 0.004); + std::remove(mfloat32path.c_str()); + std::remove(mfloat64path.c_str()); + } + + + void tearDown() + { + std::remove(mfloat32path.c_str()); + std::remove(mfloat64path.c_str()); + } + + + void test_defaults_and_attributes() + { + TS_ASSERT_DELTA(0.1, mpdfc->getGridStep(), meps); + TS_ASSERT_EQUALS(32, mpdfc->getAccumBlockSize()); + TS_ASSERT_EQUALS(true, mpdfc->getApplyRho0Background3D()); + TS_ASSERT_EQUALS(true, mpdfc->getUseCQWindow3D()); + TS_ASSERT_EQUALS(false, mpdfc->getEnableNNDelta3D()); + TS_ASSERT_EQUALS(0, mpdfc->getCalculationMode3D()); + TS_ASSERT_EQUALS(0, mpdfc->getHistogramWeightMode3D()); + TS_ASSERT_DELTA(0.0, mpdfc->getNNDelta3D(), meps); + TS_ASSERT_DELTA(0.9, mpdfc->getNNDeltaPositiveEta3D(), meps); + TS_ASSERT_DELTA(0.0, mpdfc->getDistanceDelta1_3D(), meps); + TS_ASSERT_DELTA(0.0, mpdfc->getDistanceDelta2_3D(), meps); + TS_ASSERT_EQUALS(string("*"), mpdfc->getDeltaPairA3D()); + TS_ASSERT_EQUALS(string("*"), mpdfc->getDeltaPairB3D()); + TS_ASSERT_EQUALS(1, mpdfc->getDeltaShellIndex3D()); + TS_ASSERT_DELTA( + 0.05, mpdfc->getDeltaShellTolerance3D(), meps); + TS_ASSERT_DELTA( + 1.0e-6, mpdfc->getDeltaKeyTolerance3D(), meps); + TS_ASSERT_EQUALS( + false, mpdfc->getUseADPScaleSensitivity3D()); + TS_ASSERT_DELTA(1.0, mpdfc->getADPScale3D(), meps); + TS_ASSERT_DELTA( + 1.0, mpdfc->getRho0BackgroundScale3D(), meps); + + TS_ASSERT_DELTA(0.0, + mpdfc->getDoubleAttr("enable_nn_delta3d"), meps); + TS_ASSERT_DELTA(0.0, + mpdfc->getDoubleAttr("nn_delta3d"), meps); + TS_ASSERT_DELTA(0.9, + mpdfc->getDoubleAttr("nn_delta_positive_eta3d"), meps); + TS_ASSERT_DELTA(0.0, + mpdfc->getDoubleAttr("delta1_3d"), meps); + TS_ASSERT_DELTA(0.0, + mpdfc->getDoubleAttr("delta2_3d"), meps); + TS_ASSERT_DELTA(1.0, + mpdfc->getDoubleAttr("delta_shell_index3d"), meps); + TS_ASSERT_DELTA(1.0, + mpdfc->getDoubleAttr("adp_scale3d"), meps); + + mpdfc->setDoubleAttr("calculation_mode3d", 1.0); + mpdfc->setDoubleAttr("histogram_weight_mode3d", 1.0); + mpdfc->setDoubleAttr("enable_nn_delta3d", 1.0); + mpdfc->setDoubleAttr("use_adp_scale_sensitivity3d", 1.0); + TS_ASSERT_EQUALS(1, mpdfc->getCalculationMode3D()); + TS_ASSERT_EQUALS(1, mpdfc->getHistogramWeightMode3D()); + TS_ASSERT_EQUALS(true, mpdfc->getEnableNNDelta3D()); + TS_ASSERT_EQUALS( + true, mpdfc->getUseADPScaleSensitivity3D()); + } + + + void test_invalid_configuration() + { + TS_ASSERT_THROWS(mpdfc->setGridStep(0.0), invalid_argument); + TS_ASSERT_THROWS(mpdfc->setGridStep(-0.1), invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setAccumBlockSize(0), invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setCalculationMode3D(2), invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setHistogramWeightMode3D(-1), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setDoubleAttr("calculation_mode3d", 0.5), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setDoubleAttr("histogram_weight_mode3d", 0.5), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setNNDelta3D(-0.01), invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setNNDeltaPositiveEta3D(0.0), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setNNDeltaPositiveEta3D(1.0), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setDistanceDelta1_3D(-0.01), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setDistanceDelta2_3D(-0.01), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setDeltaPairTypes3D("", "Ni"), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setDeltaShellIndex3D(0), invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setDeltaShellTolerance3D(0.0), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setDeltaKeyTolerance3D(0.0), + invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setADPScale3D(0.0), invalid_argument); + TS_ASSERT_THROWS( + mpdfc->setRho0BackgroundScale3D(-1.0), + invalid_argument); + } + + + void test_standard_grid_mode() + { + configureSmallStandardCalculation(*mpdfc); + mpdfc->eval(mstru2); + + QuantityType result = mpdfc->getThreeDPDF(); + TS_ASSERT(!result.empty()); + TS_ASSERT_EQUALS(0u, result.size() % 4); + TS_ASSERT(allSparseValuesFinite(result)); + + double positive = 0.0; + double negative = 0.0; + TS_ASSERT(findSparseValue( + result, 1.0, 0.0, 0.0, meps, positive)); + TS_ASSERT(findSparseValue( + result, -1.0, 0.0, 0.0, meps, negative)); + TS_ASSERT(positive > 0.0); + TS_ASSERT(negative > 0.0); + TS_ASSERT_DELTA(positive, negative, meps); + } + + + void test_vector_histogram_and_pre_rename_regression() + { + mpdfc->setRmax(2.0); + mpdfc->setGridStep(0.5); + mpdfc->setCalculationMode3D(1); + mpdfc->setHistogramWeightMode3D(1); + mpdfc->eval(mstru2); + + QuantityType result = mpdfc->getThreeDPDF(); + TS_ASSERT_EQUALS(8u, result.size()); + double positive = 0.0; + double negative = 0.0; + TS_ASSERT(findSparseValue( + result, 1.0, 0.0, 0.0, meps, positive)); + TS_ASSERT(findSparseValue( + result, -1.0, 0.0, 0.0, meps, negative)); + TS_ASSERT_DELTA(1.0, positive, meps); + TS_ASSERT_DELTA(1.0, negative, meps); + + QuantityType radial = mpdfc->getRadialHistogram3D(); + size_t nr = static_cast( + std::ceil(mpdfc->getRmax() / + mpdfc->getGridStep())) + 1; + TS_ASSERT_EQUALS(2 * nr, radial.size()); + TS_ASSERT_DELTA(2.0, radial[5], meps); + } + + + void test_histogram_mode_does_not_use_adps() + { + AtomicStructureAdapterPtr largeadp = + makeDimer("Ni", 0.4); + mpdfc->setRmax(2.0); + mpdfc->setGridStep(0.5); + mpdfc->setCalculationMode3D(1); + mpdfc->setHistogramWeightMode3D(1); + + mpdfc->eval(mstru2); + QuantityType small = mpdfc->getThreeDPDF(); + mpdfc->eval(largeadp); + QuantityType large = mpdfc->getThreeDPDF(); + TS_ASSERT(quantitiesClose(small, large, 0.0)); + } + + + void test_q_window_changes_grid() + { + configureSmallStandardCalculation(*mpdfc); + mpdfc->eval(mstru2); + QuantityType unwindowed = mpdfc->getThreeDPDF(); + + mpdfc->setUseCQWindow3D(true); + mpdfc->setQmin(0.0); + mpdfc->setQmax(4.0); + QuantityType windowed = mpdfc->getThreeDPDF(); + + TS_ASSERT(!windowed.empty()); + TS_ASSERT(allSparseValuesFinite(windowed)); + TS_ASSERT(!quantitiesClose(unwindowed, windowed, 1.0e-10)); + } + + + void test_rho0_background() + { + PeriodicStructureAdapterPtr periodic = makePeriodicDimer(); + configureSmallStandardCalculation(*mpdfc); + mpdfc->eval(periodic); + QuantityType raw = mpdfc->getThreeDPDF(); + TS_ASSERT(!raw.empty()); + + mpdfc->setApplyRho0Background3D(true); + mpdfc->setRho0BackgroundScale3D(0.75); + QuantityType corrected = mpdfc->getThreeDPDF(); + + double correctedvalue = 0.0; + TS_ASSERT(findSparseValue( + corrected, raw[0], raw[1], raw[2], + meps, correctedvalue)); + const double expected = + 0.75 * periodic->numberDensity(); + TS_ASSERT_DELTA( + expected, raw[3] - correctedvalue, 1.0e-10); + } + + + void test_delta_pair_filter_and_shell_covariance() + { + configureSmallStandardCalculation(*mpdfc); + mpdfc->setEnableNNDelta3D(true); + mpdfc->setDeltaPairTypes3D("O", "O"); + mpdfc->eval(mstru2); + TS_ASSERT_EQUALS(0, mpdfc->getDeltaEligiblePairCount3D()); + + mpdfc->setDeltaPairTypes3D("Ni", "Ni"); + mpdfc->eval(mstru2); + TS_ASSERT( + mpdfc->getDeltaEligiblePairCount3D() > 0); + const double bound = + mpdfc->getNNDelta3DUpperBound(); + TS_ASSERT(bound > 0.0); + QuantityType baseline = mpdfc->getThreeDPDF(); + + mpdfc->setNNDelta3D(0.5 * bound); + mpdfc->eval(mstru2); + QuantityType correlated = mpdfc->getThreeDPDF(); + TS_ASSERT(!quantitiesClose( + baseline, correlated, 1.0e-10)); + TS_ASSERT_THROWS( + mpdfc->setNNDelta3D(bound), invalid_argument); + } + + + void test_distance_delta_covariance() + { + configureSmallStandardCalculation(*mpdfc); + mpdfc->setEnableNNDelta3D(true); + mpdfc->setDeltaPairTypes3D("Ni", "Ni"); + mpdfc->eval(mstru2); + QuantityType baseline = mpdfc->getThreeDPDF(); + + mpdfc->setDistanceDelta1_3D(0.1); + mpdfc->setDistanceDelta2_3D(0.05); + mpdfc->eval(mstru2); + QuantityType correlated = mpdfc->getThreeDPDF(); + TS_ASSERT(!quantitiesClose( + baseline, correlated, 1.0e-10)); + } + + + void test_rdf_normalization_cancels_atom_scattering_scale() + { + AtomicStructureAdapterPtr carbon = + makeDimer("C", 0.004); + configureSmallStandardCalculation(*mpdfc); + mpdfc->eval(mstru2); + QuantityType nickelresult = mpdfc->getThreeDPDF(); + mpdfc->eval(carbon); + QuantityType carbonresult = mpdfc->getThreeDPDF(); + + TS_ASSERT(quantitiesClose( + nickelresult, carbonresult, 1.0e-10)); + } + + + void test_binary_export() + { + mpdfc->setRmax(1.0); + mpdfc->setGridStep(0.5); + mpdfc->setUseCQWindow3D(false); + mpdfc->setApplyRho0Background3D(false); + mpdfc->eval(mstru2); + + mpdfc->exportGrid3DBinary( + mfloat32path, true, false); + mpdfc->exportGrid3DBinary( + mfloat64path, false, false); + const size_t gridcells = 5u * 5u * 5u; + TS_ASSERT_EQUALS( + gridcells * sizeof(float), + binaryFileSize(mfloat32path)); + TS_ASSERT_EQUALS( + gridcells * sizeof(double), + binaryFileSize(mfloat64path)); + TS_ASSERT_THROWS( + mpdfc->exportGrid3DBinary( + "libdiffpy_missing_test_directory/grid.bin"), + runtime_error); + } + + + void test_serialization_round_trip() + { + mpdfc->setRmax(1.0); + mpdfc->setGridStep(0.25); + mpdfc->setAccumBlockSize(7); + mpdfc->setApplyRho0Background3D(false); + mpdfc->setUseCQWindow3D(false); + mpdfc->setCalculationMode3D(1); + mpdfc->setHistogramWeightMode3D(1); + mpdfc->setEnableNNDelta3D(true); + mpdfc->setNNDelta3D(0.02); + mpdfc->setNNDeltaPositiveEta3D(0.8); + mpdfc->setDistanceDelta1_3D(0.03); + mpdfc->setDistanceDelta2_3D(0.04); + mpdfc->setDeltaPairTypes3D("Ni", "Ni"); + mpdfc->setDeltaShellIndex3D(2); + mpdfc->setDeltaShellTolerance3D(0.06); + mpdfc->setDeltaKeyTolerance3D(1.0e-5); + mpdfc->setUseADPScaleSensitivity3D(true); + mpdfc->setADPScale3D(1.2); + mpdfc->setRho0BackgroundScale3D(0.5); + + stringstream storage(ios::in | ios::out | ios::binary); + diffpy::serialization::oarchive oa( + storage, ios::binary); + oa << mpdfc; + diffpy::serialization::iarchive ia( + storage, ios::binary); + boost::shared_ptr restored; + ia >> restored; + + TS_ASSERT_DIFFERS(restored.get(), mpdfc.get()); + TS_ASSERT_DELTA( + 0.25, restored->getGridStep(), meps); + TS_ASSERT_EQUALS(7, restored->getAccumBlockSize()); + TS_ASSERT_EQUALS( + false, restored->getApplyRho0Background3D()); + TS_ASSERT_EQUALS( + false, restored->getUseCQWindow3D()); + TS_ASSERT_EQUALS( + 1, restored->getCalculationMode3D()); + TS_ASSERT_EQUALS( + 1, restored->getHistogramWeightMode3D()); + TS_ASSERT_EQUALS( + true, restored->getEnableNNDelta3D()); + TS_ASSERT_DELTA( + 0.02, restored->getNNDelta3D(), meps); + TS_ASSERT_DELTA( + 0.8, restored->getNNDeltaPositiveEta3D(), meps); + TS_ASSERT_DELTA( + 0.03, restored->getDistanceDelta1_3D(), meps); + TS_ASSERT_DELTA( + 0.04, restored->getDistanceDelta2_3D(), meps); + TS_ASSERT_EQUALS( + string("Ni"), restored->getDeltaPairA3D()); + TS_ASSERT_EQUALS( + string("Ni"), restored->getDeltaPairB3D()); + TS_ASSERT_EQUALS( + 2, restored->getDeltaShellIndex3D()); + TS_ASSERT_DELTA( + 0.06, restored->getDeltaShellTolerance3D(), meps); + TS_ASSERT_DELTA( + 1.0e-5, restored->getDeltaKeyTolerance3D(), meps); + TS_ASSERT_EQUALS( + true, restored->getUseADPScaleSensitivity3D()); + TS_ASSERT_DELTA( + 1.2, restored->getADPScale3D(), meps); + TS_ASSERT_DELTA( + 0.5, restored->getRho0BackgroundScale3D(), meps); + + restored->eval(mstru2); + QuantityType restoredresult = + restored->getThreeDPDF(); + mpdfc->eval(mstru2); + QuantityType originalresult = + mpdfc->getThreeDPDF(); + TS_ASSERT(quantitiesClose( + originalresult, restoredresult, 0.0)); + } + +}; // class TestThreeDPDFCalculator + +// End of file From c69c32f5e1d6ffad9f63f25607f65a2f9b141a11 Mon Sep 17 00:00:00 2001 From: joeseaer <2969111832@qq.com> Date: Sun, 2 Aug 2026 15:04:30 +0800 Subject: [PATCH 4/4] MAINT: adapt 3D tests to migration smart pointers --- src/tests/TestObjCrystStructureAdapter.hpp | 3 ++- src/tests/TestThreeDPDFCalculator.hpp | 10 +++++----- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/tests/TestObjCrystStructureAdapter.hpp b/src/tests/TestObjCrystStructureAdapter.hpp index 1cacfbc9..fbe8af32 100644 --- a/src/tests/TestObjCrystStructureAdapter.hpp +++ b/src/tests/TestObjCrystStructureAdapter.hpp @@ -19,6 +19,7 @@ * *****************************************************************************/ +#include #include #include @@ -150,7 +151,7 @@ class TestObjCrystStructureAdapter : public CxxTest::TestSuite // fetchSymmetryOperations to read ObjCryst translation vectors, // including builds where ObjCryst::REAL is double. CrystalStructureAdapterPtr ni = - boost::dynamic_pointer_cast(m_ni); + std::dynamic_pointer_cast(m_ni); TS_ASSERT(ni); TS_ASSERT_EQUALS(4, ni->siteMultiplicity(0)); CrystalStructureAdapter::AtomVector equivalent = diff --git a/src/tests/TestThreeDPDFCalculator.hpp b/src/tests/TestThreeDPDFCalculator.hpp index d1d8e3d1..e4836d81 100644 --- a/src/tests/TestThreeDPDFCalculator.hpp +++ b/src/tests/TestThreeDPDFCalculator.hpp @@ -10,10 +10,10 @@ #include #include #include +#include #include #include #include -#include #include #include @@ -31,7 +31,7 @@ AtomicStructureAdapterPtr makeDimer(const string& atomtype, double uiso) { AtomicStructureAdapterPtr stru = - boost::make_shared(); + std::make_shared(); Atom atom; atom.atomtype = atomtype; atom.xyz_cartn = R3::Vector(0.0, 0.0, 0.0); @@ -49,7 +49,7 @@ PeriodicStructureAdapterPtr makePeriodicDimer() { PeriodicStructureAdapterPtr stru = - boost::make_shared(); + std::make_shared(); stru->setLatPar(10.0, 10.0, 10.0, 90.0, 90.0, 90.0); AtomicStructureAdapterPtr atoms = makeDimer("Ni", 0.004); stru->append(atoms->at(0)); @@ -136,7 +136,7 @@ class TestThreeDPDFCalculator : public CxxTest::TestSuite { private: - boost::shared_ptr mpdfc; + std::shared_ptr mpdfc; AtomicStructureAdapterPtr mstru2; double meps; const string mfloat32path = "threedpdf_test_float32.bin"; @@ -484,7 +484,7 @@ class TestThreeDPDFCalculator : public CxxTest::TestSuite oa << mpdfc; diffpy::serialization::iarchive ia( storage, ios::binary); - boost::shared_ptr restored; + std::shared_ptr restored; ia >> restored; TS_ASSERT_DIFFERS(restored.get(), mpdfc.get());