Homework 01

As always there are multiple options for the circlemat function definition. Here we show the two, that we have foreshadowed in the homework hints:

function circlemat(n)
    A = zeros(Int, n, n) # creates nxn matrix of zeros
    for i in 1:n
        for j in 1:n
            if (i == j-1 && j > 1) || (i == n && j == 1) || (i == j+1 && j < n) || (i == 1 && j == n)
                A[i,j] = 1
            end
        end
    end
    return A
end
circlemat(10)

circlemat(n) = [(i == j-1 && j > 1) || (i == n && j == 1) || (i == j+1 && j < n) || (i == 1 && j == n) ? 1 : 0 for i in 1:n, j in 1:n]
circlemat(10)
Both version give the same answer.

Extending the original polynomial to matrix valued point x = A, requires only small changes to the initialization of accumulator variable. Running directly the original code fails on MethodError, because julia cannot add a matrix to an integer.

function polynomial(a, x::AbstractMatrix) # we are limiting this function to everything that is subtype of AbstractMatrix
    accumulator = zeros(eltype(x), size(x)) # zeros of the same type and size as `x`
    for i in length(a):-1:1
        accumulator += x^(i-1) * a[i] # ! 1-based indexing for arrays
    end
    return accumulator
end

A = circlemat(10) # matrix of size 10x10
coeffs = ones(4)  # coefficients of polynomial
polynomial(coeffs, A)

The other option is to use the more abstract version that we have defined to work with generators/iterators. For example

polynomial(a, x) = sum(ia -> x^(ia[1]-1) * ia[2], enumerate(a))

works out of the box

polynomial(coeffs, A)

General notes

# no argument types specified 
# will work as a fallback when no more concrete example is found
function polynomial(a, x)
    accumulator = 0
    for i in length(a):-1:1
        accumulator += x^(i-1) * a[i]
    end
    return accumulator
end
 
# with this we are defining a special method for situations, where x is a matrix
# multiple dispatch always tries to find most specified/concrete implementation for a given argument type
function polynomial(a, x::AbstractMatrix) # equivalent to function polynomial(a, x::T) where T<:AbstractMatrix
    accumulator = zeros(size(x))
    for i in length(a):-1:1
        accumulator += x^(i-1) * a[i]
    end
    return accumulator
end

As you can see there is quite a bit of code duplication, thus in this case it is more suitable to design this differently, by taking the generic untyped function and introduce a initialization function in the place of accumulator.

function polynomial(a, x)
    accumulator = _init_zero(x)
    for i in length(a):-1:1
        accumulator += x^(i-1) * a[i]
    end
    return accumulator
end
_init_zero(x) = 0; 
_init_zero(x::AbstractMatrix) = zeros(size(x))

This is still not an ideal solution, because it does not treat the eltype of the matrix or the typeof the scalar input x. We will show in later lectures, what problems this “better” design brings. If you are interested this functionality is already built in Julia with zero and one functions.

In one of the solutions, we have seen a type of dynamic dispatch based on an if condition

function polynomial(a, x)
    if typeof(x) <: AbstractMatrix
        accumulator = zeros(size(x))
    else
        accumulator = 0
    end
 
    for i in length(a):-1:1
        accumulator += x^(i-1) * a[i]
    end
    return accumulator
end
which essentially mimics the design with the _init_zero function above. In this simple case Julia will compile away the if condition, because for each of the argument type we get a different method and different if branch as well being compiled, however this is not a good design pattern, because new type of initialization requires changes in the original method. Note that this type of dynamic dispatch is often encountered in Python, where due to the interpretation nature of the language, the condition is always evaluated.