<jupyter_start><jupyter_text>Numpy<jupyter_code>import numpy as np<jupyter_output><empty_output><jupyter_text>The core of the `numpy` package is the `array` class. Let's examine that first. We can make an array out of a sequence, like a list.<jupyter_code>d = [1, 2, 3, 4, 5]
np.array(d)<jupyter_output><empty_output><jupyter_text>data typesUnlike lists, arrays must be homogeneous, in that the data types of each element must be the same. The data type of the array is upcast to be able to represent all of the data. So, if only one element is a float, all elements will be converted to floats.<jupyter_code>d = [1, 2, 3.1415, 4, 5]
np.array(d)<jupyter_output><empty_output><jupyter_text>You can query the datatype by examining the `dtype` attribute of the array.<jupyter_code>d = [1, 2, 3.1415, 4, 5]
arr = np.array(d)
arr.dtype<jupyter_output><empty_output><jupyter_text>Array types may be defined explicity in the call<jupyter_code>arr = np.array([1, 2, 3, 4, 5], dtype='float32')
arr<jupyter_output><empty_output><jupyter_text>Complex numbers are noted with a lowercase `j` or uppercase `J`, like this<jupyter_code>cmplx = np.array([1.0+2.0j, 3.0+4.0J])
print(cmplx)
cmplx.dtype<jupyter_output>[ 1.+2.j  3.+4.j]<jupyter_text>As we have seen before, arrays are like multidimensional sequences. We can create a 2D array by supplying a list of lists as the argument.<jupyter_code>arr = np.array([[1., 2., 3.,], [4., 5., 6.]])
arr<jupyter_output><empty_output><jupyter_text>Array attributesArrays have a few other important attributes. Note attributes never have parentheses after them. Methods always do.<jupyter_code>arr.size          # The number of elements in the array
arr.shape         # The shape of the array (i.e., the size of each dimension)
arr.ndim          # The number of dimensions of the array<jupyter_output><empty_output><jupyter_text>Setting array shapeYou can set the `array.shape` attribute to change the shape of the array. This attribute does not change the elements of the array, or how it is stored in memory, just how it is seen.<jupyter_code>arr.shape = (3, 2)
arr
arr.shape = (6,)
arr<jupyter_output><empty_output><jupyter_text>Singleton dimensions add to the dimensionality of an array. The last example was a 1D array (also called a vector), the next are 2D arrays.<jupyter_code>arr.shape = (1, 6)
arr   # Note that there are *two* square brackets in the output sequence. This is a row vector.
arr.shape = (6, 1)
arr   # this is also a 2D array, like a column vector<jupyter_output><empty_output><jupyter_text>Array indexingArrays are indexed in a similar way to sequences, with `start:stop:stride` notation, except that this is used for each dimension in the array. Colons denote all the values in a particular dimension, slices indicate some particular subset of the data in that particular dimension. A common use case is to get a single row or column from a 2D array (a table of data).<jupyter_code>arr = np.arange(60).reshape(6, 10)
arr
arr[:, 4]   # the 5th column
arr[2, :]   # the 3rd row
arr[2]     # Trailing colons do not need to be explicitly typed. This is equivalent to the last example.
arr[4, 7]   # an individual element in the table<jupyter_output><empty_output><jupyter_text>--- *Exercise*> Slices can be combined in any way. Define a new array or use array `arr` and grab out every other row and the 4th column and beyond.---  Conventions concerning arrays containing spatio-temporal informationGenerally, you will want to think of arrays as representing dimensions in space and time. The conventional way to think of this is that the dimensions are $(t, z, y, x)$; missing dimensions are omitted. This will help make plotting and analysis easier. Some examples might be:    temp[:, :, :, :]      A 4D array (time, height, latitude, longitude)    press[:, :]           A 2D array (time, height)    humid[:, :]           A 2D array (latitude, longitude)  Array methodsArrays have a number of methods. Let's take a look at the `mean` method as an example.<jupyter_code>arr = np.array([[1., 2., 3.,], [4., 5., 6.]])  # reset the array to our 2x3 array.

arr.mean()        # The mean of all of the elements in the array<jupyter_output><empty_output><jupyter_text>`Mean` takes the optional argument `axis` that can be used to take the mean along a single axis of the array. Just like with indexing, the axes are reference in a zero-based system; `axis=0` means the first dimension.<jupyter_code>arr.mean(axis=0)  # The mean<jupyter_output><empty_output><jupyter_text>In this case, there are two rows in the first dimension, and `arr.mean(axis=0)` takes the average in the 'row' direction, resulting in a 1D array that is the average across the rows. --- *Exercise*> Find the mean of the array in the 'column' direction, along `axis=1`.> Use the `sum` method of the array class to get the sum of the numbers in each column. The result should be a 1D array with three elements.--- You can also use the `reshape` method to change the shape of an array.<jupyter_code>arr
arr.reshape(3, 2)<jupyter_output><empty_output><jupyter_text>You can find the mininum and maximum of an array with the `min` and `max` methods. Sometimes it is useful to find the indices of these minima and maxima. For this use `argmin` and `argmax`, like<jupyter_code>x = np.random.rand(10)
imax = x.argmax()
print(imax, x[imax], x.max())<jupyter_output>7 0.947778556457 0.947778556457<jupyter_text>Creating standard arraysThere are a few standard arrays, for example, arrays filled with zeros or ones (or empty). Here are some examples of creating arrays.<jupyter_code>o = np.ones((3, 4, 5))    # The argument is a shape, so is a tuple with the length of each dimension as an argument
b = np.ones((2, 3), dtype=np.bool)
z = np.zeros((2, 3), dtype=np.float32)

b<jupyter_output><empty_output><jupyter_text>You can also create these arrays with the same shape and datatype of the input array using `np.ones_like` and `np.zeros_like`.<jupyter_code>zb = np.zeros_like(b)
zb<jupyter_output><empty_output><jupyter_text>You can also create a diagonal array with a given vector along the diagonal. These can be offset with an optional argument `k` (default=0). This example creates a tri-diagonal array similar to that used for finite difference calculations<jupyter_code>np.diag(-2*np.ones(6)) + np.diag(np.ones(5), k=-1) + np.diag(np.ones(5), k=1)<jupyter_output><empty_output><jupyter_text>There are also a number of ways to generate sequences of numbers. - `np.arange([start,] stop [[, stride]])` Create a sequence of numbers, similar to `range` - `np.linspace(min, max, length)` Create a uniform series of specified `length` between `min` and `max`, inclusive. - `np.logspace(minpow, maxpow, length)` Create a uniform series in logspace of specified `length` between `10**minpow` and `10**maxpow`, inclusive.<jupyter_code>np.arange(10.)
np.arange(2, 10, 2)
np.linspace(2, 4, 17)<jupyter_output><empty_output><jupyter_text>You can create arrays of random numbers easily with methods in `np.random`.* `np.random.rand(d0, d1, ..., d2)`: Create an array of the given shape `d0, ..., dn` and populate it with random samples from a uniform distribution over [0,1).* `np.random.randint(low, high=None, size=None)`: Return random integers from `low` (inclusive) to `high` (exclusive). If `high` is None then return integers from [0, `low`). `size` is an int or tuple of ints to give the output shape.* `np.randon.randn(d0, d1, ..., d2)`: Create an array of the given shape `d0, ..., dn` and populate it with random samples from the "standard normal" distribution.* `np.random.random(size=None)`: Return random floats of `size` (int or tuple of ints) in the interval [0, 1).<jupyter_code>np.random.rand(2, 4)
np.random.randint(1, 50, (2, 4))<jupyter_output><empty_output><jupyter_text>--- *Exercise*> Create an array of random floats between 0 and 1 that has dimension 5 x 3. Calculate the standard deviation of the columns of the array. Then add to this a `linspace` array of the appropriate size that contains numbers between 10 and 15.---  Combining and splitting arraysGenerally, arrays can be combined with the `np.concatenate` function. The arguments are a sequence of arrays to join, and the axis along which to join them (default=0).<jupyter_code>x = np.random.rand(4, 5, 6)
y = np.random.rand(4, 5, 6)

print(np.concatenate((x, y)).shape)
print()
print(np.concatenate((x, y), axis=0).shape)
print(np.concatenate((x, y), axis=1).shape)
print(np.concatenate((x, y), axis=2).shape)<jupyter_output>(8, 5, 6)

(8, 5, 6)
(4, 10, 6)
(4, 5, 12)<jupyter_text>There are a number of convenience functions that act like concatenate for specific axes: - `np.vstack` – vertical stack (stack along axis=0) - `np.hstack` – horizontal stack (stack along axis=1) - `np.dstack` – depth stack (stack along axis=2)<jupyter_code>print(np.vstack((x, y)).shape)
print(np.hstack((x, y)).shape)
print(np.dstack((x, y)).shape)<jupyter_output>(8, 5, 6)
(4, 10, 6)
(4, 5, 12)<jupyter_text>Likewise, arrays can be split with `np.split` or `np.array_split`. There are also convenience functions to split horizontally, vertically, and with depth.<jupyter_code>x = np.random.rand(12, 2, 5)
[a.shape for a in np.split(x, 4, axis=0)]<jupyter_output><empty_output><jupyter_text>--- *Exercise*> Create an array, A, of shape (40, 50, 60). The array slices for first ten entries in the axis=1 direction of A should be filled with 1's, for the next ten filled with 2's, and on up to 5's.> Split it along axis=1 into five sections.> Concatenate two of these back together along axis 1.> What is the resulting shape of each array?  _[Advanced: can you calculate this on one line?]_---  Finding valuesThere are a number of ways to find values in an array. The simplest is always to create a boolean array, like<jupyter_code>x = np.random.rand(5, 5)
print(x)
ind = x > 0.5
print(ind)<jupyter_output>[[ 0.11873289  0.40188555  0.96426409  0.60360203  0.54526375]
 [ 0.20289134  0.98457289  0.3847002   0.00733166  0.63739985]
 [ 0.50234835  0.25930096  0.36977497  0.05787961  0.95391643]
 [ 0.02016065  0.69361985  0.5945571   0.381457    0.33079702]
 [ 0.62571904  0.18626591  0.35990447  0.65804042  0.12489787]]
[[False False  True  True  True]
 [False  True False False  True]
 [ True False False False  True]
 [False  True  True False False]
 [ True False False  True False]]<jupyter_text>The boolean array can be used as an index to other arrays. Note this will return a 1D array, no matter what dimension the origial arrays are, because there is no way to know what structure the `True` values have.<jupyter_code>x = np.random.rand(5, 5)
y = np.sin(x)

y[x > 0.5]
# or, equivalently, as two lines
idx = x > 0.5
y[idx]<jupyter_output><empty_output><jupyter_text>To get the indices of the places where the conditional is true (i.e., the locations of the `True` values in the boolean array), use the `np.where` command.<jupyter_code>x = np.random.rand(5, 5)
idx = np.where(x > 0.5)
idx<jupyter_output><empty_output><jupyter_text>Note that `np.where` always returns a tuple of indices for each dimension. This is a little strange for 1D arrays, but is done for consistency across all input values. Often, you will want to explicitly pull out the (single) array of indices from the tuple, like<jupyter_code>x = np.random.rand(10)
idx = np.where(x>0.5)[0]
print(idx)<jupyter_output><empty_output><jupyter_text>_What happens with the [0] is missing behind the call to `where`?_ --- *Exercise*> You can also use these calculated indices, or boolean matrices on the left hand side for assignment.> Create a 10x10 random array, with values between 0 and 1. Replace all of the numbers smaller than 0.5 with zero.> Do this first not using `where` and then do it using `where`.---  Array viewsThe data for an array may be stored in memory using `C` or `FORTRAN` ordered memory. Typically, there is no need to think about this, some details can be found [here](http://docs.scipy.org/doc/numpy-1.10.0/reference/internals.html).However, it is important to remember that subsets of an array can produce a different 'view' of the array that addresses the same memory as the original array. This can lead to some unexpected behaviors. One way to think of this is that assignment in Python is more like a C-pointer (i.e., a reference to a memory location) than an actual value.<jupyter_code>a = np.arange(10.0)
b = a[::2]
print(a)
print(b)
a[4] = -999   # this will modify b as well...
print(a)
print(b)
b[-1] = -888  # this will modify a as well...
print(a)
print(b)<jupyter_output>[   0.    1.    2.    3. -999.    5.    6.    7. -888.    9.]
[   0.    2. -999.    6. -888.]<jupyter_text>Normally, this will not be a problem, but if you need to make sure that a subset of an array has it's own memory, make sure you make a `copy` of the array, like<jupyter_code>a = np.arange(10.0)
b = a.copy()[::2]     # or np.copy(a)
a[4] = -999   # this will NOT modify b now
print(a)
print(b)<jupyter_output>[   0.    1.    2.    3. -999.    5.    6.    7.    8.    9.]
[ 0.  2.  4.  6.  8.]<jupyter_text>Array broadcasting(Largely taken from [SciPy docs](https://docs.scipy.org/doc/numpy-1.10.0/user/basics.broadcasting.html))Generally arrays should be the same shape for them to be multiplied together.<jupyter_code>a = np.array([1.0, 2.0, 3.0])
b = np.array([2.0, 2.0, 2.0])
a * b<jupyter_output><empty_output><jupyter_text>The term broadcasting describes how `numpy` treats arrays with different shapes during arithmetic operations. Subject to certain constraints, the smaller array is “broadcast” across the larger array so that they have compatible shapes. Broadcasting provides a means of vectorizing array operations so that looping occurs in C instead of Python. It does this without making needless copies of data and usually leads to efficient algorithm implementations.For example, the simplest broadcasting example occurs when an array and a scalar value are combined in an operation:<jupyter_code>a = np.array([1.0, 2.0, 3.0])
b = 2.0
a * b<jupyter_output><empty_output><jupyter_text>The result is equivalent to the previous example where b was an array. We can think of the scalar b being stretched during the arithmetic operation into an array with the same shape as a. The new elements in b are simply copies of the original scalar. The stretching analogy is only conceptual. NumPy is smart enough to use the original scalar value without actually making copies, so that broadcasting operations are as memory and computationally efficient as possible.The code in the second example is more efficient than that in the first because broadcasting moves less memory around during the multiplication (b is a scalar rather than an array).  General Broadcasting RulesWhen operating on two arrays, NumPy compares their shapes element-wise. It starts with the trailing dimensions, and works its way forward. Two dimensions are compatible when1. they are equal, or1. one of them is 1If these conditions are not met, a ValueError: frames are not aligned exception is thrown, indicating that the arrays have incompatible shapes. The size of the resulting array is the maximum size along each dimension of the input arrays.Arrays do not need to have the same number of dimensions. For example, if you have a 256x256x3 array of RGB values, and you want to scale each color in the image by a different value, you can multiply the image by a one-dimensional array with 3 values. Lining up the sizes of the trailing axes of these arrays according to the broadcast rules, shows that they are compatible:    Image  (3d array): 256 x 256 x 3    Scale  (1d array):             3    Result (3d array): 256 x 256 x 3When either of the dimensions compared is one, the other is used. In other words, dimensions with size 1 are stretched or “copied” to match the other.In the following example, both the A and B arrays have axes with length one that are expanded to a larger size during the broadcast operation:    A      (4d array):  8 x 1 x 6 x 1    B      (3d array):      7 x 1 x 5    Result (4d array):  8 x 7 x 6 x 5    Here are some more examples:    A      (2d array):  5 x 4    B      (1d array):      1    Result (2d array):  5 x 4    A      (2d array):  5 x 4    B      (1d array):      4    Result (2d array):  5 x 4    A      (3d array):  15 x 3 x 5    B      (3d array):  15 x 1 x 5    Result (3d array):  15 x 3 x 5    A      (3d array):  15 x 3 x 5    B      (2d array):       3 x 5    Result (3d array):  15 x 3 x 5    A      (3d array):  15 x 3 x 5    B      (2d array):       3 x 1    Result (3d array):  15 x 3 x 5    Let's create an example with arrays of random numbers.<jupyter_code>A = np.random.rand(15, 3, 5)
B = np.random.rand(3, 1)
print(A.shape, B.shape)

Result = A * B
print(Result.shape)<jupyter_output>(15, 3, 5) (3, 1)
(15, 3, 5)<jupyter_text>Here are examples of shapes that do not broadcast:    A      (1d array):  3    B      (1d array):  4  trailing dimensions do not match    A      (2d array):      2 x 1    B      (3d array):  8 x 4 x 3  second from last dimensions mismatched --- *Exercise*    a = np.random.rand(5, 7, 1, 8)    b = np.random.rand(8)    c = np.random.rand(7, 3, 8)    d = np.random.rand(5, 1, 3, 1)> Experiment with multiplying combinations of the arrays above together. Try to predict the resulting shape beforehand.--- Notice that the rules for broadcasting are based on the location of singleton dimensions. Singleton dimensions are implied forward (to the left), but not backward (to the right). So, the first example here works but not the second:    A      (2d array):  5 x 4    B      (1d array):      4    Result (2d array):  5 x 4    A      (2d array):  5 x 4    B      (1d array):  5*Compare with large set of examples above. How can the bottom example here be fixed?* This problem can be fixed by creating new singleton dimensions in arrays. This can be done by putting `np.newaxis` in the appropriate space when indexing the array. For example:<jupyter_code>A = np.random.rand(5, 4)
B = np.random.rand(5)
A*B
print(B.shape)
print(B[:,np.newaxis].shape)
(A*B[:,np.newaxis]).shape<jupyter_output><empty_output><jupyter_text>--- *Exercise*> Multiply `b = np.random.rand(8)` and `c = np.random.rand(8, 3, 7)`. What is another way you could accomplish this calculation besides using `newaxis`?--- --- *Exercise*    b = np.random.rand(2)    c = np.random.rand(2, 3)    > Concatenate arrays `b` and `c`. Along which axis would it make sense to concatenate, given the arrays dimensions? Do you need to make any changes to the arrays to get this to work?---  Flattening arrays with `a.flat` and `a.flatten()`There are two basic ways to turn any array into a 1D array. They are slightly different.`a.flatten()` returns a copy of an array, in one dimension.<jupyter_code>a = np.arange(12).reshape(3, 4)
print(a)
b = a.flatten()
print(b)<jupyter_output>[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
[ 0  1  2  3  4  5  6  7  8  9 10 11]<jupyter_text>the `flat` attribute on the other hand gives a view of the array in 1D. It looks like an iterator object (like `range` and `zip`). This allows<jupyter_code>a.flat[6] = -999
print(a)<jupyter_output>[[   0    1    2    3]
 [   4    5 -999    7]
 [   8    9   10   11]]<jupyter_text>In contrast, this does not work as expected.  _WHY?_<jupyter_code>a.flatten()[5] = -888
print(a)<jupyter_output>[[   0    1    2    3]
 [   4    5 -999    7]
 [   8    9   10   11]]<jupyter_text>Other operations can be done to the array first. For example, we can take a transpose of the array before we flatten it.<jupyter_code>a.T.flat[6] = -998
print(a)<jupyter_output>[[   0    1 -998    3]
 [   4    5 -999    7]
 [   8    9   10   11]]<jupyter_text>Here, the `T` attribute (equivalent to the `a.transpose()` method) gives a view of the array transposed (similar to MATLAB's tick notation).<jupyter_code>print(a.T)<jupyter_output>[[   0    4    8]
 [   1    5    9]
 [-998 -999   10]
 [   3    7   11]]<jupyter_text>Masked arraysMasked arrays are ways to create arrays with missing values. MATLAB&trade; uses NaNs (NaN stands for 'Not a Number'), and the NaNs are the values of the arrays at those points. This approach also works in Python. Masked arrays are preferred since they retain the masked array values, and also some plotting routines require masked arrays when plotting arrays with missing values. Masked arrays are usually created through some condition, like<jupyter_code>arr = np.random.randn(7, 8)
cond = arr > 0.1   # `cond` is True for the random values greater than 0.5

marr = np.ma.masked_where(cond, arr)

print(marr)
marr.mean(axis=0)<jupyter_output><empty_output><jupyter_text>The mask can also be supplied explicity when creating the masked array,<jupyter_code>marr = np.ma.masked_array([1, 2, 3, 4, 5], mask=[True, True, False, False, True])
marr<jupyter_output><empty_output><jupyter_text>Importing dataOne of the basic commands in `numpy` for loading in data is the `loadtxt` command. There are other ways to do this, such as the [`genfromtxt`](http://docs.scipy.org/doc/numpy-dev/user/basics.io.genfromtxt.html) command, but `loadtxt` is sufficient for most purposes, and is easy to use.<jupyter_code>data = np.loadtxt('../data/CTD.txt', comments='*')
data[:,2]    # a column of data representing temperature
data<jupyter_output><empty_output><jupyter_text>--- *Exercise*> Read in the oceanographic data file '../data/CTD.txt' into an array. You can look at the data file itself to see what variables are stored in each column.> Using this data, write a function to calculate the linear equation of state. This is an approximation of the density of water, as it depends on salinity, temperature, and some empirical constants. We will use the following form for the linear equation of state:> $\rho = 1027[1+7.6\times 10^{-4}(S-35) -1.7\times 10^{-4}(T-25)]$> where $\rho$ is the density, $S$ is the salinity, and $T$ is the temperature.> This is more free form than the homework, so you should set up all of the associated code to call the function, and write out the function yourself. Don't forget docstrings! For a check, the first value of your density array in order should equal 1021.7519981630001 and the last should equal 1028.0471353619998.---  Polynomial fittingThe basic function for fitting a polynomial (e.g., a straight line) is `np.polyfit(x, y, deg)`. There are a number of other functions that let you add (`np.polyadd`), multiply (`np.polymul`), find zeros (`np.roots`), and do other operations to polynomials.<jupyter_code>x = np.random.rand(100)
y = 5 + 3*x + 0.1*np.random.randn(100)   # A straight line with some noise

p = np.polyfit(x, y, 1)  # fit a straight line (order is 1)
print(p)  # The coefficients of the polynomial, with highest order first. (i.e,. [slope, intercept])<jupyter_output>[ 3.04891792  4.98836372]<jupyter_text>Let's plot it to make sure this makes sense:<jupyter_code>import matplotlib.pyplot as plt
%matplotlib inline

# plot data
plt.plot(x, y, '.')

# plot fitted line
plt.plot(x, p[0]*x + p[1])

plt.legend(('Data', 'Fitted line'))<jupyter_output><empty_output><jupyter_text>Once you have the fit, you can use it to find other useful things, like the value of the fitted line at $x=1$:<jupyter_code>np.polyval(p, 1)<jupyter_output><empty_output><jupyter_text>You can also use the `np.polynomial.Polynomial` class to work with polynomials. Note, these define polynomials the opposite way, with the _lowest_ order first. The Polynomial class gives an excellent example of operator overloading, and the flexibility of classes.<jupyter_code>p1 = np.polynomial.Polynomial([5, 3])         # y = 5 + 3 x
p2 = np.polynomial.Polynomial([3, 6, 8, 2])   # y = 3 + 6 x + 8 x**2 + 2 x**3<jupyter_output><empty_output><jupyter_text>You can use the Polynomial object to evaluate the value of the polynomial at various input values:<jupyter_code>print('Evaluation')
print('p1(0.0) = ', p1(0))
print('p2(5.0) = ', p2(5))<jupyter_output>Evaluation
p1(0.0) =  5.0
p2(5.0) =  483.0<jupyter_text>We can use this to make a plot to see the function:<jupyter_code>x = np.linspace(0,10)
plt.plot(x, p1(x), x, p2(x))
plt.legend(['p1', 'p2'])<jupyter_output><empty_output><jupyter_text>Other things we can do:<jupyter_code>print('Roots')
print('Roots of p2 = ', p2.roots())
print()
print('Operations')
print('p1 + p2 = ', p1 + p2)
print('p1 * p2 = ', p1 * p2)
print()
print('Calculus')
print('Derivative of p1', p1.deriv(1))
print('Integral of p2', p2.integ(4, k=[4, 3, 2, 1]))<jupyter_output>Roots
Roots of p2 =  [-3.21124270+0.j         -0.39437865-0.55818847j -0.39437865+0.55818847j]

Operations
p1 + p2 =  poly([ 8.  9.  8.  2.])
p1 * p2 =  poly([ 15.  39.  58.  34.   6.])

Calculus
Derivative of p1 poly([ 3.])
Integral of p2 poly([ 1.          2.          1.5         0.66666667  0.125       0.05
  0.02222222  0.00238095])<jupyter_text>VectorizationVectorization and array broadcasting are two big reasons that `numpy` can be efficient and fast. With these tools, you can avoid writing for loops (which are slow).The best way to do mathematical operations using `numpy` arrays is to do `vector` operations. That is, mathematical operations are defined to be element by element, and this is done much faster than looping. As a rule of thumb, you should be very concerned if your code has more than one significant `for` loop in the numerical analysis section.Here is a way to do multiply 2 big arrays using for loops, which is not how you should do it. The sum at the end is included for comparison with the subsequent approach.<jupyter_code>a = np.arange(102400.0).reshape(4, 8, 1600, 2)   # a 4D array using sequential numbers
b = np.random.rand(4, 8, 1600, 2)                # a 4D array using random numbers

li, lj, lk, lm = b.shape  # size of b in each dimension
sol = np.zeros(b.shape)
for i in range(li):
    for j in range(lj):
        for k in range(lk):
            for m in range(lm):
                sol[i,j,k,m] = a[i,j,k,m]*b[i,j,k,m]
print(sol.sum())<jupyter_output>2625932085.42<jupyter_text>The better way is to directly multiply the arrays together, taking advantage of C code that Python has in the background.<jupyter_code>sol = a * b       # element-by-element multiplication. This operation is about as fast as it can be on your computer.
print(sol.sum())<jupyter_output>2625932085.42<jupyter_text>Basic performance evaluationWe can do some very basic perfomance testing using the `%time` special function in jupyter notebooks. Lets use this to examine the time it takes to do a singular value decomposition for different sized matrices.<jupyter_code>b = np.random.randn(5000, 2000)
%time u, s, v = np.linalg.svd(b)<jupyter_output>CPU times: user 45.2 s, sys: 275 ms, total: 45.5 s
Wall time: 3.99 s<jupyter_text>`%time` runs the line once and gives the time required. However, calculation times vary depending on many things including the numbers involved and the state of your computer at the moment. In this case, the `%timeit` function can be used to perform the test a number of times to get an average calculation time.<jupyter_code>%timeit b = np.random.randn(50, 20); u, s, v = np.linalg.svd(b)<jupyter_output>The slowest run took 4.78 times longer than the fastest. This could mean that an intermediate result is being cached.
1000 loops, best of 3: 444 µs per loop<jupyter_text>For statements that are longer than a single line, the `time.time` function can be used.<jupyter_code>import time

t_start = time.time()
time.sleep(0.25)   # Do nothing for 0.25 seconds
t_stop = time.time()

print('{:6.4f} seconds have passed.'.format(t_stop-t_start))<jupyter_output>0.2506 seconds have passed.<jupyter_text>--- *Exercise*> Earlier, we discussed using array operations instead of looping because it is faster. Let's compare.> Calculate the time it takes to calculate the $a$ and $b$ arrays with dimensions [4, 8, 1600, 2] by both methods demonstrated: using a series of 4 `for` loops, one for each dimension of the arrays and using array operations. Compare the times by calculating a ratio.---  Linear algebraOne of the key elements of the `numpy` package is the `numpy.linalg` subpackage that contains a number of linear algebra functions that work efficiently on arrays.<jupyter_code>a = np.random.randn(100, 100)
e, v = np.linalg.eig(a)

b = np.random.randn(500, 200)
u, s, v = np.linalg.svd(b)<jupyter_output><empty_output><jupyter_text>Matrix multiplication is done using the `np.dot` function. In this case, matrices do _not_ need to be the same shape, but must follow the rules of matrix multiplication. E.g., the operation dot(, ) results in a 4x12 array; i.e., the inner dimensions must match (technically last and second-to-last, for arrays with more than two dimensions).<jupyter_code>x = np.random.rand(4, 5)
y = np.random.rand(5, 12)

res = np.dot(x, y)
print(res.shape)

# np.dot(y, x)  # This gives an error -- order is important.<jupyter_output>(4, 12)
StarcoderdataJupyterStructuredCleanDedup