Helper function.
199308
def sample_dynamic_data(params):
    
    model_key = os.path.basename(params['model']).replace('.txt', '')
    if 'writedir' not in params or params['writedir'] is None:
        params['writedir'] = settings.writedir + model_key + '_sim'
    if not os.path.exists(params['writedir']): os.makedirs(params['writedir'])
    readwrite.write_params(params['writedir'] + '/params.txt', params)
    # init variables
    writedir = params['writedir']
    tmax = params['tmax']
    branching = params['branching']
    noiseObs = params['noiseObs']
    noiseDyn = params['noiseDyn']
    nrRealizations = params['nrRealizations']
    step = params['step']  # step size for saving the figure

    nrSamples = 1  # how many files?
    maxRestarts = 1000
    maxNrSamples = 1

    # simple vector auto regressive process or
    # hill kinetics process simulation
    if 'krumsiek11' not in model_key:
        # create instance, set seed
        grnsim = GRNsim(model=model_key, params=params)
        nrOffEdges_list = np.zeros(nrSamples)
        for sample in range(nrSamples):
            # random topology / for a given edge density
            if 'hill' not in model_key:
                Coupl = np.array(grnsim.Coupl)
                for sampleCoupl in range(10):
                    nrOffEdges = 0
                    for gp in range(grnsim.dim):
                        for g in range(grnsim.dim):
                            # only consider off-diagonal edges
                            if g != gp:
                                Coupl[gp, g] = 0.7 if np.random.rand() < 0.4 else 0
                                nrOffEdges += 1 if Coupl[gp, g] > 0 else 0
                            else:
                                Coupl[gp, g] = 0.7
                    # check that the coupling matrix does not have eigenvalues
                    # greater than 1, which would lead to an exploding var process
                    if max(sp.linalg.eig(Coupl)[0]) < 1:
                        break
                nrOffEdges_list[sample] = nrOffEdges
                grnsim.set_coupl(Coupl)
            # init type
            real = 0
            X0 = np.random.rand(grnsim.dim)
            Xsamples = []
            for restart in range(nrRealizations+maxRestarts):
                # slightly break symmetry in initial conditions
                if 'toggleswitch' in model_key:
                    X0 = (np.array([0.8 for i in range(grnsim.dim)])
                          + 0.01*np.random.randn(grnsim.dim))
                X = grnsim.sim_model(tmax=tmax, X0=X0,
                                     noiseDyn=noiseDyn)
                # check branching
                check = True
                if branching:
                    check, Xsamples = _check_branching(X, Xsamples, restart)
                if check:
                    real += 1
                    grnsim.write_data(X[::step], dir=writedir,
                                      noiseObs=noiseObs,
                                      append=(False if restart==0 else True),
                                      branching=branching,
                                      nrRealizations=nrRealizations)
                # append some zeros
                if 'zeros' in writedir and real == 2:
                    grnsim.write_data(noiseDyn*np.random.randn(500,3), dir=writedir,
                                      noiseObs=noiseObs,
                                      append=(False if restart==0 else True),
                                      branching=branching,
                                      nrRealizations=nrRealizations)
                if real >= nrRealizations:
                    break
        if False:
            logg.info('mean nr of offdiagonal edges',nrOffEdges_list.mean(),
                      'compared to total nr',grnsim.dim*(grnsim.dim-1)/2.)

    # more complex models
    else:
        initType = 'random'

        dim = 11
        step = 5

        grnsim = GRNsim(dim=dim,initType=initType,model=model_key,params=params)
        curr_nrSamples = 0
        Xsamples = []
        for sample in range(maxNrSamples):
            # choose initial conditions such that branchings result
            if initType == 'branch':
                X0mean = grnsim.branch_init_model1(tmax)
                if X0mean is None:
                    grnsim.set_coupl()
                    continue
            real = 0
            for restart in range(nrRealizations+maxRestarts):
                if initType == 'branch':
                    # vary initial conditions around mean
                    X0 = X0mean + (0.05*np.random.rand(dim) - 0.025*np.ones(dim))
                else:
                    # generate random initial conditions within [0.3,0.7]
                    X0 = 0.4*np.random.rand(dim)+0.3
                if model_key in [5,6]:
                    X0 = np.array([0.3,0.3,0,0,0,0])
                if model_key in [7,8,9,10]:
                    X0 = 0.6*np.random.rand(dim)+0.2
                    X0[2:] = np.zeros(4)
                if 'krumsiek11' in model_key:
                    X0 = np.zeros(dim)
                    X0[grnsim.varNames['Gata2']] = 0.8
                    X0[grnsim.varNames['Pu.1']] = 0.8
                    X0[grnsim.varNames['Cebpa']] = 0.8
                    X0 += 0.001*np.random.randn(dim)
                    if False:
                        switch_gene = restart - (nrRealizations - dim)
                        if switch_gene >= dim:
                            break
                        X0[switch_gene] = 0 if X0[switch_gene] > 0.1 else 0.8
                X = grnsim.sim_model(tmax,X0=X0,
                                     noiseDyn=noiseDyn,
                                     restart=restart)
                # check branching
                check = True
                if branching:
                    check, Xsamples = _check_branching(X,Xsamples,restart)
                if check:
                    real += 1
                    grnsim.write_data(X[::step],dir=writedir,noiseObs=noiseObs,
                                     append=(False if restart==0 else True),
                                     branching=branching,
                                     nrRealizations=nrRealizations)
                if real >= nrRealizations:
                    break
    import glob
    # load the last simulation file
    filename = glob.glob(writedir + '/sim*.txt')[-1]
    logg.info('reading simulation results', filename)
    adata = readwrite._read(filename, first_column_names=True,
                            suppress_cache_warning=True)
    adata.uns['tmax_write'] = tmax/step
    return adata
747613
def ILIKE(pattern):
    """"""
    return P(lambda x: fnmatch.fnmatch(x.lower(), pattern.lower()))
78593
def Areml_eigh(self):
        """"""
        s,U = LA.eigh(self.Areml(),lower=True)
        i_pos = (s>1e-10)
        s = s[i_pos]
        U = U[:,i_pos]
        return s,U
1030201
def getvalue(x):
    """"""
    if isrepeating(x):
        raise TypeError(
            "Ambiguous call to getvalue for %r which has more than one value."
            % x)

    for value in getvalues(x):
        return value
179828
def _nargs(f) -> Optional[int]:
    '''
        
    '''
    if isinstance(f, Function):
        return f.nargs
    spec = inspect.getfullargspec(f)
    if spec.varargs is not None:
        return None
    return len(spec.args)
1374683
def _call(self, x):
        """"""

        sum_constr = abs(x.ufuncs.sum() / self.diameter - 1) <= self.sum_rtol

        nonneq_constr = x.ufuncs.greater_equal(0).asarray().all()

        if sum_constr and nonneq_constr:
            return 0
        else:
            return np.inf
886300
def _styleof(expr, styles):
    """"""
    style = dict()
    for expr_filter, sty in styles:
        if expr_filter(expr):
            style.update(sty)
    return style
702094
def model_function(self, model_name, version, func_name):
        """"""
        assert func_name in ('serializer', 'loader', 'invalidator')
        name = "%s_%s_%s" % (model_name.lower(), version, func_name)
        return getattr(self, name)
673363
def _rgetattr(obj, key):
    """"""
    for k in key.split("."):
        obj = getattr(obj, k)
    return obj
732795
def pick(seq, func, maxobj=None):
    """"""
    maxscore = None
    for obj in seq:
        score = func(obj)
        if maxscore is None or maxscore < score:
            (maxscore, maxobj) = (score, obj)
    return maxobj
975286
def _call(self, x):
        """"""
        if self.operator is None:
            return self.vector.inner(x) + self.constant
        elif self.vector is None:
            return x.inner(self.operator(x)) + self.constant
        else:
            tmp = self.operator(x)
            tmp += self.vector
            return x.inner(tmp) + self.constant
968974
def _K(self, R):
        """"""
        return self._ns * self._N / R / self._sin_alpha
1185109
def scalar(cls, value):
    """"""
    return lambda i, name: cls(i, name, (value,), lambda a: a[0], None)
460466
def _call(self, x, out=None):
        """"""
        if out is None:
            out = self.range.zero()
        else:
            out.set_zero()

        out[self.index] = x
        return out
681123
def _call(self, x, out=None):
        """"""
        if out is None:
            out = self.scalar * x
        else:
            out.lincomb(self.scalar, x)
        return out
1157550
def mu(self, value):
        """"""
        _checkParam('mu', value, self.PARAMLIMITS, self.PARAMTYPES)
        if value != self.mu:
            self._cached = {}
        self._mu = value
585929
def value(self, data):
        """"""
        value = data.get(self.name)
        if value:
            return int(value)
        return self.default