def simulation(n_obs=<NUM_LIT>, n_vars=None, alpha=None, beta=None, gamma=None, alpha_=None, t_max=None,<EOL>noise_model='<STR_LIT>', noise_level=<NUM_LIT:1>, switches=None):
from .tools.dynamical_model_utils import unspliced, spliced, vectorize<EOL>def draw_poisson(n):<EOL><INDENT>from random import uniform  <EOL>t = np.cumsum([- <NUM_LIT> * np.log(uniform(<NUM_LIT:0>, <NUM_LIT:1>)) for _ in range(n - <NUM_LIT:1>)])<EOL>return np.insert(t, <NUM_LIT:0>, <NUM_LIT:0>)  <EOL><DEDENT>def simulate_dynamics(tau, alpha, beta, gamma, u0, s0, noise_model, noise_level):<EOL><INDENT>ut = unspliced(tau, u0, alpha, beta)<EOL>st = spliced(tau, s0, u0, alpha, beta, gamma)<EOL>if noise_model is '<STR_LIT>':  <EOL><INDENT>ut += np.random.normal(scale=noise_level * np.percentile(ut, <NUM_LIT>) / <NUM_LIT:10>, size=len(ut))<EOL>st += np.random.normal(scale=noise_level * np.percentile(st, <NUM_LIT>) / <NUM_LIT:10>, size=len(st))<EOL><DEDENT>ut, st = np.clip(ut, <NUM_LIT:0>, None), np.clip(st, <NUM_LIT:0>, None)<EOL>return ut, st<EOL><DEDENT>def simulate_gillespie(alpha, beta, gamma):<EOL><INDENT>update_rule = np.array([[<NUM_LIT:1>, <NUM_LIT:0>], [-<NUM_LIT:1>, <NUM_LIT:1>], [<NUM_LIT:0>, -<NUM_LIT:1>], [<NUM_LIT:0>, <NUM_LIT:0>]])<EOL>def update(props):<EOL><INDENT>if np.sum(props) > <NUM_LIT:0>:<EOL><INDENT>props /= np.sum(props)<EOL><DEDENT>p_cumsum = props.cumsum()<EOL>p = np.random.rand()<EOL>i = <NUM_LIT:0><EOL>while p > p_cumsum[i]:<EOL><INDENT>i += <NUM_LIT:1><EOL><DEDENT>return update_rule[i]<EOL><DEDENT>u, s = np.zeros(len(alpha)), np.zeros(len(alpha))<EOL>for i, alpha_i in enumerate(alpha):<EOL><INDENT>u_, s_ = (u[i - <NUM_LIT:1>], s[i - <NUM_LIT:1>]) if i > <NUM_LIT:0> else (<NUM_LIT:0>, <NUM_LIT:0>)<EOL>du, ds = update(props=np.array([alpha_i, beta * u_, gamma * s_]))<EOL>u[i], s[i] = (u_ + du, s_ + ds)<EOL><DEDENT>return u, s<EOL><DEDENT>alpha = <NUM_LIT:5> if alpha is None else alpha<EOL>beta = <NUM_LIT> if beta is None else beta<EOL>gamma = <NUM_LIT> if gamma is None else gamma<EOL>alpha_ = <NUM_LIT:0> if alpha_ is None else alpha_<EOL>t = draw_poisson(n_obs)<EOL>if t_max is not None:<EOL><INDENT>t *= t_max / np.max(t)<EOL><DEDENT>t_max = np.max(t)<EOL>def cycle(array, n_vars=None):<EOL><INDENT>if isinstance(array, (list, tuple)):<EOL><INDENT>return array if n_vars is None else array * int(np.ceil(n_vars / len(array)))<EOL><DEDENT>else:<EOL><INDENT>return [array] if n_vars is None else [array] * n_vars<EOL><DEDENT><DEDENT>switches = cycle([<NUM_LIT>, <NUM_LIT>, <NUM_LIT:1>, <NUM_LIT>], n_vars) if switches is None else cycle(switches, n_vars)<EOL>t_ = np.array([np.max(t[t < t_i * t_max]) for t_i in switches])<EOL>noise_level = cycle(noise_level, len(switches) if n_vars is None else n_vars)<EOL>n_vars = max(len(switches), len(noise_level))<EOL>U = np.zeros(shape=(len(t), n_vars))<EOL>S = np.zeros(shape=(len(t), n_vars))<EOL>for i in range(n_vars):<EOL><INDENT>alpha_i = alpha[i] if isinstance(alpha, (tuple, list, np.ndarray)) else alpha<EOL>beta_i = beta[i] if isinstance(beta, (tuple, list, np.ndarray)) else beta<EOL>gamma_i = gamma[i] if isinstance(gamma, (tuple, list, np.ndarray)) else gamma<EOL>tau, alpha_vec, u0_vec, s0_vec = vectorize(t, t_[i], alpha_i, beta_i, gamma_i, alpha_=alpha_, u0=<NUM_LIT:0>, s0=<NUM_LIT:0>)<EOL>beta_i<EOL>if noise_model is '<STR_LIT>':<EOL><INDENT>U[:, i], S[:, i] = simulate_gillespie(alpha_vec, beta, gamma)<EOL><DEDENT>else:<EOL><INDENT>U[:, i], S[:, i] = simulate_dynamics(tau, alpha_vec, beta_i, gamma_i,<EOL>u0_vec, s0_vec, noise_model, noise_level[i])<EOL><DEDENT><DEDENT>obs = {'<STR_LIT>': t.round(<NUM_LIT:2>)}<EOL>var = {'<STR_LIT>': t_[:n_vars],<EOL>'<STR_LIT>': np.ones(n_vars) * alpha,<EOL>'<STR_LIT>': np.ones(n_vars) * beta,<EOL>'<STR_LIT>': np.ones(n_vars) * gamma,<EOL>'<STR_LIT>': np.ones(n_vars)}<EOL>layers = {'<STR_LIT>': U, '<STR_LIT>': S}<EOL>return AnnData(S, obs, var, layers=layers)<EOL>
Simulation of mRNA metabolism with transcription, splicing and degradation

    Returns
    -------
    Returns `adata` object