Julia Community ๐ŸŸฃ

Vinod V
Vinod V

Posted on

Julia Vectors and Matrices

What are Julia vectors?

Julia vectors

Julia provides:

- Builtin multi-dimensional arrays

- closer to hardware (efficiency)

- designed for scientific computation (convenience)
Enter fullscreen mode Exit fullscreen mode
    julia> v = [0, 1, 2, 3]
      4-element Vector{Int64}:     
       0                           
       1                           
       2                           
       3                           
Enter fullscreen mode Exit fullscreen mode

For example, a vector containing:

  • values of an experiment/simulation at discrete time steps

  • signal recorded by a measurement device, e.g. sound wave

  • pixels of an image, grey-level or colour

  • 3-D data measured at different X-Y-Z positions, e.g. MRI scan

  • ...

Why it is useful: Memory-efficient container that provides fast numerical
operations.

    julia> using BenchmarkTools #used to evaluate the time taken for code execution
    julia> L = 1:1000
    julia> @btime [a^2 for a in L];
  1.342 ยตs (2 allocations: 7.97 KiB)
Enter fullscreen mode Exit fullscreen mode

Julia supports multidimensional arrays

Julia Reference documentation

     julia> ? Matrix
    search: Matrix BitMatrix DenseMatrix StridedMatrix AbstractMatrix

      Matrix{T} <: AbstractMatrix{T}

      Two-dimensional dense array with elements of type T, often used to represent a mathematical matrix. Alias for Array{T,2}.

      See also fill, zeros, undef and similar for creating matrices.
Enter fullscreen mode Exit fullscreen mode
  • Looking for something:
     julia> mat<Tab><Tab> # press tab key two times
     MathConstants  Matrix
Enter fullscreen mode Exit fullscreen mode

Creating vectors

Manual construction of vectors

1-D:

    julia> a = [0, 1, 2, 3]

    julia> ndims(a)
    1
    julia> size(a)
    (4,)
    julia> length(a)
    4
Enter fullscreen mode Exit fullscreen mode

2-D, 3-D, ...

    julia> b = [0  1  2; 3 4 5]
    2ร—3 Matrix{Int64}:
     0  1  2
     3  4  5
    julia> ndims(b)
    2
    julia> size(b)
    (2, 3)
    julia> length(b)     # returns the total size = 2 x 3
    6

    julia> c = [1 3 5
               2 4 6;;;
               7 9 11
               8 10 12]
        2ร—3ร—2 Array{Int64, 3}:
        [:, :, 1] =
         1  3  5
         2  4  6

        [:, :, 2] =
         7   9  11
         8  10  12
Enter fullscreen mode Exit fullscreen mode

Exercise: Simple arrays

* Create a simple two dimensional array. First, redo the examples
  from above. And then create your own: how about odd numbers
  counting backwards on the first row, and even numbers on the second?
* Use the functions `length`, `size` on these arrays.
  How do they relate to each other? And to the ``ndims`` function on 
  the arrays?
Enter fullscreen mode Exit fullscreen mode

Functions for creating vectors

In practice, we rarely enter items one by one...
Enter fullscreen mode Exit fullscreen mode

Evenly spaced:

    julia> a = collect(1:10)
    [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]    
    julia> b = 1:2:10 # start, step ,end (inclusive), 
    julia> b
    [1, 3, 5, 7]
Enter fullscreen mode Exit fullscreen mode
  • or by number of points:
    julia> c =  range(0,1,length = 6)   # start, end, num-points
    0.0:0.2:1.0
    julia> d = LinRange(0,1,6) #less overhead, prone to floating point errors
    6-element LinRange{Float64, Int64}:
    0.0, 0.2, 0.4, 0.6, 0.8, 1.0
Enter fullscreen mode Exit fullscreen mode
  • Common arrays:
    julia> a = ones(3,3)
    3ร—3 Matrix{Float64}:
     1.0  1.0  1.0
     1.0  1.0  1.0
     1.0  1.0  1.0

    julia> b = zeros(3,3)
    3ร—3 Matrix{Float64}:
     0.0  0.0  0.0
     0.0  0.0  0.0
     0.0  0.0  0.0  

    julia> using LinearAlgebra

    julia> c = Matrix(1.0I,3,3) #identity matrix
    3ร—3 Matrix{Float64}:
     1.0  0.0  0.0
     0.0  1.0  0.0
     0.0  0.0  1.0
Enter fullscreen mode Exit fullscreen mode
  • random numbers:
  julia> using Random           

  julia> Random.seed!(123)      
  TaskLocalRNG()                

  julia> a = rand(5,) #uniform random numbers in the range [0,1]          
          5-element Vector{Float64}:    
           0.521213795535383            
           0.5868067574533484           
           0.8908786980927811           
           0.19090669902576285          
           0.5256623915420473    
   julia> b = randn(5,) #Gaussian random vector
        5-element Vector{Float64}:
          0.9809798121241488
          0.0799568295050599
          1.5491245530427917
         -1.3416092408832219
          0.41216163468296796   
Enter fullscreen mode Exit fullscreen mode

Exercise: Creating arrays using functions

  • Experiment with :, range, ones, zeros, I and diagm.
  • Create different kinds of arrays with random numbers.
  • Try setting the seed before creating an array with random values.
  • Look at the function empty. What does it do? When might this be useful?

construct 1 2 3 4 5
construct -5, -4, -3, -2, -1
construct 2 4 6 8
look what is in an empty() array
construct 15 equispaced numbers in range [0, 10]

Basic data types

You may have noticed that, in some instances, array elements are displayed with
a trailing dot (e.g. 2. vs 2). This is due to a difference in the
data-type used:

    julia> a = [1, 2, 3]
    julia> eltype(a)
    Int64
    julia> b = [1., 2., 3.]
    julia> eltype(b)
    Float64
Enter fullscreen mode Exit fullscreen mode
Different data-types allow us to store data more compactly in memory,
but most of the time we simply work with floating point numbers.
Note that, in the example above, Julia auto-detects the data-type
from the input.
Enter fullscreen mode Exit fullscreen mode

The default data type is floating point:

    julia> a = ones(3, 3)
    julia> eltype(a)
    Float64
Enter fullscreen mode Exit fullscreen mode

There are also other types:

Complex:

    julia> d = [1+2im, 3+4im, 5+6*1im]
    julia> eltype(d)
    Complex{Int64}
Enter fullscreen mode Exit fullscreen mode

Bool:

    julia> e = [true, false, false, true]
    julia> eltype(e)
    Bool
Enter fullscreen mode Exit fullscreen mode

Strings:

    julia> f = ["Bonjour", "Hello", "Hallo"]
    julia> eltype(f)
    String
Enter fullscreen mode Exit fullscreen mode

Much more:

Int8
UInt8
Int16
UInt16
Int32
UInt32
Int64
UInt64
Int128
UInt128
Float16
Float32
Float64
ComplexF32
ComplexF64
Bool
Enter fullscreen mode Exit fullscreen mode

Basic visualization

Now that we have our first data arrays, we are going to visualize them.

    julia> using Plots
    julia> x = 1:100;y = 1:100; 
    julia> plot(x, y)       # line plot    # doctest: +SKIP
Enter fullscreen mode Exit fullscreen mode
    julia> using PyPlot
    julia> PyPlot.plot(x, y)       # line plot    # doctest: +SKIP
Enter fullscreen mode Exit fullscreen mode
  • 1D plotting:
  julia> x = LinRange(0, 3, 20)
  julia> y = LinRange(0, 9, 20)
  julia> PyPlot.plot(x, y, "g^")  # ^ plot
Enter fullscreen mode Exit fullscreen mode

(img)[]

  • 2D arrays (such as images):

      julia> using PyPlot
      julia> image = rand(0:255,10,10)
      julia> PyPlot.imshow(image)
Enter fullscreen mode Exit fullscreen mode

(img)[]

Exercise: Simple visualizations

  • Plot some simple arrays: a cosine as a function of time and a 2D matrix.

Indexing and slicing

The items of an array can be accessed and assigned to the same way as
other Python sequences (e.g. lists):


    julia> a = collect(1:9)
    julia> a
    [1, 2, 3, 4, 5, 6, 7, 8, 9]
    julia> a[1], a[3], a[end]
    (1, 3, 9)
Enter fullscreen mode Exit fullscreen mode

warning

Indices begin at 1, like other Julia sequences (Fortran or Matlab).
In contrast, in C/C++/Python, indices begin at 0.

The usual Julia method for reversing a sequence is supported:

   julia> reverse(a)
   [9, 8, 7, 6, 5, 4, 3, 2, 1])
Enter fullscreen mode Exit fullscreen mode

For multidimensional arrays, indices are tuples of integers:

    julia> using LinearAlgebra

    julia> A = diagm([1,1,1,1])
    4ร—4 Matrix{Int64}:
     1  0  0  0
     0  1  0  0
     0  0  1  0
     0  0  0  1 

    julia> A[1, 1]
    1
    julia> A[2, 1] = 10 # second line, first column
    julia> A
    4ร—4 Matrix{Int64}:
      1  0  0  0
     10  1  0  0
      0  0  1  0
      0  0  0  1

    julia> A[1,:]
    4-element Vector{Int64}:
     1
     0
     0
     0

    julia> A[1:1,:]
    1ร—4 Matrix{Int64}:
     1  0  0  0
     julia> A[1]
     1
Enter fullscreen mode Exit fullscreen mode

Note:

  • In 2D, the first dimension corresponds to rows, the second to columns.
  • for multidimensional A, A[p] is interpreted by the column major concept. If p = k*m + q, A[p] = A[q+1,k+1].

Slicing: Arrays, like other Julia sequences can also be sliced:

    julia> a = collect(0:9)
    julia> a
    [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
    julia> a[2:3:9] # [start:step:end]
    [1, 4, 7]
Enter fullscreen mode Exit fullscreen mode

Note that the last index is included! :

    julia> a[1:4]
    array([0, 1, 2, 3])
Enter fullscreen mode Exit fullscreen mode

    julia> a[1:3]
    [0, 1, 2]
    julia> a[1:2:end]
    [0, 2, 4, 6, 8]
    julia> a[3:end]
    [2, 3, 4, 5, 6, 7, 8, 9]
Enter fullscreen mode Exit fullscreen mode

You can also combine assignment and slicing:

   julia> a = collect(0:9)                                               
           10-element Vector{Int64}:                                              
            0                                                                     
            1                                                                     
            2                                                                     
            3                                                                     
            4                                                                     
            5                                                                     
            6                                                                     
            7                                                                     
            8                                                                     
            9                                                              

   julia> a[5:end] .= 10                                                  
           6-element view(::Vector{Int64}, 5:10) with eltype Int64:               
            10                                                                    
            10                                                                    
            10                                                                    
            10                                                                    
            10                                                                    
            10                                                                   

   julia> b = collect(0:5)                                                
           6-element Vector{Int64}:                                               
            0                                                                     
            1                                                                     
            2                                                                     
            3                                                                     
            4                                                                     
            5                                                                     

   julia> a[5:end] .= reverse(b)                                         
           6-element view(::Vector{Int64}, 5:10) with eltype Int64:              
            5                                                                    
            4                                                                    
            3                                                                    
            2                                                                    
            1                                                                    
            0                                                                    

   julia> a                                                              
           10-element Vector{Int64}:                                             
            0                                                                    
            1                                                                    
            2                                                                    
            3                                                                    
            5                                                                    
            4                                                                    
            3                                                                    
            2                                                                    
            1                                                                    
            0                                                                    
Enter fullscreen mode Exit fullscreen mode

Exercise: Indexing and slicing

  • Try the different flavours of slicing, using start, end and step: starting from a linspace, try to obtain odd numbers counting backwards, and even numbers counting forwards.
  • Reproduce the slices in the diagram above. You may
    use the following expression to create the array:

    
        julia> A = repeat(collect(0:5),1,6)'
        6ร—6 adjoint(::Matrix{Int64}) with eltype Int64:
         0  1  2  3  4  5
         0  1  2  3  4  5
         0  1  2  3  4  5
         0  1  2  3  4  5
         0  1  2  3  4  5
         0  1  2  3  4  5
    
        julia> B = [0:10:50;]
        6-element Vector{Int64}:
          0
         10
         20
         30
         40
         50
    
        julia> A .+ B
        6ร—6 Matrix{Int64}:
          0   1   2   3   4   5
         10  11  12  13  14  15
         20  21  22  23  24  25
         30  31  32  33  34  35
         40  41  42  43  44  45
         50  51  52  53  54  55
    

Exercise: Array creation
Create the following arrays (with correct data types)

```julia [1 1 1 1;
1 1 1 1;
1 1 1 2;
1 6 1 1]

    [0. 0. 0. 0. 0.;
    2. 0. 0. 0. 0.;
    0. 3. 0. 0. 0.;
    0. 0. 4. 0. 0.;
    0. 0. 0. 5. 0.;
    0. 0. 0. 0. 6.]
Enter fullscreen mode Exit fullscreen mode


    *Hint*: Individual array elements can be accessed similarly to a list,
    e.g. ``a[1]`` or ``a[1, 2]``.

    *Hint*: Examine the help for ``diagm`` in LinearAlgebra module.

Exercise: Tiling for array creation

    Skim through the documentation for ``repeat``, and use this function
    to construct the array


```julia
[4 3 4 3 4 3;
2 1 2 1 2 1;
4 3 4 3 4 3;
2 1 2 1 2 1]
Enter fullscreen mode Exit fullscreen mode

Copies and views

A slicing operation creates a view on the original array, which is
just a way of accessing array data.

When modifying the view, the original array is modified as well:

    julia> using Statistics
    julia> x = rand(100)
    julia> @allocated cor(x[1:10],x[1:10])
    9204830
    julia> @allocated cor(view(x,1:10),view(x,11:20)) #or @allocated @views cor(x[1:10],x[1:10])
    8428206 #memory allocations are less 
Enter fullscreen mode Exit fullscreen mode

Exercise: [1, 2, 3, 4, 5] -> [1, 2, 3]
Exercise: [1, 2, 3, 4, 5] -> [4, 5]
Exercise: [1, 2, 3, 4, 5] -> [1, 3, 5]
Exercise: [1, 2, 3, 4, 5] -> [2, 4]
Exercise: create an array [1, 1, 1, 1, 0, 0, 0]
Exercise: create an array [0, 0, 0, 0, 1, 1, 1]
Exercise: create an array [0, 1, 0, 1, 0, 1, 0]
Exercise: create an array [1, 0, 1, 0, 1, 0, 1]
Exercise: create an array [1, 0, 2, 0, 3, 0, 4]

Worked example: Prime number sieve

Compute prime numbers in 0--99, with a sieve

  • Construct a shape (100,) boolean array is_prime, filled with True in the beginning:
    julia> is_prime = fill(true,100)
Enter fullscreen mode Exit fullscreen mode
  • Cross out 0 and 1 which are not primes:
   julia> is_prime[1] = false
Enter fullscreen mode Exit fullscreen mode
  • For each integer j starting from 2, cross out its higher multiples:
   julia> N_max = round(Int,sqrt(length(is_prime)))
   julia> for j in 2:N_max
                is_prime[2*j:j:end] .= false
        end
Enter fullscreen mode Exit fullscreen mode
  1. Skip ``j`` which are already known to not be primes

  2. The first number to cross out is :math:`j^2`
Enter fullscreen mode Exit fullscreen mode

Fancy indexing

Julia arrays can be indexed with slices, but also with boolean or
integer arrays (**masks**). This method is called *fancy indexing*.
It creates **copies not views**.
Enter fullscreen mode Exit fullscreen mode

Using boolean masks

    julia> a = rand(1:21,15)
    julia> a
    array([ 3, 13, 12, 10, 10, 10, 18,  4,  8,  5,  6, 11, 12, 17,  3])
    julia> (a .% 3 .== 0)
    Bool[1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1]
    julia> mask = (a .% 3 .== 0)
    julia> extract_from_a = a[mask] # or,   a[a .% 3 .== 0]
    julia> extract_from_a           # extract a sub-array with the mask
    [ 3, 12, 18,  6, 12,  3]
Enter fullscreen mode Exit fullscreen mode

Indexing with a mask can be very useful to assign a new value to a sub-array:

    julia> a[a .% 3 .== 0] .= -1
    julia> a
    [-1, 13, -1, 10, 10, 10, -1, 4, 8, 5, -1, 11, -1, 17, -1]
Enter fullscreen mode Exit fullscreen mode

Indexing with an array of integers


    julia> a = collect(0:10:90)
    julia> a
    [ 0, 10, 20, 30, 40, 50, 60, 70, 80, 90]
Enter fullscreen mode Exit fullscreen mode

Indexing can be done with an array of integers, where the same index is repeated several times:

    julia> a[[1,2,3]]
    3-element Vector{Int64}:
      0
     10
     20
Enter fullscreen mode Exit fullscreen mode

New values can be assigned with this kind of indexing:

    julia> a[[10, 8]] .= -100
    julia> a
           [0,   10,   20,   30,   40,   50,   60, -100,   80, -100]
Enter fullscreen mode Exit fullscreen mode

Adapted from https://github.com/scipy-lectures/scientific-python-lectures/blob/main/intro/numpy/array_object.rst

Top comments (1)

Collapse
 
Sloan, the sloth mascot
Comment deleted

Some comments have been hidden by the post's author - find out more