Separating two interleaved arrays into continuous sequences

Question

I'm not sure how to best explain my question in words, so I will provide a code example below. But to at least give it a try. I am solving an eigenvalue problem as a function of some external parameter, which results in two eigenvalues. Those two eigenvalues cross, as the function of the external parameter. Eigenvalue sorting then leads to the wrong classification of which eigenvalue belongs to what 'branch' of the problem. I'd like to disentangle that.

Okay, as for the example. We start from the desired result and then mess it up according to what happens in the diagonalization routine.

xs = np.arange(-np.pi, np.pi, 0.001)
fun_a = np.cos(xs) + 0.5*np.sin(xs)
fun_b = np.cos(xs) - 0.5*xs*np.sin(xs) + 0.2

plt.figure()
plt.plot(xs, fun_a)
plt.plot(xs, fun_b)

This results in two smooth branches that cross each other:

Now, what happens is that instead the eigenvalues are sorted, which I can mimmic as follows:

fun_c = np.zeros_like(fun_a)
fun_c[fun_a>fun_b] = fun_a[fun_a>fun_b]
fun_c[fun_a<=fun_b] = fun_b[fun_a<=fun_b]

fun_d = np.zeros_like(fun_a)
fun_d[fun_a>fun_b] = fun_b[fun_a>fun_b]
fun_d[fun_a<=fun_b] = fun_a[fun_a<=fun_b]

plt.figure()
plt.plot(xs, fun_c)
plt.plot(xs, fun_d)

The result is these two branches, which are no longer smooth, continuous branches

What my question boils down to is how we can go from the second case (fun_c, fun_d) to the first case (fun_a, fun_b). I suppose one could use the information in the derivative, to some extent. Using np.diff() reveals a sharp discontinuity at the point where things go wrong. But I don't immediately see how to nicely use that.

edit: I'm starting to think the best assignment would be such that each subsequent point is chosen to minimize the change in the slope. But I'm not yet sure how to do that..

Answer 1

I'm sure there are many adjustments to be made to improve performance, but here's my guess, based on second derivative thresholding:

# first derivative
deriv = np.diff(fun_c)
# second derivative
deriv2 = np.diff(deriv)
# adjust threshold to detect discontinuities
switch_points = deriv2 > 0.0002
# indices of points of intersection
touch_index = np.where(switch_points == True)[0]
# adjustment to prevent false positives
# (sometimes contigous samples like 127,128 are both detected)
duplicate_index = np.where(np.diff(touch_index) == 1)[0]
touch_index = np.delete(touch_index, duplicate_index)
# begin and ending points of sections to swap
begins = touch_index[::2]
ends = touch_index[1::2]

from itertools import zip_longest

# swap values in selected sections
tmp = fun_c.copy() # need a tmp array to swap
for begin, end in zip_longest(begins, ends, fillvalue=len(fun_c)):
    tmp[begin:end] = fun_d[begin:end]

# complete swapping, correcting fun_d
# on the indices we changed before
swapped = tmp != fun_c
fun_d[swapped] = fun_c[swapped]
fun_c = tmp

plt.plot(xs, fun_c)
plt.plot(xs, fun_d)

I found it quite dependant on the sampling_rate (might fail if it is too low).

Edit: If it wasn't for the

touch_index = np.delete(touch_index, duplicate_index)

you could entirely skip all this:

begins = touch_index[::2]
ends = touch_index[1::2]
from itertools import zip_longest
# swap values in selected sections
tmp = fun_c.copy() # need a tmp array to swap
for begin, end in zip_longest(begins, ends, fillvalue=len(fun_c)):
    tmp[begin:end] = fun_d[begin:end]

with just

np.logical_xor.accumulate(touch_index) | touch_index