JupyterのCell MagicでFortranやC++やCythonを使う

jupyter notebookではpython以外にも様々な言語が使える。

スポンサーリンク

fortran

ipythonでfortranを使えるようにするにはまず以下のコマンドを実行する。

!pip3 install -U fortran-magic

次に下記のコマンドを実行してfortranをjupyter cellにロードする。

%load_ext fortranmagic
/home/aho/.pyenv/versions/3.6.5/envs/py365/lib/python3.6/site-packages/fortranmagic.py:147: UserWarning: get_ipython_cache_dir has moved to the IPython.paths module since IPython 4.0.
  self._lib_dir = os.path.join(get_ipython_cache_dir(), 'fortran')

cell magic %%fortranでfortranが使えるようになる。

%%fortran
subroutine compute_fortran(x, y, z)
    real, intent(in) :: x(:), y(:)
    real, intent(out) :: z(size(x, 1))

    z = sin(x + y)

end subroutine compute_fortran
compute_fortran([1, 2, 3], [4, 5, 6])
array([-0.9589243,  0.6569866,  0.4121185], dtype=float32)
%%fortran --link lapack

subroutine solve(A, b, x, n)
    ! solve the matrix equation A*x=b using LAPACK
    implicit none

    real*8, dimension(n,n), intent(in) :: A
    real*8, dimension(n), intent(in) :: b
    real*8, dimension(n), intent(out) :: x

    integer :: pivot(n), ok

    integer, intent(in) :: n
    x = b

    ! find the solution using the LAPACK routine SGESV
    call DGESV(n, 1, A, n, pivot, x, n, ok)

end subroutine
A = np.array([[1, 2.5], [-3, 4]])
b = np.array([1, 2.5])

solve(A, b)
array([-0.19565217,  0.47826087])

cython

今度は、cell magic %%cythonを使えるようにする。

%load_ext Cython
%%cython
def power_of_two(int x):
    return 2.0 ** x
power_of_two(256) #2の256乗
1.157920892373162e+77
%%cython --compile-args=-fopenmp --link-args=-fopenmp --force
cimport cython
cimport openmp

import cython.parallel as cp
from cython.parallel import parallel, prange

import numpy as np
cimport numpy as np

ctypedef np.int32_t int32_t

cdef extern from "unistd.h" nogil:
    unsigned int sleep(unsigned int seconds)

@cython.boundscheck(False)
cpdef double[::1] ompfunc(int32_t size, int32_t num_threads):
    cdef int32_t idx
    cdef double[::1] result = np.empty(size, dtype=np.float64)
    with nogil, cp.parallel(num_threads=num_threads):
        for idx in prange(size, schedule='dynamic'):
            result[idx] = cp.threadid()
            sleep(1)
    return result
np.asarray(ompfunc(10, 10))
array([0., 9., 1., 2., 3., 4., 5., 6., 7., 8.])
%timeit ompfunc(10, 1)
%timeit ompfunc(10, 2)
%timeit ompfunc(10, 3)
%timeit ompfunc(10, 4)
%timeit ompfunc(10, 5)
%timeit ompfunc(10, 6)
%timeit ompfunc(10, 7)
%timeit ompfunc(10, 8)
%timeit ompfunc(10, 9)
%timeit ompfunc(10, 10)
10.1 s ± 19.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
5.07 s ± 47 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
4 s ± 3.07 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
3 s ± 2.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
2 s ± 219 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
2 s ± 2.56 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
2 s ± 217 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
2 s ± 325 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
2 s ± 257 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
1 s ± 138 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

C++

cythonを使えばC++プログラムをpythonモジュールに変換できる。

%%file fib.hpp

double fib(int n);
Overwriting fib.hpp
%%file fib.pxd

cdef extern from "fib.hpp":
    double fib(int n)
Writing fib.pxd
%%file fib.cpp

double fib(int n) {
    double a = 0, b = 1;
    for (int i=0; i<n; i++) {
        double tmp = b;
        b = a;
        a += tmp;
     }
    return a;
}
Overwriting fib.cpp
%%file setup.py
from distutils.core import setup, Extension
from Cython.Build import cythonize

ext = Extension("fib2cpp",
              sources=["fib2cpp.pyx", "fib.cpp"],
              language="c++",)

setup(name = "cython_fibcpp",
      ext_modules = cythonize(ext))
%%file fib2cpp.pyx

cimport fib

def fib(n):
    return fib.fib(n)
! python setup.py build_ext --inplace
import fib2cpp
fib2cpp.fib(100)
3.54224848179262e+20

cythonizeのtemplateはこれを使いまわせるので、自分なりの関数を自由に作れる。

参考サイトhttps://github.com/
参考サイトhttps://people.duke.edu/~ccc14/sta-663/FromCompiledToPython.html
参考サイトhttps://blogs.baruch.cuny.edu/cis3100/?p=40