From c4395ccd405280df25f0c7c1d26ea2c81d85c067 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 9 May 2025 21:22:35 -0400 Subject: [PATCH 01/15] Remove the system option for the coordinate in favor of the mapped coordinate transformation which works way better for power laws and exponentials --- exputil/OrthoFunction.cc | 40 +++++++++++++++++----------------------- include/OrthoFunction.H | 7 ++----- 2 files changed, 19 insertions(+), 28 deletions(-) diff --git a/exputil/OrthoFunction.cc b/exputil/OrthoFunction.cc index d7812a794..79f9c5951 100644 --- a/exputil/OrthoFunction.cc +++ b/exputil/OrthoFunction.cc @@ -2,8 +2,8 @@ #include OrthoFunction::OrthoFunction -(int nmax, DensFunc W, double rmin, double rmax, double scl, int dof, bool seg) : - nmax(nmax), W(W), rmin(rmin), rmax(rmax), scale(scl), dof(dof), segment(seg) +(int nmax, DensFunc W, double rmin, double rmax, double scl, int dof) : + nmax(nmax), W(W), rmin(rmin), rmax(rmax), scale(scl), dof(dof) { // Quadrature knots/weights lq = std::make_shared(knots); @@ -27,14 +27,14 @@ void OrthoFunction::generate() norm.resize(nmax+1); // Initial values - beta[0] = norm[0] = scalar_prod(0, 0); - alph[0] = scalar_prod(0, 1)/norm[0]; + beta(0) = norm(0) = scalar_prod(0, 0); + alph(0) = scalar_prod(0, 1)/norm(0); // Remaining values for (int i=1; i<=nmax; i++) { - norm[i] = scalar_prod(i, 0); - alph[i] = scalar_prod(i, 1)/norm[i]; - beta[i] = norm[i]/norm[i-1]; + norm(i) = scalar_prod(i, 0); + alph(i) = scalar_prod(i, 1)/norm[i]; + beta(i) = norm(i)/norm(i-1); } } @@ -45,10 +45,10 @@ double OrthoFunction::scalar_prod(int n, int m) for (int i=0; iknot(i); double r = x_to_r(x); - double y = 2.0*r/scale; if (segment) y = x; - Eigen::VectorXd p = poly_eval(y, n) * W(r); + Eigen::VectorXd p = poly_eval(r, n) * W(r); + ans += dx*lq->weight(i)*d_r_to_x(x)*pow(r, dof-1) * - p(n) * p(n) * pow(y, m); + p(n) * p(n) * pow(r, m); } return ans; @@ -81,15 +81,14 @@ Eigen::MatrixXd OrthoFunction::testOrtho() double x = xmin + dx*lq->knot(i); double r = x_to_r(x); double f = dx*lq->weight(i)*d_r_to_x(x)*pow(r, dof-1); - double y = 2.0*r/scale; if (segment) y = x; - // Evaluate the unnormalized polynomial - Eigen::VectorXd p = poly_eval(y, nmax) * W(r); + // Evaluate the unnnormalized polynomial + Eigen::VectorXd p = poly_eval(r, nmax) * W(r); // Compute scalar products with the normalizations for (int n1=0; n1<=nmax; n1++) { for (int n2=0; n2<=nmax; n2++) { - ret(n1, n2) += f * p(n1) * p(n2) / sqrt(norm[n1]*norm[n2]); + ret(n1, n2) += f * p(n1) * p(n2) / sqrt(norm(n1)*norm(n2)); } } } @@ -111,13 +110,12 @@ void OrthoFunction::dumpOrtho(const std::string& filename) // double x = xmin + dx*i/(number-1); double r = x_to_r(x); - double y = 2.0*r/scale; if (segment) y = x; // Evaluate polynomial and apply the normalization // - Eigen::VectorXd p = poly_eval(y, nmax) * W(r); + Eigen::VectorXd p = poly_eval(r, nmax) * W(r); fout << std::setw(16) << r; - for (int n=0; n(N); + generate(); } //! Test orthogonality From 97b81c5523c0f6a4053775e8f82df7b074fe1153 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 9 May 2025 21:23:23 -0400 Subject: [PATCH 02/15] Fix a typo in one of the docstrings --- pyEXP/CoefWrappers.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc index fc4fa42dd..498f78cfa 100644 --- a/pyEXP/CoefWrappers.cc +++ b/pyEXP/CoefWrappers.cc @@ -1646,7 +1646,7 @@ void CoefficientClasses(py::module &m) { Returns ------- numpy.ndarray - 4-dimensional numpy array containing the spherical coefficients + 4-dimensional numpy array containing the cylindrical coefficients )"); From b2e8b79ddf58211488d7131d3aa084e91cba5b66 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 9 May 2025 21:25:06 -0400 Subject: [PATCH 03/15] Fix a typo in the trapezoidal rule computation of the tapered Mestel mass; this doesn't affect any computation in EXP so far --- exputil/mestel.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exputil/mestel.cc b/exputil/mestel.cc index 926edbc24..0bbca9b69 100644 --- a/exputil/mestel.cc +++ b/exputil/mestel.cc @@ -158,7 +158,7 @@ double TaperedMestelDisk::get_mass(double R) for (int i=1; i<=N; i++) { double rr = rmin + dr*i; double cur = get_density(rr); - cum += 0.5*dr*(lst + cur) * 2.0*M_PI*rr; + cum += 0.5*dr*(lst*(rr-dr) + cur*rr) * 2.0*M_PI; lst = cur; vecR(i) = rr; From 8261ff870808d0c91b3cfbbbb2cf7b0b2d5c205e Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 9 May 2025 21:26:47 -0400 Subject: [PATCH 04/15] Remove OutVel specific keys from YAML config before forwarding the VelocityBasis constructor --- src/OutVel.cc | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/OutVel.cc b/src/OutVel.cc index 545d99b0d..f6aee96c0 100644 --- a/src/OutVel.cc +++ b/src/OutVel.cc @@ -65,6 +65,12 @@ OutVel::OutVel(const YAML::Node& conf) : Output(conf) YAML::Node node; node["parameters"] = conf; + // Remove OutVel only keys + // + for (auto s : {"nint", "nintsub", "name"}) { + if (node["parameters"][s]) node["parameters"].remove(s); + } + basis = std::make_shared(node); // Create the coefficient container based on the dimensionality. From 225011069f5e92266cf6bbb536546e4e77e4c991 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 9 May 2025 21:32:35 -0400 Subject: [PATCH 05/15] Two main changes: (1) eliminate the 'weight' field since the 'density' field sufficies'; (2) fix a shadowing error in the density file construction that did not set the class version of 'rmin' and 'rmax' from the file resulting in a poor basis construction if 'rmin' and 'rmax' were *not* set in the configuration file as well. --- expui/FieldBasis.cc | 82 ++++++++++++++++++++++++++------------------- 1 file changed, 48 insertions(+), 34 deletions(-) diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index 7c738f677..e9f13db52 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -33,6 +33,7 @@ namespace BasisClasses "ascl", "delta", "lmax", + "mmax", "nmax", "model" }; @@ -54,7 +55,7 @@ namespace BasisClasses // Allocate coefficient storage // - nfld = p.size() + 2; + nfld = p.size() + 1; allocateStore(); // Okay to register @@ -67,7 +68,7 @@ namespace BasisClasses // void FieldBasis::configure() { - nfld = 2; // Weight and density fields by default + nfld = 1; // Density field by default lmax = 4; nmax = 10; rmin = 1.0e-4; @@ -137,7 +138,8 @@ namespace BasisClasses // Compute interpolation functionoid // - double rmin = r.front(), rmax = r.back(); + rmin = r.front(); + rmax = r.back(); interp = std::make_shared(r, d); densfunc = [this](double r) @@ -168,7 +170,6 @@ namespace BasisClasses // Initialize fieldlabels // fieldLabels.clear(); - fieldLabels.push_back("weight"); fieldLabels.push_back("density"); // Debug @@ -204,9 +205,11 @@ namespace BasisClasses void FieldBasis::initialize() { - // Remove matched keys + // Check for unmatched keys // - // for (auto v : valid_keys) current_keys.erase(v); + auto unmatched = YamlCheck(conf, valid_keys); + if (unmatched.size()) + throw YamlConfigError("Basis::FieldBasis", "parameter", unmatched, __FILE__, __LINE__); // Assign values from YAML // @@ -215,6 +218,7 @@ namespace BasisClasses if (conf["model" ]) model = conf["model" ].as(); if (conf["nfld" ]) nfld = conf["nfld" ].as(); if (conf["lmax" ]) lmax = conf["lmax" ].as(); + if (conf["mmax" ]) lmax = conf["mmax" ].as(); if (conf["nmax" ]) nmax = conf["nmax" ].as(); if (conf["dof" ]) dof = conf["dof" ].as(); if (conf["rmin" ]) rmin = conf["rmin" ].as(); @@ -312,8 +316,7 @@ namespace BasisClasses double x, double y, double z, double u, double v, double w) { - constexpr std::complex I(0, 1); - constexpr double fac0 = 0.25*M_2_SQRTPI; + constexpr double fac2 = 0.5*M_2_SQRTPI/M_SQRT2; // 1/sqrt(2*pi)=0.3989422804 int tid = omp_get_thread_num(); PS3 pos{x, y, z}, vel{u, v, w}; @@ -337,18 +340,16 @@ namespace BasisClasses auto p = (*ortho)(R); - (*coefs[tid])(0, 0, 0) += mass*p(0)*fac0; - for (int m=0; m<=lmax; m++) { - std::complex P = std::exp(-I*(phi*m)); + auto P = std::exp(std::complex(0, -phi*m))*fac2; for (int n=0; n P = - std::exp(-I*(phi*m))*Ylm(l, m)*plgndr(l, m, cth) * s; - - s *= -1.0; // Flip sign for next evaluation + auto P = std::exp(std::complex(0, -phi*m)) * + Ylm(l, m)*plgndr(l, m, cth) * s; + + // Flip sign for next evaluation + s *= -1.0; for (int n=0; n(); - else coefret = std::make_shared(); - // Add the configuration data + else coefret = std::make_shared(); + + // Add the configuration data + // std::ostringstream sout; sout << node; coefret->buf = sout.str(); } @@ -419,10 +421,14 @@ namespace BasisClasses // if (dof==2) { auto p = dynamic_pointer_cast(coefret); - p->assign(store[0], nfld, lmax, nmax); + int rows = lmax + 1; + int cols = nmax; + p->assign(store[0], nfld, rows, cols); } else { auto p = dynamic_pointer_cast(coefret); - p->assign(store[0], nfld, lmax, nmax); + int rows = (lmax+1)*(lmax+2)/2; + int cols = nmax; + p->assign(store[0], nfld, rows, cols); } } @@ -431,11 +437,15 @@ namespace BasisClasses { if (dof==2) { auto cstr = std::make_shared(); - cstr->assign(store[0], nfld, lmax, nmax); + int rows = lmax + 1; + int cols = nmax; + cstr->assign(store[0], nfld, rows, cols); return cstr; } else { auto cstr = std::make_shared(); - cstr->assign(store[0], nfld, lmax, nmax); + int rows = (lmax+1)*(lmax+2)/2; + int cols = nmax; + cstr->assign(store[0], nfld, rows, cols); return cstr; } } @@ -443,18 +453,21 @@ namespace BasisClasses std::vector FieldBasis::sph_eval(double r, double costh, double phi) { - constexpr std::complex I(0, 1); + constexpr double fac2 = 0.5*M_2_SQRTPI/M_SQRT2; // 1/sqrt(2 pi) if (dof==2) { std::vector ret(nfld, 0); - auto p = (*ortho)(r*sqrt(fabs(1.0 - costh*costh))); + auto p = (*ortho)(r); for (int i=0; i0 ? 2.0 : 1.0; + auto P = std::exp(std::complex(0, -phi*m)); + + std::complex sum = 0.0; for (int n=0; n(0, -phi*m)); + double fac = m>0 ? 2.0 : 1.0; for (int n=0; n Date: Fri, 9 May 2025 21:33:12 -0400 Subject: [PATCH 06/15] Added many more tests of the orthogonal function generation and reconstruction --- utils/SL/oftest.cc | 472 +++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 435 insertions(+), 37 deletions(-) diff --git a/utils/SL/oftest.cc b/utils/SL/oftest.cc index 1426d3ae3..0f6b229df 100644 --- a/utils/SL/oftest.cc +++ b/utils/SL/oftest.cc @@ -1,6 +1,7 @@ #include #include #include +#include #include #include #include @@ -8,36 +9,71 @@ #include #include // Orthogonal basis +#include // Linear interpolation #include -// Generate a 2d exponential distribution +// Recursion relation evaluation for generalized Laguerre functions // -class genE +double generalized_laguerre(int n, double alpha, double x) { -private: + if (n == 0) { + return 1.0; + } else if (n == 1) { + return 1.0 + alpha - x; + } else { + return ((2.0 * n - 1.0 + alpha - x) * generalized_laguerre(n - 1, alpha, x) - (n - 1.0 + alpha) * generalized_laguerre(n - 2, alpha, x)) / n; + } +} + + +// Class to generate 2d distribution +// +class gen2D +{ +protected: static constexpr double TOL = 1.0e-12; std::shared_ptr gen; std::shared_ptr> uniform; - double R, a; - double f(double x) { return R - (1.0 - (1.0 + x)*exp(-x)); }; - double df(double x){ return -x*exp(-x); }; - + double a, phi, alpha; + int m; + public: - genE(double a) : a(a) { + gen2D(double a, int m, double phi, double alpha) : + a(a), m(m), phi(phi), alpha(alpha) + { unsigned seed = std::random_device{}(); gen = std::make_shared(seed); uniform = std::make_shared>(0.0, 1.0); } - std::tuple operator()() + + virtual std::tuple operator()() = 0; +}; + +// Generate a 2d exponential distribution +// +class genE : public gen2D +{ +private: + + double R; + double f(double x) { return R - (1.0 - (1.0 + x)*exp(-x)); }; + double df(double x){ return -x*exp(-x); }; + +public: + + genE(double a, int m, double phi, double alpha) : gen2D(a, m, phi, alpha) {} + + std::tuple operator()() { R = (*uniform)(*gen); double x = sqrt(R); + // Realize a radius int n=0; for (; n<100; n++) { double delF = - f(x)/df(x); @@ -45,18 +81,143 @@ class genE if (fabs(delF) (*uniform)(*gen)) break; + P = 2.0*M_PI*(*uniform)(*gen); + } + + return {x*a, P, n, k}; + } +}; + +// Generate a distribution from arrays +// +class genO : public gen2D +{ +private: + + std::unique_ptr interp; + +public: + + genO(const std::string& filename, + int m, double phi, double alpha) : gen2D(1.0, m, phi, alpha) + { + std::ifstream in (filename); + std::vector r, d; + if (in) { + std::string line; + while (std::getline(in, line)) { + std::istringstream iss(line); + double x, y; + iss >> x >> y; + if (iss) { + r.push_back(x); + d.push_back(y); + } else break; + } + } + else throw std::runtime_error("Cannot open file " + filename); + + std::vector pp(r.size(), 0.0); + for (int n=1; n(pp, r); + // Debug + std::ofstream out("genO.dat"); + if (out) { + for (int n=0; n r, std::vector d, + int m, double phi, double alpha) : gen2D(1.0, m, phi, alpha) + { + std::vector pp(r.size(), 0.0); + for (int n=1; n(pp, r); + // Debug + std::ofstream out("genO.dat"); + if (out) { + for (int n=0; n operator()() + { + double r = interp->eval((*uniform)(*gen)); + + // Realize an angle + double P = 2.0*M_PI*(*uniform)(*gen); + + int k=0; + if (m>0) { + for (; k<1000; k++) { + double fp = 0.5*(1.0 + cos((P-phi-alpha*r)*m)); + if (fp > (*uniform)(*gen)) break; + P = 2.0*M_PI*(*uniform)(*gen); + } + } + + return {r, P, 0, k}; + } +}; + +// Generate a 2d uniform distribution; a uniform disk with an edge at a. +// +class genU : public gen2D +{ +private: + +public: + + genU(double a, int m, double phi, double alpha) : gen2D(a, m, phi, alpha) {} + + std::tuple operator()() + { + double R = (*uniform)(*gen); + double x = sqrt(R); + + // Realize an angle + double P = 2.0*M_PI*(*uniform)(*gen); + int k=0; + for (; k<1000; k++) { + double fp = 0.5*(1.0 + cos((P-phi-alpha*x)*m)); + if (fp > (*uniform)(*gen)) break; + P = 2.0*M_PI*(*uniform)(*gen); + } + + return {x*a, P, 0, k}; } }; int main(int argc, char** argv) { bool logr = false; - int numr, nmax, nout, N; - double A, rmin, rmax, scale, delta; + int numr, nmax, nout, N, M, mmax, index, knots; + double A, rmin, rmax, scale, delta, pitch, phi, Rout; unsigned seed; // Will be inialized by /dev/random if // not set on the command line - std::string filename; + std::string filename, modelname, coeffile, densprof, bodyfile; // Parse command line // @@ -66,12 +227,23 @@ int main(int argc, char** argv) options.add_options() ("h,help", "Print this help message") ("logr", "Plot output grid with logarithmic spacing") + ("uniform", "Generate a uniform distribution rather than an exponential") + ("i,index", "Use an orthogonal basis of index i for realization", + cxxopts::value(index)->default_value("0")) ("A,length", "characteristic disk scale length", cxxopts::value(A)->default_value("1.0")) ("D,delta", "truncation smoothing scale", cxxopts::value(delta)->default_value("0.005")) + ("P,phi", "Position angle for M>0 distribution (degrees)", + cxxopts::value(phi)->default_value("45.0")) ("N,mc", "number of particles for Monte Carlo realization", cxxopts::value(N)->default_value("10000")) + ("M", "azimuthal order for Monte Carlo realization", + cxxopts::value(M)->default_value("0")) + ("k,knots", "Inner product knots for Stieltjes procedure", + cxxopts::value(knots)->default_value("400")) + ("mmax", "maximum azimuthal order for reconstrution", + cxxopts::value(mmax)->default_value("4")) ("nmax", "maximum number of radial harmonics in the expansion", cxxopts::value(nmax)->default_value("18")) ("r,rmin", "minimum radius for the SL grid", @@ -79,13 +251,25 @@ int main(int argc, char** argv) ("R,rmax", "maximum radius for the SL grid", cxxopts::value(rmax)->default_value("20.0")) ("s,scale", "ortho function map scale factor", - cxxopts::value(scale)->default_value("0.01")) + cxxopts::value(scale)->default_value("1.0")) + ("p,pitch", "dPhi/dR in scale length units", + cxxopts::value(pitch)->default_value("0.0")) ("nout", "number of points in the output grid per side", cxxopts::value(nout)->default_value("40")) - ("seed", "Random number seed. Default: use /dev/random", + ("rout", "half-size of output grid in scale lenghts", + cxxopts::value(Rout)->default_value("4.0")) + ("seed", "Random number seed. Default: use /dev/random", cxxopts::value(seed)) + ("model", "Input model file", + cxxopts::value(modelname)) ("o,filename", "Output filename", - cxxopts::value(filename)->default_value("testof")) + cxxopts::value(filename)->default_value("oftest")) + ("c,coeffile", "Input coefficient file", + cxxopts::value(coeffile)) + ("d,densprof", "Density profile for sampling", + cxxopts::value(densprof)) + ("b,bodyfile", "Body file for sampling", + cxxopts::value(bodyfile)) ; //=================== @@ -113,22 +297,68 @@ int main(int argc, char** argv) seed = std::random_device{}(); } + phi *= M_PI/180.0; + // Log spacing? // if (vm.count("logr")) logr = true; - // Create the weight function instance + // File or expon // - auto densfunc = [&](double r) - { - return exp(-r/A) * - 0.5*(1.0 + std::erf((rmax - 5.0*delta - r)/delta)) / (A*A); - }; - + std::function densfunc; + std::shared_ptr interp, interp2; + + // Read a model file + // + if (vm.count("model")) { + std::vector r, d; + std::ifstream in(modelname); + if (not in) throw std::runtime_error("Error opening file: " + modelname); + + std::string line; + while (std::getline(in, line)) { + auto pos = line.find_first_of("!#"); + if (pos == std::string::npos) { + std::istringstream iss(line); + double x, y; + iss >> x >> y; + if (iss) { + r.push_back(x); + d.push_back(y); + } + } + } + + // Compute interpolation functionoid + // + double rmin = r.front(), rmax = r.back(); + + interp = std::make_shared(r, d); + densfunc = [&interp](double r) + { + return interp->eval(r); + }; + } + // Create the weight function for expon + // + else { + if (vm.count("uniform")) + densfunc = [&](double r) + { + return 1.0; + }; + else + densfunc = [&](double r) + { + return exp(-r/A) * + 0.5*(1.0 + std::erf((rmax - 5.0*delta - r)/delta)) / (A*A); + }; + } + // Generate the orthogonal function instance // OrthoFunction ortho(nmax, densfunc, rmin, rmax, scale, 2); - + ortho.setKnots(knots); // Output file for grid // @@ -136,7 +366,7 @@ int main(int argc, char** argv) // Define some representative limits // - double Rmax = 4.0*A; + double Rmax = Rout*A; // Grid spacing // @@ -148,32 +378,200 @@ int main(int argc, char** argv) // Print the grid // + double tmass = 0.0; for (int i=0; i>> coef(mmax+1); + for (int m=0; m<=mmax; m++) coef[m].resize(nmax, 0.0); + + // Read a coefficient file + // + if (vm.count("coeffile")) { + std::ifstream in(coeffile); + if (in.good()) { + std::string line; + while (std::getline(in, line)) { + std::istringstream sout(line); + int n; sout >> n; + if (sout and n < nmax) { + std::vector> c; + while (sout) { + double r, v; + sout >> r >> v; + c.push_back({r, v}); + } + for (int mm=0; mm<=std::min(mmax, c.size()); mm++) { + coef[mm][n] = c[mm]; + } + } + } + } else { + throw std::runtime_error("Error opening coefficient file: " + coeffile); + } + + N = 0; + } + + if (N or vm.count("bodyfile")) { + // Distribution generator + std::unique_ptr gen; + if (vm.count("uniform")) + gen = std::make_unique(A, M, phi, pitch); + else if (vm.count("index")) { + double ratio = 1.0; + const int ngrid = 1000; + double lrmin = log(rmin), lrmax = log(rmax); + double dr = (lrmax - lrmin)/(ngrid - 1); + for (int n=0; n rr(ngrid), dd(ngrid); + for (int n=0; n(rr, dd, M, phi, pitch); + } + else if (vm.count("densprof")) + gen = std::make_unique(densprof, M, phi, pitch); + else + gen = std::make_unique(A, M, phi, pitch); + + if (vm.count("bodyfile")) { + std::ifstream in(bodyfile); + if (in) { + std::string line; + std::getline(in, line); + while (std::getline(in, line)) { + std::istringstream iss(line); + double ms, x, y; + iss >> ms >> x >> y; + auto p = ortho(sqrt(x*x + y*y)); + double P = atan2(y, x); + for (int m=0; m<=mmax; m++) { + auto azi = std::exp(std::complex(0, -P*m)); + for (int i=0; i(0, -P*m)); + for (int i=0; i tot = 0.0; + for (int m=0; m<=mmax; m++) { + std::complex c = 0.0; + for (int i=0; i(0, phi*m)); + out1 << std::setw(18) << std::real(c); + tot += c; + } + out1 << std::setw(18) << std::real(tot) << std::endl; + } + + auto p = ortho(std::abs(y)); + double phi = 0.5*M_PI; + if (y<0.0) phi = -0.5*M_PI; + + std::complex tot = 0.0; + for (int m=0; m<=mmax; m++) { + std::complex c = 0.0; + for (int i=0; i(0, phi*m)); + out << std::setw(18) << std::real(c); + tot += c; + } + out2 << std::setw(18) << y << std::setw(18) << std::real(tot) << std::endl; + } + } + + std::cout << std::endl + << "Orthogonality of the function at the grid points" + << std::endl; + + ortho.dumpOrtho("oftest.dump"); + + if (vm.count("model")==0 and vm.count("uniform")==0) { + std::ofstream out("oftest.laguerre"); + const int ngrid = 1000; + double alpha = 1.0; + double lrmin = log(rmin), lrmax = log(rmax); + double dr = (lrmax - lrmin)/(ngrid - 1); + for (int j=0; j Date: Sun, 11 May 2025 19:05:49 -0400 Subject: [PATCH 07/15] Correct dimension updates for CoefStruct for latest version --- expui/FieldBasis.cc | 20 +++++++------------- 1 file changed, 7 insertions(+), 13 deletions(-) diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index e9f13db52..1d9fa91ac 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -55,7 +55,7 @@ namespace BasisClasses // Allocate coefficient storage // - nfld = p.size() + 1; + nfld = p.size() + 1; // Add the density to the phase-space return allocateStore(); // Okay to register @@ -287,7 +287,7 @@ namespace BasisClasses // Sanity test dimensions if (nfld!=p->nfld || lmax!=p->mmax || nmax!=p->nmax) { std::ostringstream serr; - serr << "FieldBasis::set_coefs: dimension error! " + serr << "FieldBasis::set_coefs: dimension error for dof=2! " << " nfld [" << nfld << "!= " << p->nfld << "]" << " mmax [" << lmax << "!= " << p->mmax << "]" << " nmax [" << nmax << "!= " << p->nmax << "]"; @@ -302,7 +302,7 @@ namespace BasisClasses // Sanity test dimensions if (nfld!=p->nfld || lmax!=p->lmax || nmax!=p->nmax) { std::ostringstream serr; - serr << "FieldBasis::set_coefs: dimension error! " + serr << "FieldBasis::set_coefs: dimension error for dof=3! " << " nfld [" << nfld << "!= " << p->nfld << "]" << " lmax [" << lmax << "!= " << p->lmax << "]" << " nmax [" << nmax << "!= " << p->nmax << "]"; @@ -421,14 +421,10 @@ namespace BasisClasses // if (dof==2) { auto p = dynamic_pointer_cast(coefret); - int rows = lmax + 1; - int cols = nmax; - p->assign(store[0], nfld, rows, cols); + p->assign(store[0], nfld, lmax, nmax); } else { auto p = dynamic_pointer_cast(coefret); - int rows = (lmax+1)*(lmax+2)/2; - int cols = nmax; - p->assign(store[0], nfld, rows, cols); + p->assign(store[0], nfld, lmax, nmax); } } @@ -439,13 +435,11 @@ namespace BasisClasses auto cstr = std::make_shared(); int rows = lmax + 1; int cols = nmax; - cstr->assign(store[0], nfld, rows, cols); + cstr->assign(store[0], nfld, lmax, nmax); return cstr; } else { auto cstr = std::make_shared(); - int rows = (lmax+1)*(lmax+2)/2; - int cols = nmax; - cstr->assign(store[0], nfld, rows, cols); + cstr->assign(store[0], nfld, lmax, nmax); return cstr; } } From 92f6c20d4824b0b49de493e383b66855a151aa47 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 12 May 2025 10:30:00 -0400 Subject: [PATCH 08/15] Ask OrthoFunction for the correct number of basis functions --- expui/FieldBasis.cc | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index 1d9fa91ac..d8f86abe5 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -165,7 +165,11 @@ namespace BasisClasses // Generate the orthogonal function instance // - ortho = std::make_shared(nmax, densfunc, rmin, rmax, rmapping, dof); + ortho = std::make_shared + (nmax-1, densfunc, rmin, rmax, rmapping, dof); + // ^ + // | + // +--- This is the polynmial order, not the rank // Initialize fieldlabels // From 463f4ae82e649a63949f2b6ac3bdc7d70f1d89f4 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 12 May 2025 10:31:23 -0400 Subject: [PATCH 09/15] Don't use shallow copies of Eigen structures which can have intended consquences in the C++ code from the pybind11 mapping --- pyEXP/CoefWrappers.cc | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc index 498f78cfa..c4559acc5 100644 --- a/pyEXP/CoefWrappers.cc +++ b/pyEXP/CoefWrappers.cc @@ -1365,7 +1365,7 @@ void CoefficientClasses(py::module &m) { .def("getAllCoefs", [](CoefClasses::SphCoefs& A) { - auto M = A.getAllCoefs(); // Need a copy here + Eigen::Tensor, 3> M = A.getAllCoefs(); // Need a copy here py::array_t> ret = make_ndarray3>(M); return ret; }, @@ -1440,7 +1440,7 @@ void CoefficientClasses(py::module &m) { .def("getAllCoefs", [](CoefClasses::CylCoefs& A) { - auto M = A.getAllCoefs(); // Need a copy here + Eigen::Tensor, 3> M = A.getAllCoefs(); // Need a copy here py::array_t> ret = make_ndarray3>(M); return ret; }, @@ -1552,7 +1552,7 @@ void CoefficientClasses(py::module &m) { .def("getAllCoefs", [](CoefClasses::SphFldCoefs& A) { - auto M = A.getAllCoefs(); // Need a copy here + Eigen::Tensor, 4> M = A.getAllCoefs(); // Need a copy here py::array_t> ret = make_ndarray4>(M); return ret; }, @@ -1636,7 +1636,7 @@ void CoefficientClasses(py::module &m) { .def("getAllCoefs", [](CoefClasses::CylFldCoefs& A) { - auto M = A.getAllCoefs(); // Need a copy here + Eigen::Tensor, 4> M = A.getAllCoefs(); // Need a copy here py::array_t> ret = make_ndarray4>(M); return ret; }, @@ -1706,7 +1706,7 @@ void CoefficientClasses(py::module &m) { .def("getAllCoefs", [](CoefClasses::SlabCoefs& A) { - auto M = A.getAllCoefs(); // Need a copy here + Eigen::Tensor, 4> M = A.getAllCoefs(); // Need a copy here py::array_t> ret = make_ndarray4>(M); return ret; }, @@ -1800,7 +1800,7 @@ void CoefficientClasses(py::module &m) { .def("getAllCoefs", [](CoefClasses::CubeCoefs& A) { - auto M = A.getAllCoefs(); // Need a copy here + Eigen::Tensor, 4> M = A.getAllCoefs(); // Need a copy here py::array_t> ret = make_ndarray4>(M); return ret; }, @@ -1897,7 +1897,7 @@ void CoefficientClasses(py::module &m) { ------- TableData instance )", py::arg("time"), py::arg("array"), py::arg("verbose")=true) - .def("getAllCoefs", &CoefClasses::TableData::getAllCoefs, + .def("getAllCoefs", &CoefClasses::TableData::getAllCoefs, R"( Return a 2-dimensional ndarray indexed by column and time @@ -1967,9 +1967,16 @@ void CoefficientClasses(py::module &m) { ------- TrajectoryData instance )", py::arg("time"), py::arg("array"), py::arg("verbose")=true) - .def("getAllCoefs", &CoefClasses::TrajectoryData::getAllCoefs, + .def("getAllCoefs", + [](CoefClasses::TrajectoryData& A) + { + Eigen::Tensor M = A.getAllCoefs(); // Need a copy here + py::array_t ret = make_ndarray3(M); + return ret; + }, + R"( - Return a 2-dimensional ndarray indexed by column and time + Return a 3-dimensional ndarray indexed by column and time Returns ------- From 72039a84532f0cf981ce69ed56b94d8dc5f4e5bf Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 12 May 2025 10:31:43 -0400 Subject: [PATCH 10/15] Add more comments only --- include/OrthoFunction.H | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/include/OrthoFunction.H b/include/OrthoFunction.H index 7d5523ee6..e344efedf 100644 --- a/include/OrthoFunction.H +++ b/include/OrthoFunction.H @@ -78,16 +78,26 @@ private: public: - //! Constructor - OrthoFunction(int nmax, DensFunc W, double rmin, double rmax, double scale, int dof); + /** Constructor + @param norder Order of the orthogonal polynomial + @param W Density weight function + @param rmin Inner radial bound + @param rmax Outer radial bound + @param scale Mapping scale + @param dof Degrees of freedom (2 or 3) + */ + OrthoFunction(int norder, DensFunc W, double rmin, double rmax, double scale, int dof); //! Destructor virtual ~OrthoFunction() {} //! Evaluate orthogonal functions at r + //! @param r Radial coordinate + //! @return A vector of the orthogonal functions Eigen::VectorXd operator()(double r); //! Reset Legendre knots and weights for inner product + //! @param N Number of Legendre knots void setKnots(int N) { knots = N; @@ -99,6 +109,7 @@ public: Eigen::MatrixXd testOrtho(); //! Dump the orthogonal function table + //! @param filename Filename for output void dumpOrtho(const std::string& filename); }; From 4a54f93534c1f55b9b7fe721fa0984bab9775289 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 12 May 2025 11:19:47 -0400 Subject: [PATCH 11/15] Remove cruft --- expui/FieldBasis.cc | 2 -- 1 file changed, 2 deletions(-) diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index d8f86abe5..91858e768 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -437,8 +437,6 @@ namespace BasisClasses { if (dof==2) { auto cstr = std::make_shared(); - int rows = lmax + 1; - int cols = nmax; cstr->assign(store[0], nfld, lmax, nmax); return cstr; } else { From b23593cad131dd76b7081ee71312e8b7f177a185 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 12 May 2025 11:25:09 -0400 Subject: [PATCH 12/15] Use Eigen3 index operator for clarity --- exputil/OrthoFunction.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/exputil/OrthoFunction.cc b/exputil/OrthoFunction.cc index 79f9c5951..d3048fcd6 100644 --- a/exputil/OrthoFunction.cc +++ b/exputil/OrthoFunction.cc @@ -33,7 +33,7 @@ void OrthoFunction::generate() // Remaining values for (int i=1; i<=nmax; i++) { norm(i) = scalar_prod(i, 0); - alph(i) = scalar_prod(i, 1)/norm[i]; + alph(i) = scalar_prod(i, 1)/norm(i); beta(i) = norm(i)/norm(i-1); } } @@ -58,11 +58,11 @@ double OrthoFunction::scalar_prod(int n, int m) Eigen::VectorXd OrthoFunction::poly_eval(double t, int n) { Eigen::VectorXd ret(n+1); - ret[0] = 1.0; + ret(0) = 1.0; if (n) { - ret[1] = (t - alph[0])*ret[0]; + ret(1) = (t - alph(0))*ret(0); for (int j=1; j Date: Mon, 12 May 2025 11:28:50 -0400 Subject: [PATCH 13/15] Spelling in comment --- exputil/OrthoFunction.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exputil/OrthoFunction.cc b/exputil/OrthoFunction.cc index d3048fcd6..59316b1f5 100644 --- a/exputil/OrthoFunction.cc +++ b/exputil/OrthoFunction.cc @@ -82,7 +82,7 @@ Eigen::MatrixXd OrthoFunction::testOrtho() double r = x_to_r(x); double f = dx*lq->weight(i)*d_r_to_x(x)*pow(r, dof-1); - // Evaluate the unnnormalized polynomial + // Evaluate the unnormalized polynomial Eigen::VectorXd p = poly_eval(r, nmax) * W(r); // Compute scalar products with the normalizations From 755610535f2d7c7618028111fa69231212b6bf48 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Tue, 13 May 2025 10:00:16 +0100 Subject: [PATCH 14/15] comment typo; warning suppression --- expui/FieldBasis.cc | 4 ++-- utils/SL/oftest.cc | 2 ++ 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index 91858e768..ff59c2016 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -166,10 +166,10 @@ namespace BasisClasses // Generate the orthogonal function instance // ortho = std::make_shared - (nmax-1, densfunc, rmin, rmax, rmapping, dof); + (nmax-1, densfunc, rmin, rmax, rmapping, dof); // ^ // | - // +--- This is the polynmial order, not the rank + // +--- This is the polynomial order, not the rank // Initialize fieldlabels // diff --git a/utils/SL/oftest.cc b/utils/SL/oftest.cc index 0f6b229df..171cb020a 100644 --- a/utils/SL/oftest.cc +++ b/utils/SL/oftest.cc @@ -52,6 +52,8 @@ class gen2D virtual std::tuple operator()() = 0; + + virtual ~gen2D() = default; }; // Generate a 2d exponential distribution From f145ea77f9c2a935aa4b975fe6f0c78d5d214705 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Tue, 13 May 2025 10:02:27 +0100 Subject: [PATCH 15/15] version bump [skip ci] --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 7f85582c6..8dacf8e78 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.25) # Needed for CUDA, MPI, and CTest features project( EXP - VERSION "7.8.4" + VERSION "7.8.5" HOMEPAGE_URL https://github.com/EXP-code/EXP LANGUAGES C CXX Fortran)