Введение в научный Python-1. Введение в научный Python
 Скачать 6.2 Mb.
|
|
Пример. Построим график полусферы. В начале этого параграфа мы рассмотрели пример построения графика поверхности, заданной явным уравнением y x f z , Если для этого использовать формулу 1 1 , 1 1 1 2 2 y x y x z , то при некоторых значениях x,y под корнем будут стоять отрицательные числа. Поэтому задавать уравнение поверхности следует так: 1 , 0 1 , 1 2 2 2 2 2 2 y x y x y x z Но тогда придется перебирать всю сетку значений x,y и обнулять подкоренное выражение, если оно отрицательно. Вот как это можно сделать, используя логическое индексирование. import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D x=np.linspace( -1.5,1.5,31) y=np.linspace( -1.5,1.5,31) X,Y=np.meshgrid(x,y) Z=1-X**2-Y**2 Z[Z<0]=0 # обнуление отрицательных значений Z=np.sqrt(Z) fig=plt.figure() ax=Axes3D(fig) ax.plot_surface(X,Y,Z,rstride=1,cstride=1,color='cyan') plt.show() ■ 89 Если числовой массив преобразовать в логический (булев), то результат можно использовать при логическом индексировании. Для этого числовой массив можно преобразовать в логический с помощью функции bool8(массив). Она любое число, отличное от 0, преобразует в True, а нули – в False. >>> D=np.arange(5); print(D) [0 1 2 3 4] >>> C=np.bool8([1,0,-1,0,-3]); print(C) [ True False True False True] или по - другому >>> C=np.array([1,0,-1,0,-3],dtype=bool) Тогда >>> D[C]=55; print(D) [55 1 55 3 55] Или наоборот >>> D[C==False]=55; print(D) [ 0 55 2 55 4] Для получения индексов элементов, отличных от нуля, используется функция nonzero(массив). Она возвращает массив индексов. >>> a=np.array([2,0,-1,0,6,0,3]) >>> b=np.nonzero(a); b (array([0, 2, 4, 6], dtype=int64),) Эту функцию можно использовать для определения индексов элементов массива, удовлетворяющих любому условию >>> a=np.array([2,0, -1,5,6,0,3,4,7,0,3]) >>> c=np.nonzero(a>2); c (array([ 3, 4, 6, 7, 8, 10], dtype=int64),) Здесь мы получили индексы всех элементов, значения которых больше 2. Это работает потому, что в логическом массиве „ a>2 ‟ значение False интерпретируется как 0 и выражение np.nonzero(a>2) возвращает индексы True–элементов логического массива. У массивов есть метод nonzero(), который эквивалентен функции nonzero(). Функция flatnonzero() эквивалентна функции nonzero() для одномерных массивов, а с многомерными она работает как с одномерными, предварительно выравнивая их. >>> B=np.array([[3, -4,0],[-2,0,2],[0,0,-3]]) >>> np.flatnonzero(B) array([0, 1, 3, 5, 8], dtype=int64) Для массива определена операция получения среза „:‟ (двоеточие), которая позволяет выбрать диапазон его элементов. В отличие от списков, она создает ссылку, указывающую на данные внутри массива. Поэтому, изменения, сделанные через нее, видны в исходном массиве. >>> a=np.linspace(0,8,9); print(a) [ 0. 1. 2. 3. 4. 5. 6. 7. 8.] 90 >>> b=a[2:6]; b[0]=123; print(b) [123. 3. 4. 5.] >>> print(a) # изменения в массиве b привели к изменениям в массиве a [ 0. 1. 123. 3. 4. 5. 6. 7. 8.] Следующая команда также только создает новую ссылку на массив >>> c=a[:]; c[2]=-55; print(a) [ 0. 1. -55. 3. 4. 5. 6. 7. 8.] Чтобы создать копию массива следует использовать функцию array(). >>> A=np.array(a); A[2]=77; print(A) [ 0. 1. 77. 3. 4. 5. 6. 7. 8.] >>> print(a) # массив a остался без изменений [ 0. 1. -55. 3. 4. 5. 6. 7. 8.] Кроме того, для создания копии массива имеется метод copy(). >>> B=a.copy();B[2]=999; print(B) [ 0. 1. 999. 3. 4. 5. 6. 7. 8.] >>> print(a) [ 0. 1. -55. 3. 4. 5. 6. 7. 8.] Аналогично создается копия фрагмента массива. >>> b=a[2:6].copy() ; print(b) [-55. 3. 4. 5.] >>> b[0]=99; print(a) [0. 1. -55. 3. 4. 5. 6. 7. 8.] Используя срез можно переписать массив в обратном порядке. >>> b=a[len(a):0:-1]; b[1]=77; print(b) [ 8. 77. 6. 5. 4. 3. 123. 1.] >>> print(a) # изменения в массиве b приводят к изменениям в массиве a [ 0. 1. 123. 3. 4. 5. 6. 77. 8.] Можно использовать диапазон с неединичным шагом. >>> a=np.linspace(0,8,9); b=a[0:9:2]; print(b) [ 0. 2. 4. 6. 8.] Элементам среза можно присвоить скаляр или значения из массива «правильного» размера. >>> a=np.arange(1,9) >>> a[3:7]=10 ; print(a) [ 1 2 3 10 10 10 10 8] >>> a[3:7]= np.arange(-4,0,1);print(a) # срезу присваивается массив [ 1 2 3 -4 -3 -2 -1 8] >>> a=np.linspace(0,8,9,dtype=np.int16); print(a) [0 1 2 3 4 5 6 7 8] >>> a[1:9:3]=np.array([111,222,333]); print(a) # срез с шагом 3 [ 0 111 2 3 222 5 6 333 8] Присваивание срезу массива неправильного размера приводит к ошибке. >>> a[1:9:3]=np.array([111,222,333,444]) Ошибка! Доступ к столбцу или строке массива с помощью среза возвращает массив >>> a=np.array([[1,2,3],[4,5,6],[7,8,9]]) 91 >>> a[:,1] # выводим 2-й столбец array([2, 5, 8]) >>> a[2,:] # выводим 3-ю строку array([7, 8, 9]) Используя срез, можно изменять значения элементов двумерного массива. >>> a=np.array([[1,2,3],[4,5,6],[7,8,9]]) >>> a[0:2,0:2]=0;print(a) [[0 0 3] [0 0 6] [7 8 9]] Однако срез a[0:2][0:2] работает не так, как предыдущий. >>> a=np.array([[1,2,3],[4,5,6],[7,8,9]] >>> a[0:2][0:2]=0;print(a) # осторожно! [[0 0 0] [0 0 0] [7 8 9]] Разберем эту особенность для двумерных массивов. Использование «повторной» индексации a[X][Y] работает как суперпозиция двух операций индексации – вначале выполняется левая a[X]. Она возвращает одномерный массив, если X индекс, или набор одномерных массивов (т.е. двумерный массив), если X является срезом. >>> a=np.array([[1,2,3],[4,5,6],[7,8,9]] >>> a[1:3] array([[4, 5, 6], [7, 8, 9]]) Затем выполняется вторая индексация (a[X])[Y]. Если X являлся индексом, то из одномерного массива a[X] выбирается элемент или диапазон элементов Y. Если X являлся срезом, то a[X] является двумерным массивом и индекс Y выбирает соответствующую строку из этого массива. >>> a[1:3][1] array([7, 8, 9]) Если Y диапазон индексов, то для двумерного массива a[X] диапазон Y снова выбирает диапазон строк, т.е. опять возвращается двумерный массив >>> a[1:3][0:1] array([[4, 5, 6]]) Если X и Y индексы, то в соответствии с описанным алгоритмом возвращается элемент Y X a исходного массива, т.е. в этом случае a[X,Y]=a[X][Y]. Обратите также внимание на следующую особенность присваивания переменных массивам. Присваивание другого имени b=a не создает новый массив, а создает синоним его имени. >>> a=np.arange(0,5,1);print(a) [0 1 2 3 4] Следующие команды присваивают переменной b имя массива. Но этим не создается новый массив, а создается другое имя. Если изменить элементы 92 массива b, то тем самым изменяется и массив a. Изменения, сделанные через массив b, видны и в исходном массиве a. >>> b=a; b[2:4]=0;print(a) [0 1 0 0 4] Мы уже говорили, что для копирования данных одного массива в другой, нужно использовать метод copy. >>> a=arange(0,5,1);print(a) [0 1 2 3 4] >>> b=a.copy() >>> b[2:4]=0;print(b); print(a) [0 1 0 0 4] [0 1 2 3 4] Перебирать элементы массива можно в цикле. >>> A=np.linspace(0,5,6); A array([ 0., 1., 2., 3., 4., 5.]) >>> for elem in A: print(elem**2, end=" ") 0.0 1.0 4.0 9.0 16.0 25.0 А следующий цикл выполняется по строкам двумерного массива. >>> B=np.array([[3,-4,0],[-2,0,2],[0,0,-3]]) >>> for row in B: print(row) [ 3 -4 0] [-2 0 2] [ 0 0 -3] Функции работы с массивами. Модуль numpy содержит элементарные функции (дублирующие функции модуля math), аргументами которых могут быть массивы. В таком случае эти функции применяются к массивам поэлементно. Они называются универсальными функциями (ufunc). >>> type(np.sin) >>> c=np.linspace(0, np.pi/2,4) >>> print(np.sin(c)) [0. 0.5 0.8660254 1.] Методы массивов sum(),prod(),max(),min(),mean(),std() вычисляют сумму и произведение элементов массива, максимальный и минимальный элемент, среднее и среднеквадратичное отклонение. >>> b=np.array([[1,2,3],[4,5,6]]);b [[1 2 3] [4 5 6]] >>> b.sum() 21 >>> b.prod() 720 93 >>> b.max() 6 >>> b.mean() 3.5 >>> b.std() 1.707825127659933 Эти методы умеют вычислять соответствующие значения по каждой из координат массива. >>> a=np.array([[1,2,3],[4,5,6],[7,8,9]]) ; print(a) [[1 2 3] [4 5 6] [7 8 9]] >>> a.sum() # сумма всех элементов 45 >>> print(a.sum(0)) # суммирование по нулевому индексу [12 15 18] Сравните с суммой по столбцам >>> a[:,0].sum(), a[:,1].sum(),a[:,2].sum() (12, 15, 18) Аналогично >>> print(a.sum(1)) # суммирование по первому индексу [ 6 15 24] Сравните с суммой по строкам >>> a[0,:].sum(), a[1,:].sum(),a[2,:].sum() (6, 15, 24) Аналогично работают методы prod(), max(), min() и т.д. >>> a.min(1) # минимум по строкам array([1, 4, 7]) Методы argmin() и argmax() возвращают индексы минимального и максимального элементов по всему массиву или по строкам (столбцам). При определении индекса по всему массиву он вначале выравнивается построчно до одномерного массива. Затем возвращается индекс элемента в таком массиве. >>> a=np.array([[2,-1,5],[7,3,-2],[3,9,5]]); a array([[ 2, -1, 5], [ 7, 3, -2], [ 3, 9, 5]]) >>> a.argmin() # индекс минимального элемента 5 >>> a.argmax() # индекс максимального элемента 7 >>> a.argmin(0) # индексы минимальных по столбцам array([0, 0, 1], dtype=int64) >>> a.argmin(1) # индексы минимальных по строкам array([1, 2, 0], dtype=int64) 94 >>> a.argmax(0) # индексы максимальных по столбцам array([1, 2, 0], dtype=int64) >>> a.argmax(1) # индексы максимальных по строкам array([2, 0, 1], dtype=int64) Накопительной суммой массива n j a a a a ..., , ..., , , 1 0 a называется массив n j c c c c ,..., ..., , , 1 0 c , где j i i j a c 0 . Она вычисляется с помощью функции cumsum(). >>> a=np.array([-1,3,2,-4,5]) >>> np.cumsum(a) array([-1, 2, 4, 0, 5], dtype=int32) Функция numpy.diff(a, n=1,axis=-1) вычисляет конечные разности n–го порядка массива a вдоль заданной оси axis. В простейшем случае функция diff(a) вычисляет конечные разности первого порядка массива a по последней оси. Например, для одномерного массива n i i a 0 a функция возвращает вектор n i i i n i i a a a 1 1 1 d , т.е. вектор, составленный из разностей соседних элементов. >>> np.diff([1,3,8,25]) array([ 2, 5, 17]) >>> np.diff([[1,3,8,25],[3,4,5,6]]) array([[ 2, 5, 17], [ 1, 1, 1]]) >>> np.diff([[1,3,8,25],[3,4,5,6],[1,2,3,4]],axis=0) array([[ 2, 1, -3, -19], [ -2, -2, -2, -2]]) Функции round(), around() (и метод round()) округляют элементы массивов до ближайших значений, имеющих заданное количество десятичных знаков. >>> x=np.array([ -3.546, 2.215, -5.149, 4.375, 6.85, -8.4237]) >>> b=x.round(2); b array([-3.55, 2.22, -5.15, 4.38, 6.85, -8.42]) Имеется возможность вычисления кусочных выражений вида n n n n n x x x f x x x x f x x x x f x x x f x f , , , , 1 1 2 1 2 1 1 Знаки условных операций, содержащие равенства (меньше или равно, больше или равно) могут стоять в любой стороне условий. Для вычисления кусочных выражений в модуле numpy определена функция piecewise(), которая имеет следующий формат: 95 piecewise(x, список_условий, список_функций [,*args,**kw]) Здесь в качестве первого аргумента x должен стоять массив и результатом будет массив того же размера. Вот как можно вычислять функцию абсолютного значения, которая равняется x , если 0 x , и x , если 0 x >>> x=np.array([ -3.5, 2.2, -5.1, 4.3, 6.8, -8.4]) >>> np.piecewise(x, [x < 0, x >= 0], [lambda x: -x, lambda x: x]) array([ 3.5, 2.2, 5.1, 4.3, 6.8, 8.4]) Если участков «аналитичности» больше двух, то следует использовать операции логической алгебры для массивов: logical_and(), logical_or(), logical_xor(), logical_not(). >>> np.piecewise(x,[x< -3,np.logical_and(x>= -3,x<6),x>=6],[-1,1,2]) array([-1., 1., -1., 1., 2., -1.]) >>> X=np.piecewise(x,[x< -3, np.logical_and(x>= -3,x<4), np.logical_and(x>=4,x<=6), x>=6], [lambda t: -t, lambda t: t**2, lambda t: t+100, lambda t: t**0.5]) >>> X array([3.5, 4.84, 5.1, 104.3, 2.60768096, 8.4]) Имеется возможность преобразовывать Python функции для работы с аргументами массивами. Это называется векторизацией. Векторизация позволяет преобразовать скалярную функцию так, чтобы в качестве аргументов она могла принимать последовательности элементов (списки, кортежи и т.д.) или массивы, а не только скаляры. После векторизации результатом работы функции будет последовательность элементов или массив, полученные поочередным применением скалярной функции к каждому элементу аргумента. Векторизация используется для удобства, а не для ускорения вычислений. Для векторизации используется функция numpy.vectorize(pyfunc[,...]) . Ее первым аргументом является скалярная Python функция. Вот несколько примеров. >>> from math import sin >>> import numpy as np Векторизация библиотечной скалярной функции >>> mysin = np.vectorize(sin) # векторизация функции math.sin >>> x=[0,np.pi/6,np.pi/4,np.pi/2,np.pi/2] >>> mysin(x) # аргумент - список array([ 0. ,0.5 , 0.70710678, 1. , 1. ]) Векторизация Python функции пользователя. >>> def f(x): return sin(x)**2+x 96 >>> fvec = np.vectorize(f) # векторизация функции пользователя >>> fvec(x) # аргумент - список array([ 0. , 0.77359878, 1.28539816, 2.57079633, 2.57079633]) Векторизация лямбда функции одной переменной. >>> g=lambda x: x**2+x+1 >>> m=np.array([[2,5],[1,8]]) >>> gvec=np.vectorize(g) # векторизация лямбда функции >>> gvec(m) # аргумент - двумерный массив array([[ 7, 31], [ 3, 73]]) Векторизация лямбда функции двух переменных. >>> x=np.linspace(0,4,5); x array([ 0., 1., 2., 3., 4.]) >>> y=np.linspace(-2,2,5); y array([-2.,-1., 0., 1., 2.]) >>> h=lambda x,y: x**2+y**2 # лямбда функция двух аргументов >>> hvec=np.vectorize(h) >>> hvec(x,y) # аргументы – два одномерных массива array([ 4., 2., 4., 10., 20.]) Имеется большое количество других функций и методов манипулирования массивами. Кроме того, мы не описали необязательные аргументы многих функций. С ними вы можете познакомиться по справочной системе пакета numpy . Для этого следует выполнить команду help(np.имя_функции), если np – это имя модуля numpy. Подведем небольшой итог. Массивы похожи своим поведением на списки, но имеют следующие особенности: элементы массива всегда представляются объектами одного типа, например, все элементы являются вещественными числами; в момент создания массива количество его элементов должно быть известно; в порядке исключения оно может быть изменено, например, функцией resize(); объекты массивов импортируются из пакета NumPy и не являются стандартной частью Python; над массивами можно выполнять «векторные» и «матричные» операции с помощью функций пакета NumPy. 3.2 Элементы линейной алгебры Функции подмодулей numpy.linalg и scipy.linalg реализуют основные операции линейной алгебры. Они в значительной степени дублирую друг друга. Однако лучше использовать модуль scipy.linalg в силу большей эффективности его процедур. Кроме того, имеется модуль numpy.matlib, содержащий функции, которые возвращают/создают матрицы, а не массивы (в 97 NumPy тип «матрица» не совпадает с типом «массив»). В этом модуле для объектов типа «матрица» основные матричные операции реализованы в виде операций, а не функций. Здесь мы опишем только наиболее важные функции этих модулей. При этом для краткости не будем упоминать о необязательных аргументах многих функций. В тех случаях, когда они есть, но мы о них не упоминаем, в описании будут использоваться квадратные скобки, означающие необязательность, и троеточие, обозначающее пропуск в описании одного или нескольких аргументов. Например, norm(a[,ord,...]) означает, что a является обязательным аргументом, кроме того может использоваться необязательный аргумент ord и могут быть другие необязательные аргументы. Возвращаемые функциями величины обычно понятны по смыслу, однако более подробную информацию о них вам следует получить из справочной системы. Функция det(a[,...]) вычисляет определитель квадратной матрицы. >>>import numpy as np >>>import scipy.linalg as la >>>a = np.array([[1, 2],[3, 4]]) >>>la.det(a) # вычисление определителя квадратного массива -2.0 Функция norm(a[,ord,...]) вычисляет норму массива (вектора или матрицы). >>> la.norm(a) # норма Евклида 5.4772255750516612 >>> np.sqrt(1**2+2**2+3**2+4**2) 5.4772255750516612 Второй аргумент определяет вид нормы (Евклида, Фробениуса, Минковского и т.д.). Если аргумент ord опущен, то вычисляется норма Евклида. Функция inv(a[,...]) вычисляет обратную матрицу/массив. >>> a=np.array([[1,2],[3,4]]) >>> a1=la.inv(a); a1 array([[-2. , 1. ], [ 1.5, -0.5]]) Проверим, что a*a1=I. Имеем (напомним, что функция matmul появилась в пакете numpy начиная с версии 1.10) >>> np.matmul(a,a1) # или np.dot(a,a1) array([[ 1.00000000e+00, 0.00000000e+00], [ 8.88178420e-16, 1.00000000e+00]]) Здесь недиагональные элементы практически равны нулю. Неточность возникает из–за погрешности вычислений. Если в последнем матричном умножении выполнить округление с точностью до 10 цифр после десятичной точки, то получим >>> np.round_(np.dot(a,a1),10) array([[ 1., 0.], [ 0., 1.]]) |