<jupyter_start><jupyter_text>Get the first bifurcation point (without randomisations)Once you have this, you can draw a line segment bisecting the angle between the point and two targets. This will be the line about which you symmetrise to get the second bifurcation point<jupyter_code>x_fit = []
y_fit = []

bif_pt = 0.8
params1 = [1, 1, 0.5]

xs1 = xs[xs < 2.9]
ys1 = ys[xs < 2.9]
ts1 = ts[xs < 2.9]

x_sub = np.concatenate((xs1, xs1))
y_sub = np.concatenate((ys1, -ys1))
t_sub = np.concatenate((ts1, ts1))

tmin = 0
tmax = 25
for idx, t in enumerate(range(tmin, tmax)):
    window_min = t
    window_max = t + window_size

    x = x_sub[(t_sub > window_min) & (t_sub < window_max)]
    y = y_sub[(t_sub > window_min) & (t_sub < window_max)]
    tmp_img2 = density_map(x, y, stats=True)

    if idx == 0:
        tmp_img = tmp_img2
    else:
        tmp_img = np.fmax(tmp_img2, tmp_img)

x_fit = np.where(tmp_img > peak_threshold)[1] * post1_x / nbins
y_fit = (
    np.where(tmp_img > peak_threshold)[0] * (post0_y - post2_y) / nbins + post2_y
)

x_fit = x_fit
y_fit = np.abs(y_fit)
y_fit = y_fit[x_fit > bif_pt]
x_fit = x_fit[x_fit > bif_pt]

for i in range(0, 10):
    fit_params, pcov = curve_fit(
        fitfunc_vec_self, x_fit, y_fit, p0=params1, maxfev=10000
    )
    params1 = fit_params
bisector_xs = [params1[0], post2_x]
bisector_ys = [
    0,
    np.tan(np.arctan2(post2_y, post2_x - params1[0]) / 2) * (post2_x - params1[0]),
]<jupyter_output><empty_output><jupyter_text>Get the second bifurcation pointFor this, you must center the trajectories about the bifurcation point, get a new heatmap and rotate this by the angle of the bisector line<jupyter_code># center points about the first bifurcation
cxs = xs - params1[0]
cys = ys
cts = ts

cpost0_x = post0_x - params1[0]
cpost1_x = post1_x - params1[0]
cpost2_x = post2_x - params1[0]
@numba.njit(fastmath=True, parallel=True)
def parallel_rotate(xy, rmat):
    out = np.zeros(xy.shape)
    for idx in numba.prange(xy.shape[0]):
        out[idx] = np.dot(rmat[idx], xy[idx])
    return out
# clip all points to the left of and below 0 and points beyond post centers
ccxs = cxs[cxs > 0]
ccys = cys[cxs > 0]
ccts = cts[cxs > 0]
ccxs = ccxs[ccys > 0]
ccts = ccts[ccys > 0]
ccys = ccys[ccys > 0]

xy = np.concatenate((ccxs.reshape(-1, 1), ccys.reshape(-1, 1)), axis=1)
angle = np.full(ccxs.shape, np.arctan2(post2_y, post2_x - params1[0]) / 2)
rmat = np.array([[np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)]]).T
rx, ry = parallel_rotate(xy, rmat).T

blur = (101, 101)

r1 = [[0, post1_x], [post0_y, post2_y]]

tmin = 0
tmax = 25
for idx, t in enumerate(range(tmin, tmax)):
    window_min = t
    window_max = t + window_size

    x = rx[(ccts > window_min) & (ccts < window_max)]
    y = ry[(ccts > window_min) & (ccts < window_max)]
    tmp_img = density_map(x, y, stats=False)

    if idx == 0:
        tmp_img1 = tmp_img
    else:
        tmp_img1 = np.fmax(tmp_img1, tmp_img)
plt.imshow(tmp_img1, extent=[r1[0][0], r1[0][1], r1[1][0], r1[1][1]])
positive_bif2 = []

while True:
    bif_pt = 2
    params2 = [2.5, 1, 0.5]

    x_sub = np.concatenate((rx, rx))
    y_sub = np.concatenate((ry, -ry))
    np.random.shuffle(y_sub)
    t_sub = np.concatenate((ccts, ccts))

    tmin = 0
    tmax = 25
    for idx, t in enumerate(range(tmin, tmax, 2)):
        window_min = t
        window_max = t + window_size

        x = x_sub[(t_sub > window_min) & (t_sub < window_max)]
        y = y_sub[(t_sub > window_min) & (t_sub < window_max)]
        tmp_img = density_map(x, y, stats=False)

        if idx == 0:
            tmp_img1 = tmp_img
        else:
            tmp_img1 = np.fmax(tmp_img1, tmp_img)

    x_fit = np.where(tmp_img1 > peak_threshold)[1] * post1_x / nbins
    y_fit = np.where(tmp_img1 > peak_threshold)[0] * (post0_y - post2_y) / nbins + post2_y

    x_fit = x_fit
    y_fit = np.abs(y_fit)
    y_fit = y_fit[x_fit > bif_pt]
    x_fit = x_fit[x_fit > bif_pt]

    try:
        for i in range(0, 10):
            fit_params, pcov = curve_fit(
                fitfunc_vec_self, x_fit, y_fit, p0=params2, maxfev=10000
            )
            params2 = fit_params

        if (
            (fit_params[0] > 0)
            & (fit_params[0] < post0_x)
            & (fit_params[1] > 0.1)
            & (fit_params[2] > 0.1)
        ):
            positive_bif2.append(1)
        else:
            positive_bif2.append(0)
    except TypeError:
        pass

    if len(positive_bif2) >= randomisations:
        break

    if len(positive_bif2) % 10 == 0:
        print(len(positive_bif2))
np.sum(positive_bif2) / len(positive_bif2), len(positive_bif2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13540106951871658, 'fraction_numerical': 0.038716577540106954, 'mean_word_length': 2.936026936026936, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '291250', 'license': 'MIT', 'n_tokens_mistral': 2055, 'n_tokens_neox': 1934, 'n_words': 446, 'path': 'bifurcation_statistics.ipynb', 'repo_name': 'RnoB/fly-matrix', 'size': '4675'}