Tekst źródłowy strony TI/Matplotlib + numpy


Krótkie wprowadzenie do biblioteki NumPy pojawiło się już w [[TI/Programowanie_z_Pythonem/Numpy|zeszłym roku]].
==Wycinki==
<noinclude>Fragment z
http://en.wikibooks.org/wiki/Non-Programmer%27s_Tutorial_for_Python_2.0/Revenge_of_the_Strings, autorstwa [[wikibooksen:User:Jrincayc]], [[wikibooksen:User:Whiteknight]], [[wikibooksen:User:Siebengang]], [[wikibooksen:User:33rogers]].<br>
Z poprawkami!
</noinclude>

==== Indeksy ====

Indeksując napis odwołujemy się do wybranego znaku. Znaki mogą być liczone od początku napisu (pierwszy ma numer 0) lub od jego końca (ostatni ma numer –1).

<tt>
{| style="text-align:center" cellspacing="0px"
|
| style="width:2ex" |
| style="width:2ex;color:red" |0
| style="width:2ex" |
| style="width:2ex;color:red" |1
| style="width:2ex" |
| style="width:2ex;color:red" |2
| style="width:2ex" |
| style="width:2ex" | ...
| style="width:2ex" |
| style="width:2ex;color:blue" | –2
| style="width:2ex" |
| style="width:2ex;color:blue" | –1
| style="width:2ex" |
|-
|
|
| style="background:#FFFF88" |&darr;
|
| style="background:#FFFF88" |&darr;
|
| style="background:#FFFF88" |&darr;
|
| style="background:#FFFF88" |&darr;
|
| style="background:#FFFF88" |&darr;
|
| style="background:#FFFF88" |&darr;
|
|
|-
| tekst = 
|"
! style="background:#FFFF88" | S
! 
! style="background:#FFFF88" | T
!
! style="background:#FFFF88" | R
!
! style="background:#FFFF88" | I
!
! style="background:#FFFF88" | N
!
! style="background:#FFFF88" | G
|"
|-
|
|
| style="background:#FFFF88" |&uarr;
|
| style="background:#FFFF88" |
|
| style="background:#FFFF88" |
|
| style="background:#FFFF88" |
|
| style="background:#FFFF88" |
|
| style="background:#FFFF88" |&uarr;
|
|
|-
|
|
| [
|
|
|
|
|
|
|
|
|
| ]
|
|}
</tt>

Na powyższym rysunku <span style="color:red">czerwone</span> indeksy są liczone od początku napisu, a <span style="color:blue">niebieskie</span> od jego końca.

Teraz przykłady<nowiki>:</nowiki>

<tt>
{|
|-
|tekst[0] || style="width:3em;text-align:center" | &rarr; || "S"
|-
|tekst[1] || style="width:3em;text-align:center" | &rarr; || "T"
|-
|tekst[2] || style="width:3em;text-align:center" | &rarr; || "R"
|-
|tekst[&minus;2] || style="width:3em;text-align:center" | &rarr; || "N"
|-
|tekst[&minus;0] || style="width:3em;text-align:center" | &rarr; || "S"
|-
|tekst[20] || style="width:3em;text-align:center" | &rarr; || {{color|red|IndexError
|-
|}
</tt>

Nie ma indeksu &minus;0, bo w Pythonie <tt>&minus;0 == 0</tt>, więc &minus;0 również wskazuje na początek napisu. 

Numery wskazujące na miejsca poza końcem napisu lub przed jego początkiem są błędne, Python zwróci informację o błędzie, tzw. <tt>IndexError</tt>.

====Przykład &mdash; indeksy dla napisu "string"====
=====Indeksy dodatnie=====
<source lang = "pycon">
>>> tekst[0]
's'
>>> tekst[1]
't'
>>> tekst[2]
'r'
>>> tekst[3]
'i'
>>> tekst[4]
'n'
>>> tekst[5]
'g'
</source>
=====Indeksy ujemne=====
<source lang = "pycon">
>>> tekst = "string"
>>> tekst[-1]
'g'
>>> tekst[-2]
'n'
>>> tekst[-3]
'i'
>>> tekst[-5]
't'
>>> tekst[-4]
'r'
>>> tekst[-5]
't'
>>> tekst[-6]
's'
</source>
==== Wycinki ====

Pobranie wycinka sekwencji polega na wybraniu jej fragmentu &mdash; podsekwencji.
Miejsca początku i końca zakresu z którego pobieramy wycinek liczone są podobnie jak przypadku indeksowania. Istotna różnica jest taka, że wskazujemy podsekwencję, a nie pojedynczy element &mdash; i to w sposób szczególny. '''Przedział który wskazujemy jako <tt>[a:b]</tt>, jest zamknięty od strony <tt>a</tt>, a otwarty od strony <tt>b</tt>.'''
<!-- piszę: od strony a, od strony b, zamiast od lewej, od prawej, bo dla kroku ujemnego strony się zamieniają w przedziale, natomiast nie we wskazaniu i nastąpiłaby konfuzja. -->

Lewy i prawy kraniec wycinka podajemy oddzielając indeksy dwukropkiem <tt>:</tt>. Każdą z tych wartości można pominąć, w skrajnym przypadku zostawiając tylko dwukropek. Jeśli dany indeks pominiemy, zostaną użyte wartości domyślne:
* domyślne położenie początku wycinka to początek sekwencji
* domyślne położenie końca wycinka to koniec całej sekwencji

Taka konwencja ma pewne zalety, w szczególności dla dowolnego indeksu <tt>a</tt> prawdą jest, że
<source lang= python>sekwencja[:a] + sekwencja[a:] == sekwencja</source>

Konstruując wycinek możemy pobierać co ''n''-ty element napisu. Robi się to podając wartość kroku ''n'' po kolejnym dwukropku. Domyślnie ''n'' ma wartość 1.

Teraz przykłady<nowiki>:</nowiki>

<tt>
{|
|-
|tekst[1:4] || style="width:3em;text-align:center" | &rarr; || "TRI"
|-
|tekst[:5] || style="width:3em;text-align:center" | &rarr; || "STRIN"
|-
|tekst[:–1] || style="width:3em;text-align:center" | &rarr; || "STRIN"
|-
|tekst[–4:] || style="width:3em;text-align:center" | &rarr; || "RING"
|-
|tekst[:] || style="width:3em;text-align:center" | &rarr; || "STRING"  (tzw. ''pełny wycinek'', czyli kopia sekwencji)
|-
|tekst[::–1] || style="width:3em;text-align:center" | &rarr; || "GNIRTS"
|-
|tekst[&minus;10:&minus;2] || style="width:3em;text-align:center" | &rarr; || "STRI"
|}
</tt>

Numery wskazujące na miejsca poza rozmiarem napisu nie są błędne, Python wybierze tylko istniejące znaki.

== Numpy ==
=== Tablice ===
W zeszłym roku była mowa o tablicach (arrays) w numpy. Ograniczaliśmy się praktycznie do przypadku dwuwymiarowego, czyli macierzy. Warto sobie uświadomić, że obiekt array może być wielowymiarową tablicą. Przytoczymy teraz kilka przydatnych funkcji dotyczących tablic w numpy.
<source lang="python">
import numpy as np

a = np.ones((2,3,4)) # tworzy tablicę jedynek o wymiarze 2x3x4 
a = np.empty((2,3)) # tworzy pustą tablicę o wymiarze 2x3
a = np.zeros((2,3,4,5)) # twrzy tablicę zer o wymiarze 2x3x4x5
# tak, w numpy można tworzyć więcej niż 3-wymiarowe tablice!
# jeśli potrzebowalibyśmy bardzo dużo wymiarów, a nie chce nam się ich po kolei wpisywać możemy zrobić tak:
a = np.empty(tuple([10 for i in range(20)])) # 20-sto wymiarowa tablica, w której każdy wymiar ma długość 10
#np.empty przyjmuje jako argument obiekt tuple, dlatego utworzyliśmy w wygodny sposób listę wymiarów, a następnie przerabiamy ją na tuple, żeby podać jako argument
# łatwiejszym do wyobrażenia, 3-wymiarowym odpowiednikiem tego jest:
a = np.empty(tuple([10 for i in range(3)])) # 3-wymiarowa tablica, w której każdy z wymiarów ma długość 10
A = np.random.randint(10,size = (3,3)) # wytwarza tablicę losowych elementów całkowitych o wymiarach 3x3 (takie tablicę częst się przydają jako dane testowe do naszych funkcji/programów)
A = np.random.randint(10,size = (3,3,5)) # wytwarza tablicę analogiczną do powyższej, ale o wymiarze 3x3x5
wiersze,kolumny = numpy.shape(A) # co się stanie, gdy A ma wymiar mniejszy/większy niż 2?
</source>

* Jaki będzie wynik następujących instrukcji:
<source lang = "python">
a = np.arange(1,10)
</source>
<source lang = "python">
a = np.arange(10)
</source>
<source lang = "python">
a = np.arange(1., 10)
</source>
<source lang = "python">
a = np.arange(1, 10, 2.5)
</source>
<source lang = "python">
a = np.arange(1,2,10)
</source>
* Utwórz 3-wymiarową tablicę A, o wymiarach 10x5x3 taką, że A[i,j,k] = "i,j,k" (wartością każdego elementu będzie string zawierający wszystkie indeksy). Aby zainicjalizować poprawnie tablicę pamiętaj o wybraniu odpowiedniego typu zmiennych (<tt>dtype</tt>). Wypisz ją, przyjrzyj się w jakiej kolejności są wypisywane elementy, gdy napiszesz po prostu.
<source lang="python">
print A
</source>
* Utwórz 100-wymiarową tablicę zer, o każdym wymiarze = 10. Ile maksymalnie wymiarów może mieć tablica w numpy?
* Z listy [[1, 2, 3], [4, 5, 6],[7, 8, 9]] utwórz macierz.
* Zastanów się, jak zmienić swoją funkcję obliczającą splot, tak, żeby radziła sobie też z wektorami różnej długości.

===Trochę wiadomości o tablicach w <tt>numpy</tt>===
* W <tt>numpy</tt> jak i w Pythonie obowiązuje numeracja od zera &mdash; stąd pierwszy element wektora <tt>a</tt> o długości ''N'' to <tt>a[0]</tt>, a ostatni <tt>a[N-1]</tt>. 
* Elementy tablicy indeksuje się: <math>A[i_0,i_1,...,i_{n-1}]</math>, gdzie <math>i_j</math> to indeks w wymiarze ''j''. Tzn dla tablicy dwuwymiarowej: <tt>A[wiersz, kolumna]</tt>.
[[Grafika:macierz.svg|right]]
* Łączenie tablic wykonuje się za pomocą funkcji <tt>numpy.concatenate((a1, a2, ...), axis=0)</tt>. <tt>a1, a2, ...</tt> oznacza tablice, a <tt>axis</tt> oś, w której tablice zostaną połączone. Parametr ten ma domyślnie wartość 0. Tablice muszą mieć taki sam rozmiar w osi, w której będą połączone. Funkcja zwraca tablicę wynikową.
* Usuwanie elementu z tablicy (wielowymiarowej):<source lang="python">dane2 = dane[numpy.r_[:i,i+1:dane_1.shape[0]]] </source>
==Problemy==
===Proste operacje===
<source lang="python">
from numpy import *

l1 = [1,1,1]
12 = ones(3)
</source>
Jaki będzie wynik następujących operacji:
<source lang="python">l1 * 3</source>
<source lang="python">l2 * 3</source>
Jakie operacje należy wykonać na liście l1, aby otrzymać wynik analogiczny do l2*3?
<source lang="python">l1 / 3</source>
<source lang="python">l2 / 3</source>
===Plansza kółko-krzyżyk===
Zaprojektuj tablicę złużącą do gry w kółko i krzyżyk w trzech wymiarach. Kółka to zera, krzyżyki to jedynki. To znaczy zdefiniuj kilka funkcji, które będą umożliwiały grę w kółko i krzyżyk w trzech wymiarach:
* funkcję tworzącą pustą planszę
* funkcję stawiąjącą krzyżyk w wybranym miejscu(określonym współrzędnymi x,y,z)
* funkcję stawiąjącą kółko w wybranym miejscu
* funkcję wypisującą aktualny stan planszy
* Rozszerz poprzednie zadanie o wyrysowywanie trójwymiarowe stanu planszy.
* Utwórz ze zdefiniowanych funkcji moduł
===Dopasowanie metodą najmniejsztch kwadratów===
W pracy doświadczalnej wielokrotnie badamy zależność z jaką wielkość <math>y</math> zależy od wielkości <math>x</math>. 
Zwykle ta zależność fizyczna jest postaci:
<equation id="eq:ocz">
<math>y = f(x;par)</math>
</equation>
gdzie <math>f</math> &mdash; zależność fizyczna, <math>par</math> &mdash; parametry.

Żeby zobaczyć, czy funkcja <math>f(x;par)</math> dobrze oddaje realną zależność fizyczną między <math>y</math> a <math>x</math> obliczamy suma kwadratów residuów (<math>R</math>) &mdash; sumę kwadratów różnic wyników pomiaru i ich wartości oczekiwanych, wyliczonych z (<xr id="eq:ocz"/>):
<equation id="eq:ch2">
<math>R = \sum_{i=1}^N \left( y_i-\mathrm{oczekiwane}_i\right)^2= \sum_{i=1}^N \left( y_i-f(x_i;par)\right)^2</math>
</equation>
gdzie para <math>(x_i,y_i)</math> oznacza wynik piątego pomiaru wielkości <math>x</math> i <math>y</math>.

Szukając parametrów <math>par</math> (<xr id="eq:ocz"/>), wybieramy te, dla których <math>R</math> jest najmniejsze. Jest wiele różnych algorytmów, którymi można się w tym posłużyć. Można też po prostu przeskanować przestrzeń parametrów, w szczególności, kiedy wiemy, jakie są fizyczne zakresy wartości tychże; i wybrać takie, które minimalizują <math>R</math>.

Błąd dopasowania wybranego już parametru <math>\theta</math> przedstawia się następującym wzorem:
<equation id="eq:deltap">
<math>\Delta \theta = \sqrt{\frac{1}{\frac 1 2 \frac{\mathrm \partial^2 R}{\mathrm \partial \theta^2}</math>.
</equation>
====Zadanie====
Napisz moduł zawierający funkcję do obliczenia sumy kwadratów residuów (<xr id="eq:ch2"/>) i błędu dopasowania parameterów (<xr id="eq:deltap"/>).
#Załóż, że funcja opisująca relację między danymi ma postać (<xr id="eq:ocz"/>) <source lang="python">func(dane,parametry)</source>.
#Dane znajdują się w macierzy o wymiarach ''n''&times;2 (<tt>dane</tt>).
#Funkcja obliczająca (<xr id="eq:ch2"/>) powinna mieć przekazywane dane w wektorze o wymiarach (<tt>dane</tt>), nazwę funkcji charakteryzującej zależność fizyczną i służącej do obliczania wartości oczekiwanych (<tt>func</tt>) i listę parametrów funkcji (<tt>params</tt>):<source lang="python"> def chi2(dane, func, params): </source><b>Pamietaj</b> o rozpakowaniu listy z parametrami funkcji <tt>func</tt>.
#Funkcja obliczająca (<xr id="eq:deltap"/>) powinna zwracać niepewności wszystkich podanych w liście parametrów. Parametrem opcjonalnym powinno być <math>\Delta</math>, służące do obliczenia drugiej pochodnej (patrz wskazówka 1). Funkcja powinna korzystać z funkcji <tt>chi2</tt>. Nagłówek funkcji:<source lang="python">def delta_p(dane, func, params,delta=0.000001):</source>
#Metodą ''brute force'', skanując przestrzeń parametrów, sprawdź, czy do danych pasuje zależność (<xr id="eq:FR"/>) (opis danych i parametrów w tytule rysunku).
#Zobacz na wykresie, jak dobre jest Twoje dopasowanie.
#Zmodyfikuj funkcję, tak żeby mogła uwzględnić niepewności pomiarowe w obliczaniu sumy kwadratów residuów: <math>R =  \sum_{i=1}^N \left( \frac{y_i-f(x_i;par)}{\sigma_i}\right)^2</math>, gdzie <math>\sigma_i</math> ''i''-tą niepewność pomiarową <math>\mathrm{pomiar}_i</math>. Załóż, że jeżeli tablica z danymi ma trzy kolumny, w trzeciej znajdują się niepewności pomiarowe <math>\delta y</math>. Jeżeli tablica z danymi ma cztery kolumny, w czwartej znajdują się niepewności pomiarowe ''y''-ka &mdash; <math>\delta y</math>, a w trzeciej niepewności pomiarowe ''x''-a &mdash; <math>\delta x</math>.

Dane znajdują się w [[Plik:Dane.txt]]. Pierwsza kolumna przedstawia interwały <math>T</math>. Druga obliczone wartości wzmocnienia synaptycznego ''FR'', które może być opisywane następującą zależnością:
<equation id="eq:FR">
<math> FR=1+\left( 1+(1-U_\mathrm{SE}) e^{-\frac{T}{\tau_\mathrm{f}\right)^4</math>
</equation>
*<math>\tau_\mathrm{f} = \unit{5 - 100}{ms}</math>
*<math>U_\mathrm{SE} = 0.001 - 0.03</math>
W jakim zakresie <tt>Delta</tt> wynik niepewności jest stabilny?


====Wskazówka 1====
Numeryczne obliczanie pierwszej pochodnej:
::<math>\frac{\mathrm d x}{\mathrm dt} = \frac{x(t+\Delta t)-x(t)}{\Delta t}</math>
Numeryczne obliczanie drugiej pochodnej:
::<math>\frac{\mathrm d^2 x}{\mathrm dt^2} = \frac{f'(x)-f'(x-\Delta t)}{\Delta t}=\frac{x(t+\Delta t) +x(t-\Delta t) -2x(t)}{\Delta t^2}</math>
====Wskazówka 2====
Skorzystaj z modułu <tt>Numpy</tt> do wczytania danych (<tt>loadtxt</tt>) oraz obliczenia wartości funkcji wykładniczej (<tt>exp</tt>).



===Zadanie_z_chi2/rozwiązanie macierzowe===
Najpierw definiujemy funkcję służącą do obliczania FR na podstawie modelu teoretycznego:
<source lang=python>
from __future__ import division
import numpy

def FRobl(tau, Use, T):
    return 1 + (1 + (1 - Use)*numpy.exp(-T/tau))**4
</source>
Proszę zwrócić uwagę, że tak zdefiniowana funkcja działa równie dobrze dla zwykłych liczb, jak i dla tablic <tt>numpy.ndarray</tt>, o ile podane argumenty obsługują działania wykorzystywane w funkcji, czyli <tt>/</tt>,
<tt>*</tt>,  <tt>**</tt>,  <tt>-</tt>,  <tt>+</tt> i <tt>numpy.exp</tt>.<br>

Z kolei <tt>__future__.division</tt> jest importowane bo chce się upewnić, że dzielenie zawsze jest zwykłym dzieleniem, nawet jeśli funkcja zostanie wywołana jako <tt>FRobl(2, 0.03, 1)</tt>.

Następnie definiujemy funkcje do obliczania pierwszej i drugiej pochodnej:
<source lang=python>
def pochodna(f, i, *params, **opts):
    eps = opts.get('eps', 1e-6)
    params_eps = list(params)
    params_eps[i] += eps
    return (f(*params_eps) - f(*params))/eps

def pochodna2(f, i, *params, **opts):
    eps = opts.get('eps', 1e-6)
    params_eps = list(params)
    params_eps[i] += eps
    return (pochodna(f, i, *params_eps) -
            pochodna(f, i, *params)) / eps
</source>
Jeśli się pamięta wzór na pochodną, to tutaj nie ma wielkich niespodzianek.

Jak działa <tt>*params, **opts</tt> w liście parametrów?<br>
Mamy argument <tt>params</tt>, który zawierać będzie wszystkie nadmiarowe argumenty pozycyjne, oraz <tt>opts</tt> który zawierać będzie wszystkie nazwane argumenty. Dzięki temu możemy funkcję wywołać z dowolną liczbą parametrów <tt>params</tt>, które zostaną przekazane do <tt>f</tt>.

Do <tt>opts</tt> trafiają argumenty nazwane. Jeśli napiszemy <tt>eps=...</tt>, to nasz argument zostanie wydobyty ze słownika i użyty jako <tt>eps</tt>. Jeśli nie podamy argumentu <tt>eps=...</tt>, to jako <tt>eps</tt> zostanie użyty
drugi argument do <tt>get</tt>, czyli 1e-6. Jeśli podamy argumenty nazwane inne niż  <tt>eps=...</tt> to zostaną one zignorowane -- niestety nasza funkcja nie jest doskonała.

Aby obliczyć niepewność parametru, musimy mieć funkcję obliczającą rms (<tt>rmsfunc</tt>) i miejsce w którym liczymy (<tt>params</tt>) i numer parametru który nas interesuje (<tt>i</tt>).
<source lang=python>
def niepewnosc(rmsfunc, i, *params):
    return (pochodna2(rmsfunc, i, *params)/2)**-0.5
</source>

W naszym wypadku funkcja obliczająca rms może wyglądać tak:
<source lang=python>
def frms(T, FR, tau, Use):
    FR2 = FRobl(tau, Use, T)
    return ((FR2 - FR)**2).mean()**0.5
</source>
Sens jest taki, że obliczamy co daje model teoretyczny <tt>FRobl</tt> jak dla danego zestawu (<math>\tau</math>, <math>U_{SE}</math>) = (tau, Use) i pewnego <tt>T</tt> i porównujemy z wynikami doświadczalnymi <tt>FR</tt>. Otrzymujemy zestaw odchyleń, a rms jest oczywiście pierwiastkiem ze średniej kwadratów tych odchyleń.

Przechodąc do meritum, możemy albo się posłużyć danymi doświadczalnymi, albo np. wygenerować dane dla znanych ($\tau$, $U_{SE}$). To drugie podejście pozwala nam sprawdzić naszą metodę, bo wiemy dokładnie jaki wynik powinniśmy otrzymać.

Zestaw danych wygenerowanych wraz z parametrami:
<source lang=python>
def dane_symulowane(tau, Use):
    T = numpy.arange(5, 80, 5)
    return T, FRobl(tau, Use, T), tau, Use

T_opt, FR_opt, tau_opt, Use_opt = dane_symulowane(30, 0.02)
</source>

Dopasowanie wykonujemy metodą <i>brute force</i>, czyli po prostu przeszukujemy całą przestrzeń parametrów.
Implementację rozbijmy na dwie części &mdash; najpierw obliczanie rms, potem znajdywanie minimum:
<source lang=python>
def rms_brute_force(T, FR):
    tau = numpy.arange(5, 101, 1)
    Use = numpy.arange(0.001, 0.031, 0.001)

    # tau3 = tau
    # Use3 = Use[:, numpy.newaxis]
    # T3 = T[:, numpy.newaxis, numpy.newaxis]

    tau3 = tau[numpy.newaxis, numpy.newaxis, :]
    Use3 = Use[numpy.newaxis, :, numpy.newaxis]
    T3 = T[:, numpy.newaxis, numpy.newaxis]

    print('kształty tau, Use, T: ',
          tau3.shape, Use3.shape, T3.shape)

    wynik = FRobl(tau3, Use3, T3)

    FR3 = FR[:, numpy.newaxis, numpy.newaxis]
    odchylenie = wynik - FR3
    rms = (odchylenie**2).mean(axis=0)
    return tau, Use, rms
</source>
Do szybkiego przeprowadzenia obliczeń na całej siatce wykorzystujemy broadcasting. Tworzymy trzy wektory, czyli tablice o długości większej niż 1 w jednym tylko kierunku &mdash; (1, 1, M2), (1, M1, 1), (N, 1, 1). Każdy z wektorów jest skierowany w inną stronę, więc po wykonaniu działania element po elemencie (tak jak to się dzieje w <tt>FRobl</tt>), otrzymujemy macierz o wymiarze (N, M1, M2).

<source lang=python>
def dopasowanie_brute_force(T, FR):
    tau, Use, rms = rms_brute_force(T, FR)
    ind_min = numpy.unravel_index(rms.argmin(), rms.shape)
    Use_min = Use[ind_min[0]]
    tau_min = tau[ind_min[1]]
    print('minimum dopasowania to τ={} Use={}'.format(tau_min, Use_min))
    assert frms(T, FR, tau_min, Use_min) == rms.min()

    D_tau_min = niepewnosc(frms, 2, *(T, FR, tau_min, Use_min))
    D_Use_min = niepewnosc(frms, 3, *(T, FR, tau_min, Use_min))
    print('niepewności wynoszą Δτ={} ΔUse={}'.format(D_Use_min, D_tau_min))

    return tau, Use, rms, tau_min, D_tau_min, Use_min, D_Use_min
</source>
Aby znaleźć optymalne wartości parametrów, czyli takie dla których rms jest najmniejszy, wykorzystujemy 
<tt>numpy.argmin</tt>. Ta funkcja nam zwraca indeks minimalnego elementu, ale po "spłaszczeniu". Numer wiersza i kolumny uzyskujmy przez <tt>numpy.unravel_index</tt>.

Do obliczenia niepewności mamy zdefiniowaną wcześniej funkcję <tt>niepewnosc</tt>

Wynik:
<source lang=python>
tau, Use, rms, tau_min, D_tau_min, Use_min, D_use_min = dopasowanie_brute_force(T_opt, FR_opt)
</tt>

Wyniki możemy narysować przy pomocy następujących funkcji, które zostawie prawie bez komentarza:
<source lang=python>
from __future__ import unicode_literals
from matplotlib import pyplot, patches
__name__ == '__main__' and pyplot.ion()
# unicode_literals powoduje, że matplotlib wie, że wszystkie stringi to Unicode
# pyplot.ion() podowuje coś

def plot_map(tau, Use, rms):
    f = pyplot.figure()
    ax = f.add_subplot(111)
    X, Y = numpy.meshgrid(tau, Use)
    map = ax.pcolor(X, Y, rms)
    ax.set_xlabel(r'tau')
    ax.set_ylabel(r'Use')
    f.colorbar(map)
    f.show()
    return f
</source>

[[Plik:plot_map.svg]]

<source lang=python>
def plot_map2(tau, Use, rms, **kwargs):
    f = pyplot.figure()
    ax = f.add_subplot(111)
    X, Y = numpy.meshgrid(tau, Use)
    cont = ax.contour(X, Y, rms)
    ax.set_xlabel(r'tau')
    ax.set_ylabel(r'Use')
    f.colorbar(cont)
    ax.clabel(cont, inline=1, **kwargs)
    f.show()
    return f
</source>

[[Plik:plot_map2.svg]]

<source lang=python>
def plot_dane(T, FR, label1, T2, FR2, label2):
    f = pyplot.figure()
    ax = f.add_subplot(111)
    ax.set_xlabel(r'T')
    ax.set_ylabel(r'FR')
    ax.plot(T, FR, 'ro', label=label1)
    ax.plot(T2, FR2, '-g', label=label2)
    ax.legend(loc='best')
    return f

def plot_dopasowanie(T, FR, tau, Use):
    T2 = numpy.linspace(T.min() - 1, T.max() + 10, 300)
    FR2 = FRobl(tau, Use, T2)
    return plot_dane(T, FR, 'dane', T2, FR2, 'dopasowanie')
</source>

[[Plik:plot_dopasowanie.svg]]

I na koniec rysunek z zeznaczonym rms:
<source lang=python>
def plot_dopasowanie_map(tau, Use, rms,
                         tau_min, D_tau_min, Use_min, D_Use_min):
    f = plot_map(tau, Use, rms)
    ax = f.gca()

    ax.plot(tau_min, Use_min, 'xy', markersize=5)

    ellip = patches.Ellipse(xy=(tau_min, Use_min),
                            width=D_tau_min, height=D_Use_min)
    ellip.set_facecolor('none')
    ellip.set_linestyle('dashed')
    ellip.set_edgecolor('yellow')
    ax.add_artist(ellip)
    return f
</source>

[[Plik:plot_dopasowanie_map.svg]]


===Działania macierzowe===
Broadcasting (dopełnianie) jest słowem, które określa, jak numpy traktuje wektory o różnym kształcie. Przy spełnieniu pewnych warunków, mniejszy wektor jest "rozszerzony" (powiększony) tak, żeby mieć kształt korespondujący z kształtem większego wektora &mdash; odpowiednie wymiary zostają dopełnione. Dopełnienie umożliwia przeprowadzenie operacji macierzowych tak, że pętle są wykonywane w C, a nie w Pythonie, co zwykle prowadzi do bardziej efektywnych implementacji algorytmów.

[[Plik:Broadcasting_1D.png|thumb|<figure id="fig:rys_b1"/>Dodanie liczby 3 do wektora.]]
Najprostszym przypadkiem broadcastingu jest dodawanie liczby (skalara) do macierzy. W tej sytuacji, każdy element macierzy zostanie zwiększony o tęże liczbę. Przyjrzyj się następującemu przykładowi i obrazkowi <xr id="fig:rys_b1"/>:
<source lang="pycon">
>>> x = np.arange(4)
>>> x = array([0, 1, 2, 3])
>>> x + 3
array([3, 4, 5, 6])
</source>

[[Plik:Broadcasting_2D.png|thumb|<figure id="fig:rys_b2"/>Dopełnienie jednowymiarowego wektora tak, żeby można było go dodać do wektora dwuwymiarowego.]]
W przypadku więcejwymiarowym, oś, która powinna zostać dodana &mdash; w którą stronę wektor powinien zostać dopełniony, zostaje zaznaczona ''explicite'' słowem <tt>np.newaxis</tt> (albo <tt>None</tt>). Tak, jak na poniższym przykładzie i Rys. <xr id="fig:rys_b2"/>:
<source lang="pycon">
>>> a = np.arange(12).reshape((3,4))
>>> b = np.array([1,2,3])[:,np.newaxis]
>>> a + b
array([[ 1,  2,  3,  4],
       [ 6,  7,  8,  9],
       [11, 12, 13, 14]])
</source>
[[Plik:Broadcasting_3D.png|thumb|<figure id="fig:rys_b3"/>Dopełnienie dwuwymiarowego wektora tak, żeby po dodaniu wektora jednowymiarowego, utworzyć trójwymiarowy wektor.]]
Przypadek trójwymiarowy (Rys. <xr id="fig:rys_b3"/>):
<source lang="pycon">
>>> x = np.zeros((3, 5))
>>> y = np.zeros(8)
>>> (x[..., None] + y).shape
(3, 5, 8)
</source>

Można tez skorzystać z funkcji [http://docs.scipy.org/doc/numpy/reference/generated/numpy.broadcast_arrays.html <tt>np.broadcast_arrays</tt>].

====Zadanie 1====
Przed przystąpieniem do zadania polecamy zapoznać się z przykładami [http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html stąd]

====Zadanie 2==== 
Przepisz operacje zachodzące w pętlach i funkcjach w [[TI/Matplotlib_%2B_numpy#Dopasowanie_metod.C4.85_najmniejsztch_kwadrat.C3.B3w|zadaniu z dopasowaniem metodą najmniejszych kwadratów]] tak, żeby skorzystać z operacji macierzowych.

Żeby znaleźć najmniejszą wartość w macierzy z kwadratami residuów skorzystaj z metody <tt>min()</tt> strunktury <tt>ndarray</tt>. Żeby uzyskać spłaszczony indeks tego elementu skorzystaj z metody <tt>argmin()</tt> strunktury <tt>ndarray</tt> i następnie z funkcji [http://docs.scipy.org/doc/numpy/reference/generated/numpy.unravel_index.html#numpy.unravel_index <tt>numpy.unravel_index(indices, dims, order='C')</tt>], żeby znaleźć wartości parametrów, które minimalizują sumę kwadratów residuów.

===Deklarowanie własnych typów danych===
Deklarowanie własnych typów danych odbywa się za pomocą funkcji <tt>dtype</tt>.

Przykładowo (przykłady z dokumentacji):
<source lang="pycon">
>>> np.dtype([('f1', np.int16)])
dtype([('f1', '<i2')])
</source>
Daje nam rekord o polu 'f1' o typie danych <tt>int16</tt>.
<source lang="pycon">
np.dtype([('f1', np.uint), ('f2', np.int32)])
dtype([('f1', '<u4'), ('f2', '<i4')])
</source>
Daje nam dwa pola rekordu jedno typu unsigned <tt>int</tt>, drugie typu <tt>int32</tt>.

Można też zadeklarować rekord zawierający liczbę typu <tt>double</tt> i string dziesięciobajtowy:
<source lang="pycon">
np.dtype([('a','f8'),('b','S10')])
dtype([('a', '<f8'), ('b', '|S10')])
</source>

Przykład użycia:
<source lang="pycon">
>>> typ = np.dtype([("a", np.int16), ("b", "|S10")])
>>> typ
dtype([('a', '<i2'), ('b', '|S10')])
>>> a = np.array([(5, "napis")], dtype = typ)
>>> a
array([(5, 'napis')], 
      dtype=[('a', '<i2'), ('b', '|S10')])
>>> a = np.array([("n", "napis")], dtype = typ)

Traceback (most recent call last):
  File "<pyshell#171>", line 1, in <module>
    a = np.array([("n", "napis")], dtype = typ)
TypeError: expected a readable buffer object

>>> a = np.array([(5, 5)], dtype = typ)
>>> a
array([(5, '5')], 
      dtype=[('a', '<i2'), ('b', '|S10')])
>>> a = np.array([(5,)], dtype = typ)

Traceback (most recent call last):
  File "<pyshell#174>", line 1, in <module>
    a = np.array([(5,)], dtype = typ)
TypeError: expected a readable buffer object


>>> typ2
dtype([('c', '<i2')])
>>> typ2 = np.dtype([("c", (np.int16, 3))])
>>> a = np.array([((3.5,1),)], dtype = typ2)

Traceback (most recent call last):
  File "<pyshell#218>", line 1, in <module>
    a = np.array([((3.5,1),)], dtype = typ2)
TypeError: expected a readable buffer object
>>> a = np.array([((3,3,3),)], dtype = typ2)
>>> a
array([([3, 3, 3],)], 
      dtype=[('c', '<i2', 3)])

>>> typ3 = np.dtype([("c", (np.int16, 5))])
>>> a = np.array([((3,3,4,5),)], dtype = typ3)

Traceback (most recent call last):
  File "<pyshell#225>", line 1, in <module>
    a = np.array([((3,3,4,5),)], dtype = typ3)
TypeError: expected a readable buffer object
>>> a = np.array([((3,3,4,5,6),)], dtype = typ3)
>>> a
array([([3, 3, 4, 5, 6],)], 
      dtype=[('c', '<i2', 5)])
</source>
====Zadanie 1==== 
Pod [http://escher.fuw.edu.pl/~mm/KAMIL/kamil_eksp_1290627532.3900001.txt linkiem] znajduje się plik tekstowy zawierający dane o triggerach z jednego z eksperymentów. Obejrzyj dane. Zaprojektuj za pomocą funkcji <tt>dtype</tt> z modułu <tt>numpy</tt> typ danych, który opisuję dane, znajdujące się w rzędach. Następnie wczytaj je za pomocą funkcji <tt>loadtxt</tt>, pamietając o opuszczeniu pierwszej linii i zadaniu przecinków, jako separatorów między kolumnami i wypisz wektor z danymi. Program nie powinien mieć więcej niż 5 linii.

====Zadanie 2==== 
Wczytaj dane z poprzedniego zadania i spod [http://escher.fuw.edu.pl/~mm/KAMIL/kamil_eksp_1290760577.937.txt łącza], a następnie wypisz średnie czasy odpowiedzi dla bodźców typu pozytywnego, neutralnego i negatywnego