Julia Community 🟣: Vinod V

Understanding Cooley-Tukey FFT Matrix Factorization with Julia

Vinod V — Sat, 09 Dec 2023 05:34:32 +0000

#Ref : https://www.cs.cornell.edu/~bindel/class/cs5220-s10/slides/FFT.pdf
#page 17
#run(`$(Base.julia_cmd())`)

using LinearAlgebra

function getBitReversed(N)
    m = Int64(log2(N))
    return [parse(Int64,reverse(bitstring(i)[end-m+1:end]),base = 2) for i in 0:N-1]
end

function getMatrices(N)
#N = 4

    t = Int64(log2(N))
    As = []
    for q = 1:t
        L = 2^q
        Lby2 = 2^(q-1)
        r = Int64(N/L)
        I_r = Matrix(1.0I,r,r)
        tw = exp.(-2im*pi*(0:Lby2-1)/L)
        D = diagm(tw)
        B_r = [Matrix(1.0I,Lby2,Lby2) D;Matrix(1.0I,Lby2,Lby2) -D]
        A_q = kron(I_r,B_r)
        push!(As,A_q)
    end
    return As
end

function test()

    N = 16

    P = 1im*zeros(N,N)

    ind = getBitReversed(N) .+ 1

    for i in Base.OneTo(N)
        P[i,ind[i]] = 1.0
    end

    As = getMatrices(N)

    F_ct = copy(P)

    for i in Base.OneTo(length(As))
        #println("A",i)
        #display(As[i])
        F_ct = As[i]*F_ct
    end

    F_N = exp.(-2im*pi*(0:N-1)*(0:N-1)'/N)

    #=
    println("F_n")
    display(F_N)

    println("F_ct")
    display(F_ct)
    =#

    println(isapprox(F_N,F_ct))
    return F_N,As,P
end

F_N,As,P = test();

Differences between Julia vector and Python/NumPy vector

Vinod V — Mon, 20 Nov 2023 06:06:18 +0000

Let us try to understand a Julia vector [commented code]

julia> using LinearAlgebra             
julia> x = [1,2,3] #x is a Julia vector                                                                         
3-element Vector{Int64}:                                                                    
 1                                                                                          
 2                                                                                          
 3                                                                                          
julia> size(x)                                                                              
(3,)
julia> A = rand(3,3)                                                                        
3×3 Matrix{Float64}:                                                                        
 0.46205    0.907979  0.57372                                                               
 0.0121377  0.590446  0.376092                                                              
 0.672652   0.909502  0.304637
julia> A*x       #x acts as a column vecotor                                                                           
3-element Vector{Float64}:                                                                  
 3.999167421183575                                                                          
 2.3213037773354586                                                                         
 3.4055688127482697
julia> dot(A[1,:],x)    #First element of A*x is the dot product between first row of A and x                                                                    
3.999167421183575 
julia> x*A #Julia thinks x is a matrix with dimension (3,1)                                                                                 
ERROR: DimensionMismatch: matrix A has dimensions (3,1), matrix B has dimensions (3,3)

Let us try to make out the behaviour of a vector in python/numpy

>>> import numpy as np
>>> x = np.array([1,2,3])
>>> x.shape
(3,)
#x is an orientationless array 
>>> A = np.random.random((3,3)) #A is a 3 x 3 matrix
>>> A
array([[0.62256804, 0.48535832, 0.86000692],
       [0.24275571, 0.22393731, 0.5531013 ],
       [0.54999826, 0.90498169, 0.35720945]])
>>> A@x #Left mulitply x with matrix A
array([4.17330543, 2.34993421, 3.43158999]) #x acts as a column vector
>>> np.dot(A[0,:],x) #First element of A@x is the dot product between first row of A and x
4.173305432498873
>>> x@A
array([2.75807425, 3.64817801, 3.03783785])
>>> np.dot(A[:,0],x)
2.7580742486325107 #First element of x@A is the dot product between first column of A and x
>>>

Julia seems to be very particular about the orientation of julia vector while defining.
But Python treats a vector as both column and row vector.

A row vector in Julia can be made as

x = [1 2 3] #without commas in between elements
1×3 Matrix{Int64}:
 1  2  3

A row vector in Python can be made as

x = np.array([[1, 2, 3]])

A column vector in Python can be made as

x = np.array([[1], [2], [3]])

Gram–Schmidt Orthogonalization with Julia

Vinod V — Fri, 27 Oct 2023 04:53:50 +0000

Theorem : Gram–Schmidt Let $W$ be a $p$ -dimensional subspace of $R^n$ , and let $w1,w2,…,wp{w_1, w_2, \dots , w_p}$ be any basis for $W$ . Then the set of vectors ${u_1, u_2, . . . , u_p}$ is an orthogonal basis for $W$ , where

$u_1 = w_1$

$u2=w2−u1Tw2u1Tu1u1u_2 = w_2 - \frac{u_1^Tw_2}{u_1^Tu_1}u_1$

$u3=w3−u1Tw3u1Tu1u1−u2Tw3u2Tu2u2u_3 = w_3 - \frac{u_1^Tw_3}{u_1^Tu_1}u_1 - \frac{u_2^Tw_3}{u_2^Tu_2}u_2$

and where in general

$ui=wi−∑k=1k=i−1ukTwiukTukuk,2≤i≤pu_i = w_i - \sum_{k=1}^{k = i-1}\frac{u_k^Tw_i}{u_k^Tu_k}u_k, 2 \le i \le p$

By rewriting the above equations for $w_i$ ,

$wi=ui−∑k=1k=i−1ukTwiukTukuk,2≤i≤pw_i = u_i - \sum_{k=1}^{k = i-1}\frac{u_k^Tw_i}{u_k^Tu_k}u_k, 2 \le i \le p$

In matrix form

$[w1T⋮wnT]=[100⋯0t2110⋯0t31t321⋯0⋮⋮⋮⋮tn1tn2tn3⋯1][u1T⋮unT]\begin{bmatrix} w_1^T \\ \vdots \\ w_n^T\end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ t_{21} & 1 & 0 & \cdots & 0 \\ t_{31} & t_{32} & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots \\ t_{n1} & t_{n2} & t_{n3} & \cdots & 1 \end{bmatrix} \begin{bmatrix} u_1^T \\ \vdots \\ u_n^T\end{bmatrix}$

$tij=ujTwi−1/(ujTuj),∀(i>j);tij=1,∀(i=j);tij=0,∀(i<j)t_{ij} = u_j^T w_{i-1}/(u_j^T u_j), \forall (i > j); t_{ij} = 1,\forall (i = j);t_{ij} = 0, \forall (i < j)$

Given below is an implementation of Gram-Schmidt orthogonalization in Julia

using LinearAlgebra

N = 10;
A = rand(N,N); #row vectors to be orthogonalized
U = zeros(N,N);
U[1,:] = A[1,:]; 
V = zeros(N,N); # for the for loop
V[1,:] = A[1,:]; 

U[2,:] = A[2,:] - dot(U[1,:],A[2,:])*U[1,:]/dot(U[1,:],U[1,:])

U[3,:] = A[3,:] - dot(U[1,:],A[3,:])*U[1,:]/dot(U[1,:],U[1,:]) - dot(U[2,:],A[3,:])*U[2,:]/dot(U[2,:],U[2,:])

println(dot(U[1,:],U[2,:]))

println(dot(U[1,:],U[3,:]))

println(dot(U[2,:],U[3,:]))

#for loop for a matrix of size N
for i in 2:N
    V[i,:] = A[i,:]
    for j in 1:(i-1)
        V[i,:] -= dot(V[j,:],A[i,:])*V[j,:]/dot(V[j,:],V[j,:])
    end
end


println(dot(V[1,:],V[2,:]))

println(dot(V[1,:],V[3,:]))

println(dot(V[2,:],V[3,:]))

println(dot(V[N,:],V[N-1,:]))

Using row operations

using LinearAlgebra

function test()

N = 4;A = rand(N,N);B = Matrix(1.0I,N,N);U = zeros(N,N);U[1,:] = A[1,:]; V = zeros(N,N);V[1,:] = A[1,:];

row_operations = [] 

B[2,1] = -dot(U[1,:],A[2,:])/dot(U[1,:],U[1,:])

push!(row_operations,B)

X = B*A

B = Matrix(1.0I,N,N)

B[3,1] = -dot(X[1,:],A[3,:])/dot(X[1,:],X[1,:])

B[3,2] = -dot(X[2,:],A[3,:])/dot(X[2,:],X[2,:])

push!(row_operations,B)

X = B*X

println(dot(X[1,:],X[2,:]))

println(dot(X[1,:],X[3,:]))

println(dot(X[2,:],X[3,:]))

X = zeros(N,N);X = copy(A);#X[1,:] = A[1,:];
row_operations = []
for i in 2:N
    C = Matrix(1.0I,N,N)
    for j in 1:(i-1)
        C[i,j] = -dot(X[j,:],A[i,:])/dot(X[j,:],X[j,:])     
    end
    push!(row_operations,C)
    X = C*X

end

println(dot(X[1,:],X[2,:]))

println(dot(X[1,:],X[3,:]))

println(dot(X[2,:],X[3,:]))

println(dot(X[N,:],X[N-1,:]))

end

test()

Julia Vectors and Matrices

Vinod V — Tue, 10 Oct 2023 10:28:16 +0000

What are Julia vectors?

Julia vectors

Julia provides:

- Builtin multi-dimensional arrays

- closer to hardware (efficiency)

- designed for scientific computation (convenience)

    julia> v = [0, 1, 2, 3]
      4-element Vector{Int64}:     
       0                           
       1                           
       2                           
       3

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)

Julia supports multidimensional arrays

Julia Reference documentation

On the web: https://docs.julialang.org/
Interactive help:

     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.

Looking for something:

     julia> mat<Tab><Tab> # press tab key two times
     MathConstants  Matrix

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

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

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?

Functions for creating vectors

In practice, we rarely enter items one by one...

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]

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

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

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

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

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.

The default data type is floating point:

    julia> a = ones(3, 3)
    julia> eltype(a)
    Float64

There are also other types:

Complex:

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

Bool:

    julia> e = [true, false, false, true]
    julia> eltype(e)
    Bool

Strings:

    julia> f = ["Bonjour", "Hello", "Hallo"]
    julia> eltype(f)
    String

Much more:

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

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

    julia> using PyPlot
    julia> PyPlot.plot(x, y)       # line plot    # doctest: +SKIP

1D plotting:

  julia> x = LinRange(0, 3, 20)
  julia> y = LinRange(0, 9, 20)
  julia> PyPlot.plot(x, y, "g^")  # ^ plot

(img)[]

2D arrays (such as images):


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

(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)

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])

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

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]

Note that the last index is included! :

    julia> a[1:4]
    array([0, 1, 2, 3])


    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]

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

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.]



    *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]

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

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)

Cross out 0 and 1 which are not primes:

   julia> is_prime[1] = false

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

Follow-up:
- Move the above code into a script file named prime_sieve.jl
- Run it to check it works
- Use the optimization suggested in the sieve of Eratosthenes https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes

  1. Skip ``j`` which are already known to not be primes

  2. The first number to cross out is :math:`j^2`

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**.

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]

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]

Indexing with an array of integers


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

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

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]

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

Exploring Quadratic forms with Julia

Vinod V — Fri, 22 Sep 2023 09:48:19 +0000

A quadratic form is a type of polynomial in many variables where all the terms have an exponent of two.

Definition of quadratic form : A quadratic form in $n$ variables $x1,…,xnx_1,\dots,x_n$ is given by the expression $\sum_{i=1}^{i = n}\sum_{j=1}^{i = n}b_{ij}x_ix_j$ . The quadratic form $q (x)$ can also be expressed in terms of matrix vector product as $q(x) = x^T A x$ where $[x_1,\dots,x_n]^T$ and $aij=bij+bji2a_{ij} = \frac{b_{ij} + b_{ji}}{2}$ . Obviously $A$ is a symmetric matrix.

Is matrix $A$ unique for the quadratic form? It seems to be yes.

Hypothesis $H$ : Let $r(x) = x^TCx$ be a quadratic form with asymmmetric matrix $C$ . Then $r(x) = x^TBx$ , where $\frac{C+C^T}{2}$ The following Julia code tries to take out any exceptions to the hypothesis.

using Symbolics
@variables x y z
for i in 1:10
         C = rand(1:10,3,3)
         B = (A+A')/2
         q = [x y z]*C*[x y z]'
         eq = simplify(q[1],expand=true)
         q1 = [x y z]*B*[x y z]'
         eq1 = simplify(q1[1],expand=true)
         println(isequal(eq,eq1)) # always true
end

Hypothesis $H$ seems to be true.

Let $n$ be the number of variables in a quadratic form $r (x)$ . Then the number of terms $p$ in $q (x)$ is $n(n+1)2\frac{n(n+1)}{2}$ . If $p$ is given then $\sqrt(1+8p))/2)$ .

Julia function to calculate $n$ given $p$ can be f(p) = Int64(trunc((-1 + sqrt(1+8p))/2)).

A quadratic expression $q (x)$ is

Positive definite if $\forall x \in R^n$ .
Positive semi-definite if $\ge 0, \forall x \in R^n.$
Negative definite if $\forall x \in R^n$ .
Negative semi-definite if $\le 0, \forall x \in R^n.$
Indefinite if $q (x) > 0$ , for some $\in R^n$ and $q (x) < 0$ , for some $\in R^n$

Definiteness of a matrix $A$ can be characterized by eigen values $λ1,…,λn\lambda_1,\dots,\lambda_n$ of $A$ as

Positive definite if $λi>0,∀i\lambda_i > 0, \forall i$ .
Positive semi-definite if $λi≥0,∀i\lambda_i \ge 0, \forall i$ .
Negative definite if $λi<0,∀i\lambda_i < 0, \forall i$ .
Negative semi-definite if $λi≤0,∀i.\lambda_i \le 0, \forall i.$
Indefinite if $λi>0\lambda_i > 0$ , for some $i$ and $λj<0\lambda_j < 0$ , for some $\neq j$

The following code in Julia detemines the definiteness of a matrix $A$ in terms of eigenvalues.

using LinearAlgebra
function getDefiniteness(A)

    #@assert isequal(A,A') "Matrix is not symmetric"
ishermitian(A) || throw(ArgumentError("Matrix must be Hermitian"))

    if isposdef(A)
        return "Matrix is positive definite"
    end 

    if isposdef(-A)
        return "Matrix is negative definite"
    end

    #x,y = eigmin(A),eigmax(A)
    z = eigvals(A)
    x,y = z[1],z[end]

    if iszero(x)
        return "Matrix is positive semi definite"
    end

    if y > 0.0 && x < 0.0
        return "Matrix is indefinite"
    end

    if iszero(y)
        return "Matrix is negative semi definite"
    end

    #=
    if x < 0 && y < 0 #all([x, y] .< 0)
        return "Matrix is negative definite"
    end
    =#

end

A = rand(5,5)
A = (A+A')/2

A = [-1 0;0 -1.0]
getDefiniteness(A)

By diagonalizing the symmetric matrix associated with a quadratic form, it is possible to eliminate cross terms.

The Principal Axes Theorem is a mathematical concept related to quadratic forms, which states that a quadratic form $q (x)$ given by $q(x) = x^T A x$ can be diagonalized by an orthonormal set of eigenvectors of $A$ , which is called the set of principal axes. The theorem is based on the spectral theorem related to a symmetric matrix, which states that $Q^T A Q = D$ , where $D$ is a diagonal matrix containing the eigenvalues of $A$ in its diagonal and $Q$ is the eigen matrix of $A$ . By substituting $y = Q x$ , the quadratic form can be expressed as $y^T D y = \lambda_1 y_1^2 + \cdots + \lambda_n y_n^2$ , where $λi\lambda_i$ are the eigenvalues of $A$ .

The minimum and maximum of a quadratic form subjected to the constraint $x12+⋯+xn2=1x_1^2 + \cdots + x_n^2 = 1$ can be calculated by principal axes theorem.

Let $q(x) = x^T Ax, ||x|| = 1$ .

Minimum of $\lambda_1$ is the minimum of eigenvalues of $A$ .

Maximum of $\lambda_n$ is the maximum of eigenvalues of $A$ .

julia> using  LinearAlgebra

julia> A = rand(1:10,4,4)
4×4 Matrix{Int64}:
 5   6  2   7
 9   8  5   4
 2   7  5  10
 9  10  8   2

julia> A = (A+A')/2
4×4 Matrix{Float64}:
 5.0  7.5  2.0  8.0
 7.5  8.0  6.0  7.0
 2.0  6.0  5.0  9.0
 8.0  7.0  9.0  2.0

julia> minim,maxim = eigmin(A),eigmax(A)
(-7.833772750196167, 25.006216466801746)

function eigminmax(H::Hermitian) #suggested by stevenj
           T = hessenberg(H).H # real SymTridiagonal reduction
           return eigmin(T), eigmax(T)
end

julia> A = Hermitian(rand(100,100));

julia> @btime a,b = eigminmax(A)
  920.200 μs (32 allocations: 145.23 KiB)
(-5.668761039486605, 49.66411909262399)

julia> @btime x,y = eigmin(A),eigmax(A)
  1.602 ms (25 allocations: 229.81 KiB)
(-5.668761039486608, 49.66411909262401)

Discussion here : https://discourse.julialang.org/t/exploring-quadratic-forms-with-julia/104138

Boosting Julia Performance: Benchmarking Code in Function, Global Scope, and Builtin

Vinod V — Wed, 20 Sep 2023 10:22:34 +0000

Julia's code execution efficiency can vary depending on where the code is executed. The two possible locations are the local scope and the global scope.Alternatively, we can use the builtin function. To compare the efficiency of code execution in these two locations, we can use the BenchmarkTools module in Julia. For instance, consider the following code snippet that calculates the total of an array in global , local scopes and builtin. The profiling results indicate that the local scope code is more efficient in terms of execution time than the global scope code But the builtin code is most efficient (probably Julia is using the formula n(n+1)/2 to sum first n integers).

using BenchmarkTools

i = 1
total = 0

function sum1(num) #global scope
global i
global total
    for i in eachindex(num)
        total += num[i]
    end
    return total
end

function sum2(num) #local scope
    total = 0
    for i in eachindex(num)
        total += num[i]
    end
    return total
end

#a = rand(1:1000,1,10000)
a = 1:100000

println("sum1")
@btime b = sum1($a)
println("Sum2")
@btime b = sum2($a)
println("Builtin sum")
@btime b = sum($a)

sum1
  3.136 ms (199489 allocations: 3.04 MiB)
Sum2
  55.498 μs (0 allocations: 0 bytes)
Builtin sum
  0.001 ns (0 allocations: 0 bytes)
5000050000

Discussion here https://discourse.julialang.org/t/why-the-function-sum1-is-faster-than-builtin-sum/104051/12

Control Flow in Julia

Vinod V — Fri, 15 Sep 2023 06:57:03 +0000

Control Flow

In the programs we have seen till now, there has always been a series of statements faithfully executed by Julia in exact top-down order. What if you wanted to change the flow of how it works? For example, you want the program to take some decisions and do different things depending on different situations, such as printing "Good Morning" or "Good Evening" depending on the time of the day?

As you might have guessed, this is achieved using control flow statements. There are three control flow statements in Julia - if, for and while.

The `if` statement

The if statement is used to check a condition: if the condition is true, we run a block of statements (called the if-block), else we process another block of statements (called the else-block). The else clause is optional.

Example (save as if.jl):

number = 23
println("Enter an integer : ")
text = readline()
guess = parse(Int,text)

if guess == number
    # New block starts here
    println("Congratulations, you guessed it.")
    println("(but you do not win any prizes!)")
    # New block ends here
elseif guess < number
    # Another block
    println("No, it is a little higher than that")
    # You can do whatever you want in a block ...
else
    println("No, it is a little lower than that")
    # you must have guessed > number to reach here
end

println("Done")
# This last statement is always executed,
# after the if statement is executed.

Output:

Enter an integer :                     
6                                      
No, it is a little higher than that   
Done

How It Works
In this program, we take guesses from the user and check if it is the number that we have. We set the variable number to any integer we want, say 23. Then, we take the user's guess using the readline() function. Functions are just reusable pieces of programs. We'll read more about them in the next chapter.

We supply a string to the built-in readline function which prints it to the screen and waits for input from the user. Once we enter something and press 'enter' key, the readline() function returns what we entered, as a string. We then convert this string to an integer using parse and then store it in the variable guess. Actually, the Int is a data type but all you need to know right now is that you can use it to convert a string to an integer (assuming the string contains a valid integer in the text).

Next, we compare the guess of the user with the number we have chosen. If they are equal, we print a success message. Notice that we use indentation levels to tell Julia which statements belong to which block.

Then, we check if the guess is less than the number, and if so, we inform the user that they must guess a little higher than that. What we have used here is the elseif clause which actually combines two related if else-if else statements into one combined if-elseif-else statement. This makes the program easier and reduces the amount of indentation required.

The elseif and else statements must also have a colon at the end of the logical line followed by their corresponding block of statements (with proper indentation, of course)

You can have another if statement inside the if-block of an if statement and so on - this is called a nested if statement.

Remember that the elseif and else parts are optional. A minimal valid if statement is:

if true
    print("Yes, it is true")
end

After Julia has finished executing the complete if statement along with the associated elseif and else clauses, it moves on to the next statement in the block containing the if statement. In this case, it is the main block (where execution of the program starts), and the next statement is the print('Done') statement. After this, Julia sees the ends of the program and simply finishes up.

Even though this is a very simple program, I have been pointing out a lot of things that you should notice. All these are pretty straightforward (and surprisingly simple for those of you from C/C++ backgrounds). You will need to become aware of all these things initially, but after some practice you will become comfortable with them, and it will all feel 'natural' to you.

Note for C/C++ Programmers

There is no switch statement in Julia. You can use an if..elseif..else statement to do the same thing (and in some cases, use a Dict to do it quickly). Dict{K,V}() constructs a hash table with keys of type K and values of type V which will be discussed later.

The while Statement

The while statement allows you to repeatedly execute a block of statements as long as a condition is true. A while statement is an example of what is called a looping statement. A while statement can have an optional else clause.

Example (save as while.jl):

number = 23
running = true

while running == true
    text1 = readline()
    guess = parse(Int,text1)

    if guess == number
        print("Congratulations, you guessed it.")
        # this causes the while loop to stop
        running = false
    elseif guess < number
        print("No, it is a little higher than that.")
    else
        print("No, it is a little lower than that.")
    end

end

Output:

23
Congratulations, you guessed it.

How It Works

In this program, we are still playing the guessing game, but the advantage is that the user is allowed to keep guessing until he guesses correctly - there is no need to repeatedly run the program for each guess, as we have done in the previous section. This aptly demonstrates the use of the while statement.

We move the input and if statements to inside the while loop and set the variable running to true before the while loop. First, we check if the variable running is true and then proceed to execute the corresponding while-block. After this block is executed, the condition is again checked which in this case is the running variable. If it is true, we execute the while-block again, else we continue to execute the optional else-block and then continue to the next statement.

The else block is executed when the while loop condition becomes false - this may even be the first time that the condition is checked. If there is an else clause for a while loop, it is always executed unless you break out of the loop with a break statement.

The true and false are called Boolean types and you can consider them to be equivalent to the value 1 and 0 respectively.

Note for C/C++/Python Programmers

Remember that you can have an else clause for the while loop.

The `for` loop

The for..in statement is another looping statement which iterates over a sequence of objects i.e. go through each item in a sequence. We will see more about sequences in detail in later chapters. What you need to know right now is that a sequence is just an ordered collection of items.

Example (save as for.jl):

for i in 1:4:
    println(i)

Output:

How It Works

In this program, we are printing a sequence of numbers. We generate this sequence of numbers using the built-in : operator.

m:n returns a sequence of numbers starting from the first number and up to the second number. For example, 1:4 gives the sequence [1, 2, 3, 4]. By default, : takes a step count of 1. If we supply a third number to :, then that becomes the step count. For example, 1:2:7 gives [1,3,5,7]. Remember that the range extends to the second number i.e. it does include the second number.

Note that : generates only one number at a time, if you want the full list of numbers, call collect() on the :, for example, collect(1:4) will result in [1, 2, 3, 4].

The for loop then iterates over this range - for i in 1:4 is equivalent to for i = [1, 2, 3, 4] which is like assigning each number (or object) in the sequence to i, one at a time, and then executing the block of statements for each value of i. In this case, we just print the value in the block of statements.

Remember that the for..in loop works for any sequence. Here, we have a list of numbers generated by the built-in : operator, but in general we can use any kind of sequence of any kind of objects! We will explore this idea in detail in later chapters.

Note for C/C++ Programmers

In C/C++, if you want to write for (int i = 0; i < 5; i++), then in Julia you write just for i in 0:4. As you can see, the for loop is simpler, more expressive and less error prone in Julia.

The foreach statement

The foreach statement is an inbuilt function in Julia that is used to call a specified set of instructions on each element of a specified iterable. The syntax for the foreach statement is foreach(f, c...), where f is the specified set of instructions and c is the specified iterable. The function f is called element-wise, and iteration stops when any iterator is finished. The foreach statement should be used instead of map when the results of f are not needed. Below are some code examples of the foreach statement in Julia:

julia> foreach(println,[1,2,3])        
1                                      
2                                      
3                                      
  julia> tri = 1:3:7; res = Int[];

  julia> foreach(x -> push!(res, x^2), tri)

  julia> res
  3-element Vector{Int64}:
    1
   16
   49
julia> foreach((x, y) -> println(x, " with ", y), 1:3:7, 'a':'c')

The break Statement

The break statement is used to break out of a loop statement i.e. stop the execution of a looping statement, even if the loop condition has not become false or the sequence of items has not been completely iterated over.

Example (save as break.jl):

while true
    println("Enter something : ")

    s = readline()    

    if s == "quit"
        break
    end
    println("Length of the string is ", length(s))

end

Output:

Enter something : 
quit

How It Works

In this program, we repeatedly take the user's input and print the length of each input each time. We are providing a special condition to stop the program by checking if the user input is "quit". We stop the program by breaking out of the loop and reach the end of the program.

The length of the input string can be found out using the built-in length function.

Remember that the break statement can be used with the for loop as well.

The `continue` Statement

The continue statement is used to tell Julia to skip the rest of the statements in the current loop block and to continue to the next iteration of the loop.

Example (save as continue.jl):

while true
    println("Enter something : ")
    s = readline()
    if s == "quit"
        break
    end
    if length(s) < 3
        println("Too small")
        continue
    end
    println("Input is of sufficient length")    
end

Output:

Enter something :                    
2                                    
Too small
Enter something :           
quit

                             **How It Works**

In this program, we accept input from the user, but we process the input string only if it is at least 3 characters long. So, we use the built-in length function to get the length and if the length is less than 3, we skip the rest of the statements in the block by using the continue statement. Otherwise, the rest of the statements in the loop are executed, doing any kind of processing we want to do here.

Note that the continue statement works with the for loop as well.

Summary

We have seen how to use the three control flow statements - if, while and for along with their associated break and continue statements. These are some of the most commonly used parts of Julia and hence, becoming comfortable with them is essential.

Adapted from control_flow

Master Julia Dictionaries with These Easy Recipes

Vinod V — Tue, 12 Sep 2023 10:17:51 +0000

Find out what two dictionaries have in common

a = Dict("x" => 1, "y" => 2, "z" => 3)

b = Dict("w" => 10, "x" => 11, "y" => 2)

println("Common keys:", intersect(keys(a), keys(b)))
println("Keys in a not in b:", setdiff(keys(a), keys(b)))
println("(key,value) pairs in common:", intersect(collect(a), collect(b)))

Output:
Common keys:Set(["x", "y"])
Keys in a not in b:Set(["z"])
(key,value) pairs in common:["y" => 2]

Extracting a subset from a dictionary

prices = Dict(
   "ACME" => 46.23,
   "AAPL" => 612.78,
   "IBM" => 205.55,
   "HPQ" => 37.20,
   "FB" => 10.75,
   "INTL" => 30.56,
   "GS" => 322.3
)

# Make a dictionary of all prices over 100
p100 = Dict([ (key, value) for (key, value) in prices if value > 100 ])

println("All prices over 100")
println(p100)

# Make a dictionary of tech stocks
tech_shares = Set(["AAPL", "IBM", "HPQ", "MSFT","INTL"])
ptech = Dict([ (key, value) for (key, value) in prices if key in tech_shares ])

println("All techs")
println(ptech)

Output:
All prices over 100
Dict("IBM" => 205.55, "AAPL" => 612.78, "GS" => 322.3)
All techs
Dict("HPQ" => 37.2, "IBM" => 205.55, "AAPL" => 612.78, "INTL" => 30.56)

Merging dictionaries in Julia

dict1 = Dict("a" => 1, "b" => 2, "c" => 3)
dict2 = Dict("c" =>10, "e" => 5, "f" => 6)
dict3 = merge(dict1,dict2)
print(dict3)
Output:
Dict("f" => 6, "c" => 10, "e" => 5, "b" => 2, "a" => 1)

Find the N smallest or largest items

portfolio = [
   Dict("name"=> "IBM", "shares"=> 100, "price"=> 91.2),
   Dict("name"=> "AAPL", "shares"=> 50, "price"=> 643.22),
   Dict("name"=> "FB", "shares"=> 150, "price"=> 21.09),
   Dict("name"=> "HPQ", "shares"=> 35, "price"=> 31.75),
   Dict("name"=> "YHOO", "shares"=> 45, "price"=> 16.35),
   Dict("name"=> "ACME", "shares"=> 75, "price"=> 115.65)
]
N = 3
sort(portfolio, by = x -> x["price"])[1:N] #N smallest shares by price
sort(portfolio, by = x -> x["price"],rev=true)[1:N] #N largest by price

Sort a dictionary to find minimum and maximum values

prices = Dict("ACME" => 45.23, "AAPL" => 612.78, "IBM" => 205.55, "HPQ" => 37.20, "FB" => 10.75)

# Find min and max price
min_price = findmin(prices)
max_price = findmax(prices)

println("min price: ", min_price)
println("max price: ", max_price)

println("sorted prices:")
prices_sorted = sort(collect(prices), by=x->x[2])
for (name, price) in prices_sorted
    println("    ", name, " ", price)
end

Output:
min price: (10.75, "FB")
max price: (612.78, "AAPL")
sorted prices:
    FB 10.75
    HPQ 37.2
    ACME 45.23
    IBM 205.55
    AAPL 612.78

Sort a dictionary based on two keys simultaneously

rows = [Dict("fname" =>  "Mohandas", "lname" =>  "Gandhi", "uid" =>  1003),
    Dict("fname" =>  "Jawaharlal", "lname" =>  "Nehru", "uid" =>  1002),
    Dict("fname" =>  "Narasimha", "lname" =>  "Rao", "uid" =>  1001),
    Dict("fname" =>  "Arvind", "lname" =>  "Kejrival", "uid" =>  1004)]



# sort list of Dicts by two keys simultaneously
a = sort(rows, by = x -> (x["lname"], x["fname"]))

print(a)

Output:
Dict("lname" => "Gandhi", "uid" => 1003, "fname" => "Mohandas")
 Dict("lname" => "Kejrival", "uid" => 1004, "fname" => "Arvind")
 Dict("lname" => "Nehru", "uid" => 1002, "fname" => "Jawaharlal")
 Dict("lname" => "Rao", "uid" => 1001, "fname" => "Narasimha")

Operators and Expressions in Julia

Vinod V — Mon, 11 Sep 2023 05:30:41 +0000

Operators and Expressions

Most statements (logical lines) that you write will contain expressions. A simple example of an expression is 2 + 3. An expression can be broken down into operators and operands.

Operators are functionality that do something and can be represented by symbols such as + or by special keywords. Operators require some data to operate on and such data is called operands. In this case, 2 and 3 are the operands.

Operators

We will briefly take a look at the operators and their usage.

Note that you can evaluate the expressions given in the examples using the interpreter interactively. For example, to test the expression 2 + 3, use the interactive Julia interpreter prompt:

julia> 2 + 3
5
julia> 3 * 5
15
julia>

Here is a quick overview of the available operators:

+ (plus)
- Adds two objects
- 3 + 5 gives 8.
- (minus)
- Gives the subtraction of one number from the other; if the first operand is absent it is assumed to be zero.
- -5.2 gives a negative number and 50 - 24 gives 26.
* (multiply)
- Gives the multiplication of the two numbers or returns the string repeated that many times.
- 2 * 3 gives 6. "ab" ^ 3 gives "ababab".
^ (power)
- Returns x to the power of y
- 3 ^ 4 gives 81 (i.e. 3 * 3 * 3 * 3)
/ (divide)
- Divide x by y
- 13 / 3 gives 4.333333333333333

div (divide and floor)

Divide x by y and round the answer down to the nearest integer value towards zero. Note that if one of the values is a float, you'll get back a float.
div(13,3) gives 4
div(-13,3) gives -4

- `div(9,1.81)` gives `4.0`

  #from julia documentation for div
  julia> div(4, 3, RoundDown) 
  1
  julia> div(4, 3, RoundUp) 
  2
  julia> div(5, 2, RoundNearest)
  2
  julia> div(5, 2, RoundNearestTiesAway)
  3
  julia> div(-5, 2, RoundNearest)
  -2
  julia> div(-5, 2, RoundNearestTiesAway)
  -3
  julia> div(-5, 2, RoundNearestTiesUp)
  -2
  julia> div(4, 3, RoundFromZero)
  2
  julia> div(-4, 3, RoundFromZero)
  -2

// (Divide two integers or rational numbers, giving a Rational result.)
- (5//2) //(2//5) gives 25//4
% (modulo)
- Returns the remainder of the division
- 13 % 3 gives 1. -25.5 % 2.25 gives -0.75.
<< (left shift)
- Shifts the bits of the number to the left by the number of bits specified. (Each number is represented in memory by bits or binary digits i.e. 0 and 1)
- 2 << 2 gives 8. 2 is represented by 10 in bits.
- Left shifting by 2 bits gives 1000 which represents the decimal 8.
>> (right shift)
- Shifts the bits of the number to the right by the number of bits specified.
- 11 >> 1 gives 5.
- 11 is represented in bits by 1011 which when right shifted by 1 bit gives 101which is the decimal 5.
& (bit-wise AND)
- Bit-wise AND of the numbers: if both bits are 1, the result is 1. Otherwise, it's 0.
- 5 & 3 gives 1 (0101 & 0011 gives 0001)
| (bit-wise OR)
- Bitwise OR of the numbers: if both bits are 0, the result is 0. Otherwise, it's 1.
- 5 | 3 gives 7 (0101 | 0011 gives 0111)
\xor<tab> (bit-wise XOR)
- Bitwise XOR of the numbers: if both bits (1 or 0) are the same, the result is 0. Otherwise, it's 1.
- 5 \xor<tab> 3 gives 6 (O101 xor 0011 gives 0110)
~ (bit-wise invert)
- The bit-wise inversion of x is -(x+1)
- ~5 gives -6.
< (less than)
- Returns whether x is less than y. All comparison operators return true or false. Note the small case of these names.
- 5 < 3 gives false and 3 < 5 gives true.
- Comparisons can be chained arbitrarily: 3 < 5 < 7 gives true.
> (greater than)
- Returns whether x is greater than y
- 5 > 3 returns true. Comparisons can be chained arbitrarily: 3 > 5 > 7 gives false.
<= (less than or equal to)
- Returns whether x is less than or equal to y
- x = 3; y = 6; x <= y returns true
>= (greater than or equal to)
- Returns whether x is greater than or equal to y
- x = 4; y = 3; x >= 3 returns true
== (equal to)
- Compares if the objects are equal
- x = 2; y = 2; x == y returns true
- x = "str"; y = "stR"; x == y returns false
- x = 'str'; y = 'str'; x == y returns true
!= (not equal to)
- Compares if the objects are not equal
- x = 2; y = 3; x != y returns true
! (boolean NOT)
- If x is true, it returns false. If x is false, it returns true.
- x = true; !x returns false.
&& (boolean AND)
- x && y returns false if x is false, else it returns evaluation of y
- x = false; y = true; x && y returns false since x is false. In this case, Julia will not evaluate y since it knows that the left hand side of the 'and' expression is False which implies that the whole expression will be False irrespective of the other values. This is called short-circuit evaluation.
|| (boolean OR)
- If x is true, it returns True, else it returns evaluation of y
- x = true; y = false; x || y returns true. Short-circuit evaluation applies here as well.

Shortcut for math operation and assignment

It is common to run a math operation on a variable and then assign the result of the operation back to the variable, hence there is a shortcut for such expressions:

a = 2
a = a * 3

can be written as:

a = 2
a *= 3

Notice that var = var operation expression becomes var operation= expression.

Evaluation Order

If you had an expression such as 2 + 3 * 4, is the addition done first or the multiplication? Our high school maths tells us that the multiplication should be done first. This means that the multiplication operator has higher precedence than the addition operator.

The following link gives the precedence table for Julia, from the lowest precedence (least binding) to the highest precedence (most binding). This means that in a given expression, Julia will first evaluate the operators and expressions lower in the table before the ones listed higher in the table.

Mathematical Operations and Elementary Functions · The Julia Language

Changing the Order Of Evaluation

To make the expressions more readable, we can use parentheses. For example, 2 + (3 * 4) is definitely easier to understand than 2 + 3 * 4 which requires knowledge of the operator precedences. As with everything else, the parentheses should be used reasonably (do not overdo it) and should not be redundant, as in (2 + (3 * 4)).

There is an additional advantage to using parentheses - it helps us to change the order of evaluation. For example, if you want addition to be evaluated before multiplication in an expression, then you can write something like (2 + 3) * 4.

Operator precedence of an operator > can be found out by the function Base.operator_precedence(:>)

julia> Base.operator_precedence(:^)
13

Expressions

Example (save as expression.jl):

length = 5
breadth = 2

area = length * breadth
print("Area is ", area)
print("Perimeter is ", 2 * (length + breadth))

Output:

$ julia expression.jl
Area is 10
Perimeter is 14

How It Works

The length and breadth of the rectangle are stored in variables by the same name. We use these to calculate the area and perimeter of the rectangle with the help of expressions. We store the result of the expression length * breadth in the variable area and then print it using the print function. In the second case, we directly use the value of the expression 2 * (length + breadth) in the print function.

Summary

We have seen how to use operators, operands and expressions - these are the basic building blocks of any program. Next, we will see how to make use of these in our programs using statements.

Adapted from Operators and Expressions

Julia Basics

Vinod V — Sat, 09 Sep 2023 01:19:40 +0000

Basics

Just printing hello world is not enough, is it? You want to do more than that - you want to take some input, manipulate it and get something out of it. We can achieve this in Julia using constants and variables, and we'll learn some other concepts as well in this chapter.

Comments

Comments are any text to the right of the # symbol and is mainly useful as notes for the reader of the program.

For example:

print("hello world") # Note that print is a function

or:

# Note that print is a function
print("hello world")

Use as many useful comments as you can in your program to:

explain assumptions
explain important decisions
explain important details
explain problems you're trying to solve
explain problems you're trying to overcome in your program, etc.

Code tells you how, comments should tell you why.

This is useful for readers of your program so that they can easily understand what the program is doing. Remember, that person can be yourself after six months!

Literal Constants

An example of a literal constant is a number like 5, 1.23, or a string like "This is a string" or "It's a string!".

It is called a literal because it is literal - you use its value literally. The number 2 always represents itself and nothing else - it is a constant because its value cannot be changed. Hence, all these are referred to as literal constants.

Numbers

Numbers are mainly of two types - integers and floats.

An example of an integer is 2 which is just a whole number.

Examples of floating point numbers (or floats for short) are 3.23 and 52.3E-4. The E notation indicates powers of 10. In this case, 52.3E-4 means 52.3 * 10^-4^.

Strings

A string is a sequence of characters. Strings are basically just a bunch of words.

You will be using strings in almost every Julia program that you write, so pay attention to the following part.

Double Quote

You can specify strings using single quotes such as "Quote me on this".

All white space i.e. spaces and tabs, within the double quotes, are preserved as-is.

You can specify multi-line strings using double quotes. You can use single quotes freely within the double quotes. In case you want to use double quotes precede it with a "\" backslash character. An example is:

"This is a multi-line string. This is the first line.It's a text.
This is the second line.
\"What's your name?,\" I asked.
He said \"Bond, James Bond.\"
"

Strings Are Immutable

This means that once you have created a string, you cannot change it. Although this might seem like a bad thing, it really isn't. We will see why this is not a limitation in the various programs that we see later on.

The format method

Sometimes we may want to construct strings from other information. Julia uses the $ symbol to interpolate the variables directly into the string.

age = 20
name = "Swaroop"

println("$name was $age years old when he wrote this book")
println("Why is $name playing with that Julia?")

Output:

Swaroop was 20 years old when he wrote this book
Why is Swaroop playing with that Julia?

Notice that we could have achieved the same using string concatenation operator '*':

name * " is " * string(age) * " years old"

but that is much uglier and more error-prone. Second, the conversion to string would be done automatically by the $ method instead of the explicit conversion to strings needed in this case. Third, when using the $ method, we can change the message without having to deal with the variables used and vice-versa.

What Julia does in the $ method is that it substitutes each argument value into the place of the specification. There can be more detailed specifications such as:

# decimal (.) precision of 3 for float '1.23456'
using Printf
@printf("%1.3f",1.23456)

Output:

1.234

Escape Sequences

Suppose, you want to have a string which contains a double quote ("), how will you specify this string? For example, the string is "This is a "special" string". You cannot specify as such because Julia will be confused as to where the string starts and ends. So, you will have to specify that this single quote does not indicate the end of the string. This can be done with the help of what is called an escape sequence. You specify the double quote as \" : notice the backslash. Now, you can specify the string as "This is a \"special\" string".

Another useful escape sequence to know is the tab: \t. There are many more escape sequences but I have mentioned only the most useful ones here.

Variable

Using just literal constants can soon become boring - we need some way of storing any information and manipulate them as well. This is where variables come into the picture. Variables are exactly what the name implies - their value can vary, i.e., you can store anything using a variable. Variables are just parts of your computer's memory where you store some information. Unlike literal constants, you need some method of accessing these variables and hence you give them names.

Identifier Naming

Variables are examples of identifiers. Identifiers are names given to identify something. There are some rules you have to follow for naming identifiers:

Generally, variable names are written in lowercase. Underscores are used to separate different words in a variable name, but it is not advisable to use names that would require underscores. Function and macros names are in lowercase. Underscores are not used.The first character of types and modules is uppercase. The separation between words in names is done using upper camel case.The functions that modify or write to any of their arguments end with "!" symbol.

Data Types

Variables can hold values of different types called data types. The basic types are numbers and strings, which we have already discussed. In later chapters, we will see how to create our own types using struct keyword.

We will now see how to use variables along with literal constants. Save the following example and run the program.

How to write Julia programs

Henceforth, the standard procedure to save and run a Julia program is as follows:

Open a text editor (notepad,notepad++,gedit etc)
Create new file with the filename mentioned (file extension must be .jl).
Type the program code given in the example.
run the current file in julia command line as include("path to file")

Example: Using Variables And Literal Constants

Type and run the following program:

# Filename : var.jl
i = 5
print(i)
i = i + 1
print(i)

s = "This is a multi-line string.
This is the second line."
print(s)

Output:

5
6
This is a multi-line string.
This is the second line.

How It Works

Here's how this program works. First, we assign the literal constant value 5 to the variable i using the assignment operator (=). This line is called a statement because it states that something should be done and in this case, we connect the variable name i to the value 5. Next, we print the value of i using the print statement which, unsurprisingly, just prints the value of the variable to the screen.

Then we add 1 to the value stored in i and store it back. We then print it and expectedly, we get the value 6.

Similarly, we assign the literal string to the variable s and then print it.

Note for static language programmers

Variables are used by just assigning them a value. No declaration or data type definition is needed/used.

Logical And Physical Line

A physical line is what you see when you write the program. A logical line is what Python sees as a single statement. Julia implicitly assumes that each physical line corresponds to a logical line.

An example of a logical line is a statement like print('hello world') - if this was on a line by itself (as you see it in an editor), then this also corresponds to a physical line.

Implicitly, Julia encourages the use of a single statement per line which makes code more readable.

If you want to specify more than one logical line on a single physical line, then you have to explicitly specify this using a semicolon (;) which indicates the end of a logical line/statement. For example:

i = 5
print(i)

is effectively same as

i = 5;
print(i);

which is also same as

i = 5; print(i);

and same as

i = 5; print(i)

However, I strongly recommend that you stick to writing a maximum of a single logical line on each single physical line. The idea is that you should never use the semicolon. In fact, I have never used or even seen a semicolon in a Julia program.

Indentation is not important in Julia but indented code tend to be more readable.

This blog post is an adaptation of byte-of-python
/basics.md to Julia.

Filtering Data Made Easy: Tips and Tricks in Julia

Vinod V — Wed, 06 Sep 2023 01:38:58 +0000

data_list = [1, 4, -5, 11, -7, 4, 3, -1]

Values greater than zero

positives = [n for n in data_list if n > 0]
#positives = filter(>(0),data_list)
println(positives)

Values less than zero

negatives = [n for n in data_list if n < 0]
#negatives = filter(<(0),data_list)
println(negatives)

Negative values replaced with 0

neg_to_zero = [n > 0 ? n : 0 for n in data_list]
println(neg_to_zero)

Positive values replaced with 0

pos_to_zero = [n < 0 ? n : 0 for n in data_list]
println(pos_to_zero)

Filter to a range

rang = filter(x-> -3<x<1, data_list) #between -3 and 1
println(rang)

Clamp between low and high values

clamp.(data_list,-2,4)
8-element Vector{Int64}:
  1
  4
 -2
  4
 -2
  4
  3
 -1

Select values from a list based on values on some other list

addresses = [
    "5412 N Alappuzha",
    "5148 N Munnar", 
    "5800 E Kovalam",
    "2122 N Palakkad",
    "5645 N Wayanad",
    "1060 W Kumarakom",
    "4801 N Murukani",
    "1039 W Kollengode",
]

counts = [0, 3, 10, 4, 1, 7, 6, 1]

more5 = [n > 5 for n in counts]
a = addresses[more5]
a = [addresses[i] for i in eachindex(more5) if true == more5[i]]
a = [a for (a, c) in Iterators.zip(addresses, counts) if c > 5]
a = addresses[filter(i -> counts[i] > 5, 1:length(counts))]
a = addresses[counts .> 5]
println(a)

10 Julia Recipes You Can't Miss

Vinod V — Sat, 26 Aug 2023 11:18:54 +0000

a.Accept arbitray number of function arguments

function sum1(a, remaining...)
    sm = a + sum(remaining)
    println(sm)
end

sum1(1, 2, 3, 4,5)
sum1(1)
sum1(1,2,3)

b.You can search and replace a simple text pattern in a string using the replace() method.

stri = "I scream, you scream, we all scream for ice cream."
stri = replace(stri, "scream" => "cry")
println(stri)
stri = "I scream, you scream, we all scream for ice cream."
replace(stri, "scream"=>uppercase) #replace scream by SCREAM

c. Filtering a list based on some condition

temp = [56, 49, 50, 55, 53, 52, 60, 60, 42, 54]

# Only the elements that are greater than 50
filtered = [item for item in temp if item > 50]
print(filtered)

filtered = filter(>(50), temp)

filtered = filter(temp) do x
    x > 50
end

d.Unpack tuple to a list of variables

temp1 = (45, 51, 43, 49, 47)
a, b, c, d, e = temp1
print(a)
print(d)

stri = "Money"
c1, c2, c3, c4, c5 = stri

t = (a=1, b="hello", c = 3.14) #named tuple
(; a, c) = t; @show a; @show c; #a = 1,c = 3.14,b undefined

e.Dictionary comprehension in Julia

marks = Dict("Ram" => 60, "Shyam" => 33, "Manik" => 64,"Mani" => 76, "Hema" => 84)

more_than_55 = Dict(key => value for (key, value) in marks if value > 60)
println(more_than_55)

names = Set(["Ram", "Hema", "Mani","Manik"])
selected = Dict(key => value for (key, value) in marks if key in names)
println(selected)

f.Iterating in Reverse in Julia

for i in reverse(1:5)
  print(i)
end

r = reverse("Julia")
"ailuJ"

for i in 1:length(r)
  print(r[reverseind("Julia", i)])
end

g.Padding Strings in Julia

str = "Julia's syntax is good"
println(lpad(str, 40, '*'))
******************Julia's syntax is good
println(rpad(str, 40, '*'))
Julia's syntax is good******************

h.Return multiple values from a function

function retMultipleValues()
    v1,v2 = "India","Iran"
    return v1,v2
end

countries = retMultipleValues()
println(countries)

i. Creating and Writing to a Text File in Julia

fid = open("abc.txt", "x") #"x" flag ensures that the file is created only if it does not already exist. 
write(fid, "Hello Julia")
close(fid)

open("def.txt", "w") do fid #no need for close statement
    write(fid, "Hey")
end

#Read a text file line by line
for line in eachline("my_file.txt")
    print(line)
end

#most compact
foreach(println,eachline("my_file.txt"))

readlines("my_file.txt")

j. Parse date string to DateTime object

using Dates

dstr = "2023-07-03"
println(typeof(dstr))

date_time = DateTime(dstr, "yyyy-mm-dd")
println(date_time)
println(typeof(date_time))

New trick suggested by reviewers:
k. Make a 10-element vector of random named tuples

points = [(x=rand(), y=rand()) for i in 1:10]
distances = [sqrt(p.x^2 + p.y^2) for p in points if p.x > p.y]
#distances = [hypot(p.x , p.y) for p in points if p.x > p.y]

l. Unpacking of tuples of varying sizes (Julia version 1.9)

julia> a, b..., c = (1, 2, 3, 4)
(1, 2, 3, 4)

julia> a
1

julia> b
(2, 3)

julia> c
4

m. Convert a number array to string

julia> r = join([1,2,3.1], ", ")
"1.0, 2.0, 3.1"

n. Round a float to integer

julia> round(Int, 12.5) == Int(trunc(12.5))
true

o. To parse a binary string to integer

parse(Int, "1011", base=2)

p. To create a random binary string

 join(rand(['1','0'],10))

q. Convert a 1 x n matrix (row matrix) to a vector

B = [1 2 3 4]
1×4 Matrix{Int64}:
 1  2  3  4
julia> dropdims(B,dims=1)
4-element Vector{Int64}:
 1
 2
 3
 4
vec(B)
4-element Vector{Int64}:
 1
 2
 3
 4

r. Convert a row vector to a column vector

B = [1 2 3 4]
1×4 Matrix{Int64}:
 1  2  3  4
julia> permutedims(B)
4×1 Matrix{Int64}:
 1
 2
 3
 4

s. Check occurrence of a number in a vector

a = [1,2,3]
in(1,a) # true

t. The command splice! can be used to remove elements from an array by index

a = [1,2,3]
splice!(a,2)
a
[1,3]

u. Get all the permutations with replacement (Cartesian product)

colors = ['R','G']
n = 3
vec(Iterators.product(fill(colors, n)...) |> collect)
8-element Vector{Tuple{Char, Char, Char}}:      
 ('R', 'R', 'R')                                
 ('G', 'R', 'R')                                
 ('R', 'G', 'R')                                
 ('G', 'G', 'R')                                
 ('R', 'R', 'G')                                
 ('G', 'R', 'G')                                
 ('R', 'G', 'G')                                
 ('G', 'G', 'G')

v. Raw string in julia

x = raw"C:/temp"

w. Ternary operator with short circuited &&

x = 3 > 2 && 7 #x = 7
x = 2 > 3 && 7 #x = false

x. List all the files in a folder

readdir(raw"C:\temp")

y. Walk through the directory including subfolders

for (root,fold,files) in walkdir(raw"C:\temp") 
  for f in files                                                                            
    println(joinpath(root,f))  #print all full file paths (python os.walk)                                                                 
  end                                                                 
  end 
end

z. Pretty print a matrix

julia> using DelimitedFiles

julia> a = rand(2,3);

julia>  writedlm(stdout, a)
0.29272827910923105     0.7551746399784662      0.22809145708060086
0.9698016822675571      0.7551356902700419      0.740811475230424

aa. To get the indices of elements satisfying certain condition

 indices = findall(x -> x < 0, -4:10)

ab. Import a package with a different name

import LinearAlgebra as la

ac. Delete an element from a list

filter!(!=(1), [1,0,1,2,4])

ad. Convert data type of any to float in matrix

A = Matrix{Any}(rand(2,2))
B = similar(A,Float64)
B .= A

ae. Find the type of a Julia variable

a = 2
typeof(a)
a isa Int64
eltype(a)

af. Convert IEEE 754 hexadecimal string to float

hexStr = Meta.parse(string("0b", bitstring(1.0f0)))
0x3f800000
reinterpret(Float32,hexStr)
1.0f0
reinterpret(UInt8,[1.0f0])
4-element reinterpret(UInt8, ::Vector{Float32}):
 0x00
 0x00
 0x80
 0x3f

ag. To get the first half of the elements in a vector

x = rand(1024)
y =  x[1:end÷2] #y = x[1:512]

ah. To restart the julia shell

run(`$(Base.julia_cmd())`)

Julia Community 🟣: Vinod V

Understanding Cooley-Tukey FFT Matrix Factorization with Julia

Differences between Julia vector and Python/NumPy vector

Gram–Schmidt Orthogonalization with Julia

Julia Vectors and Matrices

What are Julia vectors?

Julia vectors

Basic data types

Basic visualization

Indexing and slicing

Copies and views

Fancy indexing

Exploring Quadratic forms with Julia

Boosting Julia Performance: Benchmarking Code in Function, Global Scope, and Builtin

Control Flow in Julia

Control Flow

The if statement

The while Statement

The for loop

The foreach statement

The break Statement

The continue Statement

Summary

Master Julia Dictionaries with These Easy Recipes

Find out what two dictionaries have in common

Extracting a subset from a dictionary

Merging dictionaries in Julia

Find the N smallest or largest items

Sort a dictionary to find minimum and maximum values

Sort a dictionary based on two keys simultaneously

Operators and Expressions in Julia

Operators

- div(9,1.81) gives 4.0

Shortcut for math operation and assignment

Evaluation Order

Changing the Order Of Evaluation

Expressions

Summary

Julia Basics

Basics

Comments

Literal Constants

Numbers

Strings

Double Quote

Strings Are Immutable

The format method

Escape Sequences

Variable

Identifier Naming

Data Types

How to write Julia programs

Example: Using Variables And Literal Constants

Logical And Physical Line

Filtering Data Made Easy: Tips and Tricks in Julia

Values greater than zero

Values less than zero

Negative values replaced with 0

Positive values replaced with 0

Filter to a range

Clamp between low and high values

Select values from a list based on values on some other list

10 Julia Recipes You Can't Miss

The `if` statement

The `for` loop

The `continue` Statement

- `div(9,1.81)` gives `4.0`