# Constructed data for the colocalisation guide.
# (1) Three small two-channel images (64 x 64 pixels) showing how Pearson's and Manders' coefficients differ:
#     A: channel 2 sits on the same spots as channel 1 with proportional brightness;
#     B: same spots, but brightness within the spots unrelated between channels;
#     C: half of channel 2's spots lie away from channel 1's.
# (2) Per-cell Pearson coefficients from 3 independent experiments x 2 conditions x 25 cells. There is NO true
#     treatment effect; experiments differ from each other (SD 0.08 on the Fisher z scale) and so do the two
#     conditions within an experiment by chance (SD 0.08), with cell-to-cell SD 0.15.
#     The seed (2) is the first from 1 upward for which a per-cell t-test gives p < 0.05, the situation this page
#     is about; the simulation shows how often that happens.
import numpy as np, pandas as pd
rng = np.random.default_rng(2)

def spots(rng, n=25, size=64):
    yy, xx = np.mgrid[0:size, 0:size]; img = np.zeros((size, size)); centres = rng.integers(4, size - 4, (n, 2))
    for cy, cx in centres:
        img += np.exp(-((yy - cy) ** 2 + (xx - cx) ** 2) / 6.0)
    return img, centres

base, centres = spots(rng)
bg = lambda: rng.normal(10, 2, (64, 64))
yy, xx = np.mgrid[0:64, 0:64]
ch1 = 100 * base + bg()
A2 = 80 * base + bg()
w = np.zeros((64, 64))
for cy, cx in centres:
    w += rng.uniform(0.2, 1.8) * np.exp(-((yy - cy) ** 2 + (xx - cx) ** 2) / 6.0)
B2 = 80 * w + bg()
other, _ = spots(rng, n=13)
half = np.zeros((64, 64))
for cy, cx in centres[:12]:
    half += np.exp(-((yy - cy) ** 2 + (xx - cx) ** 2) / 6.0)
C2 = 80 * (half + other) + bg()
rows = []
for name, c2 in (("A", A2), ("B", B2), ("C", C2)):
    for (y, x), v1 in np.ndenumerate(ch1):
        rows.append((name, x, y, round(v1, 2), round(c2[y, x], 2)))
pd.DataFrame(rows, columns=["image", "x", "y", "ch1", "ch2"]).to_csv("images.csv", index=False)

rows = []
for e in range(1, 4):
    ue = rng.normal(0, 0.08)
    for cond in ("control", "treated"):
        uc = rng.normal(0, 0.08)
        z = np.arctanh(0.6) + ue + uc + rng.normal(0, 0.15, 25)
        for i, zi in enumerate(z):
            rows.append((e, cond, f"e{e}{cond[0]}{i+1}", round(float(np.tanh(zi)), 4)))
pd.DataFrame(rows, columns=["experiment", "condition", "cell", "pearson"]).to_csv("cells.csv", index=False)