Pair, differences and correlations

Question

This has a couple of twists which I personally found rather awkward... All code can be downloaded in a read-to-run format from here

The problem

Consider a classic interview question: "Find given pairs with a specific difference in a given array"

Here I impose the unique conditions:

• The array is non-decreasing
• Consider combinations not permutations arr[i],arr[j] == arr[j],arr[i]
• To make the operation more interesting, calculate the correlation function defined as ((arr[i]-mean_arr)*(arr[j]-mean_arr)).mean()
• Consequently, when no pairs are found the return should be np.nan

Benchmark Naiive Example

This will scale terribly but is easy to read

import itertools
import numpy as np
diff = 0.1                 # difference
arr = np.linspace(0,1,11)  # input array example

def naiiveFn(arr, diff, tol=1e-5):
    # get all combinations of the list with itself
    all_pairs = list(itertools.product(*[list(arr)]*2))
    all_pairs = np.asarray(all_pairs)
    # select all valid pairs
    pairs = all_pairs[np.abs(np.diff(all_pairs, axis=1).ravel() - diff) < tol]
    # calculate the correlation function ((pair1_i - mean)*(pair2_i - mean)).mean()
    result = np.prod(pairs- arr.mean(), axis=1).mean()
    return result

The target is speed speed speed!!!

Test cases

In this section I provide two examples for accuracy and performance testing

Performance

A good benchmark case to challenge performance requires a lot of varied separations but must hold the non-decreasing property.

The following provides a good example, where we want performance to scale well with n

diff = 0.1; n = 2000
perf_arr_small = np.cumsum(np.around(np.random.exponential(diff, n), 1))

and when run for our basic function

%timeit naiiveFn(perf_arr, diff)
1 loop, best of 3: 1.85 s per loop

A more challenging case is n > 20000: I ran naiiveFn at n=20000, went for a coffee break and it was still going when I came back!

Accuracy

Obviously a function is totally useless if if doesn't do what it is supposed to. This will test six cases that I picked to exhibit particular behaviours these are generated by running genTestData() from below

from scipy.special import binom 
import numpy as np # it's beyond me that neither np.binom or math.binom exist

def testFn(results, test, res_pairs=None):
    """Takes in two arrays of your function results and compares to test data
    If you pass the pairs you found to res_pairs it will also neatly display those

    Required Inputs
        results :: list :: list of results from test function
        test :: list :: list of test cases to compare against

    Optional Inputs
        res_pairs :: list of lists/np.arrays :: list of pairs that were found
    """
    outcomes = ["Failed","Passed"]

    print "\nTest outcomes..."
    if res_pairs is None: res_pairs = len(test)*[None]
    for i,(r,t, pairs) in enumerate(zip(results, test, res_pairs)):
        try:
            np.testing.assert_almost_equal(r,t)
            passed = True
        except:
            passed = False
        pr = "  test:{} :: {} :: res: {:7.4f} actual: {:7.4f}".format(i+1, outcomes[passed], r, t)
        if pairs is not None: pr += " pairs: "+" ".join(["{:d}x({:3.1f} {:3.1f})".format(n,i,j) for (i,j),n in Counter(tuple(p) for p in pairs).iteritems()])
        print pr
    pass

def genTestData():
    """Generate test data
        test_cases :: list of 6 test arrays each of length 10
        test_set1  :: the four test cases with 0.1 separation
        test_set2  :: the four test cases with no separation
    """
    n = 10

    # Examples to catch most common errors
    a = np.array([0.1]*10)         # case of everything the same
    b = np.linspace(0.1, 1, 10)    # everything spaced equally
    c = np.array([0.1]*5+[0.2]*5)  # intersection of two repeating segments
    d = np.array([0.1, 0.2, 0.3] + [0.4]*5 + [0.5]*2) # a mash-up
    e = np.array([0.4]*3 + [0.5]*3 + [0.6]*3 + [5]) # series of identicals
    f = np.asarray([0, 0.2, 0.4, 0.6, 0.8] + [0.9]*3 + [1.0]*2) # no match then matches


    # a quick function used a fair bit in the case of equal incrementation
    equalSpacing = lambda seg, mean, sep: np.sum((seg[:seg.size-sep]-mean)*(seg[sep:]-mean))
    nCr52 = binom(5,2)  # ways of choosing n from r where order matters
    nCr32 = binom(3,2)
    dm = d.mean()      # both means used a lot so declaring saves space
    em = e.mean()
    fm = f.mean()

    # the test cases for 0.1 separation
    sep = 0.1
    t1a = np.nan
    t1b = equalSpacing(b, b.mean(), 1)/float(n-1)
    t1c = (0.1-c.mean())*(0.2-c.mean()) #*5**2/5**2
    t1d = (equalSpacing(d[:3], dm, 1) + (.3-dm)*(.4-dm)*5. + 5.*2*(.4-dm)*(.5-dm))/(2.+5.+2.*5.)
    t1e = ((0.4-em)*(0.5-em)*3.*3. + (0.5-em)*(0.6-em)*3*3)/(3.*3.+3.*3)
    t1f = ((0.8-fm)*(0.9-fm)*3. + (0.9-fm)*(1.0-fm)*2.*3.)/(3.+2.*3.)

    # test caess for 0 separation
    sep = 0.0
    t2a = 0.0
    t2b = np.nan
    t2c = ((0.1-c.mean())**2*nCr52 + (0.2-c.mean())**2*nCr52)/(nCr52+nCr52)
    t2d = ((0.4-dm)**2*nCr52 + (.5-dm)**2)/(nCr52+1)
    t2e = ((0.4-em)**2*nCr32 + (.5-em)**2*nCr32 + (0.6-em)**2*nCr32)/(3*nCr32)
    t2f = ((0.9-fm)**2*nCr32 + (1.0-fm)**2)/(nCr32+1)

    cases = [a,b,c,d,e, f]
    test_set1 = [t1a, t1b, t1c, t1d, t1e, t1f]
    test_set2 = [t2a, t2b, t2c, t2d, t2e, t2f]
    return cases, test_set1, test_set2

Useful

• Numbers, pairs, and a difference
• Finding number of number pairs with given difference
• Count all distinct pairs with difference equal to k

Answer 1

My Own Attempt

I am fairly confident that this can pass as an answer now to the problem and is ready for anyone to review and offer speed improvements

Performance

Here are the benchmark tests. I also ran at various other parameters to verify the scaling behaviour as shown in the figure below

diff = 0.1; n = 2000
arr = np.cumsum(np.around(np.random.exponential(diff, n), 1))
mean = arr.mean()
%timeit attempt(arr, diff, mean, n)

n = 20000
arr = np.cumsum(np.around(np.random.exponential(diff, n), 1))
mean = arr.mean()
%timeit attempt(arr, diff, mean, n)

n = 2000000
arr = np.cumsum(np.around(np.random.exponential(diff, n), 1))
mean = arr.mean()
%timeit attempt(arr, diff, mean, n)

## -- End pasted text --
100 loops, best of 3: 2.52 ms per loop
10 loops, best of 3: 25.1 ms per loop
1 loop, best of 3: 2.52 s per loop

Accuracy

Start separation of 0.1

Test outcomes...
  test:1 :: Passed :: res:     nan actual:     nan pairs: 
  test:2 :: Passed :: res:  0.0642 actual:  0.0642 pairs: 1x(0.8 0.9) 1x(0.3 0.4) 1x(0.1 0.2) 1x(0.7 0.8) 1x(0.6 0.7) 1x(0.4 0.5) 1x(0.5 0.6) 1x(0.2 0.3) 1x(0.9 1.0)
  test:3 :: Passed :: res: -0.0025 actual: -0.0025 pairs: 25x(0.1 0.2)
  test:4 :: Passed :: res:  0.0056 actual:  0.0056 pairs: 1x(0.2 0.3) 10x(0.4 0.5) 5x(0.3 0.4) 1x(0.1 0.2)
  test:5 :: Passed :: res:  0.2025 actual:  0.2025 pairs: 9x(0.4 0.5) 9x(0.5 0.6)
  test:6 :: Passed :: res:  0.0606 actual:  0.0606 pairs: 3x(0.8 0.9) 6x(0.9 1.0)

Start separation of 0.0

Test outcomes...
  test:1 :: Passed :: res:  0.0000 actual:  0.0000 pairs: 1x(0.1 0.1)
  test:2 :: Passed :: res:     nan actual:     nan pairs: 
  test:3 :: Passed :: res:  0.0025 actual:  0.0025 pairs: 1x(0.1 0.1) 1x(0.2 0.2)
  test:4 :: Passed :: res:  0.0032 actual:  0.0032 pairs: 1x(0.5 0.5) 1x(0.4 0.4)
  test:5 :: Passed :: res:  0.2092 actual:  0.2092 pairs: 1x(0.5 0.5) 1x(0.6 0.6) 1x(0.4 0.4)
  test:6 :: Passed :: res:  0.0669 actual:  0.0669 pairs: 1x(1.0 1.0) 1x(0.9 0.9)

Code

This code is a shortened version without the overhead debugging. For a longer version which generated the accuracy test see my repository for test.test_acorrMapped

def attemptShort(arr, sep, mean, n, tol=1e-7, **kwargs):
    """Shortened version for Stack Exchange"""
    if sep == 0: # fast exit for 0 separations
        # faster than np.unique as latter requires a mask over counts>1
        unique_counts = np.asarray([(v,c) for v,c in Counter(arr).iteritems() if c>1])
        if not unique_counts.size: return np.nan    # handle no unique items
        combinations = binom(unique_counts[:,1],2)  # get combinations
        return ((unique_counts[:,0]-mean)**2*combinations).sum() / combinations.sum()

    front = 1   # front "pythony-pointer-thing"
    back  = 0   # back "pythony-pointer-thing"
    bssp  = 0   # back sweep start point
    bsfp  = 0   # back sweep finish point
    ans   = 0.0 # store the answer
    count = 0   # counter for averaging
    new_front = True # the first front value is new
    while front < n:            # keep going until exhausted array
        new_front = (arr[front]-arr[front-1]>tol)  # check if front value is a new one
        back = bsfp if new_front else bssp         # this is the magical step

        diff = arr[front] - arr[back]
        if abs(diff - sep) < tol: # if equal subject to tol: pair found
            if new_front:
                bssp  = bsfp    # move sweep start point
                back  = bsfp    # and back to last front point
                bsfp  = front   # send start end point to front's position
            else:
                back  = bssp    # reset back to the sweep start point
            while back < bsfp:  # calculate the correlation function for matched pairs
                count+= 1
                ans  += (arr[front] - mean)*(arr[back] - mean)
                back += 1
        else:
            if abs(arr[bssp+1]- arr[bssp]) > tol: bsfp = front

        front +=1
    return ans/float(count) if count > 0 else np.nan # cannot calculate if no pairs

Applications

Hybrid Monte Carlo - see my repository

My implementation can be found in the function correlations.acorr.acorrMapped.

I assign a map of separations to the MCMC samples that were generated by solving Hamiltons Equations of Motion for exponentially distributed trajectories. Hence, unlike normal autocorrelations, these are separated by an exponentially fictitious length defined by the Equations of Motion.