From 6de7e42bfdd12deab21c40e3816ec9f10dbd340b Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 9 Apr 2025 16:16:49 -0400 Subject: [PATCH 1/4] Rewrite step counting logic in IntegrateOrbits to always provide ODE between interval end points --- expui/BiorthBasis.cc | 16 +++++++++++----- pyEXP/BasisWrappers.cc | 4 ++-- 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index f9237bc65..25f20ab1e 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3899,10 +3899,14 @@ namespace BasisClasses // int numT = floor( (tfinal - tinit)/h ); - // Compute output step + // Number of output steps. Will attempt to find the best stride... // - nout = std::min(numT, nout); - double H = (tfinal - tinit)/nout; + int stride = std::ceil(static_cast(numT)/static_cast(nout)); + if (stride>1) numT = nout * stride; + + // Compute the interval-matching step + // + h = (tfinal - tinit)/(numT - 1); // Return data // @@ -3934,9 +3938,11 @@ namespace BasisClasses for (int k=0; k<6; k++) ret(n, k, 0) = ps(n, k); double tnow = tinit; - for (int s=1, cnt=1; s= H*cnt-h*1.0e-8) { + if (s++ % stride == 0) { times(cnt) = tnow; for (int n=0; n nout Returns ------- From fe274b72ce1a8e850acb9e41d7504f9f766c0979 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 10 Apr 2025 08:08:53 -0400 Subject: [PATCH 2/4] Final checked version --- expui/BiorthBasis.cc | 60 +++++++++++++++++++++++++++++++++----------- 1 file changed, 46 insertions(+), 14 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 25f20ab1e..ab9945ca7 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3885,11 +3885,29 @@ namespace BasisClasses // Sanity check // + if (tfinal == tinit) { + throw std::runtime_error + ("BasisClasses::IntegrateOrbits: tinit cannot be equal to tfinal"); + } + + if (h < 0.0 and tfinal > tinit) { + std::ostringstream sout; + throw std::runtime_error + ("BasisClasses::IntegrateOrbits: tfinal must be smaller than tinit " + "when step size is negative"); + } + + if (h > 0.0 and tfinal < tinit) { + throw std::runtime_error + ("BasisClasses::IntegrateOrbits: tfinal must be larger than " + "tinit when step size is positive"); + } + if ( (tfinal - tinit)/h > static_cast(std::numeric_limits::max()) ) { - std::cout << "BasicFactor::IntegrateOrbits: step size is too small or " - << "time interval is too large.\n"; + std::cout << "BasisClasses::IntegrateOrbits: step size is too small or " + << "time interval is too large." << std::endl; // Return empty data // return {Eigen::VectorXd(), Eigen::Tensor()}; @@ -3897,16 +3915,22 @@ namespace BasisClasses // Number of steps // - int numT = floor( (tfinal - tinit)/h ); + int numT = std::ceil( (tfinal - tinit)/h ); + + // Want both end points in the output at minimum + // + numT = std::max(2, numT); + nout = std::max(2, nout); - // Number of output steps. Will attempt to find the best stride... + // Number of output steps. Compute a stride>1 if nout(numT)/static_cast(nout)); - if (stride>1) numT = nout * stride; + if (stride>1) numT = (nout-1) * stride + 1; + else nout = numT; // Compute the interval-matching step // - h = (tfinal - tinit)/(numT - 1); + h = (tfinal - tinit)/(numT-1); // Return data // @@ -3916,10 +3940,10 @@ namespace BasisClasses ret.resize(rows, 6, nout); } catch (const std::bad_alloc& e) { - std::cout << "BasicFactor::IntegrateOrbits: memory allocation failed: " + std::cout << "BasisClasses::IntegrateOrbits: memory allocation failed: " << e.what() << std::endl << "Your requested number of orbits and time steps requires " - << floor(8.0*rows*6*nout/1e9)+1 << " GB free memory" + << floor(4.0*rows*6*nout/stride/1e9)+1 << " GB free memory" << std::endl; // Return empty data @@ -3931,18 +3955,24 @@ namespace BasisClasses // Eigen::VectorXd times(nout); - // Do the work + // Assign the initial point // times(0) = tinit; for (int n=0; n Date: Thu, 10 Apr 2025 17:57:21 -0400 Subject: [PATCH 3/4] Some very minor clean up of warning output --- expui/BiorthBasis.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index ab9945ca7..0c29596f2 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3891,7 +3891,6 @@ namespace BasisClasses } if (h < 0.0 and tfinal > tinit) { - std::ostringstream sout; throw std::runtime_error ("BasisClasses::IntegrateOrbits: tfinal must be smaller than tinit " "when step size is negative"); @@ -3943,7 +3942,7 @@ namespace BasisClasses std::cout << "BasisClasses::IntegrateOrbits: memory allocation failed: " << e.what() << std::endl << "Your requested number of orbits and time steps requires " - << floor(4.0*rows*6*nout/stride/1e9)+1 << " GB free memory" + << std::floor(4.0*rows*6*nout/1e9)+1 << " GB free memory" << std::endl; // Return empty data From 804d665fafbde71a02be03b09d97c33e1580d7b5 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 17 Apr 2025 16:16:54 -0400 Subject: [PATCH 4/4] Fix round off boundary and never ignore nout --- expui/BiorthBasis.H | 2 +- expui/BiorthBasis.cc | 22 ++++++++++++++++------ pyEXP/BasisWrappers.cc | 16 +++++++++------- 3 files changed, 26 insertions(+), 14 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index e2fd16dcd..d7e239433 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -1189,7 +1189,7 @@ namespace BasisClasses std::tuple> IntegrateOrbits (double tinit, double tfinal, double h, Eigen::MatrixXd ps, std::vector bfe, - AccelFunctor F, int nout=std::numeric_limits::max()); + AccelFunctor F, int nout=0); using BiorthBasisPtr = std::shared_ptr; } diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 0c29596f2..7705c7bc1 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3914,23 +3914,33 @@ namespace BasisClasses // Number of steps // - int numT = std::ceil( (tfinal - tinit)/h ); + int numT = std::ceil( (tfinal - tinit)/h + 0.5); // Want both end points in the output at minimum // numT = std::max(2, numT); - nout = std::max(2, nout); - // Number of output steps. Compute a stride>1 if nout(numT)/static_cast(nout)); - if (stride>1) numT = (nout-1) * stride + 1; - else nout = numT; + int stride = 1; // Default stride + if (nout>0) { // User has specified output count... + nout = std::max(2, nout); + stride = std::ceil(static_cast(numT)/static_cast(nout)); + numT = (nout-1) * stride + 1; + } else { // Otherwise, use the default output number + nout = numT; // with the default stride + } // Compute the interval-matching step // h = (tfinal - tinit)/(numT-1); + // DEBUG + if (false) + std::cout << "BasisClasses::IntegrateOrbits: choosing nout=" << nout + << " numT=" << numT << " h=" << h << " stride=" << stride + << std::endl; + // Return data // Eigen::Tensor ret; diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index f792019f2..bb214cf20 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -2194,11 +2194,13 @@ void BasisFactoryClasses(py::module &m) R"( Compute particle orbits in gravitational field from the bases - Integrate a list of initial conditions from tinit to tfinal with a - step size of h using the list of basis and coefficient pairs. Every - step will be included in return unless you provide an explicit - value for 'nout', the number of desired output steps. This will - choose the 'nout' points closest to the desired time. + Integrate a list of initial conditions from 'tinit' to 'tfinal' with + a step size of 'h' using the list of basis and coefficient pairs. The + step size will be adjusted to provide uniform sampling. Every + step will be returned unless you provide an explicit value for 'nout', + the number of desired output steps. In this case, the code will + choose new set step size equal or smaller to the supplied step size + with a stride to provide exactly 'nout' output times. Parameters ---------- @@ -2216,7 +2218,7 @@ void BasisFactoryClasses(py::module &m) func : AccelFunctor the force function nout : int - the number of output intervals if (tfinal - tinit) / h > nout + the number of output points, if specified Returns ------- @@ -2225,5 +2227,5 @@ void BasisFactoryClasses(py::module &m) )", py::arg("tinit"), py::arg("tfinal"), py::arg("h"), py::arg("ps"), py::arg("basiscoef"), py::arg("func"), - py::arg("nout")=std::numeric_limits::max()); + py::arg("nout")=0); }