Basic numerics and plot#

pylab notebook just like namespace

# For interactive plot, please use 
# %pylab notebook
%pylab inline
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 3
      1 # For interactive plot, please use 
      2 # %pylab notebook
----> 3 get_ipython().run_line_magic('pylab', 'inline')

File ~/opt/anaconda3/envs/book_new/lib/python3.9/site-packages/IPython/core/interactiveshell.py:2417, in InteractiveShell.run_line_magic(self, magic_name, line, _stack_depth)
   2415     kwargs['local_ns'] = self.get_local_scope(stack_depth)
   2416 with self.builtin_trap:
-> 2417     result = fn(*args, **kwargs)
   2419 # The code below prevents the output from being displayed
   2420 # when using magics with decodator @output_can_be_silenced
   2421 # when the last Python token in the expression is a ';'.
   2422 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):

File ~/opt/anaconda3/envs/book_new/lib/python3.9/site-packages/IPython/core/magics/pylab.py:155, in PylabMagics.pylab(self, line)
    151 else:
    152     # invert no-import flag
    153     import_all = not args.no_import_all
--> 155 gui, backend, clobbered = self.shell.enable_pylab(args.gui, import_all=import_all)
    156 self._show_matplotlib_backend(args.gui, backend)
    157 print(
    158     "%pylab is deprecated, use %matplotlib inline and import the required libraries."
    159 )

File ~/opt/anaconda3/envs/book_new/lib/python3.9/site-packages/IPython/core/interactiveshell.py:3642, in InteractiveShell.enable_pylab(self, gui, import_all, welcome_message)
   3615 """Activate pylab support at runtime.
   3616 
   3617 This turns on support for matplotlib, preloads into the interactive
   (...)
   3638     This argument is ignored, no welcome message will be displayed.
   3639 """
   3640 from IPython.core.pylabtools import import_pylab
-> 3642 gui, backend = self.enable_matplotlib(gui)
   3644 # We want to prevent the loading of pylab to pollute the user's
   3645 # namespace as shown by the %who* magics, so we execute the activation
   3646 # code in an empty namespace, and we update *both* user_ns and
   3647 # user_ns_hidden with this information.
   3648 ns = {}

File ~/opt/anaconda3/envs/book_new/lib/python3.9/site-packages/IPython/core/interactiveshell.py:3588, in InteractiveShell.enable_matplotlib(self, gui)
   3567 def enable_matplotlib(self, gui=None):
   3568     """Enable interactive matplotlib and inline figure support.
   3569 
   3570     This takes the following steps:
   (...)
   3586         display figures inline.
   3587     """
-> 3588     from matplotlib_inline.backend_inline import configure_inline_support
   3590     from IPython.core import pylabtools as pt
   3591     gui, backend = pt.find_gui_and_backend(gui, self.pylab_gui_select)

File ~/opt/anaconda3/envs/book_new/lib/python3.9/site-packages/matplotlib_inline/__init__.py:1
----> 1 from . import backend_inline, config  # noqa
      2 __version__ = "0.1.6"  # noqa

File ~/opt/anaconda3/envs/book_new/lib/python3.9/site-packages/matplotlib_inline/backend_inline.py:6
      1 """A matplotlib backend for publishing figures via display_data"""
      3 # Copyright (c) IPython Development Team.
      4 # Distributed under the terms of the BSD 3-Clause License.
----> 6 import matplotlib
      7 from matplotlib import colors
      8 from matplotlib.backends import backend_agg

ModuleNotFoundError: No module named 'matplotlib'

Commands beginning with % our IPython ‘magic’ commands. This one sets up a bunch of matplotlib back end and imports numpy into the global namespace. We will do this in all of our notebooks for simplicity. It’s actually considered bad form because it puts a lot of numpy functions into the global namespace, but for exploratory work it’s very convenient.

To avoid importing everything into the global namespace, one would use instead

#%matplotlib notebook
import numpy as np
import matplotlib as mpl
#mpl.rcParams['figure.dpi']= 150
%config InlineBackend.figure_format = 'svg'
import matplotlib.pyplot as plt

which will make all of the numpy functions, such as array() and sin(), available as np.array(), np.sin(), and the plotting functions such as plot() and xlabel() as plt.plot() and plt.xlabel().

Numpy Arrays#

Numpy arrays store multidimensional arrays of objects of a fixed type. The type of an array is a dtype, which is a more refined typing system than Python provides. They are efficient maps from indices (i,j) to values. They have minimal memory overhead.

Arrays are mutable: their contents can be changed after they are created. However, their size and dtype, once created cannot be efficiently changed (requires a copy).

Arrays are good for:

  1. Representing matrices and vectors (linear algebra)

  2. Storing grids of numbers (plotting, numerical analysis)

  3. Storing data series (data analysis)

  4. Getting/changing slices (regular subarrays)

Arrays are not good for:

  1. Applications that require growing/shrinking the size.

  2. Heterogenous objects.

  3. Non-rectangular data.

Arrays are 0-indexed.

# A vector is an array with 1 index
a = array([1/sqrt(2), 0, 1/sqrt(2)])
a
array([0.70710678, 0.        , 0.70710678])
a.shape
(3,)
a.dtype
dtype('float64')
a.size
3
# We access the elements using [ ]
a[0]
0.7071067811865475
a[0] = a[0]*2
a
array([1.41421356, 0.        , 0.70710678])

We create a 2D array (that is a matrix) by passing the array() function a list of lists of numbers in the right shape.

# A matrix is an array with 2 indices
B = array( [[ 1, 0, 0],
            [ 0, 0, 1],
            [ 0, 1, 0]] )
B
array([[1, 0, 0],
       [0, 0, 1],
       [0, 1, 0]])
B.shape
(3, 3)
B.dtype
dtype('int64')
B.size
9
B[0,0]
1

Warning! There is also a type called ‘matrix’ instead of ‘array’ in numpy. This is specially for 2-index arrays but is being removed from Numpy over the next two years because it leads to bugs. Never use matrix(), only array()

Basic Linear Algebra#

There are two basic kinds of multiplication of arrays in Python:

  1. Element-wise multiplication: a*b multiplies arrays of the same shape element by element.

  2. Dot product: a@b forms a dot product of two vectors or a matrix product of two rectangular matrices.

Mathematically, for vectors,

\[ a@b = \sum_i a[i] b[i] \]

while for 2D arrays (matrices),

\[ A@B[i,j] = \sum_k A[i,k] B[k,j] \]
a = array([1,  1]) / sqrt(2)
b = array([1, -1]) / sqrt(2)
a
array([0.70710678, 0.70710678])
b
array([ 0.70710678, -0.70710678])
a*a
array([0.5, 0.5])
a@a
0.9999999999999998
# Compute the length of a
norm(a)
0.9999999999999999
sqrt(a@a)
0.9999999999999999
a*b
array([ 0.5, -0.5])
a@b
0.0

There are many, many more functions for doing linear algebra operations numerically provided by numpy and scipy. We will use some of them as we go.

Basic Plotting#

The second primary use of numpy array’s is to hold grids of numbers for analyzing and plotting. In this case, we consider a long 1D array with length N as representing the values of the x and y axis of a plot, for example.

Let’s plot a sine wave:

# create an equally spaced array of 100 numbers 
# from -2pi to 2pi 
x = np.linspace(-2*np.pi, 2*np.pi, 100)

# evaluate a function at each point in x and create a 
# corresponding array
y = 0.5*np.sin(x)
figure()
plot(x,y)
grid()
xlabel(r'$x$')
ylabel(r'$0.5 \sin(x)$')
show()
_images/55c2e2ad0b00d22306b82b6abf5ba6cf6ed64dc112c68b19b00aa12beaf5024c.svg

A few comments:

  1. The call to sin(x) is a ‘ufunc’, which automatically acts element-by-element on whatever shape array it is passed.

  2. Matplotlib supports \(\LaTeX\) style mathematical expressions in any text that it renders – just enclose in $ -signs. We use a raw (r”..”) string so that the backslashes are passed onto the \(\LaTeX\) interpreter intact rather than being interpreted as special characters.

  3. The ‘notebook’ backend (which we turned on at the beginning of the notebook) provides basic interactive plotting inside the notebook. Try it out. There are other backends which provide somewhat better interaction if you setup jupyter on your local computer.

Speed#

L = range(1000)
%timeit [i**2 for i in L]# % is a comment communicate to jupyter, not to python.
279 µs ± 4.72 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
a = arange(1000)
%timeit a**2
806 ns ± 35.5 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
mpl.rcParams['text.usetex']=True # This assigned the font style to be the MathJax style in the plot.
mpl.rcParams['font.family']='Times New Roman'
figure()
plot(a, a**2, '--')
title(r'This is a test.')
xlabel(r'$a$')
ylabel(r'$a^2$',rotation=0)# Do not rotate the y axis label.
show()
_images/0498cef04d99770efa635746e7b7c7ccc2dad315daa5dbdc5cb35451cc861508.svg

Some Common Arrays#

arange(2,10,2)
array([2, 4, 6, 8])
# compact notation for the previous
# r_ creates a *row* vector with the contents of the slice 
# notice the [ ] rather than ( )
r_[2:10:2]
array([2, 4, 6, 8])
linspace(2,10,5)
array([ 2.,  4.,  6.,  8., 10.])
# compact notation for the previous
r_[2:10:5j]
array([ 2.,  4.,  6.,  8., 10.])
ones((2,2))
array([[1., 1.],
       [1., 1.]])
zeros((3,1))
array([[0.],
       [0.],
       [0.]])
eye(3)
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])
diag([1,2,3])
array([[1, 0, 0],
       [0, 2, 0],
       [0, 0, 3]])
np.random.rand(2,2)
array([[0.80947482, 0.89261187],
       [0.021396  , 0.07445528]])
np.random.exponential(scale=3,size=(3,3))
array([[ 1.72812182, 11.29459947,  4.62183267],
       [ 2.47499787,  2.67545379,  3.34511434],
       [ 0.75359228,  3.57557972,  0.46732403]])

Array DTypes#

Numpy has its own set of data types for the elements of arrays. These dtypes are more specific than ‘int’ or ‘str’ so that you can control the number of bytes used.

  1. Integers: int16 (‘i2’), int32 (‘i4’), int64 (‘i8’), …

  2. Unsigned: uint32 (‘u4’), uint64 (‘u8’), …

  3. Float: float16 (‘f2’), float32 (‘f4’), float64 (‘f8’), …

  4. Boolean: bool

  5. Fixed Length Strings: ‘S1’, ‘S2’, …

array([1,0]).dtype
dtype('int64')
array([1.,0]).dtype
dtype('float64')
float32
numpy.float32
dtype('i4')
dtype('int32')

Basic Operations and Visualization#

x = linspace(-pi, pi, 100)
y1 = sin(x)
y2 = exp(x/pi)
figure()
plot(x,y1, label=r'$\sin(x)$')
plot(x,y2, label=r'$e^{x/\pi}$')
legend(loc='upper left')
show()
_images/428022ed399ee20b9d3546769dbfc0abce68b8518a3f0f44629fa323acba350b.svg