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
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 (Deep Learning / Quantum Chromodynamics) - 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.