%pylab inline

import numpy as np
import skimage.data as skiD
import skimage.color as clr
from time import time


img = skiD.lena()
img = clr.rgb2grey(img)

###################
p = img[100:200, 400:]

figure()
imshow(img, cmap='gray')
figure()
imshow(p, cmap='gray')

def calcHOGNaive(img, binSize):
    # calculate gradient
    dx, dy = np.gradient(img)
    
    # calculate magnitude of gradient
    mag = np.sqrt(dx ** 2 + dy ** 2).flatten()
    
    # calculate angles and round them into bins
    ## plain angle 0 .. 2pi
    # >>>> angle = np.arctan2(dy, dx) + np.pi)
    ## normalization to 0 .. binSize-0.000000000001
    # >>>> normAngle = angle / (2 * np.pi + np.spacing(10000)) * binSize)
    ## binning and flatting
    # >>>> np.floor(normAngle).flatten()    
    
    ang = (np.floor((np.arctan2(dy, dx) + np.pi) \
            / (2 * np.pi + np.spacing(10000)) * binSize)).flatten()
    
    # initialize histogram
    hog = np.zeros((binSize, ))
    
    # fill histogram
    for i in range(ang.shape[0]):
        hog[ang[i]] += mag[i]
    
    # normalize
    return hog / ang.shape[0]

hog = calcHOGNaive(p, 9)
%timeit calcHOGNaive(p, 9)

%prun -s cumulative calcHOGNaive(p, 9)

import pyTools.sandbox.lineProfExample as lPE

%load_ext line_profiler_ext
%lprun -f lPE.calcHOGNaive lPE.calcHOGNaive(p, 9)

import line_profiler
lP = line_profiler.LineProfiler(lPE.calcHOGNaive)
res = lP.runcall(lPE.calcHOGNaive, p, 9)
lP.print_stats()

#import scipy.weave as weave

def calcHOGWeaveCpp(img, binSize):
    import scipy.weave as weave
    import numpy as np
    dx, dy = np.gradient(img)
    
    mag = np.sqrt(dx ** 2 + dy ** 2).flatten()
    ang = (np.floor((np.arctan2(dy, dx) + np.pi) \
            / (2 * np.pi + np.spacing(10000)) * binSize)).flatten()
        
    hog = np.zeros((binSize, ))
    
#    for i in range(ang.shape[0]):
#        hog[ang[i]] += mag[i]
    code = """
    int pos;
    for (int i = 0; i < Nang[0]; ++i){
        pos = ang[i];
        hog[pos] += mag[i];
    }
    """    
    weave.inline(code, ['ang', 'hog', 'mag'], 
                 extra_compile_args=['-O3'], compiler='gcc')
    
    return hog / ang.shape[0]

hogCpp = calcHOGWeaveCpp(p, 9)
%timeit calcHOGWeaveCpp(p, 9)
print "numpy and cpp version are the same: ", np.allclose(hog, hogCpp)

import pyTools.imgProc.hog as fHog

dx, dy = np.gradient(p)
hogF = fHog.hog.anghist(dx, dy, 9)
%timeit dx, dy = np.gradient(p); fHog.hog.anghist(dx, dy, 9)

print "Naive and Fortran (1) equal: ", np.allclose(hog, hogF)

import pyTools.imgProc.hog as fHog

hogF2 = fHog.hog.anghist2(p, 9)
%timeit fHog.hog.anghist2(p, 9)

print "Naive and Fortran (2) equal: ", np.allclose(hog, hogF2)

from numba.decorators import jit
from numba import float64, int32

@jit(argtypes=[float64[:,:], int32])
def calcHOGNumba(img, binSize):
    # calculate gradient
    d = np.gradient(img)
    
    dx = d[0]
    dy = d[1]
    
    # calculate magnitude of gradient
    mag = np.sqrt(dx ** 2 + dy ** 2).flatten()
    
    # calculate angles and round them into bins   
    
    ang = (np.floor((np.arctan2(dy, dx) + np.pi) \
            / (2 * np.pi + np.spacing(10000)) * binSize)).flatten()
    
    # initialize histogram
    hog = np.zeros((binSize, ))
    
    # fill histogram
    for i in range(ang.shape[0]):
        hog[ang[i]] += mag[i]
    
    # normalize
    return hog / ang.shape[0]

hogN = calcHOGNumba(p, 9)
%timeit calcHOGNumba(p, 9)
print "Naive and Numba equal: ", np.allclose(hog, hogN)

from numba.decorators import jit
from numba import float64, int32

@jit(argtypes=[float64[:], int32[:], float64[:]])
def fillHist(hog, ang, mag): 
    # fill histogram
    for i in range(ang.shape[0]):
        hog[ang[i]] += mag[i]
    
    return hog
        
        
def calcHOGNumba2(img, binSize):
    # calculate gradient
    d = np.gradient(img)
    
    dx = d[0]
    dy = d[1]
    
    # calculate magnitude of gradient
    mag = np.sqrt(dx ** 2 + dy ** 2).flatten()
    
    # calculate angles and round them into bins   
    
    ang = np.int32(np.floor((np.arctan2(dy, dx) + np.pi) \
            / (2 * np.pi + np.spacing(10000)) * binSize)).flatten()
    
    # initialize histogram
    hog = np.zeros((binSize, ))
    
    hog = fillHist(hog, ang, mag)
    
    # normalize
    return hog / ang.shape[0]

hogN2 = calcHOGNumba2(p, 9)
%timeit calcHOGNumba2(p, 9)
print "Naive and Numba equal: ", np.allclose(hog, hogN2)

p = img[100:200, 400:]
dx, dy = np.gradient(p)
hog = calcHOGNaive(p, 9)
hogCpp = calcHOGWeaveCpp(p, 9)
hogF = fHog.hog.anghist(dx, dy, 9)
hogF2 = fHog.hog.anghist2(p, 9)
hogN = calcHOGNumba(p, 9)
hogN2 = calcHOGNumba2(p, 9)

print "all results are equal:", all([np.allclose(hog, hogCpp), 
                                     np.allclose(hog, hogF), 
                                     np.allclose(hog, hogF2), 
                                     np.allclose(hog, hogN), 
                                     np.allclose(hog, hogN2)])

# patch 100x112
p = img[100:200, 400:]
dx, dy = np.gradient(p)
print "naive"
%timeit -n 100 calcHOGNaive(p, 9)
print "\nweave cpp"
%timeit -n 100 calcHOGWeaveCpp(p, 9)
print "\nfortran (1)"
%timeit -n 100 dx, dy = np.gradient(p); fHog.hog.anghist(dx, dy, 9)
print "\nfortran (2)"
%timeit -n 100 fHog.hog.anghist2(p, 9)
print "\nnumba (1)"
%timeit -n 100 calcHOGNumba(p, 9)
print "\nnumba (2)"
%timeit -n 100 calcHOGNumba2(p, 9)

# over entire image
print "naive"
%timeit -n 100 calcHOGNaive(img, 9)
print "\nweave cpp"
%timeit -n 100 calcHOGWeaveCpp(img, 9)
print "\nfortran (1)"
%timeit -n 100 dx, dy = np.gradient(img); fHog.hog.anghist(dx, dy, 9)
print "\nfortran (2)"
%timeit -n 100 fHog.hog.anghist2(img, 9)
print "\nnumba (1)"
%timeit -n 100 calcHOGNumba(img, 9)
print "\nnumba (2)"
%timeit -n 100 calcHOGNumba2(img, 9)

from IPython.parallel import Client
# open clients
rc = Client()
print rc.ids

dview = rc[:]        
lbview = rc.load_balanced_view()

splits = np.split(img, 16)

results = []
for split in splits:
    results += [lbview.apply_async(calcHOGWeaveCpp, split, 9)]

from time import sleep

allReady = False
while not allReady:
    allReady = True
    for ar in results:
        allReady = allReady and ar.ready()
                
    sleep(0.01)

hogs = []

for ar in results:
    # copy data
    data = ar.get()
    hogs += [data]
    # delete data from cluster
    msgId = ar.msg_id
    #~ del lbview.results[msgId]
    del rc.results[msgId]
    del rc.metadata[msgId]
    

# close client to close socket connections
rc.close()

from PyQt4.QtCore import * 

class BaseThread(QThread):
    """ Provides infrastructure for save deleting of thread.
        Check self.exiting in heavy operations frequently and 
        abort exectution of run(), to allow the thread to get 
        deleted.
    """
    def __init__(self):
        QThread.__init__(self)
        self.exiting = False

    def __del__(self):     
        self.exiting = True
        self.wait()
        
class Worker(BaseThread):            
    def __init__(self, func, *args):
        """
        This object will exectute `func` with `args` in a
        separate thread. You can query ready() to check
        if processing finished and get() to get the result.

        Args:
            func (function pointer)
                        function will be called asyncroneously
            args (arguments)
                        arguments for func
        """
        BaseThread.__init__(self)   
        
        self.func = func
        self.args = args        
        self.result = None
        self.isReady = False
    
        self.start()
    
    def run(self):
        self.result = self.func(*self.args)
        self.isReady = True
            
    def ready(self):
        return self.isReady
    
    def get(self):
        if self.isReady:
            return self.result
        else:
            return None

w = Worker(calcHOGWeaveCpp, p, 9)
w.ready()

working = True
while working:
    isReady = w.ready()
    if isReady:
        res = w.get()
        working = False
        
        # clean up, otherwise the object will
        # stay in memory for ever (similar to
        # clients in IPython parallel)
        del w
        
print res

def test(a, *args):
    print a
    print args
    test2(*args)
    
def test2(a,b):
    print a
    print b

test(4,5,6)