<div dir="ltr">Here's some code implementing the "replace similar values with an arbitrarily chosen one" (in this case the smallest of the similar values). I didn't see any way to do this cleverly with strides, so I just did a simple loop. It's about 100 times slower in pure Python, or a bit under 10 times slower if you're willing to use a bit of Cython. Not sure if this is good enough for your purposes. I imagine you could go a bit faster if you were willing to do the lexical integration by hand (since you've already done the separate sorting of each subarray for value replacement purposes) instead of passing that off to np.lexsort.<br>
<div><div><br></div><div>Note that this approach will only work if your points are not only well-separated in space but also either well-separated or identical in each dimension as well. It's OK to have points with the same, say, x value, but if you have points that have close x values before the noise is added, then the noise can move intermediate points around in the sort order. It works well with the gridded data I used as a sample, but if you're, say, generating random points, this could be a problem:</div>
<div><br></div><div>point 1 is (1, 0, 1e-12)</div><div>point 2 is (0, 1, 0)</div><div><br></div><div>These are well separated. The algorithm will pool those z values and report 1 as coming before 2. Unless you get jitter like this:</div>
<div><br></div><div>point 1: (1, 0, 1.5e-12)</div><div>point 2: (0, 1, -0.5e-12)</div><div><br></div><div>Now they won't be pooled any more and we'll get 2 as coming before 1.</div><div><br></div><div>Anyway, here's the code:</div>
<div><br></div><div>In [1]:</div><div><br></div><div>import numpy as np</div><div>In [2]:</div><div><br></div><div>def gen_grid(n, d):</div><div> #Generate a bunch of grid points, n in each dimension of spacing d</div>
<div> vals = np.linspace(0, (n-1)*d, n)</div><div> x, y, z = np.meshgrid(vals, vals, vals)</div><div> grid = np.empty((n**3, 3))</div><div> grid[:,0] = x.flatten()</div><div> grid[:,1] = y.flatten()</div><div>
grid[:,2] = z.flatten()</div><div> return grid</div><div> </div><div>def jitter(array, epsilon=1e-12):</div><div> #Add random jitter from a uniform distribution of width epsilon</div><div> return array + np.random.random(array.shape) * epsilon - epsilon / 2</div>
<div>In [3]:</div><div><br></div><div>grid = gen_grid(4, 0.1)</div><div>print np.lexsort(grid.T)</div><div>print np.lexsort(jitter(grid.T))</div><div>[ 0 4 8 12 16 20 24 28 32 36 40 44 48 52 56 60 1 5 9 13 17 21 25 29 33</div>
<div> 37 41 45 49 53 57 61 2 6 10 14 18 22 26 30 34 38 42 46 50 54 58 62 3 7</div><div> 11 15 19 23 27 31 35 39 43 47 51 55 59 63]</div><div>[60 4 48 32 40 12 36 28 44 56 16 8 24 0 52 20 45 25 49 1 53 29 9 33 5</div>
<div> 61 41 37 17 13 21 57 22 50 18 10 2 62 58 54 6 34 26 42 38 46 14 30 3 11</div><div> 55 63 27 15 35 43 31 39 7 59 47 23 51 19]</div><div>In [4]:</div><div><br></div><div>def pool_values(A, epsilon=1e-12):</div><div>
idx = np.argsort(A)</div><div> for i in range(1, len(A)):</div><div> if A[idx[i]] - A[idx[i-1]] < epsilon:</div><div> A[idx[i]] = A[idx[i-1]]</div><div> return A</div><div> </div><div>def stable_sort(grid):</div>
<div> return np.lexsort((pool_values(grid[:,0]),</div><div> pool_values(grid[:,1]),</div><div> pool_values(grid[:,2])))</div><div>In [5]:</div><div><br></div><div>print stable_sort(grid)</div>
<div>print stable_sort(jitter(grid))</div><div>[ 0 4 8 12 16 20 24 28 32 36 40 44 48 52 56 60 1 5 9 13 17 21 25 29 33</div><div> 37 41 45 49 53 57 61 2 6 10 14 18 22 26 30 34 38 42 46 50 54 58 62 3 7</div><div> 11 15 19 23 27 31 35 39 43 47 51 55 59 63]</div>
<div>[ 0 4 8 12 16 20 24 28 32 36 40 44 48 52 56 60 1 5 9 13 17 21 25 29 33</div><div> 37 41 45 49 53 57 61 2 6 10 14 18 22 26 30 34 38 42 46 50 54 58 62 3 7</div><div> 11 15 19 23 27 31 35 39 43 47 51 55 59 63]</div>
<div>In [6]:</div><div><br></div><div>%timeit np.lexsort(jitter(grid.T))</div><div>100000 loops, best of 3: 10.4 µs per loop</div><div>In [7]:</div><div><br></div><div>%timeit stable_sort(jitter(grid))</div><div>1000 loops, best of 3: 1.39 ms per loop</div>
<div>In [8]:</div><div><br></div><div>%load_ext cythonmagic</div><div>In [12]:</div><div><br></div><div>%%cython</div><div>import numpy as np</div><div>cimport numpy as np</div><div> </div><div>cdef fast_pool_values(double[:] A, double epsilon=1e-12):</div>
<div> cdef long[:] idx = np.argsort(A)</div><div> cdef int i</div><div> for i in range(1, len(A)):</div><div> if A[idx[i]] - A[idx[i-1]] < epsilon:</div><div> A[idx[i]] = A[idx[i-1]]</div><div>
return A</div><div> </div><div>def fast_stable_sort(grid):</div><div> return np.lexsort((fast_pool_values(grid[:,0]),</div><div> fast_pool_values(grid[:,1]),</div><div> fast_pool_values(grid[:,2])))</div>
<div>In [10]:</div><div><br></div><div>%timeit np.lexsort(jitter(grid.T))</div><div>10000 loops, best of 3: 38.5 µs per loop</div><div>In [13]:</div><div><br></div><div>%timeit fast_stable_sort(jitter(grid))</div><div>1000 loops, best of 3: 309 µs per loop</div>
</div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Sun, Oct 27, 2013 at 5:41 PM, Freddie Witherden <span dir="ltr"><<a href="mailto:freddie@witherden.org" target="_blank">freddie@witherden.org</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="im">On 27/10/13 21:05, Jonathan March wrote:<br>
> If an "almost always works" solution is good enough, then sort on the<br>
> distance to some fixed random point that is in the vicinity of your N<br>
> points.<br>
<br>
</div>I had considered this. Unfortunately I need a solution which really<br>
does always work.<br>
<br>
The only pure-Python solution I can envision -- at the moment anyway --<br>
is to do some cleverness with the output of np.unique to identify<br>
similar values and replace them with an arbitrarily chosen one. This<br>
should permit the output to be passed to np.lexsort without issue.<br>
<br>
Regards, Freddie.<br>
<br>
<br>
<br>
<br>
<br>_______________________________________________<br>
NumPy-Discussion mailing list<br>
<a href="mailto:NumPy-Discussion@scipy.org">NumPy-Discussion@scipy.org</a><br>
<a href="http://mail.scipy.org/mailman/listinfo/numpy-discussion" target="_blank">http://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br>
<br></blockquote></div><br></div>