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

# # Ubrzavanje Pythona

# In[1]:


from IPython.display import Image
Image("https://raw.github.com/jrjohansson/scientific-python-lectures/master/images/optimizing-what.png")


# Katkada korištenje Numpy-a ne daje dovoljno ubrzanje, ili je teško/nespretno vektorizirati kod. Tada postoji više opcija
# 
# 1. spori dio koda (koji dio koda je spor možemo saznati profiliranjem koda, npr. pomoću IPythonove magične funkcije `%prun`) možemo napisati npr. u C-u. U Pythonu je taj proces prilično bezbolan.
# 2. spori dio koda možemo implementirati u [Cythonu](http://cython.org/), koji je proširenje Pythona
# 3. možemo ubrzati kod korištenjem specijaliziranih kompajlera koji optimiziraju strojni kod

# Pod 3., postoji cijeli niz kompajlera za Python: [PyPy](http://pypy.org/), [Nuitka](http://nuitka.net/), [shedskin](https://code.google.com/p/shedskin/), [parakeet](https://github.com/iskandr/parakeet), [psyco](http://psyco.sourceforge.net/), [Theano](http://deeplearning.net/software/theano/) i [Numba](http://numba.pydata.org/). 
# 
# Mi ćemo se ovdje pozabaviti s Cythonom i Numbom. U oba slučaja je od najveće važnosti precizirati vrste podataka (tzv. anotacija) jer se onda strojni kod može dobro optimizirati.

# ## Cython

# In[2]:


# olakšava rad u Cythonu u IPythonu
get_ipython().run_line_magic('load_ext', 'Cython')


# In[3]:


a=10
b=20


# `%%cython_inline` magična funkcija koristi `Cython.inline` da bi se kompajlirao Cython izraz. `return` služi za slanje izlaza.

# In[4]:


get_ipython().run_cell_magic('cython_inline', '', 'return a+b\n')


# `%%cython_pyximport` magična funkcija služi da se unese proizvoljan Cython kod u IPython notebook ćeliju. Taj kod se sprema u `.pyx` u radnom direktoriju i onda importira koristeći `pyximport`. Moramo specificirati ime modula (u donjem slučaju `foo`). Svi objekti modula se automatski importiraju.

# In[5]:


get_ipython().run_cell_magic('cython_pyximport', 'foo', 'def f(x):\n    return 4.0*x\n')


# In[6]:


f(10)


# U dosadašnjim primjerima Cython kod se nije razlikovao od Python koda. U sljedećem primjeru ćemo vidjeti neka od proširenja Pythona koja nudi Cython.
# 
# Magična funkcija `%%cython` je slična funkciji `%%cython_pyximport`, ali ne traži ime modula te sve datoteke sprema u privremeni direktorij.
# 
# Jedan od načina računanja broja $\pi$ je korištenjem formule
# $$ \frac{\pi}{4} = \int_0^1 \sqrt{1-x^2}\,\mathrm{d}x.$$  
# A integral možemo aproksimirati pomoću trapezne formule, što nam daje
# $$\frac{\pi}{4}\approx \frac{1}{n}\left(\frac{1}{2}+\sum_{i=1}^n\sqrt{1-\left(\frac{i}{n}\right)^2} \right) $$
# Krenimo od običnog Pythona. 

# In[7]:


from math import sqrt
def funkcija(x):
    return sqrt(1-x**2)

def integral4pi(n):
    korak = 1.0/n
    rez = (funkcija(0)+funkcija(1))/2
    for i in range(n):
        rez += funkcija(i*korak)
    return 4*rez*korak


# In[8]:


approx=integral4pi(10**7)
print ('pi={}'.format(approx))


# In[9]:


get_ipython().run_line_magic('timeit', 'integral4pi(10**7)')


# ### Cython verzija
# Ovdje je
# 
# - `cimport`: ekvivalent importa, ali možemo učitavati i iz C ili C++ biblioteka
# - `cdef`: ekvivalent za def, gdje definiramo (kao u C-u) tipove podataka; također služi za deklariranje tipa varijabli

# In[10]:


get_ipython().run_cell_magic('cython', '', '\ncimport cython\nfrom libc.math cimport sqrt\ncdef double cy_funkcija(double x):\n    return sqrt(1-x**2)\n\ndef cy_integral4pi(int n):\n    cdef double korak, rez\n    cdef int i\n    korak = 1.0/n\n    rez = (cy_funkcija(0)+cy_funkcija(1))/2\n    for i in range(n):\n        rez += cy_funkcija(i*korak)\n    return 4*rez*korak\n')


# In[11]:


cy_approx = cy_integral4pi(10**7)
print ('pi={}'.format(cy_approx))


# In[12]:


get_ipython().run_line_magic('timeit', 'cy_integral4pi(10**7)')


# ### Složeniji primjer
# Ovdje je
# 
# - `@`: Pythonova sintaksa za tzv. [dekoratore](http://en.wikipedia.org/wiki/Python_syntax_and_semantics#Decorators)
# - `nogil`: GIL je skraćenica od _global interpreter lock_, koji spriječava simultano izvršavanje koda u više niti (eng. threads), ovdje s `nogil` deklariramo da je sigurno pozvati funkciju bez GIL-a
# - `with nogil`: odgovarajući dio koda ne koristi GIL

# In[13]:


get_ipython().run_cell_magic('cython', '', 'cimport cython\nfrom libc.math cimport exp, sqrt, pow, log, erf\n@cython.cdivision(True)\ncdef double std_norm_cdf(double x) nogil:\n    return 0.5*(1+erf(x/sqrt(2.0)))\n@cython.cdivision(True)\ndef black_scholes(double s, double k, double t, double v, double rf, double div, double cp):\n    """s : početna cijena dionice\n    k : fiksirana cijena (strike price)\n    t : vrijeme \n    v : volatilnost\n    rf : bezrizična kamata\n    div : dividenda\n    cp : put/call paritet\n    """\n    cdef double d1, d2, optprice\n    with nogil:\n        d1 = (log(s/k)+(rf-div+0.5*pow(v,2))*t)/(v*sqrt(t))\n        d2 = d1 - v*sqrt(t)\n        optprice = cp*s*exp(-div*t)*std_norm_cdf(cp*d1) - cp*k*exp(-rf*t)*std_norm_cdf(cp*d2)\n    return optprice\n')


# In[14]:


black_scholes(100.0, 100.0, 1.0, 0.3, 0.03, 0.0, -1)


# In[15]:


get_ipython().run_line_magic('timeit', 'black_scholes(100.0, 100.0, 1.0, 0.3, 0.03, 0.0, -1)')


# Cython omogućava linkanje dodatnih biblioteka pri kompajliranju, u ovom slučaju standardne matematičke biblioteke

# In[33]:


get_ipython().run_cell_magic('cython', '-lm', "from libc.math cimport sin\nprint ('sin(1)=', sin(1))\n")


# ### Još jedan primjer, gdje koristimo numpy nizove
# 
# Računamo kumulativnu sumu niza

# In[34]:


import numpy as np
def py_dcumsum(a):
    b = np.empty_like(a)
    b[0] = a[0]
    for n in range(1,len(a)):
        b[n] = b[n-1]+a[n]
    return b


# In[35]:


a = np.random.rand(100000)
b = np.empty_like(a)


# In[36]:


get_ipython().run_cell_magic('cython', '', '\ncimport numpy\n\ndef cy_dcumsum2(numpy.ndarray[numpy.float64_t, ndim=1] a, numpy.ndarray[numpy.float64_t, ndim=1] b):\n    cdef int i, n = len(a)\n    b[0] = a[0]\n    for i from 1 <= i < n:\n        b[i] = b[i-1] + a[i]\n    return b\n')


# In[37]:


get_ipython().run_line_magic('timeit', 'cy_dcumsum2(a,b)')


# In[38]:


get_ipython().run_line_magic('timeit', 'py_dcumsum(a)')


# ## Numba
# 
# Numba je just-in-time kompajler koji korisiti LLVM. Korištenje numbe je zapravo dosta jednostavno.

# In[39]:


# iz numbe učitavamo jit kompajler
from numba import jit

from numpy import arange


# Funkcija `sum2d(arr)`će biti optimizirana, a cijela procedura se svela na korištenje jednog dekoratora. 

# In[40]:


# @jit je dekorator
@jit
def sum2d(arr):
    M, N = arr.shape
    result = 0.0
    for i in range(M):
        for j in range(N):
            result += arr[i,j]
    return result


# In[41]:


def py_sum2d(arr):
    M, N = arr.shape
    result = 0.0
    for i in range(M):
        for j in range(N):
            result += arr[i,j]
    return result


# In[42]:


a = arange(9).reshape(3,3)


# In[43]:


get_ipython().run_line_magic('timeit', 'sum2d(a)')


# In[44]:


get_ipython().run_line_magic('timeit', 'py_sum2d(a)')


# ### Kompliciraniji primjer

# In[45]:


import numpy
def filter2d(image, filt):
    M, N = image.shape
    Mf, Nf = filt.shape
    Mf2 = Mf // 2
    Nf2 = Nf // 2
    result = numpy.zeros_like(image)
    for i in range(Mf2, M - Mf2):
        for j in range(Nf2, N - Nf2):
            num = 0.0
            for ii in range(Mf):
                for jj in range(Nf):
                    num += (filt[Mf-1-ii, Nf-1-jj] * image[i-Mf2+ii, j-Nf2+jj])
            result[i, j] = num
    return result
from numba import double
fastfilter_2d = jit(double[:,:](double[:,:], double[:,:]))(filter2d)


# In[46]:


image = numpy.random.random((100, 100))
filt = numpy.random.random((10, 10))


# In[47]:


get_ipython().run_line_magic('timeit', 'fastfilter_2d(image, filt)')


# In[48]:


get_ipython().run_line_magic('timeit', 'filter2d(image, filt)')


# ### Paralelizacija
# Drugi način ubrzavanja izvršavanja programa je pomoću paralelizacije. To je opsežna tema sama po sebi te u nju nećemo ulaziti.
# Popularni način paralelizacije je npr. pomoću [Apache Sparka](https://spark.apache.org/) ili pomoću [Dask](http://dask.pydata.org/) biblioteke.

# In[49]:


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