Distinguishing between two "symmetrical" arrays using numpy

Question

I have two arrays of values of angles in radians. The two arrays are symmetrical with respect to a known constant angle. The arrays are shown in the figure:

The example values are following:

   one = [ 2.98153965 -1.33298928  2.94993567 -1.39909924  2.99214403  3.00138863
  3.04390642 -1.59098448 -1.65660299 -1.73146174 -1.8166248  -2.85595599
 -2.02035274 -2.64530394 -2.26451127 -2.3982946  -2.52735954 -2.17570346
 -2.77544658 -2.88566686 -1.84913768 -3.07261908 -1.66738719 -1.6029932
 -1.54596053 -1.50177363 -1.46133745 -1.42288915 -1.38241718  2.79925996
 -1.30775884 -1.27309395  2.72153718 -1.20592812 -1.18113435 -1.15029987]

two = [-1.30507254  2.9385436  -1.36415496  2.95897805 -1.43845065 -1.48295087
 -1.53346541  3.09685482 -3.11358085 -3.0466034  -2.95794156 -1.9128659
 -2.75067133 -2.13826992 -2.51567194 -2.39565127 -2.28148844 -2.65519436
 -2.05312249 -1.95523663 -2.98473857 -1.75415233  3.13322155  3.06539723
  3.00595703  2.95378704  2.90786779  2.86730208  2.831318   -1.34113191
  2.77057495  2.74479777 -1.23620286  2.70046364  2.68129889  2.66380717]

It can be seen that the values are "following" two symmetrical arctan lines, my question is how do I distinguish between the two of them, and get something like this:

I've tried several approaches but can't come up with a universal one which will work in all cases, there is often a misinterpreted section assigned to the wrong array.

Any ideas are welcome! Thanks!

Answer 1

Here is a solution that minimizes the distance between consecutive points and also the change in slope (weighted by parameter lam). Distance alone fails at the crossover point.

import numpy as np

one = list(map(float, """ 2.98153965 -1.33298928  2.94993567 -1.39909924  2.99214403  3.00138863
  3.04390642 -1.59098448 -1.65660299 -1.73146174 -1.8166248  -2.85595599
 -2.02035274 -2.64530394 -2.26451127 -2.3982946  -2.52735954 -2.17570346
 -2.77544658 -2.88566686 -1.84913768 -3.07261908 -1.66738719 -1.6029932
 -1.54596053 -1.50177363 -1.46133745 -1.42288915 -1.38241718  2.79925996
 -1.30775884 -1.27309395  2.72153718 -1.20592812 -1.18113435 -1.15029987""".split()))

two = list(map(float, """-1.30507254  2.9385436  -1.36415496  2.95897805 -1.43845065 -1.48295087
 -1.53346541  3.09685482 -3.11358085 -3.0466034  -2.95794156 -1.9128659
 -2.75067133 -2.13826992 -2.51567194 -2.39565127 -2.28148844 -2.65519436
 -2.05312249 -1.95523663 -2.98473857 -1.75415233  3.13322155  3.06539723
  3.00595703  2.95378704  2.90786779  2.86730208  2.831318   -1.34113191
  2.77057495  2.74479777 -1.23620286  2.70046364  2.68129889  2.66380717""".split()))

data = np.array([one, two])

dd = (data[[[0, 1], [1, 0]], 1:] - data[:, None, :-1] + np.pi)%(2*np.pi) - np.pi
dde2 = np.einsum('ijk,ijk->jk', dd, dd)

xovr1 = np.argmin(dde2, axis=0)
pick1 = np.r_[0, np.cumsum(xovr1) & 1]

d2d = dd[:, :, None, 1:] - dd[[[1, 0], [0, 1]], :, :-1]
d2de2 = np.r_['2', np.zeros((2, 2, 1)), np.einsum('ijkl,ijkl->jkl', d2d, d2d)]

lam = 0.5
e2 = (dde2[:, None, :] + lam * d2de2).reshape(4, -1)

xovr2 = np.argmin(e2, axis=0)>>1
pick2 = np.r_[0, np.cumsum(xovr2) & 1]

print('by position only')
print(data[pick1, np.arange(data.shape[1])])
print(data[1-pick1, np.arange(data.shape[1])])

print('by position and slope')
print(data[pick2, np.arange(data.shape[1])])
print(data[1-pick2, np.arange(data.shape[1])])


# by position only
# [ 2.98153965  2.9385436   2.94993567  2.95897805  2.99214403  3.00138863
#   3.04390642  3.09685482 -3.11358085 -3.0466034  -2.95794156 -2.85595599
#  -2.75067133 -2.64530394 -2.51567194 -2.3982946  -2.52735954 -2.65519436
#  -2.77544658 -2.88566686 -2.98473857 -3.07261908  3.13322155  3.06539723
#   3.00595703  2.95378704  2.90786779  2.86730208  2.831318    2.79925996
#   2.77057495  2.74479777  2.72153718  2.70046364  2.68129889  2.66380717]
# [-1.30507254 -1.33298928 -1.36415496 -1.39909924 -1.43845065 -1.48295087
#  -1.53346541 -1.59098448 -1.65660299 -1.73146174 -1.8166248  -1.9128659
#  -2.02035274 -2.13826992 -2.26451127 -2.39565127 -2.28148844 -2.17570346
#  -2.05312249 -1.95523663 -1.84913768 -1.75415233 -1.66738719 -1.6029932
#  -1.54596053 -1.50177363 -1.46133745 -1.42288915 -1.38241718 -1.34113191
#  -1.30775884 -1.27309395 -1.23620286 -1.20592812 -1.18113435 -1.15029987]
# by position and slope
# [ 2.98153965  2.9385436   2.94993567  2.95897805  2.99214403  3.00138863
#   3.04390642  3.09685482 -3.11358085 -3.0466034  -2.95794156 -2.85595599
#  -2.75067133 -2.64530394 -2.51567194 -2.39565127 -2.28148844 -2.17570346
#  -2.05312249 -1.95523663 -1.84913768 -1.75415233 -1.66738719 -1.6029932
#  -1.54596053 -1.50177363 -1.46133745 -1.42288915 -1.38241718 -1.34113191
#  -1.30775884 -1.27309395 -1.23620286 -1.20592812 -1.18113435 -1.15029987]
# [-1.30507254 -1.33298928 -1.36415496 -1.39909924 -1.43845065 -1.48295087
#  -1.53346541 -1.59098448 -1.65660299 -1.73146174 -1.8166248  -1.9128659
#  -2.02035274 -2.13826992 -2.26451127 -2.3982946  -2.52735954 -2.65519436
#  -2.77544658 -2.88566686 -2.98473857 -3.07261908  3.13322155  3.06539723
#   3.00595703  2.95378704  2.90786779  2.86730208  2.831318    2.79925996
#   2.77057495  2.74479777  2.72153718  2.70046364  2.68129889  2.66380717]

Answer 2

The difficulty come from the 2 pi jump, which can be resolved by :

def transform(x):
    return x+2*pi*(x<0)

This function transform the arrays in continuous ones. You must first turn your lists in ndarrays.

then :

t=arange(one.size)    
tone = transform(one)
ttwo = transform(two)

maxi=np.maximum(tone,ttwo)
subplot(211)
plot(t,tone,'o',t,ttwo,'o',maxi)

induces what to do :

i=maxi.argmin()
dicrease = np.choose(np.logical_xor(tone>ttwo,t<i),[tone,ttwo])
increase = np.choose(np.logical_xor(tone>ttwo,t<i),[ttwo,tone])
subplot(212)
plot(t,dicrease,label='dicrease')
plot(t,increase,label='increase')
legend()    

for

You can if necessary turn back in [-pi,pi[ by x -> (x + pi) % (2*pi) - pi .

EDIT

for a less ad hoc transform, I propose this other, which will probably solve more cases :

def transform2(y,gap):
    breaks=np.diff(y)**2>gap**2/2
    signs=np.sign(np.diff(y))
    offset=np.concatenate(([0],(breaks*signs).cumsum()))*gap
    return y-offset 

and an noisy example :