@codeprogrammer Scientific Python Lectures
@codeprogrammer Scientific Python Lectures
2024
SciPy Python EDITION
IP[y]:
Cython IPython
Edited by
Lectures
Pierre de Buyl
K. Jarrod Millman
Stéfan van der Walt
lectures.scientific-python.org
i
5.3 Linear algebra operations: scipy.linalg . . . . . . . . . . . . . . . . . . . . . . . . . . . 214
5.4 Interpolation: scipy.interpolate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 215
5.5 Optimization and fit: scipy.optimize . . . . . . . . . . . . . . . . . . . . . . . . . . . . 217
5.6 Statistics and random numbers: scipy.stats . . . . . . . . . . . . . . . . . . . . . . . . 222
5.7 Numerical integration: scipy.integrate . . . . . . . . . . . . . . . . . . . . . . . . . . . 224
5.8 Fast Fourier transforms: scipy.fft . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 226
5.9 Signal processing: scipy.signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 229
5.10 Image manipulation: scipy.ndimage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 231
5.11 Summary exercises on scientific computing . . . . . . . . . . . . . . . . . . . . . . . . . . 236
5.12 Full code examples for the SciPy chapter . . . . . . . . . . . . . . . . . . . . . . . . . . . 249
ii
13.3 Full code examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 446
13.4 Examples for the mathematical optimization chapter . . . . . . . . . . . . . . . . . . . . 446
13.5 Practical guide to optimization with SciPy . . . . . . . . . . . . . . . . . . . . . . . . . . 486
13.6 Special case: non-linear least-squares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 488
13.7 Optimization with constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 490
13.8 Full code examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 491
13.9 Examples for the mathematical optimization chapter . . . . . . . . . . . . . . . . . . . . 491
iii
IV About the Scientific Python Lectures 690
19 About the Scientific Python Lectures 691
19.1 Authors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 691
Index 694
iv
Scientific Python Lectures, Edition 2024.2rc0.dev0
Contents 1
Part I
2
Scientific Python Lectures, Edition 2024.2rc0.dev0
This part of the Scientific Python Lectures is a self-contained introduction to everything that is needed
to use Python for science, from the language itself, to numerical computing or plotting.
3
CHAPTER 1
Python scientific computing ecosystem
• Batteries included Rich collection of already existing bricks of classic numerical methods, plot-
ting or data processing tools. We don’t want to re-program the plotting of a curve, a Fourier
transform or a fitting algorithm. Don’t reinvent the wheel!
• Easy to learn Most scientists are not paid as programmers, neither have they been trained so.
They need to be able to draw a curve, smooth a signal, do a Fourier transform in a few minutes.
• Easy communication To keep code alive within a lab or a company it should be as readable
as a book by collaborators, students, or maybe customers. Python syntax is simple, avoiding
strange symbols or lengthy routine specifications that would divert the reader from mathematical
or scientific understanding of the code.
4
Scientific Python Lectures, Edition 2024.2rc0.dev0
• Efficient code Python numerical modules are computationally efficient. But needless to say that
a very fast code becomes useless if too much time is spent writing it. Python aims for quick
development times and quick execution times.
• Universal Python is a language used for many different problems. Learning Python avoids learning
a new software for each new problem.
Pros
• Very fast. For heavy computations, it’s difficult to outperform these languages.
Cons
• Painful usage: no interactivity during development, mandatory compilation steps,
verbose syntax, manual memory management. These are difficult languages for
non programmers.
Pros
• Very rich collection of libraries with numerous algorithms, for many different do-
mains. Fast execution because these libraries are often written in a compiled lan-
guage.
• Pleasant development environment: comprehensive and help, integrated editor, etc.
• Commercial support is available.
Cons
• Base language is quite poor and can become restrictive for advanced users.
• Not free and not everything is open sourced.
Julia
Pros
• Fast code, yet interactive and simple.
• Easily connects to Python or C.
Cons
• Ecosystem limited to numerical computing.
• Still young.
Pros
• Open-source, free, or at least cheaper than Matlab.
• Some features can be very advanced (statistics in R, etc.)
Cons
• Fewer available algorithms than in Matlab, and the language is not more advanced.
• Some software are dedicated to one domain. Ex: Gnuplot to draw curves. These
programs are very powerful, but they are restricted to a single type of usage, such
as plotting.
Python
Pros
• Very rich scientific computing libraries
• Well thought out language, allowing to write very readable and well structured code:
we “code what we think”.
• Many libraries beyond scientific computing (web server, serial port access, etc.)
• Free and open-source software, widely spread, with a vibrant community.
• A variety of powerful environments to work in, such as IPython, Spyder, Jupyter
notebooks, Pycharm, Visual Studio Code
Cons
• Not all the algorithms that can be found in more specialized software or toolboxes.
Unlike Matlab, or R, Python does not come with a pre-bundled set of modules for scientific computing.
Below are the basic building blocks that can be combined to obtain a scientific computing environment:
• NumPy: numerical computing with powerful numerical arrays objects, and routines to manip-
ulate them. https://numpy.org/
See also:
chapter on numpy
• SciPy : high-level numerical routines. Optimization, regression, interpolation, etc https://scipy.
org/
See also:
chapter on SciPy
• Matplotlib : 2-D visualization, “publication-ready” plots https://matplotlib.org/
See also:
chapter on matplotlib
Advanced interactive environments:
• IPython, an advanced Python console https://ipython.org/
• Jupyter, notebooks in the browser https://jupyter.org/
Domain-specific packages,
• pandas, statsmodels, seaborn for statistics
• sympy for symbolic computing
• scikit-image for image processing
• scikit-learn for machine learning
and many more packages not documented in the Scientific Python Lectures.
See also:
chapters on advanced topics
chapters on packages and applications
Python comes in many flavors, and there are many ways to install it. However, we recommend to install
a scientific-computing distribution, that comes readily with optimized versions of scientific modules.
Under Linux
If you have a recent distribution, most of the tools are probably packaged, and it is recommended to use
your package manager.
Other systems
There are several fully-featured scientific Python distributions:
• Anaconda
• WinPython
Interactive work to test and understand algorithms: In this section, we describe a workflow
combining interactive work and consolidation.
Python is a general-purpose language. As such, there is not one blessed environment to work in, and
not only one way of using it. Although this makes it harder for beginners to find their way, it makes it
possible for Python to be used for programs, in web servers, or embedded devices.
We recommend an interactive work with the IPython console, or its offspring, the Jupyter notebook.
They are handy to explore and understand algorithms.
Start ipython:
In [2]: print?
Signature: print(*args, sep=' ', end='\n', file=None, flush=False)
Docstring:
Prints the values to a stream, or to sys.stdout by default.
sep
string inserted between values, default a space.
end
string appended after the last value, default a newline.
file
a file-like object (stream); defaults to the current sys.stdout.
flush
whether to forcibly flush the stream.
Type: builtin_function_or_method
See also:
• IPython user manual: https://ipython.readthedocs.io/en/stable/
• Jupyter Notebook QuickStart: https://docs.jupyter.org/en/latest/start/index.html
As you move forward, it will be important to not only work interactively, but also to create and reuse
Python files. For this, a powerful code editor will get you far. Here are several good easy-to-use editors:
• Spyder: integrates an IPython console, a debugger, a profiler. . .
• PyCharm: integrates an IPython console, notebooks, a debugger. . . (freely available, but commer-
cial)
• Visual Studio Code: integrates a Python console, notebooks, a debugger, . . .
Some of these are shipped by the various scientific Python distributions, and you can find them in the
menus.
As an exercise, create a file my_file.py in a code editor, and add the following lines:
s = 'Hello world'
print(s)
Now, you can run it in IPython console or a notebook and explore the resulting variables:
In [4]: s
Out[4]: 'Hello world'
In [5]: %whos
Variable Type Data/Info
----------------------------
s str Hello world
While it is tempting to work only with scripts, that is a file full of instructions following each other,
do plan to progressively evolve the script to a set of functions:
• A script is not reusable, functions are.
• Thinking in terms of functions helps breaking the problem in small blocks.
The user manuals contain a wealth of information. Here we give a quick introduction to four useful
features: history, tab completion, magic functions, and aliases.
Command history Like a UNIX shell, the IPython console supports command history. Type up and
down to navigate previously typed commands:
In [6]: x = 10
In [7]: <UP>
In [8]: x = 10
Tab completion Tab completion, is a convenient way to explore the structure of any object you’re
dealing with. Simply type object_name.<TAB> to view the object’s attributes. Besides Python objects
and keywords, tab completion also works on file and directory names.*
In [9]: x = 10
In [10]: x.<TAB>
as_integer_ratio() conjugate() imag to_bytes()
bit_count() denominator numerator
bit_length() from_bytes() real
Magic functions The console and the notebooks support so-called magic functions by prefixing a
command with the % character. For example, the run and whos functions from the previous section are
magic functions. Note that, the setting automagic, which is enabled by default, allows you to omit the
preceding % sign. Thus, you can just type the magic function and it will work.
Other useful magic functions are:
• %cd to change the current directory.
In [11]: cd /tmp
/tmp
• %cpaste allows you to paste code, especially code from websites which has been prefixed with the
standard Python prompt (e.g. >>>) or with an ipython prompt, (e.g. in [3]):
In [12]: %cpaste
• %timeit allows you to time the execution of short snippets using the timeit module from the
standard library:
In [12]: %timeit x = 10
11.6 ns +- 0.818 ns per loop (mean +- std. dev. of 7 runs, 100,000,000 loops␣
˓→each)
See also:
Chapter on optimizing code
• %debug allows you to enter post-mortem debugging. That is to say, if the code you try to execute,
raises an exception, using %debug will enter the debugger at the point where the exception was
thrown.
In [13]: x === 10
Cell In[13], line 1
x === 10
^
SyntaxError: invalid syntax
In [14]: %debug
> /home/jarrod/.venv/lectures/lib64/python3.11/site-packages/IPython/core/
˓→compilerop.py(86)ast_parse()
87
(continues on next page)
˓→': 'exec'}
ipdb>
See also:
Chapter on debugging
Aliases Furthermore IPython ships with various aliases which emulate common UNIX command line
tools such as ls to list files, cp to copy files and rm to remove files (a full list of aliases is shown when
typing alias).
Getting help
We introduce here the Python language. Only the bare minimum necessary for getting started with
NumPy and SciPy is addressed here. To learn more about the language, consider going through the
excellent tutorial https://docs.python.org/3/tutorial. Dedicated books are also available, such as Dive
into Python 3.
Tip: Python is a programming language, as are C, Fortran, BASIC, PHP, etc. Some specific
features of Python are as follows:
• an interpreted (as opposed to compiled) language. Contrary to e.g. C or Fortran, one does not
compile Python code before executing it. In addition, Python can be used interactively: many
Python interpreters are available, from which commands and scripts can be executed.
• a free software released under an open-source license: Python can be used and distributed free
of charge, even for building commercial software.
• multi-platform: Python is available for all major operating systems, Windows, Linux/Unix,
MacOS X, most likely your mobile phone OS, etc.
• a very readable language with clear non-verbose syntax
12
Scientific Python Lectures, Edition 2024.2rc0.dev0
• a language for which a large variety of high-quality packages are available for various applications,
from web frameworks to scientific computing.
• a language very easy to interface with other languages, in particular C and C++.
• Some other features of the language are illustrated just below. For example, Python is an object-
oriented language, with dynamic typing (the same variable can contain objects of different types
during the course of a program).
See https://www.python.org/about/ for more information about distinguishing features of Python.
Tip: If you don’t have Ipython installed on your computer, other Python shells are available, such
as the plain Python shell started by typing “python” in a terminal, or the Idle interpreter. However,
we advise to use the Ipython shell because of its enhanced features, especially for interactive scientific
computing.
Tip: The message “Hello, world!” is then displayed. You just executed your first Python instruction,
congratulations!
>>> a = 3
>>> b = 2*a
>>> type(b)
<class 'int'>
>>> print(b)
6
>>> a*b
18
>>> b = 'hello'
>>> type(b)
<class 'str'>
>>> b + b
'hellohello'
>>> 2*b
'hellohello'
Tip: Two variables a and b have been defined above. Note that one does not declare the type of a
variable before assigning its value. In C, conversely, one should write:
int a = 3;
In addition, the type of a variable may change, in the sense that at one point in time it can be equal
to a value of a certain type, and a second point in time, it can be equal to a value of a different type.
b was first equal to an integer, but it became equal to a string when it was assigned the value ‘hello’.
Operations on integers (b=2*a) are coded natively in Python, and so are some operations on strings such
as additions and multiplications, which amount respectively to concatenation and repetition.
Integer
>>> 1 + 1
2
>>> a = 4
>>> type(a)
<class 'int'>
Floats
>>> c = 2.1
>>> type(c)
<class 'float'>
Complex
Booleans
>>> 3 > 4
False
>>> test = (3 > 4)
>>> test
False
>>> type(test)
<class 'bool'>
Tip: A Python shell can therefore replace your pocket calculator, with the basic arithmetic operations
+, -, *, /, % (modulo) natively implemented
>>> 7 * 3.
21.0
>>> 2**10
1024
>>> 8 % 3
2
>>> float(1)
1.0
2.2.2 Containers
Tip: Python provides many efficient types of containers, in which collections of objects can be stored.
Lists
Tip: A list is an ordered collection of objects, that may have different types. For example:
>>> colors[2]
'green'
>>> colors[-1]
'white'
>>> colors[-2]
'black'
>>> colors
['red', 'blue', 'green', 'black', 'white']
>>> colors[2:4]
['green', 'black']
Warning: Note that colors[start:stop] contains the elements with indices i such as start<=
i < stop (i ranging from start to stop-1). Therefore, colors[start:stop] has (stop - start)
elements.
Tip: For collections of numerical data that all have the same type, it is often more efficient to use the
array type provided by the numpy module. A NumPy array is a chunk of memory containing fixed-sized
items. With NumPy arrays, operations on elements can be faster because elements are regularly spaced
in memory and more operations are performed through specialized C functions instead of Python loops.
Tip: Python offers a large panel of functions to modify lists, or query them. Here are a few examples;
for more details, see https://docs.python.org/3/tutorial/datastructures.html#more-on-lists
Reverse:
>>> rcolors2
['red', 'blue', 'green', 'black', 'white']
>>> rcolors2.reverse() # in-place; reversing rcolors2 does not affect colors
>>> rcolors2
['white', 'black', 'green', 'blue', 'red']
Tip: Sort:
The notation rcolors.method() (e.g. rcolors.append(3) and colors.pop()) is our first example
of object-oriented programming (OOP). Being a list, the object rcolors owns the method function
that is called using the notation .. No further knowledge of OOP than understanding the notation .
is necessary for going through this tutorial.
Discovering methods:
Strings
This syntax error can be avoided by enclosing the string in double quotes instead of single quotes.
Alternatively, one can prepend a backslash to the second single quote. Other uses of the backslash are,
e.g., the newline character \n and the tab character \t.
Tip: Strings are collections like lists. Hence they can be indexed and sliced, using the same syntax and
rules.
Indexing:
>>> a = "hello"
>>> a[0]
'h'
>>> a[1]
'e'
>>> a[-1]
'o'
Tip: (Remember that negative indices correspond to counting from the right end.)
Slicing:
Tip: Accents and special characters can also be handled as in Python 3 strings consist of Unicode
characters.
A string is an immutable object and it is not possible to modify its contents. One may however create
new strings from the original one.
Tip: Strings have many useful methods, such as a.replace as seen above. Remember the a. object-
oriented notation and use tab completion or help(str) to search for new methods.
See also:
Python offers advanced possibilities for manipulating strings, looking for patterns or formatting. The
interested reader is referred to https://docs.python.org/3/library/stdtypes.html#string-methods and
https://docs.python.org/3/library/string.html#format-string-syntax
String formatting:
>>> 'An integer: %i ; a float: %f ; another string: %s ' % (1, 0.1, 'string') # with␣
˓→more values use tuple after %
>>> i = 102
>>> filename = 'processing_of_dataset_%d .txt' % i # no need for tuples with just␣
˓→one value after %
>>> filename
'processing_of_dataset_102.txt'
Dictionaries
Tip: It can be used to conveniently store and retrieve values associated with a name (a string for a
date, a name, etc.). See https://docs.python.org/3/tutorial/datastructures.html#dictionaries for more
information.
A dictionary can have keys (resp. values) with different types:
Tuples
Tuples are basically immutable lists. The elements of a tuple are written between parentheses, or just
separated by commas:
Things to note:
• a single object can have several names bound to it:
In [5]: a = [1, 2, 3]
In [6]: a = [1, 2, 3]
2.3.1 if/elif/else
>>> if 2**2 == 4:
... print("Obvious!")
...
Obvious!
Tip: Type the following lines in your Python interpreter, and be careful to respect the indentation
depth. The Ipython shell automatically increases the indentation depth after a colon : sign; to decrease
the indentation depth, go four spaces to the left with the Backspace key. Press the Enter key twice to
leave the logical block.
>>> a = 10
>>> if a == 1:
... print(1)
... elif a == 2:
... print(2)
... else:
... print("A lot")
...
A lot
Indentation is compulsory in scripts as well. As an exercise, re-type the previous lines with the same
indentation in a script condition.py, and execute the script with run condition.py in Ipython.
2.3.2 for/range
2.3.3 while/break/continue
>>> z = 1 + 1j
>>> while abs(z) < 100:
... z = z**2 + 1
>>> z
(-134+352j)
>>> z = 1 + 1j
>>> a = [1, 0, 2, 4]
>>> for element in a:
... if element == 0:
... continue
... print(1. / element)
1.0
0.5
0.25
if <OBJECT>
Evaluates to False:
• any number equal to zero (0, 0.0, 0+0j)
• an empty container (list, tuple, set, dictionary, . . . )
• False, None
Evaluates to True:
• everything else
a == b
Tests equality, with logics:
>>> 1 == 1.
True
a is b
Tests identity: both sides are the same object:
>>> a = 1
>>> b = 1.
>>> a == b
True
(continues on next page)
>>> a = 1
>>> b = 1
>>> a is b
True
a in b
For any collection b: b contains a
>>> b = [1, 2, 3]
>>> 2 in b
True
>>> 5 in b
False
You can iterate over any sequence (string, list, keys in a dictionary, lines in a file, . . . ):
Tip: Few languages (in particular, languages for scientific computing) allow to loop over anything but
integers/indices. With Python it is possible to loop exactly over the objects of interest without bothering
with indices you often don’t care about. This feature can often be used to make code more readable.
Warning: Not safe to modify the sequence you are iterating over.
Common task is to iterate over a sequence while keeping track of the item number.
• Could use while loop with a counter as above. Or a for loop:
Use items:
Note: The ordering of a dictionary is random, thus we use sorted() which will sort on the keys.
Instead of creating a list by means of a loop, one can make use of a list comprehension with a rather
self-explaining syntax.
Exercise
In [2]: test()
in test function
In [4]: disk_area(1.5)
Out[4]: 7.0649999999999995
2.4.3 Parameters
In [6]: double_it(3)
Out[6]: 6
In [7]: double_it()
---------------------------------------------------------------------------
(continues on next page)
In [9]: double_it()
Out[9]: 4
In [10]: double_it(3)
Out[10]: 6
Warning: Default values are evaluated when the function is defined, not when it is called. This
can be problematic when using mutable types (e.g. dictionary or list) and modifying them in the
function body, since the modifications will be persistent across invocations of the function.
Using an immutable type in a keyword argument:
In [11]: bigx = 10
In [14]: double_it()
Out[14]: 20
Using an mutable type in a keyword argument (and modifying it inside the function body):
In [15]: def add_to_dict(args={'a': 1, 'b': 2}):
....: for i in args.keys():
....: args[i] += 1
....: print(args)
....:
In [16]: add_to_dict
Out[16]: <function __main__.add_to_dict(args={'a': 1, 'b': 2})>
In [17]: add_to_dict()
{'a': 2, 'b': 3}
In [18]: add_to_dict()
{'a': 3, 'b': 4}
In [19]: add_to_dict()
{'a': 4, 'b': 5}
In [21]: rhyme = 'one fish, two fish, red fish, blue fish'.split()
In [22]: rhyme
Out[22]: ['one', 'fish,', 'two', 'fish,', 'red', 'fish,', 'blue', 'fish']
In [23]: slicer(rhyme)
Out[23]: ['one', 'fish,', 'two', 'fish,', 'red', 'fish,', 'blue', 'fish']
but it is good practice to use the same ordering as the function’s definition.
Keyword arguments are a very convenient feature for defining functions with a variable number of argu-
ments, especially when default values are to be used in most calls to the function.
Tip: Can you modify the value of a variable inside a function? Most languages (C, Java, . . . ) distinguish
“passing by value” and “passing by reference”. In Python, such a distinction is somewhat artificial, and
it is a bit subtle whether your variables are going to be modified or not. Fortunately, there exist clear
rules.
Parameters to functions are references to objects, which are passed by value. When you pass a variable
to a function, python passes the reference to the object to which the variable refers (the value). Not the
variable itself.
If the value passed in a function is immutable, the function does not modify the caller’s variable. If the
value is mutable, the function may modify the caller’s variable in-place:
>>> def try_to_modify(x, y, z):
... x = 23
... y.append(42)
... z = [99] # new reference
... print(x)
... print(y)
... print(z)
(continues on next page)
Variables declared outside the function can be referenced within the function:
In [28]: x = 5
In [30]: addx(10)
Out[30]: 15
But these “global” variables cannot be modified within the function, unless declared global in the
function.
This doesn’t work:
In [32]: setx(10)
x is 10
In [33]: x
Out[33]: 5
This works:
In [35]: setx(10)
x is 10
In [36]: x
Out[36]: 10
2.4.7 Docstrings
Documentation about what the function does and its parameters. General convention:
In [40]: funcname?
Signature: funcname(params)
Docstring:
Concise one-line sentence describing the function.
Extended summary which can contain multiple paragraphs.
File: ~/src/scientific-python-lectures/<ipython-input-13-64e466df6d64>
Type: function
In [41]: va = variable_args
2.4.9 Methods
Methods are functions attached to objects. You’ve seen these in our examples on lists, dictionaries,
strings, etc. . .
2.4.10 Exercises
Write a function that displays the n first terms of the Fibonacci sequence, defined by:
⎧
⎨ 𝑈0 = 0
𝑈1 = 1
𝑈𝑛+2 = 𝑈𝑛+1 + 𝑈𝑛
⎩
Exercise: Quicksort
function quicksort(array)
var list less, greater
if length(array) < 2
return array
select and remove a pivot value pivot from array
for each x in array
if x < pivot + 1 then append x to less
else append x to greater
return concatenate(quicksort(less), pivot, quicksort(greater))
For now, we have typed all instructions in the interpreter. For longer sets of instructions we need to
change track and write the code in text files (using a text editor), that we will call either scripts or
modules. Use your favorite text editor (provided it offers syntax highlighting for Python), or the editor
that comes with the Scientific Python Suite you may be using.
2.5.1 Scripts
Tip: Let us first write a script, that is a file with a sequence of instructions that are executed each
time the script is called. Instructions may be e.g. copied-and-pasted from the interpreter (but take care
to respect indentation rules!).
The extension for Python files is .py. Write or copy-and-paste the following lines in a file called test.py
Tip: Let us now execute the script interactively, that is inside the Ipython interpreter. This is maybe
the most common use of scripts in scientific computing.
Note: in Ipython, the syntax to execute a script is %run script.py. For example,
In [2]: message
Out[2]: 'Hello how are you?'
The script has been executed. Moreover the variables defined in the script (such as message) are now
available inside the interpreter’s namespace.
Tip: Other interpreters also offer the possibility to execute scripts (e.g., execfile in the plain Python
interpreter, etc.).
It is also possible In order to execute this script as a standalone program, by executing the script inside
a shell terminal (Linux/Mac console or cmd Windows console). For example, if we are in the same
directory as the test.py file, we can execute this in a console:
$ python test.py
Hello
how
are
you?
In file.py:
import sys
print(sys.argv)
Warning: Don’t implement option parsing yourself. Use a dedicated module such as argparse.
In [3]: import os
In [4]: os
Out[4]: <module 'os' (frozen)>
In [5]: os.listdir('.')
Out[5]:
['profile_ksttxb2h',
'.XIM-unix',
'profile_8jyd3wt0',
'systemd-private-547c90d23384455997d6600060dd5355-chrony.service-Cs0azH',
'profile_i8ku8_ln',
'profile_6c5_7yvi',
'clr-debug-pipe-589-869-in',
'profile_n4tf0uy4',
'profile_7crnvcpm',
'profile_ecjb1ld9',
'profile_m13gl1fv',
'profile_1_590dye',
'profile_v_p76c0x',
'clr-debug-pipe-1621-15129-out',
'profile_lvwe1o1z',
'profile_orrts951',
'profile_bziy8mt_',
'profile_0u5120_y',
'profile_tav_998d',
'.font-unix',
'dotnet-diagnostic-589-869-socket',
'clr-debug-pipe-589-869-out',
'profile_ujz8liq0',
'profile_rqn3apk7',
'profile_f3qz57ee',
'profile_v5h87_cl',
'profile_g8aplujp',
'systemd-private-547c90d23384455997d6600060dd5355-systemd-resolved.service-qwjgaC',
'.Test-unix',
'profile_orda1vei',
'profile_zwgjajcz',
'clr-debug-pipe-1604-14840-in',
'profile_rd6ks_2m',
'profile_hsd6fbf4',
(continues on next page)
And also:
Importing shorthands:
Warning:
from os import *
Tip: Modules are thus a good way to organize code in a hierarchical way. Actually, all the scientific
computing tools we are going to use are modules:
Tip: If we want to write larger and better organized programs (compared to simple scripts), where
some objects are defined, (variables, functions, classes) and that we want to reuse several times, we have
to create our own modules.
def print_b():
"Prints b."
print("b")
def print_a():
"Prints a."
print("a")
c = 2
d = 2
Tip: In this file, we defined two functions print_a and print_b. Suppose we want to call the print_a
function from the interpreter. We could execute the file as a script, but since we just want to have access
to the function print_a, we are rather going to import it as a module. The syntax is as follows.
In [9]: demo.print_a()
a
In [10]: demo.print_b()
b
Importing the module gives access to its objects, using the module.object syntax. Don’t forget to put
the module’s name before the object’s name, otherwise Python won’t recognize the instruction.
Introspection
In [11]: demo?
Type: module
Base Class: <type 'module'>
String Form: <module 'demo' from 'demo.py'>
Namespace: Interactive
File: /home/varoquau/Projects/Python_talks/scipy_2009_tutorial/source/
˓→demo.py
Docstring:
A demo module.
In [12]: who
demo
In [13]: whos
Variable Type Data/Info
------------------------------
demo module <module 'demo' from 'demo.py'>
In [14]: dir(demo)
Out[14]:
['__builtins__',
'__doc__',
'__file__',
'__name__',
'__package__',
'c',
'd',
'print_a',
'print_b']
In [15]: demo.<TAB>
demo.c demo.print_a demo.py
demo.d demo.print_b demo.pyc
In [17]: whos
Variable Type Data/Info
--------------------------------
demo module <module 'demo' from 'demo.py'>
print_a function <function print_a at 0xb7421534>
print_b function <function print_b at 0xb74214c4>
Tip: Sometimes we want code to be executed when a module is run directly, but not when it is imported
by another module. if __name__ == '__main__' allows us to check whether the module is being run
directly.
File demo2.py:
def print_b():
"Prints b."
print("b")
def print_a():
"Prints a."
print("a")
if __name__ == "__main__":
# print_a() is only executed when the module is run directly.
print_a()
Importing it:
Running it:
When the import mymodule statement is executed, the module mymodule is searched in a given list of
directories. This list includes a list of installation-dependent default path (e.g., /usr/lib64/python3.
11) as well as the list of directories specified by the environment variable PYTHONPATH.
The list of directories searched by Python is given by the sys.path variable
In [23]: sys.path
Out[23]:
['/home/runner/work/scientific-python-lectures/scientific-python-lectures',
'/opt/hostedtoolcache/Python/3.11.9/x64/lib/python311.zip',
'/opt/hostedtoolcache/Python/3.11.9/x64/lib/python3.11',
'/opt/hostedtoolcache/Python/3.11.9/x64/lib/python3.11/lib-dynload',
'/opt/hostedtoolcache/Python/3.11.9/x64/lib/python3.11/site-packages']
Tip: On Linux/Unix, add the following line to a file read by the shell at startup (e.g. /etc/profile,
.profile)
export PYTHONPATH=$PYTHONPATH:/home/emma/user_defined_modules
Tip:
import sys
new_path = '/home/emma/user_defined_modules'
if new_path not in sys.path:
sys.path.append(new_path)
This method is not very robust, however, because it makes the code less portable (user-dependent
path) and because you have to add the directory to your sys.path each time you want to import
from a module in this directory.
See also:
See https://docs.python.org/3/tutorial/modules.html for more information about modules.
2.5.6 Packages
A directory that contains many modules is called a package. A package is a module with submodules
(which can have submodules themselves, etc.). A special file called __init__.py (which may be empty)
tells Python that the directory is a Python package, from which modules can be imported.
$ ls
_build_utils/ fft/ _lib/ odr/ spatial/
cluster/ fftpack/ linalg/ optimize/ special/
conftest.py __init__.py linalg.pxd optimize.pxd special.pxd
constants/ integrate/ meson.build setup.py stats/
datasets/ interpolate/ misc/ signal/
_distributor_init.py io/ ndimage/ sparse/
$ cd ndimage
$ ls
_filters.py __init__.py _measurements.py morphology.py src/
filters.py _interpolation.py measurements.py _ni_docstrings.py tests/
_fourier.py interpolation.py meson.build _ni_support.py utils/
fourier.py LICENSE.txt _morphology.py setup.py
From Ipython:
In [25]: sp.__file__
Out[25]: '/opt/hostedtoolcache/Python/3.11.9/x64/lib/python3.11/site-packages/scipy/__
˓→init__.py'
In [26]: sp.version.version
Out[26]: '1.13.1'
In [27]: sp.ndimage.morphology.binary_dilation?
Signature:
sp.ndimage.morphology.binary_dilation(
input,
structure=None,
iterations=1,
mask=None,
output=None,
border_value=0,
origin=0,
brute_force=False,
)
Docstring:
Multidimensional binary dilation with the given structuring element.
...
Tip: Indenting is compulsory in Python! Every command block following a colon bears an
additional indentation level with respect to the previous line with a colon. One must therefore
indent after def f(): or while:. At the end of such logical blocks, one decreases the indentation
depth (and re-increases it if a new block is entered, etc.)
Strict respect of indentation is the price to pay for getting rid of { or ; characters that delineate
logical blocks in other languages. Improper indentation leads to errors such as
------------------------------------------------------------
IndentationError: unexpected indent (test.py, line 2)
All this indentation business can be a bit confusing in the beginning. However, with the clear
indentation, and in the absence of extra characters, the resulting code is very nice to read compared
to other languages.
• Indentation depth: Inside your text editor, you may choose to indent with any positive number
of spaces (1, 2, 3, 4, . . . ). However, it is considered good practice to indent with 4 spaces. You
may configure your editor to map the Tab key to a 4-space indentation.
• Style guidelines
Long lines: you should not write very long lines that span over more than (e.g.) 80 characters.
Long lines can be broken with the \ character
Spaces
Write well-spaced code: put whitespaces after commas, around arithmetic operators, etc.:
>>> a = 1 # yes
>>> a=1 # too cramped
A certain number of rules for writing “beautiful” code (and more importantly using the same
conventions as anybody else!) are given in the Style Guide for Python Code.
Quick read
If you want to do a first quick pass through the Scientific Python Lectures to learn the ecosystem, you
can directly skip to the next chapter: NumPy: creating and manipulating numerical data.
The remainder of this chapter is not necessary to follow the rest of the intro part. But be sure to
come back and finish this chapter later.
To be exhaustive, here are some information about input and output in Python. Since we will use the
NumPy methods to read and write files, you may skip this chapter at first reading.
We write or read strings to/from files (other types must be converted to strings). To write in a file:
In [2]: s = f.read()
In [3]: print(s)
This is a test
and another test
In [4]: f.close()
See also:
For more details: https://docs.python.org/3/tutorial/inputoutput.html
In [7]: f.close()
File modes
• Read-only: r
• Write-only: w
– Note: Create a new file or overwrite existing file.
• Append a file: a
• Read and Write: r+
• Binary mode: b
– Note: Use for binary files, especially on Windows.
Current directory:
In [1]: import os
In [2]: os.getcwd()
Out[2]: '/tmp'
List a directory:
In [3]: os.listdir(os.curdir)
Out[3]:
['profile_ksttxb2h',
'profile_crwqcy2d',
'.XIM-unix',
'profile_8jyd3wt0',
'systemd-private-547c90d23384455997d6600060dd5355-chrony.service-Cs0azH',
'profile_i8ku8_ln',
'profile_6c5_7yvi',
'clr-debug-pipe-589-869-in',
'profile_n4tf0uy4',
'profile_rwsw_k2b',
'profile_7crnvcpm',
'profile_ecjb1ld9',
'profile_m13gl1fv',
'profile_1_590dye',
'profile_v_p76c0x',
'clr-debug-pipe-1621-15129-out',
'profile_lvwe1o1z',
'profile_orrts951',
'profile_bziy8mt_',
'profile_0u5120_y',
'profile_tav_998d',
'.font-unix',
'dotnet-diagnostic-589-869-socket',
'clr-debug-pipe-589-869-out',
'profile_ujz8liq0',
'profile_rqn3apk7',
'profile_f3qz57ee',
'profile_v5h87_cl',
'profile_g8aplujp',
'systemd-private-547c90d23384455997d6600060dd5355-systemd-resolved.service-qwjgaC',
(continues on next page)
Make a directory:
In [4]: os.mkdir('junkdir')
In [9]: os.rmdir('foodir')
Delete a file:
In [12]: fp.close()
In [14]: os.remove('junk.txt')
In [17]: fp.close()
In [18]: a = os.path.abspath('junk.txt')
In [19]: a
Out[19]: '/tmp/junk.txt'
In [20]: os.path.split(a)
(continues on next page)
In [21]: os.path.dirname(a)
Out[21]: '/tmp'
In [22]: os.path.basename(a)
Out[22]: 'junk.txt'
In [23]: os.path.splitext(os.path.basename(a))
Out[23]: ('junk', '.txt')
In [24]: os.path.exists('junk.txt')
Out[24]: True
In [25]: os.path.isfile('junk.txt')
Out[25]: True
In [26]: os.path.isdir('junk.txt')
Out[26]: False
In [27]: os.path.expanduser('~/local')
Out[27]: '/home/runner/local'
In [29]: os.system('ls')
Out[29]: 0
In [30]: import sh
In [31]: com = sh.ls()
In [31]: print(com)
basic_types.rst exceptions.rst oop.rst standard_library.rst
control_flow.rst first_steps.rst python_language.rst
demo2.py functions.rst python-logo.png
demo.py io.rst reusing_code.rst
In [32]: type(com)
Out[32]: str
Walking a directory
Environment variables:
In [34]: os.environ.keys()
Out[34]: KeysView(environ({'SHELL': '/bin/bash', 'COLORTERM': 'truecolor', ...}))
In [35]: os.environ['SHELL']
Out[35]: '/bin/bash'
In [37]: glob.glob('*.txt')
Out[37]: ['junk.txt']
In [39]: sys.platform
Out[39]: 'linux'
In [40]: sys.version
Out[40]: '3.11.9 (main, May 9 2024, 14:13:20) [GCC 11.4.0]'
In [41]: sys.prefix
Out[41]: '/opt/hostedtoolcache/Python/3.11.9/x64'
In [42]: sys.argv
Out[42]:
['/opt/hostedtoolcache/Python/3.11.9/x64/lib/python3.11/site-packages/sphinx/__main__.
˓→py',
'-b',
'latex',
'-d',
'build/doctrees',
'.',
'build/latex']
sys.path is a list of strings that specifies the search path for modules. Initialized from PYTHONPATH:
In [43]: sys.path
Out[43]:
['/home/runner/work/scientific-python-lectures/scientific-python-lectures',
'/opt/hostedtoolcache/Python/3.11.9/x64/lib/python311.zip',
'/opt/hostedtoolcache/Python/3.11.9/x64/lib/python3.11',
(continues on next page)
In [48]: out
Out[48]: [1, None, 'Stan']
Exercise
path_site
It is likely that you have raised Exceptions if you have typed all the previous commands of the tutorial.
For example, you may have raised an exception if you entered a command with a typo.
Exceptions are raised by different kinds of errors arising when executing Python code. In your own code,
you may also catch errors, or define custom error types. You may want to look at the descriptions of the
the built-in Exceptions when looking for the right exception type.
2.8.1 Exceptions
In [1]: 1/0
---------------------------------------------------------------------------
ZeroDivisionError Traceback (most recent call last)
Cell In[1], line 1
----> 1 1/0
In [2]: 1 + 'e'
(continues on next page)
In [4]: d[3]
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
Cell In[4], line 1
----> 1 d[3]
KeyError: 3
In [5]: l = [1, 2, 3]
In [6]: l[4]
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
Cell In[6], line 1
----> 1 l[4]
In [7]: l.foobar
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[7], line 1
----> 1 l.foobar
As you can see, there are different types of exceptions for different errors.
try/except
In [9]: x
Out[9]: 1
try/finally
In [10]: try:
....: x = int(input('Please enter a number: '))
....: finally:
....: print('Thank you for your input')
....:
Please enter a number: a
Thank you for your input
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[10], line 2
1 try:
----> 2 x = int(input('Please enter a number: '))
3 finally:
4 print('Thank you for your input')
ValueError: invalid literal for int() with base 10: 'a'
In [14]: print_sorted('132')
132
In [16]: filter_name('Gaël')
(continues on next page)
In [17]: filter_name('Stéfan')
---------------------------------------------------------------------------
UnicodeEncodeError Traceback (most recent call last)
Cell In[17], line 1
----> 1 filter_name('Stéfan')
In [19]: x = 0
In [21]: x
Out[21]: 0.9990234375
Use exceptions to notify certain conditions are met (e.g. StopIteration) or not (e.g. custom error raising)
In the previous example, the Student class has __init__, set_age and set_major methods. Its at-
tributes are name, age and major. We can call these methods and attributes with the following notation:
classinstance.method or classinstance.attribute. The __init__ constructor is a special method
we call with: MyClass(init parameters if any).
Now, suppose we want to create a new class MasterStudent with the same methods and attributes as
the previous one, but with an additional internship attribute. We won’t copy the previous class, but
inherit from it:
The MasterStudent class inherited from the Student attributes and methods.
Thanks to classes and object-oriented programming, we can organize code with different classes corre-
sponding to different objects we encounter (an Experiment class, an Image class, a Flow class, etc.), with
their own methods and attributes. Then we can use inheritance to consider variations around a base
class and reuse code. Ex : from a Flow base class, we can create derived StokesFlow, TurbulentFlow,
PotentialFlow, etc.
Authors: Emmanuelle Gouillart, Didrik Pinte, Gaël Varoquaux, and Pauli Virtanen
This chapter gives an overview of NumPy, the core tool for performant numerical computing with Python.
Section contents
53
Scientific Python Lectures, Edition 2024.2rc0.dev0
NumPy arrays
Python objects
• high-level number objects: integers, floating point
• containers: lists (costless insertion and append), dictionaries (fast lookup)
NumPy provides
• extension package to Python for multi-dimensional arrays
• closer to hardware (efficiency)
• designed for scientific computation (convenience)
• Also known as array oriented computing
In [1]: L = range(1000)
In [3]: a = np.arange(1000)
In [5]: np.array?
Docstring:
array(object, dtype=None, *, copy=True, order='K', subok=False, ndmin=0,
like=None)
Create an array.
Parameters
----------
object : array_like
An array, any object exposing the array interface, an object whose
``__array__`` method returns an array, or any (nested) sequence.
If object is a scalar, a 0-dimensional array containing object is
returned.
dtype : data-type, optional
The desired data-type for the array. If not given, NumPy will try to use
a default ``dtype`` that can represent the values (by applying promotion
rules when necessary.)
copy : bool, optional
If true (default), then the object is copied. Otherwise, a copy will
only be made if ``__array__`` returns a copy, if obj is a nested
sequence, or if a copy is needed to satisfy any of the other
requirements (``dtype``, ``order``, etc.).
order : {'K', 'A', 'C', 'F'}, optional
Specify the memory layout of the array. If object is not an array, the
newly created array will be in C order (row major) unless 'F' is
specified, in which case it will be in Fortran order (column major).
If object is an array the following holds.
When ``copy=False`` and a copy is made for other reasons, the result is
the same as if ``copy=True``, with some exceptions for 'A', see the
Notes section. The default order is 'K'.
subok : bool, optional
If True, then sub-classes will be passed-through, otherwise
the returned array will be forced to be a base-class array (default).
ndmin : int, optional
Specifies the minimum number of dimensions that the resulting
array should have. Ones will be prepended to the shape as
needed to meet this requirement.
like : array_like, optional
Reference object to allow the creation of arrays which are not
NumPy arrays. If an array-like passed in as ``like`` supports
the ``__array_function__`` protocol, the result will be defined
(continues on next page)
.. versionadded:: 1.20.0
Returns
-------
out : ndarray
An array object satisfying the specified requirements.
See Also
--------
empty_like : Return an empty array with shape and type of input.
ones_like : Return an array of ones with shape and type of input.
zeros_like : Return an array of zeros with shape and type of input.
full_like : Return a new array with shape of input filled with value.
empty : Return a new uninitialized array.
ones : Return a new array setting values to one.
zeros : Return a new array setting values to zero.
full : Return a new array of given shape filled with value.
Notes
-----
When order is 'A' and ``object`` is an array in neither 'C' nor 'F' order,
and a copy is forced by a change in dtype, then the order of the result is
not necessarily 'C' as expected. This is likely a bug.
Examples
--------
>>> np.array([1, 2, 3])
array([1, 2, 3])
Upcasting:
Minimum dimensions 2:
Type provided:
In [6]: np.con*?
np.concatenate
np.conj
np.conjugate
np.convolve
Import conventions
• 1-D:
• 2-D, 3-D, . . . :
[[3],
[4]]])
>>> c.shape
(2, 2, 1)
• Create a simple two dimensional array. First, redo the examples from above. And then create
your own: how about odd numbers counting backwards on the first row, and even numbers on
the second?
• Use the functions len(), numpy.shape() on these arrays. How do they relate to each other?
And to the ndim attribute of the arrays?
• Evenly spaced:
• or by number of points:
• Common arrays:
You may have noticed that, in some instances, array elements are displayed with a trailing dot (e.g. 2.
vs 2). This is due to a difference in the data-type used:
Tip: Different data-types allow us to store data more compactly in memory, but most of the time we
simply work with floating point numbers. Note that, in the example above, NumPy auto-detects the
data-type from the input.
Bool
Strings
Much more
• int32
• int64
• uint32
• uint64
Now that we have our first data arrays, we are going to visualize them.
Start by launching IPython:
Or the notebook:
$ jupyter notebook
>>> %matplotlib
The inline is important for the notebook, so that plots are displayed in the notebook and not in a new
window.
Matplotlib is a 2D plotting package. We can import its functions as below:
And then use (note that you have to use show explicitly if you have not enabled interactive plots with
%matplotlib):
• 1D plotting:
See also:
More in the: matplotlib chapter
The items of an array can be accessed and assigned to the same way as other Python sequences (e.g.
lists):
>>> a = np.arange(10)
>>> a
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> a[0], a[2], a[-1]
(0, 2, 9)
Warning: Indices begin at 0, like other Python sequences (and C/C++). In contrast, in Fortran
or Matlab, indices begin at 1.
>>> a[::-1]
array([9, 8, 7, 6, 5, 4, 3, 2, 1, 0])
>>> a = np.diag(np.arange(3))
>>> a
array([[0, 0, 0],
[0, 1, 0],
[0, 0, 2]])
>>> a[1, 1]
1
>>> a[2, 1] = 10 # third line, second column
>>> a
(continues on next page)
Note:
• In 2D, the first dimension corresponds to rows, the second to columns.
• for multidimensional a, a[0] is interpreted by taking all elements in the unspecified dimensions.
>>> a[:4]
array([0, 1, 2, 3])
All three slice components are not required: by default, start is 0, end is the last and step is 1:
>>> a[1:3]
array([1, 2])
>>> a[::2]
array([0, 2, 4, 6, 8])
>>> a[3:]
array([3, 4, 5, 6, 7, 8, 9])
• Try the different flavours of slicing, using start, end and step: starting from a linspace, try to
obtain odd numbers counting backwards, and even numbers counting forwards.
• Reproduce the slices in the diagram above. You may use the following expression to create the
array:
>>> np.arange(6) + np.arange(0, 51, 10)[:, np.newaxis]
array([[ 0, 1, 2, 3, 4, 5],
[10, 11, 12, 13, 14, 15],
[20, 21, 22, 23, 24, 25],
[30, 31, 32, 33, 34, 35],
[40, 41, 42, 43, 44, 45],
[50, 51, 52, 53, 54, 55]])
Skim through the documentation for np.tile, and use this function to construct the array:
[[4, 3, 4, 3, 4, 3],
[2, 1, 2, 1, 2, 1],
[4, 3, 4, 3, 4, 3],
[2, 1, 2, 1, 2, 1]]
A slicing operation creates a view on the original array, which is just a way of accessing array data.
Thus the original array is not copied in memory. You can use np.may_share_memory() to check if two
arrays share the same memory block. Note however, that this uses heuristics and may give you false
positives.
When modifying the view, the original array is modified as well:
>>> a = np.arange(10)
>>> a
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> b = a[::2]
>>> b
array([0, 2, 4, 6, 8])
>>> np.may_share_memory(a, b)
True
>>> b[0] = 12
>>> b
array([12, 2, 4, 6, 8])
>>> a # (!)
array([12, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> a = np.arange(10)
>>> c = a[::2].copy() # force a copy
>>> c[0] = 12
>>> a
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> np.may_share_memory(a, c)
False
This behavior can be surprising at first sight. . . but it allows to save both memory and time.
• For each integer j starting from 2, cross out its higher multiples:
>>> N_max = int(np.sqrt(len(is_prime) - 1))
>>> for j in range(2, N_max + 1):
... is_prime[2*j::j] = False
Tip: NumPy arrays can be indexed with slices, but also with boolean or integer arrays (masks). This
method is called fancy indexing. It creates copies not views.
Indexing with a mask can be very useful to assign a new value to a sub-array:
>>> a[a % 3 == 0] = -1
>>> a
array([-1, 13, -1, 10, 10, 10, -1, 4, 8, 5, -1, 11, -1, 17, -1])
Indexing can be done with an array of integers, where the same index is repeated several time:
Tip: When a new array is created by indexing with an array of integers, the new array has the same
shape as the array of integers:
>>> a = np.arange(10)
>>> idx = np.array([[3, 4], [9, 7]])
>>> idx.shape
(2, 2)
>>> a[idx]
array([[3, 4],
[9, 7]])
Section contents
• Elementwise operations
• Basic reductions
• Broadcasting
• Array shape manipulation
• Sorting data
• Summary
Basic operations
With scalars:
>>> b = np.ones(4) + 1
>>> a - b
array([-1., 0., 1., 2.])
>>> a * b
array([2., 4., 6., 8.])
>>> j = np.arange(5)
>>> 2**(j + 1) - j
array([ 2, 3, 6, 13, 28])
These operations are of course much faster than if you did them in pure python:
>>> a = np.arange(10000)
>>> %timeit a + 1
10000 loops, best of 3: 24.3 us per loop
>>> l = range(10000)
>>> %timeit [i+1 for i in l]
1000 loops, best of 3: 861 us per loop
>>> c @ c
array([[3., 3., 3.],
[3., 3., 3.],
[3., 3., 3.]])
• Try simple arithmetic elementwise operations: add even elements with odd elements
• Time them against their pure python counterparts using %timeit.
• Generate:
– [2**0, 2**1, 2**2, 2**3, 2**4]
– a_j = 2^(3*j) - j
Other operations
Comparisons:
Logical operations:
Transcendental functions:
>>> a = np.arange(5)
>>> np.sin(a)
array([ 0. , 0.84147098, 0.90929743, 0.14112001, -0.7568025 ])
(continues on next page)
Shape mismatches
>>> a = np.arange(4)
>>> a + np.array([1, 2])
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
ValueError: operands could not be broadcast together with shapes (4,) (2,)
>>> a = np.arange(9).reshape(3, 3)
>>> a.T[0, 2] = 999
>>> a.T
array([[ 0, 3, 999],
[ 1, 4, 7],
[ 2, 5, 8]])
>>> a
array([[ 0, 1, 2],
[ 3, 4, 5],
[999, 7, 8]])
Computing sums
Other reductions
Logical operations:
Statistics:
Exercise: Reductions
• Given there is a sum, what other function might you expect to see?
• What is the difference between sum and cumsum?
Tip: Let us consider a simple 1D random walk process: at each time step a walker jumps right or
left with equal probability.
We are interested in finding the typical distance from the origin of a random walker after t left or
right jumps? We are going to simulate many “walkers” to find this law, and we are going to do so
using array computing tricks: we are going to create a 2D array with the “stories” (each walker has a
story) in one direction, and the time in the other:
We find a well-known result in physics: the RMS distance grows as the square root of the time!
3.2.3 Broadcasting
Let’s verify:
A useful trick:
>>> a = np.arange(0, 40, 10)
>>> a.shape
(4,)
>>> a = a[:, np.newaxis] # adds a new axis -> 2D array
>>> a.shape
(4, 1)
>>> a
array([[ 0],
[10],
[20],
[30]])
>>> a + b
array([[ 0, 1, 2],
[10, 11, 12],
[20, 21, 22],
[30, 31, 32]])
Tip: Broadcasting seems a bit magical, but it is actually quite natural to use it when we want to solve
a problem whose output data is an array with more dimensions than input data.
Let’s construct an array of distances (in miles) between cities of Route 66: Chicago, Springfield,
Saint-Louis, Tulsa, Oklahoma City, Amarillo, Santa Fe, Albuquerque, Flagstaff and Los Angeles.
>>> mileposts = np.array([0, 198, 303, 736, 871, 1175, 1475, 1544,
... 1913, 2448])
>>> distance_array = np.abs(mileposts - mileposts[:, np.newaxis])
>>> distance_array
array([[ 0, 198, 303, 736, 871, 1175, 1475, 1544, 1913, 2448],
[ 198, 0, 105, 538, 673, 977, 1277, 1346, 1715, 2250],
A lot of grid-based or network-based problems can also use broadcasting. For instance, if we want to
compute the distance from the origin of points on a 5x5 grid, we can do
Or in color:
>>> plt.pcolor(distance)
<matplotlib.collections.PolyQuadMesh object at ...>
>>> plt.colorbar()
<matplotlib.colorbar.Colorbar object at ...>
Remark : the numpy.ogrid() function allows to directly create vectors x and y of the previous example,
with two “significant dimensions”:
Tip: So, np.ogrid is very useful as soon as we have to handle computations on a grid. On the other
hand, np.mgrid directly provides matrices full of indices for cases where we can’t (or don’t want to)
benefit from broadcasting:
See also:
Broadcasting: discussion of broadcasting in the Advanced NumPy chapter.
Flattening
Reshaping
>>> a.shape
(2, 3)
>>> b = a.ravel()
>>> b = b.reshape((2, 3))
>>> b
array([[1, 2, 3],
[4, 5, 6]])
Or,
Tip:
>>> b[0, 0] = 99
>>> a
array([[99, 2, 3],
[ 4, 5, 6]])
To understand this you need to learn more about the memory layout of a NumPy array.
Adding a dimension
Indexing with the np.newaxis object allows us to add an axis to an array (you have seen this already
above in the broadcasting section):
>>> z[np.newaxis, :]
array([[1, 2, 3]])
Dimension shuffling
>>> a = np.arange(4*3*2).reshape(4, 3, 2)
>>> a.shape
(4, 3, 2)
>>> a[0, 2, 1]
5
>>> b = a.transpose(1, 2, 0)
>>> b.shape
(3, 2, 4)
>>> b[2, 1, 0]
5
>>> b[2, 1, 0] = -1
>>> a[0, 2, 1]
-1
Resizing
>>> a = np.arange(4)
>>> a.resize((8,))
>>> a
array([0, 1, 2, 3, 0, 0, 0, 0])
>>> b = a
>>> a.resize((4,))
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
ValueError: cannot resize an array that references or is referenced
by another array in this way.
Use the np.resize function or refcheck=False
• Look at the docstring for reshape, especially the notes section which has some more information
about copies and views.
• Use flatten as an alternative to ravel. What is the difference? (Hint: check which one returns
a view and which a copy)
• Experiment with transpose for dimension shuffling.
In-place sort:
>>> a.sort(axis=1)
>>> a
array([[3, 4, 5],
[1, 1, 2]])
Exercise: Sorting
• Look at the axis keyword for sort and rewrite the previous exercise.
3.2.6 Summary
• Know miscellaneous operations on arrays, such as finding the mean or max (array.max(), array.
mean()). No need to retain everything, but have the reflex to search in the documentation (online
docs, help(), lookfor())!!
• For advanced use: master the indexing with arrays of integers, as well as broadcasting. Know more
NumPy functions to handle various array operations.
Quick read
If you want to do a first quick pass through the Scientific Python Lectures to learn the ecosystem, you
can directly skip to the next chapter: Matplotlib: plotting.
The remainder of this chapter is not necessary to follow the rest of the intro part. But be sure to
come back and finish this chapter, as well as to do some more exercices.
Section contents
Casting
Forced casts:
Rounding:
Integers (signed):
int8 8 bits
int16 16 bits
int32 32 bits (same as int on 32-bit platform)
int64 64 bits (same as int on 64-bit platform)
Unsigned integers:
uint8 8 bits
uint16 16 bits
uint32 32 bits
uint64 64 bits
Floating-point numbers:
float16 16 bits
float32 32 bits
float64 64 bits (same as float)
float96 96 bits, platform-dependent (same as np.longdouble)
float128 128 bits, platform-dependent (same as np.longdouble)
>>> np.finfo(np.float32).eps
1.1920929e-07
>>> np.finfo(np.float64).eps
2.2204460492503131e-16
If you don’t know you need special data types, then you probably don’t.
Comparison on using float32 instead of float64:
• Half the size in memory and on disk
• Half the memory bandwidth required (may be a bit faster in some operations)
In [1]: a = np.zeros((int(1e6),), dtype=np.float64)
• But: bigger rounding errors — sometimes in surprising places (i.e., don’t use them unless you
really need them)
>>> samples['sensor_code']
array([b'ALFA', b'BETA', b'TAU', b'ALFA', b'ALFA', b'TAU'], dtype='|S4')
>>> samples['value']
array([0.37, 0.11, 0.13, 0.37, 0.11, 0.13])
>>> samples[0]
(b'ALFA', 1., 0.37)
Note: There are a bunch of other syntaxes for constructing structured arrays, see here and here.
• For floats one could use NaN’s, but masks work for all types:
While it is off topic in a chapter on NumPy, let’s take a moment to recall good coding practice, which
really do pay off in the long run:
Good practices
• Explicit variable names (no need of a comment to explain what is in the variable)
• Style: spaces after commas, around =, etc.
A certain number of rules for writing “beautiful” code (and, more importantly, using the same
conventions as everybody else!) are given in the Style Guide for Python Code and the Docstring
Conventions page (to manage help strings).
• Except some rare cases, variable names and comments in English.
Section contents
• Polynomials
• Loading data files
3.4.1 Polynomials
>>> t = np.linspace(0, 1, 200) # use a larger number of points for smoother plotting
>>> plt.plot(x, y, 'o', t, p(t), '-')
[<matplotlib.lines.Line2D object at ...>, <matplotlib.lines.Line2D object at ...>]
NumPy also has a more sophisticated polynomial interface, which supports e.g. the Chebyshev basis.
3𝑥2 + 2𝑥 − 1:
Example using polynomials in Chebyshev basis, for polynomials in range [-1, 1]:
Text files
Example: populations.txt:
# year hare lynx carrot
1900 30e3 4e3 48300
1901 47.2e3 6.1e3 48200
1902 70.2e3 9.8e3 41500
1903 77.4e3 35.2e3 38200
Note: If you have a complicated text file, what you can try are:
• np.genfromtxt
• Using Python’s I/O functions and e.g. regexps for parsing (Python is quite well suited for this)
Images
Using Matplotlib:
>>> plt.imshow(plt.imread('red_elephant.png'))
<matplotlib.image.AxesImage object at ...>
Other libraries:
NumPy has its own binary format, not portable but with efficient I/O:
Write a Python script that loads data from populations.txt:: and drop the last column and the first
5 rows. Save the smaller dataset to pop2.txt.
NumPy internals
If you are interested in the NumPy internals, there is a good discussion in Advanced NumPy.
[[1, 6, 11],
[2, 7, 12],
[3, 8, 13],
[4, 9, 14],
[5, 10, 15]]
and generate a new array containing its 2nd and 4th rows.
2. Divide each column of the array:
elementwise with the array b = np.array([1., 5, 10, 15, 20]). (Hint: np.newaxis).
3. Harder one: Generate a 10 x 3 array of random numbers (in range [0,1]). For each row, pick the
number closest to 0.5.
• Use abs and argmin to find the column j closest for each row.
• Use fancy indexing to extract the numbers. (Hint: a[i,j] – the array i must contain the
row numbers corresponding to stuff in j.)
Let’s do some manipulations on NumPy arrays by starting with an image of a raccoon. scipy provides
a 2D array of this image with the scipy.datasets.face function:
Here are a few images we will be able to obtain with our manipulations: use different colormaps, crop
the image, change some parts of the image.
• We will now frame the face with a black locket. For this, we
need to create a mask corresponding to the pixels we want to be black. The center of the face
is around (660, 330), so we defined the mask by this condition (y-300)**2 + (x-660)**2
>>> sy, sx = face.shape
>>> y, x = np.ogrid[0:sy, 0:sx] # x and y indices of pixels
>>> y.shape, x.shape
((768, 1), (1, 1024))
>>> centerx, centery = (660, 300) # center of the image
>>> mask = ((y - centery)**2 + (x - centerx)**2) > 230**2 # circle
then we assign the value 0 to the pixels of the image corresponding to the mask. The syntax
is extremely simple and intuitive:
>>> face[mask] = 0
>>> plt.imshow(face)
<matplotlib.image.AxesImage object at 0x...>
The data in populations.txt describes the populations of hares and lynxes (and carrots) in northern
Canada during 20 years:
>>> data = np.loadtxt('data/populations.txt')
>>> year, hares, lynxes, carrots = data.T # trick: columns to variables
Write a function f(a, b, c) that returns 𝑎𝑏 −𝑐. Form a 24x12x6 array containing its values in parameter
ranges [0,1] x [0,1] x [0,1].
Approximate the 3-d integral
∫︁ 1 ∫︁ 1 ∫︁ 1
(𝑎𝑏 − 𝑐)𝑑𝑎 𝑑𝑏 𝑑𝑐
0 0 0
over this volume with the mean. The exact result is: ln 2 − 21 ≈ 0.1931 . . . — what is your relative error?
(Hints: use elementwise operations and broadcasting. You can make np.ogrid give a number of points
in given range with np.ogrid[0:1:20j].)
Reminder Python functions:
Write a script that computes the Mandelbrot fractal. The Mandelbrot iteration:
N_max = 50
some_threshold = 50
c = x + 1j*y
z = 0
for j in range(N_max):
z = z**2 + c
• Computes the stationary distribution: the eigenvector of P.T with eigenvalue 1 (numerically: closest
to 1) => p_stationary
Remember to normalize the eigenvector — I didn’t. . .
• Checks if p_50 and p_stationary are equal to tolerance 1e-5
Toolbox: np.random, @, np.linalg.eig, reductions, abs(), argmin, comparisons, all, np.linalg.norm,
etc.
Solution: Python source file
1D plotting
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0, 3, 20)
y = np.linspace(0, 9, 20)
plt.plot(x, y)
(continues on next page)
2D plotting
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng()
image = rng.random((30, 30))
plt.imshow(image, cmap=plt.cm.hot)
plt.colorbar()
plt.show()
Distances exercise
import numpy as np
import matplotlib.pyplot as plt
Fitting to polynomial
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(27446968)
x = np.linspace(0, 1, 20)
y = np.cos(x) + 0.3 * rng.random(20)
p = np.poly1d(np.polyfit(x, y, 3))
t = np.linspace(0, 1, 200)
plt.plot(x, y, "o", t, p(t), "-")
plt.show()
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(27446968)
x = np.linspace(-1, 1, 2000)
y = np.cos(x) + 0.3 * rng.random(2000)
p = np.polynomial.Chebyshev.fit(x, y, 90)
plt.plot(x, y, "r.")
plt.plot(x, p(x), "k-", lw=3)
plt.show()
Population exercise
import numpy as np
import matplotlib.pyplot as plt
data = np.loadtxt("../../../data/populations.txt")
year, hares, lynxes, carrots = data.T
import numpy as np
import matplotlib.pyplot as plt
original figure
plt.figure()
img = plt.imread("../../../data/elephant.png")
plt.imshow(img)
plt.figure()
img_red = img[:, :, 0]
plt.imshow(img_red, cmap=plt.cm.gray)
lower resolution
plt.figure()
img_tiny = img[::6, ::6]
plt.imshow(img_tiny, interpolation="nearest")
plt.show()
Mandelbrot set
import numpy as np
import matplotlib.pyplot as plt
from numpy import newaxis
import warnings
# Mandelbrot iteration
z = c
# The code below overflows in many regions of the x-y grid, suppress
# warnings temporarily
with warnings.catch_warnings():
warnings.simplefilter("ignore")
for j in range(N_max):
z = z**2 + c
mandelbrot_set = abs(z) < some_threshold
return mandelbrot_set
Plot distance as a function of time for a random walk together with the theoretical result
import numpy as np
import matplotlib.pyplot as plt
t = np.arange(t_max)
# Steps can be -1 or 1 (note that randint excludes the upper limit)
rng = np.random.default_rng()
steps = 2 * rng.integers(0, 1 + 1, (n_stories, t_max)) - 1
Thanks
Many thanks to Bill Wing and Christoph Deil for review and corrections.
Chapter contents
• Introduction
• Simple plot
• Figures, Subplots, Axes and Ticks
• Other Types of Plots: examples and exercises
• Beyond this tutorial
• Quick references
• Full code examples
106
Scientific Python Lectures, Edition 2024.2rc0.dev0
4.1 Introduction
Tip: Matplotlib is probably the most used Python package for 2D-graphics. It provides both a quick
way to visualize data from Python and publication-quality figures in many formats. We are going to
explore matplotlib in interactive mode covering most common cases.
Tip: The Jupyter notebook and the IPython enhanced interactive Python, are tuned for the scientific-
computing workflow in Python, in combination with Matplotlib:
In [1]: %matplotlib
Jupyter notebook
In the notebook, insert, at the beginning of the notebook the following magic:
%matplotlib inline
4.1.2 pyplot
Tip: pyplot provides a procedural interface to the matplotlib object-oriented plotting library. It is
modeled closely after Matlab™. Therefore, the majority of plotting commands in pyplot have Matlab™
analogs with similar arguments. Important commands are explained with interactive examples.
Tip: In this section, we want to draw the cosine and sine functions on the same plot. Starting from
the default settings, we’ll enrich the figure step by step to make it nicer.
First step is to get the data for the sine and cosine functions:
import numpy as np
X is now a numpy array with 256 values ranging from −𝜋 to +𝜋 (included). C is the cosine (256 values)
and S is the sine (256 values).
To run the example, you can type them in an IPython interactive session:
$ ipython --matplotlib
Tip: You can also download each of the examples and run it using regular python, but you will lose
interactive data manipulation:
$ python plot_exercise_1.py
You can get source for each step by clicking on the corresponding figure.
Hint: Documentation
• plot tutorial
• plot() command
Tip: Matplotlib comes with a set of default settings that allow customizing all kinds of properties. You
can control the defaults of almost every property in matplotlib: figure size and dpi, line width, color and
style, axes, axis and grid properties, text and font properties and so on.
import numpy as np
import matplotlib.pyplot as plt
plt.plot(X, C)
plt.plot(X, S)
plt.show()
Hint: Documentation
• Customizing matplotlib
In the script below, we’ve instantiated (and commented) all the figure settings that influence the ap-
pearance of the plot.
Tip: The settings have been explicitly set to their default values, but now you can interactively play
with the values to explore their affect (see Line properties and Line styles below).
import numpy as np
import matplotlib.pyplot as plt
# Set x limits
plt.xlim(-4.0, 4.0)
# Set x ticks
plt.xticks(np.linspace(-4, 4, 9))
# Set y limits
plt.ylim(-1.0, 1.0)
# Set y ticks
(continues on next page)
Hint: Documentation
• Controlling line properties
• Line2D API
Tip: First step, we want to have the cosine in blue and the sine in red and a slightly thicker line for
both of them. We’ll also slightly alter the figure size to make it more horizontal.
...
plt.figure(figsize=(10, 6), dpi=80)
plt.plot(X, C, color="blue", linewidth=2.5, linestyle="-")
plt.plot(X, S, color="red", linewidth=2.5, linestyle="-")
...
Hint: Documentation
• xlim() command
• ylim() command
Tip: Current limits of the figure are a bit too tight and we want to make some space in order to clearly
see all data points.
...
plt.xlim(X.min() * 1.1, X.max() * 1.1)
plt.ylim(C.min() * 1.1, C.max() * 1.1)
...
Hint: Documentation
• xticks() command
• yticks() command
• Tick container
• Tick locating and formatting
Tip: Current ticks are not ideal because they do not show the interesting values (±𝜋,:math:pm pi/2)
for sine and cosine. We’ll change them such that they show only these values.
...
plt.xticks([-np.pi, -np.pi/2, 0, np.pi/2, np.pi])
plt.yticks([-1, 0, +1])
...
Hint: Documentation
• Working with text
• xticks() command
• yticks() command
• set_xticklabels()
• set_yticklabels()
Tip: Ticks are now properly placed but their label is not very explicit. We could guess that 3.142 is 𝜋
but it would be better to make it explicit. When we set tick values, we can also provide a corresponding
label in the second argument list. Note that we’ll use latex to allow for nice rendering of the label.
...
plt.xticks([-np.pi, -np.pi/2, 0, np.pi/2, np.pi],
[r'$-\pi$', r'$-\pi/2$', r'$0$', r'$+\pi/2$', r'$+\pi$'])
plt.yticks([-1, 0, +1],
[r'$-1$', r'$0$', r'$+1$'])
...
Hint: Documentation
• spines API
• Axis container
• Transformations tutorial
Tip: Spines are the lines connecting the axis tick marks and noting the boundaries of the data area.
They can be placed at arbitrary positions and until now, they were on the border of the axis. We’ll change
that since we want to have them in the middle. Since there are four of them (top/bottom/left/right),
we’ll discard the top and right by setting their color to none and we’ll move the bottom and left ones to
coordinate 0 in data space coordinates.
...
ax = plt.gca() # gca stands for 'get current axis'
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
ax.spines['bottom'].set_position(('data',0))
ax.yaxis.set_ticks_position('left')
ax.spines['left'].set_position(('data',0))
...
Hint: Documentation
• Legend guide
• legend() command
• legend API
Tip: Let’s add a legend in the upper left corner. This only requires adding the keyword argument label
(that will be used in the legend box) to the plot commands.
...
plt.plot(X, C, color="blue", linewidth=2.5, linestyle="-", label="cosine")
plt.plot(X, S, color="red", linewidth=2.5, linestyle="-", label="sine")
Hint: Documentation
• Annotating axis
• annotate() command
Tip: Let’s annotate some interesting points using the annotate command. We chose the 2𝜋/3 value
and we want to annotate both the sine and the cosine. We’ll first draw a marker on the curve as well as
a straight dotted line. Then, we’ll use the annotate command to display some text with an arrow.
...
t = 2 * np.pi / 3
plt.plot([t, t], [0, np.cos(t)], color='blue', linewidth=2.5, linestyle="--")
plt.scatter([t, ], [np.cos(t), ], 50, color='blue')
Hint: Documentation
• artist API
• set_bbox() method
Tip: The tick labels are now hardly visible because of the blue and red lines. We can make them
bigger and we can also adjust their properties such that they’ll be rendered on a semi-transparent white
background. This will allow us to see both the data and the labels.
...
for label in ax.get_xticklabels() + ax.get_yticklabels():
label.set_fontsize(16)
label.set_bbox(dict(facecolor='white', edgecolor='None', alpha=0.65))
...
A “figure” in matplotlib means the whole window in the user interface. Within this figure there can be
“subplots”.
Tip: So far we have used implicit figure and axes creation. This is handy for fast plots. We can have
more control over the display using figure, subplot, and axes explicitly. While subplot positions the plots
in a regular grid, axes allows free placement within the figure. Both can be useful depending on your
intention. We’ve already worked with figures and subplots without explicitly calling them. When we
call plot, matplotlib calls gca() to get the current axes and gca in turn calls gcf() to get the current
figure. If there is none it calls figure() to make one, strictly speaking, to make a subplot(111). Let’s
look at the details.
4.3.1 Figures
Tip: A figure is the windows in the GUI that has “Figure #” as title. Figures are numbered starting
from 1 as opposed to the normal Python way starting from 0. This is clearly MATLAB-style. There are
several parameters that determine what the figure looks like:
Tip: The defaults can be specified in the resource file and will be used most of the time. Only the
number of the figure is frequently changed.
As with other objects, you can set figure properties also setp or with the set_something methods.
When you work with the GUI you can close a figure by clicking on the x in the upper right corner. But
you can close a figure programmatically by calling close. Depending on the argument it closes (1) the
current figure (no argument), (2) a specific figure (figure number or figure instance as argument), or (3)
all figures ("all" as argument).
4.3.2 Subplots
Tip: With subplot you can arrange plots in a regular grid. You need to specify the number of rows and
columns and the number of the plot. Note that the gridspec command is a more powerful alternative.
4.3.3 Axes
Axes are very similar to subplots but allow placement of plots at any location in the figure. So if we
want to put a smaller plot inside a bigger one we do so with axes.
4.3.4 Ticks
Well formatted ticks are an important part of publishing-ready figures. Matplotlib provides a totally
configurable system for ticks. There are tick locators to specify where ticks should appear and tick
formatters to give ticks the appearance you want. Major and minor ticks can be located and formatted
independently from each other. Per default minor ticks are not shown, i.e. there is only an empty list
for them because it is as NullLocator (see below).
Tick Locators
Tick locators control the positions of the ticks. They are set as follows:
ax = plt.gca()
ax.xaxis.set_major_locator(eval(locator))
All of these locators derive from the base class matplotlib.ticker.Locator. You can make your own
locator deriving from it. Handling dates as ticks can be especially tricky. Therefore, matplotlib provides
special locators in matplotlib.dates.
Starting from the code below, try to reproduce the graphic taking care of filled areas:
n = 256
X = np.linspace(-np.pi, np.pi, n)
(continues on next page)
Starting from the code below, try to reproduce the graphic taking care of marker size, color and trans-
parency.
n = 1024
rng = np.random.default_rng()
X = rng.normal(0,1,n)
Y = rng.normal(0,1,n)
plt.scatter(X,Y)
Starting from the code below, try to reproduce the graphic by adding labels for red bars.
n = 12
X = np.arange(n)
rng = np.random.default_rng()
Y1 = (1 - X / float(n)) * rng.uniform(0.5, 1.0, n)
Y2 = (1 - X / float(n)) * rng.uniform(0.5, 1.0, n)
plt.ylim(-1.25, +1.25)
Starting from the code below, try to reproduce the graphic taking care of the colormap (see Colormaps
below).
n = 256
x = np.linspace(-3, 3, n)
y = np.linspace(-3, 3, n)
X, Y = np.meshgrid(x, y)
4.4.5 Imshow
Starting from the code below, try to reproduce the graphic taking care of colormap, image interpolation
and origin.
Hint: You need to take care of the origin of the image in the imshow command and use a colorbar()
n = 10
x = np.linspace(-3, 3, 4 * n)
y = np.linspace(-3, 3, 3 * n)
X, Y = np.meshgrid(x, y)
plt.imshow(f(X, Y))
Starting from the code below, try to reproduce the graphic taking care of colors and slices size.
rng = np.random.default_rng()
Z = rng.uniform(0, 1, 20)
plt.pie(Z)
Starting from the code below, try to reproduce the graphic taking care of colors and orientations.
n = 8
X, Y = np.mgrid[0:n, 0:n]
plt.quiver(X, Y)
4.4.8 Grids
Starting from the code below, try to reproduce the graphic taking care of line styles.
axes = plt.gca()
axes.set_xlim(0, 4)
axes.set_ylim(0, 3)
axes.set_xticklabels([])
axes.set_yticklabels([])
plt.subplot(2, 2, 1)
plt.subplot(2, 2, 3)
plt.subplot(2, 2, 4)
plt.axes([0, 0, 1, 1])
N = 20
theta = np.arange(0., 2 * np.pi, 2 * np.pi / N)
rng = np.random.default_rng()
radii = 10 * rng.random(N)
width = np.pi / 4 * rng.random(N)
bars = plt.bar(theta, radii, width=width, bottom=0.0)
4.4.11 3D Plots
fig = plt.figure()
ax = Axes3D(fig)
X = np.arange(-4, 4, 0.25)
Y = np.arange(-4, 4, 0.25)
X, Y = np.meshgrid(X, Y)
R = np.sqrt(X**2 + Y**2)
Z = np.sin(R)
4.4.12 Text
Quick read
If you want to do a first quick pass through the Scientific Python Lectures to learn the ecosystem, you
can directly skip to the next chapter: SciPy : high-level scientific computing.
The remainder of this chapter is not necessary to follow the rest of the intro part. But be sure to
come back and finish this chapter later.
Matplotlib benefits from extensive documentation as well as a large community of users and developers.
Here are some links of interest:
4.5.1 Tutorials
The code is well documented and you can quickly access a specific command from within a python
session:
Call signatures::
4.5.4 Galleries
The matplotlib gallery is also incredibly useful when you search how to render a given graphic. Each
example comes with its source.
Finally, there is a user mailing list where you can ask for help and a developers mailing list that is more
technical.
4.6.3 Markers
4.6.4 Colormaps
All colormaps can be reversed by appending _r. For instance, gray_r is the reverse of gray.
If you want to know more about colormaps, check the documentation on Colormaps in matplotlib.
The examples here are only examples relevant to the points raised in this chapter. The matplotlib
documentation comes with a much more exhaustive gallery.
Pie chart
import numpy as np
import matplotlib.pyplot as plt
n = 20
Z = np.ones(n)
Z[-1] *= 2
plt.show()
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
plt.show()
import numpy as np
import matplotlib.pyplot as plt
n = 1024
rng = np.random.default_rng()
X = rng.normal(0, 1, n)
Y = rng.normal(0, 1, n)
T = np.arctan2(Y, X)
plt.xlim(-1.5, 1.5)
plt.xticks([])
plt.ylim(-1.5, 1.5)
plt.yticks([])
plt.show()
Subplots
fig = plt.figure()
fig.subplots_adjust(bottom=0.025, left=0.025, top=0.975, right=0.975)
plt.subplot(2, 1, 1)
plt.xticks([]), plt.yticks([])
plt.subplot(2, 3, 4)
plt.xticks([])
plt.yticks([])
plt.subplot(2, 3, 5)
plt.xticks([])
plt.yticks([])
plt.subplot(2, 3, 6)
plt.xticks([])
plt.yticks([])
plt.show()
plt.figure(figsize=(6, 4))
plt.subplot(2, 1, 1)
plt.xticks([])
plt.yticks([])
plt.text(0.5, 0.5, "subplot(2,1,1)", ha="center", va="center", size=24, alpha=0.5)
plt.subplot(2, 1, 2)
plt.xticks([])
plt.yticks([])
plt.text(0.5, 0.5, "subplot(2,1,2)", ha="center", va="center", size=24, alpha=0.5)
plt.tight_layout()
plt.show()
plt.figure(figsize=(6, 4))
plt.subplot(1, 2, 1)
plt.xticks([])
plt.yticks([])
plt.text(0.5, 0.5, "subplot(1,2,1)", ha="center", va="center", size=24, alpha=0.5)
plt.subplot(1, 2, 2)
plt.xticks([])
plt.yticks([])
plt.text(0.5, 0.5, "subplot(1,2,2)", ha="center", va="center", size=24, alpha=0.5)
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
plt.show()
3D plotting
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
ax = plt.figure().add_subplot(projection="3d")
X = np.arange(-4, 4, 0.25)
Y = np.arange(-4, 4, 0.25)
X, Y = np.meshgrid(X, Y)
R = np.sqrt(X**2 + Y**2)
Z = np.sin(R)
plt.show()
Imshow elaborate
import numpy as np
import matplotlib.pyplot as plt
n = 10
x = np.linspace(-3, 3, int(3.5 * n))
y = np.linspace(-3, 3, int(3.0 * n))
X, Y = np.meshgrid(x, y)
Z = f(X, Y)
plt.xticks([])
plt.yticks([])
plt.show()
A simple example showing how to plot a vector field (quiver) with matplotlib.
import numpy as np
import matplotlib.pyplot as plt
n = 8
X, Y = np.mgrid[0:n, 0:n]
T = np.arctan2(Y - n / 2.0, X - n / 2.0)
R = 10 + np.sqrt((Y - n / 2.0) ** 2 + (X - n / 2.0) ** 2)
U, V = R * np.cos(T), R * np.sin(T)
plt.xlim(-1, n)
plt.xticks([])
plt.ylim(-1, n)
plt.yticks([])
plt.show()
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
plt.xticks([])
plt.yticks(np.arange(-1.0, 1.0, 0.2))
plt.grid()
ax = plt.gca()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
n = 256
x = np.linspace(-3, 3, n)
y = np.linspace(-3, 3, n)
X, Y = np.meshgrid(x, y)
plt.xticks([])
plt.yticks([])
plt.show()
import numpy as np
import matplotlib.pyplot as plt
N = 20
theta = np.arange(0.0, 2 * np.pi, 2 * np.pi / N)
rng = np.random.default_rng()
radii = 10 * rng.random(N)
width = np.pi / 4 * rng.random(N)
bars = plt.bar(theta, radii, width=width, bottom=0.0)
ax.set_xticklabels([])
ax.set_yticklabels([])
plt.show()
import numpy as np
import matplotlib.pyplot as plt
n = 256
X = np.linspace(-np.pi, np.pi, n)
Y = np.sin(2 * X)
plt.xlim(-np.pi, np.pi)
plt.xticks([])
plt.ylim(-2.5, 2.5)
plt.yticks([])
plt.show()
Subplot grid
plt.figure(figsize=(6, 4))
plt.subplot(2, 2, 1)
plt.xticks([])
plt.yticks([])
plt.text(0.5, 0.5, "subplot(2,2,1)", ha="center", va="center", size=20, alpha=0.5)
plt.subplot(2, 2, 2)
plt.xticks([])
plt.yticks([])
plt.text(0.5, 0.5, "subplot(2,2,2)", ha="center", va="center", size=20, alpha=0.5)
plt.subplot(2, 2, 3)
plt.xticks([])
plt.yticks([])
plt.subplot(2, 2, 4)
plt.xticks([])
plt.yticks([])
plt.text(0.5, 0.5, "subplot(2,2,4)", ha="center", va="center", size=20, alpha=0.5)
plt.tight_layout()
plt.show()
Bar plots
import numpy as np
import matplotlib.pyplot as plt
n = 12
X = np.arange(n)
rng = np.random.default_rng()
Y1 = (1 - X / float(n)) * rng.uniform(0.5, 1.0, n)
Y2 = (1 - X / float(n)) * rng.uniform(0.5, 1.0, n)
plt.xlim(-0.5, n)
plt.xticks([])
plt.ylim(-1.25, 1.25)
plt.yticks([])
plt.show()
Axes
plt.show()
Grid
ax.set_xlim(0, 4)
ax.set_ylim(0, 3)
ax.xaxis.set_major_locator(plt.MultipleLocator(1.0))
(continues on next page)
plt.show()
3D plotting
ax = plt.figure().add_subplot(projection="3d")
X, Y, Z = axes3d.get_test_data(0.05)
cset = ax.contourf(X, Y, Z)
ax.clabel(cset, fontsize=9, inline=1)
(continues on next page)
plt.xticks([])
plt.yticks([])
ax.set_zticks([])
ax.text2D(
-0.05,
1.05,
" 3D plots \n",
horizontalalignment="left",
verticalalignment="top",
bbox={"facecolor": "white", "alpha": 1.0},
family="DejaVu Sans",
size="x-large",
transform=plt.gca().transAxes,
)
ax.text2D(
-0.05,
0.975,
" Plot 2D or 3D data",
horizontalalignment="left",
verticalalignment="top",
family="DejaVu Sans",
size="medium",
transform=plt.gca().transAxes,
)
plt.show()
GridSpec
plt.figure(figsize=(6, 4))
G = gridspec.GridSpec(3, 3)
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
eqs = []
eqs.append(
r"$W^{3\beta}_{\delta_1 \rho_1 \sigma_2} = U^{3\beta}_{\delta_1 \rho_1} + \frac{1}
˓→ {8 \pi 2} \int^{\alpha_2}_{\alpha_2} d \alpha^\prime_2 \left[\frac{ U^{2\beta}_{\
˓→sigma_2}}\right]$"
)
eqs.append(
r"$\frac{d\rho}{d t} + \rho \vec{v} \cdot\nabla\vec{v} = -\nabla p + \mu\nabla^2 \
˓→vec{v} + \rho \vec{g} $"
)
eqs.append(r"$\int_{-\infty}^\infty e^{-x^2}dx=\sqrt{\pi}$")
eqs.append(r"$E = mc^2 = \sqrt{{m_0}^2c^4 + p^2c^2}$")
eqs.append(r"$F_G = G\frac{m_1m_2} {r^2}$")
rng = np.random.default_rng()
for i in range(24):
index = rng.integers(0, len(eqs))
eq = eqs[index]
size = np.random.uniform(12, 32)
x, y = np.random.uniform(0, 1, 2)
alpha = np.random.uniform(0.25, 0.75)
plt.text(
x,
y,
eq,
ha="center",
va="center",
color="#11557c",
alpha=alpha,
transform=plt.gca().transAxes,
fontsize=size,
clip_on=True,
)
plt.xticks([])
plt.yticks([])
plt.show()
Exercise 1
import numpy as np
import matplotlib.pyplot as plt
n = 256
X = np.linspace(-np.pi, np.pi, 256)
C, S = np.cos(X), np.sin(X)
plt.plot(X, C)
plt.plot(X, S)
plt.show()
Exercise 4
import numpy as np
import matplotlib.pyplot as plt
plt.show()
Exercise 3
import numpy as np
import matplotlib.pyplot as plt
plt.xlim(-4.0, 4.0)
plt.xticks(np.linspace(-4, 4, 9))
plt.ylim(-1.0, 1.0)
plt.yticks(np.linspace(-1, 1, 5))
plt.show()
Exercise 5
import numpy as np
import matplotlib.pyplot as plt
plt.show()
Exercise 6
import numpy as np
import matplotlib.pyplot as plt
plt.show()
Exercise 2
import numpy as np
import matplotlib.pyplot as plt
# Create a new figure of size 8x6 points, using 100 dots per inch
plt.figure(figsize=(8, 6), dpi=80)
# Plot cosine using blue color with a continuous line of width 1 (pixels)
plt.plot(X, C, color="blue", linewidth=1.0, linestyle="-")
# Plot sine using green color with a continuous line of width 1 (pixels)
plt.plot(X, S, color="green", linewidth=1.0, linestyle="-")
# Set x limits
plt.xlim(-4.0, 4.0)
# Set x ticks
plt.xticks(np.linspace(-4, 4, 9))
# Set y ticks
plt.yticks(np.linspace(-1, 1, 5))
Exercise 7
import numpy as np
import matplotlib.pyplot as plt
ax = plt.gca()
(continues on next page)
plt.show()
Exercise 8
import numpy as np
import matplotlib.pyplot as plt
ax = plt.gca()
ax.spines["right"].set_color("none")
ax.spines["top"].set_color("none")
ax.xaxis.set_ticks_position("bottom")
ax.spines["bottom"].set_position(("data", 0))
ax.yaxis.set_ticks_position("left")
ax.spines["left"].set_position(("data", 0))
plt.legend(loc="upper left")
plt.show()
Exercise 9
import numpy as np
import matplotlib.pyplot as plt
ax = plt.gca()
ax.spines["right"].set_color("none")
ax.spines["top"].set_color("none")
ax.xaxis.set_ticks_position("bottom")
ax.spines["bottom"].set_position(("data", 0))
ax.yaxis.set_ticks_position("left")
ax.spines["left"].set_position(("data", 0))
t = 2 * np.pi / 3
plt.plot([t, t], [0, np.cos(t)], color="blue", linewidth=1.5, linestyle="--")
plt.scatter(
(continues on next page)
plt.legend(loc="upper left")
plt.show()
Exercise
import numpy as np
import matplotlib.pyplot as plt
ax = plt.gca()
ax.spines["right"].set_color("none")
ax.spines["top"].set_color("none")
ax.xaxis.set_ticks_position("bottom")
ax.spines["bottom"].set_position(("data", 0))
ax.yaxis.set_ticks_position("left")
ax.spines["left"].set_position(("data", 0))
t = 2 * np.pi / 3
plt.plot([t, t], [0, np.cos(t)], color="blue", linewidth=1.5, linestyle="--")
plt.scatter(
[
t,
],
[
np.cos(t),
],
50,
color="blue",
)
plt.annotate(
r"$sin(\frac{2\pi}{3} )=\frac{\sqrt{3} }{2} $",
xy=(t, np.sin(t)),
xycoords="data",
xytext=(10, 30),
textcoords="offset points",
fontsize=16,
arrowprops={"arrowstyle": "->", "connectionstyle": "arc3,rad=.2"},
)
plt.show()
size = 256, 16
dpi = 72.0
figsize = size[0] / float(dpi), size[1] / float(dpi)
fig = plt.figure(figsize=figsize, dpi=dpi)
fig.patch.set_alpha(0)
plt.axes([0, 0.1, 1, 0.8], frameon=False)
plt.xlim(0, 11)
plt.xticks([])
plt.yticks([])
plt.show()
Linewidth
size = 256, 16
dpi = 72.0
figsize = size[0] / float(dpi), size[1] / float(dpi)
fig = plt.figure(figsize=figsize, dpi=dpi)
fig.patch.set_alpha(0)
plt.axes([0, 0.1, 1, 0.8], frameon=False)
plt.xlim(0, 11)
plt.ylim(0, 1)
plt.xticks([])
plt.yticks([])
plt.show()
Alpha: transparency
size = 256, 16
dpi = 72.0
figsize = size[0] / float(dpi), size[1] / float(dpi)
fig = plt.figure(figsize=figsize, dpi=dpi)
fig.patch.set_alpha(0)
plt.axes([0, 0.1, 1, 0.8], frameon=False)
plt.xlim(0, 11)
plt.xticks([])
plt.yticks([])
plt.show()
size = 128, 16
dpi = 72.0
figsize = size[0] / float(dpi), size[1] / float(dpi)
fig = plt.figure(figsize=figsize, dpi=dpi)
fig.patch.set_alpha(0)
plt.rcParams["text.antialiased"] = False
plt.text(0.5, 0.5, "Aliased", ha="center", va="center")
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.xticks([])
plt.yticks([])
plt.show()
size = 128, 16
dpi = 72.0
figsize = size[0] / float(dpi), size[1] / float(dpi)
fig = plt.figure(figsize=figsize, dpi=dpi)
fig.patch.set_alpha(0)
plt.axes([0, 0, 1, 1], frameon=False)
plt.rcParams["text.antialiased"] = True
plt.text(0.5, 0.5, "Anti-aliased", ha="center", va="center")
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.xticks([])
plt.yticks([])
plt.show()
Marker size
size = 256, 16
dpi = 72.0
figsize = size[0] / float(dpi), size[1] / float(dpi)
fig = plt.figure(figsize=figsize, dpi=dpi)
fig.patch.set_alpha(0)
plt.axes([0, 0, 1, 1], frameon=False)
plt.xlim(0, 11)
plt.xticks([])
plt.yticks([])
plt.show()
size = 256, 16
dpi = 72.0
figsize = size[0] / float(dpi), size[1] / float(dpi)
fig = plt.figure(figsize=figsize, dpi=dpi)
fig.patch.set_alpha(0)
plt.axes([0, 0, 1, 1], frameon=False)
plt.show()
Colormaps
import numpy as np
import matplotlib.pyplot as plt
plt.rc("text", usetex=False)
a = np.outer(np.arange(0, 1, 0.01), np.ones(10))
plt.figure(figsize=(10, 5))
plt.subplots_adjust(top=0.8, bottom=0.05, left=0.01, right=0.99)
maps = [m for m in plt.cm.datad if not m.endswith("_r")]
maps.sort()
l = len(maps) + 1
for i, m in enumerate(maps):
plt.subplot(1, l, i + 1)
plt.axis("off")
plt.imshow(a, aspect="auto", cmap=plt.get_cmap(m), origin="lower")
plt.title(m, rotation=90, fontsize=10, va="bottom")
plt.show()
import numpy as np
import matplotlib.pyplot as plt
size = 256, 16
dpi = 72.0
figsize = size[0] / float(dpi), size[1] / float(dpi)
fig = plt.figure(figsize=figsize, dpi=dpi)
fig.patch.set_alpha(0)
plt.axes([0, 0, 1, 1], frameon=False)
plt.xlim(0, 12)
plt.ylim(-1, 2)
plt.xticks([])
plt.yticks([])
plt.show()
import numpy as np
import matplotlib.pyplot as plt
size = 256, 16
dpi = 72.0
figsize = size[0] / float(dpi), size[1] / float(dpi)
fig = plt.figure(figsize=figsize, dpi=dpi)
fig.patch.set_alpha(0)
plt.axes([0, 0, 1, 1], frameon=False)
plt.plot(
5 + np.arange(4), np.ones(4), color="blue", linewidth=8, solid_capstyle="round"
)
plt.xlim(0, 14)
plt.xticks([])
plt.yticks([])
plt.show()
import numpy as np
import matplotlib.pyplot as plt
size = 256, 16
dpi = 72.0
figsize = size[0] / float(dpi), size[1] / float(dpi)
fig = plt.figure(figsize=figsize, dpi=dpi)
fig.patch.set_alpha(0)
plt.axes([0, 0, 1, 1], frameon=False)
rng = np.random.default_rng()
plt.xlim(0, 11)
plt.xticks([])
plt.yticks([])
(continues on next page)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
size = 256, 16
dpi = 72.0
figsize = size[0] / float(dpi), size[1] / float(dpi)
fig = plt.figure(figsize=figsize, dpi=dpi)
fig.patch.set_alpha(0)
plt.axes([0, 0, 1, 1], frameon=False)
rng = np.random.default_rng()
Dash capstyle
import numpy as np
import matplotlib.pyplot as plt
size = 256, 16
dpi = 72.0
figsize = size[0] / float(dpi), size[1] / float(dpi)
fig = plt.figure(figsize=figsize, dpi=dpi)
fig.patch.set_alpha(0)
plt.axes([0, 0, 1, 1], frameon=False)
plt.plot(
np.arange(4),
np.ones(4),
color="blue",
dashes=[15, 15],
linewidth=8,
dash_capstyle="butt",
)
plt.plot(
5 + np.arange(4),
np.ones(4),
color="blue",
dashes=[15, 15],
linewidth=8,
dash_capstyle="round",
)
plt.plot(
10 + np.arange(4),
np.ones(4),
color="blue",
dashes=[15, 15],
linewidth=8,
dash_capstyle="projecting",
)
plt.xlim(0, 14)
plt.xticks([])
plt.yticks([])
plt.show()
import numpy as np
import matplotlib.pyplot as plt
size = 256, 16
dpi = 72.0
figsize = size[0] / float(dpi), size[1] / float(dpi)
fig = plt.figure(figsize=figsize, dpi=dpi)
fig.patch.set_alpha(0)
plt.axes([0, 0, 1, 1], frameon=False)
plt.plot(
np.arange(3),
[0, 1, 0],
color="blue",
dashes=[12, 5],
linewidth=8,
dash_joinstyle="miter",
)
plt.plot(
4 + np.arange(3),
[0, 1, 0],
color="blue",
dashes=[12, 5],
linewidth=8,
dash_joinstyle="bevel",
)
plt.plot(
8 + np.arange(3),
[0, 1, 0],
color="blue",
dashes=[12, 5],
linewidth=8,
dash_joinstyle="round",
)
plt.xlim(0, 12)
plt.ylim(-1, 2)
plt.xticks([])
plt.yticks([])
plt.show()
Markers
import numpy as np
import matplotlib.pyplot as plt
markers = [
0,
1,
2,
3,
4,
5,
6,
7,
"o",
"h",
"_",
"1",
"2",
"3",
"4",
"8",
"p",
"^",
"v",
"<",
">",
"|",
(continues on next page)
n_markers = len(markers)
for i, m in enumerate(markers):
marker(m, i)
plt.show()
Linestyles
import numpy as np
import matplotlib.pyplot as plt
linestyles = [
"-",
"--",
":",
"-.",
".",
",",
"o",
"^",
"v",
"<",
">",
"s",
"+",
(continues on next page)
for i, ls in enumerate(linestyles):
linestyle(ls, i)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
def tickline():
plt.xlim(0, 10), plt.ylim(-1, 1), plt.yticks([])
ax = plt.gca()
ax.spines["right"].set_color("none")
ax.spines["left"].set_color("none")
ax.spines["top"].set_color("none")
ax.xaxis.set_ticks_position("bottom")
ax.spines["bottom"].set_position(("data", 0))
ax.yaxis.set_ticks_position("none")
ax.xaxis.set_minor_locator(plt.MultipleLocator(0.1))
ax.plot(np.arange(11), np.zeros(11))
return ax
locators = [
"plt.NullLocator()",
"plt.MultipleLocator(1.0)",
"plt.FixedLocator([0, 2, 8, 9, 10])",
"plt.IndexLocator(3, 1)",
"plt.LinearLocator(5)",
"plt.LogLocator(2, [1.0])",
"plt.AutoLocator()",
]
n_locators = len(locators)
3D plotting vignette
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
(continues on next page)
ax.text2D(
0.05,
0.93,
" 3D plots \n",
horizontalalignment="left",
verticalalignment="top",
size="xx-large",
bbox={"facecolor": "white", "alpha": 1.0},
transform=plt.gca().transAxes,
)
ax.text2D(
0.05,
0.87,
" Plot 2D or 3D data",
horizontalalignment="left",
verticalalignment="top",
size="large",
transform=plt.gca().transAxes,
)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
plt.subplot(1, 1, 1, polar=True)
N = 20
theta = np.arange(0.0, 2 * np.pi, 2 * np.pi / N)
rng = np.random.default_rng()
radii = 10 * rng.random(N)
width = np.pi / 4 * rng.random(N)
bars = plt.bar(theta, radii, width=width, bottom=0.0)
for r, bar in zip(radii, bars, strict=True):
bar.set_facecolor(plt.cm.jet(r / 10.0))
bar.set_alpha(0.5)
plt.gca().set_xticklabels([])
plt.gca().set_yticklabels([])
plt.text(
-0.2,
1.02,
" Polar Axis \n",
horizontalalignment="left",
verticalalignment="top",
size="xx-large",
bbox={"facecolor": "white", "alpha": 1.0},
transform=plt.gca().transAxes,
)
(continues on next page)
plt.text(
-0.2,
1.01,
"\n\n Plot anything using polar axis ",
horizontalalignment="left",
verticalalignment="top",
size="large",
transform=plt.gca().transAxes,
)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
n = 256
X = np.linspace(0, 2, n)
Y = np.sin(2 * np.pi * X)
(continues on next page)
ax = plt.gca()
ax.add_patch(
FancyBboxPatch(
(-0.05, 0.87),
width=0.66,
height=0.165,
clip_on=False,
boxstyle="square,pad=0",
zorder=3,
facecolor="white",
alpha=1.0,
transform=plt.gca().transAxes,
)
)
plt.text(
-0.05,
1.02,
" Regular Plot: plt.plot(...)\n",
horizontalalignment="left",
verticalalignment="top",
size="xx-large",
transform=plt.gca().transAxes,
)
plt.text(
-0.05,
1.01,
"\n\n Plot lines and/or markers ",
horizontalalignment="left",
verticalalignment="top",
size="large",
transform=plt.gca().transAxes,
)
plt.show()
ax = plt.subplot(2, 1, 1)
ax.set_xticklabels([])
ax.set_yticklabels([])
ax = plt.gca()
ax.add_patch(
FancyBboxPatch(
(-0.05, 0.72),
width=0.66,
height=0.34,
clip_on=False,
boxstyle="square,pad=0",
zorder=3,
facecolor="white",
alpha=1.0,
transform=plt.gca().transAxes,
)
)
(continues on next page)
plt.text(
-0.05,
1.02,
" Multiplot: plt.subplot(...)\n",
horizontalalignment="left",
verticalalignment="top",
size="xx-large",
transform=ax.transAxes,
)
plt.text(
-0.05,
1.01,
"\n\n Plot several plots at once ",
horizontalalignment="left",
verticalalignment="top",
size="large",
transform=ax.transAxes,
)
ax = plt.subplot(2, 2, 3)
ax.set_xticklabels([])
ax.set_yticklabels([])
ax = plt.subplot(2, 2, 4)
ax.set_xticklabels([])
ax.set_yticklabels([])
plt.show()
import numpy as np
import matplotlib.pyplot as plt
n = 5
Z = np.zeros((n, 4))
X = np.linspace(0, 2, n)
rng = np.random.default_rng()
Y = rng.random((n, 4))
plt.boxplot(Y)
plt.xticks([])
plt.yticks([])
ax = plt.gca()
ax.add_patch(
FancyBboxPatch(
(-0.05, 0.87),
width=0.66,
height=0.165,
clip_on=False,
boxstyle="square,pad=0",
zorder=3,
facecolor="white",
alpha=1.0,
transform=plt.gca().transAxes,
(continues on next page)
plt.text(
-0.05,
1.02,
" Box Plot: plt.boxplot(...)\n ",
horizontalalignment="left",
verticalalignment="top",
size="xx-large",
transform=axes.transAxes,
)
plt.text(
-0.04,
0.98,
"\n Make a box and whisker plot ",
horizontalalignment="left",
verticalalignment="top",
size="large",
transform=axes.transAxes,
)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
n = 1024
rng = np.random.default_rng()
X = rng.normal(0, 1, n)
Y = rng.normal(0, 1, n)
T = np.arctan2(Y, X)
ax = plt.gca()
ax.add_patch(
FancyBboxPatch(
(-0.05, 0.87),
width=0.66,
height=0.165,
clip_on=False,
boxstyle="square,pad=0",
(continues on next page)
plt.text(
-0.05,
1.02,
" Scatter Plot: plt.scatter(...)\n",
horizontalalignment="left",
verticalalignment="top",
size="xx-large",
transform=plt.gca().transAxes,
)
plt.text(
-0.05,
1.01,
"\n\n Make a scatter plot of x versus y ",
horizontalalignment="left",
verticalalignment="top",
size="large",
transform=plt.gca().transAxes,
)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
n = 20
X = np.ones(n)
X[-1] *= 2
plt.pie(X, explode=X * 0.05, colors=[f"{ i / float(n): f} " for i in range(n)])
fig = plt.gcf()
w, h = fig.get_figwidth(), fig.get_figheight()
r = h / float(w)
plt.xlim(-1.5, 1.5)
plt.ylim(-1.5 * r, 1.5 * r)
plt.xticks([])
plt.yticks([])
ax = plt.gca()
ax.add_patch(
FancyBboxPatch(
(-0.05, 0.87),
width=0.66,
height=0.165,
clip_on=False,
(continues on next page)
plt.text(
-0.05,
1.02,
" Pie Chart: plt.pie(...)\n",
horizontalalignment="left",
verticalalignment="top",
size="xx-large",
transform=plt.gca().transAxes,
)
plt.text(
-0.05,
1.01,
"\n\n Make a pie chart of an array ",
horizontalalignment="left",
verticalalignment="top",
size="large",
transform=plt.gca().transAxes,
)
plt.show()
Imshow demo
Demoing imshow
import numpy as np
import matplotlib.pyplot as plt
n = 10
x = np.linspace(-3, 3, 8 * n)
y = np.linspace(-3, 3, 6 * n)
X, Y = np.meshgrid(x, y)
Z = f(X, Y)
plt.imshow(Z, interpolation="nearest", cmap="bone", origin="lower")
plt.xticks([])
plt.yticks([])
ax = plt.gca()
ax.add_patch(
FancyBboxPatch(
(-0.05, 0.87),
width=0.66,
height=0.165,
clip_on=False,
(continues on next page)
plt.text(
-0.05,
1.02,
" Imshow: plt.imshow(...)\n",
horizontalalignment="left",
verticalalignment="top",
size="xx-large",
transform=plt.gca().transAxes,
)
plt.text(
-0.05,
1.01,
"\n\n Display an image to current axes ",
horizontalalignment="left",
verticalalignment="top",
family="DejaVu Sans",
size="large",
transform=plt.gca().transAxes,
)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
n = 16
X = np.arange(n)
Y1 = (1 - X / float(n)) * np.random.uniform(0.5, 1.0, n)
Y2 = (1 - X / float(n)) * np.random.uniform(0.5, 1.0, n)
plt.bar(X, Y1, facecolor="#9999ff", edgecolor="white")
plt.bar(X, -Y2, facecolor="#ff9999", edgecolor="white")
plt.xlim(-0.5, n)
plt.xticks([])
plt.ylim(-1, 1)
plt.yticks([])
ax = plt.gca()
ax.add_patch(
FancyBboxPatch(
(-0.05, 0.87),
width=0.66,
height=0.165,
clip_on=False,
boxstyle="square,pad=0",
zorder=3,
facecolor="white",
(continues on next page)
plt.text(
-0.05,
1.02,
" Bar Plot: plt.bar(...)\n",
horizontalalignment="left",
verticalalignment="top",
size="xx-large",
transform=plt.gca().transAxes,
)
plt.text(
-0.05,
1.01,
"\n\n Make a bar plot with rectangles ",
horizontalalignment="left",
verticalalignment="top",
size="large",
transform=plt.gca().transAxes,
)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
n = 8
X, Y = np.mgrid[0:n, 0:n]
T = np.arctan2(Y - n / 2.0, X - n / 2.0)
R = 10 + np.sqrt((Y - n / 2.0) ** 2 + (X - n / 2.0) ** 2)
U, V = R * np.cos(T), R * np.sin(T)
plt.quiver(X, Y, U, V, R, alpha=0.5)
plt.quiver(X, Y, U, V, edgecolor="k", facecolor="None", linewidth=0.5)
plt.xlim(-1, n)
plt.xticks([])
plt.ylim(-1, n)
plt.yticks([])
ax = plt.gca()
ax.add_patch(
FancyBboxPatch(
(-0.05, 0.87),
width=0.66,
height=0.165,
clip_on=False,
(continues on next page)
plt.text(
-0.05,
1.02,
" Quiver Plot: plt.quiver(...)\n",
horizontalalignment="left",
verticalalignment="top",
size="xx-large",
transform=plt.gca().transAxes,
)
plt.text(
-0.05,
1.01,
"\n\n Plot a 2-D field of arrows ",
horizontalalignment="left",
verticalalignment="top",
size="large",
transform=plt.gca().transAxes,
)
plt.show()
An example demoing how to plot the contours of a function, with additional layout tweaks.
/home/runner/work/scientific-python-lectures/scientific-python-lectures/intro/
˓→matplotlib/examples/pretty_plots/plot_contour_ext.py:24: UserWarning: The following␣
import numpy as np
import matplotlib.pyplot as plt
n = 256
x = np.linspace(-3, 3, n)
y = np.linspace(-3, 3, n)
X, Y = np.meshgrid(x, y)
ax = plt.gca()
ax.add_patch(
FancyBboxPatch(
(-0.05, 0.87),
width=0.66,
height=0.165,
clip_on=False,
boxstyle="square,pad=0",
zorder=3,
facecolor="white",
alpha=1.0,
transform=plt.gca().transAxes,
)
)
plt.text(
-0.05,
1.02,
" Contour Plot: plt.contour(..)\n",
horizontalalignment="left",
verticalalignment="top",
size="xx-large",
transform=plt.gca().transAxes,
)
plt.text(
-0.05,
1.01,
"\n\n Draw contour lines and filled contours ",
horizontalalignment="left",
verticalalignment="top",
size="large",
transform=plt.gca().transAxes,
)
plt.show()
Grid elaborate
axes.xaxis.set_major_locator(MultipleLocator(1.0))
axes.xaxis.set_minor_locator(MultipleLocator(0.1))
axes.yaxis.set_major_locator(MultipleLocator(1.0))
axes.yaxis.set_minor_locator(MultipleLocator(0.1))
axes.grid(which="major", axis="x", linewidth=0.75, linestyle="-", color="0.75")
axes.grid(which="minor", axis="x", linewidth=0.25, linestyle="-", color="0.75")
axes.grid(which="major", axis="y", linewidth=0.75, linestyle="-", color="0.75")
axes.grid(which="minor", axis="y", linewidth=0.25, linestyle="-", color="0.75")
axes.set_xticklabels([])
(continues on next page)
ax = plt.gca()
ax.add_patch(
FancyBboxPatch(
(-0.05, 0.87),
width=0.66,
height=0.165,
clip_on=False,
boxstyle="square,pad=0",
zorder=3,
facecolor="white",
alpha=1.0,
transform=plt.gca().transAxes,
)
)
plt.text(
-0.05,
1.02,
" Grid: plt.grid(...)\n",
horizontalalignment="left",
verticalalignment="top",
size="xx-large",
transform=axes.transAxes,
)
plt.text(
-0.05,
1.01,
"\n\n Draw ticks and grid ",
horizontalalignment="left",
verticalalignment="top",
size="large",
transform=axes.transAxes,
)
import numpy as np
import matplotlib.pyplot as plt
fig = plt.figure()
plt.xticks([])
plt.yticks([])
eqs = []
eqs.append(
r"$W^{3\beta}_{\delta_1 \rho_1 \sigma_2} = U^{3\beta}_{\delta_1 \rho_1} + \frac{1}
˓→ {8 \pi 2} \int^{\alpha_2}_{\alpha_2} d \alpha^\prime_2 \left[\frac{ U^{2\beta}_{\
˓→sigma_2}}\right]$"
)
eqs.append(
r"$\frac{d\rho}{d t} + \rho \vec{v} \cdot\nabla\vec{v} = -\nabla p + \mu\nabla^2 \
˓→vec{v} + \rho \vec{g} $"
)
eqs.append(r"$\int_{-\infty}^\infty e^{-x^2}dx=\sqrt{\pi}$")
eqs.append(r"$E = mc^2 = \sqrt{{m_0}^2c^4 + p^2c^2}$")
eqs.append(r"$F_G = G\frac{m_1m_2} {r^2}$")
rng = np.random.default_rng()
ax = plt.gca()
ax.add_patch(
FancyBboxPatch(
(-0.05, 0.87),
width=0.66,
height=0.165,
clip_on=False,
boxstyle="square,pad=0",
zorder=3,
facecolor="white",
alpha=1.0,
transform=plt.gca().transAxes,
)
)
plt.text(
-0.05,
1.02,
" Text: plt.text(...)\n",
horizontalalignment="left",
verticalalignment="top",
size="xx-large",
transform=plt.gca().transAxes,
)
plt.text(
-0.05,
1.01,
"\n\n Draw any kind of text ",
horizontalalignment="left",
verticalalignment="top",
size="large",
transform=plt.gca().transAxes,
(continues on next page)
plt.show()
Authors: Gaël Varoquaux, Adrien Chauve, Andre Espaze, Emmanuelle Gouillart, Ralf Gommers
Scipy
The scipy package contains various toolboxes dedicated to common issues in scientific computing.
Its different submodules correspond to different applications, such as interpolation, integration, opti-
mization, image processing, statistics, special functions, etc.
Tip: scipy can be compared to other standard scientific-computing libraries, such as the GSL (GNU
Scientific Library for C and C++), or Matlab’s toolboxes. scipy is the core package for scientific routines
in Python; it is meant to operate efficiently on numpy arrays, so that NumPy and SciPy work hand in
hand.
Before implementing a routine, it is worth checking if the desired data processing is not already imple-
mented in SciPy. As non-professional programmers, scientists often tend to re-invent the wheel, which
leads to buggy, non-optimal, difficult-to-share and unmaintainable code. By contrast, SciPy’s routines
are optimized and tested, and should therefore be used when possible.
Chapters contents
210
Scientific Python Lectures, Edition 2024.2rc0.dev0
Warning: This tutorial is far from an introduction to numerical computing. As enumerating the
different submodules and functions in SciPy would be very boring, we concentrate instead on a few
examples to give a general idea of how to use scipy for scientific computing.
Tip: They all depend on numpy, but are mostly independent of each other. The standard way of
importing NumPy and these SciPy modules is:
scipy.io contains functions for loading and saving data in several common formats including Matlab,
IDL, Matrix Market, and Harwell-Boeing.
Matlab files: Loading and saving:
Warning: Python / Matlab mismatch: The Matlab file format does not support 1D arrays.
>>> a = np.ones(3)
>>> a
array([1., 1., 1.])
>>> a.shape
(3,)
>>> sp.io.savemat('file.mat', {'a': a})
>>> a2 = sp.io.loadmat('file.mat')['a']
>>> a2
array([[1., 1., 1.]])
>>> a2.shape
(1, 3)
Notice that the original array was a one-dimensional array, whereas the saved and reloaded array is
a two-dimensional array with a single row.
For other formats, see the scipy.io documentation.
See also:
• Load text files: numpy.loadtxt()/numpy.savetxt()
• Clever loading of text/csv files: numpy.genfromtxt()
• Fast and efficient, but NumPy-specific, binary format: numpy.save()/numpy.load()
• Basic input/output of images in Matplotlib: matplotlib.pyplot.imread()/matplotlib.pyplot.
imsave()
• More advanced input/output of images: imageio
“Special” functions are functions commonly used in science and mathematics that are not considered to
be “elementary” functions. Examples include
• the gamma function, scipy.special.gamma(),
• the error function, scipy.special.erf(),
• Bessel functions, such as scipy.special.jv() (Bessel function of the first kind), and
• elliptic functions, such as scipy.special.ellipj() (Jacobi elliptic functions).
Other special functions are combinations of familiar elementary functions, but they offer better accuracy
or robustness than their naive implementations would.
Most of these function are computed elementwise and follow standard NumPy broadcasting rules when
the input arrays have different shapes. For example, scipy.special.xlog1py() is mathematically
equivalent to 𝑥 log(1 + 𝑦).
>>> x = 2.5
>>> y = 1e-18
>>> x * np.log(1 + y)
0.0
>>> sp.special.xlog1py(x, y)
2.5e-18
Many special functions also have “logarithmized” variants. For instance, the gamma function Γ(·) is
related to the factorial function by 𝑛! = Γ(𝑛 + 1), but it extends the domain from the positive integers
to the complex plane.
>>> x = np.arange(10)
>>> np.allclose(sp.special.gamma(x + 1), sp.special.factorial(x))
True
>>> sp.special.gamma(5) < sp.special.gamma(5.5) < sp.special.gamma(6)
True
The factorial function grows quickly, and so the gamma function overflows for moderate values of the
argument. However, sometimes only the logarithm of the gamma function is needed. In such cases, we
can compute the logarithm of the gamma function directly using scipy.special.gammaln().
Such functions can often be used when the intermediate components of a calculation would overflow
or underflow, but the final result would not. For example, suppose we wish to compute the ratio
Γ(500)/Γ(499).
>>> a = sp.special.gamma(500)
>>> b = sp.special.gamma(499)
>>> a, b
(inf, inf)
Both the numerator and denominator overflow, so performing 𝑎/𝑏 will not return the result we seek.
However, the magnitude of the result should be moderate, so the use of logarithms comes to mind.
Combining the identities log(𝑎/𝑏) = log(𝑎) − log(𝑏) and exp(log(𝑥)) = 𝑥, we get:
Similarly, suppose we wish to compute the difference log(Γ(500) − Γ(499)). For this, we use scipy.
special.logsumexp(), which computes log(exp(𝑥)+exp(𝑦)) using a numerical trick that avoids overflow.
For more information about these and many other special functions, see the documentation of scipy.
special.
scipy.linalg provides a Python interface to efficient, compiled implementations of standard linear alge-
bra operations: the BLAS (Basic Linear Algebra Subroutines) and LAPACK (Linear Algebra PACKage)
libraries.
For example, the scipy.linalg.det() function computes the determinant of a square matrix:
Mathematically, the solution of a linear system 𝐴𝑥 = 𝑏 is 𝑥 = 𝐴−1 𝑏, but explicit inversion of a matrix is
numerically unstable and should be avoided. Instead, use scipy.linalg.solve():
Linear systems with special structure can often be solved more efficiently than more general systems. For
example, systems with triangular matrices can be solved using scipy.linalg.solve_triangular():
scipy.linalg also features matrix factorizations/decompositions such as the singular value decomposi-
tion.
Many other decompositions (e.g. LU, Cholesky, QR), solvers for structured linear systems (e.g. triangu-
lar, circulant), eigenvalue problem algorithms, matrix functions (e.g. matrix exponential), and routines
for special matrix creation (e.g. block diagonal, toeplitz) are available in scipy.linalg.
On the other hand, if the data are not noisy, it may be desirable to pass exactly through each point.
The derivative and antiderivative methods of the result object can be used for differentiation and
integration. For the latter, the constant of integration is assumed to be zero, but we can “wrap” the
antiderivative to include a nonzero constant of integration.
For functions that are monotonic on an interval (e.g. sin from 𝜋/2 to 3𝜋/2), we can reverse the arguments
of make_interp_spline to interpolate the inverse function. Because the first argument is expected to
be monotonically increasing, we also reverse the order of elements in the arrays with numpy.flip().
See the summary exercise on Maximum wind speed prediction at the Sprogø station for a more advanced
spline interpolation example, and read the SciPy interpolation tutorial and the scipy.interpolate
documentation for much more information.
scipy.optimize provides algorithms for root finding, curve fitting, and more general optimization.
Warning: None of the functions in scipy.optimize that accept a guess are guaranteed to converge
for all possible guesses! (For example, try x0=1.5 in the example above, where the derivative of
the function is exactly zero.) If this occurs, try a different guess, adjust the options (like providing
a bracket as shown below), or consider whether SciPy offers a more appropriate method for the
problem.
Note that only one the root at 1.0 is found. By inspection, we can tell that there is a second root at
2.0. We can direct the function toward a particular root by changing the guess or by passing a bracket
that contains only the root we seek.
Over-constrained problems can be solved in the least-squares sense using scipy.optimize.root() with
method='lm' (Levenberg-Marquardt).
We can approximate the underlying amplitude, frequency, and phase from the data by least squares
curve fitting. To begin, we write a function that accepts the independent variable as the first argument
and all parameters to fit as separate arguments:
The temperature extremes in Alaska for each month, starting in January, are given by (in
degrees Celsius):
max: 17, 19, 21, 28, 33, 38, 37, 37, 31, 23, 19, 18
min: -62, -59, -56, -46, -32, -18, -9, -13, -25, -46, -52, -58
5.5.3 Optimization
Suppose we wish to minimize the scalar-valued function of a single variable 𝑓 (𝑥) = 𝑥2 + 10 sin(𝑥):
We can see that the function has a local minimizer near 𝑥 = 3.8 and a global minimizer near 𝑥 = −1.3,
but the precise values cannot be determined from the plot.
The most appropriate function for this purpose is scipy.optimize.minimize_scalar(). Since we know
the approximate locations of the minima, we will provide bounds that restrict the search to the vicinity
of the global minimum.
If we did not already know the approximate location of the global minimum, we could use one of
SciPy’s global minimizers, such as scipy.optimize.differential_evolution(). We are required to
pass bounds, but they do not need to be tight.
>>> bounds=[(-5, 5)] # list of lower, upper bound for each variable
>>> res = sp.optimize.differential_evolution(f, bounds=bounds)
>>> res
message: Optimization terminated successfully.
success: True
fun: -7.9458233756...
x: [-1.306e+00]
nit: 6
nfev: 111
jac: [ 9.948e-06]
For multivariate optimization, a good choice for many problems is scipy.optimize.minimize(). Sup-
pose we wish to find the minimum of a quadratic function of two variables, 𝑓 (𝑥0 , 𝑥1 ) = (𝑥0 −1)2 +(𝑥1 −2)2 .
Like scipy.optimize.root(), scipy.optimize.minimize() requires a guess x0. (Note that this is the
initial value of both variables rather than the value of the variable we happened to label 𝑥0 .)
Maximization?
This barely scratches the surface of SciPy’s optimization features, which include mixed integer linear
programming, constrained nonlinear programming, and the solution of assignment problems. For much
more information, see the documentation of scipy.optimize and the advanced chapter Mathematical
optimization: finding minima of functions.
𝑥4 2
𝑓 (𝑥, 𝑦) = (4 − 2.1𝑥2 + )𝑥 + 𝑥𝑦 + (4𝑦 2 − 4)𝑦 2
3
has multiple local minima. Find a global minimum (there is more than one, each with the
same value of the objective function) and at least one other local minimum.
Hints:
• Variables can be restricted to −2 < 𝑥 < 2 and −1 < 𝑦 < 1.
• numpy.meshgrid() and matplotlib.pyplot.imshow() can help with visualization.
• Try minimizing with scipy.optimize.minimize() with an initial guess of (𝑥, 𝑦) =
(0, 0). Does it find the global minimum, or converge to a local minimum? What about
other initial guesses?
• Try minimizing with scipy.optimize.differential_evolution().
solution
See the summary exercise on Non linear least squares curve fitting: application to point extraction in
topographical lidar data for another, more advanced example.
Consider a random variable distributed according to the standard normal. We draw a sample consisting
of 100000 observations from the random variable. The normalized histogram of the sample is an estimator
of the random variable’s probability density function (PDF):
Each of the 100+ scipy.stats distribution families is represented by an object with a __call__
method. Here, we call the scipy.stats.norm object to specify its location and scale, and it returns a
frozen distribution: a particular element of a distribution family with all parameters fixed. The frozen
distribution object has methods to compute essential functions of the particular distribution.
Suppose we knew that the sample had been drawn from a distribution belonging to the family of normal
distributions, but we did not know the particular distribution’s location (mean) and scale (standard
deviation). We perform maximum likelihood estimation of the unknown parameters using the distribution
family’s fit method:
Since we know the true parameters of the distribution from which the sample was drawn, we are not
surprised that these estimates are similar.
Generate 1000 random variates from a gamma distribution with a shape parameter of 1. Hint: the
shape parameter is passed as the first argument when freezing the distribution. Plot the histogram of
the sample, and overlay the distribution’s PDF. Estimate the shape parameter from the sample using
the fit method.
Extra: the distributions have many useful methods. Explore them using tab completion. Plot the
cumulative density function of the distribution, and compute the variance.
The sample mean is an estimator of the mean of the distribution from which the sample was drawn:
>>> np.mean(sample)
0.001576700508...
NumPy includes some of the most fundamental sample statistics (e.g. numpy.mean(), numpy.var(),
numpy.percentile()); scipy.stats includes many more. For instance, the geometric mean is a common
measure of central tendency for data that tends to be distributed over many orders of magnitude.
>>> sp.stats.gmean(2**sample)
1.0010934829...
SciPy also includes a variety of hypothesis tests that produce a sample statistic and a p-value. For
instance, suppose we wish to test the null hypothesis that sample was drawn from a normal distribution:
Here, statistic is a sample statistic that tends to be high for samples that are drawn from non-normal
distributions. pvalue is the probability of observing such a high value of the statistic for a sample
that has been drawn from a normal distribution. If the p-value is unusually small, this may be taken
as evidence that sample was not drawn from the normal distribution. Our statistic and p-value are
moderate, so the test is inconclusive.
There are many other features of scipy.stats, including circular statistics, quasi-Monte Carlo methods,
and resampling methods. For much more information, see the documentation of scipy.stats and the
advanced chapter statistics.
5.7.1 Quadrature
∫︀ 𝜋/2
Suppose we wish to compute the definite integral 0 sin(𝑡)𝑑𝑡 numerically. scipy.integrate.quad()
chooses one of several adaptive techniques depending on the parameters, and is therefore the recom-
mended first choice for integration of function of a single variable:
Other functions for numerical quadrature, including integration of multivariate functions and approxi-
mating integrals from samples, are available in scipy.integrate.
scipy.integrate also features routines for integrating Ordinary Differential Equations (ODE). For
example, scipy.integrate.solve_ivp() integrates ODEs of the form:
𝑑𝑦
= 𝑓 (𝑡, 𝑦(𝑡))
𝑑𝑡
from an initial time 𝑡0 and initial state 𝑦(𝑡 = 𝑡0 ) = 𝑡0 to a final time 𝑡𝑓 or until an event occurs (e.g. a
specified state is reached).
As an introduction, consider the initial value problem given by 𝑑𝑦 𝑑𝑡 = −2𝑦 and the initial condition
𝑦(𝑡 = 0) = 1 on the interval 𝑡 = 0 . . . 4. We begin by defining a callable that computes 𝑓 (𝑡, 𝑦(𝑡)) given
the current time and state.
Let us integrate a more complex ODE: a damped spring-mass oscillator. The position of√︀ a mass attached
to a spring obeys the 2nd order ODE 𝑦¨ + 2𝜁𝜔0 𝑦˙ + 𝜔02 𝑦 = 0 with natural frequency 𝜔0 = 𝑘/𝑚, damping
ratio 𝜁 = 𝑐/(2𝑚𝜔0 ), spring constant 𝑘, mass 𝑚, and damping coefficient 𝑐.
Before using scipy.integrate.solve_ivp(), the 2nd order ODE needs to be transformed into a system
of first-order ODEs. Note that
𝑑𝑦 𝑑𝑦˙
= 𝑦˙ = 𝑦¨ = −(2𝜁𝜔0 𝑦˙ + 𝜔02 𝑦)
𝑑𝑡 𝑑𝑡
If we define 𝑧 = [𝑧0 , 𝑧1 ] where 𝑧0 = 𝑦 and 𝑧1 = 𝑦,˙ then the first order equation:
[︂ 𝑑𝑧0 ]︂ [︂ ]︂
𝑑𝑧 𝑑𝑡 𝑧1
= 𝑑𝑧1 =
𝑑𝑡 𝑑𝑡
−(2𝜁𝜔0 𝑧1 + 𝜔02 𝑧0 )
>>> m = 0.5 # kg
>>> k = 4 # N/m
>>> c = 0.4 # N s/m
>>> zeta = c / (2 * m * np.sqrt(k/m))
>>> omega = np.sqrt(k / m)
Tip: With the option method=’LSODA’, scipy.integrate.solve_ivp() uses the LSODA (Liver-
more Solver for Ordinary Differential equations with Automatic method switching for stiff and non-stiff
problems). See the ODEPACK Fortran library for more details.
See also:
Partial Differental Equations
There is no Partial Differential Equations (PDE) solver in SciPy. Some Python packages for solving
PDE’s are available, such as fipy or SfePy.
The scipy.fft module computes fast Fourier transforms (FFTs) and offers utilities to handle them.
Some important functions are:
• scipy.fft.fft() to compute the FFT
• scipy.fft.fftfreq() to generate the sampling frequencies
• scipy.fft.ifft() to compute the inverse FFT, from frequency space to signal space
Signal FFT
As the signal comes from a real-valued function, the Fourier transform is symmetric.
The peak signal frequency can be found with freqs[power.argmax()]
Setting the Fourier component above this frequency to zero and inverting the FFT with scipy.fft.
ifft(), gives a filtered signal.
numpy.fft
NumPy also has an implementation of FFT (numpy.fft). However, the SciPy one should be preferred,
as it uses more efficient underlying implementations.
1. Examine the provided image moonlanding.png, which is heavily contaminated with periodic
noise. In this exercise, we aim to clean up the noise using the Fast Fourier Transform.
2. Load the image using matplotlib.pyplot.imread().
3. Find and use the 2-D FFT function in scipy.fft, and plot the spectrum (Fourier transform of)
the image. Do you have any trouble visualising the spectrum? If so, why?
4. The spectrum consists of high and low frequency components. The noise is contained in the
high-frequency part of the spectrum, so set some of those components to zero (use array slicing).
5. Apply the inverse Fourier transform to see the resulting image.
Solution
>>> plt.plot(t, x)
[<matplotlib.lines.Line2D object at ...>]
>>> plt.plot(t[::4], x_resampled, 'ko')
[<matplotlib.lines.Line2D object at ...>]
Tip: Notice how on the side of the window the resampling is less accurate and has a rippling effect.
This resampling is different from the interpolation provided by scipy.interpolate as it only applies to
regularly sampled data.
>>> plt.plot(t, x)
[<matplotlib.lines.Line2D object at ...>]
>>> plt.plot(t, x_detrended)
[<matplotlib.lines.Line2D object at ...>]
Filtering: For non-linear filtering, scipy.signal has filtering (median filter scipy.signal.medfilt(),
Wiener scipy.signal.wiener()), but we will discuss this in the image section.
Tip: scipy.signal also has a full-blown set of tools for the design of linear filter (finite and infinite
response filters), but this is out of the scope of this tutorial.
>>> plt.subplot(151)
<Axes: >
>>> plt.axis('off')
(-0.5, 1023.5, 767.5, -0.5)
>>> # etc.
Exercise
Tip: Mathematical morphology stems from set theory. It characterizes and transforms geometrical
structures. Binary (black and white) images, in particular, can be transformed using this theory: the
sets to be transformed are the sets of neighboring non-zero-valued pixels. The theory was also extended
to gray-valued images.
>>> el = sp.ndimage.generate_binary_structure(2, 1)
>>> el
array([[False, True, False],
[...True, True, True],
[False, True, False]])
>>> el.astype(int)
array([[0, 1, 0],
[1, 1, 1],
[0, 1, 0]])
• Erosion scipy.ndimage.binary_erosion()
• Dilation scipy.ndimage.binary_dilation()
• Opening scipy.ndimage.binary_opening()
• Closing: scipy.ndimage.binary_closing()
Exercise
An opening operation removes small structures, while a closing operation fills small holes. Such opera-
tions can therefore be used to “clean” an image.
Exercise
Check that the area of the reconstructed square is smaller than the area of the initial
square. (The opposite would occur if the closing step was performed before the opening).
For gray-valued images, eroding (resp. dilating) amounts to replacing a pixel by the minimal (resp.
maximal) value among pixels covered by the structuring element centered on the pixel of interest.
Extract the 4th connected component, and crop the array around it:
>>> sp.ndimage.find_objects(labels)[3]
(slice(30, 48, None), slice(30, 48, None))
>>> sl = sp.ndimage.find_objects(labels)[3]
>>> import matplotlib.pyplot as plt
>>> plt.imshow(sig[sl])
<matplotlib.image.AxesImage object at ...>
See the summary exercise on Image processing application: counting bubbles and unmolten grains for a
more advanced example.
The summary exercises use mainly NumPy, SciPy and Matplotlib. They provide some real-life examples
of scientific computing with Python. Now that the basics of working with NumPy and SciPy have been
introduced, the interested user is invited to try these exercises.
The exercise goal is to predict the maximum wind speed occurring every 50 years even if no measure
exists for such a period. The available data are only measured over 21 years at the Sprogø meteorological
station located in Denmark. First, the statistical steps will be given and then illustrated with functions
from the scipy.interpolate module. At the end the interested readers are invited to compute results from
raw data and in a slightly different approach.
Statistical approach
The annual maxima are supposed to fit a normal probability density function. However such function
is not going to be estimated because it gives a probability from a wind speed maxima. Finding the
maximum wind speed occurring every 50 years requires the opposite approach, the result needs to be
found from a defined probability. That is the quantile function role and the exercise goal will be to find
it. In the current model, it is supposed that the maximum wind speed occurring every 50 years is defined
as the upper 2% quantile.
By definition, the quantile function is the inverse of the cumulative distribution function. The latter
describes the probability distribution of an annual maxima. In the exercise, the cumulative probability
p_i for a given year i is defined as p_i = i/(N+1) with N = 21, the number of measured years. Thus
it will be possible to calculate the cumulative probability of every measured wind speed maxima. From
those experimental points, the scipy.interpolate module will be very useful for fitting the quantile function.
Finally the 50 years maxima is going to be evaluated from the cumulative probability of the 2% quantile.
The annual wind speeds maxima have already been computed and saved in the NumPy format in the
file examples/max-speeds.npy, thus they will be loaded by using NumPy:
Following the cumulative probability definition p_i from the previous section, the corresponding values
will be:
In this section the quantile function will be estimated by using the UnivariateSpline class which can
represent a spline from points. The default behavior is to build a spline of degree 3 and points can
have different weights according to their reliability. Variants are InterpolatedUnivariateSpline and
LSQUnivariateSpline on which errors checking is going to change. In case a 2D spline is wanted, the
BivariateSpline class family is provided. All those classes for 1D and 2D splines use the FITPACK
Fortran subroutines, that’s why a lower library access is available through the splrep and splev functions
for respectively representing and evaluating a spline. Moreover interpolation functions without the use
of FITPACK parameters are also provided for simpler use.
For the Sprogø maxima wind speeds, the UnivariateSpline will be used because a spline of degree 3
seems to correctly fit the data:
The quantile function is now going to be evaluated from the full range of probabilities:
In the current model, the maximum wind speed occurring every 50 years is defined as the upper 2%
quantile. As a result, the cumulative probability value will be:
So the storm wind speed occurring every 50 years can be guessed by:
The interested readers are now invited to make an exercise by using the wind speeds measured over 21
years. The measurement period is around 90 minutes (the original period was around 10 minutes but the
file size has been reduced for making the exercise setup easier). The data are stored in NumPy format
inside the file examples/sprog-windspeeds.npy. Do not look at the source code for the plots until you
have completed the exercise.
• The first step will be to find the annual maxima by using NumPy and plot them as a matplotlib
bar figure.
• The second step will be to use the Gumbell distribution on cumulative probabilities p_i defined as
-log( -log(p_i) ) for fitting a linear quantile function (remember that you can define the degree
of the UnivariateSpline). Plotting the annual maxima versus the Gumbell distribution should
give you the following figure.
• The last step will be to find 34.23 m/s for the maximum wind speed occurring every 50 years.
5.11.2 Non linear least squares curve fitting: application to point extraction in
topographical lidar data
The goal of this exercise is to fit a model to some data. The data used in this tutorial are lidar data
and are described in details in the following introductory paragraph. If you’re impatient and want to
practice now, please skip it and go directly to Loading and visualization.
Introduction
Lidars systems are optical rangefinders that analyze property of scattered light to measure distances.
Most of them emit a short light impulsion towards a target and record the reflected signal. This signal
is then processed to extract the distance between the lidar system and the target.
Topographical lidar systems are such systems embedded in airborne platforms. They measure distances
between the platform and the Earth, so as to deliver information on the Earth’s topography (see1 for
more details).
In this tutorial, the goal is to analyze the waveform recorded by the lidar system2 . Such a signal
contains peaks whose center and amplitude permit to compute the position and some characteristics of
the hit target. When the footprint of the laser beam is around 1m on the Earth surface, the beam can
hit multiple targets during the two-way propagation (for example the ground and the top of a tree or
building). The sum of the contributions of each target hit by the laser beam then produces a complex
signal with multiple peaks, each one containing information about one target.
One state of the art method to extract information from these data is to decompose them in a sum of
Gaussian functions where each function represents the contribution of a target hit by the laser beam.
Therefore, we use the scipy.optimize module to fit a waveform to one or a sum of Gaussian functions.
As shown below, this waveform is a 80-bin-length signal with a single peak with an amplitude of ap-
proximately 30 in the 15 nanosecond bin. Additionally, the base level of noise is approximately 3. These
values can be used in the initial solution.
1 Mallet, C. and Bretar, F. Full-Waveform Topographic Lidar: State-of-the-Art. ISPRS Journal of Photogrammetry
The signal is very simple and can be modeled as a single Gaussian function and an offset corresponding
to the background noise. To fit the signal with the function, we must:
• define the model
• propose an initial solution
• call scipy.optimize.leastsq
Model
where
• coeffs[0] is 𝐵 (noise)
• coeffs[1] is 𝐴 (amplitude)
• coeffs[2] is 𝜇 (center)
• coeffs[3] is 𝜎 (width)
Initial solution
Fit
scipy.optimize.leastsq minimizes the sum of squares of the function given as an argument. Basically,
the function to minimize is the residuals (the difference between the data and the model):
So let’s get our solution by calling scipy.optimize.leastsq() with the following arguments:
• the function to minimize
• an initial solution
• the additional arguments to pass to the function
Remark: from scipy v0.8 and above, you should rather use scipy.optimize.curve_fit() which takes
the model and the data as arguments, so you don’t need to define the residuals any more.
Going further
• Try with a more complex waveform (for instance waveform_2.npy) that contains three significant
peaks. You must adapt the model which is now a sum of Gaussian functions instead of only one
Gaussian peak.
• In some cases, writing an explicit function to compute the Jacobian is faster than letting leastsq
estimate it numerically. Create a function to compute the Jacobian of the residuals and use it as
an input for leastsq.
• When we want to detect very small peaks in the signal, or when the initial guess is too far from a
good solution, the result given by the algorithm is often not satisfying. Adding constraints to the
parameters of the model enables to overcome such limitations. An example of a priori knowledge
we can add is the sign of our variables (which are all positive).
• See the solution.
• Further exercise: compare the result of scipy.optimize.leastsq() and what you can get with
scipy.optimize.fmin_slsqp() when adding boundary constraints.
1. Open the image file MV_HFV_012.jpg and display it. Browse through the keyword arguments
in the docstring of imshow to display the image with the “right” orientation (origin in the bottom
left corner, and not the upper left corner as for standard arrays).
This Scanning Element Microscopy image shows a glass sample (light gray matrix) with some
bubbles (on black) and unmolten sand grains (dark gray). We wish to determine the fraction
of the sample covered by these three phases, and to estimate the typical size of sand grains and
bubbles, their sizes, etc.
2. Crop the image to remove the lower panel with measure information.
3. Slightly filter the image with a median filter in order to refine its histogram. Check how the
histogram changes.
4. Using the histogram of the filtered image, determine thresholds that allow to define masks for sand
pixels, glass pixels and bubble pixels. Other option (homework): write a function that determines
automatically the thresholds from the minima of the histogram.
5. Display an image in which the three phases are colored with three different colors.
6. Use mathematical morphology to clean the different phases.
7. Attribute labels to all bubbles and sand grains, and remove from the sand mask grains that are
smaller than 10 pixels. To do so, use ndimage.sum or np.bincount to compute the grain sizes.
8. Compute the mean size of bubbles.
5.11.4 Example of solution for the image processing exercise: unmolten grains in
glass
1. Open the image file MV_HFV_012.jpg and display it. Browse through the keyword arguments
in the docstring of imshow to display the image with the “right” orientation (origin in the bottom
left corner, and not the upper left corner as for standard arrays).
2. Crop the image to remove the lower panel with measure information.
3. Slightly filter the image with a median filter in order to refine its histogram. Check how the
histogram changes.
4. Using the histogram of the filtered image, determine thresholds that allow to define masks for sand
pixels, glass pixels and bubble pixels. Other option (homework): write a function that determines
automatically the thresholds from the minima of the histogram.
5. Display an image in which the three phases are colored with three different colors.
7. Attribute labels to all bubbles and sand grains, and remove from the sand mask grains that are
smaller than 10 pixels. To do so, use sp.ndimage.sum or np.bincount to compute the grain sizes.
import numpy as np
import matplotlib.pyplot as plt
def f(x):
return x**2 + 10 * np.sin(x)
x = np.arange(-5, 5, 0.1)
plt.plot(x, f(x))
import scipy as sp
plt.show()
import numpy as np
t = np.linspace(0, 5, 100)
x = np.sin(t)
Downsample it by a factor of 4
import scipy as sp
Plot
plt.figure(figsize=(5, 4))
plt.plot(t, x, label="Original signal")
plt.plot(t[::4], x_resampled, "ko", label="Resampled signal")
plt.legend(loc="best")
plt.show()
import numpy as np
t = np.linspace(0, 5, 100)
rng = np.random.default_rng()
x = t + rng.normal(size=100)
Detrend
import scipy as sp
x_detrended = sp.signal.detrend(x)
Plot
plt.figure(figsize=(5, 4))
plt.plot(t, x, label="x")
plt.plot(t, x_detrended, label="x_detrended")
plt.legend(loc="best")
plt.show()
Solve the ODE dy/dt = -2y between t = 0..4, with the initial condition y(t=0) = 1.
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
plt.figure(figsize=(4, 3))
plt.plot(res.t, res.y[0])
plt.xlabel("t")
plt.ylabel("y")
plt.title("Solution of Initial Value Problem")
plt.tight_layout()
plt.show()
Explore the normal distribution: a histogram built from samples and the PDF (probability density
function).
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
m = 0.5 # kg
k = 4 # N/m
c = 0.4 # N s/m
plt.figure(figsize=(4, 3))
plt.plot(res.t, res.y[0], label="y")
plt.plot(res.t, res.y[1], label="dy/dt")
plt.legend(loc="best")
plt.show()
import numpy as np
import matplotlib.pyplot as plt
plt.figure(figsize=(6, 4))
plt.hist(samples1, bins=bins, density=True, label="Samples 1")
plt.hist(samples2, bins=bins, density=True, label="Samples 2")
plt.legend(loc="best")
plt.show()
import numpy as np
# And plot it
import matplotlib.pyplot as plt
plt.figure(figsize=(6, 4))
plt.scatter(x_data, y_data)
import scipy as sp
print(params)
plt.figure(figsize=(6, 4))
plt.scatter(x_data, y_data, label="Data")
plt.plot(x_data, test_func(x_data, *params), label="Fitted function")
plt.legend(loc="best")
plt.show()
import numpy as np
import matplotlib.pyplot as plt
time_step = 0.01
time_vec = np.arange(0, 70, time_step)
plt.figure(figsize=(8, 5))
plt.plot(time_vec, sig)
import scipy as sp
plt.figure(figsize=(5, 4))
plt.imshow(spectrogram, aspect="auto", cmap="hot_r", origin="lower")
plt.title("Spectrogram")
plt.ylabel("Frequency band")
plt.xlabel("Time window")
plt.tight_layout()
plt.figure(figsize=(5, 4))
plt.semilogx(freqs, psd)
plt.title("PSD: power spectral density")
plt.xlabel("Frequency")
plt.ylabel("Power")
plt.tight_layout()
plt.show()
np.random.seed(0)
a = np.zeros((50, 50))
a[10:-10, 10:-10] = 1
a += 0.25 * np.random.standard_normal(a.shape)
mask = a >= 0.5
opened_mask = sp.ndimage.binary_opening(mask)
closed_mask = sp.ndimage.binary_closing(opened_mask)
# Plot
import matplotlib.pyplot as plt
plt.figure(figsize=(12, 3.5))
plt.subplot(141)
plt.imshow(a, cmap=plt.cm.gray)
plt.axis("off")
plt.title("a")
plt.subplot(142)
plt.imshow(mask, cmap=plt.cm.gray)
plt.axis("off")
plt.title("mask")
plt.subplot(143)
plt.imshow(opened_mask, cmap=plt.cm.gray)
plt.axis("off")
plt.title("opened_mask")
plt.subplot(144)
plt.imshow(closed_mask, cmap=plt.cm.gray)
plt.title("closed_mask")
plt.axis("off")
plt.show()
face = sp.datasets.face(gray=True)
plt.figure(figsize=(15, 3))
plt.subplot(151)
plt.imshow(shifted_face, cmap=plt.cm.gray)
plt.axis("off")
plt.subplot(152)
plt.imshow(shifted_face2, cmap=plt.cm.gray)
plt.axis("off")
plt.subplot(153)
plt.imshow(rotated_face, cmap=plt.cm.gray)
plt.axis("off")
plt.subplot(154)
plt.imshow(cropped_face, cmap=plt.cm.gray)
plt.axis("off")
plt.subplot(155)
plt.imshow(zoomed_face, cmap=plt.cm.gray)
plt.axis("off")
plt.show()
import numpy as np
import matplotlib.pyplot as plt
x, y = np.indices((100, 100))
sig = (
np.sin(2 * np.pi * x / 50.0)
* np.sin(2 * np.pi * y / 50.0)
(continues on next page)
plt.figure(figsize=(7, 3.5))
plt.subplot(1, 2, 1)
plt.imshow(sig)
plt.axis("off")
plt.title("sig")
plt.subplot(1, 2, 2)
plt.imshow(mask, cmap=plt.cm.gray)
plt.axis("off")
plt.title("mask")
plt.subplots_adjust(wspace=0.05, left=0.01, bottom=0.01, right=0.99, top=0.9)
import scipy as sp
labels, nb = sp.ndimage.label(mask)
plt.figure(figsize=(3.5, 3.5))
plt.imshow(labels)
plt.title("label")
plt.axis("off")
Extract the 4th connected component, and crop the array around it
sl = sp.ndimage.find_objects(labels == 4)
plt.figure(figsize=(3.5, 3.5))
plt.imshow(sig[sl[0]])
plt.title("Cropped connected component")
plt.axis("off")
plt.show()
import numpy as np
def f(x):
return x**2 + 10 * np.sin(x)
Find minima
import scipy as sp
# Global optimization
grid = (-10, 10, 0.1)
xmin_global = sp.optimize.brute(f, (grid,))
print(f"Global minima found { xmin_global} ")
# Constrain optimization
xmin_local = sp.optimize.fminbound(f, 0, 10)
print(f"Local minimum found { xmin_local} ")
Root finding
face = sp.datasets.face(gray=True)
face = face[:512, -512:] # crop out square on right
import numpy as np
noisy_face = np.copy(face).astype(float)
rng = np.random.default_rng()
noisy_face += face.std() * 0.5 * rng.standard_normal(face.shape)
blurred_face = sp.ndimage.gaussian_filter(noisy_face, sigma=3)
median_face = sp.ndimage.median_filter(noisy_face, size=5)
wiener_face = sp.signal.wiener(noisy_face, (5, 5))
plt.figure(figsize=(12, 3.5))
plt.subplot(141)
plt.imshow(noisy_face, cmap=plt.cm.gray)
plt.axis("off")
plt.title("noisy")
plt.subplot(142)
plt.imshow(blurred_face, cmap=plt.cm.gray)
plt.axis("off")
plt.title("Gaussian filter")
plt.subplot(143)
plt.imshow(median_face, cmap=plt.cm.gray)
plt.axis("off")
plt.title("median filter")
plt.subplot(144)
plt.imshow(wiener_face, cmap=plt.cm.gray)
plt.title("Wiener filter")
plt.axis("off")
plt.show()
import numpy as np
Simple visualization in 2D
plt.figure()
plt.imshow(sixhump([xg, yg]), extent=xlim + ylim, origin="lower")
plt.colorbar()
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
surf = ax.plot_surface(
xg,
yg,
sixhump([xg, yg]),
rstride=1,
cstride=1,
cmap=plt.cm.viridis,
linewidth=0,
antialiased=False,
)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("f(x, y)")
ax.set_title("Six-hump Camelback function")
Find minima
import scipy as sp
# local minimization
res_local = sp.optimize.minimize(sixhump, x0=[0, 0])
# global minimization
res_global = sp.optimize.differential_evolution(sixhump, bounds=[xlim, ylim])
plt.figure()
# Show the function in 2D
plt.imshow(sixhump([xg, yg]), extent=xlim + ylim, origin="lower")
plt.colorbar()
# Mark the minima
plt.scatter(res_local.x[0], res_local.x[1], label="local minimizer")
plt.scatter(res_global.x[0], res_global.x[1], label="global minimizer")
plt.legend()
plt.show()
Plot the power of the FFT of a signal and inverse FFT back to reconstruct a signal.
This example demonstrate scipy.fft.fft(), scipy.fft.fftfreq() and scipy.fft.ifft(). It imple-
ments a basic filter that is very suboptimal, and should not be used.
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
time_step = 0.02
period = 5.0
plt.figure(figsize=(6, 5))
plt.plot(time_vec, sig, label="Original signal")
# Find the peak frequency: we can focus on only the positive frequencies
pos_mask = np.where(sample_freq > 0)
freqs = sample_freq[pos_mask]
peak_freq = freqs[power[pos_mask].argmax()]
[]
We now remove all the high frequencies and transform back from frequencies to signal.
high_freq_fft = sig_fft.copy()
high_freq_fft[np.abs(sample_freq) > peak_freq] = 0
filtered_sig = sp.fft.ifft(high_freq_fft)
plt.figure(figsize=(6, 5))
plt.plot(time_vec, sig, label="Original signal")
(continues on next page)
plt.legend(loc="best")
/opt/hostedtoolcache/Python/3.11.9/x64/lib/python3.11/site-packages/matplotlib/cbook.
˓→py:1762: ComplexWarning: Casting complex values to real discards the imaginary part
return math.isfinite(val)
/opt/hostedtoolcache/Python/3.11.9/x64/lib/python3.11/site-packages/matplotlib/cbook.
˓→py:1398: ComplexWarning: Casting complex values to real discards the imaginary part
Note This is actually a bad way of creating a filter: such brutal cut-off in frequency space does not
control distortion on the signal.
Filters should be created using the SciPy filter design code
plt.show()
•
# Generate data
import numpy as np
rng = np.random.default_rng(27446968)
measured_time = np.linspace(0, 2 * np.pi, 20)
function = np.sin(measured_time)
noise = rng.normal(loc=0, scale=0.1, size=20)
measurements = function + noise
plt.figure(figsize=(6, 4))
plt.plot(measured_time, measurements, ".", ms=6, label="measurements")
plt.plot(interpolation_time, smooth_results, label="smoothing spline")
plt.plot(interpolation_time, np.sin(interpolation_time), "--", label="underlying curve
˓→")
plt.legend()
plt.show()
plt.legend()
plt.show()
plt.legend()
plt.show()
import numpy as np
data = np.loadtxt("../../../../data/populations.txt")
years = data[:, 0]
populations = data[:, 1:]
plt.figure()
plt.plot(years, populations * 1e-3)
plt.xlabel("Year")
plt.ylabel(r"Population number ($\cdot10^3$)")
plt.legend(["hare", "lynx", "carrot"], loc=1)
import scipy as sp
plt.figure()
plt.plot(periods, abs(ft_populations) * 1e-3, "o")
plt.xlim(0, 22)
plt.xlabel("Period")
plt.ylabel(r"Power ($\cdot10^3$)")
plt.show()
/home/runner/work/scientific-python-lectures/scientific-python-lectures/intro/scipy/
˓→examples/solutions/plot_periodicity_finder.py:39: RuntimeWarning: divide by zero␣
˓→encountered in divide
periods = 1 / frequencies
There’s probably a period of around 10 years (obvious from the plot), but for this crude a method,
there’s not enough data to say much more.
Total running time of the script: (0 minutes 0.121 seconds)
We have the min and max temperatures in Alaska for each months of the year. We would like to find a
function to describe this yearly evolution.
For this, we will fit a periodic function.
The data
import numpy as np
temp_max = np.array([17, 19, 21, 28, 33, 38, 37, 37, 31, 23, 19, 18])
temp_min = np.array([-62, -59, -56, -46, -32, -18, -9, -13, -25, -46, -52, -58])
months = np.arange(12)
plt.plot(months, temp_max, "ro")
plt.plot(months, temp_min, "bo")
plt.xlabel("Month")
plt.ylabel("Min and max temperature")
import scipy as sp
plt.figure()
plt.plot(months, temp_max, "ro")
plt.plot(days, yearly_temps(days, *res_max), "r-")
plt.plot(months, temp_min, "bo")
plt.plot(days, yearly_temps(days, *res_min), "b-")
plt.xlabel("Month")
plt.ylabel(r"Temperature ($^\circ$C)")
plt.show()
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
# read image
img = plt.imread("../../../../data/elephant.png")
plt.figure()
plt.imshow(img)
# convolve
img_ft = sp.fft.fft2(img, axes=(0, 1))
# the 'newaxis' is to match to color direction
img2_ft = kernel_ft[:, :, np.newaxis] * img_ft
img2 = sp.fft.ifft2(img2_ft, axes=(0, 1)).real
# plot output
plt.figure()
plt.imshow(img2)
The above exercise was only for didactic reasons: there exists a function in scipy that will do
this for us, and probably do a better job: scipy.signal.fftconvolve()
Note that we still have a decay to zero at the border of the image. Using scipy.ndimage.
gaussian_filter() would get rid of this artifact
plt.show()
import numpy as np
import matplotlib.pyplot as plt
im = plt.imread("../../../../data/moonlanding.png").astype(float)
plt.figure()
plt.imshow(im, plt.cm.gray)
plt.title("Original image")
import scipy as sp
im_fft = sp.fft.fft2(im)
def plot_spectrum(im_fft):
from matplotlib.colors import LogNorm
# A logarithmic colormap
(continues on next page)
plt.figure()
plot_spectrum(im_fft)
plt.title("Fourier transform")
Filter in FFT
# In the lines following, we'll make a copy of the original spectrum and
# truncate coefficients.
plt.figure()
plot_spectrum(im_fft2)
plt.title("Filtered Spectrum")
# Reconstruct the denoised image from the filtered spectrum, keep only the
# real part for display.
im_new = sp.fft.ifft2(im_fft2).real
plt.figure()
plt.imshow(im_new, plt.cm.gray)
plt.title("Reconstructed Image")
Implementing filtering directly with FFTs is tricky and time consuming. We can use the
Gaussian filter from scipy.ndimage
im_blur = sp.ndimage.gaussian_filter(im, 4)
plt.figure()
plt.imshow(im_blur, plt.cm.gray)
plt.title("Blurred image")
plt.show()
See also:
References to go further
• Some chapters of the advanced and the packages and applications parts of the SciPy lectures
• The SciPy cookbook
In [1]: help(np.van<TAB>
In [2]: help(np.vander)
Help on _ArrayFunctionDispatcher in module numpy:
The columns of the output matrix are powers of the input vector. The
order of the powers is determined by the `increasing` boolean argument.
Specifically, when `increasing` is False, the `i`-th output column is
the input vector raised element-wise to the power of ``N - i - 1``. Such
a matrix with a geometric progression in each row is named for Alexandre-
Theophile Vandermonde.
Parameters
----------
x : array_like
1-D input array.
N : int, optional
Number of columns in the output. If `N` is not specified, a square
array is returned (``N = len(x)``).
(continues on next page)
292
Scientific Python Lectures, Edition 2024.2rc0.dev0
.. versionadded:: 1.9.0
Returns
-------
out : ndarray
Vandermonde matrix. If `increasing` is False, the first column is
``x^(N-1)``, the second ``x^(N-2)`` and so forth. If `increasing` is
True, the columns are ``x^0, x^1, ..., x^(N-1)``.
See Also
--------
polynomial.polynomial.polyvander
Examples
--------
>>> x = np.array([1, 2, 3, 5])
>>> N = 3
>>> np.vander(x, N)
array([[ 1, 1, 1],
[ 4, 2, 1],
[ 9, 3, 1],
[25, 5, 1]])
>>> np.linalg.det(np.vander(x))
48.000000000000043 # may vary
>>> (5-3)*(5-2)*(5-1)*(3-2)*(3-1)*(2-1)
48
In Ipython it is not possible to open a separated window for help and documentation; however one can
always open a second Ipython shell just to display help and docstrings. . .
• Numpy’s and Scipy’s documentations can be browsed online on https://scipy.org and https://
293
Scientific Python Lectures, Edition 2024.2rc0.dev0
numpy.org. The search button is quite useful inside the reference documentation of the two
packages.
Tutorials on various topics as well as the complete API with all docstrings are found on this website.
• Numpy’s and Scipy’s documentation is enriched and updated on a regular basis by users on a wiki
https://numpy.org/doc/stable/. As a result, some docstrings are clearer or more detailed on the
wiki, and you may want to read directly the documentation on the wiki instead of the official
documentation website. Note that anyone can create an account on the wiki and write better
documentation; this is an easy way to contribute to an open-source project and improve the tools
you are using!
• The SciPy Cookbook https://scipy-cookbook.readthedocs.io gives recipes on many common prob-
lems frequently encountered, such as fitting data points, solving ODE, etc.
• Matplotlib’s website https://matplotlib.org/ features a very nice gallery with a large number of
plots, each of them shows both the source code and the resulting plot. This is very useful for
learning by example. More standard documentation is also available.
Finally, two more “technical” possibilities are useful as well:
• In Ipython, the magical function %psearch search for objects matching patterns. This is useful if,
for example, one does not know the exact name of a function.
In [4]: np.lookfor('convolution')
Search results for 'convolution'
--------------------------------
numpy.convolve
Returns the discrete, linear convolution of two one-dimensional sequences.
numpy.ma.convolve
Returns the discrete, linear convolution of two one-dimensional sequences.
numpy.polymul
Find the product of two polynomials.
numpy.bartlett
Return the Bartlett window.
numpy.correlate
Cross-correlation of two 1-dimensional sequences.
numpy.vectorize
vectorize(pyfunc=np._NoValue, otypes=None, doc=None, excluded=None,
• If everything listed above fails (and Google doesn’t have the answer). . . don’t despair! There is
a vibrant Scientific Python community. Scientific Python is present on various platform. https:
//scientific-python.org/community/
Packages like SciPy and NumPy also have their own channels. Have a look at their respective
websites to find out how to engage with users and maintainers.
294
Part II
Advanced topics
295
Scientific Python Lectures, Edition 2024.2rc0.dev0
This part of the Scientific Python Lectures is dedicated to advanced usage. It strives to educate the
proficient Python coder to be an expert and tackles various specific topics.
296
CHAPTER 7
Advanced Python Constructs
Chapter contents
297
Scientific Python Lectures, Edition 2024.2rc0.dev0
7.1.1 Iterators
Simplicity
Duplication of effort is wasteful, and replacing the various home-grown approaches with a standard
feature usually ends up making things more readable, and interoperable as well.
Guido van Rossum — Adding Optional Static Typing to Python
An iterator is an object adhering to the iterator protocol — basically this means that it has a next
method, which, when called, returns the next item in the sequence, and when there’s nothing to return,
raises the StopIteration exception.
An iterator object allows to loop just once. It holds the state (position) of a single iteration, or from the
other side, each loop over a sequence requires a single iterator object. This means that we can iterate
over the same sequence more than once concurrently. Separating the iteration logic from the sequence
allows us to have more than one way of iteration.
Calling the __iter__ method on a container to create an iterator object is the most straightforward way
to get hold of an iterator. The iter function does that for us, saving a few keystrokes.
>>> nums = [1, 2, 3] # note that ... varies: these are different objects
>>> iter(nums)
<...iterator object at ...>
>>> nums.__iter__()
<...iterator object at ...>
>>> nums.__reversed__()
<...reverseiterator object at ...>
>>> it = iter(nums)
>>> next(it)
1
>>> next(it)
2
>>> next(it)
3
>>> next(it)
Traceback (most recent call last):
(continues on next page)
When used in a loop, StopIteration is swallowed and causes the loop to finish. But with explicit
invocation, we can see that once the iterator is exhausted, accessing it raises an exception.
Using the for..in loop also uses the __iter__ method. This allows us to transparently start the iteration
over a sequence. But if we already have the iterator, we want to be able to use it in an for loop in
the same way. In order to achieve this, iterators in addition to next are also required to have a method
called __iter__ which returns the iterator (self).
Support for iteration is pervasive in Python: all sequences and unordered containers in the standard
library allow this. The concept is also stretched to other things: e.g. file objects support iteration over
lines.
The file is an iterator itself and it’s __iter__ method doesn’t create a separate object: only a single
thread of sequential access is allowed.
A second way in which iterator objects are created is through generator expressions, the basis for list
comprehensions. To increase clarity, a generator expression must always be enclosed in parentheses
or an expression. If round parentheses are used, then a generator iterator is created. If rectangular
parentheses are used, the process is short-circuited and we get a list.
The list comprehension syntax also extends to dictionary and set comprehensions. A set is cre-
ated when the generator expression is enclosed in curly braces. A dict is created when the generator
expression contains “pairs” of the form key:value:
One gotcha should be mentioned: in old Pythons the index variable (i) would leak, and in versions >=
3 this is fixed.
7.1.3 Generators
Generators
A third way to create iterator objects is to call a generator function. A generator is a function containing
the keyword yield. It must be noted that the mere presence of this keyword completely changes the
nature of the function: this yield statement doesn’t have to be invoked, or even reachable, but causes
the function to be marked as a generator. When a normal function is called, the instructions contained in
the body start to be executed. When a generator is called, the execution stops before the first instruction
in the body. An invocation of a generator function creates a generator object, adhering to the iterator
protocol. As with normal function invocations, concurrent and recursive invocations are allowed.
When next is called, the function is executed until the first yield. Each encountered yield statement
gives a value becomes the return value of next. After executing the yield statement, the execution of
this function is suspended.
Let’s go over the life of the single invocation of the generator function.
Contrary to a normal function, where executing f() would immediately cause the first print to be
executed, gen is assigned without executing any statements in the function body. Only when gen.
__next__() is invoked by next, the statements up to the first yield are executed. The second next
prints -- finish -- and execution halts on the second yield. The third next falls of the end of the
Each yield statement causes a value to be passed to the caller. This is the reason for the introduction
of generators by PEP 255. But communication in the reverse direction is also useful. One obvious way
would be some external state, either a global variable or a shared mutable object. Direct communication
is possible thanks to PEP 342. It is achieved by turning the previously boring yield statement into an
expression. When the generator resumes execution after a yield statement, the caller can call a method
on the generator object to either pass a value into the generator, which then is returned by the yield
statement, or a different method to inject an exception into the generator.
The first of the new methods is send(value), which is similar to next(), but passes value into the
generator to be used for the value of the yield expression. In fact, g.next() and g.send(None) are
equivalent.
The second of the new methods is throw(type, value=None, traceback=None) which is equivalent to:
>>> it = g()
>>> next(it)
--start--
--yielding 0--
0
>>> it.send(11)
--yield returned 11--
--yielding 1--
1
>>> it.throw(IndexError)
--yield raised IndexError()--
--yielding 2--
2
>>> it.close()
--closing--
Note: This is a preview of PEP 380 (not yet implemented, but accepted for Python 3.3).
Let’s say we are writing a generator and we want to yield a number of values generated by a second
generator, a subgenerator. If yielding of values is the only concern, this can be performed without
much difficulty using a loop such as
subgen = some_other_generator()
for v in subgen:
yield v
However, if the subgenerator is to interact properly with the caller in the case of calls to send(), throw()
and close(), things become considerably more difficult. The yield statement has to be guarded by a
try..except..finally structure similar to the one defined in the previous section to “debug” the generator
function. Such code is provided in PEP 380#id13, here it suffices to say that new syntax to properly
yield from a subgenerator is being introduced in Python 3.3:
This behaves like the explicit loop above, repeatedly yielding values from some_other_generator until
it is exhausted, but also forwards send, throw and close to the subgenerator.
7.2 Decorators
Summary
This amazing feature appeared in the language almost apologetically and with concern that it might
not be that useful.
Bruce Eckel — An Introduction to Python Decorators
Since functions and classes are objects, they can be passed around. Since they are mutable objects, they
can be modified. The act of altering a function or class object after it has been constructed but before
is is bound to its name is called decorating.
There are two things hiding behind the name “decorator” — one is the function which does the work of
decorating, i.e. performs the real work, and the other one is the expression adhering to the decorator
syntax, i.e. an at-symbol and the name of the decorating function.
Function can be decorated by using the decorator syntax for functions:
@decorator # ❷
def function(): # ❶
pass
def function(): # ❶
pass
function = decorator(function) # ❷
Decorators can be stacked — the order of application is bottom-to-top, or inside-out. The semantics
are such that the originally defined function is used as an argument for the first decorator, whatever is
returned by the first decorator is used as an argument for the second decorator, . . . , and whatever is
returned by the last decorator is attached under the name of the original function.
The decorator syntax was chosen for its readability. Since the decorator is specified before the header
of the function, it is obvious that its is not a part of the function body and its clear that it can only
operate on the whole function. Because the expression is prefixed with @ is stands out and is hard to
miss (“in your face”, according to the PEP :) ). When more than one decorator is applied, each one is
placed on a separate line in an easy to read way.
Decorators can either return the same function or class object or they can return a completely different
object. In the first case, the decorator can exploit the fact that function and class objects are mutable
and add attributes, e.g. add a docstring to a class. A decorator might do something useful even without
modifying the object, for example register the decorated class in a global registry. In the second case,
virtually anything is possible: when something different is substituted for the original function or class,
the new object can be completely different. Nevertheless, such behaviour is not the purpose of decorators:
they are intended to tweak the decorated object, not do something unpredictable. Therefore, when a
function is “decorated” by replacing it with a different function, the new function usually calls the original
function, after doing some preparatory work. Likewise, when a class is “decorated” by replacing if with
a new class, the new class is usually derived from the original class. When the purpose of the decorator
is to do something “every time”, like to log every call to a decorated function, only the second type of
decorators can be used. On the other hand, if the first type is sufficient, it is better to use it, because it
is simpler.
The only requirement on decorators is that they can be called with a single argument. This means that
decorators can be implemented as normal functions, or as classes with a __call__ method, or in theory,
even as lambda functions.
Let’s compare the function and class approaches. The decorator expression (the part after @) can be
either just a name, or a call. The bare-name approach is nice (less to type, looks cleaner, etc.), but is
only possible when no arguments are needed to customise the decorator. Decorators written as functions
can be used in those two cases:
The two trivial decorators above fall into the category of decorators which return the original function.
If they were to return a new function, an extra level of nestedness would be required. In the worst case,
three levels of nested functions.
The _wrapper function is defined to accept all positional and keyword arguments. In general we cannot
know what arguments the decorated function is supposed to accept, so the wrapper function just passes
everything to the wrapped function. One unfortunate consequence is that the apparent argument list is
misleading.
Compared to decorators defined as functions, complex decorators defined as classes are simpler. When
an object is created, the __init__ method is only allowed to return None, and the type of the created
object cannot be changed. This means that when a decorator is defined as a class, it doesn’t make
much sense to use the argument-less form: the final decorated object would just be an instance of the
decorating class, returned by the constructor call, which is not very useful. Therefore it’s enough to
discuss class-based decorators where arguments are given in the decorator expression and the decorator
__init__ method is used for decorator construction.
Contrary to normal rules (PEP 8) decorators written as classes behave more like functions and therefore
their name often starts with a lowercase letter.
In reality, it doesn’t make much sense to create a new class just to have a decorator which returns the
original function. Objects are supposed to hold state, and such decorators are more useful when the
decorator returns a new object.
A decorator like this can do pretty much anything, since it can modify the original function object and
mangle the arguments, call the original function or not, and afterwards mangle the return value.
7.2.3 Copying the docstring and other attributes of the original function
When a new function is returned by the decorator to replace the original function, an unfortunate
consequence is that the original function name, the original docstring, the original argument list are
lost. Those attributes of the original function can partially be “transplanted” to the new function
by setting __doc__ (the docstring), __module__ and __name__ (the full name of the function), and
__annotations__ (extra information about arguments and the return value of the function available in
Python 3). This can be done automatically by using functools.update_wrapper.
functools.update_wrapper(wrapper, wrapped)
>>> function
<function function at 0x...>
>>> print(function.__doc__)
extensive documentation
One important thing is missing from the list of attributes which can be copied to the replacement
function: the argument list. The default values for arguments can be modified through the __defaults__,
__kwdefaults__ attributes, but unfortunately the argument list itself cannot be set as an attribute. This
means that help(function) will display a useless argument list which will be confusing for the user of
the function. An effective but ugly way around this problem is to create the wrapper dynamically, using
eval. This can be automated by using the external decorator module. It provides support for the
decorator decorator, which takes a wrapper and turns it into a decorator which preserves the function
signature.
To sum things up, decorators should always use functools.update_wrapper or some other means of
copying function attributes.
First, it should be mentioned that there’s a number of useful decorators available in the standard library.
There are three decorators which really form a part of the language:
• classmethod causes a method to become a “class method”, which means that it can be invoked
without creating an instance of the class. When a normal method is invoked, the interpreter inserts
the instance object as the first positional parameter, self. When a class method is invoked, the
class itself is given as the first parameter, often called cls.
Class methods are still accessible through the class’ namespace, so they don’t pollute the module’s
namespace. Class methods can be used to provide alternative constructors:
class Array(object):
def __init__(self, data):
self.data = data
@classmethod
def fromfile(cls, file):
data = numpy.load(file)
return cls(data)
In this example, A.a is an read-only attribute. It is also documented: help(A) includes the
docstring for attribute a taken from the getter method. Defining a as a property allows it to be a
calculated on the fly, and has the side effect of making it read-only, because no setter is defined.
To have a setter and a getter, two methods are required, obviously:
class Rectangle(object):
def __init__(self, edge):
self.edge = edge
@property
def area(self):
"""Computed area.
@area.setter
def area(self, area):
self.edge = area ** 0.5
The way that this works, is that the property decorator replaces the getter method with a property
object. This object in turn has three methods, getter, setter, and deleter, which can be used
as decorators. Their job is to set the getter, setter and deleter of the property object (stored as
attributes fget, fset, and fdel). The getter can be set like in the example above, when creating
the object. When defining the setter, we already have the property object under area, and we add
the setter to it by using the setter method. All this happens when we are creating the class.
Afterwards, when an instance of the class has been created, the property object is special. When the
interpreter executes attribute access, assignment, or deletion, the job is delegated to the methods
of the property object.
To make everything crystal clear, let’s define a “debug” example:
Properties are a bit of a stretch for the decorator syntax. One of the premises of the decorator
syntax — that the name is not duplicated — is violated, but nothing better has been invented so
far. It is just good style to use the same name for the getter, setter, and deleter methods.
Some newer examples include:
• functools.lru_cache memoizes an arbitrary function maintaining a limited cache of argu-
ments:answer pairs (Python 3.2)
• functools.total_ordering is a class decorator which fills in missing ordering methods (__lt__,
__gt__, __le__, . . . ) based on a single available one.
Let’s say we want to print a deprecation warning on stderr on the first invocation of a function we don’t
like anymore. If we don’t want to modify the function, we can use a decorator:
class deprecated(object):
"""Print a deprecation warning once on first use of the function.
def deprecated(func):
"""Print a deprecation warning once on first use of the function.
Let’s say we have function which returns a lists of things, and this list created by running a loop. If we
don’t know how many objects will be needed, the standard way to do this is something like:
def find_answers():
answers = []
while True:
ans = look_for_next_answer()
if ans is None:
break
answers.append(ans)
return answers
This is fine, as long as the body of the loop is fairly compact. Once it becomes more complicated, as
often happens in real code, this becomes pretty unreadable. We could simplify this by using yield
statements, but then the user would have to explicitly call list(find_answers()).
We can define a decorator which constructs the list for us:
def vectorized(generator_func):
def wrapper(*args, **kwargs):
return list(generator_func(*args, **kwargs))
return functools.update_wrapper(wrapper, generator_func)
This is a class decorator which doesn’t modify the class, but just puts it in a global registry. It falls into
the category of decorators returning the original object:
class WordProcessor(object):
PLUGINS = []
def process(self, text):
for plugin in self.PLUGINS:
text = plugin().cleanup(text)
return text
@classmethod
def plugin(cls, plugin):
cls.PLUGINS.append(plugin)
(continues on next page)
@WordProcessor.plugin
class CleanMdashesExtension(object):
def cleanup(self, text):
return text.replace('—', u'\N{em dash}')
Here we use a decorator to decentralise the registration of plugins. We call our decorator with a noun,
instead of a verb, because we use it to declare that our class is a plugin for WordProcessor. Method
plugin simply appends the class to the list of plugins.
A word about the plugin itself: it replaces HTML entity for em-dash with a real Unicode em-dash
character. It exploits the unicode literal notation to insert a character by using its name in the unicode
database (“EM DASH”). If the Unicode character was inserted directly, it would be impossible to
distinguish it from an en-dash in the source of a program.
See also:
More examples and reading
• PEP 318 (function and method decorator syntax)
• PEP 3129 (class decorator syntax)
• https://wiki.python.org/moin/PythonDecoratorLibrary
• https://docs.python.org/dev/library/functools.html
• https://pypi.org/project/decorator
• Bruce Eckel
– Decorators I: Introduction to Python Decorators
– Python Decorators II: Decorator Arguments
– Python Decorators III: A Decorator-Based Build System
A context manager is an object with __enter__ and __exit__ methods which can be used in the with
statement:
var = manager.__enter__()
try:
do_something(var)
finally:
manager.__exit__()
In other words, the context manager protocol defined in PEP 343 permits the extraction of the boring
part of a try..except..finally structure into a separate class leaving only the interesting do_something
block.
1. The __enter__ method is called first. It can return a value which will be assigned to var. The
as-part is optional: if it isn’t present, the value returned by __enter__ is simply ignored.
2. The block of code underneath with is executed. Just like with try clauses, it can either execute
successfully to the end, or it can break, continue or return, or it can throw an exception. Either
way, after the block is finished, the __exit__ method is called. If an exception was thrown,
the information about the exception is passed to __exit__, which is described below in the next
subsection. In the normal case, exceptions can be ignored, just like in a finally clause, and will
be rethrown after __exit__ is finished.
Let’s say we want to make sure that a file is closed immediately after we are done writing to it:
Here we have made sure that the f.close() is called when the with block is exited. Since closing files is
such a common operation, the support for this is already present in the file class. It has an __exit__
method which calls close and can be used as a context manager itself:
The common use for try..finally is releasing resources. Various different cases are implemented
similarly: in the __enter__ phase the resource is acquired, in the __exit__ phase it is released, and the
exception, if thrown, is propagated. As with files, there’s often a natural operation to perform after the
object has been used and it is most convenient to have the support built in. With each release, Python
provides support in more places:
• all file-like objects:
– file ➥ automatically closed
– fileinput, tempfile
– bz2.BZ2File, gzip.GzipFile, tarfile.TarFile, zipfile.ZipFile
– ftplib, nntplib ➥ close connection
• locks
– multiprocessing.RLock ➥ lock and unlock
– multiprocessing.Semaphore
– memoryview ➥ automatically release
• decimal.localcontext ➥ modify precision of computations temporarily
• _winreg.PyHKEY ➥ open and close hive key
• warnings.catch_warnings ➥ kill warnings temporarily
• contextlib.closing ➥ the same as the example above, call close
• parallel programming
– concurrent.futures.ThreadPoolExecutor ➥ invoke in parallel then kill thread pool
– concurrent.futures.ProcessPoolExecutor ➥ invoke in parallel then kill process pool
– nogil ➥ solve the GIL problem temporarily (cython only :( )
When an exception is thrown in the with-block, it is passed as arguments to __exit__. Three arguments
are used, the same as returned by sys.exc_info(): type, value, traceback. When no exception is thrown,
None is used for all three arguments. The context manager can “swallow” the exception by returning a
true value from __exit__. Exceptions can be easily ignored, because if __exit__ doesn’t use return
and just falls of the end, None is returned, a false value, and therefore the exception is rethrown after
__exit__ is finished.
The ability to catch exceptions opens interesting possibilities. A classic example comes from unit-tests
— we want to make sure that some code throws the right kind of exception:
class assert_raises(object):
# based on pytest and unittest.TestCase
def __init__(self, type):
self.type = type
def __enter__(self):
pass
def __exit__(self, type, value, traceback):
if type is None:
raise AssertionError('exception expected')
if issubclass(type, self.type):
return True # swallow the expected exception
raise AssertionError('wrong exception type')
with assert_raises(KeyError):
{}['foo']
When discussing generators, it was said that we prefer generators to iterators implemented as classes
because they are shorter, sweeter, and the state is stored as local, not instance, variables. On the other
hand, as described in Bidirectional communication, the flow of data between the generator and its caller
can be bidirectional. This includes exceptions, which can be thrown into the generator. We would like to
implement context managers as special generator functions. In fact, the generator protocol was designed
to support this use case.
@contextlib.contextmanager
def some_generator(<arguments>):
<setup>
try:
yield <value>
finally:
<cleanup>
The contextlib.contextmanager helper takes a generator and turns it into a context manager. The
generator has to obey some rules which are enforced by the wrapper function — most importantly it
must yield exactly once. The part before the yield is executed from __enter__, the block of code
protected by the context manager is executed when the generator is suspended in yield, and the rest
is executed in __exit__. If an exception is thrown, the interpreter hands it to the wrapper through
__exit__ arguments, and the wrapper function then throws it at the point of the yield statement.
Through the use of generators, the context manager is shorter and simpler.
Let’s rewrite the closing example as a generator:
@contextlib.contextmanager
def closing(obj):
(continues on next page)
@contextlib.contextmanager
def assert_raises(type):
try:
yield
except type:
return
except Exception as value:
raise AssertionError('wrong exception type')
else:
raise AssertionError('exception expected')
Prerequisites
• NumPy
• Cython
• Pillow (Python imaging library, used in a couple of examples)
Chapter contents
• Life of ndarray
– It’s. . .
315
Scientific Python Lectures, Edition 2024.2rc0.dev0
– Block of memory
– Data types
– Indexing scheme: strides
– Findings in dissection
• Universal functions
– What they are?
– Exercise: building an ufunc from scratch
– Solution: building an ufunc from scratch
– Generalized ufuncs
• Interoperability features
– Sharing multidimensional, typed data
– The old buffer protocol
– The old buffer protocol
– Array interface protocol
• Array siblings: chararray, maskedarray
– chararray: vectorized string operations
– masked_array missing data
– recarray: purely convenience
• Summary
• Contributing to NumPy/SciPy
– Why
– Reporting bugs
– Contributing to documentation
– Contributing features
– How to help, in general
8.1.1 It’s. . .
ndarray =
block of memory + indexing scheme + data type descriptor
• raw data
• how to locate an element
• how to interpret an element
/* Block of memory */
char *data;
/* Indexing scheme */
int nd;
npy_intp *dimensions;
npy_intp *strides;
/* Other stuff */
PyObject *base;
int flags;
PyObject *weakreflist;
} PyArrayObject;
>>> x.__array_interface__['data'][0]
64803824
>>> x.__array_interface__
{'data': (..., False), 'strides': None, 'descr': [('', '<i4')], 'typestr': '<i4',
˓→'shape': (3,), 'version': 3}
>>> x = b'1234'
x is a string (in Python 3 a bytes), we can represent its data as an array of ints:
>>> y.flags
C_CONTIGUOUS : True
F_CONTIGUOUS : True
OWNDATA : False
WRITEABLE : False
ALIGNED : True
WRITEBACKIFCOPY : False
The owndata and writeable flags indicate status of the memory block.
See also:
array interface
The descriptor
>>> np.dtype(int).type
<class 'numpy.int64'>
>>> np.dtype(int).itemsize
8
>>> np.dtype(int).byteorder
'='
chunk_id "RIFF"
chunk_size 4-byte unsigned little-endian integer
format "WAVE"
fmt_id "fmt "
fmt_size 4-byte unsigned little-endian integer
audio_fmt 2-byte unsigned little-endian integer
num_channels 2-byte unsigned little-endian integer
sample_rate 4-byte unsigned little-endian integer
byte_rate 4-byte unsigned little-endian integer
block_align 2-byte unsigned little-endian integer
bits_per_sample 2-byte unsigned little-endian integer
data_id "data"
data_size 4-byte unsigned little-endian integer
See also:
wavreader.py
>>> wav_header_dtype['format']
dtype('S4')
>>> wav_header_dtype.fields
mappingproxy({'chunk_id': (dtype('S4'), 0), 'chunk_size': (dtype('uint32'), 4),
˓→'format': (dtype('S4'), 8), 'fmt_id': (dtype('S4'), 12), 'fmt_size': (dtype('uint32
>>> wav_header_dtype.fields['format']
(dtype('S4'), 8)
• The first element is the sub-dtype in the structured data, corresponding to the name format
• The second one is its offset (in bytes) from the beginning of the item
Exercise
Mini-exercise, make a “sparse” dtype by using offsets, and only some of the fields:
>>> wav_header_dtype = np.dtype(dict(
... names=['format', 'sample_rate', 'data_id'],
... offsets=[offset_1, offset_2, offset_3], # counted from start of structure in␣
˓→bytes
and use that to read the sample rate, and data_id (as sub-array).
>>> wav_header['sample_rate']
array([16000], dtype=uint32)
>>> wav_header['data_id']
array([[['d', 'a'],
['t', 'a']]],
dtype='|S1')
>>> wav_header.shape
(1,)
>>> wav_header['data_id'].shape
(1, 2, 2)
Note: There are existing modules such as wavfile, audiolab, etc. for loading sound data. . .
casting
• on assignment
• on array construction
• on arithmetic
• etc.
• and manually: .astype(dtype)
data re-interpretation
• manually: .view(dtype)
Casting
Re-interpretation / viewing
– 4 of uint8, OR,
– 4 of int8, OR,
– 2 of int16, OR,
– 1 of int32, OR,
– 1 of float32, OR,
– ...
How to switch from one to another?
1. Switch the dtype:
>>> y = x.view("<i4")
>>> y
array([67305985], dtype=int32)
>>> 0x04030201
67305985
Note:
• .view() makes views, does not copy (or alter) the memory block
• only changes the dtype (and adjusts array shape):
>>> x[1] = 5
>>> y
array([328193], dtype=int32)
>>> y.base is x
True
See also:
view-colors.py
You have RGBA data in an array:
where the last three dimensions are the R, B, and G, and alpha channels.
How to make a (10, 10) structured array with field names ‘r’, ‘g’, ‘b’, ‘a’ without copying data?
>>> y = ...
Solution
We ask NumPy to interpret these bytes as elements of dtype int16—each of which occupies two
bytes in memory. Therefore, 0x01 0x03 becomes the first uint16 and 0x02 0x04 the second.
You may then expect to see 0x0103 (259, when converting from hexadecimal to decimal) as the first
result. But your computer likely stores most significant bytes first, and as such reads the number as
0x0301 or 769 (go on and type 0x0301 into your Python terminal to verify).
We can do the same on a copy of y (why doesn’t it work on y directly?):
>>> y.copy().view(np.int16)
array([[ 513],
[1027]], dtype=int16)
Can you explain these numbers, 513 and 1027, as well as the output shape of the resulting array?
Main point
The question:
>>> x.strides
(3, 1)
>>> byte_offset = 3 * 1 + 1 * 2 # to find x[1, 2]
>>> x.flat[byte_offset]
6
>>> x[1, 2]
6
simple, flexible
Note: The Python built-in bytes returns bytes in C-order by default which can cause confusion when
trying to inspect memory layout. We use numpy.ndarray.tobytes() with order=A instead, which
preserves the C or F ordering of the bytes in memory.
𝑠𝐹
𝑗 = 𝑑1 𝑑2 ...𝑑𝑗−1 × itemsize
Transposition does not affect the memory layout of the data, only strides
>>> x.strides
(2, 1)
>>> y.strides
(1, 2)
>>> x.tobytes('A')
b'\x01\x02\x03\x04'
>>> y.tobytes('A')
b'\x01\x03\x02\x04'
• Everything can be represented by changing only shape, strides, and possibly adjusting the data
pointer!
• Never makes copies of the data
>>> y = x[2:]
>>> y.__array_interface__['data'][0] - x.__array_interface__['data'][0]
8
(continues on next page)
But: not all reshaping operations can be represented by playing with strides:
>>> bytes(a.data)
b'\x00\x01\x02\x03\x04\x05'
>>> b
array([[0, 2, 4],
[1, 3, 5]], dtype=int8)
>>> c = b.reshape(3*2)
>>> c
array([0, 2, 4, 1, 3, 5], dtype=int8)
Here, there is no way to represent the array c given one stride and the block of memory for a. Therefore,
the reshape operation needs to make a copy here.
Stride manipulation
Warning: as_strided does not check that you stay inside the memory block bounds. . .
See also:
stride-fakedims.py
Exercise
Spoiler
Stride can also be 0 :
Broadcasting
• Doing something useful with it: outer product of [1, 2, 3, 4] and [5, 6, 7]
>>> x2 * y2
array([[ 5, 10, 15, 20],
[ 6, 12, 18, 24],
[ 7, 14, 21, 28]], dtype=int16)
See also:
stride-diagonals.py
Challenge
• Pick diagonal entries of the matrix: (assume C memory order):
However,
>>> y.flags.owndata
False
See also:
stride-diagonals.py
Challenge
Compute the tensor trace:
>>> x = np.arange(5*5*5*5).reshape(5, 5, 5, 5)
>>> s = 0
>>> for i in range(5):
... for j in range(5):
... s += x[j, i, j, i]
Solution
In [1]: x = np.zeros((20000,))
In [2]: y = np.zeros((20000*67,))[::67]
Parts of an Ufunc
1. Provided by user
char types[3]
Making it easier
𝑧 ← 𝑧2 + 𝑐
where 𝑐 = 𝑥 + 𝑖𝑦 is a complex number. This iteration is repeated – if 𝑧 stays finite no matter how long
the iteration runs, 𝑐 belongs to the Mandelbrot set.
• Make ufunc called mandel(z0, c) that computes:
z = z0
for k in range(iterations):
z = z*z + c
say, 100 iterations or until z.real**2 + z.imag**2 > 1000. Use it to determine which c are in
the Mandelbrot set.
• Our function is a simple one, so make use of the PyUFunc_* helpers.
• Write it in Cython
See also:
mandel.pyx, mandelplot.py
#
# Fix the parts marked by TODO
#
#
# Compile this file by (Cython >= 0.12 required because of the complex vars)
#
# cython mandel.pyx
# python setup.py build_ext -i
#
# and try it out with, in this directory,
#
# >>> import mandel
# >>> mandel.mandel(0, 1 + 2j)
#
#
#
# Some points of note:
#
# - It's *NOT* allowed to call any Python functions here.
#
# The Ufunc loop runs with the Python Global Interpreter Lock released.
# Hence, the ``nogil``.
#
# - And so all local variables must be declared with ``cdef``
#
# - Note also that this function receives *pointers* to the data
#
#
# TODO: write the Mandelbrot iteration for one point here,
# as you would write it in Python.
#
# Say, use 100 as the maximum number of iterations, and 1000
# as the cutoff for z.real**2 + z.imag**2.
#
import_array()
import_ufunc()
#
# Reminder: some pre-made Ufunc loops:
#
# ================ =======================================================
# ``PyUfunc_f_f`` ``float elementwise_func(float input_1)``
# ``PyUfunc_ff_f`` ``float elementwise_func(float input_1, float input_2)``
# ``PyUfunc_d_d`` ``double elementwise_func(double input_1)``
# ``PyUfunc_dd_d`` ``double elementwise_func(double input_1, double input_2)``
# ``PyUfunc_D_D`` ``elementwise_func(complex_double *input, complex_double* complex_
˓→double)``
# ================ =======================================================
#
# The full list is above.
#
#
# Type codes:
#
# NPY_BOOL, NPY_BYTE, NPY_UBYTE, NPY_SHORT, NPY_USHORT, NPY_INT, NPY_UINT,
# NPY_LONG, NPY_ULONG, NPY_LONGLONG, NPY_ULONGLONG, NPY_FLOAT, NPY_DOUBLE,
# NPY_LONGDOUBLE, NPY_CFLOAT, NPY_CDOUBLE, NPY_CLONGDOUBLE, NPY_DATETIME,
# NPY_TIMEDELTA, NPY_OBJECT, NPY_STRING, NPY_UNICODE, NPY_VOID
#
mandel = PyUFunc_FromFuncAndData(
loop_func,
elementwise_funcs,
input_output_types,
1, # number of supported input types
TODO, # number of input args
TODO, # number of output args
0, # `identity` element, never mind this
"mandel", # function name
"mandel(z, c) -> computes z*z + c", # docstring
0 # unused
)
Type codes:
NPY_BOOL, NPY_BYTE, NPY_UBYTE, NPY_SHORT, NPY_USHORT, NPY_INT, NPY_UINT,
NPY_LONG, NPY_ULONG, NPY_LONGLONG, NPY_ULONGLONG, NPY_FLOAT, NPY_DOUBLE,
NPY_LONGDOUBLE, NPY_CFLOAT, NPY_CDOUBLE, NPY_CLONGDOUBLE, NPY_DATETIME,
NPY_TIMEDELTA, NPY_OBJECT, NPY_STRING, NPY_UNICODE, NPY_VOID
#
# Some points of note:
(continues on next page)
# Straightforward iteration
for k in range(100):
z = z*z + c
if z.real**2 + z.imag**2 > 1000:
break
import_array()
import_ufunc()
loop_func[0] = PyUFunc_DD_D
input_output_types[0] = NPY_CDOUBLE
input_output_types[1] = NPY_CDOUBLE
(continues on next page)
elementwise_funcs[0] = <void*>mandel_single_point
mandel = PyUFunc_FromFuncAndData(
loop_func,
elementwise_funcs,
input_output_types,
1, # number of supported input types
2, # number of input args
1, # number of output args
0, # `identity` element, never mind this
"mandel", # function name
"mandel(z, c) -> computes iterated z*z + c", # docstring
0 # unused
)
"""
Plot Mandelbrot
================
"""
import numpy as np
import mandel
https://github.com/cython/cython/wiki/MarkLodato-CreatingUfuncs
loop_funcs[0] = PyUFunc_DD_D
input_output_types[0] = NPY_CDOUBLE
input_output_types[1] = NPY_CDOUBLE
input_output_types[2] = NPY_CDOUBLE
elementwise_funcs[0] = <void*>mandel_single_point
loop_funcs[1] = PyUFunc_FF_F
input_output_types[3] = NPY_CFLOAT
input_output_types[4] = NPY_CFLOAT
input_output_types[5] = NPY_CFLOAT
elementwise_funcs[1] = <void*>mandel_single_point_singleprec
mandel = PyUFunc_FromFuncAndData(
loop_func,
elementwise_funcs,
input_output_types,
2, # number of supported input types <----------------
2, # number of input args
1, # number of output args
0, # `identity` element, never mind this
"mandel", # function name
"mandel(z, c) -> computes iterated z*z + c", # docstring
0 # unused
)
ufunc
output = elementwise_function(input)
Both output and input can be a single array element only.
generalized ufunc
output and input can be arrays with a fixed number of dimensions
For example, matrix trace (sum of diag elements):
(n, n) -> ()
Matrix product:
Status in NumPy
• matrix multiplication this way could be useful for operating on many small matrices at once
• Also see tensordot and einsum
int i;
input_1 += steps[0];
input_2 += steps[1];
output += steps[2];
}
}
Suppose you
1. Write a library than handles (multidimensional) binary data,
2. Want to make it easy to manipulate the data with NumPy, or whatever other library,
3. . . . but would not like to have NumPy as a dependency.
Currently, 3 solutions:
1. the “old” buffer interface
2. the array interface
3. the “new” buffer interface (PEP 3118)
Q:
Check what happens if data is now modified, and img saved again.
"""
From buffer
============
Show how to exchange data between numpy and a library that only knows
the buffer interface.
"""
import numpy as np
import Image
#
# Modify the original data, and save again.
#
# It turns out that PIL, which knows next to nothing about NumPy,
# happily shares the same data.
#
x[:, :, 1] = 254
img.save("test2.png")
• Multidimensional buffers
• Data type information present
• NumPy-specific approach; slowly deprecated (but not going away)
• Not integrated in Python otherwise
See also:
Documentation: https://numpy.org/doc/stable/reference/arrays.interface.html
::
Note: .view() has a second meaning: it can make an ndarray an instance of a specialized ndarray
subclass
Masked arrays are arrays that may have missing or invalid entries.
For example, suppose we have an array where the fourth entry is invalid:
>>> mx.mean()
2.75
>>> np.mean(mx)
2.75
Warning: Not all NumPy functions respect masks, for instance np.dot, so check the return types.
>>> mx[1] = 9
>>> x
array([ 1, 9, 3, -99, 5])
The mask
>>> mx[1] = 9
>>> mx
masked_array(data=[1, 9, 3, --, 5],
mask=[False, False, False, True, False],
fill_value=999999)
>>> mx.mask
array([False, False, False, True, False])
The masked entries can be filled with a given value to get an usual array back:
>>> x2 = mx.filled(-1)
>>> x2
array([ 1, 9, 3, -1, 5])
Domain-aware functions
Note: Streamlined and more seamless support for dealing with missing data in arrays is making its
way into NumPy 1.7. Stay tuned!
Canadian rangers were distracted when counting hares and lynxes in 1903-1910 and 1917-1918, and
got the numbers are wrong. (Carrot farmers stayed alert, though.) Compute the mean populations
over time, ignoring the invalid numbers.
>>> data = np.loadtxt('data/populations.txt')
>>> populations = np.ma.masked_array(data[:,1:])
>>> year = data[:, 0]
>>> populations.mean(axis=0)
masked_array(data=[40472.72727272727, 18627.272727272728, 42400.0],
mask=[False, False, False],
fill_value=1e+20)
>>> populations.std(axis=0)
masked_array(data=[21087.656489006717, 15625.799814240254, 3322.5062255844787],
mask=[False, False, False],
fill_value=1e+20)
>>> arr = np.array([('a', 1), ('b', 2)], dtype=[('x', 'S1'), ('y', int)])
>>> arr2 = arr.view(np.recarray)
>>> arr2.x
array([b'a', b'b'], dtype='|S1')
>>> arr2.y
array([1, 2])
8.5 Summary
8.6.1 Why
• “There’s a bug?”
• “I don’t understand what this is supposed to do?”
• “I have this fancy code. Would you like to have it?”
• “I’d like to help! What can I do?”
>>> rng.permutation(12)
array([ 2, 6, 4, 1, 8, 11, 10, 5, 9, 3, 7, 0])
>>> rng.permutation(12.) #doctest: +SKIP
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "_generator.pyx", line 4844, in numpy.random._generator.Generator.permutation
numpy.exceptions.AxisError: axis 0 is out of bounds for array of dimension 0
I'm using NumPy 1.4.1, built from the official tarball, on Windows
64 with Visual studio 2008, on Python.org 64-bit Python.
>>> print(np.__version__)
1...
>>> print(np.__file__)
/...
1. Documentation editor
• https://numpy.org/doc/stable/
• Registration
– Register an account
– Subscribe to scipy-dev mailing list (subscribers-only)
– Problem with mailing lists: you get mail
∗ But: you can turn mail delivery off
∗ “change your subscription options”, at the bottom of
https://mail.python.org/mailman3/lists/scipy-dev.python.org/
– Send a mail @ scipy-dev mailing list; ask for activation:
Hi,
Cheers,
N. N.
Prerequisites
• NumPy
• IPython
• nosetests
• pyflakes
• gdb for the C-debugging part.
Chapter contents
• Avoiding bugs
– Coding best practices to avoid getting in trouble
– pyflakes: fast static analysis
• Debugging workflow
• Using the Python debugger
– Invoking the debugger
349
Scientific Python Lectures, Edition 2024.2rc0.dev0
Brian Kernighan
“Everyone knows that debugging is twice as hard as writing a program in the first place. So if you’re
as clever as you can be when you write it, how will you ever debug it?”
• In TextMate
Menu: TextMate -> Preferences -> Advanced -> Shell variables, add a shell variable:
TM_PYCHECKER = /Library/Frameworks/Python.framework/Versions/Current/bin/pyflakes
autocmd FileType python let &mp = 'echo "*** running % ***" ; pyflakes %'
autocmd FileType tex,mp,rst,python imap <Esc>[15~ <C-O>:make!^M
autocmd FileType tex,mp,rst,python map <Esc>[15~ :make!^M
autocmd FileType tex,mp,rst,python set autowrite
(define-minor-mode pyflakes-mode
"Toggle pyflakes mode.
With no argument, this command toggles the mode.
Non-null prefix argument turns on the mode.
Null prefix argument turns off the mode."
;; The initial value.
nil
;; The indicator for the mode line.
" Pyflakes"
;; The minor mode bindings.
'( ([f5] . pyflakes-thisfile) )
)
• In vim
– Use the pyflakes.vim plugin:
1. download the zip file from https://www.vim.org/scripts/script.php?script_id=2441
2. extract the files in ~/.vim/ftplugin/python
3. make sure your vimrc has filetype plugin indent on
– Alternatively: use the syntastic plugin. This can be configured to use flake8 too and also
handles on-the-fly checking for many other languages.
• In emacs
Use the flymake mode with pyflakes, documented on https://www.emacswiki.org/emacs/FlyMake
and included in Emacs 26 and more recent. To activate it, use M-x (meta-key then x) and enter
flymake-mode at the prompt. To enable it automatically when opening a Python file, add the
following line to your .emacs file:
If you do have a non trivial bug, this is when debugging strategies kick in. There is no silver bullet. Yet,
strategies help:
For debugging a given problem, the favorable situation is when the problem is
isolated in a small number of lines of code, outside framework or application
code, with short modify-run-fail cycles
1. Make it fail reliably. Find a test case that makes the code fail every time.
2. Divide and Conquer. Once you have a failing test case, isolate the failing code.
• Which module.
• Which function.
• Which line of code.
=> isolate a small reproducible failure: a test case
3. Change one thing at a time and re-run the failing test case.
4. Use the debugger to understand what is going wrong.
5. Take notes and be patient. It may take a while.
Note: Once you have gone through this process: isolated a tight piece of code reproducing the bug
and fix the bug using this piece of code, add the corresponding code to your test suite.
The python debugger, pdb: https://docs.python.org/3/library/pdb.html, allows you to inspect your code
interactively.
Specifically it allows you to:
• View the source code.
• Walk up and down the call stack.
• Inspect values of variables.
• Modify values of variables.
• Set breakpoints.
Yes, print statements do work as a debugging tool. However to inspect runtime, it is often more
efficient to use the debugger.
Postmortem
4 def index_error():
5 lst = list("foobar")
----> 6 print(lst[len(lst)])
In [2]: %debug
> /home/jarrod/src/scientific-python-lectures/advanced/debugging/index_error.
˓→py(6)index_error()
4 def index_error():
(continues on next page)
ipdb> list
1 """Small snippet to raise an IndexError."""
2
3
4 def index_error():
5 lst = list("foobar")
----> 6 print(lst[len(lst)])
7
8
9 if __name__ == "__main__":
10 index_error()
ipdb> len(lst)
6
ipdb> print(lst[len(lst) - 1])
r
ipdb> quit
In some situations you cannot use IPython, for instance to debug a script that wants to be called from
the command line. In this case, you can call the script with python -m pdb script.py:
$ python -m pdb index_error.py
> /home/jarrod/src/scientific-python-lectures/advanced/debugging/index_error.py(1)
˓→<module>()
index_error()
File "/home/jarrod/src/scientific-python-lectures/advanced/debugging/index_error.
˓→py", line 6, in index_error
print(lst[len(lst)])
~~~^^^^^^^^^^
IndexError: list index out of range
Uncaught exception. Entering post mortem debugging
Running 'cont' or 'step' will restart the program
> /home/jarrod/src/scientific-python-lectures/advanced/debugging/index_error.
˓→py(6)index_error()
-> print(lst[len(lst)])
(Pdb)
Step-by-step execution
Situation: You believe a bug exists in a module but are not sure where.
For instance we are trying to debug wiener_filtering.py. Indeed the code runs, but the filtering does
not work well.
• Run the script in IPython with the debugger using %run -d wiener_filtering.py :
ipdb> n
> /home/jarrod/src/scientific-python-lectures/advanced/debugging/wiener_
˓→filtering.py(3)<module>()
ipdb> b 29
Breakpoint 1 at /home/jarrod/src/scientific-python-lectures/advanced/debugging/
˓→wiener_filtering.py:29
ipdb> c
> /home/jarrod/src/scientific-python-lectures/advanced/debugging/wiener_
˓→filtering.py(29)iterated_wiener()
• Step into code with n(ext) and s(tep): next jumps to the next statement in the current execution
context, while step will go across execution contexts, i.e. enable exploring inside function calls:
ipdb> s
> /home/jarrod/src/scientific-python-lectures/advanced/debugging/wiener_
˓→filtering.py(30)iterated_wiener()
28 """
1 29 noisy_img = noisy_img
---> 30 denoised_img = local_mean(noisy_img, size=size)
31 l_var = local_var(noisy_img, size=size)
32 for i in range(3):
(continues on next page)
ipdb> n
> /home/jarrod/src/scientific-python-lectures/advanced/debugging/wiener_
˓→filtering.py(31)iterated_wiener()
1 29 noisy_img = noisy_img
30 denoised_img = local_mean(noisy_img, size=size)
---> 31 l_var = local_var(noisy_img, size=size)
32 for i in range(3):
33 res = noisy_img - denoised_img
ipdb> print(l_var)
[[2571 2782 3474 ... 3008 2922 3141]
[2105 708 475 ... 469 354 2884]
[1697 420 645 ... 273 236 2517]
...
[2437 345 432 ... 413 387 4188]
[2598 179 247 ... 367 441 3909]
[2808 2525 3117 ... 4413 4454 4385]]
ipdb> print(l_var.min())
0
Oh dear, nothing but integers, and 0 variation. Here is our bug, we are doing integer arithmetic.
When we run the wiener_filtering.py file, the following warnings are raised:
In [2]: %run wiener_filtering.py
/home/jarrod/src/scientific-python-lectures/advanced/debugging/wiener_filtering.
˓→py:35: RuntimeWarning: divide by zero encountered in divide
We can turn these warnings in exception, which enables us to do post-mortem debugging on them,
and find our problem more quickly:
In [3]: np.seterr(all='raise')
Out[3]: {'divide': 'warn', 'over': 'warn', 'under': 'ignore', 'invalid': 'warn'}
Warning: When running nosetests, the output is captured, and thus it seems that the debugger
does not work. Simply run the nosetests with the -s flag.
• pudb is a good semi-graphical debugger with a text user interface in the console.
• The Visual Studio Code integrated development environment includes a debugging mode.
• The Mu editor is a simple Python editor that includes a debugging mode.
ipdb> help
Undocumented commands:
======================
interact
If you have a segmentation fault, you cannot debug it with pdb, as it crashes the Python interpreter
before it can drop in the debugger. Similarly, if you have a bug in C code embedded in Python, pdb is
useless. For this we turn to the gnu debugger, gdb, available on Linux.
Before we start with gdb, let us add a few Python-specific tools to it. For this we add a few macros to
our ~/.gdbinit. The optimal choice of macro depends on your Python version and your gdb version. I
have added a simplified version in gdbinit, but feel free to read DebuggingWithGdb.
To debug with gdb the Python script segfault.py, we can run the script in gdb as follows
$ gdb python
...
(gdb) run segfault.py
Starting program: /usr/bin/python segfault.py
[Thread debugging using libthread_db enabled]
We get a segfault, and gdb captures it for post-mortem debugging in the C level stack (not the Python
call stack). We can debug the C call stack using gdb’s commands:
(gdb) up
#1 0x004af4f5 in _copy_from_same_shape (dest=<value optimized out>,
src=<value optimized out>, myfunc=0x496780 <_strided_byte_copy>,
swap=0)
at numpy/core/src/multiarray/ctors.c:748
748 myfunc(dit->dataptr, dest->strides[maxaxis],
As you can see, right now, we are in the C code of numpy. We would like to know what is the Python
code that triggers this segfault, so we go up the stack until we hit the Python execution loop:
(gdb) up
#8 0x080ddd23 in call_function (f=
Frame 0x85371ec, for file /home/varoquau/usr/lib/python2.6/site-packages/numpy/
˓→core/arrayprint.py, line 156, in _leading_trailing (a=<numpy.ndarray at remote␣
at ../Python/ceval.c:3750
3750 ../Python/ceval.c: No such file or directory.
in ../Python/ceval.c
(gdb) up
#9 PyEval_EvalFrameEx (f=
Frame 0x85371ec, for file /home/varoquau/usr/lib/python2.6/site-packages/numpy/
˓→core/arrayprint.py, line 156, in _leading_trailing (a=<numpy.ndarray at remote␣
at ../Python/ceval.c:2412
2412 in ../Python/ceval.c
(gdb)
Once we are in the Python execution loop, we can use our special Python helper function. For instance
we can find the corresponding Python code:
(gdb) pyframe
/home/varoquau/usr/lib/python2.6/site-packages/numpy/core/arrayprint.py (158): _
˓→leading_trailing
(gdb)
This is numpy code, we need to go up until we find code that we have written:
(gdb) up
...
(gdb) up
#34 0x080dc97a in PyEval_EvalFrameEx (f=
Frame 0x82f064c, for file segfault.py, line 11, in print_big_array (small_array=
˓→<numpy.ndarray at remote 0x853ecf0>, big_array=<numpy.ndarray at remote 0x853ed20>),
˓→ throwflag=0) at ../Python/ceval.c:1630
def make_big_array(small_array):
big_array = stride_tricks.as_strided(
small_array, shape=(int(2e6), int(2e6)), strides=(32, 32)
)
return big_array
Thus the segfault happens when printing big_array[-10:]. The reason is simply that big_array has
been allocated with its end outside the program memory.
Note: For a list of Python-specific commands defined in the gdbinit, read the source of this file.
Wrap up exercise
The following script is well documented and hopefully legible. It seeks to answer a problem of actual
interest for numerical computing, but it does not work. . . Can you debug it?
Python source code: to_debug.py
Donald Knuth
Prerequisites
• line_profiler
Chapters contents
• Optimization workflow
• Profiling Python code
– Timeit
– Profiler
– Line-profiler
• Making code go faster
– Algorithmic optimization
∗ Example of the SVD
361
Scientific Python Lectures, Edition 2024.2rc0.dev0
10.2.1 Timeit
In [2]: a = np.arange(1000)
In [3]: %timeit a ** 2
867 ns +- 2.86 ns per loop (mean +- std. dev. of 7 runs, 1,000,000 loops each)
In [5]: %timeit a * a
978 ns +- 13.9 ns per loop (mean +- std. dev. of 7 runs, 1,000,000 loops each)
Note: For long running calls, using %time instead of %timeit; it is less precise but faster
10.2.2 Profiler
Useful when you have a large program to profile, for example the following file:
# For this example to run, you also need the 'ica.py' file
import numpy as np
import scipy as sp
if __name__ == "__main__":
test()
Note: This is a combination of two unsupervised learning techniques, principal component analysis
(PCA) and independent component analysis (ICA). PCA is a technique for dimensionality reduction,
i.e. an algorithm to explain the observed variance in your data using less dimensions. ICA is a source
separation technique, for example to unmix multiple signals that have been recorded through multiple
sensors. Doing a PCA first and then an ICA can be useful if you have more sensors than signals. For
more information see: the FastICA example from scikits-learn.
To run it, you also need to download the ica module. In IPython we can time the script:
Clearly the svd (in decomp.py) is what takes most of our time, a.k.a. the bottleneck. We have to find a
way to make this step go faster, or to avoid this step (algorithmic optimization). Spending time on the
rest of the code is useless.
Similar profiling can be done outside of IPython, simply calling the built-in Python profilers cProfile
and profile.
$ python -m cProfile -o demo.prof demo.py
Using the -o switch will output the profiler results to the file demo.prof to view with an external tool.
This can be useful if you wish to process the profiler output with a visualization tool.
10.2.3 Line-profiler
The profiler tells us which function takes most of the time, but not where it is called.
For this, we use the line_profiler: in the source file, we decorate a few functions that we want to inspect
with @profile (no need to import it)
@profile
def test():
rng = np.random.default_rng()
data = rng.random((5000, 100))
u, s, v = linalg.svd(data)
pca = u[:, :10] @ data
results = fastica(pca.T, whiten=False)
Then we run the script using the kernprof command, with switches -l, --line-by-line and -v,
--view to use the line-by-line profiler and view the results in addition to saving them:
$ kernprof -l -v demo.py
The SVD is taking all the time. We need to optimise this line.
Once we have identified the bottlenecks, we need to make the corresponding code go faster.
The first thing to look for is algorithmic optimization: are there ways to compute less, or better?
For a high-level view of the problem, a good understanding of the maths behind the algorithm helps.
However, it is not uncommon to find simple changes, like moving computation or memory allocation
outside a for loop, that bring in big gains.
In both examples above, the SVD - Singular Value Decomposition - is what takes most of the time.
Indeed, the computational cost of this algorithm is roughly 𝑛3 in the size of the input matrix.
However, in both of these example, we are not using all the output of the SVD, but only the first few rows
of its first return argument. If we use the svd implementation of SciPy, we can ask for an incomplete
version of the SVD. Note that implementations of linear algebra in SciPy are richer then those in NumPy
and should be preferred.
def test():
rng = np.random.default_rng()
(continues on next page)
Real incomplete SVDs, e.g. computing only the first 10 eigenvectors, can be computed with arpack,
available in scipy.sparse.linalg.eigsh.
For certain algorithms, many of the bottlenecks will be linear algebra computations. In this case,
using the right function to solve the right problem is key. For instance, an eigenvalue problem with
a symmetric matrix is easier to solve than with a general matrix. Also, most often, you can avoid
inverting a matrix and use a less costly (and more numerically stable) operation.
Know your computational linear algebra. When in doubt, explore scipy.linalg, and use %timeit to
try out different alternatives on your data.
A complete discussion on advanced use of NumPy is found in chapter Advanced NumPy, or in the article
The NumPy array: a structure for efficient numerical computation by van der Walt et al. Here we discuss
only some commonly encountered tricks to make code faster.
• Vectorizing for loops
Find tricks to avoid for loops using NumPy arrays. For this, masks and indices arrays can be
useful.
• Broadcasting
Use broadcasting to do operations on arrays as small as possible before combining them.
• In place operations
In [18]: a = np.zeros(1e7)
note: we need global a in the timeit so that it work, as it is assigning to a, and thus considers it
as a local variable.
• Be easy on the memory: use views, and not copies
Copying big arrays is as costly as making simple numerical operations on them:
In [21]: a = np.zeros(1e7)
In [23]: %timeit a + 1
10 loops, best of 3: 112 ms per loop
In [27]: c.strides
Out[27]: (80000, 8)
This is the reason why Fortran ordering or C ordering may make a big difference on operations:
In [32]: c = np.ascontiguousarray(a.T)
In [33]: %timeit b @ c
8.04 ms +- 23.8 us per loop (mean +- std. dev. of 7 runs, 100 loops each)
Note that copying the data to work around this effect may not be worth it:
Using numexpr can be useful to automatically optimize code for such effects.
• Use compiled code
The last resort, once you are sure that all the high-level optimizations have been explored, is to
transfer the hot spots, i.e. the few lines or functions in which most of the time is spent, to compiled
code. For compiled code, the preferred option is to use Cython: it is easy to transform exiting
Python code in compiled code, and with a good use of the NumPy support yields efficient code on
NumPy arrays, for instance by unrolling loops.
Warning: For all the above: profile and time your choices. Don’t base your optimization on
theoretical considerations.
• If you need to profile memory usage, you could try the memory_profiler
• If you need to profile down into C extensions, you could try using gperftools from Python with
yep.
• If you would like to track performance of your code across time, i.e. as you make new commits to
your repository, you could try: asv
• If you need some interactive visualization why not try RunSnakeRun
11.1 Introduction
369
Scientific Python Lectures, Edition 2024.2rc0.dev0
11.1.4 Prerequisites
• numpy
• scipy
• matplotlib (optional)
• ipython (the enhancements come handy)
Examples
offset: row
2: 9
1: --10------
0: 1 . 11 .
-1: 5 2 . 12
-2: . 6 3 .
-3: . . 7 4
---------8
• matrix-vector multiplication
Examples