Practical 3 continued: Mathematical Functions, Aliasing and Copying

element by element, m * n
[[ 5 12]
 [21 32]]
matrix product, m @ n
[[19 22]
 [43 50]]
the top left of the product is 1*5 + 2*7 = 19

The universal functions

A universal function, or ufunc, is a mathematical function that NumPy applies to every element. There is one for everything the mathematics paper needs.

import numpy as np

a = np.array([1.0, 4.0, 9.0, 16.0])
angles = np.array([0.0, np.pi / 6, np.pi / 4, np.pi / 2])

print("sqrt     ", np.sqrt(a))
print("square   ", np.square(a))
print("exp      ", np.round(np.exp(np.array([0.0, 1.0, 2.0])), 4))
print("log      ", np.round(np.log(a), 4))
print("log10    ", np.log10(np.array([1.0, 10.0, 100.0, 1000.0])))
print("sin      ", np.round(np.sin(angles), 4))
print("cos      ", np.round(np.cos(angles), 4))
print("degrees  ", np.degrees(angles))
print("abs      ", np.abs(np.array([-3, 4, -5])))
print("floor    ", np.floor(np.array([1.2, 1.8, -1.2])))
print("ceil     ", np.ceil(np.array([1.2, 1.8, -1.2])))
print("round    ", np.round(np.array([1.24, 1.25, 1.26]), 1))
print("power    ", np.power(np.array([2, 3, 4]), 3))

sqrt      [1. 2. 3. 4.]
square    [  1.  16.  81. 256.]
exp       [1.     2.7183 7.3891]
log       [0.     1.3863 2.1972 2.7726]
log10     [0. 1. 2. 3.]
sin       [0.     0.5    0.7071 1.    ]
cos       [1.     0.866  0.7071 0.    ]
degrees   [ 0. 30. 45. 90.]
abs       [3 4 5]
floor     [ 1.  1. -2.]
ceil      [ 2.  2. -1.]
round     [1.2 1.2 1.3]
power     [ 8 27 64]

np.round(x, 4) is used above so the output fits the page and does not depend on how many digits a version chooses to print. In your own journal print them unrounded once, so you have seen the full values.

Look at the round row: 1.24 gave 1.2, 1.26 gave 1.3, and 1.25 gave 1.2, not 1.3. That is not a bug. NumPy, like Python's own round, rounds a value sitting exactly halfway to the nearest even last digit, so 1.25 goes down to 1.2 and 2.25 would go up to 2.2 as well. It is called round half to even, and it exists because always rounding halves upward makes a long column of figures drift high. Expect it, and do not treat it as an error in your program.

sin 30 degrees is 0.5, and the second value in the sin row is exactly that, which is the check that the angles were given in radians. NumPy's trigonometric functions take radians, never degrees. np.radians() converts if your input is in degrees.

Aggregate functions, and the axis

An aggregate reduces many numbers to one. On a two dimensional array you can reduce down the columns, along the rows, or over everything.

import numpy as np

marks = np.array([[78, 65, 80],
                  [55, 72, 60],
                  [90, 88, 95],
                  [40, 51, 45]])

print(marks)
print("total of everything      ", marks.sum())
print("mean of everything       ", round(marks.mean(), 4))
print("axis=0, down the columns ", marks.sum(axis=0))
print("axis=1, along the rows   ", marks.sum(axis=1))
print("each student's mean      ", marks.mean(axis=1))
print("each subject's highest   ", marks.max(axis=0))
print("each subject's lowest    ", marks.min(axis=0))
print("standard deviation       ", round(marks.std(), 4))
print("row 2 is the best student, index", marks.sum(axis=1).argmax())

[[78 65 80]
 [55 72 60]
 [90 88 95]
 [40 51 45]]
total of everything       819
mean of everything        68.25
axis=0, down the columns  [263 276 280]
axis=1, along the rows    [223 187 273 136]
each student's mean       [74.33333333 62.33333333 91.         45.33333333]
each subject's highest    [90 88 95]
each subject's lowest     [40 51 45]
standard deviation        17.5979
row 2 is the best student, index 2

axis=0 goes down, axis=1 goes across. That is the fact to memorise, and here is the way to keep it straight: the axis you name is the one that disappears. marks has shape (4, 3); marks.sum(axis=0) removes the 4 and leaves 3 numbers, one per column.

Check the first column by hand: 78 + 55 + 90 + 40 = 263, and that is the first number in the axis=0 row.

Broadcasting

Two arrays of different shapes can still be combined, if one of them can be stretched to fit. That stretching is called broadcasting, and it is what makes a * 2 work: the 2 is broadcast to every element.

The rule, in one line: compare the shapes from the right, and a dimension of 1, or a missing dimension, is stretched.

import numpy as np

marks = np.array([[78, 65, 80],
                  [55, 72, 60],
                  [90, 88, 95]])
weights = np.array([2, 1, 1])

print("marks shape  ", marks.shape)
print("weights shape", weights.shape)
print("weighted marks, each column scaled")
print(marks * weights)
print("add 5 to everybody")
print(marks + 5)
print("subtract each subject's mean from its column")
print(np.round(marks - marks.mean(axis=0), 4))

marks shape   (3, 3)
weights shape (3,)
weighted marks, each column scaled
[[156  65  80]
 [110  72  60]
 [180  88  95]]
add 5 to everybody
[[ 83  70  85]
 [ 60  77  65]
 [ 95  93 100]]
subtract each subject's mean from its column
[[  3.6667 -10.       1.6667]
 [-19.3333  -3.     -18.3333]
 [ 15.6667  13.      16.6667]]

And when the shapes cannot be made to fit, NumPy says so rather than guessing:

import numpy as np

a = np.array([[1, 2, 3], [4, 5, 6]])
b = np.array([1, 2])
print(a + b)

ValueError: operands could not be broadcast together with shapes (2,3) (2,)

Read that message: it names both shapes. (2, 3) and (2,) compared from the right gives 3 against 2, and neither is 1, so there is nothing to stretch. b of shape (3,) would have worked.

Part two: aliasing, views and copying

This is MU's third bullet and it is the part of this practical that gets asked about. There are three cases, not two, and they have to be told apart by writing through each one and looking at the original.

An alias is a second name

import numpy as np

a = np.array([1, 2, 3, 4, 5])
b = a

print("a is", a)
print("b is", b)
print("are they the same object?", b is a)
print("same id?", id(a) == id(b))

b[0] = 999
print("after b[0] = 999")
print("a is", a)
print("b is", b)

a is [1 2 3 4 5]
b is [1 2 3 4 5]
are they the same object? True
same id? True
after b[0] = 999
a is [999   2   3   4   5]
b is [999   2   3   4   5]

There is one array in that program and two names for it. b = a copies nothing at all; it writes the same reference into a second name. So of course changing one changes the other, and b is a proves they are the same object.

This is not special to NumPy. A Python list behaves identically, and it is the source of one of the nastiest bugs a beginner writes:

first = [1, 2, 3]
second = first
second.append(4)
print("first ", first)
print("second", second)
print("the same object?", second is first)

first  [1, 2, 3, 4]
second [1, 2, 3, 4]
the same object? True

A view is a window

import numpy as np

a = np.array([1, 2, 3, 4, 5, 6])
window = a[1:4]
reshaped = a.reshape(2, 3)

print("a       ", a)
print("window  ", window, " is it a view?", window.base is not None)
print("reshaped\n", reshaped, " is it a view?", reshaped.base is not None)

window[0] = 999
print("after window[0] = 999, a is", a)

reshaped[1, 2] = 777
print("after reshaped[1, 2] = 777, a is", a)
print("are they the same object? ", window is a, reshaped is a)

a        [1 2 3 4 5 6]
window   [2 3 4]  is it a view? True
reshaped
 [[1 2 3]
 [4 5 6]]  is it a view? True
after window[0] = 999, a is [  1 999   3   4   5   6]
after reshaped[1, 2] = 777, a is [  1 999   3   4   5 777]
are they the same object?  False False

Note the last line: a view is not the same object as the array, so is says False, and yet writing through it still changes the array. That is exactly why is is not the test here. base is.

A copy is a new array

import numpy as np

a = np.array([1, 2, 3, 4, 5])

shallow_name = a
sliced_view = a[:]
real_copy = a.copy()
also_a_copy = np.array(a)

for label, other in [("b = a          ", shallow_name),
                     ("a[:]           ", sliced_view),
                     ("a.copy()       ", real_copy),
                     ("np.array(a)    ", also_a_copy)]:
    print(f"{label} same object {other is a!s:>5}   owns its data {other.base is None!s:>5}")

real_copy[0] = 999
also_a_copy[1] = 888
print("after writing 999 and 888 into the two copies, a is", a)

sliced_view[2] = 777
print("after writing 777 into a[:], a is         ", a)

b = a           same object  True   owns its data  True
a[:]            same object False   owns its data False
a.copy()        same object False   owns its data  True
np.array(a)     same object False   owns its data  True
after writing 999 and 888 into the two copies, a is [1 2 3 4 5]
after writing 777 into a[:], a is          [  1   2 777   4   5]

There is the whole exercise in one output. a.copy() and np.array(a) own their data and leave a alone. a[:] does not, and this is the line that catches people. On a Python list a[:] is the classic way to make a copy. On a NumPy array it is a view. A student who carries the list habit across corrupts their data silently.

Deep copy, and when you actually need it

a.copy() is enough for any array of numbers. It is not enough for an array or a list holding other containers, because copying the outer one still leaves the inner ones shared.

from copy import deepcopy

original = [[1, 2], [3, 4]]
shallow = list(original)
deep = deepcopy(original)

shallow[0][0] = 999
print("after writing through the shallow copy:", original)

deep[1][1] = 888
print("after writing through the deep copy   :", original)
print("the deep one changed only itself      :", deep)

after writing through the shallow copy: [[999, 2], [3, 4]]
after writing through the deep copy   : [[999, 2], [3, 4]]
the deep one changed only itself      : [[1, 2], [3, 888]]

list(original) made a new outer list, but its two slots still refer to the same two inner lists, so writing into shallow[0][0] reached the original. deepcopy copied all the way down. For NumPy numeric arrays this never arises, because the numbers are stored in the array itself rather than referred to.

The one line summary to say at the table

import numpy as np

a = np.arange(6)

for label, other in [("alias  b = a   ", a),
                     ("view   a[1:4]  ", a[1:4]),
                     ("copy   a.copy()", a.copy())]:
    print(f"{label}  is a: {other is a!s:>5}"
          f"   shares memory: {np.shares_memory(other, a)!s:>5}")

alias  b = a     is a:  True   shares memory:  True
view   a[1:4]    is a: False   shares memory:  True
copy   a.copy()  is a: False   shares memory: False

is a true means an alias. False but sharing memory means a view, and false and sharing nothing means a copy. Two tests, three answers.

One caution about base, because [Practical 4: NumPy Slicing, Basic and Advanced Indexing] meets it head on. base says whether an array owns its own buffer, which for a plain slice or a plain copy is the same question as "is this a view of that array". For a subscript that mixes a slice with a list of positions it is not the same question, and base is then not None even though nothing is shared. Where the answer matters, ask np.shares_memory(x, a).

Procedure

1. Save as practical3b.py. Apply +, -, *, /, //, % and ** to an array and

to two arrays, and print each result.

1. Show that list 2 repeats and array 2 multiplies.
2. Apply sqrt, exp, log, sin, cos, abs, floor, ceil and power, and check

that sin of pi over 6 is 0.5.

1. Build a two dimensional array of marks. Take sum, mean, max and min over

everything, then with axis=0 and axis=1. Check one column total by hand.

1. Show broadcasting with a row of weights, and show the shape error when it cannot work.
2. Save as practical3c.py. Make an alias, a view and a copy of one array. For each, print

whether it is the same object and whether it shares memory with the original, then write through it and print the original.

1. Show that a[:] copies a list but views an array.

Result

Every operator and nine mathematical functions were applied element by element; sin of pi over 6 came out 0.5, confirming radians. Column totals matched the hand check 78 + 55 + 90 + 40 = 263. Broadcasting worked for shapes (3, 3) and (3,) and raised ValueError for (2, 3) and (2,). Writing through the alias and through the view both changed the original; writing through a.copy() and np.array(a) did not; and a[:] behaved as a view on the array and as a copy on the list.

Where marks are lost

• Using a loop to apply a function to an array. The whole point is that you do not.
• Calling a * b matrix multiplication. It is element by element. @ is the product.
• Giving degrees to np.sin. It takes radians.
• Getting axis backwards. The axis you name is the one that disappears.
• Calling a view a copy. a[1:4] is a view, and writing to it changes a.
• Using a[:] to copy an array out of list habit. On an array it is a view.
• Testing with is to tell a view from a copy. A view is a different object; use

np.shares_memory.

• Showing only the code. Aliasing is only proved by printing the original after the

write.

For the journal

Two entries under practical 3. For the mathematical functions: the aim, a table of the operators with one example each, the ufunc program and its output, and the axis table of marks with one column total checked by hand. One sentence on broadcasting: shapes are compared from the right and a dimension of 1 is stretched.

For aliasing and copying: the aim, the three row table of alias, view and copy, then the program that writes through each one, and the original printed after each write. Record np.shares_memory for each of the three, because that is the line that proves the table. The conclusion in one sentence: b = a makes a second name, a slice makes a window, and only copy() makes a new array, so the first two change the original and the third does not.

Quick revision

• An operation on an array applies to every element. No loop.
• list 2 repeats the list; array 2 multiplies each element.
• Between two arrays every operator is element by element. a * b is NOT the matrix

product; a @ b is.

• ufuncs: sqrt, square, exp, log, log10, sin, cos, abs, floor, ceil,

round, power. Trigonometry takes radians; np.radians() converts.

• Aggregates: sum, mean, min, max, std, argmax.

axis=0 goes down the columns and axis=1 along the rows, and the axis named is the one that disappears.

• Broadcasting: compare shapes from the right; a 1 or a missing dimension stretches.

Otherwise ValueError, and the message names both shapes.

• Alias b = a: one array, two names. b is a is True.
• View a[1:4], a.reshape(), a.T: different object, same memory.
• Copy a.copy() or np.array(a): new memory.
• The test is np.shares_memory(x, a). x.base is a quick indicator only: on a mixed

subscript such as m[1:, [0, 2]] the base is a temporary, so it is not None while nothing is shared.

• Writing through an alias or a view changes the original. Through a copy it does not.
• a[:] copies a list and views an array.
• deepcopy is only needed for containers holding containers.

Questions you should be able to answer

1. Why does an array not need a loop to add 1 to every element? Because arithmetic on a NumPy array is applied element by element by the library itself, in compiled code. That is called vectorising.

2. What does [1, 2] 3 give, and what does np.array([1, 2]) 3 give? [1, 2, 1, 2, 1, 2], because multiplying a list repeats it, and [3 6], because multiplying an array scales every element.

3. What is the difference between a b and a @ b? a b multiplies element by element. a @ b is the matrix product.

4. marks has shape (4, 3). What does marks.sum(axis=0) give? An array of shape (3,): one total per column, that is, per subject. The axis named is the one that disappears.

5. State the broadcasting rule. Compare the shapes from the right. A dimension that is 1, or absent, is stretched to match. If two dimensions differ and neither is 1, it is an error.

6. What are the three ways a second name can relate to an array? An alias, which is the same object; a view, a different object sharing the same memory; and a copy, a different object with its own memory.

7. How do you tell a view from a copy at the keyboard? np.shares_memory(x, a), which is True exactly when writing to one is seen in the other. x.base is a quick indicator of whether x owns its buffer, but on a mixed subscript the base is a temporary array, so it can be not None while nothing at all is shared.

8. Why is b = a never a copy? Because assignment binds a name to the object that is already there. Nothing is duplicated, so both names reach one array.

9. A student writes backup = a[:] and then changes backup. What happens? The array a changes too, because on a NumPy array a[:] is a view. On a Python list the same line would have made a real copy, which is why the mistake is so easy.

10. When do you need deepcopy rather than copy? Only when the container holds other containers. copy duplicates the outer one and leaves the inner ones shared. An array of numbers never needs it.