# Data-generating model shared by make_data.py and simulation.py. Python 3.13, NumPy.
# Normalised log2 expression for 2,000 genes in 4 control and 4 treated samples (after a bulk RNA-seq normalisation).
import numpy as np
GENES, N = 2000, 4
def true_lfc():
    # 300 changed genes with |log2 fold change| evenly spread from 0.25 to 2 (alternating up and down); 1,700 unchanged
    t = np.zeros(GENES); t[:300] = np.linspace(0.25, 2.0, 300) * np.tile([1, -1], 150); return t
def study(rng, lfc=None):
    lfc = true_lfc() if lfc is None else lfc
    mean = rng.normal(7, 2, GENES); sd = np.sqrt(rng.gamma(4, 0.25 ** 2 / 4, GENES))   # gene-wise SD around 0.25 (log2)
    ctl = mean + rng.normal(0, 1, (N, GENES)) * sd
    trt = mean + lfc + rng.normal(0, 1, (N, GENES)) * sd
    return ctl, trt