Last active
August 19, 2022 19:06
-
-
Save pushpendre/595f1531c7995adfa161270377e14aed to your computer and use it in GitHub Desktop.
A high-level design for exposing the symplectic integrator in boost
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This gist contains code for this scipy issue | |
https://github.com/scipy/scipy/issues/12690#issuecomment-1214276702 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"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" | |
} | |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"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" | |
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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 > ×) | |
: 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; | |
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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 | |
); | |
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
""" 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() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"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 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment