Assistant Professor Justin Dressel
Faculty of Mathematics, Physics, and Computation
Schmid College of Science and Technology
Task 1:
Why Julia?
_
_ _ _(_)_ | A fresh approach to technical computing
(_) | (_) (_) | Documentation: https://docs.julialang.org
_ _ _| |_ __ _ | Type "?help" for help.
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 0.6.0 (2017-06-19 13:05 UTC)
_/ |\__'_|_|_|\__'_| | Official http://julialang.org/ release
|__/ | x86_64-pc-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.
julia> x = linspace(-2,2,10000)
julia> x .^ 2
julia> x .^ 2
julia> x .^ 3
julia> typeof(x)
StepRangeLen{Float64,...}
julia> typeof(x .^ 2)
Array{Float64,1}
(define x as a range of values)
(compiles x here to square each element,
Note: .^ applies ^ across the whole array)
(subsequent operations that use the compiled x are significantly faster)
(StepRangeLen 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)
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.)
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 REPL[9]:1
square(x) at REPL[6]:1
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
Official Documentation : https://docs.julialang.org/en/stable/
Textbook Introduction: https://lectures.quantecon.org/jl/learning_julia.html
Cheat Sheet : https://brakmic.github.io/Julia-Cheat-Sheet/
Cheat Sheet : https://github.com/QuantEcon/QuantEcon.cheatsheet/blob/master/julia-cheatsheet.rst
Task 2 :
Task 3 :
"""
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
Task 4 :
(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)
# Evaluate j0 on single complex points. Note: im is the imaginary unit in Julia
print( j0( complex(1.1, 0.3) ) ) # Recommended construction for complex numbers
print( j0( 1.1 + 0.3im ) ) # Equivalent result, but constructs z in 2 steps
# Evaluate j0 across an array - the . notation automatically vectorizes any function
a = linspace(complex(0.1,0.3), complex(1.5,0.3), 100)
print( j0.(a) )
Create the following cell and execute it. Discuss what it is doing in your notebook.
Note in particular the line showing that the convention f.(a) for a function f vectorizes the function f across an entire array of arguments
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 cplane = complex_plane()
@time j0p = j0.(cplane)
Discuss exactly how this code works. What's the difference between the comma and the semicolon in the list of arguments?
Create the following cell and run it:
mutable struct ComplexPlane
x :: LinSpace{Float64}
y :: LinSpace{Float64}
z :: Array{Number,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:
cplane = ComplexPlane(xpoints=200,ypoints=200);
typeof(cplane)
How does a mutable struct compare to a Python class?
How is the constructor like the __init__ method of a Python class?
print(typeof(cplane.x))
cplane.z = j0.(cplane.z)
# 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.build("PyPlot) # Check that everything is built properly
# This should only need to be done once to verify it is installed properly
For plotting, we can call Python's matplotlib directly from within Julia
c = -1.037 + 0.17im # Set starting point of julia set
j(z) = juliamap(c, z) # Create julia map
cplane = ComplexPlane() # Create 2000x2000 point complex plane
jp = j.(cplane.z); # Apply julia map to entire plane
using PyPlot # Load PyPlot package into the current namespace
figure(1)
xlabel("Re(z)")
ylabel("Im(z)")
title("Julia Set: c = " * string(c))
pcolormesh(cplane.x, cplane.y, jp, cmap=PyPlot.cm_get_cmap("hot"))
savefig("julia.png") # Also 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 ~5 years old. However, it is very rapidly becoming a mature and promising scientific programming language.
It will reach version 1.0 in 2018.
My prediction is that julia will become a major language for data science and scientific computing.
Evidence: https://lectures.quantecon.org/jl
Task 5 :