Skip to main content
deleted 41 characters in body
Source Link

Hybrid Monte Carlo (Deep Learning / Quantum Chromodynamics) - see my repository

Hybrid Monte Carlo (Deep Learning / Quantum Chromodynamics) - see my repository

Hybrid Monte Carlo - see my repository

added 221 characters in body
Source Link

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

enter image description here

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

enter image description here

Source Link

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.