Table des matières
The objective of this lab is to give an introduction to visualization with Matplotlib and scientific computing with Numpy. We'll implement a simple function with different programming methods, measure the performances of each and visualize the results.
Classic programming methods
Recursion
Develop a function fibonacci_rec(n) that computes the value of the n-th term of the Fibonacci sequence. This function must be recursive. The expected results for the first values are: {$ 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181... $}
ANSWER ELEMENTS
def fibonacci_rec(n):
if n == 0:
return 0
if n == 1:
return 1
else:
return fibonacci_rec(n-1)+fibonacci_rec(n-2)
Time analysis
The timeit function of the package bearing the same name is used to measure the execution time of Python code. In the following example, timeit returns the total time of 100 executions of the Python code to smooth the results.
Execution times may vary from machine to machine; you should expect an execution time of a few thousandths of a second. Try it.
from timeit import timeit timeit(lambda : fibonacci_rec(10), number=100)
We defined a function measure() that determines the execution time of a function for a sequence of values.
When we measure the performances of fibonacci() with different values of {$n$}, we should observe a high execution time for values
beyond {$n = 20$}.
def measure(func, values, number=100):
m = []
for i in values :
m.append(timeit(lambda : func(i), number=number))
return m
measure(fibonacci_rec, [1, 2, 3, 5, 10, 20, 25])
Iteration
Develop a function fibonacci_iter(n) that computes the value of the n-th term of the Fibonacci sequence. This function must be iterative.
Remark for advanced students : do not use memoization.
Measure the execution time of fibonacci_iter(n). It should be significantly faster than the recursive version.
ANSWER ELEMENTS
def fibonacci_iter(n):
if n == 0:
return 0
else:
f = [0, 1]
for i in range(2, n + 1):
f.append(f[-1] + f[-2])
return f[n]
Visualization
We'll now plot these values. We use the pyplot module from the matplotlib library.
The following instructions plot a series of points defined by two lists of coordinates x and y.
import matplotlib.pyplot as plt plt.plot(xvalues, yvalues, label="name for legend") plt.legend() plt.show()
The graph should always contain a title and a legend, have named axes, the scale of the y-axis must be logarithmic, the scale of the x-axis must start at 0. Here are some useful functions.
plt.yscale("log")
plt.xlabel("Name")
plt.ylabel("Name")
plt.title("Name")
Be inclusive, don't forget color-blind people (about 1 in 20 people) and choose line styles that can be told apart without using colors. Line styles can be set with the following codes: ":", "-", "-." and "--".
plt.plot(xvalues, yvalues, "--", label="Dashed line")
Put it all together
Define a function main() that measures that execution time of fibonacci_rec() and fibonacci_iter(), and
plots the performance curves.
ANSWER ELEMENTS
import matplotlib.pyplot as plt
def main():
l = list(range(1, 20, 1))
frec = measure(fibonacci_rec, l)
fiter = measure(fibonacci_iter, l)
plt.yscale("log")
plt.title("Execution time of fibonacci functions")
plt.xlabel("Fibonacci number")
plt.ylabel("Time in second")
plt.plot(l, frec, "-.", label="fibonacci_rec")
plt.plot(l, fiter, "--", label="fibonacci_iter")
plt.legend()
plt.show()
main()
Alternative method
An alternative method to compute the Fibonacci sequence consists in using matrices. A matricial form of the Fibonacci sequence is given below:
{$$ \left( \begin{matrix} F_1 \\ F_2 \end{matrix} \right) = \left( \begin{matrix} 0 & 1 \\ 1 & 1 \end{matrix} \right) * \left( \begin{matrix} F_0 \\ F_1 \end{matrix} \right) $$}
and in general:
{$$ \left( \begin{matrix} F_n \\ F_{n+1} \end{matrix} \right) = \left( \begin{matrix} 0 & 1 \\ 1 & 1 \end{matrix} \right)^n * \left( \begin{matrix} F_0 \\ F_1 \end{matrix} \right) $$}
We'll use numpy to compute operations on matrices.
Here are some simple steps to help you.
Creation of matrices and vectors.
import numpy as np mat1 = np.array([[5,6],[10,3]]) vec = np.array([[9],[2]])
Access to one element.
mat1[1,0] # => 10
Some term-to-term operations.
mat2 = mat1 + mat1 # => mat2[i,j] = mat1[i, j] + mat1[i,j] mat3 = mat1 * 2 # => mat3[i,j] = mat1[i,j] * 2 ...
Some matrix operations.
mat4 = np.matmul(mat1, mat1) # matrix product mat5 = mat1 @ mat1 # abbreviated product syntax mat6 = np.linalg.matrix_power(mat1, 2) # mat1**2
Take a moment to familiarize yourself with numpy matrix operations in a Python shell.
WARNING: you may already have used np.matrix instead of np.array. They work slightly differently. In particular, operations like M ** 2 or M * vec are available for np.matrix; it all seems very practical. HOWEVER, the numpy documentation is very clear on this point, the np.matrix object will be removed from the library soon, so you have to stop using them: "It is no longer recommended to use this class, even for linear algebra. Instead use regular arrays. The class may be removed in the future".
Matrix
Develop a function fibonacci_mat(n) that computes the value of the n-th term of the Fibonacci sequence. This function must use numpy functions.
Include this function in your performance analysis.
ANSWER ELEMENTS
import numpy as np
def fibonacci_mat(n):
matrix = np.array([[0,1],[1,1]])
vec = np.array([[0],[1]])
return (np.linalg.matrix_power(matrix, n-1) @ vec)[1, 0]
def main():
l = list(range(1, 20, 1))
frec = measure(fibonacci_rec, l)
fiter = measure(fibonacci_iter, l)
fmat = measure(fibonacci_mat, l)
plt.yscale("log")
plt.title("Execution time of fibonacci functions")
plt.xlabel("Fibonacci number")
plt.ylabel("Time in second")
plt.plot(l, frec, "-.", label="fibonacci_rec")
plt.plot(l, fiter, "--", label="fibonacci_iter")
plt.plot(l, fmat, "-", label="fibonacci_mat")
plt.legend()
plt.show()
main()
Large values of {$n$}
Your recursive function is probably much slower than the other two variants. The execution time seems to increase exponentially with the value of {$n$}. And it does. Test your functions with larger values of {$n$}. Start with range(1, 20, 2), then range(1, 30, 3) and so on until the execution time of fibonacci_rec becomes too long.
Test the three functions.
ANSWER ELEMENTS
def main():
l = list(range(0, 20, 1))
l2 = list(range(0, 20, 1)) + list(range(20, 1000, 10))
frec = measure(fibonacci_rec, l)
fiter = measure(fibonacci_iter, l2)
fmat = measure(fibonacci_mat, l2)
plt.yscale("log")
plt.title("Execution time of fibonacci functions")
plt.xlabel("Fibbnacci number")
plt.ylabel("Time in second")
plt.plot(l, frec, "-.", label="fibonacci_rec")
plt.plot(l2, fiter, "--", label="fibonacci_iter")
plt.plot(l2, fmat, "-", label="fibonacci_mat")
plt.legend()
plt.show()
main()
The devil is in the detail
Compare the results of fibonacci_iter(93) and fibonacci_mat(93). Where does the problem come from?
ANSWER ELEMENTS
In Python integers haven't a fixed size, the only limit is the available memory. The interpreter manages the overflows and, if necessary, a greater number of bits is used.
Integer values in Numpy have a fixed size. By default an integer is a signed 64 bits value. Overflows are not detected. Conclusion: Numpy is a precious library but you need to be careful when you use it.
Well done
You have completed the first part of this lab. You studied the basics of visualization with matplotlib and, for the first time, touched with hand the simplicity of matrix calculations with numpy.
If you still have time, you may want to check the next following questions.
Optional questions
The problem with the recursive function comes from the exponentially large number of recursive calls. The function may be called many times with the same input value. It is possible to work around this problem through memoization. The rationale of memoization is to store the intermediate results in a dictionary and make a recursive call only if the result of the call is not in the dictionary yet.
Memoization
Modify the implementation of your recursive function to use memoization and compare its performances with the two other functions.
You can find out more on memoization here.
ANSWER ELEMENTS
Warning: when the default value of a parameter is a mutable object (dictionary, list, etc.), this object is preserved between calls. Hence, in the following example, the memoization dictionary is completed only once. This biases the comparisons, as the timeit module averages the execution time over several calls.
def fibonacci_rec(n, memo={}):
if n == 0:
return 0
if n == 1:
return 1
else:
if n in memo:
return memo[n]
else :
f = fibonacci_rec(n-1, memo=memo) + fibonacci_rec(n-2, memo=memo)
memo[n] = f
return f
Two solutions are possible. The first is to use an auxiliary function to explicitly initialize the dictionary. The second is to set the default value to None.
def fibonacci_rec(n, memo=None):
if n == 0:
return 0
if n == 1:
return 1
else:
if memo is None:
memo = {}
if n in memo:
return memo[n]
else :
f = fibonacci_rec(n-1, memo=memo) + fibonacci_rec(n-2, memo=memo)
memo[n] = f
return f
Very large values of {$n$}
Up to what value of {$n$} can you compute results? It depends on your machine. Increase the measurement range, ex. from {$0$} to {$20,000$} with a step of {$1,000$}. One of the functions should generate an error message for too large values of {$n$}.
In your opinion
- Why does
fibonacci_mat(n)have an execution time which increases for small values of {$n$} and then seems to be constant? - What does
fibonacci_rec(n)error message means for large values of {$n$}?
ANSWER ELEMENTS
- To compute the power of a matrix {$M^n$},
numpyuses a "exponentiation by squaring" algorithm. The number of operations performed is proportional to {$\log_2(n)$}, to observe this increase we would need values of {$n$} much more important. The notion associated with this question is the complexity of algorithms, we will come back to this during the ST2 algorithm course. - Recursive calls are accumulated on the stack of the processor. The stack has a limited memory size. When the stack size is exceeded the processor can no longer continue its normal execution, an exception is raised. You'll take a closer look at this notion if you'll follow the third year major "Computer systems architecture".
Try the matrix computation your way
Since matrix computation with Numpy generates overflows, you can implement your own method using lists of native Python integer numbers. You'll need a function to compute the product of matrices (matrices will be represented as lists of lists), and a fast exponentiation function.
ANSWER ELEMENTS
def prod(m1, m2):
result = [[sum(a * b for a, b in zip(X_row, Y_col)) for Y_col in zip(*m2)] for X_row in m1]
return result
def fast_expo(m, e):
result = [[1, 0], [0, 1]]
while e > 0:
if e % 2 == 1:
result = prod(result, m)
m = prod(m, m)
e = e // 2
return result
def fibonacci_mat(n):
matrix = [[0, 1], [1, 1]]
vec = [[0], [1]]
e = fast_expo(matrix, n-1)
return prod(e, vec)[1][0]