Assistant Professor Justin Dressel
Faculty of Mathematics, Physics, and Computation
Schmid College of Science and Technology
Task 1:
Why Julia?
~$ 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.
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)
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 none:1
square(x) at none: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
Quick Summary : http://learnxinyminutes.com/docs/julia/
Cheat Sheet : https://dl.dropboxusercontent.com/u/8252984/julia.html
Official Documentation : http://docs.julialang.org/en/release-0.4/manual/introduction/
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)
# 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?
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?
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)
# 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 :