Matematikai Algoritmusok és Felfedezések I.

13. Előadás: Scipy

2021 május 10.

Tudományos csomagok (Python Scientific stack)

csomag
NumPy Hatékony N-dimenziós tömb
SciPy Numerikus számítások
Matplotlib Grafikonok és rajzok
IPython (Jupyter) Interaktív notebook
SymPy Szimbolikus számítások
Pandas Adatbányászat
In [1]:
%matplotlib inline
import sympy as sym
import math
import numpy as np
import matplotlib.pyplot as plt  
import scipy
 

# grafikonok stílusának beállítása
plt.rcParams['figure.figsize'] = (20, 10)

Numerikus módszerek

A numerikus számítások célja, hogy a matematikai számításokat végezzünk hatékonyan.

Előnyök:

  • Gyors
  • Hasznos
  • Sok helyzetben alakalmazható

Hátrányok:

  • Csak közelítő eredményt ad
  • Emiatt figyelni kell a hiábkra
  • A matematikai intuíciónkat hátráltatja elméleti kérdésekben (pl felismernéd, hogy ez melyik algebrai szám közelítése? 2.732050807568877)

SciPy

Részletes tutorial: https://docs.scipy.org/doc/scipy/reference/tutorial/

A háttérben sokszor Fortran fut.

Alcsomagok:

scipy.cluster Vector quantization / Kmeans
scipy.constants Physical and mathematical constants
scipy.fftpack Fourier transform
scipy.integrate Integration routines
scipy.interpolate Interpolation
scipy.io Data input and output
scipy.linalg Linear algebra routines
scipy.ndimage n-dimensional image package
scipy.odr Orthogonal distance regression
scipy.optimize Optimization
scipy.signal Signal processing
scipy.sparse Sparse matrices
scipy.spatial Spatial data structures and algorithms
scipy.special Any special mathematical functions
scipy.stats Statistics

Lineáris algebra

A numpy lineáris algebra csomagját bővíti ki. Pár függvény hatékonyabb a háttérben.

LU felbontás

Adott $A$ mátrixot írjunk fel $A=LU$ alakban, ahol $L$ alsó, $U$ pedig felső háromszög mátrix.

In [2]:
from scipy import linalg, optimize
In [4]:
from scipy.linalg import lu
A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
p, l, u = lu(A)
print(A,p,l,u, sep="\n")
print(A - p @ l @ u)
np.allclose(A - p @ l @ u, np.zeros((4, 4)))
[[2 5 8 7]
 [5 2 2 8]
 [7 5 6 6]
 [5 4 4 8]]
[[0. 1. 0. 0.]
 [0. 0. 0. 1.]
 [1. 0. 0. 0.]
 [0. 0. 1. 0.]]
[[ 1.          0.          0.          0.        ]
 [ 0.28571429  1.          0.          0.        ]
 [ 0.71428571  0.12        1.          0.        ]
 [ 0.71428571 -0.44       -0.46153846  1.        ]]
[[ 7.          5.          6.          6.        ]
 [ 0.          3.57142857  6.28571429  5.28571429]
 [ 0.          0.         -1.04        3.08      ]
 [ 0.          0.          0.          7.46153846]]
[[ 0.0000000e+00  0.0000000e+00  0.0000000e+00  0.0000000e+00]
 [ 8.8817842e-16  0.0000000e+00 -4.4408921e-16  0.0000000e+00]
 [ 0.0000000e+00  0.0000000e+00  0.0000000e+00  0.0000000e+00]
 [ 8.8817842e-16  0.0000000e+00  0.0000000e+00  0.0000000e+00]]
Out[4]:
True

QR felbontás

Adott $A$ mátrixot írjunk fel $A=QR$ alakban, ahol $Q$ ortogonális, $R$ pedig felső háromszög mátrix.

In [5]:
A = np.random.randn(4, 4)
print(A)
q, r = linalg.qr(A)
print(q,r,sep='\n') 
[[-1.44840301  0.67738334  1.43656256 -1.78688189]
 [ 0.41986004  1.58434968  0.97647205 -1.36177752]
 [-1.36852932 -1.15752736  0.69033396 -0.78003637]
 [ 0.22152299 -0.29598561 -1.43680908  0.16168519]]
[[-0.70707678 -0.54274783  0.3340504  -0.30639444]
 [ 0.20496594 -0.72731261 -0.16693366  0.63335494]
 [-0.66808429  0.38019553 -0.26065223  0.58410201]
 [ 0.10814239  0.17858476  0.89028304  0.40473301]]
[[ 2.0484381   0.58709294 -1.43219775  1.52295975]
 [ 0.         -2.01291109 -1.48402198  1.69257241]
 [ 0.          0.         -1.14222561 -0.0223183 ]
 [ 0.          0.          0.         -0.70517932]]
In [6]:
q @ q.T
Out[6]:
array([[ 1.00000000e+00, -1.19749865e-17, -1.03130015e-16,
         1.50165469e-16],
       [-1.19749865e-17,  1.00000000e+00,  5.06667007e-17,
        -2.36702558e-16],
       [-1.03130015e-16,  5.06667007e-17,  1.00000000e+00,
        -1.74345312e-16],
       [ 1.50165469e-16, -2.36702558e-16, -1.74345312e-16,
         1.00000000e+00]])
In [7]:
q @ r
Out[7]:
array([[-1.44840301,  0.67738334,  1.43656256, -1.78688189],
       [ 0.41986004,  1.58434968,  0.97647205, -1.36177752],
       [-1.36852932, -1.15752736,  0.69033396, -0.78003637],
       [ 0.22152299, -0.29598561, -1.43680908,  0.16168519]])

Lineáris egyenlet megoldás

Oldjuk meg az $Ax=y$ egyeneletet.

In [8]:
A=np.array([[1,2],[3,4]])
In [9]:
y = np.array([[5.], [7.]])
x=np.linalg.solve(A, y)
x
Out[9]:
array([[-3.],
       [ 4.]])
In [10]:
A @ x
Out[10]:
array([[5.],
       [7.]])

Integrálás

Számístuk ki $\int_a^b f$ értékét.

A scipy.integrate.quad(f, a, b) integrálja $a$-tól $b$-ig az $f$ függvényt. Két értékkel tér vissza, az inregrálás eredményével és egy hibahatárral.

In [11]:
import scipy.integrate
f= lambda x:2*x
i = scipy.integrate.quad(f, 0, 1)
print(i)
(1.0, 1.1102230246251565e-14)
In [12]:
import scipy.integrate

def f(x):
    return 2*x

i = scipy.integrate.quad(f, 0, 1)
print(i)
(1.0, 1.1102230246251565e-14)
In [13]:
f= lambda x:np.exp(-x**2)
i = scipy.integrate.quad(f, 0, 1)
print(i)
(0.7468241328124271, 8.291413475940725e-15)
In [14]:
## szimbolikus számítással
xval = sym.Symbol('x')
sym.integrate(sym.exp(-xval ** 2),  xval)
Out[14]:
$\displaystyle \frac{\sqrt{\pi} \operatorname{erf}{\left(x \right)}}{2}$
In [15]:
((sym.erf(1)-sym.erf(0))*np.pi**(1/2)/2).evalf() 
Out[15]:
$\displaystyle 0.746824132812427$

Próbáljunk ki egy ronda függvényt

$$\int_0^1 e^{\sin(\log(x^2+x)-x)}dx$$
In [ ]:
#sym.integrate(sym.exp(sym.sin(-sym.log(xval ** 2+xval)-xval)),  xval)
In [16]:
f= lambda x:np.exp(np.sin(-np.log(x ** 2+x)-x)) 
i = scipy.integrate.quad(f, 0, 1)
print(i)
(1.1008683826459191, 9.950209101106111e-09)

Több változóban is tudunk integrálni a scipy.integrate.dblquad(func, a, b, gfun, hfun) paranccsal. A gfun és hfun parancsok adják meg a belső integrál határait.

$$\int_{0}^{1/2} \left( \int_{0}^{\sqrt{1-4y^2}} 16xy \:dx\right)dy$$
In [17]:
f = lambda x, y : 16*x*y
g = lambda x : 0
h = lambda y : math.sqrt(1-4*y**2)
i = scipy.integrate.dblquad(f, 0, 0.5, g, h)
print(i)
(0.5, 1.7092350012594845e-14)

Interpoláció

A fő gondolat:

  • Adott két adat tömb, $x_g$ és $y_g$, gondolhatunk rájuk úgy, hogy $y_g$ a mérési eredmény $x_g$ pontokban.
  • Egy köztes $(x,y)$ adatpontot szeretnénk megbecsülni.
  • Ha adott $x_g,y_g$ és $x$, akkor mi legyen $y$?

Módszerek

  • Illeszünk egy görbét a pontokra és értékeljük ki a görbét $x$-nél
  • Vegyük a legközelebbi adat pontot
  • Lineárisan interpoláljunk
  • Magasabb rendű polinommal interpoláljunk

Mi a különbség görbe illesztés és interpoláció között?

  • Általában az interpoláció lokális, csak néhány környező adatpontot használ
  • A görbe illesztés az összes pontra nézve próbál optimális görbét találni
  • Ezért a görbe illesztés nem fog átmenni az adatpontokon általában
  • Ezért a görbe illesztés akkor hasznos, ha az adat "zajos" és mi a zajtól szeretnénk megszabadulni vagy modelt szeretnénk felállítani
In [18]:
from scipy import interpolate
 
xg = np.linspace(0, 4, 12)
yg = np.cos(xg**2/3+4)
print (xg,'\n\n',yg)
[0.         0.36363636 0.72727273 1.09090909 1.45454545 1.81818182
 2.18181818 2.54545455 2.90909091 3.27272727 3.63636364 4.        ] 

 [-0.65364362 -0.61966189 -0.51077021 -0.31047698 -0.00715476  0.37976236
  0.76715099  0.99239518  0.85886263  0.27994201 -0.52586509 -0.99582185]
In [20]:
figure, axis = plt.subplots(1, 1,figsize=(8,8))
axis.plot(xg, yg,'o')
plt.show()
In [22]:
f1 = interpolate.interp1d(xg, yg,kind = 'linear')
In [23]:
xnew = np.linspace(0, 4,100)

figure, axis = plt.subplots(1, 1,figsize=(8,8))
axis.plot(xg, yg, 'o')
axis.plot( xnew, f1(xnew), '-')
axis.legend(['data', 'linear'], loc = 'best')

plt.show()
In [24]:
f2 = interpolate.interp1d(xg, yg, kind = 'cubic')
xnew = np.linspace(0, 4,100)
#plt.plot(x, y, 'o', xnew, f(xnew), '-', xnew, f2(xnew), '--')
figure, axis = plt.subplots(1, 1,figsize=(8,8))
axis.plot(xg, yg, 'o')
#axis.plot( xnew, f1(xnew), '-')
axis.plot( xnew, f2(xnew), '--')
 
axis.legend(['data', 'cubic'], loc = 'best')

plt.show()

Legkisebb négyzetek módszere

In [25]:
f = np.poly1d([5, 1])

x = np.linspace(0, 10, 30)
y = f(x) + 6*np.random.normal(size=len(x))
xn = np.linspace(0, 10, 200)

figure, axis = plt.subplots(1, 1,figsize=(6,6))
axis.plot(x, y, 'or')
plt.show()
In [26]:
a = np.vstack([x, np.ones(len(x))]).T
np.dot(np.linalg.inv(np.dot(a.T, a)), np.dot(a.T, y))
Out[26]:
array([4.82023058, 0.0604941 ])
In [27]:
np.linalg.lstsq(a, y,rcond=None)[0]
Out[27]:
array([4.82023058, 0.0604941 ])
In [30]:
np.polyfit(x, y, 1)
Out[30]:
array([4.82023058, 0.0604941 ])
In [31]:
m, c = np.polyfit(x, y, 1)
yn = np.polyval([m, c], xn)

figure, axis = plt.subplots(1, 1,figsize=(6,6))
axis.plot(x, y, 'or')
axis.plot(xn, yn)
plt.show()

Gyök keresés

Keressük a $2x+3cos(x)$ függvény egy gyökét.

A scipy.optimize.root(fn,x0) függvényt használhatjuk. A hátterben különböző iteratív módszerek futnak, amik folyamatosan javítanak a megoldáson. A második paraméter a kezdeti tipp a gyökre.

In [32]:
from scipy.optimize import root

def func(x):
    return 2*x + 3 * np.cos(x)

sol = root(func,0)
print(sol)
    fjac: array([[-1.]])
     fun: array([0.])
 message: 'The solution converged.'
    nfev: 9
     qtf: array([-3.65345532e-10])
       r: array([-4.37742425])
  status: 1
 success: True
       x: array([-0.91485648])
In [33]:
func(-0.91485648) 
Out[33]:
-6.797106655298535e-09
In [34]:
sym.solveset( xval*2 + 2 * sym.cos(xval), xval) # a szimolikus megoldó csak széttárja a kezét.
Out[34]:
$\displaystyle \left\{x\; \middle|\; x \in \mathbb{C} \wedge 2 x + 2 \cos{\left(x \right)} = 0 \right\}$

Optimalizáció

Adott egy célfüggvény, amit bizonyos feltélek mellett szeretnék minimalizálni (maximalizálni). Például:

$$\min(x_0^3+x_1^2-x_2)$$

Feltéve hogy: $$x_0^2+x_1^2=10 $$ $$x_0x_2\le4 $$ $$0\le x_0,x_1,x_2\le100$$

Legyen a kiindulási tipp (2,3,4)

In [35]:
from scipy.optimize import minimize
In [36]:
# cél és feltételek
def objective(x):
    return x[0]**3+x[1]**2-x[2] 

def constraint1(x):
    return x[0]**2+x[1]**2-10.0

def constraint2(x):
    return 4-x[0]*x[2]
In [37]:
#Kezdeti tipp
x0=np.array([2,3,4]) 

#Kezdeti érték
print("Cél függvény értéke a kiindulási pontban:",objective(x0))
Cél függvény értéke a kiindulási pontban: 13
In [38]:
# optimalizálás
b = (0.0,100.0)
bnds = (b, b, b)
con1 = {'type': 'eq', 'fun': constraint1}
con2 = {'type': 'ineq', 'fun': constraint2}
cons = ([con1,con2])
solution = minimize(objective,x0,method='SLSQP',bounds=bnds,constraints=cons)
x = solution.x

# show final objective
print('Cél függvény értéke: ' + str(objective(x)))

# print solution
print('Optimális értékek:')
print('x1 = ' + str(x[0]))
print('x2 = ' + str(x[1]))
print('x3 = ' + str(x[2]))
Cél függvény értéke: -89.9999999999999
Optimális értékek:
x1 = 0.0
x2 = 3.162277660168381
x3 = 99.99999999999991
In [39]:
#Ellenőrzés
print(constraint1(x))
print(constraint2(x))
8.881784197001252e-15
4.0

Közönséges differenciálegyenletek

Oldjuk meg az $f'(x)=-kf(x)$ egyenletet.

Használhatjuk az scipy.integrate.odeint függvényt.

In [40]:
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# függvény, ami visszatért a derivált értékével
def model(f,x):
    k = 0.3
    deri = -k * f
    return deri

# kezdeti érték
f0 = 5

# x értékek, ahol ki akarjuk számolni a függvényt
x = np.linspace(0,20)

# egyenelet megoldása
f = odeint(model,f0,x)

# plot results
plt.plot(x,f)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.show()

Geometriai algoritmusok

  • Alapvetően nem ajánlott, mivel a Scipy ezen a területen keveset tud

Delaunay háromszögelés

Adott egy ponthalmaz. Kössük össze azokat a pontpárokat, akikhez létezik olyan körlap, ami csak őket tartalmazza a ponthalmazból.

In [41]:
from scipy.spatial import Delaunay
points = np.random.rand(30, 2) # 30 random points in 2-D
tri = Delaunay(points)
import matplotlib.pyplot as plt
plt.triplot(points[:,0], points[:,1], tri.simplices.copy())
plt.plot(points[:,0], points[:,1], 'o')
plt.show()

Konvex burok

Adott egy ponthalmaz, számítsuk ki a konvex burkát.

In [42]:
from scipy.spatial import ConvexHull
points = np.random.rand(30, 2) # 30 random points in 2-D
hull = ConvexHull(points)
plt.plot(points[:,0], points[:,1], 'o')
for simplex in hull.simplices:
    plt.plot(points[simplex,0], points[simplex,1], 'k-')
plt.show()

Köszönöm a figyelmet!