5
\$\begingroup\$

I've developed some code in Python to apply a harmonic analysis of timeseries (for satellite imagery data). It's based on this, but then I would like to optimize the performance. So instead of a single timeseries as input I've an array of 10000 timeseries as input.

Even though I've tried to make use of some nice functions in NumPy that deals with big multidimensional arrays, I'm sure I just have touched the surface as it comes to the actual capabilities.

import numpy as np

# Computing diagonal for each row of a 2d array. See: http://stackoverflow.com/q/27214027/2459096
def makediag3d(M):
    b = np.zeros((M.shape[0], M.shape[1]*M.shape[1]))
    b[:, ::M.shape[1]+1] = M
    return b.reshape(M.shape[0], M.shape[1], M.shape[1])

# Function to apply the Harmonic analysis of time series applied to arrays    
def HANTS(ni,y,nf=3,HiLo='Hi',low=0.,high=255,fet=5,delta=0.1):
    """
    ni    = nr. of images (total number of actual samples of the time series)
    nb    = length of the base period, measured in virtual samples 
            (days, dekads, months, etc.)        
    nf    = number of frequencies to be considered above the zero frequency
    y     = array of input sample values (e.g. NDVI values)
    ts    = array of size ni of time sample indicators 
            (indicates virtual sample number relative to the base period); 
            numbers in array ts maybe greater than nb
            If no aux file is used (no time samples), we assume ts(i)= i, 
            where i=1, ..., ni         
    HiLo  = 2-character string indicating rejection of high or low outliers
            select from 'Hi', 'Lo' or 'None'    
    low   = valid range minimum
    high  = valid range maximum (values outside the valid range are rejeced
            right away)    
    fet   = fit error tolerance (points deviating more than fet from curve 
            fit are rejected)
    dod   = degree of overdeterminedness (iteration stops if number of 
            points reaches the minimum required for curve fitting, plus 
            dod). This is a safety measure            
    delta = small positive number (e.g. 0.1) to suppress high amplitudes                  
    """

    # define some parameters
    nb= ni # 
    ts=np.arange(ni)
    dod=1# (2*nf-1)

    # create empty arrays to fill
    mat = np.zeros(shape=(min(2*nf+1,ni),ni))

    yr = np.zeros(shape=(y.shape[0],ni))

    # check which setting to set for outlier filtering
    if HiLo == 'Hi':
        sHiLo = -1
    elif HiLo == 'Lo':
        sHiLo = 1
    else:
        sHiLo = 0

    # initiate parameters
    nr = min(2*nf+1,ni) # number of 2*+1 frequecies, or number of input images    
    noutmax = ni-nr-dod # number of input images - number of 2*+1 frequencies - degree of overdeterminedness        
    mat[0,:] = 1    
    ang = 2*np.pi*np.arange(nb)/nb
    cs = np.cos(ang)
    sn = np.sin(ang)        

    # create some standard sinus and cosinus functions and put in matrix
    i = np.arange(1,nf+1)
    for j in np.arange(ni):    
        index = np.mod(i*ts[j],nb)
        mat[2 * i-1,j] =cs.take(index)
        mat[2 * i ,j] = sn.take(index)

    # repeat the mat array over the number of arrays in y
    # and create arrays with ones with shape y where high and low values are set to 0
    mat = np.tile(mat[None].T, (1,y.shape[0])).T 
    p = np.ones_like(y) 
    p[(low >= y) | (y > high)] = 0
    nout = np.sum(p==0,axis=-1) # count the outliers for each timeseries      

    # prepare for while loop
    ready = np.zeros((y.shape[0]), dtype=bool) # all timeseries set to false
    a = np.arange(ni)
    it = np.nditer(a)

    while ((not it.finished) & (not ready.all())):

        #print '--------*-*-*-*',it.value, '*-*-*-*--------'        
        # multipy outliers with timeseries
        za = np.einsum('ijk,ik->ij', mat, p*y)

        # multiply mat with the multiplication of multiply diagonal of p with transpose of mat
        diag = makediag3d(p)
        A = np.einsum('ajk,aki->aji',mat,np.einsum('aij,jka->ajk', diag, mat.T))         
        # add delta to supress high amplitudes but not for [0,0]
        A = A + np.tile(np.diag(np.ones(nr))[None].T,(1,y.shape[0])).T * delta        
        A[:,0,0] = A[:,0,0] - delta 

        # solve linear matrix equation and define reconstructed timeseries
        zr = np.linalg.solve(A, za)
        yr = np.einsum('ijk,kj->ki',mat.T, zr)

        # calculate error and sort err by index
        err = p * (sHiLo * (yr - y))
        rankVec = np.argsort(err, axis=1, )

        # select maximum error and compute new ready status        
        maxerr = np.diag(err.take(rankVec[:,ni-1], axis=-1))
        ready = (maxerr <= fet) | (nout == noutmax) 

        # if ready is still false 
        if (not all(ready)):
            i = ni # i is number of input images
            j = rankVec.take(i-1, axis=-1)

            p.T[j.T, np.indices(j.shape)] = p.T[j.T, np.indices(j.shape)]*ready.astype(int)#*check
            nout += 1
            i -= 1

        it.iternext()
    return yr

# Compute semi-random time series array with numb standing for number of timeseries
def array_in(numb):
    y= np.array([5.0,  2.0,  10.0, 12.0, 18.0, 23.0, 27.0, 40.0, 60.0, 70.0, 90.0,160.0,190.0, 
                 210.0,104.0,90.0,170.0, 50.0,120.0, 60.0, 40.0, 30.0, 28.0, 24.0, 15.0, 10.0 ])
    y = np.tile(y[None].T, (1,numb)).T
    kl = (np.random.randint(2, size=(numb, 26)) * 
          np.random.randint(2, size=(numb, 26)) + 1)
    kl[kl==2] = 0
    y=y*kl
    return y

By applying this on an array of 10000 I get the following output:

y = array_in(10000)
%timeit HANTS(ni=26, y=y, nf=3, HiLo='Lo')
1 loops, best of 3: 10.5 s per loop

Which gives a possible output like this:

and some possible output

Even though it works I assume it's all in all a little bit on the slow side.

I've tested this code in both IPython and Python version 2.7 and 3.4 with NumPy 1.8 and 1.9. Other versions seems to have problems.

\$\endgroup\$
4
\$\begingroup\$

Well, I'm going to take a break and may not get back to it too soon, but I created https://github.com/pckujawa/harmonic_analysis_of_time_series and an accompanying "journal" to show how I profile and (will) optimize the code. For now, you can just look at each commit in order for my advice :). (BTW, I really recommend High Performance Python - O'Reilly Media - it's good all around, but it hits especially on numpy and scientific computing.)

I used to know numpy pretty well, but I've forgotten a lot. Any experts out there could probably suggest better functions to achieve what you want. But without being a domain expert, it's hard to grok this code.

Some highlights

Profiling shows that ndarray.take is the bottleneck. This is happening twice inside a loop and is where about 64% of the time is spent. Next are calls to einsum which take about 20% of the time. So if we can find different ways to express those, we'll see big wins.

*** PROFILER RESULTS ***
HANTS (hants/__init__.py:12)
function called 1 times

         1384 function calls in 11.479 seconds

   Ordered by: cumulative time, internal time, call count
   List reduced from 48 to 40 due to restriction <40>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    1.312    1.312   11.479   11.479 __init__.py:12(HANTS)
       83    7.402    0.089    7.402    0.089 {method 'take' of 'numpy.ndarray' objects}
       64    2.293    0.036    2.293    0.036 {numpy.core.multiarray.einsum}
       16    0.175    0.011    0.257    0.016 __init__.py:6(makediag3d)
       16    0.092    0.006    0.096    0.006 linalg.py:296(solve)
\$\endgroup\$
  • \$\begingroup\$ Very nice overview with the profiler. As ndarray.take is called 83 times I expect it's not only called twice inside this loop. And following wesmckinney.com/blog/?p=215 I guess that .take actually should be able to optimize the performance if it can replace the fancy indexing (line 108, 119, 126 in github.com/pckujawa/harmonic_analysis_of_time_series/blob/…), which I expect being categorized under ndarray.take as well (explaining the high number of ncalls) \$\endgroup\$ – Mattijn Dec 7 '14 at 14:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.