hans

hans

【Python】2. Numpy - "Scientific Computing with Python"


In standard installations of Python, a list can be used to store a group of values, which can act as an array. However, since the elements of a list can be any object, what is stored in the list are pointers to those objects. Thus, to store a simple [1,2,3], three pointers and three integer objects are needed. This structure is clearly wasteful in terms of memory and CPU computation time for numerical operations.
Additionally, Python provides an array module, where array objects differ from lists in that they directly store values, similar to one-dimensional arrays in C. However, since it does not support multi-dimensional arrays and lacks various operation functions, it is also not suitable for numerical computations.
The birth of NumPy addresses these shortcomings, providing two basic objects: ndarray (N-dimensional array object) and ufunc (universal function object). ndarray (hereafter referred to as array) is a multi-dimensional array that stores a single data type, while ufunc is a function that can process arrays.#

Here, I would like to explain the differences between lists, tuples, and arrays:
List: Elements can be any type of object, must be initialized, and can be modified after initialization. When all element objects are the same, it can be used as an array.

>>> a = ['a', "hello", 1, 2.2]
>>> a
['a', 'hello', 1, 2.2]
>>> a[1] = "HELLO"
>>> a
['a', 'HELLO', 1, 2.2]

Tuple: Elements can be any type of object, must be initialized, and cannot be modified after initialization.

>>> b = ('a', "hello", 1, 2.2)
>>> b
('a', 'hello', 1, 2.2)
>>> b[2] = "HELLO"
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: 'tuple' object does not support item assignment

Array: Elements can only be of one type of object, initialization is not mandatory, and can be modified after initialization.

>>> c = np.array([1,2,3])
>>> d = np.array(['a','b','c'])
>>> e = np.array([1,2,'c','d'])
>>> e
array(['1', '2', 'c', 'd'],
      dtype='|S21')

As you can see, in array e, 1 and 2 are defaulted to string types to match the subsequent 'c' and 'd'.
If you observe carefully, you can see that tuples use parentheses (), while lists and arrays use square brackets [].#

1. Importing the Library#

import numpy as np

Import the numpy library and use the abbreviation np thereafter.

2. Creation#

a and b are one-dimensional arrays, while c is a multi-dimensional array.

>>> a = np.array([1, 2, 3, 4])
>>> b = np.array((5, 6, 7, 8))
>>> c = np.array([[1, 2, 3, 4],[4, 5, 6, 7], [7, 8, 9, 10]])
>>> b
array([5, 6, 7, 8])
>>> c
array([[1, 2, 3, 4],
       [4, 5, 6, 7],
       [7, 8, 9, 10]])
>>> c.dtype
dtype('int32')

The size of an array can be obtained through its shape attribute:

>>> a.shape
(4,)
>>> c.shape
(3, 4)

Array a's shape has only one element, so it is a one-dimensional array. Array c's shape has two elements, so it is a two-dimensional array, where the length of the 0th axis (3 horizontal axes) is 3, and the length of the 1st axis (4 vertical axes) is 4. You can also change the length of each axis of the array by modifying the shape attribute while keeping the number of array elements unchanged. The following example changes the shape of array c to (4,3). Note that changing from (3,4) to (4,3) is not transposing the array, but simply changing the size of each axis; the positions of the array elements in memory do not change:

>>> c.shape = 4,3
>>> c
array([[ 1,  2,  3],
       [ 4,  4,  5],
       [ 6,  7,  7],
       [ 8,  9, 10]])

Looking at the above, array a is a one-dimensional array, and a.shape has only one number, which indicates the number of elements in the array. We can also express a using the multi-dimensional array expression method. Note the difference:

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

When an axis has an element of -1, the length of that axis will be automatically calculated based on the number of array elements, so the following program changes the shape of array c to (2,6):

>>> c.shape = 2,-1
>>> c
array([[ 1,  2,  3,  4,  4,  5],
       [ 6,  7,  7,  8,  9, 10]])

Using the reshape method of the array, you can create a new array with changed dimensions, while the shape of the original array remains unchanged:

>>> d = a.reshape((2,2))
>>> d
array([[1, 2],
       [3, 4]])
>>> a
array([1, 2, 3, 4])

Arrays a and d actually share the same data storage memory area, so modifying any element of one array will simultaneously modify the contents of the other array:

>>> a[1] = 100 # Change the first element of array a to 100
>>> d # Note that the 2 in array d has also changed
array([[  1, 100],
       [  3,   4]])

The element type of an array can be obtained through the dtype attribute. In the above example, all elements of the parameter sequence are integers, so the created array's element type is also integer and is a 32-bit long integer. You can specify the element type when creating it through the dtype parameter:

>>> np.array([[1, 2, 3, 4],[4, 5, 6, 7], [7, 8, 9, 10]], dtype=np.float)
array([[  1.,   2.,   3.,   4.],
       [  4.,   5.,   6.,   7.],
       [  7.,   8.,   9.,  10.]])
>>> np.array([[1, 2, 3, 4],[4, 5, 6, 7], [7, 8, 9, 10]], dtype=np.complex)
array([[  1.+0.j,   2.+0.j,   3.+0.j,   4.+0.j],
       [  4.+0.j,   5.+0.j,   6.+0.j,   7.+0.j],
       [  7.+0.j,   8.+0.j,   9.+0.j,  10.+0.j]])

The above examples all create arrays using the array function. Below are three functions specifically designed to create arrays.
The arange function is similar to Python's range function, creating a one-dimensional array by specifying the start value, end value, and step size. Note that the array does not include the end value:

>>> np.arange(0,1,0.1)
array([ 0. ,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9])

The linspace function creates a one-dimensional array by specifying the start value, end value, and number of elements, and you can specify whether to include the end value using the endpoint keyword; the default setting is to include the end value:

>>> np.linspace(0, 1, 12)
array([ 0.        ,  0.09090909,  0.18181818,  0.27272727,  0.36363636,
        0.45454545,  0.54545455,  0.63636364,  0.72727273,  0.81818182,
        0.90909091,  1.        ])

The logspace function is similar to linspace, but it creates a geometric sequence. The following example generates a geometric sequence from 1 (10^0) to 100 (10^2) with 20 elements:

>>> np.logspace(0, 2, 20)
array([   1.        ,    1.27427499,    1.62377674,    2.06913808,
          2.6366509 ,    3.35981829,    4.2813324 ,    5.45559478,
          6.95192796,    8.8586679 ,   11.28837892,   14.38449888,
         18.32980711,   23.35721469,   29.76351442,   37.92690191,
         48.32930239,   61.58482111,   78.47599704,  100.        ])

Additionally, functions like frombuffer, fromstring, fromfile can create arrays from byte sequences. Below is an example using fromstring:

>>> s = "abcdefgh"

Python strings are actually byte sequences, with each character occupying one byte. Therefore, if we create an 8-bit integer array from string s, the resulting array is exactly the ASCII encoding of each character in the string:

>>> np.fromstring(s, dtype=np.int8)
array([ 97,  98,  99, 100, 101, 102, 103, 104], dtype=int8)

If we create a 16-bit integer array from string s, then two adjacent bytes represent one integer, treating byte 98 and byte 97 as one 16-bit integer, its value is 98*256+97 = 25185. It can be seen that the data is stored in memory in little-endian format (least significant byte first).

>>> np.fromstring(s, dtype=np.int16)
array([25185, 25699, 26213, 26727], dtype=int16)
>>> 98*256+97
25185

We can write a Python function that converts array indices into corresponding values in the array, and then use this function to create an array:

>>> def func(i): # Define function, function name func, parameter i
...   return i%4+1 # Return the remainder of i divided by 4 plus one
...
>>> np.fromfunction(func, (10,))
array([ 1.,  2.,  3.,  4.,  1.,  2.,  3.,  4.,  1.,  2.])

The second parameter n of fromfunction can be understood as having n values starting from 0, assigned to the preceding function.
In this example, i first takes 0, and after taking the remainder and adding one, returns the value 1.
Then i takes 1, and after taking the remainder and adding one, returns the value 2.
And so on.
Finally, when i takes 9, it returns the value 2 after taking the remainder and adding one.#

The first parameter of the fromfunction function is the function that calculates each array element, and the second parameter is the size of the array (shape). Since it supports multi-dimensional arrays, the second parameter must be a sequence; in this example, we create a one-dimensional array with 10 elements using (10,).

The following example creates a two-dimensional array representing the multiplication table, where each element a[i, j] in the output array a equals func2(i, j):

>>> def func2(i, j): # Define function, function name func, parameters i, j
...     return (i+1) * (j+1) # Function return value
...
>>> a = np.fromfunction(func2, (9,9))
>>> a
array([[  1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.],
       [  2.,   4.,   6.,   8.,  10.,  12.,  14.,  16.,  18.],
       [  3.,   6.,   9.,  12.,  15.,  18.,  21.,  24.,  27.],
       [  4.,   8.,  12.,  16.,  20.,  24.,  28.,  32.,  36.],
       [  5.,  10.,  15.,  20.,  25.,  30.,  35.,  40.,  45.],
       [  6.,  12.,  18.,  24.,  30.,  36.,  42.,  48.,  54.],
       [  7.,  14.,  21.,  28.,  35.,  42.,  49.,  56.,  63.],
       [  8.,  16.,  24.,  32.,  40.,  48.,  56.,  64.,  72.],
       [  9.,  18.,  27.,  36.,  45.,  54.,  63.,  72.,  81.]]) # Can be understood as the product of the indices at each point

3. Accessing Elements#

The method of accessing array elements is the same as the standard method in Python:

>>> a = np.arange(10)
>>> a[5]    # Using an integer as an index can retrieve a specific element in the array
5
>>> a[3:5]  # Using a range as an index retrieves a slice of the array, including a[3] but excluding a[5]
array([3, 4])
>>> a[:5]   # Omitting the starting index means starting from a[0]
array([0, 1, 2, 3, 4])
>>> a[:-1]  # Indices can use negative numbers, indicating counting backwards from the last element
array([0, 1, 2, 3, 4, 5, 6, 7, 8])
>>> a[2:4] = 100,101    # Indices can also be used to modify element values
>>> a
array([  0,   1, 100, 101,   4,   5,   6,   7,   8,   9])
>>> a[1:-1:2]   # The third parameter in the range indicates the step size; 2 means taking every other element
array([  1, 101,   5,   7])
>>> a[::-1] # Omitting the starting and ending indices with a step size of -1 reverses the entire array
array([  9,   8,   7,   6,   5,   4, 101, 100,   1,   0])
>>> a[5:1:-2] # When the step size is negative, the starting index must be greater than the ending index
array([  5, 101])

Unlike Python's list sequences, the new array obtained by accessing elements through index ranges is a view of the original array. It shares the same data space with the original array:

>>> b = a[3:7] # A new array b is created through the index range, which shares the same data space with a
>>> b
array([101,   4,   5,   6])
>>> b[2] = -10 # Modifying the 2nd element of b to -10
>>> b
array([101,   4, -10,   6])
>>> a # The 5th element of a has also been modified to -10
array([  0,   1, 100, 101,   4, -10,   6,   7,   8,   9])

In addition to accessing elements using index ranges, NumPy also provides two advanced methods for accessing elements.
Using Integer Sequences
When accessing array elements using an integer sequence, each element in the integer sequence is used as an index. The integer sequence can be a list or an array. The array obtained by using an integer sequence as an index does not share data space with the original array.

>>> x = np.arange(10,1,-1)
>>> x
array([10,  9,  8,  7,  6,  5,  4,  3,  2])
>>> x[[3, 3, 1, 8]] # Retrieve four elements from x at indices 3, 3, 1, 8, forming a new array
array([7, 7, 9, 2])
>>> b = x[np.array([3,3,-3,8])]  # Indices can be negative
>>> b[2] = 100
>>> b
array([7, 7, 100, 2])
>>> x   # Since b and x do not share data space, the values in x remain unchanged
array([10,  9,  8,  7,  6,  5,  4,  3,  2])
>>> x[[3,5,1]] = -1, -2, -3 # Integer sequence indices can also be used to modify element values
>>> x
array([10, -3,  8, -1,  6, -2,  4,  3,  2])

Using Boolean Arrays
When using a boolean array b as an index to access elements in array x, all elements in x corresponding to True indices in b will be collected. The array obtained by using a boolean array as an index does not share data space with the original array. Note that this method only corresponds to boolean arrays and cannot use boolean lists. (Do you remember the differences between lists, tuples, and arrays?)

>>> x = np.arange(5,0,-1)
>>> x
array([5, 4, 3, 2, 1])
>>> x[np.array([True, False, True, False, False])]
>>> # In the boolean array, indices 0 and 2 are True, so we retrieve elements at indices 0 and 2 from x
array([5, 3])
>>> x[[True, False, True, False, False]]
>>> # If it is a boolean list, it treats True as 1 and False as 0, retrieving elements from x in integer sequence manner
array([4, 5, 4, 5, 5])
>>> x[np.array([True, False, True, True])]
>>> # If the boolean array is not long enough, the missing parts are treated as False
array([5, 3, 2])
>>> x[np.array([True, False, True, True])] = -1, -2, -3
>>> # Boolean array indices can also be used to modify elements
>>> x
array([-1,  4, -2, -3,  1])

Boolean arrays are generally not generated manually but are produced using boolean operations with ufunc functions. For more information on ufunc functions, please refer to the ufunc operations section.

>>> x = np.random.rand(10) # Generate an array of length 10 with random values between 0 and 1
>>> x
array([ 0.72223939,  0.921226  ,  0.7770805 ,  0.2055047 ,  0.17567449,
        0.95799412,  0.12015178,  0.7627083 ,  0.43260184,  0.91379859])
>>> x>0.5
>>> # Each element in array x is compared with 0.5, resulting in a boolean array where True indicates that the corresponding value in x is greater than 0.5
array([ True,  True,  True, False, False,  True, False,  True, False,  True], dtype=bool)
>>> x[x>0.5]
>>> # Using the boolean array returned by x>0.5 to collect elements from x, thus obtaining an array of all elements in x greater than 0.5
array([ 0.72223939,  0.921226  ,  0.7770805 ,  0.95799412,  0.7627083 ,
        0.91379859])

4. Multi-dimensional Arrays#

Accessing multi-dimensional arrays is similar to one-dimensional arrays. Since multi-dimensional arrays have multiple axes, their indices need to be represented by multiple values. NumPy uses tuples as array indices. As shown in the figure, a is a 6x6 array, with colors distinguishing various indices and their corresponding selection areas.
1668631158233.jpg


The statement a[2::2,::2] can be broken down as follows: first, "2:" indicates starting from index 2 of the horizontal axis (0th axis) to the end, then ":2" indicates taking every other row. In fact, it retrieves rows 2 and 4. The first ":" indicates all columns (1st axis), and ":2" indicates taking every other column, which means taking columns 0, 2, and 4.

For example, the index of 24 is (2,4), where 2 is the index on the 0th axis, and 4 is the index on the 1st axis.


Array a is an addition table, with values on the vertical axis being 0, 10, 20, 30, 40, 50; and values on the horizontal axis being 0, 1, 2, 3, 4, 5. Each element on the vertical axis is summed with each element on the horizontal axis.

>>> np.arange(0, 60, 10).reshape(-1, 1) + np.arange(0, 6)
array([[ 0,  1,  2,  3,  4,  5],
       [10, 11, 12, 13, 14, 15],
       [20, 21, 22, 23, 24, 25],
       [30, 31, 32, 33, 34, 35],
       [40, 41, 42, 43, 44, 45],
       [50, 51, 52, 53, 54, 55]])

Multi-dimensional arrays can also use integer sequences and boolean arrays for access.
1668631278744.jpg
a[(0,1,2,3,4),(1,2,3,4,5)] :
Used to access array indices, which is still a tuple with two elements, where each element is an integer sequence corresponding to the 0th and 1st axes of the array. It retrieves two integers from the corresponding positions of the two sequences to form the index:
a[0,1], a[1,2], ..., a[4,5].
a[3:, [0, 2, 5]] : The 0th axis in the index is a range, selecting all rows after the 3rd row; the 1st axis is an integer sequence, selecting columns 0, 2, and 5.
a[mask, 2] : The 0th axis in the index is a boolean array, selecting rows 0, 2, and 5; the 1st axis is an integer, selecting column 2.

5. Array Structures#

Suppose we need to define a structured array where each element has fields name, age, and weight. In NumPy, it can be defined as follows:

import numpy as np
persontype = np.dtype({
    'names':['name', 'age', 'weight'],
    'formats':['S32','i', 'f']})
a = np.array([("Zhang",32,75.5),("Wang",24,65.2)],
    dtype=persontype)

We first create a dtype object persontype, describing the structure type's various fields through its dictionary parameter. The dictionary has two keywords: names and formats. Each keyword corresponds to a list. names defines the name of each field in the structure, while formats defines the type of each field:
S32 : String type of 32 bytes; since the size of each element in the structure must be fixed, the length of the string needs to be specified.
i : 32-bit integer type, equivalent to np.int32.
f : 32-bit single-precision floating-point type, equivalent to np.float32.
Then we call the array function to create the array, specifying the element type of the created array as the structure persontype through the keyword parameter dtype=persontype. After running the above program, we can execute the following statement in IPython to check the element type of array a:

>>> a.dtype
dtype([('name', '|S32'), ('age', '<i4'), ('weight', '<f4')])

The access method for structured arrays is the same as for general arrays; elements can be obtained through indices. Note that the values of the elements look like tuples, but they are actually structures:

>>> a[0]
('Zhang', 32, 75.5)
>>> a[0].dtype
dtype([('name', '|S32'), ('age', '<i4'), ('weight', '<f4')])

a[0] is a structured element, and it shares memory data with array a, so we can change the corresponding field in the original array by modifying its fields:

>>> c = a[1]
>>> c["name"] = "Li"
>>> a[1]["name"]
"Li"

Structures can be accessed like dictionaries by using string indices to get their corresponding field values:

>>> a[0]["name"]
'Zhang'

We can not only obtain a specific field of a structured element but also directly obtain the fields of a structured array, which returns a view of the original array, so modifying b[0] will change a[0]["age"]:

>>> b=a[:]["age"] # or a["age"]
>>> b
array([32, 24])
>>> b[0] = 40
>>> a[0]["age"]
40

By calling a.tostring or a.tofile methods, we can directly output the binary form of array a:

>>> a.tofile("test.bin")

6. ufunc Operations#

ufunc is short for universal function, which is a function that can operate on each element of an array. Many built-in ufunc functions in NumPy are implemented at the C language level, so their computation speed is very fast. Let's look at an example:

>>> x = np.linspace(0, 2*np.pi, 10)
# Perform sine calculations on each element in array x, returning a new array of the same size
>>> y = np.sin(x)
>>> y
array([  0.00000000e+00,   6.42787610e-01,   9.84807753e-01,
         8.66025404e-01,   3.42020143e-01,  -3.42020143e-01,
        -8.66025404e-01,  -9.84807753e-01,  -6.42787610e-01,
        -2.44921271e-16])

First, we use linspace to generate 10 evenly spaced numbers from 0 to 2*PI, then pass them to the sin function. Since np.sin is a ufunc function, it computes the sine value for each element in x, then returns the result and assigns it to y. After the computation, the values in x remain unchanged; a new array is created to store the results. If we want the results of the sin function to directly overwrite the array x, we can pass the array to be overwritten as the second parameter to the ufunc function. For example:

>>> t = np.sin(x,x)
>>> x
array([  0.00000000e+00,   6.42787610e-01,   9.84807753e-01,
         8.66025404e-01,   3.42020143e-01,  -3.42020143e-01,
        -8.66025404e-01,  -9.84807753e-01,  -6.42787610e-01,
        -2.44921271e-16])
>>> id(t) == id(x)
True

The second parameter of the sin function is also x, so what it does is compute the sine value for each value in x and save the results in the corresponding positions in x. At this point, the return value of the function is still the entire computation result, but it is x, so the ids of the two variables are the same (variables t and x point to the same memory area).
The following small program compares the computation speed of numpy.math and Python's standard library math.sin:

import time
import math
import numpy as np

x = [i * 0.001 for i in xrange(1000000)]
start = time.clock()
for i, t in enumerate(x):
    x[i] = math.sin(t)
print "math.sin:", time.clock() - start

x = [i * 0.001 for i in xrange(1000000)]
x = np.array(x)
start = time.clock()
np.sin(x,x)
print "numpy.sin:", time.clock() - start

# Output
# math.sin: 1.15426932753
# numpy.sin: 0.0882399858083

Calculating the sine value 1 million times, numpy.sin is more than 10 times faster than math.sin. This is due to the loop calculations of numpy.sin being at the C language level. numpy.sin also supports computing the sine of a single value, such as: numpy.sin(0.5). However, it is worth noting that for single number calculations, math.sin is much faster than numpy.sin. Let's look at the following test program:

x = [i * 0.001 for i in xrange(1000000)]
start = time.clock()
for i, t in enumerate(x):
    x[i] = np.sin(t)
print "numpy.sin loop:", time.clock() - start

# Output
# numpy.sin loop: 5.72166965355

Please note that the computation speed of numpy.sin is only 1/5 that of math.sin. This is because numpy.sin is designed to support both array and single value calculations, making its internal implementation in C more complex than that of math.sin. If we also perform loops at the Python level, the differences will become apparent. Additionally, the types of numbers returned by numpy.sin and math.sin differ; math.sin returns Python's standard float type, while numpy.sin returns a numpy.float64 type:

>>> type(math.sin(0.5))
<type 'float'>
>>> type(np.sin(0.5))
<type 'numpy.float64'>

Through the above example, we understand how to use the math library and numpy library's mathematical functions most efficiently.
NumPy has many ufunc functions that provide various calculations. In addition to single-input functions like sin, there are also many multi-input functions, with the add function being one of the most commonly used examples. Let's look at an example:

>>> a = np.arange(0,4)
>>> a
array([0, 1, 2, 3])
>>> b = np.arange(1,5)
>>> b
array([1, 2, 3, 4])
>>> np.add(a,b)
array([1, 3, 5, 7])
>>> np.add(a,b,a)
array([1, 3, 5, 7])
>>> a
array([1, 3, 5, 7])

The add function returns a new array, where each element is the sum of the corresponding elements of the two parameter arrays. It accepts a third parameter to specify the array where the computation result should be written; if specified, the add function will not produce a new array.
Due to Python's operator overloading feature, adding two arrays can be simply written as a+b, while np.add(a,b,a) can be represented as a+=b. Below is a list of operators for arrays and their corresponding ufunc functions; note that the meaning of the division operator "/" varies depending on whether future.division is activated.
y = x1 + x2: add(x1, x2 [, y])
y = x1 - x2: subtract(x1, x2 [, y])
y = x1 * x2: multiply (x1, x2 [, y])
y = x1 / x2: divide (x1, x2 [, y]), if the elements of both arrays are integers, then integer division is performed.
y = x1 / x2: true divide (x1, x2 [, y]), always returns the exact quotient.
y = x1 // x2: floor divide (x1, x2 [, y]), always rounds down the return value.
y = -x: negative(x [,y])
y = x1**x2: power(x1, x2 [, y])
y = x1 % x2: remainder(x1, x2 [, y]), mod(x1, x2, [, y])
Array objects support these operators, greatly simplifying the writing of expressions. However, be aware that if your expressions are very complex and the arrays to be operated on are very large, it may reduce the efficiency of the program due to the generation of a large number of intermediate results. For example, suppose arrays a, b, and c are used in the expression x=a*b+c, it is equivalent to:

t = a * b
x = t + c
del t

This means that an array t is needed to store the multiplication result, and then the final result array x is produced. We can manually decompose an expression into x = a*b; x += c to reduce one memory allocation.

By combining standard ufunc function calls, various expressions for array calculations can be achieved. However, sometimes these expressions are not easy to write, while functions that compute for each element can be easily implemented in Python. In such cases, we can use the frompyfunc function to convert a function that computes a single element into a ufunc function. This allows us to conveniently use the generated ufunc function for array calculations. Let's look at an example.

We want to describe a triangular wave using a piecewise function, as shown in the figure:
1668631406858.jpg
We can easily write the following function to compute the y-coordinate of a point on the triangular wave based on the figure:

def triangle_wave(x, c, c0, hc):
    x = x - int(x) # The period of the triangular wave is 1, so we only take the decimal part of the x-coordinate for calculation
    if x >= c: r = 0.0
    elif x < c0: r = x / c0 * hc
    else: r = (c-x) / (c-c0) * hc
    return r

Clearly, the triangle_wave function can only compute single values and cannot directly process arrays.
Another method is:

def triangle_func(c, c0, hc):
    def trifunc(x):
        x = x - int(x) # The period of the triangular wave is 1, so we only take the decimal part of the x-coordinate for calculation
        if x >= c: r = 0.0
        elif x < c0: r = x / c0 * hc
        else: r = (c-x) / (c-c0) * hc
        return r

    # Create a ufunc function using trifunc, which can directly compute arrays, but the result obtained through this function
    # is an Object array, requiring type conversion
    return np.frompyfunc(trifunc, 1, 1)

y2 = triangle_func(0.6, 0.4, 1.0)(x)

We wrap the three parameters of the triangular wave in the triangle_func function, defining a function trifunc for computing the triangular wave inside it. The trifunc function will use the parameters of triangle_func during its call. Finally, triangle_func returns the result converted by frompyfunc.

It is worth noting that the function obtained using frompyfunc computes array elements of type object because the frompyfunc function cannot guarantee that the data types returned by Python functions are completely consistent. Therefore, it is necessary to convert it again using y2.astype(np.float64) to convert it into a double-precision floating-point array.

7. Broadcasting#

When we use ufunc functions to compute two arrays, the ufunc function computes the corresponding elements of these two arrays, so it requires both arrays to have the same size (same shape). If the shapes of the two arrays are different, broadcasting will be performed as follows:

  1. All input arrays will align with the longest shape among them, with insufficient parts in the shape being padded with 1 at the front.
  2. The output array's shape is the maximum value of the lengths of the corresponding axes of the input arrays.
  3. If the length of an axis of an input array matches the length of the corresponding axis of the output array or is 1, that array can be used for computation; otherwise, an error will occur.
  4. When the length of an axis of an input array is 1, the first value along that axis will be used for computation.
    Understanding these four rules may be a bit tricky, so let's look at a practical example.

First, create a two-dimensional array a with shape (6,1):

>>> a = np.arange(0, 60, 10).reshape(-1, 1)
>>> a
array([[ 0], [10], [20], [30], [40], [50]])
>>> a.shape
(6, 1)
Next, create a one-dimensional array b with shape (5,):
>>> b = np.arange(0, 5)
>>> b
array([0, 1, 2, 3, 4])
>>> b.shape
(5,)

Calculating the sum of a and b results in an addition table, which is equivalent to calculating the sum of all element groups in a and b, resulting in an array with shape (6,5):

>>> c = a + b
>>> c
array([[ 0,  1,  2,  3,  4],
       [10, 11, 12, 13, 14],
       [20, 21, 22, 23, 24],
       [30, 31, 32, 33, 34],
       [40, 41, 42, 43, 44],
       [50, 51, 52, 53, 54]])
>>> c.shape
(6, 5)

Since the lengths of a and b's shapes (i.e., the ndim attribute) are different, according to rule 1, we need to align b's shape with a's by padding 1 at the front, resulting in (1,5). This is equivalent to performing the following calculation:

>>> b.shape=1,5
>>> b
array([[0, 1, 2, 3, 4]])

Thus, the shapes of the two input arrays for the addition operation are (6,1) and (1,5). According to rule 2, the output array's lengths for each axis are the maximum lengths of the input arrays' axes, so the output array's shape is (6,5).

Since the length of b's 0th axis is 1 and a's 0th axis length is 6, to allow them to be added along the 0th axis, b's length along the 0th axis needs to be expanded to 6, which is equivalent to:

>>> b = b.repeat(6,axis=0)
>>> b
array([[0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4],
       [0, 1, 2, 3, 4]])

Since a's 1st axis length is 1 and b's 1st axis length is 5, to allow them to be added along the 1st axis, a's length along the 1st axis needs to be expanded to 5, which is equivalent to:

>>> a = a.repeat(5, axis=1)
>>> a
array([[ 0,  0,  0,  0,  0],
       [10, 10, 10, 10, 10],
       [20, 20, 20, 20, 20],
       [30, 30, 30, 30, 30],
       [40, 40, 40, 40, 40],
       [50, 50, 50, 50, 50]])

After the above processing, a and b can be added element-wise.

Of course, when numpy performs the a+b operation, it does not actually expand the length-1 axes using the repeat function; doing so would be too wasteful of space.

Since this broadcasting computation is very common, numpy provides a quick way to generate arrays like a and b above: the ogrid object:

>>> x,y = np.ogrid[0:5,0:5]
>>> x
array([[0],
       [1],
       [2],
       [3],
       [4]])
>>> y
array([[0, 1, 2, 3, 4]])

ogrid is a very interesting object; it behaves like a multi-dimensional array, allowing access using slice tuples as indices, returning a set of arrays that can be used for broadcasting calculations. The slice indices have two forms:

  1. Start value: End value: Step size, similar to np.arange(start value, end value, step size).

  2. Start value: End value: Length j, when the third parameter is a complex number, it indicates the length of the returned array, similar to np.linspace(start value, end value, length):

    x, y = np.ogrid[0:1:4j, 0:1:3j]
    x
    array([[ 0. ],
    [ 0.33333333],
    [ 0.66666667],
    [ 1. ]])
    y
    array([[ 0. , 0.5, 1. ]])

Using the return values of ogrid, I can easily calculate the values at each point on the x, y grid surface, or the values at each point on the x, y, z grid volume. Below is a program to plot the 3D surface x * exp(x2 - y2):

import numpy as np
from enthought.mayavi import mlab

x, y = np.ogrid[-2:2:20j, -2:2:20j]
z = x * np.exp( - x**2 - y**2)

pl = mlab.surf(x, y, z, warp_scale="auto")
mlab.axes(xlabel='x', ylabel='y', zlabel='z')
mlab.outline(pl)

This program uses the mayavi mlab library to quickly plot the 3D surface shown in the figure.
1668631471880.jpg

8. Methods of ufunc#

The reduce method is similar to Python's reduce function; it operates on the array along the axis axis, equivalent to inserting the operator into all subarrays or elements along the axis.
For example:

>>> np.add.reduce([1,2,3]) # 1 + 2 + 3
6
>>> np.add.reduce([[1,2,3],[4,5,6]], axis=1) # 1,4 + 2,5 + 3,6
array([ 6, 15])

The accumulate method is similar to the reduce method, except that it returns an array with the same shape as the input array, saving all intermediate computation results:

>>> np.add.accumulate([1,2,3])
array([1, 3, 6])
>>> np.add.accumulate([[1,2,3],[4,5,6]], axis=1)
array([[ 1,  3,  6],
       [ 4,  9, 15]])

The outer method, .outer(a,b) method:

>>> np.multiply.outer([1,2,3,4,5],[2,3,4])
array([[ 2,  3,  4],
       [ 4,  6,  8],
       [ 6,  9, 12],
       [ 8, 12, 16],
       [10, 15, 20]])

It can be seen that the result calculated by the outer method is as follows:

#    2, 3, 4
# 1
# 2
# 3
# 4
# 5

If these two arrays are computed step by step according to the equivalent program, it will be found that the multiplication table is ultimately calculated through broadcasting.

9. Matrix Operations#

The product of matrices can be calculated using the dot function. For two-dimensional arrays, it computes the matrix product; for one-dimensional arrays, it computes their dot product. When one-dimensional arrays need to be treated as column vectors or row vectors for matrix operations, it is recommended to first use the reshape function to convert the one-dimensional array into a two-dimensional array:

>>> a = array([1, 2, 3])
>>> a.reshape((-1,1))
array([[1],
       [2],
       [3]])
>>> a.reshape((1,-1))
array([[1, 2, 3]])

dot: For two one-dimensional arrays, it calculates the sum of the products of the corresponding elements (mathematically known as the inner product); in fact, it is the dot product!

10. File Storage#

NumPy provides various file operation functions to facilitate the storage and retrieval of array contents. The file storage format is divided into two categories: binary and text. The binary format files are further divided into NumPy's specific formatted binary type and unformatted type.

Using the array method function tofile can conveniently write the data in the array to a file in binary format. The output data from tofile has no format, so when reading it back with numpy.fromfile, you need to format the data yourself:

>>> a = np.arange(0,12)
>>> a.shape = 3,4
>>> a
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
>>> a.tofile("a.bin")
>>> b = np.fromfile("a.bin", dtype=np.float) # Read data in as float type
>>> b # The read data is incorrect
array([  2.12199579e-314,   6.36598737e-314,   1.06099790e-313,
         1.48539705e-313,   1.90979621e-313,   2.33419537e-313])
>>> a.dtype # Check a's dtype
dtype('int32')
>>> b = np.fromfile("a.bin", dtype=np.int32) # Read data in as int32 type
>>> b # The data is one-dimensional
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])
>>> b.shape = 3, 4 # Modify b's shape according to a's shape
>>> b # This time it is correct
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])

From the above example, it can be seen that setting the correct dtype and shape during reading is necessary to ensure data consistency. Additionally, the tofile function outputs data in C language format regardless of whether the array's arrangement order is in C language format or Fortran language format.

Moreover, if the fromfile and tofile functions are called with the sep keyword parameter specified, the arrays will be input and output in text format.

The numpy.load and numpy.save functions save data in NumPy's specific binary type, and these two functions automatically handle element types and shape information, making it much easier to read and write arrays. However, the files output by numpy.save are difficult to read in programs written in other languages:

>>> np.save("a.npy", a)
>>> c = np.load( "a.npy" )
>>> c
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])

If you want to save multiple arrays into one file, you can use the numpy.savez function. The first parameter of the savez function is the filename, and the subsequent parameters are the arrays to be saved. You can also use keyword parameters to name the arrays; non-keyword parameters will automatically be named arr_0, arr_1, ... . The savez function outputs a compressed file (with the .npz extension), where each file is an npy file saved by the save function, and the filenames correspond to the array names. The load function automatically recognizes npz files and returns a dictionary-like object, allowing you to access the contents of the arrays using their names as keys:

>>> a = np.array([[1,2,3],[4,5,6]])
>>> b = np.arange(0, 1.0, 0.1)
>>> c = np.sin(b)
>>> np.savez("result.npz", a, b, sin_array = c)
>>> r = np.load("result.npz")
>>> r["arr_0"] # Array a
array([[1, 2, 3],
       [4, 5, 6]])
>>> r["arr_1"] # Array b
array([ 0. ,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9])
>>> r["sin_array"] # Array c
array([ 0.        ,  0.09983342,  0.19866933,  0.29552021,  0.38941834,
        0.47942554,  0.56464247,  0.64421769,  0.71735609,  0.78332691])

If you open the result.npz file with decompression software, you will find three files: arr_0.npy, arr_1.npy, and sin_array.npy, each containing the contents of arrays a, b, and c respectively.

Using numpy.savetxt and numpy.loadtxt, you can read and write 1D and 2D arrays:

>>> a = np.arange(0,12,0.5).reshape(4,-1)
>>> np.savetxt("a.txt", a) # By default, save data in '%.18e' format, separated by spaces
>>> np.loadtxt("a.txt")
array([[  0. ,   0.5,   1. ,   1.5,   2. ,   2.5],
       [  3. ,   3.5,   4. ,   4.5,   5. ,   5.5],
       [  6. ,   6.5,   7. ,   7.5,   8. ,   8.5],
       [  9. ,   9.5,  10. ,  10.5,  11. ,  11.5]])
>>> np.savetxt("a.txt", a, fmt="%d", delimiter=",") # Change to save as integers, separated by commas
>>> np.loadtxt("a.txt",delimiter=",") # When reading, specify comma separation
array([[  0.,   0.,   1.,   1.,   2.,   2.],
       [  3.,   3.,   4.,   4.,   5.,   5.],
       [  6.,   6.,   7.,   7.,   8.,   8.],
       [  9.,   9.,  10.,  10.,  11.,  11.]])

This article is referenced from: http://old.sebug.net/paper/books/scipydoc/numpy_intro.html

Loading...
Ownership of this post data is guaranteed by blockchain and smart contracts to the creator alone.