# Shared model for daily EAE-style clinical scores (0 to 5 in steps of 0.5), days 7 to 28.
# Each mouse: develops disease with probability P; onset day ~ 10 + Poisson(2) (+ DELAY with treatment);
# the score rises over 3 days to a peak (normal around PEAK, SD 0.8, limited to 1-4.5), then partly recovers by
# 0.05 per day; daily noise SD 0.25; rounded to the nearest 0.5.
import numpy as np
DAYS = np.arange(7, 29)
def simulate(rng, n, p, peak, delay=0):
    s = np.zeros((n, len(DAYS)))
    for i in range(n):
        if rng.random() >= p: continue
        onset = 10 + rng.poisson(2) + delay; pk = np.clip(rng.normal(peak, 0.8), 1, 4.5)
        for j, day in enumerate(DAYS):
            if day < onset: continue
            k = day - onset
            v = pk * min(1, (k + 1) / 3) - 0.05 * max(0, k - 3) + rng.normal(0, 0.25)
            s[i, j] = np.clip(np.round(v * 2) / 2, 0.5, 5)
    return s