import numpy as np
import matplotlib.pyplot as plt
Vectorized calculations
We will rather loosely be using the term "vectorized" to mean array-at-a-time computations. True (hardware-level) vectorization actually runs a small set of multiple calculations at the same time on your device - this may or may not happen in our array computations - and it isn't important at the moment.
a = np.array([1, 2, 3, 4])
b = np.array([5, 6, 7, 8])
# Classic loops:
c1 = np.empty_like(a)
for i in range(len(a)):
c1[i] = a[i] + b[i]
# Array-at-a-time, or "vectorized":
c2 = a + b
print(c1)
print(c2)
Let's try in two diminsions:
Note: since I'm explicitly looking at the output values, I'm okay to use square matricies. It's often a nice check to use non-square matricies if you have the choice; it catches missmatches in diminsions!
a = np.array([[1, 2], [3, 4]])
b = np.array([[5, 6], [7, 8]])
# Classic loops:
c1 = np.empty_like(a)
for i in range(a.shape[0]):
for j in range(a.shape[1]):
c1[i, j] = a[i, j] + b[i, j]
# Array-at-a-time, or "vectorized":
c2 = a + b
print(c1)
print(c2)
np.array([1, 2]) + 1
That expands to multiple diminsions like this:
(a, b, ..., z) OP (a, b, ..., z) -> element wise
(a, b, ..., z) OP (a, 1, ..., z) -> length 1 instead of b -> duplicated b times
(a, b, ..., z) OP (b, ..., z) -> prepend a diminsion of length 1 -> see aboveThe third rule explains scalars, but is a bit problematic. Try to use the second rule when possible instead.
print(np.ones([2, 3]) + np.ones([2, 3]))
print(np.ones([2, 3]) + np.ones([1, 3]))
print(np.ones([2, 3]) + np.ones([3]))
Memory layout
Computer memory is linear. It might be fragmented a bit, but it is not 2D or 3D. (and if it was, it would be even harder to manage). So, all objects more than 1D have to be laid out in some fasion. There are three choices for 2D arrays:
- Pointers to 1D arrays. Requires look ups to traverse - bad for performance and memory.
- Row major C order. Default in Python. Similar to the way we write.
- Column major F order (from Fortran). Default in Matlab, optional in Python.
Law: whenever there are two ways to do something, it will always be done both ways.
A = np.array([[1, 2], [3, 4]], order="C")
print(A)
print(A.flatten("K"))
A = np.array([[1, 2], [3, 4]], order="F")
print(A)
print(A.flatten("K"))
This matters when we make loops - how do we want the computer to look up the items in the list?
A = np.array([[1, 2], [3, 4]])
print("Continious:", end=" ")
for i in range(2):
print(A[0, i], end=" ")
print()
print("Not continious:", end=" ")
for i in range(2):
print(A[i, 0], end=" ")
Numpy hides the details whenever it can. Even if you use Fortran style arrays, everything will still "just work"; unless you try hard, even "flatten" works the same way - it turns it into a C style flat array even if that's not how it's laid out in memory (add "K" to fix). Any method or function that depends on the order will have an order= argument. UFuncts will loop in the best order automatically. And Transposes are basically no-ops - they just change from C style to F style and vise versa!
print(np.isfortran(A))
print(np.isfortran(A.T))
Note: we'll discuss why A.T is the transpose of a matrix in a minute.
np.ones([3]).flatten("K")
Adding a 1 length dimension is basically a "no-op" - the memory of the array does not change, just the interpretation of it.
So why is the auto up-scaling broadcasting (rule 3 above) irritating to use?
Given these two brodcasting expressions:
$$ \left[ \begin{matrix} 1 & 2 \\ 3 & 4 \end{matrix} \right] \times \left[ \begin{matrix} 1 & 2 \end{matrix} \right] $$$$ \left[ \begin{matrix} 1 & 2 \\ 3 & 4 \end{matrix} \right] \times \left[ \begin{matrix} 1 \\ 2 \end{matrix} \right] $$What happens if the second one is just a 1D array?
a = np.array([[1, 2], [3, 4]])
print(a * np.array([[1, 2]]))
print(a * np.array([[1], [2]]))
print(a * np.array([1, 2]))
A = np.array([[1, 2], [3, 4]])
A.T
Aside:
Python objects (since Python 2.2, at least) provide properties - those are methods that look at at like members!
So the following:
value = object.member
can secretly call a function! (Think object.member(), though it's technically calling object.member.__get__().) A similar trick exists for setting, as well.
Why would you use them?
- You can turn a real member into a property as your class gets more complex
- You can avoid some
()(great for chaining) - You can take advantage of the nice set syntax
When would you use them?
- Not for long calculations - a user does not expect a program to hang when accessing or setting a member
- Only for simple properties that do not have any options
You can access the transpose of a matrix:
a.T
Since a 1D transpose does nothing, the following expression does not do what you might naively expect:
# print(a * np.array([1,2]).T)
The solution? Make this an explicit row or column vector. Generally, when doing linear algebra, you should be using 2D arrays. This does not mean that Python is wrong for having 1D arrays - it just adds flexibility.
Quick way to make a column vector:
print(np.array([[1, 2]]).T)
Interesting note: The transpose actually just converts an array between C ordering and F ordering!
A.flags
A.T.flags
We have points in 3D space with a weight W. What is the weighted mean?
state = np.random.RandomState(42)
x = state.rand(10)
y = state.rand(10)
z = state.rand(10)
w = state.rand(10)
xsum = 0
ysum = 0
zsum = 0
wsum = 0
for i in range(len(x)):
xsum += x[i] * w[i]
ysum += y[i] * w[i]
zsum += z[i] * w[i]
wsum += w[i]
xsum /= wsum
ysum /= wsum
zsum /= wsum
fig, axs = plt.subplots(1, 3, figsize=(14, 4))
axs[0].scatter(x, y, s=w * 1000)
axs[0].plot(xsum, ysum, "rx")
axs[1].scatter(x, z, s=w * 1000)
axs[1].plot(xsum, zsum, "rx")
axs[2].scatter(y, z, s=w * 1000)
axs[2].plot(ysum, zsum, "rx")
plt.show()
# state = np.random.RandomState(42)
# M = state.rand(4, 10)
M = np.stack([x, y, z])
M.shape
Wsum = np.sum(w)
Msum = np.sum(M * w, axis=1) / Wsum
Msum.shape
print(Msum == np.array([xsum, ysum, zsum]))
M1 = np.stack([x, y, z, w])
state = np.random.RandomState(42)
M2 = state.rand(4, 10)
print(M1 == M2)
We can also do this with a 4x10 array.
Wsum2 = np.sum(M2[3])
Msum2 = np.sum(M2[:3] * M2[3], axis=1) / Wsum
Msum2.shape
Msum2 == np.array([xsum, ysum, zsum])
This (as we've seen) is nicer to generate, and keeps our number of variables down, but has the issue that it becomes less clear to read. We can use one of two tricks:
- Set
w = M2[3]andM=[:3]. These will still be the same memory space (no copy is made - they are just "views"), so this does not come at a penalty. - Use a record array, or Pandas. We'll cover this (Pandas, at least) later. If the types of the columns do not match, this is the only way to do it.