Scientific Computing Python 3
User Manual: scientific-computing-python-3
Open the PDF directly: View PDF 
.
Page Count: 322
| Download | |
| Open PDF In Browser | View PDF |
Scientific Computing with
Python 3
An example-rich, comprehensive guide for all of your Python
computational needs
Claus Führer
Jan Erik Solem
Olivier Verdier
BIRMINGHAM - MUMBAI
Scientific Computing with Python 3
Copyright © 2016 Packt Publishing
All rights reserved. No part of this book may be reproduced, stored in a retrieval system, or
transmitted in any form or by any means, without the prior written permission of the
publisher, except in the case of brief quotations embedded in critical articles or reviews.
Every effort has been made in the preparation of this book to ensure the accuracy of the
information presented. However, the information contained in this book is sold without
warranty, either express or implied. Neither the authors, nor Packt Publishing, and its
dealers and distributors will be held liable for any damages caused or alleged to be caused
directly or indirectly by this book.
Packt Publishing has endeavored to provide trademark information about all of the
companies and products mentioned in this book by the appropriate use of capitals.
However, Packt Publishing cannot guarantee the accuracy of this information.
First published: December 2016
Production reference: 1141216
Published by Packt Publishing Ltd.
Livery Place
35 Livery Street
Birmingham
B3 2PB, UK.
ISBN 978-1-78646-351-7
www.packtpub.com
Credits
Authors
Copy Editor
Claus Führer
Vikrant Phadkay
Jan Erik Solem
Sneha Singh
Olivier Verdier
Reviewers
Project Coordinator
Helmut Podhaisky
Nidhi Joshi
Commissioning Editor
Proofreader
Veena Pagare
Safis Editing
Acquisition Editor
Indexer
Sonali Vernekar
Mariammal Chettiyar
Content Development Editor
Graphics
Aishwarya Pandere
Disha Haria
Technical Editor
Production Coordinator
Karan Thakkar
Arvindkumar Gupta
About the Authors
Claus Führer is a professor of scientific computations at Lund University, Sweden. He has
an extensive teaching record that includes intensive programming courses in numerical
analysis and engineering mathematics across various levels in many different countries and
teaching environments. Claus also develops numerical software in research collaboration
with industry and received Lund University’s Faculty of Engineering Best Teacher Award
in 2016.
Jan Erik Solem is a Python enthusiast, former associate professor, and currently the CEO of
Mapillary, a street imagery computer vision company. He has previously worked as a face
recognition expert, founder and CTO of Polar Rose, and computer vision team leader at
Apple. Jan is a World Economic Forum technology pioneer and won the Best Nordic Thesis
Award 2005-2006 for his dissertation on image analysis and pattern recognition. He is also
the author of "Programming Computer Vision with Python" (O'Reilly 2012).
Olivier Verdier began using Python for scientific computing back in 2007 and received a
PhD in mathematics from Lund University in 2009. He has held post-doctoral positions in
Cologne, Trondheim, Bergen, and Umeå and is now an associate professor of mathematics
at Bergen University College, Norway.
About the Reviewer
Helmut Podhaisky works in the Institute of Mathematics at the Martin Luther University in
Halle-Wittenberg, where he teaches mathematics and scientific computing. He has coauthored a book on numerical methods for ordinary differential equations as well as several
research papers on numerical methods. For work and fun, he uses Python, Fortran, Octave,
Mathematica, and Haskell.
www.PacktPub.com
For support files and downloads related to your book, please visit www.PacktPub.com.
Did you know that Packt offers eBook versions of every book published, with PDF and
ePub files available? You can upgrade to the eBook version at www.PacktPub.com and as a
print book customer, you are entitled to a discount on the eBook copy. Get in touch with us
at service@packtpub.com for more details.
At www.PacktPub.com, you can also read a collection of free technical articles, sign up for a
range of free newsletters and receive exclusive discounts and offers on Packt books and
eBooks.
https://www.packtpub.com/mapt
Get the most in-demand software skills with Mapt. Mapt gives you full access to all Packt
books and video courses, as well as industry-leading tools to help you plan your personal
development and advance your career.
Why subscribe?
Fully searchable across every book published by Packt
Copy and paste, print, and bookmark content
On demand and accessible via a web browser
Acknowledgement
We want to acknowledge the competent and helpful comments and suggestions by Helmut
Podhaisky, Halle University, Germany. To have such a partner in the process of writing a
book is a big luck and chance for the authors.
We would also like to express our gratitude towards the reviewers of the first edition of this
book, [7], Linda Kann, KTH Stockholm, Hans Petter Langtangen, Simula Research
Laboratory, and Alf Inge Wang, NTNU Trondheim.
A book has to be tested in teaching. And here we had fantastic partners: the teaching
assistants from the course "Beräkningsprogramering med Python" during the years and the
colleagues involved in teaching: Najmeh Abiri, Christian Andersson, Dara Maghdid, Peter
Meisrimel, Fatemeh Mohammadi, Azahar Monge, Anna-Maria Persson, Alexandros
Sopasakis, Tony Stillfjord, Lund University. Najmeh Abiri also tested most of the Jupyter
notebook material which you find on the book's webpage.
A book has not only to be written, it has to be published, and in this process Aishwarya
Pandere and Karan Thakkar, PACKT Publishing, were always constructive, friendly and
helpful partners bridging different time zones and different text processing tools. Thanks.
Claus Führer, Jan-Erik Solem, Olivier Verdier Lund, Bergen 2016
Table of Contents
Preface
Chapter 1: Getting Started
Installation and configuration instructions
Installation
Anaconda
Configuration
Python Shell
Executing scripts
Getting Help
Jupyter – Python notebook
Program and program flow
Comments
Line joining
Basic types
Numbers
Strings
Variables
Lists
Operations on lists
Boolean expressions
Repeating statements with loops
Repeating a task
Break and else
Conditional statements
Encapsulating code with functions
Scripts and modules
Simple modules – collecting functions
Using modules and namespaces
Interpreter
Summary
Chapter 2: Variables and Basic Types
Variables
Numeric types
Integers
1
10
11
11
12
13
13
14
14
14
15
16
16
17
17
17
18
18
19
19
20
20
21
21
22
23
23
24
24
25
26
26
27
28
Plain integers
Floating point numbers
Floating point representation
Infinite and not a number
Underflow – Machine Epsilon
Other float types in NumPy
Complex numbers
Complex Numbers in Mathematics
The j notation
Real and imaginary parts
Booleans
Boolean operators
Boolean casting
Automatic Boolean casting
Return values of and and or
Boolean and integer
Strings
Operations on strings and string methods
String formatting
Summary
Exercises
Chapter 3: Container Types
28
29
29
30
31
32
33
33
33
34
36
36
37
38
38
39
39
41
42
43
43
47
Lists
Slicing
Strides
Altering lists
Belonging to a list
List methods
In–place operations
Merging lists – zip
List comprehension
Arrays
Tuples
Dictionaries
Creating and altering dictionaries
Looping over dictionaries
Sets
Container conversions
Type checking
Summary
[ ii ]
47
48
50
51
51
52
52
53
53
54
56
56
57
58
58
60
61
62
Exercises
62
Chapter 4: Linear Algebra – Arrays
65
Overview of the array type
Vectors and matrices
Indexing and slices
Linear algebra operations
Solving a linear system
Mathematical preliminaries
Arrays as functions
Operations are elementwise
Shape and number of dimensions
The dot operations
The array type
Array properties
Creating arrays from lists
Accessing array entries
Basic array slicing
Altering an array using slices
Functions to construct arrays
Accessing and changing the shape
The shape function
Number of dimensions
Reshape
Transpose
Stacking
Stacking vectors
Functions acting on arrays
Universal functions
Built-in universal functions
Create universal functions
Array functions
Linear algebra methods in SciPy
Solving several linear equation systems with LU
Solving a least square problem with SVD
More methods
Summary
Exercises
Chapter 5: Advanced Array Concepts
Array views and copies
Array views
65
65
67
67
68
69
69
69
70
71
73
73
74
75
75
77
77
78
78
79
79
81
82
82
83
83
83
84
86
87
87
89
90
91
91
94
94
94
[ iii ]
Slices as views
Transpose and reshape as views
Array copy
Comparing arrays
Boolean arrays
Checking for equality
Boolean operations on arrays
Array indexing
Indexing with Boolean arrays
Using where
Performance and Vectorization
Vectorization
Broadcasting
Mathematical view
Constant functions
Functions of several variables
General mechanism
Conventions
Broadcasting arrays
The broadcasting problem
Shape mismatch
Typical examples
Rescale rows
Rescale columns
Functions of two variables
Sparse matrices
Sparse matrix formats
Compressed sparse row
Compressed Sparse Column
Row-based linked list format
Altering and slicing matrices in LIL format
Generating sparse matrices
Sparse matrix methods
Summary
Chapter 6: Plotting
95
95
96
96
96
97
98
99
99
100
101
102
103
103
104
105
105
107
107
107
109
110
110
110
110
112
113
113
115
115
116
116
117
118
119
Basic plotting
Formatting
Meshgrid and contours
Images and contours
Matplotlib objects
The axes object
119
125
128
132
134
135
[ iv ]
Modifying line properties
Annotations
Filling areas between curves
Ticks and ticklabels
Making 3D plots
Making movies from plots
Summary
Exercises
136
137
138
140
141
145
146
147
Chapter 7: Functions
149
Basics
Parameters and arguments
Passing arguments – by position and by keyword
Changing arguments
Access to variables defined outside the local namespace
Default arguments
Beware of mutable default arguments
Variable number of arguments
Return values
Recursive functions
Function documentation
Functions are objects
Partial application
Using Closures
Anonymous functions – the lambda keyword
The lambda construction is always replaceable
Functions as decorators
Summary
Exercises
Chapter 8: Classes
149
150
150
151
152
153
154
154
156
157
159
160
160
161
161
162
163
164
165
167
Introduction to classes
Class syntax
The __init__ method
Attributes and methods
Special methods
Reverse operations
Attributes that depend on each other
The property function
Bound and unbound methods
Class attributes
[v]
168
169
169
170
172
174
176
177
178
178
Class methods
Subclassing and inheritance
Encapsulation
Classes as decorators
Summary
Exercises
179
181
184
185
188
188
Chapter 9: Iterating
190
The for statement
Controlling the flow inside the loop
Iterators
Generators
Iterators are disposable
Iterator tools
Generators of recursive sequences
Arithmetic geometric mean
Convergence acceleration
List filling patterns
List filling with the append method
List from iterators
Storing generated values
When iterators behave as lists
Generator expression
Zipping iterators
Iterator objects
Infinite iterations
The while loop
Recursion
Summary
Exercises
Chapter 10: Error Handling
190
191
192
193
194
194
196
196
198
200
200
200
201
202
202
203
204
205
205
206
207
207
210
What are exceptions?
Basic principles
Raising exceptions
Catching exceptions
User-defined exceptions
Context managers — the with statement
Finding Errors: Debugging
Bugs
The stack
[ vi ]
210
212
212
213
215
216
218
218
218
The Python debugger
Overview – debug commands
Debugging in IPython
Summary
219
221
222
223
Chapter 11: Namespaces, Scopes, and Modules
Namespace
Scope of a variable
Modules
Introduction
Modules in IPython
224
224
225
227
227
229
229
229
230
230
The IPython magic command
The variable __name__
Some useful modules
Summary
Chapter 12: Input and Output
231
File handling
Interacting with files
Files are iterable
File modes
NumPy methods
savetxt
loadtxt
Pickling
Shelves
Reading and writing Matlab data files
Reading and writing images
Summary
Chapter 13: Testing
231
231
233
233
234
234
234
235
236
237
237
238
239
Manual testing
Automatic testing
Testing the bisection algorithm
Using unittest package
Test setUp and tearDown methods
Parameterizing tests
Assertion tools
Float comparisons
Unit and functional tests
Debugging
[ vii ]
239
240
241
243
244
245
247
247
249
250
Test discovery
Measuring execution time
Timing with a magic function
Timing with the Python module timeit
Timing with a context manager
Summary
Exercises
Chapter 14: Comprehensive Examples
Polynomials
Theoretical background
Tasks
The polynomial class
Newton polynomial
Spectral clustering
Solving initial value problems
Summary
Exercises
250
250
251
252
253
254
254
256
256
256
258
259
263
265
269
273
273
Chapter 15: Symbolic Computations - SymPy
274
What are symbolic computations?
Elaborating an example in SymPy
Basic elements of SymPy
Symbols – the basis of all formulas
Numbers
Functions
274
276
278
278
279
279
280
281
282
283
284
285
287
289
290
292
292
294
Undefined functions
Elementary Functions
Lambda – functions
Symbolic Linear Algebra
Symbolic matrices
Examples for Linear Algebra Methods in SymPy
Substitutions
Evaluating symbolic expressions
Example: A study on the convergence order of Newton's Method
Converting a symbolic expression into a numeric function
A study on the parameter dependency of polynomial coefficients
Summary
Appendix: References
Index
295
298
[ viii ]
Preface
Python can be used for more than just general-purpose programming. It is a free, open
source language and environment that has tremendous potential for use within the domain
of scientific computing. This book presents Python in tight connection with mathematical
applications and demonstrates how to use various concepts in Python for computing
purposes, including examples with the latest version of Python 3. Python is an effective tool
to use when coupling scientific computing and mathematics and this book will teach you
how to use it for linear algebra, arrays, plotting, iterating, functions, polynomials, and much
more.
What this book covers
Chapter 1, Getting Started, addresses the main language elements of Python without going
into detail. Here we make a brief tour through all. It is a good starting point for those who
want to start directly. It is a quick reference for those readers who want in a later chapter
understand an example which uses might use constructs like functions before functions
were explained in deep .
Chapter 2, Variables and Basic Types, presents the most important and basic types in Python.
Float is the more important datatype in scientific computing together with the special
numbers nan and inf. Booleans, integers, complex, and strings are other basic datatypes,
which will be used throughout this book.
Chapter 3, Container Types, explains how to work with container types, mainly lists.
Dictionaries and tuples will be explained as well as indexing and looping, through
container objects. Occasionally, one uses even sets as a special container type.
Chapter 4, Linear Algebra, works with the most important objects in linear algebra--vectors
and matrices. This book chooses NumPy array as the central tool for describing matrices
and even higher order tensors. Arrays have many advanced features and allows also for
universal functions acting on matrices or vectors elementwise. The book emphasizes on
array indexing, slices, and the dot product as the basic operation in most computing tasks.
Some linear algebra examples are worked out to demonstrate the use of SciPy's submodule
linalg.
Preface
Chapter 5, Advanced Array Concepts, explains some more advanced aspects of arrays. The
difference between array copies and views is explained extensively as views make
programs using arrays very fast but are often a source for errors, which are hard to debug.
The use of Boolean arrays to write effective, compact, and readable code is shown and
demonstrated. Finally, the technique of array broadcasting-- a unique feature of NumPy
arrays -- is explained by its analogy to operations performed on functions.
Chapter 6, Plotting, shows how to make plots, mainly classical x/yplots but also 3D plots
and histograms. Scientific computing requires good tools for visualizing the results.
Python's module matplotlib is introduced starting from the handy plotting commands in its
submodule pyplot. Finetuning and modifying plots becomes possible by creating graphical
objects such as axes. We show how attributes of these objects can be changed and
annotations can be made.
Chapter 7, Functions, form the fundamental building block in programming, which is
probably nearest to underlying mathematical concepts. Function definition and function
calls are explained as the different ways to set function arguments. Anonymous lambda
functions are introduced and used in various examples throughout the book.
Chapter 8, Classes, defines objects as instances of classes, which we provide with methods
and attributes. In mathematics, class attributes often depend on each other, which requires
special programming techniques for setter and getter functions. Basic mathematical
operations such as + can be defined for special mathematic datatypes. Inheritance and
abstraction are mathematical concepts which are reflected by object oriented programming.
We demonstrate the use of inheritance by a simple solver class for ordinary differential
equations.
Chapter 9, Iterating, presents iteration using loops and iterators. There is now a chapter in
this book without loops and iterations, but here we come to principles of iterators and
create own generator objects. In this chapter, you learn why a generator can be exhausted
and how infinite loops can be programmed. Python's module itertools is a useful
companion for this chapter.
Chapter 10, Error Handling, covers errors and exceptions and how to find and fix them. An
error or an exception is an event, which breaks the execution of a program unit. This
chapter shows what to do then, that is, how an exception can be handled. You learn to
define your own exception classes and how to provide valuable information, which can be
used for catching these exceptions. Error handling is more than printing an error message.
[2]
Preface
Chapter 11, Namespaces, Scopes and Modules, covers Python modules. What are local and
global variables? When is a variable known and when is it unknown to a program unit?
This is discussed in this chapter. A variable can be passed to a function by a parameter list
or tacitly injected by making use of its scope. When should this technique be applied and
when not? This chapter tries to give an answer to this central question.
Chapter 12, Input and Output, covers some options for handling data files. Data files are
used for storing and providing data for a given problem, often large scale measurements.
This chapter describes how this data can be accessed and modified using different formats.
Chapter 13, Testing, focuses on testing for scientific programming. The key tool is unittest,
which allows for automatic testing and parametrized tests. By considering the classical
bisection algorithm in numerical mathematics, we exemplify different steps to design
meaningful tests, which as a side effect also deliver a documentation of the use of a piece of
code. Careful testing provides test protocols which can be later helpful when debugging a
complex code often written by many different programmers.
Chapter 14, Comprehensive Examples, presents some comprehensive and longer examples
together with a brief introduction to the theoretical background and their complete
implementation. These examples make use of all constructs shown in the book so far and
put them in a larger and more complex context. They are open for extensions by the reader.
Chapter 15, Symbolic Computations - SymPy, speaks about symbolic computations. Scientific
computing is mainly numeric computations with inexact data and approximative results.
This is contrasted by symbolic computations often formal manipulation, which aims for
exact solutions in a closed form expression. In this last chapter of the book, we introduce
this technique in Python, which is often used for deriving and verifying theoretically
mathematical models and numerical results. We emphasize on high precision floating point
evaluation of symbolic expressions.
What you need for this book
You would need Pyhon3.5 or higher, SciPy, NumPy, Matplotlib, IPython shell (we
recommend strongly to install Python and its packages through Anaconda). The examples
of the book do not have any special hardware requirements on memory and graphics.
[3]
Preface
Who this book is for
This book is the outcome of a course on Python for scientific computing which is taught at
Lund University since 2008. The course expanded over the years, and condensed versions
of the material were taught at universities in Cologne, Trondheim, Stavanger, Soran,
Lappeenranta and also in computation oriented companies.
Our belief is that Python and its surrounding scientific computing ecosystem — SciPy,
NumPY and matplotlib — represent a tremendous progress in scientific computing
environment. Python and the aforementioned libraries are free and open source. What’s
more, is a modern language featuring all the bells and whistles that this adjective entails:
object oriented programming, testing, advanced shell with IPython, etc. When writing this
book we had two groups of readers in mind:
The reader who chooses Python as his or her first programming language will
use this book in a teacher-led course. The book guides into the different topics
and offers background reading and experimenting. A teacher typically selects
and orders the material from this book in such a way, that it fits to the specific
learning outcomes of an introductory course.
The reader who already has some experience in programming, and some taste for
scientific computing or mathematics will use this book as a companion when
diving into the world of Scipy and Numpy. Programming in Python can be quite
different from programming in MATLAB, say. The book wants to point out the
"pythonic" way of programming, which makes programming a pleasure.
Our goal is to explain the steps to get started with Python in the context of scientific
computing. The book may be read either from the first page to the last, or by picking the
bits that seem most interesting. Needless to say, as improving one’s programming skills
requires considerable practice, it is highly advisable to experiment and play with the
examples and the exercises in the book.
We hope that the readers will enjoy programming with Python, SciPy, NumPY and
matplotlib as much as we do.
[4]
Preface
Python vs Other Languages
When it comes to deciding what language to use for a book on scientific computing many
factors come in to play. The learning threshold of the language itself is important for
newcomers, here scripting languages usually provide the best options. A wide range of
modules for numerical computing is necessary, preferably with a strong developer
community. If these core modules are built on a well-tested, optimized foundation of fast
libraries like e.g. LAPACK, even better. Finally, if the language is also usable in a wider
setting and a wider range of applications, the chance of the reader using the skills learned
from this book outside an academic setting is greater. Therefore the choice of Python was a
natural one.
In short, Python is
free and open source
a scripting language, meaning that it is interpreted
a modern language (object oriented, exception handling, dynamic typing etc.)
concise, easy to read and quick to learn
full of freely available libraries, in particular scientific ones (linear algebra,
visualization tools, plotting, image analysis, differential equations solving,
symbolic computations, statistics etc.)
useful in a wider setting: scientific computing, scripting, web sites, text parsing,
etc.
widely used in industrial applications
There are other alternatives to Python. Some of them and the differences to Python are
listed here.
Java, C++ : Object oriented, compiled languages. More verbose and low level compared to
Python. Few scientific libraries.
C, FORTRAN : Low level compiled languages. Both languages are extensively used in
scientific computing, where computational time matters. Nowadays these languages are
often combined with Python wrappers.
PHP, Ruby, other interpreted languages. PHP is web oriented. Ruby is as flexible as Python
but has few scientific libraries.
[5]
Preface
MATLAB, Scilab, Octave : MATLAB is a tool for matrix computation that evolved for
scientific computing. The scientific library is huge. The language features are not as
developed as those of Python. Neither free nor open source. SciLab and Octave are open
source tools which are syntactically similar to MATLAB.
Haskell : Haskell is a modern functional language and follows different programming
paradigms than Python. There are some common constructions like list comprehension.
Haskell is rarely used in scientific computing. See also [12].
Other Python literature
Here we give some hints to literature on Python which can serve as complementary sources
or as texts for parallel reading. Most introductory books on Python are devoted to teach this
language as a general purpose tool. One excellent example which we want to mention here
explicitly is [19]. It explains the language by simple examples, e.g. object oriented
programming is explained by organizing a pizza bakery.
There are very few books dedicated to Python directed towards scientific computing and
engineering. Among these few books we would like to mention the two books by
Langtangen which combine scientific computing with the modern "pythonic" view on
programming, [16,17].
This "pythonic" view is also the guiding line of our way of teaching programming of
numerical algorithms. We try to show how many well-established concepts and
constructions in computer science can be applied to problems within scientific computing.
The pizza-bakery example is replaced by Lagrange polynomials, generators become time
stepping methods for ODEs, and so on.
Finally we have to mention the nearly infinite amount of literature on the web. The web
was also a big source of knowledge when preparing this book. Literature from the web
often covers things that are new, but can also be totally outdated. The web also presents
solutions and interpretations which might contradict each other. We strongly recommend
to use the web as additional source, but we consider a "traditional" textbook with the web
resources "edited" as the better entry point to a rich new world.
[6]
Preface
Conventions
In this book, you will find a number of text styles that distinguish between different kinds
of information. Here are some examples of these styles and an explanation of their meaning.
Code words in text, database table names, folder names, filenames, file extensions,
pathnames, and user input are shown as follows: "install additional packages with conda
install within your virtual environment"
A block of code is set as follows:
from scipy import *
from matplotlib.pyplot import *
Any command-line input or output is written as follows:
jupyter notebook
New terms and important words are shown in bold. Words that you see on the screen, for
example, in menus or dialog boxes, appear in the text like this: "The Jupyter notebook is a
fantastic tool for demonstrating your work."
Warnings or important notes appear in a box like this.
Tips and tricks appear like this.
Reader feedback
Feedback from our readers is always welcome. Let us know what you think about this
book-what you liked or disliked. Reader feedback is important for us as it helps us develop
titles that you will really get the most out of. To send us general feedback, simply email feedback@packtpub.com, and mention the book's title in the subject of your
message. If there is a topic that you have expertise in and you are interested in either
writing or contributing to a book, see our author guide at www.packtpub.com/authors.
[7]
Preface
Customer support
Now that you are the proud owner of a Packt book, we have a number of things to help you
to get the most from your purchase.
Downloading the example code
You can download the example code files for this book from your account at http://www.p
acktpub.com. If you purchased this book elsewhere, you can visit http://www.packtpub.c
om/supportand register to have the files e-mailed directly to you.
You can download the code files by following these steps:
1.
2.
3.
4.
5.
6.
7.
Log in or register to our website using your e-mail address and password.
Hover the mouse pointer on the SUPPORT tab at the top.
Click on Code Downloads & Errata.
Enter the name of the book in the Search box.
Select the book for which you're looking to download the code files.
Choose from the drop-down menu where you purchased this book from.
Click on Code Download.
Once the file is downloaded, please make sure that you unzip or extract the folder using the
latest version of:
WinRAR / 7-Zip for Windows
Zipeg / iZip / UnRarX for Mac
7-Zip / PeaZip for Linux
The code bundle for the book is also hosted on GitHub at https://github.com/PacktPubl
ishing/Scientific-Computing-with-Python-3. We also have other code bundles from
our rich catalog of books and videos available at https://github.com/PacktPublishing/.
Check them out!
Downloading the color images of this book
We also provide you with a PDF file that has color images of the screenshots/diagrams used
in this book. The color images will help you better understand the changes in the output.
You can download this file from https://www.packtpub.com/sites/default/files/down
loads/ScientificComputingwithPython3_ColorImages.pdf.
[8]
Preface
Errata
Although we have taken every care to ensure the accuracy of our content, mistakes do
happen. If you find a mistake in one of our books-maybe a mistake in the text or the codewe would be grateful if you could report this to us. By doing so, you can save other readers
from frustration and help us improve subsequent versions of this book. If you find any
errata, please report them by visiting http://www.packtpub.com/submit-errata, selecting
your book, clicking on the Errata Submission Form link, and entering the details of your
errata. Once your errata are verified, your submission will be accepted and the errata will
be uploaded to our website or added to any list of existing errata under the Errata section of
that title.
To view the previously submitted errata, go to https://www.packtpub.com/books/conten
t/supportand enter the name of the book in the search field. The required information will
appear under the Errata section.
Piracy
Piracy of copyrighted material on the Internet is an ongoing problem across all media. At
Packt, we take the protection of our copyright and licenses very seriously. If you come
across any illegal copies of our works in any form on the Internet, please provide us with
the location address or website name immediately so that we can pursue a remedy.
Please contact us at copyright@packtpub.com with a link to the suspected pirated
material.
We appreciate your help in protecting our authors and our ability to bring you valuable
content.
Questions
If you have a problem with any aspect of this book, you can contact us
at questions@packtpub.com, and we will do our best to address the problem.
[9]
1
Getting Started
In this chapter, we will give a brief overview of the principal syntactical elements of Python.
Readers who have just started learning programming are guided through the book in this
chapter. Every topic is presented here in a how-to way and will be explained later in the
book in a deeper conceptual manner and will also be enriched with many applications and
extensions.
Readers who are already familiar with another programming language will come across, in
this chapter, the Python way of doing classical language constructs. It offers them a quick
start to Python programming.
Both types of readers are encouraged to take this chapter as a brief guideline when
zigzagging through the book. However, before we start we have to make sure that
everything is in place and you have the correct version of Python installed together with the
main modules for Scientific Computing and tools, such as a good editor and a shell, which
helps in code developing and testing.
Read the following section, even if you already have access to a computer with Python
installed. You might want to adjust things to have a working environment conforming to
the presentation in this book.
Getting Started
Installation and configuration instructions
Before diving into the subject of the book you should have all the relevant tools installed on
your computer. We will give you some advice and recommend tools that you might want to
use. We only describe public domain and free tools.
Installation
There are currently two major versions of Python; the 2.x branch and the new 3.x branch.
There are language incompatibilities between these branches and one has to be aware of
which one to use. This book is based on the 3.x branch, considering the language up to
release 3.5.
For this book you need to install the following:
The interpreter: Python 3.5 (or later)
The modules for scientific computing: SciPy with NumPy
The module for graphical representation of mathematical results: matplotlib
The shell: IPython
A Python related editor: Spyder (refer to the following Figure 1.1, Spyder), Geany
The installation of these is eased by the so-called distribution packages. We recommend that
you use Anaconda. The default screen of Spyder consists of an editor window on left, a
console window in the lower right corner which gives access to an IPython shell and a help
window in the upper right corner as shown in the following figure:
[ 11 ]
Getting Started
Figure 1.1: The default screen of Spyder consists of an editor window on left, a
console window in the lower right corner which gives access to an IPython shell
and a help window in the upper right corner.
Anaconda
Even if you have Python pre-installed on your computer, we recommend that you create
your personal Python environment that allows you to work without the risk of accidentally
affecting the software on which your computer's functionality might depend. With a virtual
environment, such as Anaconda, you are free to change language versions and install
packages without the unintended side-effects.
If the worst happens and you screw things up totally, just delete the Anaconda directory
and start again. Running the Anaconda installer will install Python, a Python development
environment and editor (Spyder), the shell IPython, and the most important packages for
numerical computations, for example SciPy, NumPy, and matplotlib.
[ 12 ]
Getting Started
You can install additional packages with conda install within your virtual environment
created by Anaconda (refer for official documentation from[2]) .
Configuration
Most Python codes will be collected in files. We recommend that you use the following
header in all your Python files:
from scipy import *
from matplotlib.pyplot import *
With this, you make sure that all standard modules and functions used in this book, such as
SciPy, are imported. Without this step, most of the examples in the book would raise errors.
Many editors, such as Spyder, provide the possibility to create a template for your files.
Look for this feature and put the preceding header into a template.
Python Shell
The Python shell is good but not optimal for interactive scripting. We therefore recommend
using IPython instead (refer to [26] for the official documentation). IPython can be started
in different ways:
In a terminal shell by running the following command: ipython
By directly clicking on an icon called Jupyter QT Console
When working with Spyder you should use an IPython console (refer to Figure
1.1, Spyder).
[ 13 ]
Getting Started
Executing scripts
You often want to execute the contents of a file. Depending on the location of the file on
your computer, it is necessary to navigate to the correct location before executing the
contents of a file.
Use the command cd in IPython in order to move to the directory where your file
is located.
To execute the contents of a file named myfile.py, just run the following
command in the IPython shell
run myfile
Getting Help
Here are some tips on how to use IPython:
To get help on an object, just type ? after the object's name and then return.
Use the arrow keys to reuse the last executed commands.
You may use the Tab key for completion (that is, you write the first letter of a
variable or method and IPython shows you a menu with all the possible
completions).
Use Ctrl+D to quit.
Use IPython's magic functions. You can find a list and explanations by applying
%magic on the command prompt.
You can find out more about IPython in its online documentation, [15].
Jupyter – Python notebook
The Jupyter notebook is a fantastic tool for demonstrating your work. Students might want
to use it to make and document homework and exercises and teachers can prepare lectures
with it, even slides and web pages.
[ 14 ]
Getting Started
If you have installed Python via Anaconda, you already have everything for Jupyter in
place. You can invoke the notebook by running the following command in the terminal
window:
jupyter notebook
A browser window will open and you can interact with Python through your web browser.
Program and program flow
A program is a sequence of statements that are executed in a top-down order. This linear
execution order has some important exceptions:
There might be a conditional execution of alternative groups of statements
(blocks), which we refer to as branching.
There are blocks that are executed repetitively, which is
called looping (refer to the following Figure 1.2, Program flow).
There are function calls that are references to another piece of code, which is
executed before the main program flow is resumed. A function call breaks the
linear execution and pauses the execution of a program unit while it passes the
control to another unit–a function. When this gets completed, its control is
returned to the calling unit.
Figure 1.2: Program flow
[ 15 ]
Getting Started
Python uses a special syntax to mark blocks of statements: a keyword, a colon, and an
indented sequence of statements, which belong to the block (refer to the following Figure
1.3, Block command).
Figure 1.3: Block command
Comments
If a line in a program contains the symbol #, everything following on the same line is
considered as a comment:
# This is a comment of the following statement
a = 3 # ... which might get a further comment here
Line joining
A backslash \ at the end of the line marks the next line as a continuation line, that is, explicit
line joining. If the line ends before all the parentheses are closed, the following line will
automatically be recognized as a continuation line, that is, implicit line joining.
[ 16 ]
Getting Started
Basic types
Let's go over the basic data types that you will encounter in Python.
Numbers
A number may be an integer, a real number, or a complex number. The usual operations
are:
addition and subtraction, + and multiplication and division, * and /
power, **
Here is an example:
2 ** (2 + 2) # 16
1j ** 2 # -1
1. + 3.0j
The symbol for complex numbers
j is a symbol to denote the imaginary part of a complex number.
It is a syntactic element and should not be confused with multiplication by
a variable. More on complex numbers can be found in section Numeric
Types of Chapter 2, Variables and Basic Types.
Strings
Strings are sequences of characters, enclosed by simple or double quotes:
'valid string'
"string with double quotes"
"you shouldn't forget comments"
'these are double quotes: ".." '
You can also use triple quotes for strings that have multiple lines:
"""This is
a long,
long string"""
[ 17 ]
Getting Started
Variables
A variable is a reference to an object. An object may have several references. One uses the
assignment operator = to assign a value to a variable:
x =
y =
del
del
[3,
x #
x #
y #
4] # a list object is created
this object now has two labels: x and y
we delete one of the labels
both labels are removed: the object is deleted
The value of a variable can be displayed by the print function:
x = [3, 4] # a list object is created
print(x)
Lists
Lists are a very useful construction and one of the basic types in Python. A Python list is an
ordered list of objects enclosed by square brackets. One can access the elements of a list
using zero-based indexes inside square brackets:
L1 = [5, 6]
L1[0] # 5
L1[1] # 6
L1[2] # raises IndexError
L2 = ['a', 1, [3, 4]]
L2[0] # 'a'
L2[2][0] # 3
L2[-1] # last element: [3,4]
L2[-2] # second to last: 1
Indexing of the elements starts at zero. One can put objects of any type inside a list, even
other lists. Some basic list functions are as follows:
list(range(n))} creates a list with n elements, starting with zero:
print(list(range(5))) # returns [0, 1, 2, 3, 4]
len gives the length of a list:
len(['a', 1, 2, 34]) # returns 4
[ 18 ]
Getting Started
append is used to append an element to a list:
L = ['a', 'b', 'c']
L[-1] # 'c'
L.append('d')
L # L is now ['a', 'b', 'c', 'd']
L[-1] # 'd'
Operations on lists
The operator + concatenates two lists:
L1 = [1, 2]
L2 = [3, 4]
L = L1 + L2 # [1, 2, 3, 4]
As one might expect, multiplying a list with an integer concatenates the list with
itself several times:
n*L is equivalent to making n additions.
L = [1, 2]
3 * L # [1, 2, 1, 2, 1, 2]
Boolean expressions
A Boolean expression is an expression that may have the value True or False. Some
common operators that yield conditional expressions are as follow:
Equal, ==
Not equal, !=
Less than, Less than or equal to, < , <=
Greater than, Greater than or equal to, > , >=
[ 19 ]
Getting Started
One combines different Boolean values with or and and.
The keyword not , gives the logical negation of the expression that follows. Comparisons
can be chained so that, for example, x < y < z is equivalent to x < y and y < z. The
difference is that y is only evaluated once in the first example.
In both cases, z is not evaluated at all when the first condition, x < y, evaluates to False:
2 >= 4 # False
2 < 3 < 4 # True
2 < 3 and 3 < 2 # False
2 != 3 < 4 or False # True
2 <= 2 and 2 >= 2 # True
not 2 == 3 # True
not False or True and False # True!
Precedence rules
The <, >, <=, >=, !=, and == operators have higher precedence than not.
The operators and, or have the lowest precedence. Operators with higher
precedence rules are evaluated before those with lower.
Repeating statements with loops
Loops are used to repetitively execute a sequence of statements while changing a variable
from iteration to iteration. This variable is called the index variable. It is successively
assigned to the elements of a list, (refer to Chapter 9, Iterating):
L = [1, 2, 10]
for s in L:
print(s * 2) # output: 2 4 20
The part to be repeated in the for loop has to be properly indented:
for elt in my_list:
do_something
something_else
print("loop finished") # outside the for block
Repeating a task
One typical use of a for loop is to repeat a certain task a fixed number of times:
n = 30
for iteration in range(n):
do_something # this gets executed n times
[ 20 ]
Getting Started
Break and else
The for statement has two important keywords: break and else. break quits the for loop
even if the list we are iterating is not exhausted:
for x in x_values:
if x > threshold:
break
print(x)
The finalizing else checks whether the for loop was broken with the break keyword. If it
was not broken, the block following the else keyword is executed:
for x in x_values:
if x > threshold:
break
else:
print("all the x are below the threshold")
Conditional statements
This section covers how to use conditions for branching, breaking, or otherwise controlling
your code. A conditional statement delimits a block that will be executed if the condition is
true. An optional block, started with the keyword else will be executed if the condition is
not fulfilled (refer to Figure 1.3, Block command diagram). We demonstrate this by printing
|x|, the absolute value of x:
The Python equivalent is as follows:
x = ...
if x >= 0:
print(x)
else:
print(-x)
Any object can be tested for the truth value, for use in an if or while statement. The rules
for how the truth values are obtained are explained in section Boolean of Chapter 2,
Variables and Basic Types.
[ 21 ]
Getting Started
Encapsulating code with functions
Functions are useful for gathering similar pieces of code in one place. Consider the
following mathematical function:
The Python equivalent is as follows:
def f(x):
return 2*x + 1
In Figure 1.4 Anatomy of a function the elements of a function block are explained.
The keyword def tells Python we are defining a function.
f is the name of the function.
x is the argument, or input of the function.
What is after return is called the output of the function.
Figure 1.4: Anatomy of a function
Once the function is defined, it can be called using the following code:
f(2) # 5
f(1) # 3
[ 22 ]
Getting Started
Scripts and modules
A collection of statements in a file (which usually has a py extension), is called a script.
Suppose we put the contents of the following code into a file named smartscript.py:
def f(x):
return 2*x + 1
z = []
for x in range(10):
if f(x) > pi:
z.append(x)
else:
z.append(-1)
print(z)
In a Python or IPython shell, such a script can then be executed with the exec command
after opening and reading the file. Written as a one-liner it reads:
exec(open('smartscript.py').read())
The IPython shell provides the magic command %run as a handy alternative way to execute
a script:
%run smartscript
Simple modules – collecting functions
Often one collects functions in a script. This creates a module with additional Python
functionality. To demonstrate this, we create a module by collecting functions in a single
file, for example smartfunctions.py:
def f(x):
return 2*x + 1
def g(x):
return x**2 + 4*x - 5
def h(x):
return 1/f(x)
These functions can now be used by any external script or directly in the IPython
environment.
Functions within the module can depend on each other.
Grouping functions with a common theme or purpose gives modules that can be
shared and used by others.
[ 23 ]
Getting Started
Again, the command exec(open('smartfunctions.py').read()) makes these
functions available to your IPython shell (note that there is also the IPython magic function
run). In Python terminology, one says that they are put into the actual namespace.
Using modules and namespaces
Alternatively, the modules can be imported by the command import. It creates a
named namespace. The command from puts the functions into the general namespace:
import smartfunctions
print(smartfunctions.f(2))
# 5
from smartfunctions import g
print(g(1)) # 0
#import just this function
from smartfunctions import *
print(h(2)*f(2))
#import all
# 1.0
Import
The commands import and from import the functions only once into the
respective namespace. Changing the functions after the import has no
effect for the current Python session. More on modules can be found in
section Modules of Chapter 11, Namespaces, Scopes and Modules.
Interpreter
The Python interpreter executes the following steps:
First, run the syntax.
Then execute the code line by line.
Code inside a function or class declaration is not executed (but checked for
syntax).
def f(x):
return y**2
a = 3
# here both a and f are defined
[ 24 ]
Getting Started
You can run the preceding program because there are no syntactical errors. You get an error
only when you call the function f.
f(2) # error, y is not defined
Summary
In this chapter, we briefly addressed the main language elements of Python without going
into detail.
You should now be able to start playing with small pieces of code and to test different
program constructs. All this is intended as an appetizer for the following chapters in which
we will give you the details, examples, exercises, and more background information.
[ 25 ]
2
Variables and Basic Types
In this chapter, we will present the most important and basic types in Python. What is a
type? It is a set consisting of data content, its representation, and all possible operations.
Later in this book, we will make this definition much more precise, when we introduce the
concepts of a class in Chapter 8, Classes.
Variables
Variables are references to Python objects. They are created by assignments, for example:
a = 1
diameter = 3.
height = 5.
cylinder = [diameter, height] # reference to a list
Variables take names that consist of any combination of capital and small letters, the
underscore _ , and digits. A variable name must not start with a digit. Note that variable
names are case sensitive. A good naming of variables is an essential part of documenting
your work, so we recommend that you use descriptive variable names.
Variables and Basic Types
Python has some reserved keywords, which cannot be used as variable names (refer to
following table, Table 2.1). An attempt to use such a keyword as variable name would raise
a syntax error.
Table 2.1: Reserved Python keywords.
As opposed to other programming languages, variables require no type declaration. You
can create several variables with a multiple assignment statement:
a = b = c = 1
# a, b and c get the same value 1
Variables can also be altered after their definition:
a = 1
a = a + 1 # a gets the value 2
a = 3 * a
# a gets the value 6
The last two statements can be written by combining the two operations with an assignment
directly by using increment operators:
a += 1
a *= 3
# same as a = a + 1
# same as a = 3 * a
Numeric types
At some point, you will have to work with numbers, so we start by considering different
forms of numeric types in Python. In mathematics, we distinguish between natural
numbers (ℕ), integers (ℤ), rational numbers (ℚ), real numbers (ℝ) and complex numbers
(ℂ). These are infinite sets of numbers. Operations differ between these sets and may even
not be defined. For example, the usual division of two numbers in ℤ might not result in an
integer — it is not defined on ℤ.
[ 27 ]
Variables and Basic Types
In Python, like many other computer languages, we have numeric types:
The numeric type int, which is at least theoretically the entire ℤ
The numeric type float, which is a finite subset of ℝ and
The numeric type complex, which is a finite subset of ℂ
Finite sets have a smallest and a largest number and there is a minimum spacing between
two numbers; refer to the section on Floating Point Representation for further details.
Integers
The simplest numerical type is the integer type.
Plain integers
The statement k = 3 assigns the variable k to an integer.
Applying an operation of the type +, -, or * to integers returns an integer. The division
operator, //, returns an integer, while / may return a float:
6 // 2
7 // 2
7 / 2
# 3
# 3
# 3.5
The set of integers in Python is unbounded; there is no largest integer. The limitation here is
the computer’s memory rather than any fixed value given by the language.
If the division operator (/) in the example returns 3, you might not have
installed the correct Python version.
[ 28 ]
Variables and Basic Types
Floating point numbers
If you execute the statement a = 3.0 in Python, you create a floating-point number
(Python type: float). These numbers form a subset of rational numbers, ℚ.
Alternatively the constant could have been given in exponent notation as a = 30.0e-1 or
simply a = 30.e-1. The symbol e separates the exponent from the mantissa, and the
expression reads in mathematical notation a = 30.0 × 10−1. The name floating-point number
refers to the internal representation of these numbers and reflects the floating position of the
decimal point when considering numbers over a wide range.
Applying the elementary mathematical operations +, -, *, and / to two floating-point
numbers or to an integer and a floating-point number returns a floating-point number.
Operations between floating-point numbers rarely return the exact result expected from
rational number operations:
0.4 - 0.3 # returns 0.10000000000000003
This facts matters, when comparing floating point numbers:
0.4 - 0.3 == 0.1 # returns False
Floating point representation
Internally, floating-point numbers are represented by four quantities: the sign, the mantissa,
the exponent sign, and the exponent:
with β ϵ ℕ and x0≠ 0, 0 ≤ xi≤ β
x0…xt-1 is called the mantissa, β the basis and e the exponent |e| ≤ U . t is called the mantissa
length. The condition x0 ≠ 0 makes the representation unique and saves, in the binary case (β
= 2), one bit.
There exist two-floating point zeros +0 and -0, both represented by the mantissa 0.
On a typical Intel processor, β = 2 . To represent a number in the float type 64 bits are
used, namely 2 bits for the signs, t = 52 bits for the mantissa and 10 bits for the exponent
10
|e|. The upper bound U for the exponent is consequently 2 -1 = 1023.
[ 29 ]
Variables and Basic Types
With this data the smallest positive representable number is
flmin = 1.0 × 2-1023≈ 10-308 and the largest is flmax = 1.111…1 × 21023≈ 10308.
Note that floating-point numbers are not equally spaced in [0, flmax]. There is in particular a
gap at zero (refer to [29]). The distance between 0 and the first positive number is 2-1023,
while the distance between the first and the second is smaller by a factor 2-52≈ 2.2 × 10-16. This
effect, caused by the normalization x0 ≠ 0, is visualized in Figure 2.1.
This gap is filled equidistantly with subnormal floating-point numbers to which such a
result is rounded. Subnormal floating-point numbers have the smallest possible exponent
and do not follow the convention that the leading digit x0 has to differ from zero; refer to
[13].
Infinite and not a number
There are in total
floating-point numbers. Sometimes a numerical
algorithm computes floating-point numbers outside this range.
This generates number over- or underflow. In SciPy the special floating-point number inf
is assigned to overflow results:
exp(1000.) # inf
a = inf
3 - a
# -inf
3 + a
# inf
Working with inf may lead to mathematically undefined results. This is indicated in
Python by assigning the result another special floating-point number, nan. This stands for
not-a-number, that is, an undefined result of a mathematical operation:
a + a # inf
a - a # nan
a / a # nan
[ 30 ]
Variables and Basic Types
There are special rules for operations with nan and inf. For instance, nan compared to
anything (even to itself) always returns False:
x
x
x
x
= nan
< 0 # False
> 0 # False
== x # False
See Exercise 4 for some surprising consequences of the fact that nan is never equal to itself.
The float inf behaves much more as expected:
0 < inf
#
inf <= inf #
inf == inf #
-inf < inf #
inf - inf
#
exp(-inf)
#
exp(1 / inf)
True
True
True
True
nan
0
# 1
One way to check for nan and inf is to use the isnan and isinf functions. Often, one
wants to react directly when a variable gets the value nan or inf. This can be achieved by
using the NumPy command seterr. The following command
seterr(all = 'raise')
would raise an error if a calculation were to return one of those values.
Underflow – Machine Epsilon
Underflow occurs when an operation results in a rational number that falls into the gap at
zero; refer to Figure 2.1.
Figure 2.1: The floating point gap at zero, here t = 3, U = 1
The machine epsilon or rounding unit is the largest number ε such that float(1.0 + ε) = 1.0.
[ 31 ]
Variables and Basic Types
Note that ε ≈ β1-t/2 = 1.1102 × 10-16 on most of today’s computers. The value that is valid on
the actual machine you are running your code on is accessible using the following
command:
import sys
sys.float_info.epsilon # 2.220446049250313e-16 (something like that)
The variable sys.float_info contains more information about the internal
representation of the float type on your machine.
The function float converts other types to a floating-point number—if possible. This
function is especially useful when converting an appropriate string to a number:
a = float('1.356')
Other float types in NumPy
NumPy also provides other float types, known from other programming languages as
double-precision and single-precision numbers, namely float64 and float32:
a = pi
a1 = float64(a)
a2 = float32(a)
a - a1
a - a2
#
#
#
#
#
returns
returns
returns
returns
returns
3.141592653589793
3.1415926535897931
3.1415927
0.0
-8.7422780126189537e-08
The second last line demonstrates that a and a1 do not differ in accuracy. In the first two
lines, they only differ in the way they are displayed. The real difference in accuracy is
between a and its single-precision counterpart, a2.
The NumPy function finfo can be used to display information on these floating-point
types:
f32 = finfo(float32)
f32.precision
# 6 (decimal digits)
f64 = finfo(float64)
f64.precision
# 15 (decimal digits)
f = finfo(float)
f.precision
# 15 (decimal digits)
f64.max
# 1.7976931348623157e+308 (largest number)
f32.max
# 3.4028235e+38 (largest number)
help(finfo)
# Check for more options
[ 32 ]
Variables and Basic Types
Complex numbers
Complex numbers are an extension of the real numbers frequently used in many scientific
and engineering fields.
Complex Numbers in Mathematics
Complex numbers consist of two floating-point numbers, the real part a of the number and
its imaginary part b. In mathematics, a complex number is written as z=a+bi, where i defined
by i2 = -1 is the imaginary unit. The conjugate complex counterpart of z is
.
If the real part a is zero, the number is called an imaginary number.
The j notation
In Python, imaginary numbers are characterized by suffixing a floating-point number with
the letter j, for example, z = 5.2j. A complex number is formed by the sum of a floatingpoint number and an imaginary number, for example, z = 3.5 + 5.2j.
While in mathematics the imaginary part is expressed as a product of a real number b with
the imaginary unit i, the Python way of expressing an imaginary number is not a product: j
is just a suffix to indicate that the number is imaginary.
This is demonstrated by the following small experiment:
b
z
z
z
=
=
=
=
5.2
bj
# returns a NameError
b*j # returns a NameError
b*1j # is correct
The method conjugate returns the conjugate of z:
z = 3.2 + 5.2j
z.conjugate() # returns (3.2-5.2j)
[ 33 ]
Variables and Basic Types
Real and imaginary parts
One may access the real and imaginary parts of a complex number z using the real and
imag attributes. Those attributes are read-only:
z = 1j
z.real
z.imag
z.imag = 2
# 0.0
# 1.0
# AttributeError: readonly attribute
It is not possible to convert a complex number to a real number:
z = 1 + 0j
z == 1
# True
float(z)
# TypeError
Interestingly, the real and imag attributes as well as the conjugate method work just as
well for complex arrays (Chapter 4, Linear Algebra – Arrays). We demonstrate this by
computing the Nth roots of unity which are
, that is, the N solutions
of the equation
:
N = 10
# the following vector contains the Nth roots of unity:
unity_roots = array([exp(1j*2*pi*k/N) for k in range(N)])
# access all the real or imaginary parts with real or imag:
axes(aspect='equal')
plot(unity_roots.real, unity_roots.imag, 'o')
allclose(unity_roots**N, 1) # True
[ 34 ]
Variables and Basic Types
The resulting figure (Figure 2.2) shows the roots of unity together with the unit circle. (For
more details on how to make plots, refer Chapter 6, Plotting.)
Figure 2.2: Roots of unity together with a unit circle
It is of course possible to mix the previous methods, as illustrated by the following
examples:
z = 3.2+5.2j
(z + z.conjugate()) / 2.
# returns (3.2+0j)
((z + z.conjugate()) / 2.).real
# returns 3.2
(z - z.conjugate()) / 2.
# returns 5.2j
((z - z.conjugate()) / 2.).imag
# returns 5.2
sqrt(z * z.conjugate())
# returns (6.1057350089894991+0j)
[ 35 ]
Variables and Basic Types
Booleans
Boolean is a datatype named after George Boole (1815-1864). A Boolean variable can take
only two values, True or False. The main use of this type is in logical expressions. Here are
some examples:
a = True
b = 30 > 45
# b gets the value False
Boolean expressions are often used in conjunction with the if statement:
if x > 0:
print("positive")
else:
print("nonpositive)
Boolean operators
Boolean operations are performed using the and, or, and not keywords in Python:
True and False # False
False or True # True
(30 > 45) or (27 < 30) # True
not True # False
not (3 > 4) # True
The operators follow some precedence rules (refer to section Executing scripts in Chapter 1,
Getting started) which would make the parentheses in the third line and in the last obsolete
(it is a good practice to use them anyway to increase the readability of your code). Note that
the and operator is implicitly chained in the following Boolean expressions:
a < b < c
a == b == c
# same as: a < b and b < c
# same as: a == b and b == c
[ 36 ]
Variables and Basic Types
Rules of Conversion to Booleans:
Table 2.2 : Rule of conversion to Boolean
Boolean casting
Most objects Python may be converted to Booleans; this is called Boolean casting. The built-in
function bool performs that conversion. Note that most objects are cast to True, except 0,
the empty tuple, the empty list, the empty string, or the empty array. These are all cast to
False.
It is not possible to cast arrays into Booleans unless they contain no or only one element; this
is explained further in Chapter 5, Advanced Array Concepts. The previous table contains
summarized rules for Boolean casting. Some usage examples:
bool([])
# False
bool(0)
# False
bool(' ')
# True
bool('')
# False
bool('hello')
# True
bool(1.2)
# True
bool(array([1]))
# True
bool(array([1,2]))
# Exception raised!
[ 37 ]
Variables and Basic Types
Automatic Boolean casting
Using an if statement with a non-Boolean type will cast it to a Boolean. In other words, the
following two statements are always equivalent:
if a:
...
if bool(a): # exactly the same as above
...
A typical example is testing whether a list is empty:
# L is a list
if L:
print("list not empty")
else:
print("list is empty")
An empty array, list, or tuple will return False. You can also use a variable in the if
statement, for example, an integer:
# n is an integer
if n % 2:
print("n is odd")
else:
print("n is even")
Note that we used % for the modulo operation, which returns the remainder of an integer
division. In this case, it returns 0 or 1 as the remainder after modulo 2.
In this last example, the values 0 or 1 are cast to bool. Boolean operators or,and , and not
will also implicitly convert some of their arguments to a Boolean.
Return values of and and or
Note that the operatorsand and or do not necessarily produce Boolean values. The
expression x and y is equivalent to:
def and_as_function(x,y):
if not x:
return x
else:
return y
[ 38 ]
Variables and Basic Types
And the expression x or y is equivalent to:
def or_as_function(x,y):
if x:
return x
else:
return y
Interestingly, this means that when executing the statement True or x, the variable x need
not even be defined! The same holds for False and x.
Note that, unlike their counterparts in mathematical logic, these operators are no longer
commutative in Python. Indeed, the following expressions are not equivalent:
[1] or 'a' # produces [1]
'a' or [1] # produces 'a'
Boolean and integer
In fact, Booleans and integers are the same. The only difference is in the string
representation of 0 and 1 which is in the case of Booleans False and True respectively. This
allows constructions like this (for the format method refer section on string formatting):
def print_ispositive(x):
possibilities = ['nonpositive', 'positive']
return "x is {}".format(possibilities[x>0])
We note for readers already familiar with the concept of subclasses, that the type bool is a
subclass of the type int (refer to Chapter 8, Classes). Indeed, all four inquiries
isinstance(True, bool), isinstance(False, bool), isinstance(True, int),
and isinstance(False, int) return the value True (refer to section Type Checking in
Chapter 3, Container Types).
Even rarely used statements such as True+13 are syntactically correct.
Strings
The type string is a type used for text:
name = 'Johan Carlsson'
child = "Åsa is Johan Carlsson's daughter"
book = """Aunt Julia
and the Scriptwriter"""
[ 39 ]
Variables and Basic Types
A string is enclosed either by single or double quotes. If a string contains several lines, it has
to be enclosed by three double quotes """ or three single quotes '''.
Strings can be indexed with simple indexes or slices (refer to Chapter 3, Container Types, for
a comprehensive explanation on slices):
book[-1] # returns 'r'
book[-12:] # returns 'Scriptwriter'
Strings are immutable; that is, items cannot be altered. They share this property with tuples.
The command book[1] = 'a' returns:
TypeError: 'str' object does not support item assignment
The string '\n' is used to insert a line break and 't' inserts a horizontal tabulator (TAB)
into the string to align several lines:
print('Temperature:\t20\tC\nPressure:\t5\tPa')
These strings are examples of escape sequences. Escape sequences always start with a
backslash, \ . A multi line string automatically includes escape sequences:
a="""
A multiline
example"""
a # returns '\nA multiline \nexample'
A special escape sequence is "", which represents the backslash symbol in text:
latexfontsize="\\tiny"
The same can be achieved by using a raw string instead:
latexfs=r"\tiny"
# returns "\tiny"
latexfontsize == latexfs # returns True
Note that in raw strings, the backslash remains in the string and is used to escape some
special characters:
r"\"\"
r"\\"
r"\"
# returns '\\"'
# returns '\\\\'
# returns an error
[ 40 ]
Variables and Basic Types
Operations on strings and string methods
Addition of strings means concatenation:
last_name = 'Carlsson'
first_name = 'Johanna'
full_name = first_name + ' ' + last_name
# returns 'Johanna Carlsson'
Multiplication is just repeated addition:
game = 2 * 'Yo' # returns 'YoYo'
When strings are compared, lexicographical order applies and the uppercase form precedes
the lowercase form of the same letter:
'Anna' > 'Arvi' # returns false
'ANNA' < 'anna' # returns true
'10B' < '11A'
# returns true
Among the variety of string methods, we will mention here only the most important ones:
Splitting a string: This method generates a list from a string by using a single or
multiple blanks as separators. Alternatively, an argument can be given by
specifying a particular string as a separator:
text = 'quod erat
demonstrandum'
text.split() # returns ['quod', 'erat', 'demonstrandum']
table = 'Johan;Carlsson;19890327'
table.split(';') # returns ['Johan','Carlsson','19890327']
king = 'CarlXVIGustaf'
king.split('XVI') # returns ['Carl','Gustaf']
Joining a list to a string: This is the reverse operation of splitting:
sep = ';'
sep.join(['Johan','Carlsson','19890327'])
# returns 'Johan;Carlsson;19890327'
Searching in a string: This method returns the first index in the string, where a
given search substring starts:
birthday = '20101210'
birthday.find('10') # returns 2
If the search string is not found, the return value of the method is -1 .
[ 41 ]
Variables and Basic Types
String formatting
String formatting is done using the format method:
course_code = "NUMA21"
print("This course's name is {}".format(course_code))
# This course's name is NUMA21
The function format is a string method; it scans the string for the occurrence of
placeholders, which are enclosed by curly brackets. These placeholders are replaced in a
way specified by the argument of the format method. How they are replaced depends on
the format specification defined in each {} pair. Format specifications are indicated by a
colon, ":", as their prefix.
The format method offers a range of possibilities to customize the formatting of objects
depending on their types. Of particular use in scientific computing are the formatting
specifiers for the float type. One may choose either the standard with {:f} or the
exponential notation with {:e}:
quantity = 33.45
print("{:f}".format(quantity)) # 33.450000
print("{:1.1f}".format(quantity)) # 33.5
print("{:.2e}".format(quantity)) # 3.35e+01
The format specifiers allow to specify the rounding precision (digits following the decimal
point in the representation). Also the total number of symbols including leading blanks to
represent the number can be set.
In this example, the name of the object that gets its value inserted is given as an argument to
the format method. The first {} pair is replaced by the first argument and the following
pairs by the subsequent arguments. Alternatively, it may also be convenient to use the keyvalue syntax:
print("{name} {value:.1f}".format(name="quantity",value=quantity))
# prints "quantity 33.5"
Here, two values are processed, a string name without a format specifier and a float
value that is printed in fixed point notation with one digit after the decimal point. (We
refer to the complete reference documentation for more details on string formatting, [34]).
[ 42 ]
Variables and Basic Types
Braces in the string
Sometimes, a string might contain a pair of curly braces, which should not
be considered as placeholders for a format method. In that case, double
braces are used:
r"we {} in LaTeX \begin{{equation}}".format('like')
This returns the following string: 'we like in LaTeX
\\begin{equation}'.
Summary
In this chapter, you met the basic data types in Python and saw the corresponding syntax
elements. We will work mostly with numeric types such as integers, floats and complex.
Booleans are needed for setting conditions, and by using strings, we often communicate
results and messages.
Exercises
Ex. 1 → Check whether x = 2.3 is a zero of the function:
Ex. 2 → According to de Moivre's formula, the following holds:
Choose numbers n and x and verify that formula in Python.
Ex. 3 → Complex numbers. Verify Euler's formula in the same way:
[ 43 ]
Variables and Basic Types
Ex. 4 → Suppose we are trying to check the convergence of a diverging sequence (here the
sequence is defined by the recursive relation un+1 = 2un and u0 = 1.0):
u = 1.0 # you have to use a float here!
uold = 10.
for iteration in range(2000):
if not abs(u-uold) > 1.e-8:
print('Convergence')
break # sequence has converged
uold = u
u = 2*u
else:
print('No convergence')
1. Since the sequence does not converge, the code should print the No
convergence message. Execute it to see what happens.
2. What happens if you replace the line:
if not abs(u-uold) > 1.e-8
with:
if abs(u-uold) < 1.e-8
It should give exactly the same result, shouldn't it? Run the code again to see
what happens.
3. What happens if you replace u=1.0 by u=1 (without decimal point). Run the code
to check your predictions.
4. Explain the unexpected behavior of this code. The key to understand what
happens is that inf evaluates to nan, and the comparison of nan with anything
else is returns always the value False .
Ex. 5 → An implication C = (A ⇒ B) is a Boolean expression that is defined as
C is True if A is False or A and B are both True
C is False otherwise
Write a Python function implication(A, B).
[ 44 ]
Variables and Basic Types
Ex. 6 → This exercise is to train Boolean operations. Two binary digits (bits) are added by
using a logical device called a half adder. It produces a carry bit (the digit of the next higher
value) and the sum as defined by the following table, and half adder circuit.
p
q
sum
carry
1
1
0
1
1
0
1
0
0
1
1
0
0
0
0
0
Definition of the half adder operation
Figure 2.3: A half adder circuit
[ 45 ]
Variables and Basic Types
A full adder consists of two half adders and sums up two bits and an additional carry bit on
the input (refer to the following figure):
Figure 2.4: A full adder circuit
Write a function that implements a half adder and another that implements a full adder.
Test these functions.
[ 46 ]
3
Container Types
Container types are used to group objects together. The main difference between the
different container types is the way individual elements are accessed and how operations
are defined.
Lists
A list is, as the name hints, a list of objects of any kind:
L = ['a' 20.0, 5]
M = [3,['a', -3.0, 5]]
The individual objects are enumerated by assigning each element an index. The first
element in the list gets index 0. This zero-based indexing is frequently used in mathematical
notation. Consider the usual indexing of coefficients of a polynomial.
The index allows us to access the following objects:
L[1] # returns 20.0
L[0] # returns 'a'
M[1] # returns ['a',-3.0,5]
M[1][2] # returns 5
The bracket notation here corresponds to the use of subscripts in mathematical formulas. L
is a simple list, while M itself contains a list so that one needs two indexes to access an
element of the inner list.
A list containing subsequent integers can easily be generated by the command range:
L=list(range(4)) # generates a list with four elements: [0, 1, 2 ,3]
Container Types
A more general use is to provide this command with start, stop, and step parameters:
L=list(range(17,29,4)) # generates [17, 21, 25]
The command len returns the length of the list:
len(L) # returns 3
Slicing
Slicing a list between i and j creates a new list containing the elements starting at index
i and ending just before j.
For slicing, a range of indexes has to be given. L[i:j] means create a list by taking all
elements from L starting at L[i] until L[j-1]. In other words, the new list is obtained by
removing the first i elements from L and taking the next j-i elements (for j > i ≥ 0). See
the following figure (Figure 3.1) for more examples:
Figure 3.1: Some typical slicing situations
[ 48 ]
Container Types
Here, L[i:] means remove the i first elements, L[:i] means take only the first i elements,
and similarly, L[:-i] means remove the last i elements, and L[-i:] means take only the
last i elements. This may be combined in L[i:-j] to remove the first i and the last j
elements:
L = ['C', 'l', 'o', 'u', 'd', 's']
L[1:5] # remove one element and take four from there:
# returns ['l', 'o', 'u', 'd']
One may omit the first or last bound of the slicing:
L = ['C', 'l', 'o', 'u','d', 's']
L[1:] # ['l', 'o', 'u', 'd','s']
L[:5] # ['C', 'l', 'o','u','d']
L[:] # the entire list
Python allows the use of negative indexes for counting from the right. In particular, the
element L[-1] is the last element in the list L.
Some list indexing descriptions:
L[i:] amounts to taking all elements except the i first ones
L[:i] amounts to taking the first i elements
L[-i:] amounts to taking the last i elements
L[:-i] amounts to taking all elements except the i last ones
Here is an example:
L = ['C', 'l', 'o', 'u', 'd', 's']
L[-2:] # ['d', 's']
L[:-2] # ['C', 'l', 'o','u']
Omitting one index in the range corresponds to half-open intervals in ℝ. The half-open
interval (∞, a) means, take all numbers strictly lower than a; this is similar to the syntax
L[:j].
Out of-bound slices
Notice that you never get index errors with out-of-bound slices. Possibly,
you may obtain empty lists.
[ 49 ]
Container Types
Here is an example:
L = list(range(4)) # [0, 1, 2, 3]
L[4] # IndexError: list index out of range
L[1:100] # same as L[1:]
L[-100:-1] # same as L[:-1]
L[-100:100] # same as L[:]
L[5:0] # empty list []
L[-2:2] # empty list []
Be careful when using variables in indexing that may become negative, since it changes the
slice completely. This might lead to unexpected results:
a = [1,2,3]
for iteration in range(4):
print(sum(a[0:iteration-1]))
The result is 3, 0, 1, 3 while one expects 0, 0, 1, 3.
Strides
When computing slices, one may also specify a stride, which is the length of the step from
one index to the other. The default stride is one. Here is an example:
L = list(range(100))
L[:10:2] # [0, 2, 4, 6, 8]
L[::20] # [0, 20, 40, 60, 80]
L[10:20:3] # [10, 13, 16, 19]
Note that the stride may also be negative:
L[20:10:-3] # [20, 17, 14, 11]
It is also possible to create a new list that is reversed, using a negative stride (find about
reverse method in section In-place operations):
L = [1, 2, 3]
R = L[::-1] # L is not modified
R # [3, 2, 1]
[ 50 ]
Container Types
Altering lists
Typical operations on lists are insertion and deletion of elements and list concatenation.
With the slicing notation, list insertion and deletion become obvious; deletion is just
replacing a part of a list by an empty list []:
L = ['a', 1, 2, 3, 4]
L[2:3] = [] # ['a', 1, 3, 4]
L[3:] = [] # ['a', 1, 3]
Insertion means replacing an empty slice with the list to be inserted:
L[1:1] = [1000, 2000] # ['a', 1000, 2000, 1, 3]
Two lists are concatenated by the plus operator + :
L = [1, -17]
M = [-23.5, 18.3, 5.0]
L + M # gives [1, -17, 23.5, 18.3, 5.0]
Concatenating a list n times with itself motivates the use of the multiplication operator *:
n = 3
n * [1.,17,3] # gives [1., 17, 3, 1., 17, 3, 1., 17, 3]
[0] * 5 # gives [0,0,0,0,0]
There is no arithmetic operations on list, such as elementwise summation or division. For
such operations we use arrays (refer to section Array).
Belonging to a list
One may use the keywords in and not in to determine whether an element belongs to a
list or not which is similar to and in mathematics:
L = ['a', 1, 'b', 2]
'a' in L # True
3 in L # False
4 not in L # True
[ 51 ]
Container Types
List methods
Some useful methods of the list type are collected in the following Table 3.1:
Command
Action
list.append(x)
Add x to the end of the list.
list.expand(L)
Expand the list by the elements of the list L.
list.insert(i,x) Insert x at positioni.
list.remove(x)
Remove the first item from the list whose value is x.
list.count(x)
The number of times x appears in the list.
list.sort()
Sort the items of the list, in place.
list.reverse()
Reverse the elements of the list, in place.
list.pop()
Remove the last element of the list, in place.
Table 3.1: Methods of the datatype list
There are two ways list methods can act:
They can directly alter the list, that is, in-place operations.
They produce a new object.
In–place operations
All methods that result in a list are in-place operating methods, for example, reverse:
L = [1, 2, 3]
L.reverse() # the list
L is now reversed
L # [3, 2, 1]
Be aware of in-place operations. One might be tempted to write:
L=[3, 4, 4, 5]
newL = L.sort()
This is correct Python. But it results in a possibly unintended alternation of L in a variable
newL having the value None. The reason is that sort operates in-place.
[ 52 ]
Container Types
Here we demonstrate in-place operating methods:
L = [0, 1, 2, 3, 4]
L.append(5) # [0, 1, 2,
L.reverse() # [5, 4, 3,
L.sort() # [0, 1, 2, 3,
L.remove(0) # [1, 2, 3,
L.pop() # [1, 2, 3, 4]
L.pop() # [1, 2, 3]
L.extend(['a','b','c'])
3,
2,
4,
4,
4, 5]
1, 0]
5]
5]
# [1, 2, 3, 'a', 'b', 'c']
L is altered. The count method is an example of a method that generates a new object:
L.count(2) # returns 1
Merging lists – zip
A particularly useful function for lists is zip. It can be used to merge two given lists into a
new list by pairing the elements of the original lists. The result is a list of tuples (refer
section Tuples for more information):
ind = [0,1,2,3,4]
color = ["red", "green", "blue", "alpha"]
list(zip(color,ind)) # gives [('red', 0), ('green', 1),
('blue', 2), ('alpha', 3)]
This example also demonstrates what happens if the lists have different lengths. The length
of the zipped list is the shorter of the two input lists.
zip creates a special iterable object that can be turned into a list by applying the list
function, as in the preceding example. Refer to section Iterators in Chapter 9, Iterating, for
more details on iterable objects.
List comprehension
A convenient way to build up lists is by using the list comprehension construct, possibly
with a condition inside. The syntax of a list comprehension is:
[ for in ]
[ 53 ]
Container Types
or more generally:
[ for in if ]
Here is an example:
L = [2, 3, 10, 1, 5]
L2 = [x*2 for x in L] # [4, 6, 20, 2, 10]
L3 = [x*2 for x in L if 4 < x <= 10] # [20, 10]
It is possible to have several for loops inside a list comprehension:
M = [[1,2,3],[4,5,6]]
flat = [M[i][j] for i in range(2) for j in range(3)]
# returns [1, 2, 3, 4, 5, 6]
This is of particular interest when dealing with arrays.
Set notation
List comprehension is closely related to the mathematical notation for sets.
Compare:
and L2 = [2*x for x in L].
One big difference though, is that lists are ordered while sets aren't (Refer,
section Sets for more information).
Arrays
The NumPy package offers arrays, which are container structures for manipulating vectors,
matrices, or even higher order tensors in mathematics. In this section, we point out the
similarities between arrays and lists. But arrays deserve a broader presentation, which will
be given in Chapter 4, Linear Algebra – Arrays, and Chapter 5, Advanced Array Concepts.
Arrays are constructed from lists by the function array :
v = array([1.,2.,3.])
A = array([[1.,2.,3.],[4.,5.,6.]])
To access an element of a vector, we need one index, while an element of a matrix is
addressed by two indexes:
v[2]
A[1,2]
# returns 3.0
# returns 6.0
[ 54 ]
Container Types
At first glance, arrays are similar to lists, but be aware that they are different in a
fundamental way, which can be explained by the following points:
Access to array data corresponds to that of lists, using square brackets and slices.
They may also be used to alter the array:
M = array([[1.,2.],[3.,4.]])
v = array([1., 2., 3.])
v[0] # 1
v[:2] # array([1.,2.])
M[0,1] # 2
v[:2] = [10, 20] # v is now array([10., 20., 3.])
The number of elements in a vector, or the number of rows of a matrix, is
obtained by the function len :
len(v) # 3
Arrays store only elements of the same numeric type (usually float or complex
but also int). Refer to section Array properties in Chapter 4, Liner Algebra –
Arrays, for more information.
The operations +, *, /, and - are all elementwise. The dot function and, in
Python versions ≥ 3.5, the infix operator @ are used for the scalar product and the
corresponding matrix operations.
Unlike lists, there is no append method for arrays. Nevertheless, there are special
methods to construct arrays by stacking smaller size arrays (Refer to section
Stacking in Chapter 4, Linear Algebra – Arrays, for more information.). A related
point is that arrays are not elastic as lists; one cannot use slices to change their
length.
Vector slices are views; that is, they may be used to modify the original array.
Refer to section Array views and copies in Chapter 5, Advanced Array Concepts, for
more information.
[ 55 ]
Container Types
Tuples
A tuple is an immutable list. Immutable means that it cannot be modified. A tuple is just a
comma-separated sequence of objects (a list without brackets). To increase readability, one
often encloses a tuple in a pair of parentheses:
my_tuple = 1, 2, 3
# our first tuple
my_tuple = (1, 2, 3)
# the same
my_tuple = 1, 2, 3,
# again the same
len(my_tuple) # 3, same as for lists
my_tuple[0] = 'a'
# error! tuples are immutable
The comma indicates that the object is a tuple:
singleton = 1,
len(singleton)
# note the comma
# 1
Tuples are useful when a group of values goes together; for example, they are used to
return multiple values from functions (refer to section Returns Values in Chapter 7,
Functions. One may assign several variables at once by unpacking a list or tuple:
a, b = 0, 1 #
a, b = [0, 1]
(a, b) = 0, 1
[a,b] = [0,1]
a
#
#
#
gets 0 and b gets 1
exactly the same effect
same
same thing
The swapping trick
Use packing and unpacking to swap the contents of two variables:
a, b = b, a
To summarize:
Tuples are nothing other than immutable lists with a notation without brackets.
In most cases, lists may be used instead of tuples.
The notation without parentheses is convenient but dangerous. You should use
parentheses when you are not sure:
a, b = b, a # the swap trick; equivalent to:
(a, b) = (b, a)
# but
1, 2 == 3, 4 # returns (1, False, 4)
(1, 2) == (3, 4) # returns False
[ 56 ]
Container Types
Dictionaries
Lists, tuples, and arrays are ordered sets of objects. The individual objects are inserted,
accessed, and processed according to their place in the list. On the other hand, dictionaries
are unordered sets of pairs. One accesses dictionary data by keys.
Creating and altering dictionaries
For example, we may create a dictionary containing the data of a rigid body in mechanics,
as follows:
truck_wheel = {'name':'wheel','mass':5.7,
'Ix':20.0,'Iy':1.,'Iz':17.,
'center of mass':[0.,0.,0.]}
A key/data pair is indicated by a colon, :. These pairs are comma separated and listed
inside a pair of curly brackets, {}.
Individual elements are accessed by their keys:
truck_wheel['name']
truck_wheel['mass']
# returns 'wheel'
# returns 5.7
New objects are added to the dictionary by creating a new key:
truck_wheel['Ixy'] = 0.0
Dictionaries are also used to provide parameters to a function (refer to section Parameters
and arguments in Chapter 7, Functions, for further information). Keys in a dictionary can be,
among others, strings, functions, tuples with immutable elements, and classes. Keys cannot
be lists or arrays. The command dict generates a dictionary from a list with key/value
pairs:
truck_wheel = dict([('name','wheel'),('mass',5.7),('Ix',20.0),
('Iy',1.), ('Iz',17.),
('center of mass',[0.,0.,0.])])
The zip function may come in handy in this context (refer to section Merging List).
[ 57 ]
Container Types
Looping over dictionaries
There are mainly three ways to loop over dictionaries:
By keys:
for key in truck_wheel.keys():
print(key) # prints (in any order) 'Ix', 'Iy', 'name',...
or equivalently:
for key in truck_wheel:
print(key) # prints (in any order) 'Ix', 'Iy', 'name',...
By value:
for value in truck_wheel.value():
print(value)
# prints (in any order) 1.0, 20.0, 17.0, 'wheel', ...
By item, that is, key/value pairs:
for item in truck_wheel.items():
print(item)
# prints (in any order) ('Iy', 1.0), ('Ix, 20.0),...
Please, refer to section Shelves in Chapter 12, Input and Output, for a special dictionary
object for file access.
Sets
Sets are containers that share properties and operations with sets in mathematics. A
mathematical set is a collection of distinct objects. Here are some mathematical set
expressions:
[ 58 ]
Container Types
And their Python counterparts:
A
B
C
D
E
5
= {1,2,3,4}
= {5}
= A.union(B)
# returns set([1,2,3,4,5])
= A.intersection(C)
# returns set([1,2,3,4])
= C.difference(A)
# returns set([5])
in C
# returns True
Sets contain an element only once, corresponding to the aforementioned definition:
A = {1,2,3,3,3}
B = {1,2,3}
A == B # returns True
And a set is unordered; that is, the order of the elements in the set is not defined:
A = {1,2,3}
B = {1,3,2}
A == B # returns True
Sets in Python can contain all kinds of hashable objects, that is, numeric objects, strings, and
Booleans.
There are union and intersection methods:
A={1,2,3,4}
A.union({5})
A.intersection({2,4,6}) # returns set([2, 4])
Also, sets can be compared using the methods issubset and issuperset :
{2,4}.issubset({1,2,3,4,5}) # returns True
{1,2,3,4,5}.issuperset({2,4}) # returns True
Empty set
An empty set is defined in Python by empty_set=set([]) and not by {},
which would define an empty dictionary!
[ 59 ]
Container Types
Container conversions
We summarize in the following Table 3.2 the most important properties of the container
types presented so far. Arrays will be treated in Chapter 4, Linear Algebra – Arrays.
Table 3.2 : Container Types
As you can see in the previous table, there is a difference in accessing container elements,
and sets and dictionaries are not ordered.
Due to the different properties of the various container types, we frequently convert one
type to another:
[ 60 ]
Container Types
Type checking
The direct way to see the type of a variable is to use the type command:
label = 'local error'
type(label) # returns str
x = [1, 2] # list
type(x) # returns list
However, if you want to test for a variable to be of a certain type, you should use
isinstance (instead of comparing the types with type):
isinstance(x, list) # True
The reason for using isinstance becomes apparent after having read Chapter 8,
Classes, and in particular the concept of subclassing and inheritance in section Subclassing
and Inheritance in Chapter 8, Classes. In short, often different types share some common
properties with some basic type. The classical example is the type bool, which is derived by
subclassing from the more general type int. In this situation, we see how the command
isinstance can be used in a more general way:
test = True
isinstance(test, bool) # True
isinstance(test, int) # True
type(test) == int # False
type(test) == bool # True
So, in order to make sure that the variable test is as good as an integer (the particular type
may be irrelevant), you should check that it is an instance of integer:
if isinstance(test, int):
print("The variable is an integer")
Type checking
Python is not a typed language. What that means is that objects are
identified by what they can do rather than what they are. For instance, if
you have a string manipulating function that acts on an object by using
the len method, then your function will probably be useful for any objects
implementing that method.
So far, we have come across different types: float, int, bool, complex, list, tuple,
module, function, str, dict, and array.
[ 61 ]
Container Types
Summary
In this chapter, you learned to work with container types, mainly lists. It is important to
know how to fill these containers and how to access their content. We saw that there is
access by position or by keyword.
We will meet the important concept of slicing again in the next chapter on arrays. These are
containers specially designed for mathematical operations.
Exercises
Ex. 1 → Execute the following statements:
L = [1, 2]
L3 = 3*L
1. What is the content of L3?
2. Try to predict the outcome of the following commands:
L3[0]
L3[-1]
L3[10]
3. What does the following command do?
L4 = [k**2 for k in L3]
4. Concatenate L3 and L4 to a new list L5.
Ex. 2 → Use the range command and a list comprehension to generate a list with 100
equidistantly spaced values between 0 and 1.
Ex. 3 → Assume that the following signal is stored in a list:
L = [0,1,2,1,0,-1,-2,-1,0]
[ 62 ]
Container Types
What is the outcome of:
L[0]
L[-1]
L[:-1]
L + L[1:-1] + L
L[2:2] = [-3]
L[3:4] = []
L[2:5] = [-5]
Do this exercise by inspection only, that is, without using your Python Shell.
Ex. 4 → Consider the Python statements:
L = [n-m/2 for n in range(m)]
ans = 1 + L[0] + L[-1]
and assume that the variable m has been previously assigned an integer value. What is the
value of ans? Answer this question without executing the statements in Python.
Ex. 5 → Consider the recursion formula:
with n = 0,…, 1000, h= 1/1000, and a = -0.5.
0
1. Create a list u. Store in its first three elements e , eha, and e2ha. These represent the
starting values u0, u1, and u2 in the given formula. Build up the complete list from
the recursion formula.
2. Construct a second list, td, in which you store the values nh, with n = 0, …, 1000.
Plot td versus u (refer section Basic plotting in Chapter 6, Plotting, for more
information). Make a second plot in which you plot the difference, that is, |eat –
un|, where tn represents the values inside the vector td . Set axis labels and a title.
n
The recursion is a multistep formula to solve the differential equation u' = au with the initial
value u(0) = u0 = 1. un approximates u(nh) = eanhu0.
[ 63 ]
Container Types
Ex. 6 → Let A and B be sets. The set (A \ B) ∪ (B \ A) is called the symmetric difference of
the two sets. Write a function that performs this operation. Compare your results to the
result of the command:
A.symmetric_difference(B).
Ex. 7 → Verify in Python the statement that the empty set is a subset of any set.
Ex. 8 → Study other operations on sets. You find a complete list of those by using the
command completion feature of IPython. In particular, study the update and
intersection_update methods. What is the difference between intersection and
intersection_update?
[ 64 ]
4
Linear Algebra – Arrays
Linear algebra is one of the essential building blocks of computational mathematics. The
objects of linear algebra are vectors and matrices. The package NumPy includes all the
necessary tools to manipulate those objects.
The first task is to build matrices and vectors, or to alter existing ones by slicing. The other
main task is the dot operation, which embodies most of the linear algebra operations (scalar
product, matrix-vector product, and matrix-matrix product). Finally, various methods are
available to solve linear problems.
Overview of the array type
For the impatient, here is how to use arrays in a nutshell. Be aware though that the behavior
of arrays may be surprising at first, so we encourage you to read on after this introductory
section.
Vectors and matrices
Creating vectors is as simple as using the function array to convert a list to an array:
v = array([1.,2.,3.])
The object v is now a vector that behaves much like a vector in linear algebra. We have
already emphasized the differences with the list object in Python (refer to section Arrays in
Chapter 3, Containers Type). Here are some illustrations of the basic linear algebra
operations on vectors:
# two vectors with three components
v1 = array([1., 2., 3.])
Linear Algebra – Arrays
v2 = array([2, 0, 1.])
# scalar multiplications/divisions
2*v1 # array([2., 4., 6.])
v1/2 # array([0.5, 1., 1.5])
# linear combinations
3*v1 # array([ 3., 6., 9.])
3*v1 + 2*v2 # array([ 7., 6., 11.])
# norm
from scipy.linalg import norm
norm(v1) # 3.7416573867739413
# scalar product
dot(v1, v2) # 5.
v1 @ v2 # 5 ; alternative formulation
Note that all basic arithmetic operations are performed elementwise:
# elementwise operations:
v1 * v2 # array([2., 0., 3.])
v2 / v1 # array([2.,0.,.333333])
v1 - v2 # array([-1., 2., 2.])
v1 + v2 # array([ 3., 2., 4.])
Some functions act elementwise on arrays as well:
cos(v1) # cosine, elementwise: array([ 0.5403,
-0.4161, -0.9899])
This subject will be covered in the section Functions Acting on Arrays.
A matrix is created in a similar way to a vector, but from a list of lists instead:
M = array([[1.,2],[0.,1]])
Vectors are no column – and no row matrices
The n vector, an n × 1, and a 1 × n matrix are three different objects even if
they contain the same data.
To create a row matrix containing the same data as the vector v = array([1., 2., 1.]),
we do this:
R = array([[1.,2.,1.]]) # notice the double brackets:
# this is a matrix
shape(R)
# (1,3): this is a row matrix
[ 66 ]
Linear Algebra – Arrays
The corresponding column matrix is obtained by the method reshape:
C = array([1., 2., 1.]).reshape(3, 1)
shape(C) # (3,1): this is a column matrix
Indexing and slices
Indexing and slicing are similar to that of a list. The main difference is that there may be
several indexes or slices when the array is a matrix. The subject will be covered in depth in
section Array indexing; here, we just give some illustrating examples of indexing and slicing:
v = array([1., 2., 3])
M = array([[1., 2],[3., 4]])
v[0] # works as for lists
v[1:] # array([2., 3.])
M[0, 0] # 1.
M[1:] # returns the matrix array([[3., 4]])
M[1] # returns the vector array([3., 4.])
# access
v[0] # 1.
v[0] = 10
# slices
v[:2] # array([10., 2.])
v[:2] = [0, 1] # now v == array([0., 1., 3.])
v[:2] = [1, 2, 3] # error!
Linear algebra operations
The essential operator that performs most of the usual operations of linear algebra is the
Python function dot. It is used for matrix-vector multiplications:
dot(M, v) # matrix vector multiplication; returns a vector
M @ v # alternative formulation
It may be used to compute a scalar product between two vectors:
dot(v, w) # scalar product; the result is a scalar
v @ w # alternative formulation
[ 67 ]
Linear Algebra – Arrays
Lastly, it is used to compute matrix-matrix products:
dot(M, N) # results in a matrix
M @ N # alternative formulation
Solving a linear system
If A is a matrix and b is a vector, you can solve the linear equation:
Using the solve method, which has this syntax:
from scipy.linalg import solve
x = solve(A, b)
For example, we want to solve:
Here is the solution for the preceding equation:
from scipy.linalg import solve
A = array([[1., 2.], [3., 4.]])
b = array([1., 4.])
x = solve(A, b)
allclose(dot(A, x), b) # True
allclose(A @ x, b) # alternative formulation
The command allclose is used here to compare two vectors. If they are close enough to
each other, this command returns True. Optionally a tolerance value can be set. For more
methods related to linear equations systems, refer to section Linear algebra methods in SciPy.
[ 68 ]
Linear Algebra – Arrays
Mathematical preliminaries
In order to understand how arrays work in NumPy, it is useful to understand the
mathematical parallel between accessing tensor (matrix and vector) elements by indexes
and evaluating mathematical functions by providing arguments. We also cover in this
section the generalization of the dot product as a reduction operator.
Arrays as functions
Arrays may be considered from several different points of view. We believe that the most
fruitful one in order to understand arrays is that of functions of several variables.
For instance, selecting a component of a given vector in ℝn may just be considered a
function from the set of ℕn to ℝ, where we define the set:
Here the set ℕn has n elements. The Python function range generates ℕn.
Selecting an element of a given matrix, on the other hand, is a function of two parameters,
taking its value in ℝ. Picking a particular element of an m × n matrix may thus be
considered a function from ℕm × ℕn to ℝ.
Operations are elementwise
NumPy arrays are essentially treated as mathematical functions. This is in particular true
for operations. Consider two functions, f and g, defined on the same domain and taking real
values. The product f g of those two functions is defined as the pointwise product, that is:
[ 69 ]
Linear Algebra – Arrays
Note that this construction is possible for any operation between two functions. For an
arbitrary operation defined on two scalars, which we denote here by , we could define
as follows:
This innocuous remark allows us to understand NumPy's stance on operations; all
operations are elementwise in arrays. For instance, the product between two matrices m and
n is defined, as with functions, as follows:
Shape and number of dimensions
There is a clear distinction between a:
Scalar: Function with no arguments
Vector: Function with one argument
Matrix: Function with two arguments
Higher order tensor: Function with more than two arguments
In what follows, the number of dimensions is the number of arguments of a function. The
shape corresponds essentially to the domain of definition of a function.
For instance, a vector of size n is a function from the set ℕn to ℝ. As a result, its domain of
definition is ℕn. Its shape is defined as the singleton (n,). Similarly, a matrix of size m × n is a
function defined on ℕm × ℕm. The corresponding shape is simply the pair (m, n). The shape of
an array is obtained by the numpy.shape function, and the number of dimensions by the
numpy.ndim function.
[ 70 ]
Linear Algebra – Arrays
The dot operations
Treating arrays as functions, although very powerful, completely neglects the linear algebra
structures we are familiar with, that is, matrix-vector and matrix-matrix
operations. Fortunately, these linear algebra operations may all be written in a similar
unified form:
The vector-vector operation:
The matrix-vector operation:
The matrix-matrix operation:
The vector-matrix operation:
The essential mathematical concept is that of reduction. For a matrix-vector operation, the
reduction is given by:
In general, a reduction operation defined between two tensors T and U of respective
number of dimensions m and n may be defined as:
[ 71 ]
Linear Algebra – Arrays
Clearly, the shapes of the tensors must be compatible for that operation to make any sense.
This requirement is familiar for matrix-matrix multiplication. The multiplication M N of
matrices M and N only makes sense if the number of columns of M equals the number of
rows of N.
Another consequence of the reduction operation is that it produces a new tensor with m + n
– 2 dimensions. In the following table, we gather the output of the reduction operation for
the familiar cases involving matrices and vectors:
Table 4.1: Output of the reduction operation for the familiar cases involving matrices and
vectors
In Python, all reduction operations are performed using the dot function:
angle = pi/3
M = array([[cos(angle), -sin(angle)],
[sin(angle), cos(angle)]])
v = array([1., 0.])
y = dot(M, v)
As in mathematical textbooks, also in modern Python (Version 3.5 and higher), the dot
product is sometimes preferred to be written in its operator form, dot(M, v), or by using
the more handy infix notation, M @ v. From now on we stick to the operator form; you can
modify the examples if the other form is preferred.
Elementwise versus matrix multiplication
The multiplication operator * is always elementwise. It has nothing to do
with the dot operation. Even if A is a matrix and v is a vector, A*v is still a
legal operation.
The matrix-vector multiplication is performed using the dot function.
Refer to section Broadcasting of Chapter 5, Advanced Array Concepts, for
more information.
[ 72 ]
Linear Algebra – Arrays
The array type
The objects used to manipulate vectors, matrices, and more general tensors in NumPy are
called arrays. In this section, we examine their essential properties, how to create them, and
how to access their information.
Array properties
Arrays are essentially characterized by three properties, which is given in the following
table (Table 4.2):
Name
Description
shape
It describes how the data should be interpreted, as a vector, a matrix or as a higher
order tensor, and it gives the corresponding dimension. It is accessed with the
shape attribute.
dtype
It gives the type of the underlying data (float, complex, integer, and so on).
strides This attribute specifies in which order the data should be read. For instance, a
matrix could be stored in memory contiguously column by column (the FORTRAN
convention), or row by row (the C convention). The attribute is a tuple with the
numbers of bytes that have to be skipped in memory to reach the next row and the
number of bytes to be skipped to reach the next column. The strides attribute
even allows for a more flexible interpretation of the data in memory, which is what
makes array views possible.
Table 4.2 : Properties of Arrays
Consider the following array:
A = array([[1, 2, 3], [3, 4, 6]])
A.shape
# (2, 3)
A.dtype
# dtype('int64')
A.strides # (24, 8)
Its elements have type 'int64'; that is, they use 64 bits or 8 bytes in memory. The complete
array is stored in memory row-wise. The distance from A[0, 0] to the first element in the
next row A[1,0] is thus 24 bytes (three matrix elements) in memory. Correspondingly, the
distance in memory between A[0,0] and A[0,1] is 8 bytes (one matrix element). These
values are stored in the attribute strides .
[ 73 ]
Linear Algebra – Arrays
Creating arrays from lists
The general syntax to create an array is the function array . The syntax to create a real
vector would be:
V = array([1., 2., 1.], dtype=float)
To create a complex vector with the same data:
V = array([1., 2., 1.], dtype=complex)
When no type is specified, the type is guessed. The array function chooses the type that
allows storing of all the specified values:
V = array([1, 2]) # [1, 2] is a list of integers
V.dtype # int
V = array([1., 2]) # [1., 2] mix float/integer
V.dtype # float
V = array([1. + 0j, 2.]) # mix float/complex
V.dtype # complex
Silent type conversion
NumPy silently casts floats into integers, which might give unexpected results:
a = array([1, 2, 3])
a[0] = 0.5
a # now: array([0, 2, 3])
The same often unexpected array type casting happens from complex to float.
Array and Python parentheses
As we have noticed in section Program and program flow in Chapter 1, Getting Started,
Python allows a line break when some opening brace or parenthesis is not closed. This
allows a convenient syntax for array creation, which makes it more pleasing to the human
eye:
# the identity matrix in 2D
Id = array([[1., 0.], [0., 1.]])
# Python allows this:
Id = array([[1., 0.],
[0., 1.]])
# which is more readable
[ 74 ]
Linear Algebra – Arrays
Accessing array entries
Array entries are accessed by indexes. In contrast to vector coefficients two indexes are
needed to access matrix coefficients. These are given in one pair of brackets. This
distinguishes the array syntax from a list of lists. There, two pairs of brackets are needed to
access elements.
M = array([[1., 2.],[3., 4.]])
M[0, 0] # first row, first column: 1.
M[-1, 0] # last row, first column: 3.
Basic array slicing
Slices are similar to those of lists except that there might now be in more than one
dimension:
M[i,:] is a vector filled by the row i of M.
M[:,j] is a vector filled by the column i of M.
M[2:4,:] is a slice of 2:4 on the rows only.
M[2:4,1:4] is a slice on rows and columns.
The result of matrix slicing is given in the following figure (Figure 4.1):
Figure 4.1: The result of matrix slicing
[ 75 ]
Linear Algebra – Arrays
Omitting a dimension
If you omit an index or a slice, NumPy assumes you are taking rows only.
M[3] is a vector that is a view on the third row of M and M[1:3] is a
matrix that is a view on the second and third rows of M.
Changing the elements of a slice affects the entire array:
v = array([1., 2., 3.])
v1 = v[:2] # v1 is array([1., 2.])
v1[0] = 0. # if v1 is changed ...
v # ... v is changed too: array([0., 2., 3.])
General slicing rules are given in the following table (Table 4.3):
Table 4.3: General Slicing Rules
The results of slicing operations for an array M of shape (4, 4) are given in the following table
(Table 4.4):
Table 4.4: Result of slicing operation for an array M of shape (4,4)
[ 76 ]
Linear Algebra – Arrays
Altering an array using slices
You may alter an array using slices or by direct access. The following changes only one
element in a 5 × 3 matrix M:
M[1, 3] = 2.0 # scalar
But we may change one full row of the matrix:
M[2, :] = [1., 2., 3.] # vector
We may also replace a full submatrix:
M[1:3, :] = array([[1., 2., 3.],[-1.,-2., -3.]])
There is a distinction between a column matrix and a vector. The
following assignment with a column matrix returns no error
M[1:4, 2:3] = array([[1.],[0.],[-1.0]])
while the assignment with a vector returns a Value Error
M[1:4, 2:3] = array([1., 0., -1.0]) # error
The general slicing rules are shown in Table 4.2. The matrices and vectors in the preceding
examples must have the right size to fit into matrix M. You may also make use of the
broadcasting rules (for more information, refer to section Broadcasting of Chapter 5,
Advanced Array Concepts) to determine the allowed size of the replacement arrays. If the
replacement array does not have the right shape, a ValueError exception will be raised.
Functions to construct arrays
The usual way to set up an array is via a list. But there are also a couple of convenient
methods for generating special arrays, which are given in the following table (Table 4.5):
Methods
Shape Generates
zeros((n,m))
(n,m)
Matrix filled with zeros
ones((n,m))
(n,m)
Matrix filled with ones
diag(v,k)
(n,n)
(Sub-, super-) diagonal matrix from a vector v
random.rand(n,m)
(n,m)
Matrix filled with uniformly distributed random numbers in
(0,1)
[ 77 ]
Linear Algebra – Arrays
arange(n)
(n,)
First n integers
linspace(a,b,n)
(n,)
Vector with n equispaced points between a and b
Table 4.5: Commands to create arrays
These commands may take additional arguments. In particular, the commands zeros,
ones, and arange take dtype as an optional argument. The default type is float, except
for arange. There are also methods such as zeros_like and ones_like, which are slight
variants of the preceding ones. For instance, the zeros_like(A) method is equivalent to
zeros(shape(A)).
Here is the identity function, which constructs an identity matrix of a given size:
I = identity(3)
The command is identical to:
I = array([[ 1., 0., 0.],
[ 0., 1., 0.],
[ 0., 0., 1.]])
Accessing and changing the shape
The number of dimensions is what distinguishes a vector from a matrix. The shape is what
distinguishes vectors of different sizes, or matrices of different sizes. In this section, we
examine how to obtain and change the shape of an array.
The shape function
The shape of a matrix is the tuple of its dimensions. The shape of an n × m matrix is the
tuple (n, m). It can be obtained by the shape function:
M = identity(3)
shape(M) # (3, 3)
For a vector, the shape is a singleton containing the length of that vector:
v = array([1., 2., 1., 4.])
shape(v) # (4,) <- singleton (1-tuple)
[ 78 ]
Linear Algebra – Arrays
An alternative is to use the array attribute shape, which gives the same result:
M = array([[1.,2.]])
shape(M) # (1,2)
M.shape # (1,2)
However, the advantage of using shape as a function is that this function may be used on
scalars and lists as well. This may come in handy when code is supposed to work with both
scalars and arrays:
shape(1.) # ()
shape([1,2]) # (2,)
shape([[1,2]]) # (1,2)
Number of dimensions
The number of dimensions of an array is obtained with the function numpy.ndim or using
the array attribute ndarray.ndim :
ndim(A) # 2
A.ndim # 2
Note that the number of dimensions, given by the function ndim, of a tensor T (a vector,
matrix, or higher order tensor) is always equal to the length of its shape:
T = zeros((2,2,3)) # tensor of shape (2,2,3); three dimensions
ndim(T) # 3
len(shape(T)) # 3
Reshape
The method reshape gives a new view of the array, with a new shape, without copying the
data:
v = array([0,1,2,3,4,5])
M = v.reshape(2,3)
shape(M) # returns (2,3)
M[0,0] = 10 # now v[0] is 10
[ 79 ]
Linear Algebra – Arrays
Reshape does not copy
Reshape does not create a new array. It rather gives a new view on the
existing array. In the preceding example, changing one element of M
would automatically result in a change in the corresponding element in v.
When this behavior is not acceptable, you need to copy the data.
The various effects of the reshape method on an array defined by arange(6) are given in
the following figure :
Figure 4.2: The various effects of the reshape method on an array defined by arange(6)
If one tries to reshape an array with a shape that does not multiply to the original shape, an
error is raised:
ValueError: total size of new array must be unchanged.
[ 80 ]
Linear Algebra – Arrays
Sometimes, it is convenient to specify only one shape parameter and let Python determine
the other in such a way that it multiplies to the original shape. This is done by setting the
free shape parameter -1:
v = array([1, 2, 3, 4, 5, 6, 7, 8])
M = v.reshape(2, -1)
shape(M) # returns (2, 4)
M = v.reshape(-1, 2)
shape(M) # returns (4,2 )
M = v.reshape(3,- 1) # returns error
Transpose
A special form of reshaping is transposing. It just switches the two shape elements of the
matrix. The transpose of a matrix A is a matrix B such that:
Which is resolved in the following way:
A = ...
shape(A) # 3,4
B = A.T # A transpose
shape(B) # 4,3
Transpose does not copy
Transposition is very similar to reshaping. In particular, it does not copy
the data either and just returns a view on the same array:
A= array([[ 1., 2.],[ 3., 4.]])
B=A.T
A[1,1]=5.
B[1,1] # 5
Transposing a vector makes no sense since vectors are tensors of one dimension, that is,
functions of one variable. NumPy will, however, comply and return exactly the same object:
v = array([1., 2., 3.])
v.T # exactly the same vector!
[ 81 ]
Linear Algebra – Arrays
What you have in mind when you want to transpose a vector is probably to create a row or
column matrix. This is done using reshape:
v.reshape(-1, 1) # column matrix containing v
v.reshape(1, -1) # row matrix containing v
Stacking
The universal method to build matrices from a couple of (matching) submatrices is
concatenate. Its syntax is:
concatenate((a1, a2, ...), axis = 0)
This command stacks the submatrices vertically (on top of each other) when axis=0 is
specified. With the axis=1 argument, they are stacked horizontally, and this generalizes
according to arrays with more dimensions. This function is called by several convenient
functions, as follows:
hstack: Used to stack matrices horizontally
vstack: Used to stack matrices vertically
columnstack: Used to stack vectors in columns
Stacking vectors
One may stack vectors row-wise or column-wise using vstack and column_stack, as
illustrated in the following figure:
[ 82 ]
Linear Algebra – Arrays
hstack would produce the concatenation of v1 and v2.
Let us consider the symplectic permutation as an example for vector stacking:
We have a vector of size 2n. We want to perform a symplectic transformation of a vector
with an even number of components, that is, exchange the first half with the second half of
the vector with sign change:
This operation is resolved in Python as follows:
# v is supposed to have an even length.
def symp(v):
n = len(v) // 2 # use the integer division //
return hstack([v[-n:], -v[:n]])
Functions acting on arrays
There are different types of functions acting on arrays. Some act elementwise, and they
return an array of the same shape. Those are called universal functions. Other array
functions return an array of a different shape.
Universal functions
Universal functions are functions that act elementwise on arrays. They thus have an output
array that has the same shape as the input array. These functions allow us to compute the
result of a scalar function on a whole array at once.
Built-in universal functions
A typical example is the cos function (the one provided by NumPy):
cos(pi) # -1
cos(array([[0, pi/2, pi]])) # array([[1, 0, -1]])
[ 83 ]
Linear Algebra – Arrays
Note that universal functions work on arrays in a componentwise manner. This is also true
for operators, such as multiplication or exponent:
2 * array([2, 4]) # array([4, 8])
array([1, 2]) * array([1, 8]) # array([1, 16])
array([1, 2])**2 # array([1, 4])
2**array([1, 2]) # array([1, 4])
array([1, 2])**array([1, 2]) # array([1, 4])
Create universal functions
Your function will automatically be universal if you use only universal functions in it. If,
however, your function uses functions that are not universal, you might get scalar results,
or even an error, when trying to apply them on an array:
def const(x):
return 1
const(array([0, 2])) # returns 1 instead of array([1, 1])
Another example is the following:
def heaviside(x):
if x >= 0:
return 1.
else:
return 0.
heaviside(array([-1, 2])) # error
The expected behaviour would be that the heaviside function applied to a vector [a, b]
would return [heaviside(a), heaviside(b)]. Alas, this does not work because the
function always returns a scalar, no matter the size of the input argument. Besides, using
the function with an array input would raise an exception. The NumPy function
vectorize allows us to quickly solve this problem:
vheaviside = vectorize(heaviside)
vheaviside(array([-1, 2])) # array([0, 1]) as expected
[ 84 ]
Linear Algebra – Arrays
A typical application of this method is its use when plotting a function:
xvals = linspace(-1, 1, 100)
plot(xvals, vectorize(heaviside)(xvals))
axis([-1.5, 1.5, -0.5, 1.5])
The following graph shows the heaviside function:
The vectorize function does not improve performance. It provides only
a convenient way to quickly transform a function, so that it operates
elementwise on list and arrays.
[ 85 ]
Linear Algebra – Arrays
Array functions
There are a number of functions acting on arrays that do not act componentwise. Examples
of such functions are max, min, and sum. These functions may operate on the whole matrix,
row-wise, or column-wise. When no argument is provided, they act on the whole matrix.
Suppose A is the following matrix:
The sum function acting on that matrix returns a scalar:
sum(A) # 36
The command has an optional parameter, axis . It allows us to choose along which axis to
perform the operation. For instance, if the axis is 0, it means that the sum should be
computed along the first axis. The sum along axis 0 of an array of shape (m, n) will be a
vector of length n.
Suppose we compute the sum of A along the axis 0:
sum(A, axis=0) # array([ 6, 8, 10, 12])
This amounts to computing the sum on the columns:
The result is a vector:
Now suppose we compute the sum along the axis 1:
A.sum(axis=1) # array([10, 26])
[ 86 ]
Linear Algebra – Arrays
This amounts to computing the sum on the rows:
The result is a vector:
Linear algebra methods in SciPy
SciPy offers a large range of methods from numerical linear algebra in its scipy.linalg
module. Many of these methods are Python wrapping programs from LAPACK, a collection
of well-approved FORTRAN subroutines used to solve linear equation systems and
eigenvalue problems. Linear algebra methods are the core of any method in scientific
computing, and the fact that SciPy uses wrappers instead of pure Python code makes these
central methods extremely fast. We present in detail here how two linear algebra problems
are solved with SciPy to give you a flavour of this module.
Solving several linear equation systems with LU
Let A be an n × n matrix and b1, b2, …, bk be a sequence of n-vectors. We consider the
problem to find n vectors xi such that:
We assume that the vectors bi are not known simultaneously. In particular, it is quite a
common situation that the ith problem has to be solved before bi+1 becomes available.
LU factorization is a way to organize the classical Gauss elimination method in such a way
that the computation is done in two steps:
A factorization step of the matrix A to get matrices in triangular form
A relatively cheap backward and forward elimination step that works on the bi's
and benefits from the more time-consuming factorization step
[ 87 ]
Linear Algebra – Arrays
The method also uses the fact that if P is a permutation matrix such that PA is the original
matrix with its rows permuted.
The two systems
have the same solution.
LU factorization finds a permutation matrix P, a lower triangular matrix L, and an upper
triangular matrix U such that:
.
Such a factorization always exists. Furthermore, L can be determined in such a way that Lii =
1. Thus, the essential data from L that has to be stored is Lij with i > j. Consequently, L and U
can be stored together in an n × n array, while the information about the permutation matrix
P just requires an n integer vector – the pivot vector.
In SciPy, there are two methods to compute the LU factorization. The standard one is
scipy.linalg.lu, which returns the three matrices L, U, and P. The other method
islu_factor. That is the method we describe here, because it will be conveniently used
later in combination with lu_solve:
import scipy.linalg as sl
[LU,piv] = sl.lu_factor(A)
Here, the A matrix is factorized and an array with the information about L and U is returned,
together with the pivot vector. With this information, the system can be solved by
performing row interchanges of the vectors bi according to the information stored in the
pivot vector, backward substitution using U, and finally forward substitution using L. This
is bundled in Python, in the lu_solve method. The following code snippet shows how the
system Axi = bi is solved once the LU factorization is performed and its results stored in the
tuple (LU, piv):
import scipy.linalg as sl
xi = sl.lu_solve((LU, piv), bi)
[ 88 ]
Linear Algebra – Arrays
Solving a least square problem with SVD
A linear equation system Ax = b, with A being an m × n matrix and m > n, is called an
overdetermined linear system. In general it has no classical solution and one seeks a vector
x* ℝn with the property:
Here,
denotes the Euclidean vector norm
.
This problem is called a least square problem. A stable method to solve it is based on
factorizing A = UΣVT, with U being a m × m orthogonal matrix, V a n × n orthogonal matrix,
and Σ = (σij) an m × n matrix with the property σij = 0 for all i ≠j. This factorization is called a
singular value decomposition (SVD).
We write,
with a diagonal n × n matrix Σ1. If we assume that A has full rank, then Σ1 is invertible and
it can be shown that,
. If we split U = [U1 U2] with U1 being an m × n submatrix,
then the preceding equation can be simplified to
.
SciPy provides a function called svd, which we use to solve this task:
import scipy.linalg as sl
[U1, Sigma_1, VT] = sl.svd(A, full_matrices = False,
compute_uv = True)
xast = dot(VT.T, dot(U1.T, b) / Sigma_1)
r = dot(A, xast) - b # computes the residual
nr = sl.norm(r, 2) # computes the Euclidean norm of r
[ 89 ]
Linear Algebra – Arrays
The keyword full_matrices says that only the portion U1 of U needs to be computed. As
one often uses svd to compute only singular values, σii, we have to explicitly demand the
computation of U and V by using the keyword compute_uv. The SciPy function
scipy.linalg.lstsq solves the least square problem similarly by using a singular value
decomposition.
More methods
In the examples so far, you met a couple of methods for computational tasks in linear
algebra, for example, solve. Most common methods are available after the
command import scipy.linalg as sl is executed. We refer to their documentation for
further reference. Some linear algebra functions of the scipy.linalg module are given in
the following table (Table 4.6):
Methods
Description (matrix methods)
sl.det
Determinant of a matrix
sl.eig
Eigenvalues and eigenvectors of a matrix
sl.inv
Matrix inverse
sl.pinv
Matrix pseudoinverse
sl.norm
Matrix or vector norm
sl.svd
Singular value decomposition
sl.lu
LU decomposition
sl.qr
QR decomposition
sl.cholesky
Cholesky decomposition
sl.solve
Solution of a general or symmetric linear system: Ax = b
sl.solve.banded The same for banded matrices
sl.lstsq
Least squares solution
Table 4.6: Linear algebra functions of the scipy.linalg module
Execute import scipy.linalg as sl first.
[ 90 ]
Linear Algebra – Arrays
Summary
In this chapter, we worked with the most important objects in linear algebra – vectors and
matrices. For this, we learned how to define arrays and we met important array methods. A
smaller section demonstrated how to use modules from scipy.linalg to solve central
tasks in linear algebra.
Exercises
Ex. 1 → Consider a 4 × 3 matrix M:
1. Construct this matrix in Python using the function array .
2. Construct the same matrix using the function arange followed by a suitable
reshape.
3. What is the result of the expression M[2,:] ? What is the result of the similar
expression M[2:]?
Ex. 2 → Given a vector x, construct in Python the following matrix:
[ 91 ]
Linear Algebra – Arrays
Here, xi are the components of the vector x (numbered from zero). Given a vector y, solve in
Python the linear equation system Va = y. Let the components of a be denoted by ai, i = 0, …,
5. Write a function poly, which has a and z as input and which computes the polynomial:
Plot this polynomial and depict in the same plot the points (xi, yi) as small stars. Try your
code with the vectors:
x = (0.0, 0.5, 1.o, 1.5, 2.0, 2.5)
y = (-2.0, 0.5, -2.0, 1.0, -0.5, 1.0)
Ex. 3 → The matrix V in Ex. 2 is called a Vandermonde matrix. It can be set up in Python
directly by the command vander. Evaluating a polynomial defined by a coefficient vector
can be done with the Python command polyval. Repeat Ex. 2 by using these commands.
Ex. 4 → Let u be a one dimensional array. Construct another array ξ with values ξi = (u1 + ui+1
+ ui+2)/3. In statistics, this array is called the moving average of u. In approximation theory, it
plays the role as the Greville abscissae of cubic splines. Try to avoid the use of for loops in
your script.
Ex. 5 →
1.
2.
3.
4.
Construct from the matrix V given in Ex. 2 a matrix A by deleting V's first column.
Form the matrix B = (AT A)-1 AT.
Compute c = B y with y from Ex. 2.
Use c and polyval to plot the polynomial defined by c. Plot in the same picture
again the points (xi, yi).
Ex. 6 → Ex. 5 describes the least square method. Repeat that exercise but use SciPy's
scipy.linalg.lstsq method instead.
[ 92 ]
Linear Algebra – Arrays
Ex. 7 → Let v be a vector written in its coordinate form as a 3 × 1 matrix
T
[1 -1 1] . Construct the projection matrices:
.
Show experimentally that v is an eigenvector for both matrices P and Q. What are the
corresponding eigenvalues?
Ex. 8 → In numerical linear algebra the m × m matrix A with the property
is used as an example for an extreme growth-factor, when performing LU factorization.
Set up this matrix in Python for various m, compute its LU factorization using the
command scipy.linalg.lu and derive experimentally a statement about the growth
factor
in relation to m.
[ 93 ]
5
Advanced Array Concepts
In this chapter, we will explain some more advanced aspects of arrays. First, we will cover
the notion of an array view, followed by Boolean arrays and how to compare arrays. We
briefly describe indexing and vectorization, explain sparse arrays, and some special topics
such as broadcasting.
Array views and copies
In order to control precisely how memory is used, NumPy offers the concept of view of an
array. Views are smaller arrays that share the same data as a larger array. This works just
like a reference to one single object (refer to section Basic Types in Chapter 1, Getting
Started).
Array views
The simplest example of a view is given by a slice of an array:
M = array([[1.,2.],[3.,4.]])
v = M[0,:] # first row of M
The preceding slice is a view of M. It shares the same data as M. Modifying v will modify M as
well:
v[-1] = 0.
v # array([[1.,0.]])
M # array([[1.,0.],[3.,4.]]) # M is modified as well
Advanced Array Concepts
It is possible to access the object that owns the data using the array attribute base:
v.base # array([[1.,0.],[3.,4.]])
v.base is M # True
If an array owns its data, the attribute base is none :
M.base # None
Slices as views
There are precise rules on which slices will return views and which ones will return copies.
Only basic slices (mainly index expressions with :) return views, whereas any advanced
selections (such as slicing with a Boolean) will return a copy of the data. For instance, it is
possible to create new matrices by indexing with lists (or arrays):
a = arange(4) # array([0.,1.,2.,3.])
b = a[[2,3]] # the index is a list [2,3]
b # array([2.,3.])
b.base is None # True, the data was copied
c = a[1:3]
c.base is None # False, this is just a view
In the preceding example, the array b is not a view, whereas the array c, obtained with a
simpler slice, is a view.
There is an especially simple slice of an array that returns a view of the whole array:
N = M[:] # this is a view of the whole array M
Transpose and reshape as views
Some other important operations return views. For instance, transpose returns a view:
M = random.random_sample((3,3))
N = M.T
N.base is M # True
The same applies for all reshaping operations:
v = arange(10)
C = v.reshape(-1,1) # column matrix
C.base is v # True
[ 95 ]
Advanced Array Concepts
Array copy
Sometimes it is necessary to explicitly request that the data be copied. This is simply
achieved with the NumPy function called array:
M = array([[1.,2.],[3.,4.]])
N = array(M.T) # copy of M.T
We may verify that the data has indeed been copied:
N.base is None # True
Comparing arrays
Comparing two arrays is not as simple as it may seem. Consider the following code, which
is intended to check whether two matrices are close to each other:
A = array([0.,0.])
B = array([0.,0.])
if abs(B-A) < 1e-10: # an exception is raised here
print("The two arrays are close enough")
This code raises the exception when the if statement is executed:
ValueError: The truth value of an array with more than one element is
ambiguous. Use a.any() or a.all()
In this section, we explain why this is so and how to remedy this state of affairs.
Boolean arrays
Boolean arrays are useful for advanced array indexing (refer to section Indexing with Boolean
arrays). A Boolean array is simply an array for which the entries have the type bool:
A = array([True,False]) # Boolean array
A.dtype # dtype('bool')
[ 96 ]
Advanced Array Concepts
Any comparison operator acting on arrays will create a Boolean array instead of a simple
Boolean:
M = array([[2, 3],
[1, 4]])
M > 2 # array([[False, True],
# [False, True]])
M == 0 # array([[False, False],
# [False, False]])
N = array([[2, 3],
[0, 0]])
M == N # array([[True, True],
# [False, False]])
...
Note that because array comparison creates Boolean arrays, one cannot use array
comparison directly in conditional statements, for example, if statements. the solution is to
use the methods all and any:
A = array([[1,2],[3,4]])
B = array([[1,2],[3,3]])
A == B # creates array([[True, True], [True, False]])
(A == B).all() # False
(A != B).any() # True
if (abs(B-A) < 1e-10).all():
print("The two arrays are close enough")
Checking for equality
Checking for equality of two float arrays is not straight forward, because two floats may be
very close without being equal. In NumPy, it is possible to check for equality with
allclose. This function checks for equality of two arrays up to a given precision:
data = random.rand(2)*1e-3
small_error = random.rand(2)*1e-16
data == data + small_error # False
allclose(data, data + small_error, rtol=1.e-5, atol=1.e-8)
# True
The tolerance is given in terms of a relative tolerance bound, rtol, and an absolute error
bound, atol. The command allclose is a short form of: (abs(A-B) <
atol+rtol*abs(B)).all().
[ 97 ]
Advanced Array Concepts
Note that allclose can be also applied to scalars:
data = 1e-3
error = 1e-16
data == data + error # False
allclose(data, data + error, rtol=1.e-5, atol=1.e-8)
#True
Boolean operations on arrays
You cannot use and, or, and not on Boolean arrays. Indeed, those operators force the
casting from array to Boolean, which is not permitted. Instead, we can use the operators
given in the following table (Table 5.1) for componentwise logical operations on Boolean
arrays:
Logic operator
Replacement for Boolean arrays
A and B
A & B
A or B
A | B
not A
~ A
Table 5.1 Logical operators and, or and not do not work with arrays.
A = array([True, True, False, False])
B = array([True, False, True, False])
A and B # error!
A & B # array([True, False, False, False])
A | B # array([True, True, True, False])
~A # array([False, False, True, True])
Here is an example usage of logical operators with Boolean arrays:
Suppose that we have a sequence of data that is marred with some measurement error.
Suppose further that we run a regression and it gives us a deviation for each value. We wish
to obtain all the exceptional values and all the values with little deviation that are lower
than a given threshold:
data = linspace(1,100,100) # data
deviation = random.normal(size=100) # the deviations
#don't forget the parentheses in next statement!
exceptional = data[(deviation<-0.5)|(deviation>0.5)]
exceptional = data[abs(deviation)>0.5] # same result
small = data[(abs(deviation)<0.1)&(data<5.)] # small deviation and data
[ 98 ]
Advanced Array Concepts
Array indexing
We have already seen that one may index arrays by combinations of slices and integers, this
is the basic slicing technique. There are, however, many more possibilities, which allow for
a variety of ways to access and modify array elements.
Indexing with Boolean arrays
It is often useful to access and modify only parts of an array, depending on its value. For
instance, one might want to access all the positive elements of an array. This turns out to be
possible using Boolean arrays, which act like masks to select only some elements of an
array. The result of such an indexing is always a vector. For instance, consider the following
example:
B = array([[True, False],
[False, True]])
M = array([[2, 3],
[1, 4]])
M[B] # array([2,4]), a vector
In fact, the M[B] call is equivalent to M.flatten()[B]. One may then replace the resulting
vector by another vector. For instance, one may replace all the elements by zero (refer to
section Broadcasting for more information):
M[B] = 0
M # [[0, 3], [1, 0]]
Or one may replace all the selected values by others:
M[B] = 10, 20
M # [[10, 3], [1, 20]]
By combining the creation of Boolean arrays (M > 2), smart indexing (indexing with
Boolean array), and broadcasting, one may use the following elegant syntax:
M[M>2] = 0
# all the elements > 2 are replaced by 0
The expression broadcasting here refers to the tacit conversion of the scalar 0 to a vector of
an appropriate shape.
[ 99 ]
Advanced Array Concepts
Using where
The command where gives a useful construct that can take a Boolean array as a condition
and either return the indexes of the array elements satisfying the condition or return
different values depending on the values in the Boolean array.
The basic structure is:
where(condition, a, b)
This will return values from a when the condition is True and values from b when it is
False.
For instances consider, a Heaviside function:
The following code implements a Heaviside function:
def H(x):
return where(x < 0, 0, 1)
x = linspace(-1,1,11) # [-1. -0.8 -0.6 -0.4 -0.2 0. 0.2 0.4 0.6 0.8 1. ]
print(H(x))
# [0 0 0 0 0 1 1 1 1 1 1]
The second and third arguments can be either arrays of the same size as the condition (the
Boolean array) or scalars. We give two more example to demonstrated how to manipulate
elements from an array or a scalar depending on a condition:
x = linspace(-4,4,5)
# [-4. -2. 0. 2. 4.]
print(where(x > 0, sqrt(x), 0))
# [ 0.+0.j 0.+0.j 0.+0.j 1.41421356+0.j 2.+0.j ]
print(where(x > 0, 1, -1)) # [-1 -1 -1 1 1]
If the second and third arguments are omitted, then a tuple containing the indexes of the
elements satisfying the condition is returned.
[ 100 ]
Advanced Array Concepts
For example consider the use of where with only one argument in the following code:
a = arange(9)
b = a.reshape((3,3))
print(where(a > 5))
# (array([6, 7, 8]),)
print(where(b > 5))
# (array([2, 2, 2]), array([0, 1, 2]))
Performance and Vectorization
When it comes to performance of your Python code, it often boils down to the difference
between interpreted code and compiled code. Python is an interpreted programming
language and basic Python code is executed directly without any intermediate compilation
to machine code. With a compiled language, the code needs to be translated to machine
instructions before execution.
The benefits of an interpreted language are many but interpreted code cannot compete with
compiled code for speed. To make your code faster, you can write some parts in a compiled
language like FORTRAN, C, or C++. This is what NumPy and SciPy do.
For this reason, it is best to use functions in NumPy and SciPy over interpreted versions
whenever possible. NumPy array operations such as matrix multiplication, matrix-vector
multiplication, matrix factorization, scalar products, and so on are much faster than any
pure Python equivalent. Consider the simple case of scalar products. The scalar product is
much slower than the compiled NumPy function, dot(a,b) (more than 100 times slower
for arrays with about 100 elements):
def my_prod(a,b):
val = 0
for aa,bb in zip(a,b):
val += aa*bb
return val
Measuring the speed of your functions is an important aspect of scientific computing. Refer
to section Measuring execution time in Chapter 13, Testing, for details on measuring
execution times.
[ 101 ]
Advanced Array Concepts
Vectorization
To improve performance, one has to vectorize the code often. Replacing for loops and
other slower parts of the code with NumPy slicing, operations, and functions can give
significant improvements. For example, the simple addition of a scalar to a vector by
iterating over the elements is very slow:
for i in range(len(v)):
w[i] = v[i] + 5
where using NumPy's addition is much faster:
w = v + 5
Using NumPy slicing can also give significant speed improvements over iterating with for
loops. To demonstrate this let us consider forming the average of neighbors in a twodimensional array:
def my_avg(A):
m,n = A.shape
B = A.copy()
for i in range(1,m-1):
for j in range(1,n-1):
B[i,j] = (A[i-1,j] + A[i+1,j] + A[i,j-1] + A[i,j+1])/4
return B
def slicing_avg(A):
A[1:-1,1:-1] = (A[:-2,1:-1] + A[2:,1:-1] +
A[1:-1,:-2] + A[1:-1,2:])/4
return A
These functions both assign each element the average of its four neighbors. The second
version, using slicing, is much faster.
Besides replacing for loops and other slower constructions with NumPy functions, there is
a useful function called vectorize, refer to section Functions acting on arrays in Chapter 4,
Linear Algebra – Arrays. This will take a function and create a vectorized version that applies
the function on all elements of an array using functions wherever possible.
[ 102 ]
Advanced Array Concepts
Consider the following example for vectorizing a function:
def my_func(x):
y = x**3 - 2*x + 5
if y>0.5:
return y-0.5
else:
return 0
Applying this by iterating over an array is very slow:
for i in range(len(v)):
v[i] = my_func(v[i])
Instead, use vectorize to create a new function, like this:
my_vecfunc = vectorize(my_func)
This function can then be applied to the array directly:
v = my_vecfunc(v)
The vectorized option is much faster (around 10 times faster with arrays of length 100).
Broadcasting
Broadcasting in NumPy denotes the ability to guess a common, compatible shape between
two arrays. For instance, when adding a vector (one-dimensional array) and a scalar (zerodimensional array), the scalar is extended to a vector, in order to allow for the addition. The
general mechanism is called broadcasting. We will first review that mechanism from a
mathematical point of view, and then proceed to give the precise rules for broadcasting in
NumPy.
Mathematical view
Broadcasting is often performed in mathematics, mainly implicitly. Examples are
expressions such as f(x) + C or f(x) + g(y). We'll give an explicit description of that technique
in this section.
[ 103 ]
Advanced Array Concepts
We have in mind the very close relationship between functions and NumPy arrays, as
described in section Mathematical preliminaries of Chapter 4, Linear Algebra – Arrays.
Constant functions
One of the most common examples of broadcasting is the addition of a function and a
constant; if C is a scalar, one often writes:
This is an abuse of notation since one should not be able to add functions and constants.
Constants are however implicitly broadcast to functions. The broadcast version of the
constant C is the function defined by:
Now it makes sense to add two functions together:
We are not being pedantic for the sake of it, but because a similar situation may arise for
arrays, as in the following code:
vector = arange(4) # array([0.,1.,2.,3.])
vector + 1.
# array([1.,2.,3.,4.])
In this example, everything happens as if the scalar 1. had been converted to an array of
the same length as vector, that is, array([1.,1.,1.,1.]), and then added to vector.
This example is exceedingly simple, so we proceed to show less obvious situations.
[ 104 ]
Advanced Array Concepts
Functions of several variables
A more intricate example of broadcasting arises when building functions of several
variables. Suppose, for instance, that we were given two functions of one variable, f and g,
and that we want to construct a new function F according to the formula:
This is clearly a valid mathematical definition. We would like to express this definition as
the sum of two functions in two variables defined as
and now we may simply write:
The situation is similar to that arising when adding a column matrix and a row matrix:
C = arange(2).reshape(-1,1) # column
R = arange(2).reshape(1,-1) # row
C + R
# valid addition: array([[0.,1.],[1.,2.]])
This is especially useful when sampling functions of two variables, as shown in section
Typical examples.
General mechanism
We have seen how to add a function and a scalar and how to build a function of two
variables from two functions of one variable. Let us now focus on the general mechanism
that makes this possible. The general mechanism consists of two steps: reshaping and
extending.
[ 105 ]
Advanced Array Concepts
First, the function g is reshaped to a function that takes two arguments. One of these
arguments is a dummy argument, which we take to be zero, as a convention:
Mathematically, the domain of definition of is now
in a way similar to:
Then the function f is reshaped
Now both and take two arguments, although one of them is always zero. We proceed to
the next step, extending. It is the same step that converted a constant into a constant
function (refer to the constant function example).
The function is extended to:
The function is extended to:
Now the function of two variables F, which was sloppily defined by F(x,y) = f(x) + g(y), may
be defined without reference to its arguments:
For example, let us describe the preceding mechanism for constants. A constant is a scalar,
that is, a function of zero arguments. The reshaping step is thus to define the function of one
(empty) variable:
Now the extension step proceeds simply by:
[ 106 ]
Advanced Array Concepts
Conventions
The last ingredient is the convention on how to add the extra arguments to a function, that
is, how the reshaping is automatically performed. By convention, a function is
automatically reshaped by adding zeros on the left.
For example, if a function g of two arguments has to be reshaped to three arguments, the
new function would be defined by:
Broadcasting arrays
We now repeat the observation that arrays are merely functions of several variables (refer to
section Mathematical preliminaries in Chapter 4, Linear Algebra – Arrays). Array broadcasting
thus follows exactly the same procedure as explained above for mathematical functions.
Broadcasting is done automatically in NumPy.
In the following figure (Figure 5.1), we show what happens when adding a matrix of shape
(4, 3) to a matrix of size (1, 3). The second matrix is of the shape (4, 3):
Figure 5.1: Broadcasting between a matrix and a vector.
The broadcasting problem
When NumPy is given two arrays with different shapes, and is asked to perform an
operation that would require the two shapes to be the same, both arrays are broadcast to a
common shape.
[ 107 ]
Advanced Array Concepts
Suppose the two arrays have shapes s1 and s2. This broadcasting is performed in two steps:
1. If the shape s1 is shorter than the shape s2 then ones are added on the left of the
shape s1. This is a reshaping.
2. When the shapes have the same length, the array is extended to match the shape
s2 (if possible).
Suppose we want to add a vector of shape (3, ) to a matrix of shape (4, 3). The vector needs
be broadcast. The first operation is a reshaping; the shape of the vector is converted from (3,
) to (1, 3). The second operation is an extension; the shape is converted from (1, 3) to (4, 3).
For instance, suppose a vector of size n is to be broadcast to the shape (m, n):
1. v is automatically reshaped to (1, n).
2. v is extended to (m, n).
To demonstrate this we consider a matrix defined by:
M = array([[11, 12, 13, 14],
[21, 22, 23, 24],
[31, 32, 33, 34]])
and vector given by:
v = array([100, 200, 300, 400])
Now we may add M and v directly:
M + v # works directly
The result is this matrix:
[ 108 ]
Advanced Array Concepts
Shape mismatch
It is not possible to automatically broadcast a vector v of length n to the shape (n,m). This
is illustrated in the following figure:
The broadcasting will fail, because the shape (n,) may not be automatically broadcast to
the shape (m, n). The solution is to manually reshape v to the shape (n,1). The
broadcasting will now work as usual (by extension only):
M + v.reshape(-1,1)
Here is another example, define a matrix by:
M = array([[11, 12, 13, 14],
[21, 22, 23, 24],
[31, 32, 33, 34]])
and a vector by:
v = array([100, 200, 300])
Now automatic broadcasting will fail, because automatic reshaping does not work:
M + v # shape mismatch error
The solution is thus to take care of the reshaping manually. What we want in that case is to
add 1 on the right, that is, transform the vector into a column matrix. The broadcasting then
works directly:
M + v.reshape(-1,1)
[ 109 ]
Advanced Array Concepts
For the shape parameter -1, refer to section Accessing and changing the shape of Chapter 4,
Linear Algebra - Arrays. The result is this matrix:
Typical examples
Let us examine some typical examples where broadcasting may come in handy.
Rescale rows
Suppose M is an n × m matrix, and we want to multiply each row by a coefficient. The
coefficients are stored in a vector coeff with n components. In that case, automatic
reshaping will not work, and we have to execute:
rescaled = M*coeff.reshape(-1,1)
Rescale columns
The setup is the same here, but we would like to rescale each column with a coefficient
stored in a vector coeff of length m. In this case, automatic reshaping will work:
rescaled = M*coeff
Obviously, we may also do the reshaping manually and achieve the same result with:
rescaled = M*coeff.reshape(1,-1)
Functions of two variables
Suppose u and v are vectors and we want to form the matrix W with elements wij = ui + vj.
This would correspond to the function F(x, y) = x + y. The matrix W is merely defined by:
W=u.reshape(-1,1) + v
[ 110 ]
Advanced Array Concepts
If the vectors u and v are [0, 1] and [0, 1, 2] respectively, the result is:
More generally, suppose that we want to sample the function w(x, y) := cos(x) + sin(2y).
Supposing that the vectors x and y are defined, the matrix w of sampled values is obtained
by:
w = cos(x).reshape(-1,1) + sin(2*y)
Note that this is very frequently used in combination with ogrid. The vectors obtained
from ogrid are already conveniently shaped for broadcasting. This allows for the following
elegant sampling of the function cos(x) + sin(2y):
x,y = ogrid[0:1:3j,0:1:3j]
# x,y are vectors with the contents of linspace(0,1,3)
w = cos(x) + sin(2*y)
The syntax of ogrid needs some explanation. First, ogrid is no function. It is an instance of
a class with a __getitem__ method (refer to section Attributes in Chapter 8, Classes). That
is why it is used with brackets instead of parentheses.
The two commands are equivalent:
x,y = ogrid[0:1:3j, 0:1:3j]
x,y = ogrid.__getitem__((slice(0, 1, 3j),slice(0, 1, 3j)))
The stride parameter in the preceding example is a complex number. This is to indicate that
it is the number of steps instead of the step size. The rules for the stride parameter might be
confusing at first glance:
If the stride is a real number, then it defines the size of the steps between start
and stop and stop is not included in the list.
If the stride is a complex number s, then the integer part of s.imag defines the
number of steps between start and stop and stop is included in the list.
[ 111 ]
Advanced Array Concepts
Another example for the output of ogrid is a tuple with two arrays, which can be used for
broadcasting:
x,y = ogrid[0:1:3j, 0:1:3j]
gives:
array([[
[
[
array([[
0. ],
0.5],
1. ]])
0. , 0.5,
1. ]])
which is equivalent to:
x,y = ogrid[0:1.5:.5, 0:1.5:.5]
Sparse matrices
Matrices with a small number of nonzero entries are called sparse matrices. Sparse matrices
occur, for example, in scientific computing when describing discrete differential operators
in the context of numerically solving partial differential equations.
Sparse matrices often have large dimensions, sometimes so large that the entire matrix
(with zero entries) would not even fit in the available memory. This is one motivation for a
special type for sparse matrices. Another motivation is better performance of operations
where zero matrix entries can be avoided.
There are only a very limited number of algorithms for general, unstructured sparse
matrices in linear algebra. Most of them are iterative in nature and based on efficient
implementations of matrix-vector multiplication for sparse matrices.
Examples for sparse matrices are diagonal or banded matrices. The simple pattern of these
matrices allows straightforward storing strategies; the principal diagonal and the sub- and
super-diagonals are stored in 1D arrays. Conversion from a sparse representation to the
classical array type and vice-versa can be done by the command diag.
[ 112 ]
Advanced Array Concepts
In general, there is not such a simple structure and the description of sparse matrices
requires special techniques and standards. Here we present a row and a column oriented
type for sparse matrices, both available through the module scipy.sparse .
Figure 5.2: A stiffness matrix from a finite element model of an elastic plate. The pixels denote nonzero entries
in the 1250 × 1250 matrix
Sparse matrix formats
The scipy.sparse module provides many different storing formats from sparse matrices.
We describe here only the most important ones: CSR, CSC, and LIL. The LIL format should
be used for generating and altering sparse matrices; CSR and CSC are efficient formats for
matrix-matrix and matrix-vector operations.
Compressed sparse row
The compressed sparse row format (CSR) uses three arrays: data, indptr, and indices:
The 1D array data stores all the nonzero values in order. It has as many elements
as there are nonzero elements, often denoted by the variable nnz.
[ 113 ]
Advanced Array Concepts
The 1D array indptr contains integers such that indptr[i] is the index of the
element in data, which is the first nonzero element of row i. If the entire row i is
zero, then indptr[i]==indptr[i+1]. If the original matrix has m rows, then
len(indptr)==m+1.
The 1D array indices contains the column index information in such a way that
indices[indptr[i]:indptr[i+1]] is an integer array with the column
indexes of the nonzero elements in row i. Obviously,
len(indices)==len(data)==nnz.
Let's see an example:
The CSR format of the matrix:
is given by the three arrays:
data = (1. 2. 3. 4.)
indptr = (0 2 2 3 5)
indices = (0 2 0 0 3)
The module scipy.sparse provides a type, csr_matrix, with a constructor, which can be
used in the following ways:
With a 2D array as argument
With a matrix in one of the other sparse formats in scipy.sparse
With a shape argument, (m,n), to generate a zero matrix in CSR format
By a 1D array for the data and an integer array ij with the shape
(2,len(data)) such that ij[0,k] is the row index and ij[1,k] is the column
index of data[k] of the matrix
The three arguments, data, indptr, and indices, can be given to the
constructor directly
The first two options are there for conversion purposes while the last two directly define the
sparse matrix.
[ 114 ]
Advanced Array Concepts
Consider the above example in python look like:
import scipy.sparse as sp
A = array([[1,0,2,0],[0,0,0,0],[3.,0.,0.,0.],[1.,0.,0.,4.]])
AS = sp.csr_matrix(A)
Among others, the following attributes are provided:
AS.data
AS.indptr
AS.indices
AS.nnz
#
#
#
#
returns
returns
returns
returns
array([ 1., 2., 3., 1., 4.])
array([0, 2, 2, 3, 5])
array([0, 2, 0, 0, 3])
5
Compressed Sparse Column
The CSR format has a column oriented twin – the compressed sparse column (CSC) format.
The only difference in it compared to the CSR format is the definition of the indptr and
indices arrays, which are now column-related. The type for the CSC format is
csc_matrix and its use corresponds to csr_matrix, explained previously in this section.
Continuing the same example in CSC format:
import scipy.sparse as sp
A = array([[1,0,2,0],[0,0,0,0],[3.,0.,0.,0.],[1.,0.,0.,4.]])
AS = sp.csc_matrix(A)
AS.data
# returns array([ 1., 3., 1., 2., 4.])
AS.indptr
# returns array([0, 3, 3, 4, 5])
AS.indices
# returns array([0, 2, 3, 0, 3])
AS.nnz
# returns 5
Row-based linked list format
The linked list sparse format stores the nonzero matrix entries rowwise in a list data such
that data[k] is a list of the nonzero entries in row k. If all entries in that row are 0, it
contains an empty list.
A second list, rows, contains at position k a list of column indexes of the nonzero elements
in row k. Here is an example in Row-Based linked List Format (LIL) format:
import scipy.sparse as sp
A = array([[1,0,2,0],[0,0,0,0], [3.,0.,0.,0.], [1.,0.,0.,4.]])
AS = sp.lil_matrix(A)
AS.data
# returns array([[1.0, 2.0], [], [3.0], [1.0, 4.0]],
dtype=object)
[ 115 ]
Advanced Array Concepts
AS.rows
AS.nnz
# returns array([[0, 2], [], [0], [0, 3]], dtype=object)
# returns 5
Altering and slicing matrices in LIL format
The LIL format is the one best suited for slicing, that is, extracting submatrices in LIL
format, and for changing the sparsity pattern by inserting nonzero elements. Slicing is
demonstrated by the next example:
BS = AS[1:3,0:2]
BS.data
# returns array([[], [3.0]], dtype=object)
BS.rows
# returns array([[], [0]], dtype=object)
Insertion of a new nonzero element automatically updates the attributes:
AS[0,1] = 17
AS.data # returns array([[1.0, 17.0, 2.0], [], [3.0], [1.0, 4.0]])
AS.rows
# returns array([[0, 1, 2], [], [0], [0, 3]])
AS.nnz
# returns 6
These operations are discouraged in the other sparse matrix formats as they are extremely
inefficient.
Generating sparse matrices
The NumPy commands eye, identity, diag, and rand have their sparse counterparts.
They take an additional argument; it specifies the sparse matrix format of the resulting
matrix.
The following commands generate the identity matrix but in different sparse matrix
formats:
import scipy.sparse as sp
sp.eye(20,20,format = 'lil')
sp.spdiags(ones((20,)),0,20,20, format = 'csr')
sp.identity(20,format ='csc')
The sp.rand command takes an additional argument describing the density of the
generated random matrix. A dense matrix has density 1 while a zero matrix has density 0:
import scipy.sparse as sp
AS=sp.rand(20,200,density=0.1,format=’csr’)
AS.nnz # returns 400
[ 116 ]
Advanced Array Concepts
There is no direct correspondence to the NumPy command zeroes. Matrices completely
filled with zeros are generated by instantiating the corresponding type with the shape
parameters as constructor parameters:
import scipy.sparse as sp
Z=sp.csr_matrix((20,200))
Z.nnz
# returns 0
Sparse matrix methods
There are methods to convert one sparse type into another or into an array:
AS.toarray # converts sparse formats to a numpy array
AS.tocsr
AS.tocsc
AS.tolil
The type of a sparse matrix can be inspected by the methods issparse , isspmatrix_lil,
isspmatrix_csr, and isspmatrix_csc.
Elementwise operations +, *, /, and ** on sparse matrices are defined as for NumPy arrays.
Regardless of the sparse matrix format of the operands, the result is always a csr_matrix.
Applying elementwise operating functions to sparse matrices requires first transforming
them to either CSR or CSC format and applying the functions to their data attribute, as
demonstrated by the next example.
The elementwise sine of a sparse matrix can be defined by an operation on its data
attribute:
import scipy.sparse as sp
def sparse_sin(A):
if not (sp.isspmatrix_csr(A) or sp.isspmatrix_csc(A)):
A = A.tocsr()
A.data = sin(A.data)
return A
[ 117 ]
Advanced Array Concepts
For matrix-matrix or matrix-vector multiplications, there is a sparse matrix method, dot. It
returns either a csr_matrix or a 1D NumPy array:
import scipy.sparse as sp
A = array([[1,0,2,0],[0,0,0,0],[3.,0.,0.,0.],[1.,0.,0.,4.]])
AS = sp.csr_matrix(A)
b = array([1,2,3,4])
c = AS.dot(b)
# returns array([ 7., 0., 3., 17.])
C = AS.dot(AS)
# returns csr_matrix
d = dot(AS,b)
# does not return the expected result!
Avoid using NumPy's command dot on sparse matrices, as this might
lead to unexpected results. Use the command dot from scipy.sparse
instead.
Other linear algebra operations such as system solving, least squares, eigenvalues, and
singular values are provided by the scipy.sparse.linalg module.
Summary
The concept of views is one of the important topics you should have learned from this
chapter. Missing this topic will give you a hard time when debugging your code. Boolean
arrays occur at various places throughout this book. They are handy and compact tools for
avoiding lengthy if constructions and loops when working with arrays. In nearly all large
computational projects, sparse matrices become an issue. You saw how these are handled
and which related methods are available.
[ 118 ]
6
Plotting
Plotting in Python can be done with the pyplot part of the matplotlib module. With
matplotlib you can create high-quality figures and graphics and also plot and visualize your
results. Matplotlib is open source and freely available software, [21]. The matplotlib
website also contains excellent documentation with examples, [35]. In this section, we will
show you how to use the most common features. The examples in the upcoming sections
assume that you have imported the module as:
from matplotlib.pyplot import *
In case you want to use the plotting commands in IPython, it is recommended that you run
the magic command %matplotlib directly after starting the IPython shell. This prepares
IPython for interactive plotting.
Basic plotting
The standard plotting function is plot. Calling plot(x,y) creates a figure window with a
plot of y as a function of x. The input arguments are arrays (or lists) of equal length. It is also
possible to use plot(y), in which case the values in y will be plotted against their index,
that is, plot(y) is a short form of plot(range(len(y)),y).
Here is an example that shows how to plot sin(x) for x ϵ [-2π, 2π] using 200 sample points
and sets markers at every fourth point:
# plot sin(x) for some interval
x = linspace(-2*pi,2*pi,200)
plot(x,sin(x))
# plot marker for every 4th point
samples = x[::4]
Plotting
plot(samples,sin(samples),'r*')
# add title and grid lines
title('Function sin(x) and some points plotted')
grid()
The result is shown in the following figure (Figure 6.1):
Figure 6.1: A plot of the function sin(x) with grid lines shown.
[ 120 ]
Plotting
As you can see, the standard plot is a solid blue curve. Each axis gets automatically scaled
to fit the values but can also be set manually. Color and plot options can be given after the
first two input arguments. Here, r* indicates red star-shaped markers. Formatting is
covered in more detail in the next section. The title command puts a title text string above
the plot area.
Calling plot multiple times will overlay the plots in the same window. To get a new clean
figure window, use figure(). The figure command might contain an integer, for
example, figure(2), which can be used to switch between figure windows. If there is no
figure window with that number, a new one is created, otherwise, the window is activated
for plotting and all subsequent plotting commands apply to that window.
Multiple plots can be explained using the legend function, along with adding labels to each
plot call. The following example fits polynomials to a set of points using the commands
polyfit and polyval, and plots the result with a legend:
# —Polyfit example—
x = range(5)
y = [1,2,1,3,5]
p2 = polyfit(x,y,2)
p4 = polyfit(x,y,4)
# plot the polynomials and points
xx = linspace(-1,5,200)
plot(xx, polyval(p2, xx), label='fitting polynomial of degree 2')
plot(xx, polyval(p4, xx),
label='interpolating polynomial of degree 4')
plot(x,y,'*')
# set the axis and legend
axis([-1,5,0,6])
legend(loc='upper left', fontsize='small')
[ 121 ]
Plotting
Here you can also see how to manually set the range of the axis using
axis([xmin,xmax,ymin,ymax]). The legend command takes optional arguments on
placement and formatting; in this case the legend is put in the upper-left corner and typeset
with a small font size, as shown in the following figure (Figure 6.2).
Figure 6.2: Two polynomials fitted to the same points.
As final examples for basic plotting, we demonstrate how to do scatter plots and
logarithmic plots in two dimensions.
[ 122 ]
Plotting
Example of 2D point scatter plot:
# create random 2D points
import numpy
x1 = 2*numpy.random.standard_normal((2,100))
x2 = 0.8*numpy.random.standard_normal((2,100)) + array([[6],[2]])
plot(x1[0],x1[1],'*')
plot(x2[0],x2[1],'r*')
title('2D scatter plot')
Figure 6.3(a): An example of a scatter plot
[ 123 ]
Plotting
The following code is an example of a logarithmic plot using loglog:
# log both x and y axis
x = linspace(0,10,200)
loglog(x,2*x**2, label = 'quadratic polynomial',
linestyle = '-', linewidth = 3)
loglog(x,4*x**4, label = '4th degree polynomial',
linestyle = '-.', linewidth = 3)
loglog(x,5*exp(x), label = 'exponential function', linewidth = 3)
title('Logarithmic plots')
legend(loc = 'best')
Figure 6.3(b): An example of a plot with logarithmic x and y axis
The examples shown in the preceding figure (Figure 6.3(a) and Figure 6.3(b)) used some
parameters of plot and loglog which allow special formatting. In the next section, we will
explain the parameters in more detail.
[ 124 ]
Plotting
Formatting
The appearance of figures and plots can be styled and customized to look how you want
them to look. Some important variables are linewidth, which controls the thickness of plot
lines; xlabel, ylabel, which set the axis labels, color for plot colors, and transparent
for transparency. This section will tell you how to use some of them. The following is an
example with more keywords:
k = 0.2
x = [sin(2*n*k) for n in range(20)]
plot(x, color='green', linestyle='dashed', marker='o',
markerfacecolor='blue', markersize=12, linewidth=6)
There are short commands that can be used if you only need basic style changes, for
example, setting the color and line style. The following table (Table 6.1) shows some
examples of these formatting commands. You may use either the short string syntax
plot(...,'ro-'), or the more explicit syntax plot(..., marker='o', color='r',
linestyle='-').
Table 6.1: Some common plot formatting arguments
[ 125 ]
Plotting
To set the color to green with the 'o' marker we write:
plot(x,'go')
To plot histograms instead of regular plots, the hist command is used:
# random vector with normal distribution
sigma, mu = 2, 10
x = sigma*numpy.random.standard_normal(10000)+mu
hist(x,50,normed=1)
z = linspace(0,20,200)
plot(z, (1/sqrt(2*pi*sigma**2))*exp(-(z-mu)**2/(2*sigma**2)),'g')
# title with LaTeX formatting
title('Histogram with '.format(mu,sigma))
Figure 6.4 normal distribution with 50 bins and a green curve indicating the true
distribution
The resulting plot looks similar to the preceding figure (Figure 6.4). The title, and any other
text, can be formatted using LaTeX to show mathematical formulas. LaTeX formatting is
enclosed within a pair of $ signs. Also, note the string formatting done using the format
method, refer to section Strings in Chapter 2, Variables and Basic Types.
[ 126 ]
Plotting
Sometimes the brackets for the string formatting interfere with LaTeX bracket
environments. If this occurs, replace the LaTeX bracket with a double bracket, for example,
x_{1} should be replaced with x_{{1}}. The text might contain sequences that overlap
with string escape sequences, for example, \tau will be interpreted as the tab character \t.
An easy workaround is to add r before the string, for example r'\tau'; this makes it a
raw string.
Placing several plots in one figure window can be done using the subplot command.
Consider the following example, which iteratively averages out the noise on a sine curve.
def avg(x):
""" simple running average """
return (roll(x,1) + x + roll(x,-1)) / 3
# sine function with noise
x = linspace(-2*pi, 2*pi,200)
y = sin(x) + 0.4*rand(200)
# make successive subplots
for iteration in range(3):
subplot(3, 1, iteration + 1)
plot(x,y, label = '{:d} average{}'.format(iteration, 's' if iteration >
1 else ''))
yticks([])
legend(loc = 'lower left', frameon = False)
y = avg(y) #apply running average
subplots_adjust(hspace = 0.7)
Figure 6.5: An example of plotting several times in the same figure window.
[ 127 ]
Plotting
The function avg uses a roll call to shift all values of the array. subplot takes three
arguments: the number of vertical plots, the number of horizontal plots, and an index
indicating which location to plot in (counted row-wise). Note that we used the
subplots_adjust command to add extra space to adjust the distance between both the
subplots.
A useful command is savefig which lets you save a figure as an image (this can also be
done from the figure window). Many image and file formats are supported by this
command, they are specified by the filename's extension as:
savefig('test.pdf')
# save to pdf
savefig('test.svg')
# save to svg (editable format)
or
You can place the image against a non-white background, for example, a webpage. For this,
the transparent parameter can be set to make the figure's background transparent:
savefig('test.pdf', transparent=True)
If you intend to embed a figure into a LaTeX document, it is recommended that you reduce
the surrounding white space by setting the figure's bounding box tight around the drawing,
as shown here:
savefig('test.pdf', bbox_inches='tight')
Meshgrid and contours
A common task is a graphical representation of a scalar function over a rectangle:
For this, first we have to generate a grid on the rectangle [a,b] x [c,d]. This is done using the
meshgrid command:
n = ... # number of discretization points along the x-axis
m = ... # number of discretization points along the x-axis
X,Y = meshgrid(linspace(a,b,n), linspace(c,d,m))
[ 128 ]
Plotting
X and Y are arrays with (n,m) shape such that
contains the coordinates of the grid point
as
shown in the next figure (Figure 6.6):
Figure 6.6: A rectangle discretized by meshgrid
A rectangle discretized by meshgrid will be used to visualize the behavior of an iteration.
Bur first we will use it to plot level curves of a function. This is done by the command
contour.
As an example we choose Rosenbrock's banana function:
It is used to challenge optimization methods. The function values descend towards a
banana-shaped valley, which itself decreases slowly towards the function’s global minimum
at (1, 1).
First we display the level curves using contour.
rosenbrockfunction = lambda x,y: (1-x)**2+100*(y-x**2)**2
X,Y = meshgrid(linspace(-.5,2.,100), linspace(-1.5,4.,100))
Z = rosenbrockfunction(X,Y)
contour(X,Y,Z,logspace(-0.5,3.5,20,base=10),cmap='gray')
title('Rosenbrock Function: ')
xlabel('x')
ylabel('y')
[ 129 ]
Plotting
This plots the level curve at the levels given by the fourth parameter and uses the colormap
0.5
3
gray. Furthermore, we used logarithmically spaced steps from 10 to 10 using the function
logscale to define the levels, as shown in the next figure.
Figure 6.7: A contour plot of Rosenbrock function
In the preceding example, an anonymous function indicated by the keyword lambda is
used to keep the code compact. Anonymous functions are explained in section Anonymous
functions – the lambda keyword in Chapter 7, Functions, Anonymous functions. If levels are not
given as arguments to contour, the function chooses appropriate levels by itself .
The contourf function performs the same function as contour but fills the plot with colors
according to different levels. Contour plots are ideal for visualizing the behavior of a
numerical method. We illustrate this here by showing the iterations of an optimization
method.
[ 130 ]
Plotting
We continue the preceding example and depict the steps towards the minimum of the
Rosenbrock function generated by Powell's method, [27], which we will apply to find the
minimum of the Rosenbrock function:
import scipy.optimize as so
rosenbrockfunction = lambda x,y: (1-x)**2+100*(y-x**2)**2
X,Y=meshgrid(linspace(-.5,2.,100),linspace(-1.5,4.,100))
Z=rosenbrockfunction(X,Y)
cs=contour(X,Y,Z,logspace(0,3.5,7,base=10),cmap='gray')
rosen=lambda x: rosenbrockfunction(x[0],x[1])
solution, iterates = so.fmin_powell(rosen,x0=array([0,-0.7]),retall=True)
x,y=zip(*iterates)
plot(x,y,'ko') # plot black bullets
plot(x,y,'k:',linewidth=1) # plot black dotted lines
title("Steps of Powell's method to compute a minimum")
clabel(cs)
The iterative method fmin_powell applies Powell's method to find a minimum. It is
started by a given start value of x0 and reports all iterates when the option retall=True is
given. After sixteen iterations, the solution x=0, y=0 was found. The iterations are depicted
as bullets in the following contour plot (Figure 6.8).
Figure 6.8: A contour plot of Rosenbrock function with a search path of an optimization
method
[ 131 ]
Plotting
contour also creates a contour set object that we assigned to the variable cs. This is then
used by clabel to annotate the levels of the corresponding function values, as shown in the
preceding figure (Figure 6.8).
Images and contours
Let us take a look at some examples of visualizing arrays as images. The following function
will create a matrix of color values for the Mandelbrot fractal. Here we consider a fixed
point iteration, that depends on a complex parameter c:
Depending on the choice of this parameter it may or may not create a bounded sequence of
complex values zn.
For every value of c, we check if zn exceeds a prescribed bound. If it remains below the
bound within maxit iterations, we assume the sequence to be bounded.
Note how, in the following piece of code,meshgrid is used to generate a matrix of complex
parameter values c:
def mandelbrot(h,w, maxit=20):
X,Y = meshgrid(linspace(-2, 0.8, w), linspace(-1.4, 1.4, h))
c = X + Y*1j
z = c
exceeds = zeros(z.shape, dtype=bool)
for iteration in range(maxit):
z = z**2 + c
exceeded = abs(z) > 4
exceeds_now = exceeded & (logical_not(exceeds))
exceeds[exceeds_now] = True
z[exceeded] = 2 # limit the values to avoid overflow
return exceeds
imshow(mandelbrot(400,400),cmap='gray')
axis('off')
[ 132 ]
Plotting
The command imshow displays the matrix as an image. The selected color map shows the
regions where the sequence appeared unbounded in white and others in black. Here we
used axis('off') to turn off the axis as this might be not so useful for images.
Figure 6.9: An example of using imshow to visualize a matrix as an image.
By default, imshow uses interpolation to make the images look nicer. This is clearly seen
when the matrices are small. The next figure shows the difference between using:
imshow(mandelbrot(40,40),cmap='gray')
and
imshow(mandelbrot(40,40), interpolation='nearest', cmap='gray')
[ 133 ]
Plotting
In the second example, pixel values are just replicated.
Figure 6.10: The difference between using the linear interpolation of imshow
compared to using nearest neighbor interpolation
For more details on working and plotting with images using Python refer to [30].
Matplotlib objects
Till now, we have used the pyplot module of matplotlib. This module makes it easy for us
to use the most important plot commands directly. Mostly, we are interested in creating a
figure and display it immediately. Sometimes, though, we want to generate a figure that
should be modified later by changing some of its attributes. This requires us to work with
graphical objects in an object-oriented way. In this section, we will present some basic steps
to modify figures. For a more sophisticated object oriented approach to plotting in Python,
you have to leave pyplot and have to dive directly into matplotlib with its extensive
documentation.
[ 134 ]
Plotting
The axes object
When creating a plot that should be modified later, we need references to a figure and an
axes object. For this we have to create a figure first and then define some axes and their
location in the figure. And we should not forget to assign these objects to a variable:
fig = figure()
ax = subplot(111)
A figure can have several axes objects depending on the use of subplot. In a second step
plots are associated with a given axes object:
fig = figure(1)
ax = subplot(111)
x = linspace(0,2*pi,100)
# We set up a function that modulates the amplitude of the sin function
amod_sin = lambda x: (1.-0.1*sin(25*x))*sin(x)
# and plot both...
ax.plot(x,sin(x),label = 'sin')
ax.plot(x, amod_sin(x), label = 'modsin')
Here we used an anonymous function indicated by the lambda keyword . We will explain
this construct later in section Anonymous functions – the lambda keyword in Chapter 7,
Functions. In fact, these two plot commands fill the list ax.lines with two Lines2D objects:
ax.lines #[, ]
It is a good practice to use labels so that we can later identify objects in an easy way:
for il,line in enumerate(ax.lines):
if line.get_label() == 'sin':
break
We set up now things in a way that allows further modifications. The figure we got so far is
shown in preceding figure (Figure 6.11, left).
[ 135 ]
Plotting
Modifying line properties
We just identified a particular line object by its label. It is an element of the list ax.lines
list with the index il . All its properties are collected in a dictionary
dict_keys(['marker', 'markeredgewidth', 'data', 'clip_box',
'solid_capstyle', 'clip_on', 'rasterized', 'dash_capstyle', 'path',
'ydata', 'markeredgecolor', 'xdata', 'label', 'alpha', 'linestyle',
'antialiased', 'snap', 'transform', 'url',
'transformed_clip_path_and_affine', 'clip_path', 'path_effects',
'animated', 'contains', 'fillstyle', 'sketch_params', 'xydata',
'drawstyle', 'markersize', 'linewidth', 'figure', 'markerfacecolor',
'pickradius', 'agg_filter', 'dash_joinstyle', 'color', 'solid_joinstyle',
'picker', 'markevery', 'axes', 'children', 'gid', 'zorder', 'visible',
'markerfacecoloralt'])
which can be obtained by the command:
ax.lines[il].properties()
They can be changed by corresponding setter methods. Let us change the line style of the
sine – curve:
ax.lines[il].set_linestyle('-.')
ax.lines[il].set_linewidth(2)
We can even modify the data, as shown:
ydata=ax.lines[il].get_ydata()
ydata[-1]=-0.5
ax.lines[il].set_ydata(ydata)
[ 136 ]
Plotting
The result is shown in the next figure(Figure 6.11, right):
Figure 6.11: The amplitude modulated sine-function (left) and a curve with the
last data point corrupted (right).
Annotations
One useful axes method is annotate. It sets an annotation at a given position and points,
with an arrow, to another position in the drawing. The arrow can be given properties in a
dictionary:
annot1=ax.annotate('amplitude modulated\n curve', (2.1,1.0),(3.2,0.5),
arrowprops={'width':2,'color':'k',
'connectionstyle':'arc3,rad=+0.5',
'shrink':0.05},
verticalalignment='bottom', horizontalalignment='left',fontsize=15,
bbox={'facecolor':'gray', 'alpha':0.1, 'pad':10})
annot2=ax.annotate('corrupted data', (6.3,-0.5),(6.1,-1.1),
arrowprops={'width':0.5,'color':'k','shrink':0.1},
horizontalalignment='center', fontsize=12)
In the first annotation example above, the arrow points to a point with the coordinates (2.1,
1.0) and the left bottom coordinate of the text is (3.2, 0.5). If not otherwise specified, the
coordinates are given in the convenient data-coordinate system, which refers to the data
used to generate the plots.
[ 137 ]
Plotting
Furthermore, we demonstrated a couple of arrow properties specified by the arrowprop
dictionary. You can scale the arrow by the shrink key. The setting 'shrink':0.05
reduces the arrow size by 5% to keep a distance to the curve it points to. You can let the
arrow follow a spline arc or give it other shapes using the connectionstyle key.
Text properties or even a bounding box around the text can be made by extra keyword
arguments to the annotate method, refer to the following figure (Figure 6.12, left):
Experimenting with annotations requires sometimes to remove attempts that we would like
to reject. Therefore we assigned the annotate object to a variable, which allows us to remove
the annotation by its remove method:
annot1.remove()
Filling areas between curves
Filling is an ideal tool to highlight differences between curves, such as noise on top of
expected data, approximations versus exact functions, and so on.
Filling is done by the axis method
ax.fill_between(x,y1,y2)
For the next figure we used:
axf = ax.fill_between(x, sin(x), amod_sin(x), facecolor='gray')
where is a very convenient parameter that needs a Boolean array to specify the additional
filling conditions.
axf = ax.fill_between(x, sin(x), amod_sin(x),where=amod_sin(x)-sin(x) > 0,
facecolor=’gray’)
The Boolean array which selects the regions to fill is amod_sin(x)-sin(x) > 0.
[ 138 ]
Plotting
The next figure shows the curve with both variants of filling areas:
Figure 6.12: The amplitude modulated sin-function with annotations and filled
areas(left) and a modified figure with only partially filled areas by using the
where parameter (right).
If you test these commands yourself, do not forget to remove the complete filling before you
try out the partial filling, otherwise you will not see any change:
axf.remove()
Related filling commands are fill and fill_betweenx.
[ 139 ]
Plotting
Ticks and ticklabels
Figures in talks, posters, and publications look much nicer if they are not overloaded with
unnecessary information. You want to direct the spectator to those parts that contain the
message. In our example, we clean up the picture by removing ticks from the x-axis and yaxis and by introducing problem related tick labels:
Figure 6.13: The completed example of the amplitude modulated sine – function with
annotations and filled areas and modified ticks and tick labels.
ax.set_xticks(array([0,pi/2,pi,3/2*pi,2*pi]))
ax.set_xticklabels(('$0$','$\pi/2$','$\pi$','$3/2 \pi$','$2
\pi$'),fontsize=18)
ax.set_yticks(array([-1.,0.,1]))
ax.set_yticklabels(('$-1$','$0$','$1$'),fontsize=18)
[ 140 ]
Plotting
Note that we used LaTeX formatting in the strings to represent Greek letters, to set formulas
correctly, and to use a LaTeX font. It is also a good practice to increase the font size so that
the resulting figure can be scaled down into a text document without affecting the
readability of the axes. The final result of this guiding example is shown in the previous
figure (Figure 6.13).
Making 3D plots
There are some useful matplotlib tool kits and modules that can be used for a variety of
special purposes. In this section, we describe a method for producing 3D-plots.
The mplot3d toolkit provides 3D plotting of points, lines, contours, surfaces, and all other
basic components as well as 3D rotation and scaling. Making a 3D plot is done by adding
the keyword projection='3d' to the axes object as shown in the following example:
from mpl_toolkits.mplot3d import axes3d
fig = figure()
ax = fig.gca(projection='3d')
# plot points in 3D
class1 = 0.6 * random.standard_normal((200,3))
ax.plot(class1[:,0],class1[:,1],class1[:,2],'o')
class2 = 1.2 * random.standard_normal((200,3)) + array([5,4,0])
ax.plot(class2[:,0],class2[:,1],class2[:,2],'o')
class3 = 0.3 * random.standard_normal((200,3)) + array([0,3,2])
ax.plot(class3[:,0],class3[:,1],class3[:,2],'o')
[ 141 ]
Plotting
As you can see, you need to import the axes3D type from mplot3d. The resulting plot
displays the scattered 3D-data which can be seen in the following figure (Figure 6.14)
Figure 6.14: Plotting 3D data using mplot3d toolkit
Plotting surfaces is just as easy. The following example uses the built-in function
get_test_data to create a sample data for plotting a surface. Consider the following
example of a surface plot with transparency.
X,Y,Z = axes3d.get_test_data(0.05)
fig = figure()
ax = fig.gca(projection='3d')
# surface plot with transparency 0.5
ax.plot_surface(X,Y,Z,alpha=0.5)
[ 142 ]
Plotting
The alpha value sets the transparency. The surface plot is shown in the following figure
(Figure 6.15).
Figure 6.15: Example for plotting a surface mesh with three 2D projections.
You can also plot contours in any of the coordinate projections as in the next example.
fig = figure()
ax = fig.gca(projection = '3d')
ax.plot_wireframe(X,Y,Z,rstride = 5,cstride = 5)
# plot contour projection on each
ax.contour(X,Y,Z, zdir='z',offset
ax.contour(X,Y,Z, zdir='x',offset
ax.contour(X,Y,Z, zdir='y',offset
axis plane
= -100)
= -40)
= 40)
# set axis limits
ax.set_xlim3d(-40,40)
ax.set_ylim3d(-40,40)
ax.set_zlim3d(-100,100)
# set labels
ax.set_xlabel('X axis')
ax.set_ylabel('Y axis')
ax.set_zlabel('Z axis')
[ 143 ]
Plotting
Note the commands for setting the axis limits. The standard matplotlib commands for
setting the axis limits are axis([-40, 40, -40, 40]), this works fine for 2D plots.
However, axis([-40,40,-40,40,-40,40]) does not work. For 3D plots you need to use
the object oriented version of the commands, ax.set_xlim3d(-40,40) and likewise. The
same goes for labeling the axis; note the commands for setting the labels. For 2D plots you
can do xlabel(’X axis’) and ylabel(’Y axis’) but there is no zlabel command.
Instead, in 3D plots you need to use ax.set_xlabel(’X axis’) and likewise, as shown in
the preceding example.
The resulting figure from this code is the following
There are many options for formatting the appearance of the plots, including color and
transparency of surfaces. The mplot3d documentation website, [23], has the details.
[ 144 ]
Plotting
Making movies from plots
If you have data that evolves, you might want to save it as a movie besides showing it in a
figure window, similar to the savefig command. One way to do this is with the visvis
module available at visvis (refer to [37] for more information).
Here is a simple example of evolving a circle using an implicit representation. Let the circle
be represented by the zero level,
, of a function
. Alternatively, consider
the disk
inside the zero set. If the value of f decreases at a rate v then the circle
will move outward with rate
.
This can be implemented as:
import visvis.vvmovie as vv
# create initial function values
x = linspace(-255,255,511)
X,Y = meshgrid(x,x)
f = sqrt(X*X+Y*Y) - 40 #radius 40
# evolve and store in a list
imlist = []
for iteration in range(200):
imlist.append((f>0)*255)
f -= 1 # move outwards one pixel
vv.images2swf.writeSwf('circle_evolution.swf',imlist)
The result is a Flash movie (*.swf file) of a growing black circle, as shown in the next figure
(Figure 6.16):
Figure 6.16: An example of evolving a circle
[ 145 ]
Plotting
In this example, a list of arrays was used to create the movie. The visvis module can also
save an GIF animation and on certain platforms an AVI animation (*.gif and *.avi files), and
there is also the possibility to capturing movie frames directly from the figure window.
These options, however, require some more packages to be installed on your system (for
example, PyOpenGL and PIL, the Python Imaging Library). See the documentation on the
visvis webpage for more details.
Another option is to use savefig to create images, one for each frame.
# create initial function values
x = linspace(-255,255,511)
X,Y = meshgrid(x,x)
f = sqrt(X*X+Y*Y) - 40 #radius 40
for iteration in range(200):
imshow((f>0)*255)
gray()
axis('off')
savefig('circle_evolution_{:d}.png'.format(iteration))
f -= 1
These images can then be combined using a standard video editing software, for example,
Mencoder or ImageMagick. This approach has the advantage that you can make highresolution videos by saving high-resolution images.
Summary
A graphical representation is the most compact form in which to present mathematical
results or the behavior of an algorithm. This chapter provided you with the basic tools for
plotting and introduced you to a more sophisticated way to work with graphical objects,
such as figures, axes, and lines in an object-oriented way.
In this chapter, you learned how to make plots, not only classical x/y-plots but also 3D-plots
and histograms. We gave you an appetizer on making films. You also saw how to modify
plots considering them to be graphical objects with related methods and attributes which
can be set, deleted, or modified.
[ 146 ]
Plotting
Exercises
Ex. 1 → Write a function that plots an ellipse given its center coordinates (x,y), the half axis a
and b rotation angle θ.
Ex. 2 → Write a short program that takes a 2D array, e.g., the preceding Mandelbrot contour
image, and iteratively replace each value by the average of its neighbors. Update a contour
plot of the array in a figure window to animate the evolution of the contours. Explain the
behavior.
Ex. 3 → Consider an N × N matrix or image with integer values. The mapping
is an example of a mapping of a toroidal square grid of points onto itself. This has the
interesting property that it distorts the image by shearing and then moving the pieces
outside the image back using the modulu function mod. Applied iteratively, this results in
randomizing the image in a way that eventually returns the original. Implement the
following sequence:
and save out the first N steps to files or plot them in a figure window.
As an example image, you can use the classic 512 × 512 Lena test image from scipy.misc.
from scipy.misc import lena
I = lena()
The result should look like this:
…
0
1
…
128
…
256
[ 147 ]
511
512
Plotting
Compute the x and y mappings and use array indexing (refer to section
Array Indexing in Chapter 5, Advance Array Concepts) to copy the pixel
values.
Ex. 4 → Reading and plotting on images. SciPy comes with the imread function (in the
scipy.misc module) for reading images, (refer to section Reading and Writing Images in
Chapter 12, Input and output). Write a short program that reads an image from file and plots
the image contour at a given gray level value overlaid on the original image.
You can get a gray level version of the image by averaging the color
channels like this: mean(im,axis=2)
Ex. 5 → Image edges. The zero crossings of the 2D Laplacian are a good indication of image
edges. Modify the program in the previous exercise to use the gaussian_laplace or
laplace function from the scipy.ndimage module to compute the 2D Laplacian and
overlay the edges on top of the image.
Ex. 6 → Reformulate the Mandelbrod fractal example (see section Images and Contours) by
using orgid instead of meshgrid, see also the explanation ogrid in Function of two variables
in Chapter 5, Advanced Array Concepts. What is the difference between orgid, mgrid, and
meshgrid?
[ 148 ]
7
Functions
This chapter introduces functions, a fundamental building block in programming. We show
how to define them, how to handle input and output, how to properly use them, and how
to treat them as objects.
Basics
In mathematics, a function is written as a map that uniquely assigns an element y from the
range R to every element x from the domain D.
This is expressed by f : D → R
Alternatively, when considering particular elements x and y, one writes f : x → y
Here, f is called the name of the function and f(x) is its value when applied to x. Here, x is
sometimes called the argument of f. Let's first look at an example before considering
functions in Python.
For example, D = ℝ x ℝ and y = f(x1, x2) = x1– x2. This function maps two real numbers to
their difference.
In mathematics, functions can have numbers, vectors, matrices, and even other functions as
arguments. Here is an example of a function with mixed arguments:
.
Functions
In this case, a number is returned. When working with functions, we have to distinguish
between two different steps:
The definition of the function
The evaluation of the function, that is, the computation of f(x) for a given value of
x
The first step is done once, while the second can be performed many times for various
arguments. Functions in programming languages follow the same concept and apply it to a
wide range of types of input arguments, for example, strings, lists, or any object. We
demonstrate a definition of the function by considering the given example again:
def subtract(x1, x2):
return x1 - x2
The keyword def indicates that we are going to define a function. subtract is the
function’s name and x1 and x2 are its parameters. The colon indicates that we are using a
block command and the value that is returned by the function follows the return keyword.
Now, we can evaluate this function. This function is called with its parameters replaced by
input arguments:
r = subtract(5.0, 4.3)
The result 0.7 is computed and assigned to the r variable.
Parameters and arguments
When defining a function, its input variables are called the parameters of the function. The
input used when executing the function is called its argument.
Passing arguments – by position and by keyword
We will consider the previous example again, where the function takes two parameters,
namely x1 and x2.
[ 150 ]
Functions
Their names serve to distinguish the two numbers, which in this case cannot be
interchanged without altering the result. The first parameter defines the number from
which the second parameter is subtracted. When subtract is called, every parameter is
replaced by an argument. Only the order of the arguments matters; the arguments can be
any objects. For instance, we may call the following:
z = 3
e = subtract(5,z)
Besides this standard way of calling a function, which is by passing the arguments by
position, it might sometimes be convenient to pass arguments using keywords. The names
of the parameters are the keywords; consider the following instance:
z = 3
e = subtract(x2 = z, x1 = 5)
Here, the arguments are assigned to the parameters by name and not by position in the call.
Both ways of calling a function can be combined so that the arguments given by position
come first and the arguments given by keyword follow last. We show this by using
the function plot, which was described in Chapter 6, Plotting:
plot(xp, yp, linewidth = 2,label = 'y-values')
Changing arguments
The purpose of parameters is to provide the function with the necessary input data.
Changing the value of the parameter inside the function normally has no effect on its value
outside the function:
def subtract(x1, x2):
z = x1 - x2
x2 = 50.
return z
a = 20.
b = subtract(10, a)
# returns -10
a
# still has the value 20
This applies to all immutable arguments, such as strings, numbers, and tuples. The situation
is different if mutable arguments, such as lists or dictionaries, are changed.
[ 151 ]
Functions
For example, passing mutable input arguments to a function and changing them inside the
function can change them outside the function too:
def subtract(x):
z = x[0] - x[1]
x[1] = 50.
return z
a = [10,20]
b = subtract(a)
# returns -10
a
# is now [10, 50.0]
Such a function misuses its arguments to return results. We strongly dissuade you from
such constructions and recommend that you do not change input arguments inside the
function (for more information refer to Default Arguments section).
Access to variables defined outside the local
namespace
Python allows functions to access variables defined in any of its enclosing program units.
These are called global variables, in contrast to local variables. The latter are only accessible
within the function. For example, consider the following code:
import numpy as np # here the variable np is defined
def sqrt(x):
return np.sqrt(x) # we use np inside the function
This feature should not be abused. The following code is an example of what not to do:
a = 3
def multiply(x):
return a * x # bad style: access to the variable a defined outside
When changing the variable a the function, multiply tacitly changes its behavior:
a=3
multiply(4)
a=4
multiply(4)
# returns 12
# returns 16
[ 152 ]
Functions
It is much better in that case to provide the variable as a parameter through the argument
list:
def multiply(x, a):
return a * x
Global variables can be useful when working with closures. Namespaces and scopes are
discussed more extensively in Chapter 11, Namespaces, Scopes, and Modules.
Default arguments
Some functions can have many parameters, and among them some might only be of interest
in nonstandard situations. It would be practical if arguments could automatically be set to
standard (default) values. We demonstrate the use of default arguments by looking at the
command norm in the scipy.linalg module. It computes various norms of matrices and
vectors.
The following snippet calls for computing the Frobenius norm of the 3 × 3 identity matrix
are equivalent (more on matrix norms can be found in [10]):
import scipy.linalg as sl
sl.norm(identity(3))
sl.norm(identity(3), ord = 'fro')
sl.norm(identity(3), 'fro')
Note that in the first call, no information about the ord keyword is given. How does Python
know that it should compute the Frobenius norm and not another norm, for example, the
Euclidean 2-norm?
The answer to the previous question is the use of default values. A default value is a value
already given by the function definition. If the function is called without providing this
argument, Python uses the value that the programmer provided when the function was
defined.
Suppose we call the function subtract only one argument; we get an error message:
TypeError: subtract() takes exactly 2 arguments (1 given)
To allow the omission of the argument x2, the definition of the function has to provide a
default value, for example:
def subtract(x1, x2 = 0):
return x1 - x2
[ 153 ]
Functions
To summarize, arguments can be given as positional arguments and keyword arguments.
All positional arguments have to be given first. You do not need to provide all keyword
arguments as long as those omitted arguments have default values in the function
definition.
Beware of mutable default arguments
The default arguments are set upon function definition. Changing mutable arguments
inside a function has a side effect when working with default values, for example:
def my_list(x1, x2 = []):
x2.append(x1)
return x2
my_list(1) # returns [1]
my_list(2) # returns [1,2]
Variable number of arguments
Lists and dictionaries may be used to define or call functions with a variable number of
arguments. Let's define a list and a dictionary as follows:
data = [[1,2],[3,4]]
style = dict({'linewidth':3,'marker':'o','color':'green'})
Then we can call the plot function using starred (*) arguments:
plot(*data,**style)
A variable name prefixed by * , such as *data in the preceding example, means that a list
that gets unpacked in the function call is provided. In this way, a list generates positional
arguments. Similarly, a variable name prefixed by **, such as **style in the example,
unpacks a dictionary to keyword arguments. Refer to the following figure (Figure 7.1):
[ 154 ]
Functions
Figure 7.1: Starred arguments in function calls
You might also want to use the reverse process, where all given positional arguments are
packed into a list and all keyword arguments are packed into a dictionary when passed to a
function.
In the function definition, this is indicated by parameters prefixed by * and **,
respectively. You will often find the *args and **kwargs parameters in code
documentation, refer Figure 7.2.
Figure 7.2: Starred arguments in function definitions
[ 155 ]
Functions
Return values
A function in Python always returns a single object. If a function has to return more than
one object, these are packed and returned as a single tuple object.
For instance, the following function takes a complex number z and returns its polar
coordinate representation as magnitude r and angle according to Euler’s formula:
And the Python counterpart would be this:
def complex_to_polar(z):
r = sqrt(z.real ** 2 + z.imag ** 2)
phi = arctan2(z.imag, z.real)
return (r,phi) # here the return object is formed
Here, we used the sqrt(x) NumPy function for the square root of a number x and
-1
arctan2(x,y) for the expression tan (x/y).
Let us try our function:
z =
a =
r =
phi
3 + 5j # here we define a complex number
complex_to_polar(z)
a[0]
= a[1]
The last three statements can be written more elegantly in a single line:
r,phi = complex_to_polar(z)
We can test our function by calling polar_to_comp; refer to Exercise 1.
If a function has no return statement, it returns the value None. There are many cases
where a function does not need to return any value. This could be because the variables
passed to a function may be subject to modification. Consider, for instance, the following
function:
def append_to_list(L, x):
L.append(x)
[ 156 ]
Functions
The preceding function does not return anything because it modifies one of the objects that
is given as an argument. We mentioned in Parameters and Arguments section why this is
useful. There are many methods that behave in the same way. To mention the list methods
only, the append, extend, reverse, and sort methods do not return anything (that is,
they return None ). When an object is modified by a method in this way, the modification is
called in-place. It is difficult to know whether a method changes an object, except by
looking at the code or the documentation.
Another reason for a function, or a method, not to return anything is when it prints out a
message or writes to a file.
The execution stops at the first occurring return statement. Lines after that statement are
dead code, which will never be executed:
def function_with_dead_code(x):
return 2 * x
y = x ** 2 # these two lines ...
return y
# ... are never executed!
Recursive functions
In mathematics, many functions are defined recursively. In this section, we will show how
this concept can be used even when programming a function. This makes the relation of the
program to its mathematical counterpart very clear, which may ease the readability of the
program.
Nevertheless, we recommend that you use this programming technique with care,
especially within scientific computing. In most applications, the more straightforward
iterative approach is more efficient. This will become immediately clear from the following
example.
Chebyshev polynomials are defined by a three-term recursion:
[ 157 ]
Functions
Such a recursion needs to be initialized, that is, T0(x) =1, T1(x) = x.
In Python, this three term recursion can be realized by the following function definition:
def chebyshev(n, x):
if n == 0:
return 1.
elif n == 1:
return x
else:
return 2. * x * chebyshev(n - 1, x)
- chebyshev(n - 2 ,x)
The function is then called like this:
chebyshev(5, 0.52) # returns 0.39616645119999994
This example also illustrates the risk of dramatically wasting computation time. The
number of function evaluations increases exponentially with the recursion level and most of
these evaluations are just duplicates of previous computations. While it might be tempting
to use recursive programs for demonstrating the strong relation between code and
mathematical definition, a production code will avoid this programming technique (refer to
Exercise 6). We also refer to a technique called memoization (refer to [22] for more details)
that combines recursive programming with a caching technique to save replicated function
evaluations.
A recursive function usually has a level parameter. In the previous example, it is n. It is
used to control the function's two main parts:
The base case, here, the first two if branches
The recursive body, in which the function itself is called once or several times
with smaller level parameters.
The number of levels passed by an execution of a recursive function is called the recursion
depth. This quantity should not be too large; otherwise the computation might no longer be
effective and in the ultimate case, the following error will be raised:
RuntimeError: maximum recursion depth exceeded
The maximal recursion depth depends on the memory of the computer you use. This error
also occurs when the initialization step is missing in the function definition. We encourage
the use of recursive programs for very small recursion depths (for more information, refer
to the section Infinite Iteration of Chapter 9, Iterating.
[ 158 ]
Functions
Function documentation
You should document your functions using a string at the beginning. This is
called docstring:
def newton(f, x0):
"""
Newton's method for computing a zero of a function
on input:
f (function) given function f(x)
x0 (float) initial guess
on return:
y (float) the approximated zero of f
"""
...
When calling help(newton), you get this docstring displayed together with the call of this
function:
Help on function newton in module __main__:
newton(f, x0)
Newton's method for computing a zero of a function
on input:
f (function) given function f(x)
x0 (float) initial guess
on return:
y (float) the approximated zero of f
The docstring is internally saved as an attribute, __doc__, of the given function. In the
example, it's newton.__doc__. The minimal information you should provide in a docstring
is the purpose of the function and the description of the input and output objects. There are
tools to automatically generate full code documentation by collecting all docstrings in your
program (for more information refer to [32]).
[ 159 ]
Functions
Functions are objects
Functions are objects, like everything else in Python. One may pass functions as arguments,
change their names, or delete them. For example:
def square(x):
"""
Return the square of x
"""
return x ** 2
square(4) # 16
sq = square # now sq is the same as square
sq(4) # 16
del square # square doesn't exist anymore
print(newton(sq, .2)) # passing as argument
Passing functions as arguments is very common when applying algorithms in scientific
computing. The functions fsolve in scipy.optimize for computing a zero of a given
function or quad in scipy.integrate for computing integrals are typical examples.
A function itself can have a different number of arguments with differing types. So, when
passing your function f to another function g as argument, make sure that f has exactly the
form described in the docstring of g.
The docstring of fsolve gives information about its func parameter:
func -- A Python function or method which takes at least one
(possibly vector) argument.
Partial application
Let's start with an example of a function with two variables.
The function
can be viewed as a function in two variables. Often one
considers ω not as a free variable but as a fixed parameter of a family of functions fω:
This interpretation reduces a function in two variables to a function in one variable, t, given
a fixed parameter value ω. The process of defining a new function by fixing (freezing) one
or several parameters of a function is called partial application.
[ 160 ]
Functions
Partial applications are easily created using the Python module functools, which provides
a function called partial for precisely this purpose. We illustrate this by constructing a
function that returns a sine for a given frequency:
import functools
def sin_omega(t, freq):
return sin(2 * pi * freq * t)
def make_sine(frequency):
return functools.partial(sin_omega, freq = frequency)
Using Closures
Using the view that functions are objects, partial applications can be realized by writing a
function, which itself returns a new function, with a reduced number of input arguments.
For instance, the function could be defined as follows:
def make_sine(freq):
"Make a sine function with frequency freq"
def mysine(t):
return sin_omega(t, freq)
return mysine
In this example the inner function mysine has access to the variable freq; it is neither a
local variable of this function nor is it passed to it via the argument list. Python allows such
a construction (refer to Namespaces section in Chapter 11, Namespaces, Scopes, and Modules).
Anonymous functions – the lambda
keyword
The keyword lambda is used in Python to define anonymous functions, that is; functions
without a name and described by a single expression. You might just want to perform an
operation on a function that can be expressed by a simple expression without naming this
function and without defining this function by a lengthy def block.
The name lambda originates from a special branch of calculus and
mathematical logic, the -calculus.
[ 161 ]
Functions
For instance, to compute the following expression, we may use SciPy’s function quad,
which requires the function to be integrated as its first argument and the integration bounds
as the next two arguments:
Here, the function to integrate is just a simple one-liner and we use the lambda keyword to
define it:
import scipy.integrate as si
si.quad(lambda x: x ** 2 + 5, 0, 1)
The syntax is as follows:
lambda parameter_list: expression
The definition of the lambda function can only consist of a single expression and in
particular, cannot contain loops. lambda functions are, just like other functions, objects and
can be assigned to variables:
parabola = lambda x: x ** 2 + 5
parabola(3) # gives 14
The lambda construction is always replaceable
It is important to note that lambda construction is only syntactic sugar in Python. Any
lambda construction may be replaced by an explicit function definition:
parabola = lambda x: x**2+5
# the following code is equivalent
def parabola(x):
return x ** 2 + 5
The main reason to use a construction is for very simple functions, when a full function
definition would be too cumbersome.
[ 162 ]
Functions
lambda functions provide a third way to make closures as we demonstrate by continuing
with the previous example:
We use the sin_omega function to compute the integral of the sine function for various
frequencies:
import scipy.integrate as si
for iteration in range(3):
print(si.quad(lambda x: sin_omega(x, iteration*pi), 0, pi/2.) )
Functions as decorators
In the partial application section, we saw how a function can be used to modify another
function. A decorator is a syntax element in Python that conveniently allows us to alter the
behavior of a function without changing the definition of the function itself. Let us start
with the following situation:
Assume that we have a function that determines the degree of sparsity of a matrix:
def how_sparse(A):
return len(A.reshape(-1).nonzero()[0])
This function returns an error if it is not called with an array object as input. More precisely,
it will not work with an object that does not implement the reshape method. For instance,
the how_sparse function will not work with a list, because lists have no reshape method.
The following helper function modifies any function with one input parameter so that it
tries to make a type conversion to an array:
def cast2array(f):
def new_function(obj):
fA = f(array(obj))
return fA
return new_function
[ 163 ]
Functions
Thus, the modified function how_sparse = cast2array(how_sparse) can be applied to
any object that can be cast to an array. The same functionality is achieved if the definition of
how_sparse is decorated with the type conversion function. It is recommend also to
consider the functools.wraps (refer to [8] for more details):
@cast2array
def how_sparse(A):
return len(A.reshape(-1).nonzero()[0])
To define a decorator, you need a callable object such as a function that modifies the
definition of the function to be decorated. The main purposes are:
To increase code readability by separating parts from a function that do not
directly serve its functionality (for example, memoizing)
To put common preamble and epilogue parts of a family of similar functions in a
common place (for example, type checking)
To be able to easily switch off and on additional functionalities of a function (for
example, test prints, tracing)
Summary
Functions are not only the ideal tools for making your program modular, but they also
reflect mathematical thinking. You learned the syntax of function definitions and to
distinguish between defining and calling a function.
We considered functions as objects that can be modified by other functions. When working
with functions, it is important to be familiar with the notion of the scope of a variable and
how information is passed into a function by parameters.
Sometimes, it is convenient to define functions on the fly with so-called anonymous
functions. For this, we introduced the keyword lambda.
[ 164 ]
Functions
Exercises
Ex 1 → Write a function polar_to_comp, which takes two arguments r and and returns
the complex number
Use the NumPy function exp for the exponential function.
Ex 2 → In the description of the Python module functools, (refer to [8] for more detail on
functools) you find the following Python function:
def partial(func, *args, **keywords):
def newfunc(*fargs, **fkeywords):
newkeywords = keywords.copy()
newkeywords.update(fkeywords)
return func(*(args + fargs), **newkeywords)
newfunc.func = func
newfunc.args = args
newfunc.keywords = keywords
return newfunc
Explain and test this function.
Ex 3 → Write a decorator for the function how_sparse, which cleans the input matrix A by
setting the elements that are less than 1.e-16 to zero (consider example in section Function as
decorators).
Ex 4 → A continuous function f with f(a)f(b) < 0 changes its sign in the interval [a, b] and has
at least one root (zero) in this interval. Such a root can be found with the bisection method.
This method starts from the given interval. Then it investigates the sign changes in the
subintervals,
If the sign changes in the first subinterval b is redefined to be:
Otherwise, it is redefined in the same manner to:
[ 165 ]
Functions
And the process is repeated until the b-a is less than a given tolerance.
Implement this method as a function that takes as arguments:
– the function f
– the initial interval [a, b]
– the tolerance
This function bisec should return the final interval and its midpoint.
Test the method with the function arctan and also with the polynomial f(x) = 3 x2
-5 in the interval [1.1, 1.4], and alternatively in [1.3, 1.4].
Ex. 5 → The greatest common divisor of two integers can be computed with Euclid’s
algorithm described by the following recursion:
Write a function that computes the greatest common divisor of two integers. Write another
function that computes the least common multiple of these numbers using the relation:
Ex. 6 → Study the recursive implementation of Chebyshev polynomials, consider the
example in section Recursive Function. Rewrite the program in a non-recursive way and
study computation time versus polynomial degree (see also the timeit module).
[ 166 ]
8
Classes
In mathematics, when we write sin, we refer to a mathematical object for which we know
many methods from elementary calculus. For example:
We might want to evaluate sin x at x=0.5, that is, compute sin(0.5), which returns
a real number
We might want to compute its derivative, which gives us another mathematical
object, cos
We might want to compute the first three coefficients of its Taylor polynomial
These methods may be applied not only to sin but also to other sufficiently smooth
functions. There are, however, other mathematical objects (for example, the number 5) for
which these methods make no sense. Objects that have the same methods are grouped
together in abstract classes, for example, functions. Every statement and every method
that can be applied to functions applies in particular to sin or cos. Other examples for such
classes might be a rational number, for which a denominator and numerator method exist;
an interval, which has a left and right boundary method; an infinite sequence, for which we
can ask whether it has a limit, and so on.
In this case, sin is called an instance of the class. The mathematical phrase Let g be a function
… is, in this context, called instantiation. Here, g is the name of the function; one of many
attributes that can be assigned to it. Another attribute might be its domain.
The mathematical object p(x) = 2x2– 5 is just like the sine function. Every function method
applies to p, but we can also define special methods for p. We might, for instance, ask for p’s
coefficients. These methods can be used to define the class of polynomials. As polynomials
are functions, they additionally inherit all methods of the function class.
Classes
In mathematics, we often use the same operator symbol for completely different operations.
For instance, in 5+4 and sin + cos, the operator symbol + has different meanings. By using
the same symbol, one tries to express the similarities of mathematical operations. We have
introduced these terms from object-oriented programming by applying them to
mathematical examples:
Classes
Instance and instantiation
Inheritance
Methods
Attributes
Operator overloading
In this chapter, we will show how these concepts are used in Python.
Introduction to classes
We will illustrate the concept of classes with an example of rational numbers, that is,
numbers of the form q= qN ⁄ qD, where qN and qD are integers.
Figure 8.1: An example of a class declaration
[ 168 ]
Classes
We use rational numbers here only as an example for the class concept. For future work in
Python with rational numbers use the fractions module (refer to [6]).
Class syntax
The definition of a class is made by a block command with the class keyword, the name of
the class, and some statements in the block (refer to Figure 8.1):
class RationalNumber:
pass
An instance of this class (or in other words, an object of the type RationalNumber) is
created by
r = RationalNumber()
and a query type(a) returns the answer, . If we
want to investigate whether an object is an instance of this class, we can use this:
if isinstance(a, RationalNumber):
print('Indeed it belongs to the class RationalNumber')
So far we've generated an object of the RationalNumber type, which has no data yet.
Furthermore, there are no methods defined to perform operations with these objects. This
will be the subject of the next sections.
The __init__ method
Now we provide our example class with some attributes; that is, we give it defining data. In
our case, this data will be the values of the denominator and the numerator. To this end, we
have to define a method, __init__, used to initialize the class with these values:
class RationalNumber:
def __init__(self, numerator, denominator):
self.numerator = numerator
self.denominator = denominator
[ 169 ]
Classes
Before we explain the special __init__ function, which we added to the class, we
demonstrate the instantiation of a RationalNumber object:
q = RationalNumber(10, 20)
# Defines a new object
q.numerator
# returns 10
q.denominator
# returns 20
A new object of type RationalNumber is created by using the class name as if it was a
function. This statement does two things:
It first creates an empty object, q.
Then it applies the __init__ function to it; that is, q.__init__(10, 20) is
executed.
The first parameter of __init__ refers to the new object itself. On function call, this first
parameter is replaced by the object’s instance. This applies to all methods of the class and
not only to the special method __init__. The special role of this first parameter is reflected
by the convention to name it self. In the previous example, the __init__ function defines
two attributes of the new object, numerator and denominator.
Attributes and methods
One of the main reasons for working with classes is that objects can be grouped together
and bound to a common object. We saw this already when looking at rational numbers;
denominator and numerator are two objects which we bound to an instance of the
RationalNumber class. They are called attributes of the instance. The fact that an object is
an attribute of a class instance becomes apparent from the way they are referenced, which
we have used tacitly before: