пятница, 12 мая 2017 г.

Математика в Python

Некоторые библиотеки python используемые в научных вычислениях, визуализации и обработке данных:
  • numpy - работа с матрицами и многомерными массивами
  • scipy - научные и инженерные вычисления, использует numpy;
  • matplotlib - построение графиков и диаграмм;
  • seaborn - визуализация статистических данных, эстетичнее чем matplotlib; 
  • mpld3 - использование D3.js для построения интерактивных matplotlib графиков в окне браузера;
  • pandas - анализ данных: статистики, регрессия, визуализация и т.п.

Для вычислений на лету, экспериментов и вообще работы с научными данными рекомендуется использовать Jupiter Notebook - интерпретатор Python (и не только) прямо в браузере. Интерфейс похож на оный математических пакетов. Попробовать прямо сейчас.

Линейная алгебра

import numpy as np 

Создание матрицы

A = np.array([ [1,2,3],
               [3,2,1],
               [2,1,3] ] ) 
Данная матрица это двумерный массив (список), только в обёртке numpy.

Операции.

Умножение на скаляр
A * 3.14
> array([[ 3.14,  6.28,  9.42],
         [ 9.42,  6.28,  3.14],
         [ 6.28,  3.14,  9.42]])
 
Умножение на вектор
A * [2,3,4]
Причём здесь вектор не обязательно должен быть типом numpy.

Сложение матриц 
A+B

Умножение матриц
A @ B
или
np.dot(A,B)

Обратная матрица
np.linalg.inv(A)

Решение СЛАУ
import numpy
numpy.set_printoptions(precision=4, suppress=True)
A = numpy.matrix([
    [1, 1, 3],
    [1, 2, 2],
    [1, 3, 1]])
B = numpy.matrix([[1], [2], [3]])
X = numpy.matrix([[0], [0], [0]])
X = numpy.linalg.solve(A, B)
 
Стоит обратить внимание на то, что вектор-столбец определяется как матрица из одного столбца, а не как список.

Интерполяция

from scipy.interpolate import interp1d
 
X = [1,2,3,4,5,6,7]
Y = [1,4,9,16,25,36,49]
func = interp1d(X, Y, kind='cubic')
func(5.3)
 
>>> array(28.090000000000003)
 
 
Минимизация функции
from scipy.optimize import minimize
def funcXX(x): 
   return (x-3)*(x-3) - 5
x0 = 100 # начальное значение
minimize(funcXX, x0)
 
Вывод:
nfev: 12
jac: array([ 0.])
message: 'Optimization terminated successfully.'
fun: -5.0
success: True
x: array([ 3.])
njev: 4
hess_inv: array([[ 0.5]])
status: 0
Минимум функции: f(3) = -5


Для функции с дополнительными (не оптимизируемыми) параметрами:
def funcXX(x, c1=0, c2=0, c3=0): 
   return (x-c1)*(x-c2) - c3
x0 = 100
minimize(funcXX,100, args=(3,3,5))

Оптимизация по нескольким параметрам (по списку):
def funcXX(x): 
    return ( 1 - x[0] )**2 + 100 * ( x[1] - x[0]**2 )**2
X0 = [10,10]
minimize(funcXX, X0)

Вывод:
     nfev: 292
      jac: array([  1.38046018e-06,  -6.34189745e-07])
  message: 'Optimization terminated successfully.'
      fun: 1.956291539083183e-11
  success: True
        x: array([ 0.99999558,  0.99999115])
     njev: 73
 hess_inv: array([[ 0.50636053,  1.01234176],
       [ 1.01234176,  2.02894044]])
   status: 0
Параметром method указывается конкретный метод оптимизации из числа возможных:
'Nelder-Mead', 'Powell', 'CG', 'BFGS', 'Newton-CG', 'Anneal', 'L-BFGS-B', 'TNC',  'COBYLA',  'SLSQP', 'dogleg', 'trust-ncg'. По-умолчанию применяется BFGS.

Чиленное интегрирование и дифферецирование [Integration (scipy.integrate)]

from scipy.integrate import quad 
# интеграл от sin(x) в пределах от 0 до pi/2
quad(lambda x: sin(x): 0, pi/2) 
>>> (0.9999999999999999, 1.1102230246251564e-14)

Результат выполнения функции - кортеж, первый элемент которого - значение интеграла.

Значение производной в точке 2 (dx=1e-6 - шаг для вычисления производной численным методом):
from scipy.misc import derivative 
derivative(lambda t: -3/(t+2), 2, dx=1e-6) 
 
0.18750000002620837

График функции [подробнее, ещё
Чтобы  

Путём задания функции:
from pylab import *
t = arange(0.0, 2.0, 0.01)
s = sin(2*pi*t)
plot(t, s)
xlabel('time (s)')
ylabel('voltage (mV)')
# добавить горизонтальную штриховую линию красного цвета
plt.axhline(y=-sh.S,ls ="dashed",color="red")
# можно указать размер шрифта в заголовке
title('About as simple as it gets, folks', fontsize=16)
grid(True)
savefig("test.png")
show()

График по точкам:
import matplotlib.pyplot as plt
plt.plot([1, 2, 3, 4], [1, 4, 9, 16], 'ro-', linewidth=2.0)
# включить легенду и подобрать для неё оптимальное пложение
pl.legend(loc='best')
plt.show()
 
'ro-' - красные круглые маркеры, соединенные сплошной линией
***
Закрасить область между прямой y = 0.9 и кривой y = Ffunc(); между абсциссами x1 и x2.
fill_between(np.linspace(x1,x2,100), 0.9, Ffunc(np.linspace(x1,x2,100)), color = 'yellow', alpha=0.3)
 
Кривая представляется массивом значений y. Интервал абсцисс задаётся массивом той же размерности.

***
Добавить подпись (греческие и проч. символы) в LaTeX нотации:
plt.ylabel("$\\sigma_r$ [MPa]")

***
Нарисовать наклонные подписи под осью [SO]
plt.xticks(rotation=30)

***
Text rendering With LaTeX
http://matplotlib.org/users/mathtext.html



dx, dy = 0.25, 0.25
# generate 2 2d grids for the x & y bounds
y, x = np.mgrid[slice(-30, 30 + dy, dy), slice(-30, 30 + dx, dx)]
z = (1 - x / 2. + x ** 5 + y ** 3) * np.exp(-x ** 2 - y ** 2)
z = z[:-1, :-1]
z_min, z_max = -np.abs(z).max(), np.abs(z).max()
plt.clf()
plt.pcolor(x, y, z, cmap=cm.gray, vmin=z_min, vmax=z_max)
plt.title('F(u0,ssr0)')
# set the limits of the plot to the limits of the data
plt.axis([x.min(), x.max(), y.min(), y.max()])
plt.colorbar()
plt.show() 



Теплокарта на основе набора значений X,Y [SO]
Чем больше раз встречается то или иное значение тем "горячее" будет точка.
import numpy as np
import numpy.random
import matplotlib.pyplot as plt
 
# Generate some test data
x = np.random.randn(8873)
y = np.random.randn(8873)
 
heatmap, xedges, yedges = np.histogram2d(x, y, bins=(512, 512))
extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]
 
plt.clf()
plt.imshow(heatmap, extent=extent)
plt.show()


***
Сделать картинки красивыми поможет пакет seaborn. Эту надстройку над matplotlib можно использовать непосредственно. Однако, пакет удобен прежде всего для визуализации различных данных, а не исключительно математических графиков. Поэтому при необходимости используются только стиль, а остальные операции выполняются с помощью непосредственно matplotlib.
import seaborn as sb
sb.set_style("darkgrid")  # установить стиль seaborn
Методом set можно изменить другие параметры рисования, например ращмер шрифта
sb.set(font_scale=1.55)

***

Создавать одновременно красивые и интерактивные диаграммы и графики можно с помощью пакета plotly. Удобнее всего работать с этим пакетом, если оболочка Питона запущена прямо в браузере ибо созданные диаграммы представляют собой автоматически сгенерированные html страницы с включением кода на JS. 

Графы

Библиотеки: NetworkX,

Работа с массивами данных

Сохранить массивы в cvs файл
import numpy
a = numpy.asarray([ [1,2,3], [4,5,6], [7,8,9] ])
numpy.savetxt("foo.csv",  numpy.transpose(a), delimiter=" ", fmt="%1.4f", header = "val1, val2, val3")
foo.cvs:
# val1, val2, val3
1.0000 4.0000 7.0000
2.0000 5.0000 8.0000
3.0000 6.0000 9.0000
Загрузка данных из cvs файла.
numpy.genfromtxt
from numpy import genfromtxt
my_data = genfromtxt('my_file.csv', delimiter=',', unpack=True).transpose()
Функция genfromtxt возвращет данные как они записаны в вайле (см. выше). Чтобы получить вместо списка данных по строками (одна строка - один массив), список данных по столбцам, нужно перевернуть (.transpose()) или передать в функцию параметр unpack=True.
>>> my data
 
array([[ 1.,  2.,  3.],
           [ 4.,  5.,  6.],
           [ 7.,  8.,  9.]])
При необходимости можно сразу записать данных в отдельные массивы. Их число должно совпадать с числом полей.
v1, v2, v3 = genfromtxt('data/test2', unpack=True) 
Функция numpy.genfromtxt выполняет ту же работу, но с дополнительными возможностами вроде автоматического заполнения пустых полей.

Алсо

Pandas - анализ данных

Ссылки

1 комментарий:

  1. Отличная статья! Долго искал информацию по созданию матриц и вычислений над ними и только здесь вся информация понятна и корректна. Спасибо.

    ОтветитьУдалить