NumPy Ufunc operation

The method of Ufunc operation broadcast Ufunc References

Ufunc Operations

Ufunc is the abbreviation for universal function, which is an operation of each element of an array. Many of the Ufunc functions built into the NumPy are implemented at the C language level, so they are calculated very quickly. Let's take a look at an example:

In [$]: x = np.linspace (0, 2*np.pi,) in

[%]: y = np.sin (x) in

[si]: y
out[54]: 
Array ([  0.00000000 e+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.44929360e-16])

First we use Linspace to produce a 10 number of equidistant from 0 to 2*pi, and then pass it to the SIN function, because np.sin is a ufunc function, so it asks for the sine of each element in X, then returns the result and assigns it to Y. After the calculation, the value in X does not change, but a new array is created to save the result. If we want to overwrite the result of the SIN function directly on the array x, we can pass the overridden array as the second argument to the Ufunc function. For example::

In [$]: T = np.sin (x,x)

in [+]: x
out[56]: 
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.44929360e-16] in

[()]: id (t) = = ID (x)
out[57]: True

The second parameter of the SIN function is also x, so what it does is to get the sine of each value in X and save the result to the corresponding position in X. The return value of the function is still the result of the entire calculation, except that it is X, so the IDs of the two variables are the same (the variable t and the variable x point to the same memory region).

Use%timeit to compare the calculation speed of Numpy.sin and Math.sin in Python3

In [to]: import NumPy as NP in                      
[[): Import NumPy as NP in

[+]: Import math in [MB

]: x = [i * 0.001 for I In range (1000000)] in

[?]:%timeit [Math.sin (i) for i-X]
183 ms±4.2 ms per loop (mean±std. Dev. 7 run  S, 1 loop each) in

[$]: x = Np.array (x) in

[$]:%timeit np.sin (x)
39.6 ms±1.14 ms per loop (MEAN±STD. Dev. of 7 runs, loops each)

This is the result of the operation on an Ubuntu virtual machine, you can see that the average operation speed is 4 times times faster numpy.sin than Math.sin, and the speed variance is much smaller than that of Math.sin.

Calculate 1 million times sine on my computer, Numpy.sin is 10 times times faster than Math.sin. This is conducive to numpy.sin in the C language level of the loop calculation. Numpy.sin also supports the sine of a single value, for example: Numpy.sin (0.5). It is worth noting, however, that the calculation of a single number is much faster than the Numpy.sin, let's look at the following test program: Math.sin

In [to]: x = [i * 0.001 for I in range (1000000)] in

[-]:%timeit [Np.sin (i) for i-X]
929 ms±28.6 ms per loop (mean±std. Dev. of 7 runs, 1 loop each)

Please note that the Numpy.sin calculation speed is only 1/5 of Math.sin. This is because Numpy.sin's internal implementation of the C language is much more complex than math.sin in order to support both array and single value calculations, and we can see the difference if we also loop at the Python level. In addition, the type of number returned by Numpy.sin differs from the type returned by Math.sin, Math.sin returns the standard float type of Python, and Numpy.sin returns a Numpy.float64 type:

In [to]: type (Math.sin (0.5))
out[70]: float in

[]: Type (np.sin (0.5))
out[71]: Numpy.float64

There are many ufunc functions in numpy that provide us with a wide range of computations. In addition to sin's single input function, there are many more input functions, the Add function is one of the most common examples. Let's take a look at an example:

in [[to]: a = Np.arange (0, 4) in
[a]: a = Np.arange (0, 4) in [an]

: a
out[73]: Array ([0, 1, 2, 3]) in

[74]: b = Np.arange (1, 5) in

[$]: b
out[75]: Array ([1, 2, 3, 4]) in

[[]: Np.add (A, B, a)
out[76]: Array ([1, 3 , 5, 7]) in

[+]: a
out[77]: Array ([1, 3, 5, 7])

The Add function returns a new array in which each element of the array is the sum of the corresponding elements of the two-parameter array. It accepts the 3rd parameter to specify the array to write the result to, and if so, the Add function no longer produces a new array.

Because of Python's operator overloading capabilities, the addition of two arrays can be simply written as A+b, while Np.add (A,b,a) can be represented by a+=b. The following is a list of the operators of the array and their corresponding ufunc functions, noting that the meaning of the division "/" differs depending on whether the __future__.division is activated.

Operation	corresponding function
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 two arrays are integers, then divide by integer
y = x1/x2	True Divide (x1, x2 [, y]), always returns the exact quotient
y = x1//x2	Floor Divide (x1, x2 [, y]), always rounding 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, which greatly simplifies the writing of formulas, but note that if your formula is complex and the array you want to compute is large, you can reduce the efficiency of the program by generating a large number of intermediate results. For example, suppose a B C three array is computed using the formula X=a*b+c, then it is equivalent to:

t = A * b
x = t + C
del t

That is, you need to generate an array T to hold the result of the multiplication, and then produce the final result array x. We can decompose an equation into x = A*b by hand; x + = c to reduce the memory allocation one time.

By combining the call of standard Ufunc function, we can realize the array calculation of various formulas. But sometimes this formula is not easy to write, and for each element of the calculation function is easy to implement in Python, then you can use the Frompyfunc function to calculate a single element of the function to convert to the Ufunc function. This makes it easy to use the Ufunc function to calculate the array. Let's take a look at an example.

We want to use a piecewise function to describe the triangle wave, which looks like this in Figure 2.4:

Write the following function to compute the y-coordinate of a triangle wave, as shown in the above figure:

def triangle_wave (x, C, C0, HC):
    x = X-int (x) # The trigonometric wave has a period of 1, so only the decimal part of the x-coordinate is calculated to compute
    if x >= c:r = 0.0
    elif x < C0:R = x/c0 * HC
    else:r = (c-x)/(C-C0) * HC return
    R

Obviously the Triangle_wave function can only compute a single numeric value and cannot be processed directly by the array. We can use the following method to first work out a list and then use the array function to convert the list to an array: comprehension

x = Np.linspace (0, 2, 1000)
y = Np.array ([Triangle_wave (T, 0.6, 0.4, 1.0) for T in X])

This approach requires calling functions with the list containment syntax each time, which can be cumbersome for multidimensional arrays. Let's look at how to use the Frompyfunc function to solve this problem:

Triangle_ufunc = Np.frompyfunc (lambda x:triangle_wave (x, 0.6, 0.4, 1.0), 1, 1)
y2 = triangle_ufunc (x)

The invocation format for Frompyfunc is Frompyfunc (func, Nin, Nout), where Func is the function that computes a single element, Nin is the number of input parameters for this function, and Nout is the number of return values for this function. Although the Triangle_wave function has 4 parameters, because of the latter three C, C0, HC in the entire calculation of the value is fixed, so the resulting ufunc function is actually only one parameter. To satisfy this condition, we use a lambda function to wrap the parameters of the Triangle_wave one at a time. So the function that passes in the Frompyfunc has only one argument. In this way, the efficiency is not too high, there is another way:

def triangle_func (c, C0, HC):
    def trifunc (x):
        x = X-int (x) # triangle Wave has a period of 1, so only the decimal part of the x-coordinate is calculated
        if x >= c:r = 0.0
        Elif x < C0:R = X/C0 * HC
        else:r = (c-x)/(C-C0) * HC return
        R

    # Creates a Trifunc function with the Ufunc function, which can be counted directly on the array Calculate, but through this function
    # computed is an object array, which requires type conversion return
    Np.frompyfunc (Trifunc, 1, 1)

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

We use the function Triangle_func to package the three parameters of the triangular wave, and define a function of the trigonometric wave in its interior trifunc,trifunc function is evaluated using TRIANGLE_FUNC parameters when called. Finally Triangle_func returns the result with the Frompyfunc conversion.

It is noteworthy that the type of the array element computed using the Frompyfunc function is object, because the Frompyfunc function does not guarantee that the data types returned by the Python function are exactly the same. You also need to Y2.astype (Np.float64) again to convert it to a double-precision floating-point number group. Broadcast

When we use the UFUNC function to compute two arrays, the Ufunc function computes the corresponding elements of the two arrays, so it requires that the two arrays have the same size (shape). If the shape of the two arrays is different, the following broadcast (broadcasting) processing is performed: Let all input arrays be aligned to the longest array of shape, and the insufficient parts of the shape are padded by the front plus 1. The shape of the output array is the maximum value on each axis of the input array shape if an axis of the input array is the same length as the corresponding axis of the output array, or if it is 1 o'clock, the array can be used to compute, otherwise an error

When the length of an axis of an input array is 1 o'clock, the first set of values on this axis is used along this axis

The above 4 rules may be more difficult to understand, let's look at a practical example.

First create a two-dimensional array of a, whose shape is (6,1):

In [m]: a = Np.arange (0, a). Reshape ( -1, 1)

in [m]: a
out[79]: 
Array ([[0],
       [ten], [M
       ],
       [ ], [+]
       ,
       [6]]) in

[[]: A.shape
out[80]: (1)

Then create a one-dimensional array of B, whose shape is (5,):

In [i]: b = Np.arange (0, 5) in

[to]: B
out[82]: Array ([0, 1, 2, 3, 4]) in

[]: B.shape
out[83]: (5,)

Calculate the sum of A and B, and get an addition table, which is equivalent to the sum of all the elements in the a,b and a shape (6,5) array:

In [$]: c = a + b in

[[]: C
out[85]: 
Array ([[0,  1,  2,  3,  4],
       [A, A, 14],
  [20, [
       50, 51, 52, 53, 54]], [K, M, M], [K, K, K, O];]

Because the shape length of a and B (that is, the Ndim property) is different, according to Rule 1, the shape of B needs to be aligned to a, so the shape of B is preceded by a 1, which is padded (1,5). Equivalent to doing the following calculations:

in [[to]: b.shape=1,5 in

[to]: B
out[87]: Array ([[[0, 1, 2, 3, 4]]) in

[]: B.ndim
out[88]: 2

The shape of the two input arrays of the addition operation is (6,1) and (1,5), according to rule 2, the length of each axis of the output array is the maximum value of the length of each axis of the input array, and the shape of the output array is (6,5).

Because the length on the No. 0 axis of B is 1, and the length on the No. 0 axis of a is 6, so in order for them to add on the No. 0 axis, the length of B on the No. 0 axis needs to be extended to 6, which is equivalent to:

In [m]: b = b.repeat (6,axis=0)

in [m]: b
out[90]: 
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]]

Because the 1th axis of a is 1 and the first axis length is 5, so in order for them to add on the 1th axis, you need to extend the length of a on the 1th axis to 5, which is equivalent to:

In [$]: a = a.repeat (5, Axis=1)

in [$]: a
out[92]: 
Array ([[0,  0,  0,  0,  0],
       [ [
       50,
       50, 50, 50, [A, A, M], [A, A, a], [a], [A, A, a], [a]. , 50]])

After the above processing, a and B can be added to the corresponding elements of the operation.

Of course, NumPy does not actually extend an axis with a length of 1 to a repeat function when performing a a+b operation, which would be a waste of space.

Because this kind of broadcast computing is very common, NumPy provides a quick way to produce an array such as the above a,b: Ogrid object:

In [$]: X,y = np.ogrid[0:5,0:5]

in [m]: x
out[94]: 
Array ([[0],
       [1], [
       2], [
       3],
       [4]]) In

[to]: Y
OUT[95]: Array ([[0, 1, 2, 3, 4]])

in [%]: X.ndim, Y.ndim
out[96]: (2, 2)

Ogrid is an interesting object that, like a multidimensional array, is accessed using a slice element as a subscript, returning a set of arrays that can be used to broadcast calculations. Its slices are labeled in two ways:

Start value: End value: Step, and Np.arange (start value, end value, step size) Similar

Start value: End value: Length J, when the third argument is an imaginary number, it represents the length of the returned array, similar to the 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.]])

Ogrid Why not function

According to Python syntax, slicing syntax separated by colons can only be used in brackets, and if Ogrid is a function, then these slices must be created using the slice function, which obviously increases the length of the code.

With the return value of Ogrid, I can easily compute the values of the dots on the X, y, and z grids. The following is a program to draw a three-dimensional surface X * EXP (x**2-y**2):

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 draw the 3D surface as shown in Figure 2.5, and the relevant content about Mlab will be covered in future chapters.

the method of Ufunc

The Ufunc function itself has some methods that are valid only for the UFUNC function of two input one output, and the other Ufunc objects will throw valueerror exceptions when they are called by these methods.

The reduce method, like the python reduce function, operates on an axis axis to manipulate the array, which is equivalent to inserting an operator into all of the array or elements along the axis axis.

<op>.reduce (array=, axis=0, Dtype=none)

For example:

in [[]: Np.add.reduce ([1,2,3]) # 1 + 2 + 3
out[97]: 6 in [[]

: Np.add.reduce ([[[1,2,3],[4,5,6]], Axis=1) # 1,4 + 2  , 5 + 3,6
out[98]: Array ([6,]) in

[$]: Np.add.accumulate ([1,2,3])
out[99]: Array ([1, 3, 6]) in

[100]: Np.add.accumulate ([[[1,2,3],[4,5,6]], Axis=1)
out[100]: 
Array ([[1,  3,  6],
       [4,  9, 15]])

The Reduceat method calculates the results of multiple groups of reduce, specifying the start and end positions of a series of reduce with the indices parameter. Reduceat's calculations are somewhat special, let's explain this by an example:

In [$]: a = Np.array ([1,2,3,4]) in

[/]: result = Np.add.reduceat (a,indices=[0,1,0,2,0,3,0]) in

[105]: Result
out[105]: Array ([1,  2,  3,  3,  6,  4, 10])

The reduce function is called to calculate a value for each element in indices, so the end result is the same length as the indices. The result array, with the exception of the last element, is calculated as follows:

If indices[i] < indices[i+1]:
    result[i] = Np.reduce (a[indices[i]:indices[i+1])
else:
    result[i] = a [Indices[i]]

And the last element is calculated as follows:

Np.reduce (A[indices[-1]:])

So in the example above, each element of the result is calculated as follows:

1:a[0] = 1
2:a[1] = 2
3:a[0] + a[1] = 1 + 2
3:a[2] = 3
6:a[0] + a[1] + a[2] =  1 + 2 + 3 = 6< C5/>4:A[3] = 4
10:a[0] + a[1] + a[2] + a[4] = 1+2+3+4 = 10

Outer method, the calculation of the <op>.outer (A,b) method is equivalent to the following program:

>>> A.shape + + (1,) *b.ndim
>>> <op> (a,b)
>>> a = A.squeeze ()

The function of the squeeze is to eliminate the axis of length 1 in array A. If you don't understand this equivalent program, let's look at an example:

>>> Np.multiply.outer ([1,2,3,4,5],[2,3,4])
Array ([[2,  3,  4],
       [4,  6,  8],
       [6,  9,],
       [8, a],
       [10, 15, 20]]

You can see that the results computed by the outer method are the following multiplication tables:

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

If the two arrays are computed step-by-step by the equivalent program, the multiplication table will eventually be calculated by broadcast. References

Doing scientific calculations with Python