Moved from GSL Runge-Kutta to Boost Bulirsch-Stoer

Makefile

@@ -1,5 +1,6 @@
 OPTIMIZATION ?= 3
 CXXFLAGS     += -O$(OPTIMIZATION)
+INC          += -I/home/meiron/boost_1_72_0
 LIB          += -lgsl
 
 EXECUTABLE   ?= main

main.cpp

@@ -1,8 +1,7 @@
 #include <algorithm>
+#include <array>
+#include <boost/numeric/odeint.hpp>
 #include <fstream>
-#include <gsl/gsl_errno.h>
-#include <gsl/gsl_math.h>
-#include <gsl/gsl_odeiv2.h>
 #include <gsl/gsl_spline.h>
 #include <iostream>
 #include <numeric>
@@ -13,12 +12,17 @@
 #include "loadtxt.h"
 
 
-extern "C" const int gsl_success() { return GSL_SUCCESS; } // It's zero, but just for clarity sake.
 
 // Our units are {kiloparsec, solar mass, gigayear}
 constexpr double G = 4.498317481097514e-06;
 
 class Interp {
+    // This is a wrapper around GSL spline interpolation. I tried to use
+    // boost::math::interpolators but as of version 1.72 none were suitable:
+    // barycentric_rational is the one suitalbe for non-uniform sampling but it
+    // is very slow. I also tried to resample the data uniformly using
+    // barycentric rational interpolation and then using cardinal_cubic_b_spline
+    // on the result, but was still slower than GSL.
 public:
     Interp(std::vector<double>& x, std::vector<double>& y)
     {
@@ -63,12 +67,12 @@ public:
         auto& t_data      = data[0];
         auto& halo_m_data = data[1];
         auto& halo_b_data = data[2];
-        std::transform(t_data.begin(), t_data.end(), t_data.begin(), [](const double& x){ return x-2.145; });
+        std::transform(begin(t_data), end(t_data), begin(t_data), [t0=t_data[0]](const double& x){ return x-t0; });
         std::transform(halo_b_data.begin(), halo_b_data.end(), halo_b_data.begin(), [](const double& x){ return x*0.7664209365408798; });
         interp_halo_m = Interp(t_data, halo_m_data);
         interp_halo_b = Interp(t_data, halo_b_data);
     }
-    int func(double t, const double y[], double f[], void *params)
+    void func(const std::array<double, 6> &y, std::array<double, 6> &f, const double t)
     {
         double halo_m = interp_halo_m(t);
         double halo_b = interp_halo_b(t);
@@ -76,56 +80,86 @@ public:
         f[0] = y[3]; // vx -> x'
         f[1] = y[4]; // vy -> y'
         f[2] = y[5]; // vz -> z'
-        plummer.calc_acceleration(y, f+3); // a -> v'
+        plummer.calc_acceleration(y.data(), f.data()+3); // a -> v'
-        return GSL_SUCCESS;
     }
 
 private:
     Interp interp_halo_m;
     Interp interp_halo_b;
 } galaxy("file.dat");
-// Not very nice to have it as a global variable but GSL will have problem otherwise.
-
-int jac(double t, const double y[], double *dfdy, double dfdt[], void *params) { return GSL_SUCCESS; }
-
-inline int func(double t, const double y[], double f[], void *params)
-{
-    return galaxy.func(t, y, f, params);
-}
 
 extern "C"
-int integrate(const double y0[], const double t_max, const double step_size, double y[])
+int integrate(const double y0[], const double t_max, const double stride_size, double y[])
 {
-    double t = 0;
+    using namespace boost::numeric::odeint;
-    constexpr double h = 1./4096.;
+    using Coordinates = std::array<double, 6>;
-
+    auto stepper = bulirsch_stoer<Coordinates>(1E-7, 0);
-    if (step_size/h - (int)(step_size/h) != 0) throw std::runtime_error("step_size must be a multiple of h");
+    auto function_wrapper = [](const Coordinates &x, Coordinates &dxdt, const double t) { return galaxy.func(x, dxdt, t); };
-    constexpr double epsabs = 1e-7;
+    const int stride_count = t_max / stride_size;
-    constexpr double epsrel = 0;
+    Coordinates y_current;
-    const gsl_odeiv2_step_type *T = gsl_odeiv2_step_rk8pd;
+    std::copy(y0, y0+6, begin(y_current));
-    gsl_odeiv2_step *s    = gsl_odeiv2_step_alloc(T, 6);
-    gsl_odeiv2_evolve *e  = gsl_odeiv2_evolve_alloc(6);
-    gsl_odeiv2_control *c = gsl_odeiv2_control_y_new(epsabs, 0);
-    gsl_odeiv2_system sys = {func, jac, 6, nullptr};
-    gsl_odeiv2_driver *d  = gsl_odeiv2_driver_alloc_y_new(&sys, T, h, epsabs, epsrel);
-
-    int step = 0;
-    const int step_max = t_max / step_size;
     std::copy(y0, y0+6, y);
-    for (int step=0; step<step_max; step++) {
+    double t = 0;
-        std::copy(y+step*6, y+(step+1)*6, y+(step+1)*6);
+    const double h = 1./4096.;
-        int status = gsl_odeiv2_driver_apply(d, &t, (step+1)*step_size, y+(step+1)*6);
+    for (int i=0; i<stride_count; i++) {
-        if (status != GSL_SUCCESS) return status;
+        integrate_adaptive(stepper, function_wrapper, y_current, t, t+stride_size, h);
+        // NOTE h here is just a recommended initial step size for the stepper,
+        // the actual step is adapted as needed. Since the integration is
+        // interrupted here in order for the data to be saved, the result
+        // somewhat depends on stride_size.
+        std::copy(begin(y_current), end(y_current), y+(i+1)*6);
+        t += stride_size;
     }
-    return GSL_SUCCESS;
 }
 
+double circular_energy(double(*pot)(double), double L)
+{
+    auto effective_potential = [&](double r) { return pot(r)+0.5*(L*L)/(r*r); };
+    return pot(3);
+}
+
+/*def circular_energy(Phi, L, PhiPrime=None, InitialGuess=1.0, GetRadius=False):
+    EffectivePotential = lambda r: Phi(r) + 0.5*(L/r)**2
+    EffectivePotentialPrime = None
+    if not PhiPrime is None: EffectivePotentialPrime = lambda r: PhiPrime(r) - L**2/r**3
+    Minimization = scipy.optimize.minimize(EffectivePotential, InitialGuess, method='BFGS', jac=EffectivePotentialPrime, tol=1.0E-08)
+    # Minimization = scipy.optimize.minimize(EffectivePotential, InitialGuess, method='L-BFGS-B', bounds=[(0,None)], jac=EffectivePotentialPrime)
+    # There is some risk that the negative value of r will be found as the minimum. To avoid this we can either set boundary conditions which might complicate the the numerical procedure, or just accept this and return the absolute value, which we can do here.
+    if GetRadius: return Minimization.fun, abs(Minimization.x[0])
+    return Minimization.fun*/
+
 int main()
 {
-    std::cout << "bye" << std::endl;
+    std::cout << "Hi" << std::endl;
+
+/*std::vector<double> x(500);
+std::vector<double> y(500);
+for (int i=0; i<500; i++) {
+    x[i] = i;
+    y[i] = 2*i;
+}
+// populate x, y, then:
+boost::math::barycentric_rational<double> interpolant(std::move(x), std::move(y));
+
+
+std::cout << interpolant(5.5) << std::endl;*/
+
+
+
     double y[12];
     double y0[] = {80,0,0,0,80,0};
-    for (int i=0; i<30000; i++)
+    for (int i=0; i<8000; i++)
         integrate(y0, 10, 10, y);
+
+/*double y0[] = {80,0,0,0,80,0};
+double y[12];
+integrate(y0, 10, 10, y);
+for (int i=0; i<12; i++) std::cout << y[i] << std::endl;*/
+
+    //std::vector<double> w = {8, 0.12, 0.04, -0.08, 200, -0.001};
+
+
+
+    std::cout << "Bye" << std::endl;
     return 0;
 }

plot.py

@@ -1,10 +1,13 @@
+#pylint: disable=W0401,W0614,W0622
+#All pylint codes: http://pylint.pycqa.org/en/latest/technical_reference/features.html
+
 import numpy as np
 import ctypes, subprocess
 
 recompile   = True
 
 if recompile:
-    p = subprocess.Popen('g++ -shared -o libmain.so -fPIC loadtxt.cpp main.cpp -lgsl'.split())
+    p = subprocess.Popen('g++ -I/home/meiron/boost_1_72_0 -shared -o libmain.so -fPIC loadtxt.cpp main.cpp -lgsl'.split())
     p.wait()
     if p.returncode != 0: raise RuntimeError(p.returncode)
 
@@ -18,12 +21,21 @@ def integrate(y0, t_max, step_size=None):
     y = np.array(y).reshape(size,6)
     return np.array(y), status
 
-gsl_success = libmain.gsl_success()
+def integrate_boost(y0, t_max, step_size=None):
+    y0 = (ctypes.c_double*6)(*y0)
+    if step_size is None: step_size = t_max
+    size = int(t_max // step_size) + 1
+    y = (ctypes.c_double*size*6)()
+    libmain.integrate_boost(y0, ctypes.c_double(t_max), ctypes.c_double(step_size), y)
+    y = np.array(y).reshape(size,6)
+    return np.array(y), 0
+
+#gsl_success = libmain.gsl_success()
 
 from pylab import *
 t_max = 10
 ic = [80,0,0,0,80,0]
-res = integrate(ic, t_max, step_size=1/4096)
+res = integrate(ic, t_max, step_size=32/4096)
 x, y, z, vx, vy, vz = res[0].T
 zzz = x[-1]
 plot(x,y)
@@ -32,5 +44,5 @@ res = integrate(ic, t_max)
 x, y, z, vx, vy, vz = res[0].T
 plot(x,y,'o')
 print(zzz - x[-1])
-# gca().set_aspect('equal')
+gca().set_aspect('equal')
 show()