Somewhat like ctypes
and similar libraries in Python, Julia has a built-in ccall
feature to call functions in external compiled (C ABI) libraries.
ccall(:printf, Cint, (Ptr{Uint8},), "Hello, world!")
13
Hello, world!
The format is ccall(function name, return type, argument types, arguments...)
.
You can also call functions in arbitrary shared libraries / DLLs:
mysin(x) = ccall((:sin,"libm"), Cdouble, (Cdouble,), x)
mysin (generic function with 1 method)
mysin(3.0) - sin(3.0)
0.0
mysin(3) # note that Julia automatically converts types as necessary
0.1411200080598672
Unlike Python, however, Julia's speed means that it is perfectly fine to call C functions operating on small data, like single numbers — you don't have to "vectorize" on the C side first, and you can instead vectorize on the Julia side.
@vectorize_1arg Real mysin
mysin (generic function with 4 methods)
mysin([1,2,3,4])
4-element Array{Float64,1}: 0.841471 0.909297 0.14112 -0.756802
code_native(mysin, (Float64,))
.section __TEXT,__text,regular,pure_instructions Filename: In[4] Source line: 1 push RBP mov RBP, RSP movabs RAX, 140735653238720 Source line: 1 call RAX pop RBP ret
Thanks to a package called PyCall
, Julia call arbitrary Python functions by calling directly into CPython's libpython
:
using PyCall
@pyimport math as pymath
pymath.cos(3)
-0.9899924966004454
Let's break this down. When you run pymath.sin(3)
, Julia:
3
into the corresponding Python object via PyObject(3)
.pymath.sin
via the libpython
routine PyObject_Call
.In terms of lower-level steps, it is doing:
three = PyObject(3) # calls PyInt_FromSsize_t in CPython library
PyObject 3
One slight annoyance is that Julia doesn't (yet) let you override .
, so foo.bar
in Python generally becomes foo[:bar]
in Julia (or foo["bar"]
if you want to leave the result as an unconverted Python object). This will change in a future Julia release.
mathmodule = pyimport("math") # calls PyImport_AddModule in CPython
sinfunc = mathmodule["sin"] # calls PyObject_GetAttrString
PyObject <built-in function sin>
returnval = pycall(sinfunc, PyObject, three) # calls PyObject_Call
PyObject 0.1411200080598672
convert(Float64, returnval) # if we know the type we want, we can specify it
# calls PyFloat_AsDouble in CPython
0.1411200080598672
convert(PyAny, returnval) # if we don't know the type we want, PyAny will detect it
0.1411200080598672
PyCall allows large arrays and dictionaries to be passed without making a copy.
For example, Julia arrays are wrapped by NumPy arrays with shared data.
A = rand(3,5)
3x5 Array{Float64,2}: 0.388691 0.485167 0.913292 0.373852 0.0255558 0.299324 0.881018 0.924125 0.349984 0.88512 0.575673 0.277566 0.0499966 0.0258964 0.157567
Apy = PyObject(A)
PyObject array([[ 0.38869052, 0.48516746, 0.91329182, 0.37385161, 0.02555581], [ 0.29932371, 0.88101823, 0.92412522, 0.34998437, 0.88511959], [ 0.57567336, 0.27756567, 0.04999658, 0.02589644, 0.15756664]])
A[1,1] = 17
Apy
PyObject array([[ 17. , 0.48516746, 0.91329182, 0.37385161, 0.02555581], [ 0.29932371, 0.88101823, 0.92412522, 0.34998437, 0.88511959], [ 0.57567336, 0.27756567, 0.04999658, 0.02589644, 0.15756664]])
@pyimport numpy as np
x = [-100, 39, 59, 55, 20]
np.irr(x)
0.2809484211599609
Warning: imported binding for transpose overwritten in module __anon__
By default, PyCall makes a copy of arrays when converting from Python back to Julia:
np.cumsum(x)
5-element Array{Int64,1}: -100 -61 -2 53 73
But we can specify a copy-free return if desired by calling the lower-level pycall
function and specifying the desired return type as PyArray
:
xsum = pycall(np.pymember("cumsum"), PyArray, x)
5-element Int64 PyArray: -100 -61 -2 53 73
The resulting NumPy array can be passed to other Julia routines without making a copy:
mean(xsum)
-7.4
There is also a PyVector
type that you can use to wrap Python lists without making a copy:
syspath = pyimport("sys")["path"]
PyObject ['//anaconda/lib/python27.zip', '//anaconda/lib/python2.7', '//anaconda/lib/python2.7/plat-darwin', '//anaconda/lib/python2.7/plat-mac', '//anaconda/lib/python2.7/plat-mac/lib-scriptpackages', '//anaconda/lib/python2.7/lib-tk', '//anaconda/lib/python2.7/lib-old', '//anaconda/lib/python2.7/lib-dynload', '//anaconda/lib/python2.7/site-packages', '//anaconda/lib/python2.7/site-packages/PIL', '//anaconda/lib/python2.7/site-packages/runipy-0.1.0-py2.7.egg', '//anaconda/lib/python2.7/site-packages/setuptools-3.6-py2.7.egg']
syspath_julia = PyVector(syspath)
12-element Any PyVector: "//anaconda/lib/python27.zip" "//anaconda/lib/python2.7" "//anaconda/lib/python2.7/plat-darwin" "//anaconda/lib/python2.7/plat-mac" "//anaconda/lib/python2.7/plat-mac/lib-scriptpackages" "//anaconda/lib/python2.7/lib-tk" "//anaconda/lib/python2.7/lib-old" "//anaconda/lib/python2.7/lib-dynload" "//anaconda/lib/python2.7/site-packages" "//anaconda/lib/python2.7/site-packages/PIL" "//anaconda/lib/python2.7/site-packages/runipy-0.1.0-py2.7.egg" "//anaconda/lib/python2.7/site-packages/setuptools-3.6-py2.7.egg"
There is also a PyDict
type that you can use to share a dictionary between Julia and Python.
d = PyDict()
PyDict{PyAny,PyAny} with 0 entries
d["hello"] = 7
d[23] = "goodbye"
d
PyDict{PyAny,PyAny} with 2 entries: "hello" => 7 23 => "goodbye"
For fun, we'll use pyeval
to pass d
as a local variable dict
to an arbitrary string of Python code that we want to evaluate, in this case a list comprehension in Python:
pyeval("[x for x in dict.keys()]", dict=d)
2-element Array{Any,1}: "hello" 23
Arbitrary Julia functions can be passed to Python. They get converted into callable Python objects of a custom class, whose __call__
method executes the Julia code:
foo(x) = x + 1
pyfoo = PyObject(foo)
PyObject <PyCall.jlwrap foo>
pycall(pyfoo, PyAny, 17)
18
This is extremely useful for calling functions for optimization, root-finding, etcetera, from SciPy. For example, let's solve a transcendental equation to find a root of $f(x) = \cos(x) - x$:
@pyimport scipy.optimize as so
function f(x)
println(" calling f($x)")
cos(x) - x
end
so.newton(f, 1.2)
calling f(1.2) calling f(1.2002199999999998) calling f(0.7664554749111869) calling f(0.7412167885608414) calling f(0.7390978176492645) calling f(0.7390851391787693)
0.7390851332151773
There is a bit of magic going on in passing Julia functions to Python. To define a new Python type from the CPython API, we create a C PyTypeObject
struct, and we need to stick a C function pointer into its tp_call
slot to give it a __call__
method.
A C function pointer is just the address of the compiled machine instructions, and since Julia has these instructions from its JIT compiler it can give you the address of the instructions using an intrinsic called cfunction
, e.g.
cfunction(f, Float64, (Float64,))
Ptr{Void} @0x000000011dd12210
This ability to get C function pointers from Julia functions is the key to calling any C API expecting a callback routine, not just Python. See the blog post: Passing Julia Callback Functions to C.
To get Matplotlib working in IJulia, we had to do a bit more work, similar to what IPython had to do to get its pylab
option working. For GUI windows, we had to implement the GUI event loop for the Python GUI toolkit(s) (PyQt, Gtk, wx) in Julia. For inline plots, we had to monkey-patch Matplotlib to intercept its drawing commands queue the figure for rendering as an image to be sent to the front-end. All of this is done by the PyPlot
Julia module:
using PyPlot
x = linspace(0,2π,1000)
plot(x, sin(3x + cos(5x)), "b--")
title("a funny plot")
PyObject <matplotlib.text.Text object at 0x1220d5a50>
y = linspace(0,2π,50)
surf(y, y, sin(y) .* cos(y)')
PyObject <mpl_toolkits.mplot3d.art3d.Poly3DCollection object at 0x12375a390>
Really, the whole Matplotlib API is available for use. It has everything you might want (in 2d, at least), if you dig long enough through the manuals:
xkcd()
fig = figure(figsize=(10,5))
ax = axes()
p = plot(x,sin(3x + cos(5x)))
ax[:set_xlim]([0.0,6])
annotate("A little\nsketchy",xy=[0.98,.001],arrowprops=["arrowstyle"=>"->"],xytext=[1.3,-0.3])
xticks([])
yticks([])
xlabel("TIME")
ylabel("PRODUCTIVITY")
title("An xkcd-style plot in Julia")
ax[:spines]["top"][:set_color]("none") # Remove the top axis boundary
ax[:spines]["right"][:set_color]("none") # Remove the right axis boundary
On Wednesday, Stéfan van der Walt showed some interesting demos using scikit-image. It is straightforward to try a couple of his examples from Julia:
@pyimport skimage.data as skdata
image = skdata.camera()
512x512 Array{Uint8,2}: 0x9c 0x9d 0xa0 0x9f 0x9e 0x9c … 0x97 0x96 0x97 0x98 0x98 0x98 0x9c 0x9d 0x9f 0x9e 0x9e 0x9c 0x97 0x96 0x97 0x98 0x98 0x98 0x9e 0x9d 0x9c 0x9c 0x9d 0x9d 0x99 0x98 0x99 0x98 0x98 0x98 0xa0 0x9d 0x9a 0x9a 0x9c 0x9d 0x9a 0x9a 0x9b 0x99 0x98 0x98 0x9e 0x9d 0x9c 0x9c 0x9d 0x9c 0x98 0x98 0x99 0x99 0x98 0x98 0x9c 0x9d 0x9f 0x9f 0x9f 0x9c … 0x97 0x97 0x98 0x99 0x99 0x99 0x9e 0x9d 0x9c 0x9c 0x9d 0x9c 0x98 0x98 0x99 0x99 0x99 0x99 0xa0 0x9d 0x9a 0x9a 0x9c 0x9d 0x9a 0x9a 0x9a 0x99 0x99 0x99 0x9e 0x9b 0x99 0x99 0x9b 0x9d 0x99 0x99 0x98 0x98 0x98 0x98 0x9c 0x9a 0x99 0x99 0x9b 0x9d 0x99 0x98 0x97 0x97 0x98 0x98 0x9b 0x9a 0x9a 0x99 0x9a 0x9c … 0x97 0x96 0x96 0x96 0x98 0x98 0x9b 0x9b 0x9b 0x9a 0x9b 0x9c 0x96 0x96 0x96 0x97 0x98 0x98 0x9b 0x9a 0x99 0x99 0x9c 0x9c 0x96 0x96 0x97 0x97 0x96 0x96 ⋮ ⋮ ⋱ ⋮ 0x82 0x7d 0x77 0x77 0x78 0x6d … 0x66 0x5f 0x5a 0x5d 0x61 0x62 0x71 0x76 0x7b 0x7b 0x77 0x63 0x59 0x56 0x54 0x5c 0x64 0x66 0x65 0x6e 0x76 0x76 0x72 0x60 0x5b 0x61 0x67 0x67 0x66 0x67 0x62 0x6a 0x72 0x72 0x6f 0x64 0x66 0x72 0x7d 0x75 0x6b 0x69 0x6f 0x70 0x71 0x73 0x75 0x70 0x75 0x7b 0x7e 0x78 0x72 0x70 0x7e 0x79 0x74 0x77 0x7c 0x7b … 0x84 0x82 0x7d 0x7a 0x78 0x77 0x79 0x7b 0x7d 0x7d 0x7b 0x73 0x86 0x86 0x84 0x7c 0x75 0x73 0x72 0x7c 0x85 0x82 0x7a 0x6e 0x86 0x88 0x8a 0x7d 0x71 0x6e 0x74 0x7b 0x82 0x80 0x7e 0x81 0x85 0x86 0x86 0x7b 0x70 0x6d 0x79 0x7b 0x7e 0x7e 0x83 0x96 0x85 0x85 0x82 0x79 0x71 0x6f 0x79 0x7b 0x7e 0x7e 0x83 0x98 … 0x85 0x85 0x82 0x79 0x71 0x6f 0x79 0x7b 0x7e 0x7e 0x82 0x96 0x85 0x85 0x82 0x79 0x71 0x6f
@pyimport skimage.exposure as exposure
subplot(1,3,1)
imshow(image, cmap="jet")
subplot(1,3,2)
imshow(image, cmap="gray")
subplot(1,3,3)
imshow(exposure.equalize_adapthist(image), cmap="gray")
PyObject <matplotlib.image.AxesImage object at 0x12426ae90>
@pyimport skimage.filter as filters
threshold = filters.threshold_otsu(image)
println("Otsu threshold = $threshold")
imshow(image .> threshold, cmap="gray")
Otsu threshold = 87
PyObject <matplotlib.image.AxesImage object at 0x125fc6b10>