View All Pages

Julia instructions

Help for Julia users in SF2524 & SF3580

This page can be modified both by students and teachers.

You are welcome to improve it.

Something missing? Try contacting eringh at kth se

Introduction

The material treated in the tutorial isfound in the folder SF2524/Files/tutorial.

Using Julia at KTH

Julia should already be available on the (Ubuntu?) computers at KTH, at least the computers administered by the Math department. If you want to use another version you can simply download the precompiled binaries to a folder and run that Julia-executable. To add it to your path you can:

  • Either, add a folder called "bin" to your home folder and add a symbolic link (with a name of you chocie = the command you type to start Julia) in the folder pointing to your Julia executable. Then you modify your path by adding "export PATH=$PATH:afs/kth.se/home/<PATH TO YOUR HOME FOLDER>/bin" to the file "<HOME>/.dotfiles/.bashrc".
  • Or, create an "alias" (alias name = the command you type to start Julia) by adding "alias <ALIAS NAME>=<PATH TO YOUR JULIA EXECUTABLE>" to the file "<HOME>/.dotfiles/.bashrc".

For more on the Ubuntu environment at KTH, see this page.

Common mistakes / problems

If you initiate vectors use concrete types and not abstract types. Instead of creating vectors of the abstract type Complex, use the concrete ComplexF64. The following example illustrates the cost of implicit conversions.

using BenchmarkTools, Random
n=1000;
x=zeros(ComplexF64,n);
A=randn(n,n);
@btime A*x;
@show typeof(A*x);
x=zeros(Complex,n);
@btime A*x;
@show typeof(A*x);

Running the code produces the output below, which tells us that the operations with concrete typing is in this case about 100 times faster.

551.700 μs (1 allocation: 15.75 KiB)
typeof(A * x) = Array{Complex{Float64},1}  
39.440 ms (3005001 allocations: 76.42 MiB)
typeof(A * x) = Array{Complex,1}

The example also shows how to use BenchmarkTools Links to an external site. to test performance of code. However, note that it is in general better to encapsulate the code in a function when testing performance.

 

Regarding HW1:

Supporting files are vailable in files/skeleton_codes_hw1_julia.

Loading matlab files

MAT Links to an external site. is the Julia package to interface for  Matlab data files

julia> ]
(v1.0) pkg> add MAT
using MAT

Then, for loading the data Bwedge.mat of the homework, just run

B=matread("Bwedge.mat")["B"];

Loading matrices from the "gallery"

In Julia we have access to most of the test matrices of the matlab gallery Links to an external site.  by simply installing the package MatrixDepot Links to an external site.. The following example shows how to load the matrix of hw1 second exercise. For linux users: you may need to install Zlib 1.2.9.

julia> ]
(v1.0) pkg> add MatrixDepot

Then, the following code will load the matrix "Wathen", referenced in Homework 1.

using MatrixDepot, Random
nn=10;
Random.seed!(0);
A=matrixdepot("wathen",nn,nn);

Loading the wathen matrix

The MatrixDepot package above is quite large and has complicated dependencies which cannot always been satisfied. Problems reported mostly on windows.  Here is the code to create the wathen matrix. Use as wathen(10,10). This code is copied from https://github.com/JuliaLinearAlgebra/MatrixDepot.jl/blob/51634d1c507e5300aebb7647f2e482dc892c6c06/src/higham.jl#L1521 Links to an external site.

 

using LinearAlgebra, SparseArrays;


function wathen(::Type{T}, nx::Integer, ny::Integer) where T
    e1 = T[6 -6 2 -8;-6 32 -6 20;2 -6 6 -6;-8 20 -6 32]
    e2 = T[3 -8 2 -6;-8 16 -8 20;2 -8 3 -8;-6 20 -8 16]
    e3 = [e1 e2; e2' e1]/45
    n = 3 * nx * ny + 2 * nx + 2 * ny + 1
    ntriplets = nx * ny * 64
    Irow = zeros(Int, ntriplets)
    Jrow = zeros(Int, ntriplets)
    Xrow = zeros(T, ntriplets)
    ntriplets = 0
    rho = 100 * rand(nx, ny)
    node = zeros(T, 8)

    for j = 1:ny
        for i = 1:nx

            node[1] = 3 * j * nx + 2 * i + 2 * j + 1
            node[2] = node[1] - 1
            node[3] = node[2] - 1
            node[4] = (3 * j - 1) * nx + 2 * j + i - 1
            node[5] = (3 * j - 3) * nx + 2 * j + 2 * i - 3
            node[6] = node[5] + 1
            node[7] = node[5] + 2
            node[8] = node[4] + 1

            em = convert(T, rho[i,j]) * e3

            for krow = 1:8
                for kcol = 1:8
                    ntriplets +=  1
                    Irow[ntriplets] = node[krow]
                    Jrow[ntriplets] = node[kcol]
                    Xrow[ntriplets] = em[krow, kcol]
                end
            end

        end
    end
    return sparse(Irow, Jrow, Xrow, n, n)
end
wathen(nx,ny) = wathen(Float64,nx,ny)

 

 

Loading and manipulating videos in Julia

This show how you can load an image reshape it into a vector. Packages not yet installed can be installed in the same way as shown above. (For Linux users: you may need to install Zlib 1.2.9.)(For Windows users: Avoid non-english characters and spaces in file- and folder-names (issue Links to an external site.).)

using ImageView, Images, ImageIO, ImageMagick, Printf, FileIO, LinearAlgebra, Arpack
basefilename="./market_snapshots"
# Load the image into an image object, here the first frame
img=load(abspath(basefilename*"_0001.jpg")); 
# Visualize it
imshow(img)
# Determine the image size
sz=size(img); szv=sz[1]*sz[2]; 
# Reshape the image to a vector:
R=float(red.(img));
G=float(green.(img));
B=float(blue.(img));
v=vcat(vec(R), vec(G), vec(B))

One way to do video processing is to separate the frames of the video into separate image files. This has already been done for you in the homework so you can just download the zip-file with images.

(Not necessary for homework: If you want to do the conversion yourself on your own recorded video or laspalmas.mkv or plattan.mkv here, you can generate snapshot images in ubuntu with the commands: export frames_per_second=0.1; ffmpeg -i laspalmas.mkv -r $frames_per_second laspalmas_%04d.png)

The frames can be loaded into a matrix by appropriate for loop


m=5  # Number of frames to load
A=zeros(size(v,1),m) # Matrix to store all the frames in
for k=1:m
    fname=@sprintf("%s_%04d.jpg",basefilename,k);
    img=load(abspath(fname));
    println(fname)
    R=float(red.(img));
    G=float(green.(img));
    B=float(blue.(img));
    v=vcat(vec(R), vec(G), vec(B));
    A[:,k]=v;
end

The process to reshape an image into a vector can be reversed (and you need it for the homework).

v=A[:,3]; # Third column of A = third frame
vv=reshape(v,szv,3);
R=reshape(vv[:,1],sz[1],sz[2]);
G=reshape(vv[:,2],sz[1],sz[2]);
B=reshape(vv[:,3],sz[1],sz[2]);
newimg=RGB.(R,G,B);
imshow(newimg);

 

Regarding HW2:

Minimizing functions

In the hw2 you will need to minimize a function. In Julia you can use the package optim.jl Links to an external site., instead of Matlab's fminsearch Links to an external site.. Read a bit of the documentation in order to have a basic understanding of the functionalities. We provide an example below where we minimize the Rayleigh quotient of a symmetric matrix and check that the minimum value is an eigenvalue and the minimum point is an eigenvector (known from the theory).

using Optim
A=rand(5,5); A=A+A';
r(x) = (x'*A*x)/(x'*x);
x0 = rand(5);
result = optimize(r, x0, BFGS());
x=result.minimizer;
λ=result.minimum;
display(norm(A*x-λ*x)/norm(x))

 

Regarding HW3:

Supporting files are vailable in files/skeleton_codes_hw3_julia.