Section 11.1 Testing Arguments
Let’s return to the factorial function that we saw in
Section 4.8:
fact(n::Integer) = n==0 ? 1 : n*fact(n-1)
We saw before that if we call this with a non-integer number like 3.2 that we will get an error, because
3.2 does not match the argument type of
Integer but what if we include a negative number? If you run
fact(-5) you will see:
StackOverflowError:
Stacktrace:
[1] fact(n::Int64) (repeats 2 times)
@ Main ~/code/sci-comp-book/Julia-output/adv-functions.ipynb:1
This occurs because when a negative number is put into the
fact function that since
\(n\) is not 0, it computes
n*fact(n-1) which will evaluate
fact again at a more negative number. Subsequent calls to
fact with more negative numbers would never end except that recursive functions are evaluated on what is called the stack and there is a maximum number of items that can be put on the stack and thus this error occurs. We can alert the user to an invalid argument with the following:
function fact(n::Integer)
n>=0 || throw(ArgumentError("The argument must be a nonnegative integer."))
n==0 ? 1 : n*fact(n-1)
end
The first line of this evaluates
n >= 0. If that is false, the statement after the
|| is evaluated and an error is thrown. Alternatively, if you enter a positive number for
n, then the first line
n>=0 evaluates as
true and then shortcircuits, which means as soon as the first part of the compound or statement hits true, it stops evaluating.
An astute reader will ask, "Why not just use an if statement?" And of course the answer is, "you may". That is the first line of the function (line 2 in the block above) can be replaced with:
if n < 0
throw(ArgumentError("The argument must be a nonnegative integer."))
end
and the same result will occur. We’ll use the style of the one-line version because it is shorter and is common in the Julia packages, however everyone’s style is their own.
If we use this newer version of
fact by calling
fact(-5), we get
ArgumentError: The argument must be a nonnegative integer.
which is a much more helpful error than the Stack Overflow error.
Check Your Understanding 11.1.
In
Chapter 9, the
findAllFactors function doesn’t make any sense if the number is less than 1. Add a line to the function that throws an appropriate error.
Section 11.2 Optional arguments
Let’s return to Newton’s method, which we wrote before as
using ForwardDiff
function newton(f::Function, x0::Number)
for _ = 1:10
dx = f(x0)/ForwardDiff.derivative(f, x0)
x0 -= dx
abs(dx) < 1e-6 && return x0
end
x0
end
Notice that we hard-coded the stopping criteria and the max number of steps. This is not good practice and will use a default value as an
optional argument.
Let’s define the tolerance (
tol) and the maximum number of steps (
max_steps) in the following way:
function newton(f::Function, x0::Number, tol::Real = 1e-6, max_steps::Integer = 10)
for _ = 1:max_steps
dx = f(x0)/ForwardDiff.derivative(f, x0)
x0 -= dx
abs(dx) < tol && return x0
end
x0
end
Note that when entering this, it says there are 3 methods with the name
newton. This is because Julia will build three different function signatures. One with 2 arguments (and both
tol and
max_steps use the default values), one with 3 arguments (and
max_steps using it’s default value) and one with all 4 arguments. Also remember that
1e-6 means
\(10^{-6}\text{,}\) which is the default tolerance.
This seems more robust in that we can now call Newton’s method with different values of tolerance and steps. So:
returns
2.236067977915804. But if we use a lower tolerance:
newton(x -> x^2-5,2,1e-3)
returns
2.236111111111111, which is slightly different than before. Probably the number of steps of Newton was a bit different.
If we want to change the number of steps, however, we need to include the tolerance as well like:
newton(x -> x^2-5,2,1e-6,5)
which results in
2.236111111111111. We will see later an alternative way to handle these parameters, called
keyword parameters which requires names of parameters, not the order, which is helpful as the number of arguments/parameters for a function gets large.
This is quite an improvement, however, the arguments
tol and
max_steps optional parameters should both be positive, so we will add checks on these as
function newton(f::Function, x0::Number, tol::Real = 1e-6, max_steps::Integer = 10)
tol > 0 || throw(ArgumentError("The parameter tol must be positive"))
max_steps > 0 || throw(ArgumentError("The parameter max_steps must be positive"))
for _ = 1:max_steps
dx = f(x0)/ForwardDiff.derivative(f, x0)
x0 -= dx
abs(dx) < tol && return x0
end
x0
end
and including negative numbers for either of these two now throws an error.
Check Your Understanding 11.2.
(a)
Recall that the positive root of
\(f(x)=x^2-10\) is
\(\sqrt{10}\text{.}\) Use the function above to find
\(\sqrt{10}\) with a tolerance of
\(10^{-10}\text{.}\)
(b)
Find the absolute and relative error. Use
sqrt(10) as the actual value and the result from the
newton function as the Hint: recall this from
Chapter 10.
(c)
Put in arguments for the tolerance and maximum number of steps that will throw an error. Check that this actually occurs.
Section 11.3 Handling Special Cases
root=newton(x -> x^2+1,2)
returns
2.4008803928468465, but if we evaluate the function at the root, with
f(root) the result is
6.764226660756428, which is definitely not 0 (or very close), so this doesn’t appear to be a root. Recall this occurred because the function
\(f(x)=x^2+1\) doesn’t have a root. If we temporarily print out the values of
x1 within the loop, we’ll see that the x values bounce all around and then just stops. It’s not clear, but what happens here is that the max number of steps is reached, but you are not alerted.
Note that if it never reaches the case the
dx < tol, then it doesn’t return from within the loop. When this case happens, the for loop ends, therefore if we replace line 9 with
throw(ErrorException("The maximum number of steps: $max_steps was reached without convergence"))
and then rerunning
newton(x -> x^2+1,2) gives the error
The maximum number of steps: 10 was reached without convergence which explains to the user that a root was not reached.
You, the astute reader, probably noticed that we had an alternative way to handle this case in
Chapter 12, in which we created a struct called
Root to handle all of the information from Newton’s method. Either using the
Root structure of the method of throwing an exception like above is fine as you may need to know how you will use this code.
Section 11.5 Parametric Methods
Let’s look at writing a function that finds the maximum of some number of real values. There is a built-in function for this, but going through this will illustrate a point.
If we want to find the maximum of two reals, then
findMax(x::Real,y::Real) = x > y ? x : y
and if we want any number of reals then as we saw in
Section 4.6, we can use variable arguments (varargs) like:
function findMax(nums::Real...)
local max = -Inf
for num in nums
if num > max
max = num
end
end
max
end
Another
findMax function that would be helpful is in an array of
Reals (mainly because we will need to do something different for
Complex arrays). This will look a lot like the one above:
function findMax(arr::Vector{Real})
local max = -Inf
for num in arr
if num > max
max = num
end
end
max
end
Testing it with the array
x=collect(1:10) using
findMax(x) we get the error:
MethodError: no method matching findMax(::Vector{Int64})
The function `findMax` exists, but no method is defined for this combination of argument types.
and this is because even though
x is of type
Vector{Int64} that doesn’t fit the signature
Vector{Real}.
One way to get around this is to write a function for each number type we may want to find the max of (
Int128,
Int64, ...,
BigInt,
Float64,
Float32,
Float16) and all of the flavors of
Rational but if you notice every function will be nearly identical and this would be painful and a lot of code copying. Also, this is a nightmare for maintenance, in that if there was a bug in one it would need to be tracked down in each one.
Instead, we will use something in Julia called a
parametric function which allows us to write a single function that works on a set of different types. To do this with the
findMax function, we will write:
function findMax(arr::Vector{<:Real})
local max = -Inf
for num in arr
if num > max
max = num
end
end
max
end
and this will now match a vector of any type that is a subtype of
Real (including integers, floats and rational) and then
findMax(x) returns
10. Also if we have
findMax(collect(1.0:10.0))
findMax(collect(big(1):big(10)))
returns
10 (as a
BigInt).
The function signature for the
findMax function is a shorthand notation for
function findMax(arr::Vector{T}) where T <: Real
which will match
Vector{T} for any subtypeof
T being
Real. In the case above, we didn’t need to know the type of vector, however consider a function that calculates the product of all of the elements of a real vector. The following does this:
function findProduct(arr::Vector{T}) where T <: Number
local prod = one(T)
for x in arr
prod *= x
end
prod
end
and the type is needed to call the
one function which return the number
1 for the type
T. (Try calling
one(Float64), one(Int), one(Rational{Int64}) to determine that this does what is intended.)
Check Your Understanding 11.4.
Create a
median method on an array of reals. It should be parametric with either of the following two sigatures:
function median(arr::Vector{T}) where T <: Real
end
function median(arr::Vector{<Real})
end
depending on if you need to know the type
T or not. Recall that the median of a set of sorted numbers is either the middle number if the length is odd or the mean of the middle two numbers if the length is even.
Try the function with each of the following
(a)
(b)
arr2=collect(big(1):big(10))
(c)
arr3=collect(3.0:4.0:56.0)
(d)
arr4=collect(1:1//3:7//3)
Section 11.6 Variable Scope
If you are using the command line Julia (also called the REPL) or in a notebook, then entering
x=6 will create
x as a global variable, even without saying so explicitly. If we consider the factorial function above:
function fact(n::Integer)
prod=1
for i=1:n
prod *= i
end
prod
end
then recall that
n is the function argument and the variable
prod is actually declared locally implicitly and it is good form to say it is local by
function fact(n::Integer)
local prod=1
for i=1:n
prod *= i
end
prod
end
The variable
prod is not defined or visible outside the function. If we type
2*prod outside the function for example, an error will occur saying that
prod is not defined. It is actually a good thing that variables inside a function is only visible inside the function (that is has local scope), because this encapsulates the code inside a function in that the function doesn’t affect anything outside of it.
Let’s look at local a global variables in a bit more detail, by the following function whose only purpose is to understand global and local variables.
x=3
function f()
local x=2
x
end
f()
then evaluating the function with
f() this will return
2 and typing
x outside the function just to see it’s (to see it’s value) returns
3, indicating that there are two different
x variables.
If we indeed need a variable inside of a function to use a global variable, we can use the
global keyword. If instead we have
x=3
function g()
global x
x += 1
end
then executing
g() results in
4 showing that it is using the global value. Entering
x after evaluating the function, also will return
4. Although it is rare to need to use a global variable, we will use this technique to determine the number of function evaluations to test a function in
Section 8.8.
Section 11.7 Function Arguments
As described in the
Julia documentation on functions, arguments follow what is called "pass-by-sharing", which means that values are not copied when they are passed to functions. Instead, the argument is a local variable that is passed the value when the function is called. Consider
function g(x::Number)
2*x
end
and if we make a function call like
g(4), then within the function
x is a variable that is given the value
4. We can update the variable. For example.
function g(x::Number)
x *= 2
x += 7
x
end
where recall that the
*= and
+= updates the value on the left. When called as
it returns
15 and the value of
x remains as
4. This is because the
x within the function is now a local function that starts with the value 4.
You may recall some functions seem to update the value of its arguments. For example,
A = [4,5,2,3,1]
sort!(A)
A
returns the sorted array
[1,2,3,4,5]. So how does this work, you may ask? Let’s write our own function that may do this.
function change3!(A::Vector{T}) where T <: Number
A[3] *= 2
end
which takes a vector of numbers and double the 3rd element of the vector. We used the convention in Julia that if the arguments are changed we append a
! to the name. If we enter
A = [1,2,3,4]
change3!(A)
which returns 6 (recall that since the last line of the function is the double of the 3rd element, we get 6). And if we type
A to view it again, we get the vector
[1,2,6,4] indicating that indeed the vector
A has been updated.
The reason that
A is updated with this function whereas the
x above is not is the way that array variables differ from scalar variables (like numbers and strings). The assignment
A=[1,2,3,4], means that
A is a reference to the array instead of the array itself. Thus, when you call
change3!(A) the variable
A within the function takes on the reference to the original function. Thus, when values are updated within the function, then the array is actually updated. When
x is updated in the function above, one can think that this is a copy of the passed in value, which doesn’t leave the function.
