8) Functions

Last time:

C pointers and memory pools

Today:

Using Julia
Introduction to floating point arithmetic
2.1. Plotting
Machine Epsilon

1. Julia¶

To recap from first class, Julia is a relatively new programming language. Think of it as MATLAB done right, open source, and fast. It’s nominally general-purpose, but mostly for numerical/scientific/statistical computing. There are great learning resources. We’ll introduce concepts and language features as we go.

# The last line of a cell is output by default
x = 3
y = 4

4

println("$x + $y = $(x + y)") # string formatting/interpolation
4;  # trailing semicolon suppresses output

3 + 4 = 7

@show x + y
x * y

x + y = 7

12

1.1. Numbers¶

3, 3.0, 3.0f0, big(3.0) # integers, double precision, single precision, and convert to a maximum precision representation

(3, 3.0, 3.0f0, 3.0)

typeof(3), typeof(3.0), typeof(3.0f0), typeof(big(3.0))

(Int64, Float64, Float32, BigFloat)

# automatic promotion
@show 3 + 3.0
@show 3.0 + 3.0f0
@show 3 + 3.0f0;

3 + 3.0 = 6.0
3.0 + 3.0f0 =

6.0
3 + 3.0f0 = 6.0f0

# floating and integer division
@show 4 / 2
@show -3 ÷ 2; # type `\div` and press TAB

4 / 2 = 2.0
-3 ÷ 2 = -1

1.2 Arrays¶

[1, 2, 3]

3-element Vector{Int64}:
 1
 2
 3

# explicit typing
Float64[1,2,3]

3-element Vector{Float64}:
 1.0
 2.0
 3.0

# promotion rules similar to arithmetic
[1,2,3.] + [1,2,3]

3-element Vector{Float64}:
 2.0
 4.0
 6.0

x = [10., 20, 30]
x[2] # one-based indexing

20.0

x[2] = 3.5

3.5

# multi-dimensional array
A = [10 20 30; 40 50 60]

2×3 Matrix{Int64}:
 10  20  30
 40  50  60

Compare with the Python notation:

A = np.array([[10, 20, 30], [40, 50, 60]])

1.3 Functions¶

In Julia, a function is an object that maps a tuple of argument values to a return value.
- A tuple is a built-in, fixed-length data structure that can hold any values, but cannot be modified (it is immutable). Components of tuples can optionally be named, in which case a named tuple is constructed.
The value returned by a function is the value of the last expression evaluated, which, by default, is the last expression in the body of the function definition. You don’t need to have a return statement necessarily, but it helps readability and in cases in which the control flow of the program is not linear (e.g., if you have for, while, or branching of conditional cases).


# define a function called 'f' that only takes two arguments
function f(x, y)
    sqrt(x*x + y*y)
end

f (generic function with 1 method)

What is the difference between a method and a function in Julia? Check this resource.
A definition of one possible behavior for a function is called a method. A function can then have multiple methods (or definitions).

# define a function called 'f' that takes three arguments, one of which is optional and has a default value
function f(x, y; z=0)
    sqrt(x*x + y*y) * z
end

# define another function also called 'f' that only takes two arguments
function f(x, y)
    sqrt(x*x + y*y)
end

f (generic function with 1 method)

# if I invoke f this way, which one will the Julia compiler use?
f(3, 4, z=5)

25.0

# if instead I invoke f this way, which one will the Julia compiler use?
f(2, 3)

3.605551275463989

In other words, optional arguments are tied to a function, not to any specific method of that function. It depends on the types of the optional arguments which method is invoked.
Keyword arguments behave quite differently from ordinary positional arguments. In particular, they do not participate in method dispatch. Methods are dispatched based only on positional arguments, with keyword arguments processed after the matching method is identified.

We have just seen an example of one of the most powerful features of Julia: dynamic (multiple) dispatch:

Functions can have multiple definitions as long each definition restricts the type of the parameters differently. It is the type of the parameters that define which “definition” (or “method” in Julia terminology) will be called.
This is the way Julia mimics polymorphism that some compiled languages (e.g., C++) have.

There are also compact “assignment form” function definitions:

g(x, y) = sqrt(x^2 + y^2)
g(3, 4)

5.0

In the assignment form above, the body of the function must be a single expression, although it can be a compound expression.

Anonymous functions are usually functions so short-lived that they do not need a name. This is done with an arrow notation. Example:

((x, y) -> sqrt(x^2 + y^2))(3, 4)

5.0

Argument Passing Behavior¶

Julia function arguments follow a convention sometimes called “pass-by-sharing”, which means that values are not copied when they are passed to functions. Function arguments themselves act as new variable bindings (new “names” that can refer to values), much like assignments argument_name = argument_value, so that the objects they refer to are identical to the passed values. Modifications to mutable values (such as Arrays) made within a function will be visible to the caller. (This is the same behavior found in Scheme, most Lisps, Python, Ruby and Perl, among other dynamic languages.)

As a common convention in Julia (not a syntactic requirement), a function that modifies any of its arguments would typically be named with a !, e.g., f!(x, y) rather than f(x, y), as a visual reminder at the call site that at least one of the arguments (often the first one) is being mutated.

To read more about Julia functions, check the Julia documentation.

Argument-type declarations¶

You can declare the types of function arguments by appending ::TypeName to the argument name. Example:

function f(x::Float64, y::Float64)
    sqrt(x*x + y*y)
end

f (generic function with 2 methods)

Note: Argument-type declarations normally have no impact on performance: regardless of what argument types (if any) are declared, Julia compiles a specialized version of the function for the actual argument types passed by the caller.

The most common reasons to declare argument types in Julia are, instead: dispatch, correctness, and clarity.

1.4 Loops¶

# ranges
1:50000000

1:50000000

collect(1:5)

5-element Vector{Int64}:
 1
 2
 3
 4
 5

x = 0
for n in 1:50000000
    x += 1/n^2
end
@show x
x - pi^2/6 # Basel problem (solved by Euler; generalized by Riemann's zeta function)

x = 1.6449340467988642

-2.0049362170482254e-8

# list comprehensions
sum([1/n^2 for n in 1:1000])

1.643934566681561

2. Floating point arithmetic¶

0.1 + 0.2
(1 + 1.2e-16) - 1

2.220446049250313e-16

using Pkg
Pkg.add("Plots")

using Plots
default(linewidth=4)
plot(x -> (1 + x) - 1, xlims=(-1e-15, 1e-15),
xlabel="x", ylabel="y")
plot!(x -> x)

     Project No packages added to or removed from `~/.julia/environments/v1.13/Project.toml`
    Manifest No packages added to or removed from `~/.julia/environments/v1.13/Manifest.toml`

3. Machine epsilon¶

We approximate real numbers with floating point arithmetic, which can only represent discrete values. In particular, there exists a largest number, which we call $\epsilon_{\text{machine}}$ , such that

1 \oplus x = 1 \quad \text{for all}\ \lvert x \rvert < \epsilon_{\text{machine}}.

(1)

The notation $\oplus, \ominus, \odot, \oslash$ represent the elementary operation carried out in floating point arithmetic.

eps = 1
while 1 + eps != 1
    eps = eps / 2
end
eps

1.1102230246251565e-16

eps = 1.f0
while 1 + eps != 1
    eps = eps / 2
end
eps

5.9604645f-8

3.1 Approximating `exp`¶

e^x = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \dotsb

(2)

Suppose we want to compute $f(x) = e^x - 1$ for small values of $x$ .

f1(x) = exp(x) - 1
y1 = f1(1e-8)

9.99999993922529e-9

f2(x) = x + x^2/2 + x^3/6
y2 = f2(1e-8)

1.000000005e-8

Which answer is more accurate?

@show (y1 - y2)        # Absolute difference
@show (y1 - y2) / y2;  # Relative difference

y1 - y2 = -1.1077470910720506e-16
(y1 - y2) / y2 = -1.1077470855333152e-8

0 / 0

NaN

1. Julia¶

1.1. Numbers¶

1.2 Arrays¶

1.3 Functions¶

Argument Passing Behavior¶

Argument-type declarations¶

1.4 Loops¶

2. Floating point arithmetic¶

3. Machine epsilon¶

3.1 Approximating exp¶

3.1 Approximating `exp`¶