A high-level design for exposing the symplectic integrator in boost

.AA_README_12690_The_case_for_Symplectic_Solvers_in_ODE_Solver

This gist contains code for this scipy issue
https://github.com/scipy/scipy/issues/12690#issuecomment-1214276702

.vscode-launch.json

{
    "configurations": [
        {
            "name": "C/C++: g++-12 build and debug active file",
            "type": "cppdbg",
            "request": "launch",
            "program": "${fileDirname}/${fileBasenameNoExtension}",
            "args": [],
            "stopAtEntry": false,
            "cwd": "${fileDirname}",
            "environment": [],
            "externalConsole": false,
            "MIMode": "lldb",
            "preLaunchTask": "C/C++: g++-12 build active file"
        }
    ],
    "version": "2.0.0"
}

.vscode-tasks.json

{
    "tasks": [
        {
            "type": "cppbuild",
            "label": "C/C++: g++-12 build active file",
            "command": "/usr/local/Cellar/gcc/12.1.0/bin/g++-12",
            "args": [
                "-fdiagnostics-color=always",
                "-g",
                "${file}",
                "-I", 
                "/usr/local/anaconda3/envs/py310/lib/python3.10/site-packages/scipy/_lib/",
                "-Wno-deprecated-declarations",
                "-o",
                "${fileDirname}/${fileBasenameNoExtension}"
            ],
            "options": {
                "cwd": "${fileDirname}"
            },
            "problemMatcher": [
                "$gcc"
            ],
            "group": {
                "kind": "build",
                "isDefault": true
            },
            "detail": "compiler: /usr/local/Cellar/gcc/12.1.0/bin/g++-12"
        }
    ],
    "version": "2.0.0"
}

_symplectic.hpp

#include <iostream>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>
// #include <boost/numeric/odeint/stepper/symplectic_rkn_sb3a_mclachlan.hpp>
using namespace boost::numeric::odeint;
typedef std::vector< double > pq_t;
typedef pq_t (*f_pq2dpq_t)(pq_t);

typedef std::pair<pq_t, pq_t> state_t;
typedef symplectic_rkn_sb3a_mclachlan< pq_t > stepper_type;
typedef std::tuple<
    std::vector<pq_t>, 
    std::vector<pq_t>, 
    std::vector<double>> result;

struct D_wrapper {
    const f_pq2dpq_t f_pq2dpq;

    D_wrapper (const f_pq2dpq_t f) : f_pq2dpq(f) {}

    void operator()(const pq_t &pq, pq_t &dpqdt) const {
        pq_t tmp = f_pq2dpq(pq);
        for (size_t i = 0; i < pq.size(); ++i) dpqdt[i] = tmp[i];
    }
};

struct Observer
{
    std::vector< pq_t >& m_p;
    std::vector< pq_t >& m_q;
    std::vector< double >& m_times;

    Observer(std::vector< pq_t > &p, std::vector< pq_t > &q, std::vector< double > &times)
    : m_p(p) , m_q(q), m_times(times) { }

    void operator()(const state_t &x , double t) const
    {
        std::cout << x.first[0] << ',' << x.second[0] << ',' << t << '\n';
        m_q.push_back( x.first );
        m_p.push_back( x.second );
        m_times.push_back( t );
    }
};

result solve_ivp_symplectic (
    double dt, double t0, double tn,
    f_pq2dpq_t q2dp, f_pq2dpq_t p2dq,
    pq_t q0, pq_t p0) 
{
    std::vector<pq_t> m_p;
    std::vector<pq_t> m_q;
    std::vector<double> m_t;
    Observer observer(m_p, m_q, m_t);

    integrate_const(
        stepper_type() ,
        std::make_pair( D_wrapper( q2dp ) , D_wrapper( p2dq ) ) ,
        std::make_pair( boost::ref( q0 ) , boost::ref( p0 ) ) ,
        t0 , 
        tn , 
        dt ,  
        observer);
    result res = std::make_tuple(m_q, m_p, m_t);
    return res;
}

setup.sh

get_boost_dir () {
    dirname $(python -c 'from scipy._lib._boost_utils import _boost_dir ; print(_boost_dir())')
}
first_time () {
    mkdir symplectic_ode_boost
    cd symplectic_ode_boost
    conda create -n py310 python=3.10
    conda activate py310
    which python
    pip install scipy
    python --version # Python 3.10.4
    python -c 'import scipy; print(scipy.__version__)' # 1.9.0
    git clone --depth 1 https://github.com/scipy/boost-headers-only
    ### These patches seem to already be applied
    # cd boost-headers-only
    # for f in patches/*.patch ; do git apply $f ; done
    boost_dir=$(get_boost_dir)
    ln -s $PWD/boost-headers-only/boost $boost_dir
}
# cd symplectic_ode_boost
conda activate py310
boost_dir=$(get_boost_dir)
echo $boost_dir

g++-12  -I $boost_dir -Wno-deprecated-declarations  symplectic_demo.cpp
./a.out > symplectic.log
python symplectic_demo.py

symplectic_demo.cpp

#include "_symplectic.hpp"
int main(){
    double dt = 1;
    double t0 = 0;
    double tn = 1000;
    pq_t q0 {2,};
    pq_t p0 {-1,};
    f_pq2dpq_t q2dp = [](pq_t q){
        pq_t tmp(q.size());
        for (size_t i = 0; i < q.size(); ++i) tmp[i] = -q[i];
        return tmp;};
    f_pq2dpq_t p2dq = [](pq_t p){
        pq_t tmp(p.size()); 
        for (size_t i = 0; i < p.size(); ++i) tmp[i] = p[i];
        return tmp;};
    result res = solve_ivp_symplectic(
        dt, t0, tn, q2dp, p2dq, q0, p0
    );
}

symplectic_demo.py

""" Consider a system with 2 scalar generalized coordinates that
interact as follows

dx1dt = x2
dx2dt = -x1

The solution is
x1(t) = r cos(a + t)
x2(t) = r sin(a + t)

Where r and a are determined by the initial value. 
Now Let's use a symplectic integrator to solve this sytem.
"""
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from matplotlib import animation, rc

def dx(t, x):
    return np.array([x[1], -x[0]])

x0 = np.array([2, -1])
# True solution: r = sqrt(5)
# a = acos(2/sqrt(5)) radians
t_span = (0, 1000)
t_eval = np.linspace(t_span[0], t_span[1], 1000)
# method in 
for method in 'RK45 RK23 DOP853 LSODA'.split():
    solution = solve_ivp(dx, t_span, x0, t_eval=t_eval, method=method)
    t = solution.t
    x = solution.y
    plt.plot(x[0], x[1])
    plt.title(method)
    plt.savefig(f'/tmp/{method}.pdf')
    plt.close()

# from _symplectic import solve_ivp_symplectic
# solution = solve_ivp_symplectic(dx, t_span, x0)
# t = solution.t
# x = solution.y
method = 'symplectic'
x = np.array([[float(e) for e in row.strip().split(',')] 
    for row in open('symplectic.log')]).T
plt.plot(x[0], x[1])
plt.title(method)
plt.savefig(f'/tmp/{method}.pdf')
plt.close()

vscode-c_cpp_properties.json

{
    "configurations": [
        {
            "name": "Mac",
            "includePath": [
                "${workspaceFolder}/**",
                "/usr/local/anaconda3/envs/py310/lib/python3.10/site-packages/scipy/_lib/**"
            ],
            "defines": [],
            "macFrameworkPath": [
                "/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/System/Library/Frameworks"
            ],
            "compilerPath": "/usr/local/Cellar/gcc/12.1.0/bin/g++-12",
            "cStandard": "c17",
            "cppStandard": "c++17",
            "intelliSenseMode": "macos-gcc-x64"
        }
    ],
    "version": 4
}