Timing relative to C (lower is better). http://julialang.org
"""approximate π^2/6 with the first n terms in an infinite series."""
function pi_sum(n)
sum = 0.0
for k in 1:n
sum += 1.0 / (k*k)
end
return sum
end
pi_sum (generic function with 1 method)
sqrt(6 * pi_sum(10000000.7))
3.1415925580959025
#define
in C but more powerful)One of the best things about coming to Julia from Python is that the languages are quite similar in semantics. Specifically, the way variables are assigned and passed to functions is identical. While you have to remember the surface syntax differences, you don't have to re-learn how to think about your code.
a = [1.0, 2.0, 3.0, 4.0] # some array
4-element Array{Float64,1}: 1.0 2.0 3.0 4.0
b = a # assign the name "b" to the same array that 'a' is pointing to.
b[1] = 5.0 # modify the first element in that array
b
4-element Array{Float64,1}: 5.0 2.0 3.0 4.0
a # change is reflected in a
4-element Array{Float64,1}: 5.0 2.0 3.0 4.0
# define a function that modifies an array
function double!(x)
for i=1:length(x)
x[i] *= 2.0
end
end
double! (generic function with 1 method)
a = [1.0, 2.0, 3.0, 4.0]
4-element Array{Float64,1}: 1.0 2.0 3.0 4.0
double!(a)
println(a) # modification is reflected to caller,
# because there was only ever one array!
[2.0,4.0,6.0,8.0]
Summary: Your hard work learning Python can transfer well to Julia.
For more on how both languages treat names and values, http://nedbatchelder.com/text/names1.html is a great reference.
In Python, most of us are heavy users of numpy, which provides a ndarray
class for homogenous arrays. On the other hand, we also have Python's built-in list
type, which are heterogeneous 1-d arrays. It can sometimes be awkward dealing with two types that have such overlapping functionality.
x = [1, 'two', 3.0] # heterogeneous
y = np.array([1.0, 2.0, 3.0]) # homogeneous
(I end up starting a lot of functions with x = np.asarray(x)
.)
# equivalent of Python list or ndarray with dtype='object'
a = [1.0, 2, "three", 4+0im]
4-element Array{Any,1}: 1.0 2 "three" 4+0im
typeof(a) # a is an array of heterogenous objects
Array{Any,1}
map(typeof, a)
4-element Array{Any,1}: Float64 Int64 ASCIIString Complex{Int64}
# equivalent of Python ndarray with dtype=float64
b = [1.0, 2.0, 3.0, 4.0]
typeof(b)
Array{Float64,1}
# array only takes up 4 * 8 bytes
sizeof(b)
32
You can't do this (efficiently) in NumPy.
immutable Point
x::Float64
y::Float64
end
x = [Point(1., 2.), Point(3., 4.), Point(5., 6.)]
3-element Array{Point,1}: Point(1.0,2.0) Point(3.0,4.0) Point(5.0,6.0)
sizeof(x) # points are stored efficiently in-line
48
For performance in Python, you'd have to do something like
class Points(object):
"""A container for two arrays giving x and y coordinates."""
def __init__(self, x, y):
self.x = x
self.y = y
def add_offset(self, x_offset, y_offset):
self.x += x_offset
self.y += y_offset
# ... other methods that operate element-wise
What you really want is a Point
object, but if you write classes that way in Python, performance will suffer.
# two 200 x 200 matricies
n = 200
A = rand(n, n)
B = rand(n, n);
f(A, B) = 2A + 3B + 4A.*A # function we want to optimize
f (generic function with 1 method)
using TimeIt
@timeit f(A, B);
1000 loops, best of 3: 315.26 µs per loop
We get similar performance in Python:
In [5]: n = 200
In [6]: from numpy.random import rand
In [7]: A = rand(n, n);
In [8]: B = rand(n, n);
In [9]: %timeit 2*A + 3*B + 4*A*A
1000 loops, best of 3: 354 µs per loop
But if needed to optimize this further, we'd have to reach for a specialized tool such as cython, numba, ...
function f2(A, B)
length(A) == length(B) || error("array length mismatch")
C = similar(A, promote_type(eltype(A),eltype(B)))
@inbounds for i=1:length(C)
C[i] = 2A[i] + 3B[i] + 4A[i]*A[i]
end
return C
end
f2 (generic function with 1 method)
@timeit f2(A, B);
10000 loops, best of 3: 65.61 µs per loop
function f3!(A, B, C)
length(A) == length(B) == length(C) || error("array length mismatch")
@inbounds @simd for i=1:length(C)
C[i] = 2A[i] + 3B[i] + 4A[i]*A[i]
end
end
f3! (generic function with 1 method)
C = similar(A, promote_type(eltype(A),eltype(B)))
@timeit f3!(A, B, C);
10000 loops, best of 3: 48.96 µs per loop
As an example, we have the sine function built-in:
sin(1.0)
0.8414709848078965
But what if we want to call a different implementation of the sine function?
ccall((:sin, "libopenlibm"), Float64, (Float64,), 1.0)
0.8414709848078965
In general, if you've compiled code to a shared library, libname.so
, you can call it with:
ccall((:function_name, "libname"),
ReturnType,
(Arg1Type, Arg2Type,...),
arg1, arg2, ...)
Packages and managing dependencies are super important. Julia's Pkg
is declarative (like conda).
Pkg.add("Cosmology")
would add "Cosmology" to the requirements file:
;cat ~/.julia/v0.4/REQUIRE
IJulia ForwardDiff Requests HTTPClient DocOpt Example Gadfly Winston Logging SIUnits PyCall GeneralizedSampling JLD
Julia figures out dependencies and installs the optimal version of every package to satisfy dependencies minimally.
In the browser: http://juliabox.com
At NERSC/Cori (experimental support):
module load julia
PyCall
is pretty good.)