Merge pull request #125 from EXP-code/tfinalFix

Tfinal fix

Author: michael-petersen

expui/BiorthBasis.H

@@ -1189,7 +1189,7 @@ namespace BasisClasses
   std::tuple<Eigen::VectorXd, Eigen::Tensor<float, 3>>
   IntegrateOrbits (double tinit, double tfinal, double h,
 		   Eigen::MatrixXd ps, std::vector<BasisCoef> bfe,
-		   AccelFunctor F, int nout=std::numeric_limits<int>::max());
+		   AccelFunctor F, int nout=0);
 
   using BiorthBasisPtr = std::shared_ptr<BiorthBasis>;
 }

expui/BiorthBasis.cc

@@ -3885,24 +3885,61 @@ 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) {
+      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<double>(std::numeric_limits<int>::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<float, 3>()};
       }
 
     // Number of steps
     //
-    int numT = floor( (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);
 
-    // Compute output step
+    // Number of output steps
     //
-    nout = std::min<int>(numT, nout);
-    double H = (tfinal - tinit)/nout;
+    int stride = 1;		// Default stride
+    if (nout>0) {		// User has specified output count...
+      nout = std::max(2, nout);
+      stride = std::ceil(static_cast<double>(numT)/static_cast<double>(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
     //
@@ -3912,10 +3949,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"
+		<< std::floor(4.0*rows*6*nout/1e9)+1 << " GB free memory"
 		<< std::endl;
 
       // Return empty data
@@ -3927,27 +3964,37 @@ namespace BasisClasses
     //
     Eigen::VectorXd times(nout);
 
-    // Do the work
+    // Assign the initial point
     //
     times(0) = tinit;
     for (int n=0; n<rows; n++)
       for (int k=0; k<6; k++) ret(n, k, 0) = ps(n, k);
 
+    // Sign of h
+    int sgn = (0 < h) - (h < 0);
+
+    // Set the counters
     double tnow = tinit;
-    for (int s=1, cnt=1; s<numT; s++) {
+    int s = 0, cnt = 1;
+
+    // Do the integration using stride for output
+    while (s++ < numT) {
+      if ( (tfinal - tnow)*sgn < h*sgn) h = tfinal - tnow;
       std::tie(tnow, ps) = OneStep(tnow, h, ps, accel, bfe, F);
-      if (tnow >= H*cnt-h*1.0e-8) {
+      if (cnt < nout and s % stride == 0) {
 	times(cnt) = tnow;
 	for (int n=0; n<rows; n++)
 	  for (int k=0; k<6; k++) ret(n, k, cnt) = ps(n, k);
 	cnt += 1;
       }
     }
 
+    // Corrects round off at end point
+    //
     times(nout-1) = tnow;
     for (int n=0; n<rows; n++)
       for (int k=0; k<6; k++) ret(n, k, nout-1) = ps(n, k);
-    
+
     return {times, ret};
   }
 

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 closed 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
+            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<int>::max());
+	py::arg("nout")=0);
 }