import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
X = np.random.random((1000, 3))

def parse_cap(s):
    return cap.stdout.split(' ')[-4:-2]
Benchs = ['python\nloop', 'numpy\nbroadc.', 'sklearn', 'fortran/\nf2py', 'scipy', 'cython', 'numba','parakeet']
Benchmarks = dict((bench,0) for bench in Benchs)

np.show_config()

%%capture cap

def pairwise_numpy(X):
    return np.sqrt(((X[:, None, :] - X) ** 2).sum(-1))
%timeit pairwise_numpy(X)

Benchmarks['numpy\nbroadc.'] = parse_cap(cap)
print cap.stdout

def pairwise_python(X):
    M = X.shape[0]
    N = X.shape[1]
    D = np.empty((M, M), dtype=np.float)
    for i in range(M):
        for j in range(M):
            d = 0.0
            for k in range(N):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = np.sqrt(d)
    return D

%%capture cap
%timeit pairwise_python(X)

Benchmarks['python\nloop'] = parse_cap(cap)
print cap.stdout

%%capture cap

import numba
from numba import double, jit,autojit

@jit
def pairwise_numba(X):
    M = X.shape[0]
    N = X.shape[1]
    D = np.empty((M, M))#, dtype=np.float)
    for i in range(M):
        for j in range(M):
            d = 0.0
            for k in range(N):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = np.sqrt(d)
    return D

%timeit pairwise_numba(X)

Benchmarks['numba'] = parse_cap(cap)
print numba.__version__, cap.stdout

%load_ext cythonmagic

%%cython

import numpy as np
cimport cython
from libc.math cimport sqrt

@cython.boundscheck(False)
@cython.wraparound(False)
def pairwise_cython(double[:, ::1] X):
    cdef int M = X.shape[0]
    cdef int N = X.shape[1]
    cdef double tmp, d
    cdef double[:, ::1] D = np.empty((M, M), dtype=np.float64)
    for i in range(M):
        for j in range(M):
            d = 0.0
            for k in range(N):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = sqrt(d)
    return np.asarray(D)

%%capture cap
%timeit pairwise_cython(X)

Benchmarks['cython'] = parse_cap(cap)
print cap.stdout

%%capture cap

from parakeet import jit

pairwise_parakeet = jit(pairwise_python)

%timeit pairwise_parakeet(X)

Benchmarks['parakeet'] = parse_cap(cap)
print cap.stdout

%%file pairwise_fort.f

      subroutine pairwise_fort(X,D,m,n)
          integer :: n,m
          double precision, intent(in) :: X(m,n)
          double precision, intent(out) :: D(m,m) 
          integer :: i,j,k
          double precision :: r 
          do i = 1,m 
              do j = 1,m 
                  r = 0
                  do k = 1,n 
                      r = r + (X(i,k) - X(j,k)) * (X(i,k) - X(j,k)) 
                  end do 
                  D(i,j) = sqrt(r) 
              end do 
          end do 
      end subroutine pairwise_fort

!f2py -c --help-fcompiler

!f2py -c --help-compiler

# Compile the Fortran with f2py.
# We'll direct the output into /dev/null so it doesn't fill the screen
#!f2py -c --compiler=intel pairwise_fort.f -m pairwise_fort > /dev/null
!f2py -c --fcompiler=gnu95 pairwise_fort.f -m pairwise_fort > /dev/null

%%capture cap

from pairwise_fort import pairwise_fort
XF = np.asarray(X, order='F')
%timeit pairwise_fort(XF)

Benchmarks['fortran/\nf2py'] = parse_cap(cap)
print cap.stdout

%%capture cap

from scipy.spatial.distance import cdist
%timeit cdist(X, X)

Benchmarks['scipy'] = parse_cap(cap)
print cap.stdout

%%capture cap

from sklearn.metrics import euclidean_distances
%timeit euclidean_distances(X, X)

Benchmarks['sklearn'] = parse_cap(cap)
print cap.stdout

%matplotlib inline

labels = Benchmarks.keys()
timings = [np.timedelta64(int(float(value[0])),str(value[1]))/np.timedelta64(1, 's') for value in Benchmarks.values()]
x = np.arange(len(labels))
ax = plt.axes(xticks=x, yscale='log')
ax.bar(x - 0.3, timings, width=0.6, alpha=0.4, bottom=1E-6)
ax.grid()
ax.set_xlim(-0.5, len(labels) - 0.5)
ax.set_ylim(0.5*min(timings), 1.1*max(timings))
ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda i, loc: labels[int(i)]))
ax.set_ylabel('time (s)')
ax.set_title("Pairwise Distance Timings")

labels = [val for val in Benchmarks.iterkeys() if val != 'python\nloop']
timings = [np.timedelta64(int(float(value[0])),str(value[1]))/np.timedelta64(1, 's')\
           for k, value in Benchmarks.items() if k != 'python\nloop']
x = np.arange(len(labels))
ax = plt.axes(xticks=x, yscale='log')
ax.bar(x - 0.3, timings, width=0.6, alpha=0.4, bottom=1E-6)
ax.grid()
ax.set_xlim(-0.5, len(labels) - 0.5)
ax.set_ylim(0.5*min(timings), 1.1*max(timings))
ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda i, loc: labels[int(i)]))
ax.set_ylabel('time (s)')
ax.set_title("Pairwise Distance Timings")

#Restart notebook!

S0 = 100
K = 100
T = 1
vol = 0.3
r = 0.03
q = 0
optype = -1

from numpy import log, sqrt, zeros, exp, arange
from scipy.special import erf

def parse_cap(s):
    return cap.stdout.split(' ')[-4:-2]

benchs = ['python','numba','cython','matlab']
BS_bench = dict(zip(benchs,zeros(4)))

def std_norm_cdf(x):
    return 0.5*(1+erf(x/sqrt(2.0)))

def black_scholes(s, k, t, v, rf, div, cp):
    """Price an option using the Black-Scholes model.
    
    s : initial stock price
    k : strike price
    t : expiration time
    v : volatility
    rf : risk-free rate
    div : dividend
    cp : +1/-1 for call/put
    """
    d1 = (log(s/k)+(rf-div+0.5*pow(v,2))*t)/(v*sqrt(t))
    d2 = d1 - v*sqrt(t)
    optprice = cp*s*exp(-div*t)*std_norm_cdf(cp*d1) - \
            cp*k*exp(-rf*t)*std_norm_cdf(cp*d2)
    return optprice

%%capture cap
%timeit black_scholes(S0, K, T, vol, r, q, optype)

BS_bench['python'] = parse_cap(cap)
print cap.stdout

from numba import double
from numba.decorators import jit, autojit
import math

@autojit
def std_norm_cdf_numba(x):
    return 0.5*(1+math.erf(x/math.sqrt(2.0)))

@autojit
def black_scholes_numba(s, k, t, v, rf, div, cp):
    """Price an option using the Black-Scholes model.
    
    s : initial stock price
    k : strike price
    t : expiration time
    v : volatility
    rf : risk-free rate
    div : dividend
    cp : +1/-1 for call/put
    """
    d1 = (math.log(s/k)+(rf-div+0.5*v*v)*t)/(v*math.sqrt(t))
    d2 = d1 - v*math.sqrt(t)
    optprice = cp*s*math.exp(-div*t)*std_norm_cdf_numba(cp*d1) - \
            cp*k*math.exp(-rf*t)*std_norm_cdf_numba(cp*d2)
    return optprice

%%capture cap
%timeit black_scholes_numba(S0, K, T, vol, r, q, optype)

BS_bench['numba'] = parse_cap(cap)
print cap.stdout

%load_ext cythonmagic

%%cython
cimport cython
from libc.math cimport exp, sqrt, pow, log, erf

@cython.cdivision(True)
cdef double std_norm_cdf(double x) nogil:
    return 0.5*(1+erf(x/sqrt(2.0)))

@cython.cdivision(True)
def black_scholes(double s, double k, double t, double v,
                 double rf, double div, double cp):
    """Price an option using the Black-Scholes model.
    
    s : initial stock price
    k : strike price
    t : expiration time
    v : volatility
    rf : risk-free rate
    div : dividend
    cp : +1/-1 for call/put
    """
    cdef double d1, d2, optprice
    with nogil:
        d1 = (log(s/k)+(rf-div+0.5*pow(v,2))*t)/(v*sqrt(t))
        d2 = d1 - v*sqrt(t)
        optprice = cp*s*exp(-div*t)*std_norm_cdf(cp*d1) - \
            cp*k*exp(-rf*t)*std_norm_cdf(cp*d2)
    return optprice

%%capture cap
%timeit black_scholes(S0, K, T, vol, r, q, optype)

BS_bench['cython'] = parse_cap(cap)
print cap.stdout

!cat ../std_norm_cdf.m

!cat ../black_scholes.m

%matplotlib inline

BS_bench['matlab'] = [u'55',u'us']
labels= BS_bench.keys()
timings = [np.timedelta64(int(float(value[0])),str(value[1]))/np.timedelta64(1, 's') for value in BS_bench.values()]
x = arange(len(labels))
ax = plt.axes(xticks=x, yscale='log')
ax.bar(x - 0.3, timings, width=0.6, alpha=0.4, bottom=1E-6)
ax.grid()
ax.set_xlim(-0.5, len(labels) - 0.5)
ax.set_ylim(min(timings), 1.1*max(timings))
ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda i, loc: labels[int(i)]))
ax.set_ylabel('time ($s$)')
ax.set_title("Black-Scholes Timings")

%load_ext version_information
%version_information numpy,numba,parakeet,cython,matplotlib,scipy,sklearn

from IPython.core.display import HTML
def css_styling():
    styles = open("./styles/custom.css", "r").read()
    return HTML(styles)
css_styling()