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

This kind of works, and covers for both cases, when the maxi and mini functions meet, and when there is a "jump". Still, there is a constant which needs adjusting (in my case 0.8) to cover for the slight variations in mini or maxi.

one = one + 2*np.pi*(one < 0)
two = two + 2*np.pi*(two < 0)

plt.figure(figsize=(8, 8))
plt.title(run)

t = np.arange(one.size)

maxi = np.maximum(one, two)
mini = np.minimum(one, two)

breaks = np.array([0])
ct = 0

first = np.array([])
second = np.array([])

ch = [maxi, mini]

for i in range(1, maxi.size):
    if maxi[i] < mini[i-1] or mini[i] > maxi[i-1]:
        breaks = np.append(breaks, [i])
        ct += 1
        first = np.append(first, ch[0][breaks[ct-1]:breaks[ct]])
        second = np.append(second, ch[1][breaks[ct-1]:breaks[ct]])
        ch = ch[::-1]
    elif i != maxi.size-1:
        if maxi[i] - mini[i] < 0.8*min(maxi[i-1] - mini[i-1], maxi[i+1] - mini[i+1]):
            breaks = np.append(breaks, [i])
            ct += 1
            first = np.append(first, ch[0][breaks[ct - 1]:breaks[ct]])
            second = np.append(second, ch[1][breaks[ct - 1]:breaks[ct]])
            ch = ch[::-1]

first = np.append(first, ch[0][breaks[ct]:])
second = np.append(second, ch[1][breaks[ct]:])

plt.plot(maxi, 'r.', label='maxi')
plt.plot(mini, 'b.', label='mini')
plt.plot(first, label='first')
plt.plot(second, label='second')
plt.legend(prop={'size': 8})
plt.show()

So it ends up looking like this:

or this:

These were the true value datasets, it is only going to start being fun when I start using the noisy outputs from an estimation algorithm.

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 :

Answer 3

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]