Julia and IJulia Overview

Assistant Professor Justin Dressel

Faculty of Mathematics, Physics, and Computation

Schmid College of Science and Technology

 

Julia

Task 1:

  • Log into your Sage Cloud account, open Terminal.
  • Run the command "julia"
    You should be greeted with the following julia prompt:

Why Julia?

  • Easy to learn, but well-designed; new, and rapidly growing
  • Philosophy:  code should be as easy to write and read as Python and MATLAB, but execute (almost) as efficiently as FORTRAN and C
  • Rapid growth in the past few years, with great promise
  • As a just-in-time (JIT) compiled language using the LLVM, it is a nice middle ground between interpreted (Python) and compiled (C, Fortran)
~$ julia
               _
   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: http://docs.julialang.org
   _ _   _| |_  __ _   |  Type "?help" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.4.0 (2015-10-08 06:20 UTC)
 _/ |\__'_|_|_|\__'_|  |  Official http://julialang.org release
|__/                   |  x86_64-linux-gnu
 
julia>

The julia prompt acts like an interpreter, but runs every line of code through a compiler before executing that compiled code.  This makes execution significantly faster after the compilation step than an interpreted language like Python. 

JIT Compilation

julia> x = linspace(-2,2,10000)

julia> x .^ 2

julia> x .^ 2

julia> x .^ 3

julia> typeof(x)                                                                                                                                                                                                              
LinSpace{Float64}
 
julia> typeof(x .^ 2)                                                                                                                                                                                                         
Array{Float64,1}

(define x as a LinSpace)

(compiles x here to square each element)

(subsequent operations that use the compiled x are significantly faster)

(LinSpace in Julia is like a generator in Python)

(Array in Julia is like a NumPy Array in Python)

Try running these!

(Note: syntax deliberately similar to MATLAB)

Type Inference

Structurally, julia is primarily a functional language with a multiple dispatch inferential type system

julia> square(x) = x .^ 2                                                                                                                                                                                                     
square (generic function with 1 method)
 
julia> square(3)
 
julia> square(1:10)

julia> square(s :: AbstractString) = s * s                                                                                                                                                                                    
square (generic function with 2 methods)
 
julia> square("foo")

julia> square(3)




julia> square2(x :: Real) = x .^ 2                                                                                                                                                                                            
square2 (generic function with 1 method)
 
julia> square2(3)                                                                                                                                                                                                             
 
julia> square2(1:10)                                                                                                                                                                                                          
ERROR: MethodError: `square2` has no method matching square2(::UnitRange{Int64})

(The type of x is inferred to be the most general possible type that supports the element-wise power operation .^)

(Here the type of x is explicitly specified to a Real number only. More specificity can mean better compiler optimization, but restricts generality.  Use sparingly)

(Here the type of x is restricted to a string.  When called, square will pick the most restricted type-match to execute first.)

Multiple Dispatch

Notably, julia is not object-oriented like Python.

The multiple-dispatch type system is much more flexible.

abstract Number
abstract Real     <: Number
abstract AbstractFloat <: Real
abstract Integer  <: Real
abstract Signed   <: Integer
abstract Unsigned <: Integer

Abstract Types themselves have a hierarchy, like classes in Python.  Here <: means "is a subtype of"

Methods are not attached to a single type as with Python classes.  Instead, all functions become "methods" by specifying the argument types they act on.

julia> methods(square)                                                                                                                                                                                                        
# 2 methods for generic function "square":
square(s::AbstractString) at none:1
square(x) at none:1

Composite Types

julia> type Foo                                                                                                                                                                                                               
         x :: Int                                                                                                                                                                                                             
         y :: Float64                                                                                                                                                                                                         
       end                                                                                                                                                                                                                    
 
julia> f = Foo(3, 2.0)                                                                                                                                                                                                        
Foo(3,2.0)                                                                                                                                                                                                                    
                                                                                                                                                                                                                              
julia> f.x                                                                                                                                                                                                                    
3
 
julia> f.y                                                                                                                                                                                                                    
2.0
 
julia> g = Foo(3,4)                                                                                                                                                                                                           
Foo(3,4.0)
 
julia> g.x                                                                                                                                                                                                                    
3
 
julia> g.y                                                                                                                                                                                                                    
4.0

Defining structured data types (known as structs in C) is simple. Such composite types behave in the same way as Python object attributes

Specifying the type automatically coerces the type in constructor (if it can) to ensure type safety and optimized compilation

How to Learn a Language Quickly

  1. Identify basic syntax
  2. Identify primary data types
  3. Identify primary control structures
  4. Identify primary organizational structures
  5. Find and emulate idiomatic examples
  6. Don't be afraid to play until you understand

Task 2 :

  • Open up these links as tabs in your browser for reference

Task 3 :

  • Create a new Jupyter Notebook:  learnjulia.ipynb
  • In the Kernel menu, go to Change kernel => Julia
    ​​    (You are now using Julia inside a Jupyter notebook!)
  • Create a Markdown cell at the top, add a title, then a subtitle with your name(s).
  • By investigating the links on the previous slide, playing with the interpreter, and discussion, answer the following questions in your notebook:
    1. What is the julia syntax for the following?
          comments, variable declarations, printing output
    2. What are the basic julia data types?  How do you define and use them?
    3. What basic control structures are available in julia? (for, while, etc.)
    4. What are some notable differences between Python and julia?  What are some notable similarities?
  • For each of these questions, create an example cell of working julia code that illustrates the answer, as well as descriptive Markdown cells.  After you are done, there should be no extra fluff in the notebook, and it should run without errors.
  • Try writing a function that generates the first 100 Fibonacci numbers in julia.
    (The start of this sequence is 1, 1, 2, 3, 5, 8, 13, ..., obtained by successive pair-wise sums.)

Jupyter Notebooks

"""
juliamap(c,z; maxiter) :
  Implement the iteration algorithm for a Julia Set.

**Returns:** integer number of iterations, or zero
if the iteration never diverges.

 - c : complex constant definining the set
 - z : complex number being iterated
 - maxiter : maximum iteration number, defaults to 100
"""
function juliamap(c, z; maxiter=100)
    for n = 1:maxiter
        z = z^2 + c
        if abs(z) > 2
            return n
        end
    end
    return 0
end

@doc juliamap

Julia Sets in Julia

Task 4 :

  • Create a new julia notebook: juliainjulia.ipynb
  • Add the following cell

(Note the different docstring convention from Python, which renders as Markdown.)

(The @doc macro finds and prints the docstring.)

# Specialize juliamap to c=0
j0(z) = juliamap(0,z)

# Vectorize j0 over arrays of Complex numbers
@vectorize_1arg Complex j0

# List the available methods for j0 for different types
methods(j0)

Create the following cell and execute it.  Discuss what it is doing in your notebook.

::AbstractArray{T<:Complex{T<:Real},N}

Note this type:

An AbstractArray of dimension N containing any type T that is a subtype of Complex, which is composed of parts of any (different) type T that is a subtype of Real.

Why does j0 now include a method that handles this type?

Complex Plane:  Take 1

Create the following cell and run it:

# Create a complex plane
function complex_plane(xmin=-2, xmax=2, ymin=-2, ymax=2; xpoints=2000, ypoints=2000)
    # y is a column vector
    y = linspace(ymin, ymax, ypoints)

    # x uses a transpose, yielding a row vector
    x = linspace(xmin, xmax, xpoints)'

    # z uses broadcasted addition and multiplication to create a plane
    z = x .+ y.*im;

    # The final line of a block is treated as the return value, in the absence
    # of an explicit return statement
end

Create the following cell and run it:

# The vectorized function can be applied directly to the plane
@time cp = complex_plane()
@time j0p = j0(cp)

Discuss exactly how this code works.  What's the difference between the comma and the semicolon in the list of arguments?

Complex Plane:  Take 2

Create the following cell and run it:

immutable ComplexPlane
    x :: LinSpace{Float64}
    y :: LinSpace{Float64}
    z :: Array{Complex{Float64},2}
    
    function ComplexPlane(xmin=-2, xmax=2, ymin=-2, ymax=2;
                            xpoints=2000, ypoints=2000)
        x = linspace(xmin, xmax, xpoints)
        y = linspace(ymin, ymax, ypoints)
        z = x' .+ y.*im
        new(x,y,z)
    end
end

Create the following cells and run them:

cp = ComplexPlane(xpoints=200,ypoints=200);
typeof(cp)

What is being done here? What is "immutable" vs. "type"?

  How does this method compare to a Python class?

print(typeof(cp.x))
j0(cp.z)

Plotting with PyPlot

# Run the following in a terminal julia interpreter, not the notebook!
julia> Pkg.add("PyPlot")    #  This downloads the package via git and installs it
julia> Pkg.status()         #  Check that everything is installed properly

For plotting, we can call Python's matplotlib directly from within Julia

using PyPlot                         #  Load package into the current namespace
c = -1.037 + 0.17im                  #  Set starting point of julia set
j(z) = juliamap(c, z)                #  Create julia map
@vectorize_1arg Complex j            #  Vectorize julia map
cp = ComplexPlane()                  #  Create 2000x2000 point complex plane
jp = j(cp.z);                        #  Apply julia map to entire plane

xlabel("Re(z)")
ylabel("Im(z)")
title("Julia Set: c = " * string(c))
pcolormesh(cp.x, cp.y, jp, cmap=PyPlot.cm_get_cmap("hot"))
display(gcf())                       #  Get Current Figure and display in notebook

# alternatively, instead of displaying in notebook, replace previous line with
# savefig("julia.png")               #  Output figure to png file

Julia can conveniently install packages automatically via git and the web.

NOTE

Julia first appeared in 2012, so it is only ~3 years old.  However, it is very rapidly becoming a mature and promising scientific programming language.

 

It is worth keeping track of how the language keeps maturing in the near future.

Task 5 :

  • Change the seed value (c) in the code for generating Julia Sets to see what happens.  Find values of c that you like.
  • Plot 5 of your favorite Julia Sets in your notebook.
  • Commit the notebook in GitHub.