my_var = 10
println(my_var)
my_var_sqrt = sqrt(my_var)
println(my_var_sqrt)

arr = [2, 3, 5, 7, 11]
println(arr*2)

function add_one(x)
    return x + 1
end

for i in 1:10
    if add_one(i) == 6
        println("i+1 = 6!")
    else
        println("i+1 != 6, it is ", add_one(i))
    end
end

my_arr = [1, 2, 3]
println(my_arr[1])
println(my_arr[0])  # Gives error, no such index.

arr = [1, 2, 3]
println(sqrt.(arr))
println(arr .^ 2)  # Note that ^ is the power operator in Julia (** in Python).

# The following causes errors, because sqrt without dot is a scalar operator.
println(sqrt(arr))
println(arr ^ 2)

function demonstrate_types()
    x::Int = 1  # Define x as an integer
    x = x + 0.5  # Causes error, as 1.5 is not an integer.
end

demonstrate_types()

# Courtesy of Julia's documentation.

# Define an object Point, with x- and y-coordinates of some unspecified type T.
struct Point{T}
           x::T
           y::T
end

# Define the norm of point, for all points where the type T is a subtype of Real (Int, Float, etc.).
function norm(p::Point{<:Real})
    sqrt(p.x^2 + p.y^2)
end

# Firstly, we must include the linear algebra module.
using LinearAlgebra

A = [1 2 3; 0 1 4; 5 6 0]
b = [1; 2; 3]
#=
Note: We here use Julia's syntax for 
multiline comments, '#= ...  =#'.

    1 2 3
A = 0 1 4
    5 6 0
,
    1
b = 2
    3

Solve Ax=b.
=#
x = A \ b
println("x: ", x)
# Verify
println("Ax = ", A * x)

# Construct a symmetric tridiagonal matrix
diagonal = [1, 2, 3, 4]
off_diagonal = [1, 0, 1]
A = SymTridiagonal(diagonal, off_diagonal)
#=
    1 1 0 0
A = 1 2 0 0
    0 0 3 1
    0 0 1 4
=#

# eigen calls the correct method by looking at the type of the input, here LinearAlgebra.SymTridiagonal.
eigen(A)

# Create some grid
N = 1000
grid = rand(N, N)  # Populate our grid with random numbers between 0 and 1.

function energy(grid)
    energy = 0
    # This notation gives us two for-loops,
    # an outer loop over i and
    # an inner loop over j.
    for i in 2:N-1, j in 2:N-1
        right = grid[i+1, j]
        left = grid[i-1, j]
        up = grid[i, j+1]
        down = grid[i, j-1]
    
        # Nearest neighbor interaction
        nn_interaction = grid[i,j] * (right + left + up + down)
        energy += nn_interaction
    end
    return energy
end

# Measure the time it takes to execute the function
@time energy(grid)