**Unformatted text preview: **Crossing Language Barriers with , SciPy, and thon Steven G. Johnson MIT Applied MathemaBcs Where I’m coming from… [ google “Steven Johnson MIT” ] ComputaBonal soPware you may know… … mainly C/C++ libraries & soPware … … oPen with Python interfaces … (& Matlab & Scheme & …) jdj.mit.edu/nlopt Nanophotonics jdj.mit.edu/meep & other EM simulators… erf(z) (and erfc, erﬁ, …) in SciPy 0.12+ jdj.mit.edu/book Confession: I’ve used Python’s internal C API more than I’ve coded in Python… A new programming language? Viral Shah Jeﬀ Bezanson Alan Edelman julialang.org Stefan Karpinski [ 17+ developers with 100+ commits ] [begun 2009, “0.1” in 2013, ~20k commits] First reacBon: You’re doomed. [ usual fate of all new languages ] … subsequently: … probably doomed … s7ll might be doomed but, in the meanBme, I’m having fun with it… … and it solves a real problem with technical compuBng in high-‐level languages. The “Two-‐Language” Problem Want a high-‐level language that you can work with interacBvely = easy development, prototyping, exploraBon ⇒ dynamically typed language Plenty to choose from: Python, Matlab / Octave, R, Scilab, … (& some of us even like Scheme / Guile) Historically, can’t write performance-‐criBcal code (“inner loops”) in these languages… have to switch to C/Fortran/… (staBc). [ e.g. SciPy git master is ~70% C/C++/Fortran] Workable, but Python → Python+C = a huge jump in complexity. Just vectorize your code? = rely on mature external libraries, operaBng on large blocks of data, for performance-‐criBcal code Good advice! But… • Someone has to write those libraries. • Eventually that person may be you. — some problems are impossible or just very awkward to vectorize. Dynamic languages don’t have to be slow. Lots of progress in JIT compilers, driven by web applicaBons. & excellent free/open-‐source JIT via LLVM. Javascript in modern browsers achieves C-‐like speed. Many other eﬀorts to speed up dynamic languages, e.g. PyPy, Numba / Cython (really 2nd lower-‐level language embedded in Python). What if a dynamic language were designed for JIT from the beginning, with the goal of being as high-‐
level as possible while staying within 2× C speed? (and it’s easier to call SciPy from Julia than from PyPy) Today A brief introducBon to Julia, its key features, and how it gets performance. How Julia leverages Python and IPython to lessen the “infrastructure problem” of new languages Cme permiDng: How tools can ﬂow in the other direcBon too… Julia [ julialang.org ] Dynamically typed MulBple dispatch: a generalizaBon of OO Metaprogramming (code that writes code) Direct calling of C and Python funcBons CorouBnes, asychronous I/O Most of Julia (70%+) Designed for Unicode is wriuen in Julia. Distributed-‐memory parallelism User-‐deﬁned types == builBn types … extensible promoBon and conversion rules, etc. • Large built-‐in library: regex, linear algebra, special funcBons, integraBon, etcetera… • git-‐based package manager •
•
•
•
•
•
•
• goto live Julia REPL demo… This command line is so 1990… [e.g. GNU Readline, 1989] vs. Come back in 5 years once you implement something more modern… [ ipython.org ] (roughly) How IPython Notebooks Work programmer enter Python code ZeroMQ HTTP code Web browser: Notebook Display/InteracBon HTM
L
+ JS IPython Web Server code mulB-‐ media results IPython kernel How IJulia Notebooks Work programmer enter Julia code HTTP code Web browser: Notebook Display/InteracBon HTM
L
+ JS ZeroMQ IPython Jupyter Web Server code mulB-‐ media results IJulia kernel goto live IJulia notebook demo… Why is Julia fast? Why is Julia fast? Julia performance on syntheBc benchmarks [ loops, recursion, etc., implemented in most straighxorward style ] What about real problems, compared to highly opBmized code? Special FuncBons in Julia Special funcBons s(x): classic case that cannot be vectorized well … switch between various polynomials depending on x Many of Julia’s special funcBons come from the usual C/Fortran libraries, but some are wriuen in pure Julia code. Pure Julia erﬁnv(x) [ = erf–1(x) ] 3–4× faster than Matlab’s and 23× faster than SciPy’s (Fortran Cephes). Pure Julia polygamma(m, z) [ = (m+1)th derivaBve of the ln Γ funcBon ] ~ 2× faster than SciPy’s (C/Fortran) for real z … and unlike SciPy’s, same code supports complex argument z Julia code can actually be faster than typical “opBmized” C/Fortran code, by using techniques [metaprogramming/
code generaBon] that are hard in a low-‐level language. Pure-‐Julia FFT performance double-precision complex, 1d transforms
powers of two
14000
13000 (FFTW, MKL: “unfair” factor of ~2 from manual SIMD) 12000
11000
10000
speed (mflops) intel-mkl-dfti in-place
intel-mkl-dfti out-of-place
fftw3 out-of-place
fftw3 in-place
fftw3-no-simd out-of-place
fftw3-no-simd in-place
dfftpack
emayer
julia
bloodworth
cross
cwplib
esrfft 9000
8000
7000
6000
5000 FFTW Julia 4000 already comparable to FFTPACK [ probably some tweaks to inlining will make it beuer ] w/o SI
M D 3000
2000
1000
262144 131072 65536 32768 16384 8192 4096 2048 1024 512 256 128 64 32 16 8 4 2 0 FFTW 1.0-‐like code generaBon + recursion in Julia ~ 1/3 lines of code compared to FFTPACK, more funcBonality Why is Julia fast? Why is Julia fast? Why can Julia be fast? (You can write slow code in any language, of course.) … and couldn’t Python do the same thing? Type Inference To generate fast code for a funcBon f(x,y), the compiler needs to be able to infer the types of variables in f, map them to hardware types (registers) where possible, and call specialized code paths for those types (e.g. you want to inline +, but this depends on types). At compile-‐Bme, the compiler generally only knows types of x,y, not values, and it needs to be able to cascade this informaBon to infer types throughout f and in any funcBons called by f. Julia and its standard library are designed so type inference is possible for code following straighxorward rules. … someBmes this requires subtle choices that would be painful to retroﬁt onto an exisBng language. Type Stability Key to predictable, understandable type inference: • the type of funcBon’s return value should only depend on the types of its arguments A counter-‐example in Matlab and GNU Octave: sqrt(1) == 1.0 — real ﬂoaBng-‐point sqrt(–1) == 0.0+1.0i — complex ﬂoaBng-‐point Hence, any non-‐vector code that calls sqrt(x) in Matlab cannot be compiled to fast code even if x is known to be real scalar — anything “touched” by the sqrt(x) is “poisoned” with an unknown type — unless the compiler knows x ≥ 0. Beuer to throw an excepBon for sqrt(–1), requiring sqrt(-‐1+0i). Type Stability Key to predictable, understandable type-‐inference: • the type of funcBon’s return value should only depend on the types of its arguments Common counter-‐examples in Python Typical idiom: foo(x) returns y, or None if [excepBonal condiBon] [e.g. numpy.ma.notmasked_edges, scipy.constants.ﬁnd, …] Beuer: Throw an excepBon. Type Stability Key to predictable, understandable type-‐inference: • the type of funcBon’s return value should only depend on the types of its arguments A counter-‐example in Python integer arithmeBc Integer arithmeBc in Python automaBcally uses bignums to prevent overﬂow. Unless the compiler can detect that overﬂow is impossible [which may be detectable someBmes!], integers can’t be compiled to integer registers & hw arithmeBc. Julia tradeoﬀ: default integers are 64-‐bit, overﬂow possible … use explicit BigInt type if you are doing number theory. goto live IJulia notebook demo… Julia: fun, fast, and you don’t lose your Python stuﬀ. New languages are always a risk… …but maybe not doomed? Acknowledgements julialang.org Julia core team: Jeﬀ Bezanson (MIT) Stefan Karpinski (MIT) Viral Shah …(17+ developers with 100+ commits)… Prof. Alan Edelman (MIT) & Shashi Gowda (GSoC) ipython.org Fernando Perez Mauhias Bussonier Min RK & Jake Bolewski (pyjulia) ...

View
Full Document

- Spring '16