#!/usr/bin/env python
# coding: utf-8

# Numpy
# =====

# `numpy` je paket (modul) za (efikasno) numeričko računanje u Pythonu. Naglasak je na efikasnom računanju s nizovima, vektorima i matricama, uključivo višedimenzionalne stukture. Napisan je u C-u i Fortanu te koristi BLAS biblioteku.

# In[21]:


from numpy import *


# ## Kreiranje nizova pomoću `numpy` modula

# In[22]:


v = array([1,2,3,4])
v


# In[23]:


M = array([[1, 2], [3, 4]])
M


# In[24]:


type(v), type(M)


# In[25]:


v.shape


# In[26]:


M.shape


# In[27]:


v.size, M.size


# Možemo koristiti i funkcije `numpy.shape`, `numpy.size`

# In[28]:


shape(M)


# In[29]:


size(M)


# Koja je razlika između `numpy.ndarray` tipa i standardnih lista u Pythonu?
# 
# - liste u Pythonu mogu sadržavati bilo kakve vrste objekata, to nije slučaj s `numpy.ndarray`
# - `numpy.ndarray` nisu dinamički objekti: pri kreiranju im je određen tip
# - za `numpy.ndarray` implementirane su razne efikasne metode važne u numerici
# - de facto sva računanja se odvijaju u C-u i Fortranu pomoću BLAS rutina
# 
# `dtype` (data type) nam daje informaciju o tipu podataka u nizu:

# In[30]:


M.dtype


# Kako je *M* statički objekt, ne možemo napraviti ovo:

# In[31]:


M[0,0] = "hello"


# Naravno, ovo je ok:

# In[32]:


M[0,0]=5


# `dtype` se može eksplicitno zadati:

# In[33]:


M = array([[1, 2], [3, 4]], dtype=complex)
M


# Tipično `dtype` su: `int`, `float`, `complex`, `bool`, `object`, itd.
# 
# Ali možemo biti i eksplicitni vezano za veličinu registra: `int64`, `int16`, `float128`, `complex128`.

# ### Funkcije koje generiraju nizove

# In[34]:


x = arange(0, 10, 1) # argumenti: početak, kraj, korak

x # 10 nije u nizu!


# In[35]:


x = arange(-1, 1, 0.1)
x


# In[36]:


# ovdje su i početak i kraj uključeni!
linspace(0, 10, 25)


# In[37]:


logspace(0, 10, 10, base=e)


# In[38]:


x, y = mgrid[0:5, 0:5] # slično kao meshgrid u MATLAB-u


# In[39]:


x


# In[40]:


y


# In[41]:


from numpy import random


# In[42]:


# uniformna distribucija na [0,1]
random.rand(5,5)


# In[43]:


# standardna normalna distribucija
random.randn(5,5)


# In[44]:


# dijagonalna matrica
diag([1,2,3])


# In[45]:


# matrica sa sporednom dijagonalom
diag([1,2,3], k=1)


# In[46]:


zeros((3,3))


# In[47]:


ones((3,3))


# ### Učitavanje podataka

# Često učitavamo podatke iz datoteka (lokalno ili s weba). Važni formati su cvs (comma-separated values) i tsv (tab-separated values).

# In[48]:


get_ipython().system('head tpt-europe.csv')


# In[49]:


data = genfromtxt('tpt-europe.csv')


# In[50]:


data.shape, data.dtype


# Uz `numpy.savetxt` možemo napraviti i obrnuto.

# In[51]:


M = random.rand(3,3)
M


# In[52]:


savetxt("random-matrix.csv", M)


# In[53]:


get_ipython().system('cat random-matrix.csv')


# In[54]:


savetxt("random-matrix.csv", M, fmt='%.5f') # s fmt specificiramo format

get_ipython().system('cat random-matrix.csv')


# Postoji i interni format za `numpy` nizove:

# In[55]:


save("random-matrix.npy", M)

get_ipython().system('file random-matrix.npy')


# In[56]:


load("random-matrix.npy")


# In[57]:


M.itemsize # byte-ovi po elementu


# In[58]:


M.nbytes


# In[59]:


M.ndim


# ## Rad s nizovima

# Indeksiranje funkcionira standardno.

# In[60]:


v[0]


# In[61]:


M[1,1]


# In[62]:


M


# In[63]:


M[1]


# Naravno, možemo koristiti i `:` operator: 

# In[64]:


M[1,:] # redak 1


# In[65]:


M[:,1] # stupac 1


# In[66]:


M[1,:] = 0
M[:,2] = -1


# In[67]:


M


# In[68]:


A = array([1,2,3,4,5])
A


# In[69]:


A[1:3]


# In[70]:


A[1:3] = [-2,-3]
A


# In[71]:


A[::]


# In[72]:


A[::2]


# In[73]:


A[:3]


# In[74]:


A[3:]


# S negativnim indeksima računamo od kraja niza:

# In[75]:


A = array([1,2,3,4,5])


# In[76]:


A[-1] # zadnji element niza


# In[77]:


A[-3:] # zadnja tri elementa


# Naravno, iste operacije imamo i za višedimenzionalne nizove.

# In[78]:


A = array([[n+m*10 for n in range(5)] for m in range(5)])
A


# In[79]:


A[1:4, 1:4]


# In[80]:


A[::2, ::2]


# In[81]:


indeksi_redaka = [1, 2, 3]
A[indeksi_redaka]


# In[82]:


indeksi_stupaca = [1, 2, -1]
A[indeksi_redaka, indeksi_stupaca]


# Možemo koristiti i tzv. maske: ako je maska `numpy` niz tipa `bool`, tada se izabiru oni elementi koji u maski odgovaraju vrijednosti `True`.

# In[83]:


B = array([n for n in range(5)])
B


# In[84]:


maska = array([True, False, True, False, False])
B[maska]


# In[85]:


maska = array([1,0,1,0,0], dtype=bool)
B[maska]


# Zanimljiviji primjer:

# In[86]:


x = arange(0, 10, 0.5)
x


# In[87]:


maska = (5 < x) * (x < 7.5)
maska


# In[88]:


x[maska]


# ## Funkcije na nizovima

# In[89]:


indeksi = where(maska)
indeksi


# In[90]:


x[indeksi]


# In[91]:


print(A)
diag(A)


# In[92]:


diag(A, -1)


# In[93]:


v2 = arange(-3,3)
v2


# In[94]:


indeksi_redaka = [1, 3, 5]
v2[indeksi_redaka]


# In[95]:


v2.take(indeksi_redaka)


# U sljedećem primjeru `take` djeluje na listu, a izlaz je `array`: 

# In[96]:


take([-3, -2, -1,  0,  1,  2], indeksi_redaka)


# Funkcija `choose`:

# In[97]:


koji = [1, 0, 1, 0]
izbori = [[-1,-2,-3,-4], [5,4,3,2]]

choose(koji, izbori)


# Što radi ova funkcija?

# ## Vektorizacija koda
# 
# Što je više operacija s nizovima, to će kod generalno biti brži.

# In[98]:


v1 = arange(0, 5)


# In[99]:


v1 * 2


# In[100]:


v1 + 2


# In[101]:


print(A)
A * 2, A + 2


# Defaultne operacije na nizovima su **uvijek** definirane po elementima.

# In[102]:


A * A


# In[103]:


v1 * v1


# In[104]:


A.shape, v1.shape


# In[122]:


print(A,v1)
A * v1


# Kako doći do standardnog umnoška?

# In[123]:


dot(A, A)


# In[124]:


A @ A # nova operacija definirana u Python-u 3.5+


# In[125]:


matmul(A,A) # @ je zapravo pokrata za matmul, dot i matmul nisu iste operacije (poklapaju se na 1D i 2D nizovima)


# In[126]:


dot(A, v1)


# In[127]:


A @ v1


# In[128]:


v1 @ v1 # analogno dot(v1, v1)


# Matrice mogu biti i višedimenzionalne

# In[129]:


a = random.rand(8,13,13)
b = random.rand(8,13,13)
matmul(a, b).shape


# Postoji i tip `matrix`. Kod njega operacije `+, -, *` se ponašaju onako kako smo navikli.

# In[130]:


M = matrix(A)
v = matrix(v1).T # da bi dobili stupčasti vektor


# In[131]:


v


# In[132]:


M*M


# In[133]:


M*v


# In[134]:


# skalarni produkt
v.T * v


# In[135]:


v + M*v


# Naravno, dimenzije trebaju biti kompatibilne.

# In[136]:


v = matrix([1,2,3,4,5,6]).T


# In[137]:


shape(M), shape(v)


# In[138]:


M * v


# Još neke funkcije: `inner`, `outer`, `cross`, `kron`, `tensordot`.

# In[139]:


C = matrix([[1j, 2j], [3j, 4j]])
C


# In[140]:


conjugate(C)


# Adjungiranje:

# In[141]:


C.H


# Za izvlačenje realnog, odnosno imaginarnog dijela: `real` i `imag`:

# In[142]:


real(C) # isto što i C.real


# In[143]:


imag(C) # isto što i C.imag


# In[144]:


angle(C+1) # u MATLAB-u je to funkcija arg, dakle argument (faza) kompleksnog broja


# In[145]:


abs(C)


# In[146]:


from numpy.linalg import inv, det
inv(C) # isto što i C.I 


# In[147]:


C.I * C


# In[148]:


det(C)


# In[149]:


det(C.I)


# ## Izvlačenje osnovih informacija iz nizova

# In[150]:


# u stockholm_td_adj.dat su podaci o vremenu za Stockholm
dataStockholm = genfromtxt('stockholm_td_adj.dat')
dataStockholm.shape


# In[151]:


# temperatura se nalazi u 4. stupcu (znači stupcu broj 3)
mean(dataStockholm[:,3])


# Prosječna dnevna temperatura u Stockholmu u zadnjiih 200 godina je bila 6.2 C.

# In[152]:


std(dataStockholm[:,3]), var(dataStockholm[:,3])


# In[153]:


dataStockholm[:,3].min()


# In[154]:


dataStockholm[:,3].max()


# In[155]:


d = arange(0, 10)
d


# In[156]:


sum(d)


# In[157]:


prod(d+1)


# In[158]:


# kumulativa suma
cumsum(d)


# In[159]:


# kumulativan produkt
cumprod(d+1)


# In[160]:


# isto što i: diag(A).sum()
trace(A)


# Naravno, sve ove operacije možemo raditi na dijelovima nizova.

# In[161]:


get_ipython().system('head -n 3 stockholm_td_adj.dat')


# Format je: godina, mjesec, dan, prosječna dnevna temperatura, najniža, najviša, lokacija.
# 
# Recimo da nas zanimaju samo temperature u veljači.

# In[162]:


# mjeseci su 1.,..., 12.
unique(dataStockholm[:,1])


# In[163]:


maska_velj = dataStockholm[:,1] == 2


# In[164]:


mean(dataStockholm[maska_velj,3])


# Sada nije problem doći do histograma za prosječne mjesečne temperature u par redaka.

# In[165]:


mjeseci = arange(1,13)
mjeseci_prosjek = [mean(dataStockholm[dataStockholm[:,1] == mjesec, 3]) for mjesec in mjeseci]

from pylab import *
get_ipython().run_line_magic('matplotlib', 'inline')
fig, ax = subplots()
ax.bar(mjeseci, mjeseci_prosjek)
ax.set_xlabel("Mjesec")
ax.set_ylabel("Prosj. mj. temp.");


# ## Rad s višedimenzionalnim podacima

# In[166]:


m = rand(3,3)
m


# In[167]:


m.max()


# In[168]:


# max u svakom stupcu
m.max(axis=0)


# In[169]:


# max u svakom retku
m.max(axis=1)


# Oblik niza se može promijeniti bez da se dira memorija, dakle mogu se primijenjivati i na veliku količinu podataka.

# In[170]:


A


# In[171]:


n, m = A.shape


# In[172]:


B = A.reshape((1,n*m))
B


# In[173]:


B[0,0:5] = 5 # promijenili smo B
B


# In[174]:


A # a time smo promijenili i A


# Funkcija `flatten` radi kopiju.

# In[175]:


B = A.flatten()
B


# In[176]:


B[0:5] = 10
B


# In[177]:


A # A je sad ostao isti


# In[178]:


v = array([1,2,3])


# In[179]:


shape(v)


# In[180]:


# pretvorimo v u matricu
v[:, newaxis]


# In[181]:


v[:,newaxis].shape


# In[182]:


v[newaxis,:].shape


# In[183]:


a = array([[1, 2], [3, 4]])


# In[184]:


# ponovimo svaki element tri puta
repeat(a, 3)


# In[185]:


tile(a, 3)


# In[186]:


b = array([[5, 6]])


# In[187]:


concatenate((a, b), axis=0)


# In[188]:


concatenate((a, b.T), axis=1)


# In[189]:


vstack((a,b))


# In[190]:


hstack((a,b.T))


# ### Kopiranje nizova

# In[191]:


A = array([[1, 2], [3, 4]])
A


# In[192]:


# B je isto što i A (bez kopiranja podataka)
B = A


# Ako želimo napraviti novu kopiju, koristimo funkciju `copy`:

# In[193]:


B = copy(A)


# In[194]:


v = array([1,2,3,4])

for element in v:
    print (element)


# In[195]:


M = array([[1,2], [3,4]])

for row in M:
    print ("redak {}".format(row))
    
    for element in row:
        print (element)


# Funkcija `enumerate` nam daje i element i njegov indeks:

# In[196]:


for row_idx, row in enumerate(M):
    print ("indeks retka {} redak {}".format(row_idx, row))
    
    for col_idx, element in enumerate(row):
        print ("col_idx {} element {}".format(col_idx, element))
        M[row_idx, col_idx] = element ** 2


# ## Vektorizacija funkcija

# In[197]:


def Theta(x):
    """
    Sklarna verzija step funkcije.
    """
    if x >= 0:
        return 1
    else:
        return 0


# In[198]:


Theta(array([-3,-2,-1,0,1,2,3]))


# In[199]:


Theta_vec = vectorize(Theta)


# In[200]:


Theta_vec(array([-3,-2,-1,0,1,2,3]))


# To smo mogli napraviti i ručno.

# In[201]:


def Theta(x):
    """
    Vektorska verzija step funkcije.
    """
    return 1 * (x >= 0)


# In[202]:


Theta(array([-3,-2,-1,0,1,2,3]))


# In[203]:


# radi naravno i za skalare
Theta(-1.2), Theta(2.6)


# In[204]:


M


# In[205]:


if (M > 5).any():
    print ("barem jedan element iz M je veći od 5")
else:
    print ("svi elementi iz M su manji ili jednaki od 5")


# In[206]:


if (M > 5).all():
    print ("svi elementi iz M su veći od 5")
else:
    print ("barem jedan element je manji ili jednak od 5")


# Eksplicitno pretvaranje podataka. Uvijek stvara novi niz.

# In[207]:


M.dtype


# In[208]:


M2 = M.astype(float)
M2


# In[209]:


M2.dtype


# In[210]:


M3 = M.astype(bool)
M3


# In[211]:


from verzije import *
from IPython.display import HTML
HTML(print_sysinfo()+info_packages('numpy,matplotlib'))


# # Zadaci za vježbu
# 
# - Matrica $A$ reda 5 je dobijena s `A=random.randn(5,5)`. Izvucite $2\times 2$ matricu $B$ iz gornjeg desnog kuta matrice $A$. Pomnožite matricu $B$ s nekom drugom $2\times 2$ matricom te ubacite tu matricu u gornji lijevi kut matrice $A$.
# - Dan je kod
# ```
# x = linspace(0, 1, 3)
# # y = 2*x + 1:
# y=x; y*=2; y+=1
# # z = 4*x - 4:
# z=x; z*=4; z-=4 
# print (x, y, z)
# ```
# Izvršite ovaj kod. Objasnite zašto `x`, `y` i `z` imaju iste vrijednosti.

# - Vektorizirajte kod iz 3. zadatka za vježbu iz prvog predavanja. Vektorizirana funkcija treba za argumente primati funkciju `f` i broj `n`
#  te dva niza `a,b` a vraćati niz brojeva koji odgovaraju aproksimaciji integrala $\int_{a_i}^{b_i} f(x) dx$ trapeznom formulom.

# - Dana je funkcija
# $$
# f(x)=\frac{n}{n+1}\begin{cases}
# (1/2)^{1+1/n}-(1/2-x)^{1+1/n}, & 0\le x\le 1/2,\\
# (1/2)^{1+1/n}-(x-1/2)^{1+1/n}, & 1/2 < x\le 1.
# \end{cases}
# $$
# Ovdje je $n$ realan broj, $0 < n \le 1$. 
# Napišite vektoriziranu funkciju za računanje $f(x)$ na skupu $m$ ekvidistribuiranih brojeva između 0 i 1 (dakle, bez korištenja petlji).
# -  Neka su $x_1,\ldots,x_n$ uniformno distribuirani slučajni brojevi između $a$ i $b$. Tada se $\int_a^b f(x)\mathrm{d}\,x$ može aproksimirati s 
# $\frac{b-a}{n} \sum_{i=1}^n f(x_i)$, postupak koji se obično zove _Monte Carlo integracija_.
# Napišite funkciju koja kao ulazne varijable prima $a$, $b$, $f$ te $n$ s predefiniranom vrijednošću 1000, a vraća rezultat Monte Carlo integracije.
# I u ovom zadatku se ne smiju koristiti petlje.

# In[ ]: