From c2051e9596acd1a7d820261c73f21fe037422800 Mon Sep 17 00:00:00 2001
From: Yohai Meiron <sahmes@gmail.com>
Date: Tue, 21 Apr 2020 20:05:56 -0400
Subject: [PATCH] Moved from GSL Runge-Kutta to Boost Bulirsch-Stoer

---
 Makefile |   1 +
 main.cpp | 112 ++++++++++++++++++++++++++++++++++++-------------------
 plot.py  |  20 ++++++++--
 3 files changed, 90 insertions(+), 43 deletions(-)

diff --git a/Makefile b/Makefile
index 7e5a87f..02addc9 100644
--- a/Makefile
+++ b/Makefile
@@ -1,5 +1,6 @@
 OPTIMIZATION ?= 3
 CXXFLAGS     += -O$(OPTIMIZATION)
+INC          += -I/home/meiron/boost_1_72_0
 LIB          += -lgsl
 
 EXECUTABLE   ?= main
diff --git a/main.cpp b/main.cpp
index e2e5624..b7b1f85 100644
--- a/main.cpp
+++ b/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'
-        return GSL_SUCCESS;
+        plummer.calc_acceleration(y.data(), f.data()+3); // a -> v'
     }
 
 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;
-    constexpr double h = 1./4096.;
-
-    if (step_size/h - (int)(step_size/h) != 0) throw std::runtime_error("step_size must be a multiple of h");
-    constexpr double epsabs = 1e-7;
-    constexpr double epsrel = 0;
-    const gsl_odeiv2_step_type *T = gsl_odeiv2_step_rk8pd;
-    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;
+    using namespace boost::numeric::odeint;
+    using Coordinates = std::array<double, 6>;
+    auto stepper = bulirsch_stoer<Coordinates>(1E-7, 0);
+    auto function_wrapper = [](const Coordinates &x, Coordinates &dxdt, const double t) { return galaxy.func(x, dxdt, t); };
+    const int stride_count = t_max / stride_size;
+    Coordinates y_current;
+    std::copy(y0, y0+6, begin(y_current));
     std::copy(y0, y0+6, y);
-    for (int step=0; step<step_max; step++) {
-        std::copy(y+step*6, y+(step+1)*6, y+(step+1)*6);
-        int status = gsl_odeiv2_driver_apply(d, &t, (step+1)*step_size, y+(step+1)*6);
-        if (status != GSL_SUCCESS) return status;
+    double t = 0;
+    const double h = 1./4096.;
+    for (int i=0; i<stride_count; i++) {
+        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;
 }
\ No newline at end of file
diff --git a/plot.py b/plot.py
index dcfa01d..d805035 100644
--- a/plot.py
+++ b/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()
\ No newline at end of file