INTRODUCTION TO JULIA

In [ ]:

    

# Here is another feature - a built in package manager
installed_packages = Pkg.installed()
# len tab for autocomple
length(installed_packages)


    

In [ ]:

    

haskey(installed_packages, "IJulia")


    

In [ ]:

    

# Some unicode
a = 1:10;
11 ∈ a


    

In [ ]:

    

b = 4:20;
a ∩ b


    

In [ ]:

    

# Explore: Docstring contains the query
apropos("eigenvalues")


    

In [ ]:

    

"""Compute the Collatz chain for number n."""
function collatz_chain(n)
    k = 1
    while n > 1
        n = isodd(n) ? 3n+1 : n >> 1
        k += 1
        # println(n)
    end
    k
end

"""Which of the number [1, stop) has the longest Collatz chain."""
function solve_euler(stop)
    n, N, N_max = 1, 0, 0  
    while n < stop
        value = collatz_chain(n)
        if value > N_max
            N = n
            N_max = value
        end
        n += 1
    end
    (N, N_max)
end


    

In [ ]:

    

# Let's time it
N = 1000000
t0 = tic()
answer = solve_euler(N)
t1 = toc();
answer, t1


    

In [ ]:

    

# Julia code is a translation of the definition
"""Color Julia set. Code from https://www.youtube.com/watch?v=PsjANO10KgM"""
function julia(z, c)
    for n in 1:80
        if abs2(z) > 4
            return n-1
        end
        z = z*z + c
    end
    return 80
end

"""Color Julia set. Prealoc version"""
function julia(x, y, c)
    J = zeros(length(x), length(y))
    index = 0
    for r in y, i in x
        index += 1
        J[index] = julia(complex(r, i), c);
    end
    J
end


    

In [ ]:

    

# We explore a few points
cs = (complex(-0.06, 0.67), complex(0.279, 0), complex(-0.4, 0.6), complex(0.285, 0.01))

x = collect(1:-0.002:-1);
y = collect(-1.5:0.002:1.5);

Js = []
# Evaluate fractal generation
t0 = tic()
for c in cs
    #push!(Js, [julia(complex(r, i), c) for i in x, r in y]);
    push!(Js, julia(x, y, c));
end
t1 = toc()

@printf("Generated in %.4f s\n", t1);
println("Image size $(size(Js[1], 1))x$(size(Js[1], 2))")


    

In [ ]:

    

# Here we use python interopt(more on that later) to plot the fractals
using PyPlot

for J in Js
    figure(figsize=(12, 8))
    imshow(J, cmap="viridis", extent=[-1.5, 1.5, -1, 1])
end
show()


    

In [ ]:

    

# Read the content on the directory and keep only the files
files = filter(isfile, readdir("."))


    

In [ ]:

    

# How many of the files are Julia files
n_jl_files = count(f -> endswith(f, ".jl"), files)
println("$(pwd()) contains $(n_jl_files) files.")


    

In [ ]:

    

# Get those files and print them
jl_files = files[find(f -> endswith(f, ".jl"), files)]
for (i, file) in enumerate(jl_files)
    println("Julia file $(file)")
end


    

In [ ]:

    

# We have seen Array container. Let's show Sets
jl_files2 = append!(jl_files, jl_files)
jl_set = Set(jl_files2)
println("$(pwd()) contains $(length(jl_set)) files.")


    

In [ ]:

    

# Reading from file. How many words in the MIT licence?
LICENCE = open("LICENSE")
content = lowercase(readall(LICENCE))
words = split(content, r"\W", keep=false)
close(LICENCE)
length(words)


    

In [ ]:

    

# We could have done this with bash. Show shell interopt
cmd = pipeline(`cat LICENSE`, `wc -w`)
run(cmd)  # 172 vs 171, less than 1% error is okay:)


    

In [ ]:

    

# Count them with multiplicity
unique = Dict()
for word in words
    unique[word] = get(unique, word, 0) + 1
end


    

In [ ]:

    

# Most common ones
for (k, v) in sort(collect(unique), by=last, rev=true)
    println("$(k) $(v)")
end


    

In [ ]:

    

# : is => in the dictionary construction
dict = Dict(1 => "Miro", 2 => "Kuchta")


    

In [ ]:

    

# in checks for (key, value) pair
haskey(dict, 1), (1 => "Miro") ∈ dict


    

In [ ]:

    

# iteration is over (key, values) pairs
for (key, value) in dict
    println("$key => $value")  # "%.4g" % (...)
end


    

In [ ]:

    

# String concatenation is done with *. With + the operation should be commutative, which is not the case here
"Julia " * " Language"


    

In [ ]:

    

# / is for floating point division. Integer division is done by div
/(7, 4), div(7, 4)


    

In [ ]:

    

# Indexing starts at one. The last element in range is included. 0..last not included has same length as 1..last included
a = [1, 2, 3]
try
    a[0]
catch BoundsError
    println("No no, a[1]=", a[1], " ", a[end] == a[3] == 3)
end

collect(1:10)[10] == 10


    

In [ ]:

    

# complex numbers
complex(1, 1) == 1 + 1im


    

In [ ]:

    

# Convenience
x = 4
y = 2x  # no *
y == 8


    

In [ ]:

    

# More convenience. Save one "for"
[(i, j) for i in 1:4, j in "abcd"]


    

In [ ]:

    

# The arrays are column majored. In particular vector is a column vector
a = [1, 2, 3]


    

In [ ]:

    

# Contrast this with Python
A = reshape(collect(1:25), (5, 5))


    

In [ ]:

    

# In numpy you get to iterate over rows of matrix (for vectorization). In Julia this is not default. 
entries = [e for e in A]  #join(map(string, collect(A)), ", ")
# Moverover, columns are contiguous so accessing rows is not efficient.
# Finaly you iterate over values not pointers
for a in A
    a += 1
end
A


    

In [ ]:

    

# Column majored linear indexing
A[8] == 8


    

In [ ]:

    

# Multiline lambdas
map(f -> begin
    name, ext = splitext(f)
    uppercase(name)*uppercase(ext)^2
    end, readdir("."))


    

In [ ]:

    

# One lined named function definition
f(x, y) = 2x + 4y
f(1, 0) == 2


    

In [ ]:

    

# Implicit return statements
function f(name)
    if isequal(name, "Miro")
        return "slav"
    end
    "!"
end

(f("Miro"), f("You")) == ("slav", "!")


    

In [ ]:

    

# Variable positional arguments
mysum(args...) = join(map(string, args), "+")
mysum(10, 20)


    

In [ ]:

    

# Keyword arguments. Not a dictionary
function foo(args...; kwargs...)
    kwargs = Dict([k => v for (k, v) in kwargs])
    ans = join(map(string, args), "?")
    if haskey(kwargs, :x)  # Ignore :x for a while. It comes from the 'strange' kwargs array. Symbol
        ans = ans * "?" * string(kwargs[:x])
    end
    ans 
end
foo(10, 20, x=40)


    

In [ ]:

    

using PyCall


    

In [ ]:

    

# Let's load numpy and have it compute for us the eigenvalues of some matrix
@pyimport numpy.linalg as npla
A = diagm([1, 2, 3, 4])          # Julia matrix
eigv, eigw = npla.eig(A)         # Passed to python no copying
eigv                             # The result is Julia Array


    

In [ ]:

    

# Let's have numpy make the matrix and then Julia takes over for the actual work
@pyimport numpy as np
Apy = np.diag([1, 2, 3, 4])  # Note that this is Julia array passed to Numpy which comes back as a Julia matrix
eigvals(Apy)


    

In [ ]:

    

# Julia has native packages for plotting, but why give up matplotlib
x = linspace(-1, 1, 100)
y = sin(2π*x)
figure()
plot(x, y)
xlabel("\$x\$")  # Note that dollars have to be escaped for they have special role in Julia's string (interpolation)
show()


    

In [ ]:

    

#Your code is parsed, type annotated(types can be inferred). 
simplef(x, y) = x+y
@code_warntype simplef(1, 2)   # Note x, y have types here


    

In [ ]:

    

# When f is called with specific values JIT compiles the function. With type annotation the compiler can reason about
# the code and generate efficient code.
@code_llvm simplef(1, 2)        # Note @julia_simplef_FOO(i64, i64)


    

In [ ]:

    

# If called with values of same type. We get a cached function
@code_llvm simplef(1, 2)        # Note @julia_simplef_FOO(i64, i64)


    

In [ ]:

    

# Different value types will invoke JIT again
@code_llvm simplef(1, 2.)        # @julia_simplef_BAR(i64, double)


    

In [ ]:

    

map(typeof, (1, 2., [1, 2.]))


    

In [ ]:

    

# This reads as is lhs an instance of rhs
1::Number, 1::Real, 1::Int


    

In [ ]:

    

# But 
try
    1::Int8
catch TypeError
    println("Type of 1 is $(typeof(1)) and that is not a subtype of Int8")
end


    

In [ ]:

    

depth = 0
function printtypes(T)
    global depth
    if isleaftype(T)
        println("-"^(2*depth), "$T")
    else
        println("-"^depth, "$T")
        depth += 1
        map(printtypes, subtypes(T))
        depth -= 1
    end
end

printtypes(Real);


    

In [ ]:

    

# + is a function +(x, y) = ... . Which one is called based on args
@which +(1, 3)


    

In [ ]:

    

@which 0x1+0x3


    

In [ ]:

    

@which 1+2.


    

In [ ]:

    

"""Polynomial with value c[i]*sin(i*pi*x)"""
immutable SinePolynomial{T<:Real}
    coefs::Vector{T}
end


    

In [ ]:

    

# If not specified otherwise our parent is Any - the very tip of the type tree
super(SinePolynomial)


    

In [ ]:

    

# Our type is composite i.e. has fields
fieldnames(SinePolynomial)


    

In [ ]:

    

import Base: +, -, call
# Add () for evaluating the polynomial
call{T}(p::SinePolynomial{T}, x::Real) = sum([c*sum(k*π*x) for (k, c) in enumerate(p.coefs)])

# We can add/subtract polynomials. This is simple because of othogonality
function +{T, S}(p::SinePolynomial{T}, q::SinePolynomial{S})
    np, nq = length(p.coefs), length(q.coefs)
    n = max(np, nq)
    p = [p.coefs; zeros(T, n-np)]
    q = [q.coefs; zeros(S, n-nq)]
    SinePolynomial(p+q)
end

function -{T, S}(p::SinePolynomial{T}, q::SinePolynomial{S})
    np, nq = length(p.coefs), length(q.coefs)
    n = max(np, nq)
    p = [p.coefs; zeros(T, n-np)]
    q = [q.coefs; zeros(S, n-nq)]
    SinePolynomial(p-q)
end

# See that + has been added to lookup table for dispatching
methods(+)


    

In [ ]:

    

import Base.*
# We better prevent multiplications
*{T, S}(p::SinePolynomial{T}, q::SinePolynomial{S}) = throw(ArgumentError("No multiplication until CosinePolynomial"))

p = SinePolynomial([1, 2, 3])
q = SinePolynomial([2, -1])

try
    p*q
catch ArgumentError
    println("Is it okay? ", applicable(*, (typeof(p), typeof(q))))
    nothing
end


    

In [ ]:

    

Base.call{T}(p::SinePolynomial{T}, X) = [p(x) for x in X]

x = linspace(-1, 1, 100)
length(p(x))


    

In [ ]:

    

@which p(x)         # See which one was celled


    

In [ ]:

    

# Type INstability of s
function unstable(n)
    s = 0          # S is Float
    for i in 1:n
        s += s/i   # S is Float
    end
end

# Type stability of s
function stable(n)
    s = 0.0         # S is Float
    for i in 1:n
        s += s/i   # S is Float
    end
end


    

In [ ]:

    

@time stable(1000000)


    

In [ ]:

    

@time unstable(1000000)


    

In [ ]:

    

println("The speed difference is huge ", 0.063284/0.000002, ". And so is the difference in the number of lines of generated code")


    

In [ ]:

    

# Run @code_warntype here as well to see the (un)certainty
@code_llvm stable(3)


    

In [ ]:

    

@code_llvm unstable(3)


    

In [ ]:

    

# Example setenv via libc
function my_setenv(var::AbstractString, value::AbstractString)
    ccall((:setenv, "libc"), Cint, (Ptr{UInt8}, Ptr{UInt8}, Cint), var, value, 1)
end


    

In [ ]:

    

my_setenv("TODAY", "TUESDAY");


    

In [ ]:

    

# We will write getenv to check it
function my_getenv(var::AbstractString)
    bytestring(ccall((:getenv, "libc"), Ptr{UInt8}, (Ptr{UInt8}, ), var))
end


    

In [ ]:

    

my_getenv("TODAY")


    

In [ ]:

    

# Polynomial evaluation using GNU scientific library. Parametrized functions, restricted types
"""Evaluate polynomial c[1]*x^0 + c[2]*x^1 ... at x"""
function polyval{T<:Real, S<:Real}(c::Vector{T}, x::S)
    c = convert(Vector{Float64}, c)
    n = length(c)
    value = ccall((:gsl_poly_eval, "libgsl"), Cdouble, (Ptr{Cdouble}, Cint, Cdouble), c, n, x)
    Float64(value)
end

"""Vectorized version"""
function polyval{T, S}(c::Vector{T}, X::Vector{S})
    [polyval(c, x) for x in X]
end


    

In [ ]:

    

# Plot some random polynomials
x = collect(linspace(-1, 1, 100)) # 1im will

figure()
for i in 1:4
    c = rand(5)
    y = polyval(c, x)
    plot(x, y)
end
show()


    

In [ ]:

    

# A basis building block is symbol
typeof(:x)


    

In [ ]:

    

# Symbols can be combined into expressions
f_def = quote
    cheer(name) = name * "!"
end

# Meta.show_sexpr(f_def)
# f_def
dump(f_def)


    

In [ ]:

    

# And expression can be manipulated. Change the name to yell and !!
f_def.args[2].args[1].args[1] = :yell;
f_def.args[2].args[2].args[2].args[3] = "!!!";
f_def


    

In [ ]:

    

# Expressions can be evaluated. There will be no cheer.
eval(f_def);
isdefined(:cheer)


    

In [ ]:

    

# However we can yell
yell("Fire")


    

In [ ]:

    

# When is macro run?
macro makeuppercase(f)   
    info("I was called.")
    return Expr(:call, :uppercase, f)
end

# When the code is parsed the macro executes and changes body of f 
function transform(text::AbstractString)
    @makeuppercase text
end


    

In [ ]:

    

# It does not exacute when the function is run.
transform("ApuhuAs")


    

In [ ]:

    

# We try to write the same code but evaluation of the polynomial will be done via call to external routine.
function _horner(x, p...)
    c = p[end]
    for i in length(p)-1:-1:1
        c = p[i] + c*x
    end
    c
end

function _erfinv(x::Float64)
    a = abs(x)
    if a >= 1.0
        if x == 1.0
            return Inf
        elseif x == -1.0
            return -Inf
        end
        throw(DomainError())
    elseif a <= 0.75 # Table 17 in Blair et al.
        t = x*x - 0.5625
        return x * _horner(t, 0.16030_49558_44066_229311e2,
                             -0.90784_95926_29603_26650e2,
                              0.18644_91486_16209_87391e3,
                             -0.16900_14273_46423_82420e3,
                              0.65454_66284_79448_7048e2,
                             -0.86421_30115_87247_794e1,
                              0.17605_87821_39059_0) /
                   _horner(t, 0.14780_64707_15138_316110e2,
                             -0.91374_16702_42603_13936e2,
                              0.21015_79048_62053_17714e3,
                             -0.22210_25412_18551_32366e3,
                              0.10760_45391_60551_23830e3,
                             -0.20601_07303_28265_443e2,
                              0.1e1)
    elseif a <= 0.9375 # Table 37 in Blair et al.
        t = x*x - 0.87890625
        return x * _horner(t, -0.15238_92634_40726_128e-1,
                               0.34445_56924_13612_5216,
                              -0.29344_39867_25424_78687e1,
                               0.11763_50570_52178_27302e2,
                              -0.22655_29282_31011_04193e2,
                               0.19121_33439_65803_30163e2,
                              -0.54789_27619_59831_8769e1,
                               0.23751_66890_24448) /
                   _horner(t, -0.10846_51696_02059_954e-1,
                               0.26106_28885_84307_8511,
                              -0.24068_31810_43937_57995e1,
                               0.10695_12997_33870_14469e2,
                              -0.23716_71552_15965_81025e2,
                               0.24640_15894_39172_84883e2,
                              -0.10014_37634_97830_70835e2,
                               0.1e1)
    else # Table 57 in Blair et al.
        t = 1.0 / sqrt(-log(1.0 - a))
        return _horner(t, 0.10501_31152_37334_38116e-3,
                          0.10532_61131_42333_38164_25e-1,
                          0.26987_80273_62432_83544_516,
                          0.23268_69578_89196_90806_414e1,
                          0.71678_54794_91079_96810_001e1,
                          0.85475_61182_21678_27825_185e1,
                          0.68738_08807_35438_39802_913e1,
                          0.36270_02483_09587_08930_02e1,
                          0.88606_27392_96515_46814_9) /
              (copysign(t, x) *
               _horner(t, 0.10501_26668_70303_37690e-3,
                          0.10532_86230_09333_27531_11e-1,
                          0.27019_86237_37515_54845_553,
                          0.23501_43639_79702_53259_123e1,
                          0.76078_02878_58012_77064_351e1,
                          0.11181_58610_40569_07827_3451e2,
                          0.11948_78791_84353_96667_8438e2,
                          0.81922_40974_72699_07893_913e1,
                          0.40993_87907_63680_15361_45e1,
                          0.1e1))
    end
end


    

In [ ]:

    

# Make sure it works
erfinv(0.002) - _erfinv(0.002)


    

In [ ]:

    

# Is there a difference in performance?
@time [erfinv(x) for x in -1:0.00001:1];


    

In [ ]:

    

# _erfinv is about 2x slower. Note the difference in memory
@time [_erfinv(x) for x in -1:0.00001:1];


    

In [ ]:

    

# We make the polynomial degree available to the type system
"x^n"
type Monomial{n}
end


    

In [ ]:

    

# We generate code based on n
@generated function evaluate{n}(::Monomial{n}, x::Real)
    info("I was called")
    coefs = zeros(Real, n+1)
    coefs[n+1] = 1
    quote
        return Base.Math.@horner x $(coefs...)
    end
end


    

In [ ]:

    

p = Monomial{2}()
evaluate(p, 1.0)  # Print message


    

In [ ]:

    

evaluate(p, 2.0)  # No message


    

In [ ]:

    

evaluate(Monomial{4}(), 3.0)  # New type -> message


    

In [ ]:

    

using PyCall
@pyimport scipy.fftpack as pyfft
# Show that FFT works
v = rand(20)

fft_jl = FFTW.r2r(v, [FFTW.REDFT10])
fft_py = pyfft.dct(v)

norm(fft_jl - fft_py)


    

In [ ]:

    

# Native plotting with Gadfly
using Gadfly
plot(y=fft_jl)


    

In [ ]:

    

# Linear algebra. The favourite Poisson problem on UnitInterval. Convergence of CG1 FEM
A = 0  # Keep A outside for loop
for ncells in [2^i for i in 8:14]
    h = 1/ncells
    x = Float64[h*i for i in 0:ncells]
    F = sin(2π*x)
    U = sin(2π*x)/4/π^2
    
    # Based on type of A specialized solver will be called
    A = SymTridiagonal(fill(2, ncells+1), fill(-1, ncells))/h
    b = h*F

    u = A\b
    e = norm(u - U)
    println(e)
end

# Which factorize was dispatched
@which factorize(A)


    

In [ ]:

    

# And now the eigenvalues. Convergence of the smallest eig of -u'' = lambda* u in (0, 1) with homog. Dirichlet bcs
for ncells in [2^i for i in 3:11]
    h = 1/ncells
    x = Float64[h*i for i in 0:ncells]

    A = SymTridiagonal(fill(2, ncells+1), fill(-1, ncells))/h
    M = SymTridiagonal([1.; fill(4h/6, ncells-1); 1.], fill(h/6, ncells))
    
    eigw, eigv = eig(full(A), full(M))                          # Dense
    # eigw = first(eigs(sparse(A), sparse(M), nev=6, which=:SM))  # Sparse

    lmin = minimum(eigw)
    println("System size $(length(eigw)), error in λ=$(abs(lmin-π^2))")
end