def parameter_sweep(model, sigma, ns):
    generate_equations(model)
    logp = [numpy.log10(p.value) for p in model.parameters]
    ts = numpy.linspace(0, 20*3600, 20*60)
    solver = Solver(model, ts)
    pf.set_fig_params()
    plt.figure(figsize=(1.8, 1), dpi=300)
    for i in range(ns):
        psample = sample_params(logp, 0.05)
        res = solver.run(param_values=psample)
        signal = res.observables['p53_active']
        plt.plot(signal, color=(0.7, 0.7, 0.7), alpha=0.3)
    # Highlighted
    colors = ['g', 'y', 'c']
    for c in colors:
        psample = sample_params(logp, 0.05)
        res = solver.run(param_values=psample)
        signal = res.observables['p53_active']
        plt.plot(signal, c)

    # Nominal
    solver = Solver(model, ts)
    res = solver.run()
    signal = res.observables['p53_active']
    plt.plot(signal, 'r')

    plt.xticks([])
    plt.xlabel('Time (a.u.)', fontsize=7)
    plt.ylabel('Active p53', fontsize=7)
    plt.yticks([])
    plt.ylim(ymin=0)
    pf.format_axis(plt.gca())
    plt.savefig(model.name + '_sample.pdf')

    def __init__(self, solver, values_to_sample, objective_function,
                 observable):
        self._model = solver.model
        self._logger = get_logger(__name__, model=self._model)
        self._logger.info('%s created for observable %s' % (
            self.__class__.__name__, observable))
        generate_equations(self._model)
        self._ic_params_of_interest_cache = None
        self._values_to_sample = values_to_sample
        if solver is None or not isinstance(solver,
                                            pysb.simulator.base.Simulator):
            raise(TypeError, "solver must be a pysb.simulator object")
        self._solver = solver
        self.objective_function = objective_function
        self.observable = observable

        # Outputs
        self.b_matrix = []
        self.b_prime_matrix = []
        self.p_prime_matrix = np.zeros(self._size_of_matrix)
        self.p_matrix = np.zeros(self._size_of_matrix)

        self._original_initial_conditions = np.zeros(len(self._model.species))
        self._index_of_species_of_interest = self._create_index_of_species()
        self.simulation_initials = self._setup_simulations()
        # Stores the objective function value for the original unperturbed
        # model
        self._objective_fn_standard = None

def test_simres_observable():
    """ Test on demand observable evaluation """
    models = [tyson_oscillator.model, robertson.model,
              expression_observables.model, earm_1_3.model,
              bax_pore_sequential.model, bax_pore.model,
              bngwiki_egfr_simple.model]
    for model in models:
        generate_equations(model)
        spm = SpeciesPatternMatcher(model)
        for obs in model.observables:
            dyn_obs = spm.match(pattern=obs.reaction_pattern, index=True,
                                counts=True)

            # Need to sort by species numerical order for comparison purposes
            dyn_obs = collections.OrderedDict(sorted(dyn_obs.items()))

            dyn_species = list(dyn_obs.keys())

            if obs.match == 'species':
                dyn_coeffs = [1] * len(dyn_obs)
            else:
                dyn_coeffs = list(dyn_obs.values())

            assert dyn_species == obs.species
            assert dyn_coeffs == obs.coefficients

    def __init__(self, model, inputs=None, stop_time=1000,
                 outside_name_map=None):
        self.model = model
        generate_equations(model)

        self.inputs = inputs if inputs else []
        self.stop_time = stop_time
        self.outside_name_map = outside_name_map if outside_name_map else {}

        self.dt = numpy.array(10.0)
        self.units = 'seconds'
        self.sim = None
        self.attributes = copy.copy(default_attributes)
        self.species_name_map = {}
        for idx, species in enumerate(self.model.species):
            monomer = species.monomer_patterns[0].monomer
            self.species_name_map[monomer.name] = idx
        self.input_vars = self._get_input_vars()
        # These attributes are related to the simulation state
        self.state = numpy.array([100.0 for s in self.species_name_map.keys()])
        self.time = numpy.array(0.0)
        self.status = 'start'
        self.time_course = [(self.time, self.state)]
        # EMELI needs a DONE attribute
        self.DONE = False

def stoichiometry_matrix(model):
    generate_equations(model)
    sm = np.zeros((len(model.species), len(model.reactions)))
    for i_s, sp in enumerate(model.species):
        for i_r, r in enumerate(model.reactions):
            sm[i_s][i_r] = r['products'].count(i_s) - r['reactants'].count(i_s)
    return sm

    def setUp(self):
        Monomer('A', ['a'])
        Monomer('B', ['b'])

        Parameter('ksynthA', 100)
        Parameter('ksynthB', 100)
        Parameter('kbindAB', 100)

        Parameter('A_init', 0)
        Parameter('B_init', 0)

        Initial(A(a=None), A_init)
        Initial(B(b=None), B_init)

        Observable("A_free", A(a=None))
        Observable("B_free", B(b=None))
        Observable("AB_complex", A(a=1) % B(b=1))

        Rule('A_synth', None >> A(a=None), ksynthA)
        Rule('B_synth', None >> B(b=None), ksynthB)
        Rule('AB_bind', A(a=None) + B(b=None) >> A(a=1) % B(b=1), kbindAB)

        self.model = model
        generate_equations(self.model)

        # Convenience shortcut for accessing model monomer objects
        self.mon = lambda m: self.model.monomers[m]

        # This timespan is chosen to be enough to trigger a Jacobian evaluation
        # on the various solvers.
        self.time = np.linspace(0, 1)
        self.sim = BngSimulator(self.model, tspan=self.time)

def stoichimetry_matrix_passengers(model, pruned_system):
    generate_equations(model)
    sm = np.zeros((len(pruned_system.keys()), len(model.reactions)))
    for i_s, pa in enumerate(pruned_system.keys()):
        for i_r, r in enumerate(model.reactions):
            if r['rate'] in pruned_system[pa].as_coefficients_dict().keys():
                sm[i_s][i_r] = r['products'].count(pa) - r['reactants'].count(pa)
    return sm

def test_non_python_name_phos():
    st = Phosphorylation(Agent('14-3-3'), Agent('BRAF kinase'))
    pa = PysbAssembler([st])
    pa.make_model()
    names = [m.name for m in pa.model.monomers]
    assert 'BRAF_kinase' in names
    assert 'p14_3_3' in names
    bng.generate_equations(pa.model)

def create_model_files(model, model_name, directory=None):
    curr_dir = os.getcwd()
    if not model.odes:
        generate_equations(model)
    
    if directory != None:
        os.chdir(directory)
    
    file_basename = model_name+'_model'
    dill.dump(model, open(file_basename+'.p', 'wb'))
    
    os.chdir(curr_dir)

def build_all_models():
    mek_seq_rand = ['seq', 'random']
    mkp_seq_rand = ['seq', 'random']
    erk_dimerization = [False, 'any', 'uT', 'phos']
    mkp_activation = [False, True]

    model_combinations = itertools.product(*(mek_seq_rand, mkp_seq_rand,
                                             erk_dimerization, mkp_activation))

    models = {}
    for mc in model_combinations:
        model = build_model(*mc)
        generate_equations(model)
        models[model.name] = model
    return models

    def __init__(self, solver, values_to_sample, objective_function,
                 observable, sens_type='initials', sample_list=None):
        if not isinstance(solver, pysb.simulator.base.Simulator):
            raise TypeError("solver must be a pysb.simulator object")
        self._model = solver.model
        self._logger = get_logger(__name__, model=self._model)
        self._logger.info('%s created for observable %s' % (
            self.__class__.__name__, observable))
        generate_equations(self._model)
        self._values_to_sample = values_to_sample
        self._solver = solver
        self.objective_function = objective_function
        self._observable = observable
        self._sens_type = sens_type
        if self._sens_type not in ('params', 'initials', 'all'):
            if sample_list is None:
                raise ValueError("Please provide 'sens_type' or 'sample_list'")
        if sample_list is not None:
            _valid_options = [i.name for i in self._model.parameters]
            for i in sample_list:
                if i not in _valid_options:
                    raise ValueError("{} not in model.parameters".format(i))
            self.index = sample_list
        elif self._sens_type == 'params':
            self.index = [i.name for i in self._model.parameters_rules()]
        elif self._sens_type == 'initials':
            self.index = [i[1].name for i in self._model.initial_conditions]
        elif self._sens_type == 'all':
            self.index = [i.name for i in self._model.parameters]

        self.orig_vals = [i.value for i in self._model.parameters]
        self.index_of_param = {i.name: n for n, i in
                               enumerate(self._model.parameters)}
        self._n_sam = len(self._values_to_sample)
        self._n_species = len(self.index)
        self._nm = self._n_species * self._n_sam
        self._size_of_matrix = self._nm ** 2
        self._shape_of_matrix = self._nm, self._nm

        # Outputs
        self.b_matrix = []
        self.b_prime_matrix = []
        self.params_to_run = self._setup_simulations()
        self.p_prime_matrix = np.zeros(self._size_of_matrix)
        self.p_matrix = np.zeros(self._size_of_matrix)
        # Stores the objective function value for the original unperturbed
        # model
        self._objective_fn_standard = None

def print_model_stats(model):
    """ provides stats for the models

    :param model:
    """
    generate_equations(model)
    print("Information about model {0}".format(model.name))
    print("Number of rules {0}".format(len(model.rules)))
    print("Number of parameters {0}".format(len(model.parameters)))
    print(
    "Number of parameter rules {0}".format(len(model.parameters_rules())))
    print("Number of reactions {0}".format(len(model.reactions)))
    print(
    "Number of initial conditions {0}".format(len(model.initial_conditions)))
    print("Number of species {0}".format(len(model.species)))
    print('{}'.format('-' * 24))

 def __init__(self, model, tspan=None, initials=None,
              param_values=None, verbose=False, **kwargs):
     super(StochKitSimulator, self).__init__(model,
                                             tspan=tspan,
                                             initials=initials,
                                             param_values=param_values,
                                             verbose=verbose,
                                             **kwargs)
     self.cleanup = kwargs.pop('cleanup', True)
     if kwargs:
         raise ValueError('Unknown keyword argument(s): {}'.format(
             ', '.join(kwargs.keys())
         ))
     self._outdir = None
     generate_equations(self._model,
                        cleanup=self.cleanup,
                        verbose=self.verbose)

def convert_odes(model, p_name_map, s_name_map_by_pattern):
    """Substitutes species and parameter names using the given name mappings.

    Parameters
    ----------
    model : pysb.core.Model
        The model to be converted.
    p_name_map : dict
        A dict where the keys are the parameter names of the PySB model and the
        values are the parameter names of the original model, e.g.
        {'bind_BidT_Bcl2_kf': 'ka_tBid_Bcl2'}
    s_name_map : dict
        A dict where the keys are the string representations of the species
        from the PySB model, generated by calling str(species), and the values
        are the species names in the original model, e.g.
        {'Bid(bf=None, state=T)': 'Act'}

    Returns
    -------
    A list of strings, with one entry for each ODE in the model. Each ODE
    is represented as a string, e.g. "d[Act]/dt = ..."
    """

    generate_equations(model)

    # Get the index of each species
    s_name_map_by_num = {}
    for i, s in enumerate(model.species):
        for key in s_name_map_by_pattern:
            if key == str(s):
                s_name_map_by_num['s%d' % i] = s_name_map_by_pattern[key]

    name_map = {}
    name_map.update(p_name_map)
    name_map.update(s_name_map_by_num)

    # Substitute new names into the ODEs
    ode_list = {} 
    for i, ode in enumerate(model.odes):
        new_ode = ode.subs(name_map)
        ode_species = s_name_map_by_pattern[str(model.species[i])]
        ode_list[ode_species] = str(new_ode)
        #new_ode = 'd[%s]/dt = %s' % (s_name_map_by_num['s%d' % i], str(new_ode))
        #ode_list.append(new_ode)

    return ode_list

def run_one_model(model, data, prior_vals, ns, pool=None):
    # Vector of estimated parameters
    pe = [p for p in model.parameters if p.name.startswith('k')]

    # Generate model equations
    generate_equations(model)
    Solver._use_inline = True
    sol = Solver(model, numpy.linspace(0,10,10))
    sol.run()

    # Number of temperatures, dimensions and walkers
    ntemps = 20
    ndim = len(pe)
    blocksize = 48
    nwalkers = get_num_walkers(ndim, blocksize)
    print 'Running %d walkers at %d temperatures for %d steps.' %\
          (nwalkers, ntemps, ns)

    sampler = emcee.PTSampler(ntemps, nwalkers, ndim, likelihood, prior,
         threads=1, pool=pool, betas=None, a=2.0, Tmax=None,
         loglargs=[model, data], logpargs=[model, prior_vals],
         loglkwargs={}, logpkwargs={})

    # Random initial parameters for walkers
    p0 = numpy.ones((ntemps, nwalkers, ndim))
    for i in range(ntemps):
        for j in range(nwalkers):
            for k, pp in enumerate(pe):
                p0[i, j, k] = prior_vals.vals[j][pp.name]
    print p0
    # Run sampler
    fname = scratch_path + 'chain_%s.dat' % model.name
    step = 0
    for result in sampler.sample(p0, iterations=ns, storechain=True):
        print '---'
        position = result[0]
        with open(fname, 'a') as fh:
            for w in range(nwalkers):
                for t in range(ntemps):
                    pos_str = '\t'.join(['%f' % p for p in position[t][w]])
                    fh.write('%d\t%d\t%d\t%s\n' % (step, w, t, pos_str))
        step += 1
    return sampler

def compare_odes(model, p_name_map, s_index_map, m_ode_list):
    """Compare the PySB model ODEs to the original MATLAB ODEs."""
    s_name_map = [('s%d' % i, 'x(%d)' % j) for i, j in enumerate(s_index_map)]
    name_map = {}
    name_map.update(p_name_map)
    name_map.update(s_name_map)
    generate_equations(model)
    ode_list = []
    result_list = []
    for i, ode in enumerate(model.odes):
        new_ode = ode.subs(name_map)
        old_ode = m_ode_list[s_index_map[i] - 1]
        result = new_ode == old_ode
        result_list.append(result)
        if not result:
            print "Mismatch for species " + str(i) + ": "
            print "Pysb ODE   : %s" % str(new_ode)
            print "Matlab ODE : %s" % str(old_ode)

    return result_list

def check_all_species_generated(model):
    """Check the lists of rules and triggering species match BNG"""
    generate_equations(model)
    spm = SpeciesPatternMatcher(model)

    # Get the rules and species firing them using the pattern matcher
    species_produced = dict()
    for rule, rp_species in spm.rule_firing_species().items():
        species_produced[rule.name] = set()
        for sp_list in rp_species:
            for sp in sp_list:
                species_produced[rule.name].add(
                    model.get_species_index(sp))

    # Get the equivalent dictionary of {rule name: [reactant species]} from BNG
    species_bng = collections.defaultdict(set)
    for rxn in model.reactions:
        for rule in rxn['rule']:
            species_bng[rule].update(rxn['reactants'])

    # Our output should match BNG
    assert species_produced == species_bng

def test_species_pattern_matcher():
    # See also SpeciesPatternMatcher doctests

    # Check that SpeciesPatternMatcher raises exception if model has no species
    model = robertson.model
    model.reset_equations()
    assert_raises(Exception, SpeciesPatternMatcher, model)

    model = bax_pore.model
    generate_equations(model)
    spm = SpeciesPatternMatcher(model)
    BAX = model.monomers['BAX']
    sp_sets = spm.species_fired_by_reactant_pattern(
        as_reaction_pattern(BAX(t1=None, t2=None))
    )
    assert len(sp_sets) == 1
    assert len(sp_sets[0]) == 2

    sp_sets = spm.species_fired_by_reactant_pattern(
        as_reaction_pattern(BAX(t1=WILD, t2=ANY))
    )
    assert len(sp_sets) == 1
    assert len(sp_sets[0]) == 10

def run_cupsoda(tspan, model, simulations, vol, name, mem):
    tpb = 32

    solver = CupSodaSimulator(model, tspan, atol=ATOL, rtol=RTOL, gpu=0,
                              max_steps=mxstep, verbose=False, vol=vol,
                              obs_species_only=True,
                              memory_usage=mem_dict[mem])

    cols = ['model', 'nsims', 'tpb', 'mem', 'cupsodatime', 'cupsoda_io_time',
            'pythontime', 'rtol', 'atol', 'mxsteps', 't_end', 'n_steps', 'vol',
            'card']
    generate_equations(model)
    nominal_values = np.array([p.value for p in model.parameters])
    all_output = []
    for num_particles in simulations:
        # set the number of blocks to make it 16 threads per block
        n_blocks = int(np.ceil(1. * num_particles / tpb))
        solver.n_blocks = n_blocks

        # create a matrix of initial conditions
        c_matrix = np.zeros((num_particles, len(nominal_values)))
        c_matrix[:, :] = nominal_values
        y0 = np.zeros((num_particles, len(model.species)))
        for ic in model.initial_conditions:
            for j in range(len(model.species)):
                if str(ic[0]) == str(model.species[j]):
                    y0[:, j] = ic[1].value
                    break

        # setup a unique log output file to extract timing from
        log_file = 'logfile_{}.log'.format(num_particles)
        # remove it if it already exists (happens if rerunning simulation)
        if os.path.exists(log_file):
            os.remove(log_file)

        setup_logger(logging.INFO, file_output=log_file, console_output=True)
        start_time = time.time()

        # run the simulations
        x = solver.run(param_values=c_matrix, initials=y0)

        end_time = time.time()

        # create an emtpy list to put all data for this run
        out_list = list()
        out_list.append(name)
        out_list.append(num_particles)
        out_list.append(str(tpb))
        out_list.append(mem)
        cupsoda_time = 'error'
        total_time = 'error'
        with open(log_file, 'r') as f:
            for line in f:
                if 'reported time' in line:
                    good_line = line.split(':')[-1].split()[0]
                    cupsoda_time = float(good_line)
                if 'I/O time' in line:
                    good_line = line.split(':')[-1].split()[0]
                    total_time = float(good_line)
            f.close()

        out_list.append(cupsoda_time)
        out_list.append(total_time)
        out_list.append(end_time - start_time)
        out_list.append(RTOL)
        out_list.append(ATOL)
        out_list.append(len(tspan))
        out_list.append(mxstep)
        out_list.append(np.max(tspan))
        out_list.append(vol)
        out_list.append(card)
        all_output.append(out_list)
        print(" ".join(str(i) for i in out_list))
        print('out==')
        print(x.observables[0][0])
        print(x.observables[0][-1])
        print('==out\n')

    df = pd.DataFrame(all_output, columns=cols)
    print(df)
    df.to_csv('{}_cupsoda_timings_{}.csv'.format(name, mem), index=False)
    return df

    # Add tBid initial condition to mito only model
    Bid = m_mito.all_components()["Bid"]
    try:
        p = Parameter("tBid_0", 1, _export=False)
        m_mito.initial(Bid(state="T", bf=None), p)
        m_mito.add_component(p)
    except Exception:
        pass  # Duplicate initial condition for tBid

    if model in ["chen_2007_febs_direct", "howells"]:
        Bad = m_mito.all_components()["Bad"]
        p = Parameter("mBad_0", 1, _export=False)
        m_mito.initial(Bad(state="M", bf=None, serine="U"), p)
        m_mito.add_component(p)

    generate_equations(m_mito)
    generate_equations(m_full)
    print(
        "%-24s %11d %10d %12d %11d %10d %12d"
        % (
            model,
            len(m_mito.rules),
            len(m_mito.odes),
            len(m_mito.parameters),
            len(m_full.rules),
            len(m_full.odes),
            len(m_full.parameters),
        )
    )

# Full models