def bistability_analysis():
    f2_range = linspace(0, 0.4, 41)
    t = linspace(0, 50000, 1000)
    ion()
    ss_aBax_vals_up = []
    ss_aBax_vals_down = []

    for f2 in f2_range:
        model.parameters['Bid_0'].value = f2 * 1e-1
        bax_total = 2e-1

        # Do "up" portion of hysteresis plot
        model.parameters['aBax_0'].value = 0
        model.parameters['cBax_0'].value = bax_total
        x = odesolve(model, t)
        figure('up')
        plot(t, x['aBax_']/bax_total)
        ss_aBax_vals_up.append(x['aBax_'][-1]/bax_total)

        # Do "down" portion of hysteresis plot
        model.parameters['aBax_0'].value = bax_total
        model.parameters['cBax_0'].value = 0
        x = odesolve(model, t)
        figure('down')
        plot(t, x['aBax_']/bax_total)
        ss_aBax_vals_down.append(x['aBax_'][-1]/bax_total)

    figure()
    plot(f2_range, ss_aBax_vals_up, 'r')
    plot(f2_range, ss_aBax_vals_down, 'g')

def time_dep_sensitivity(model, t, results = None, delta=0.1, parameters=[]):
	if not parameters:
		return
	all_param_results = {}

	if results is None:
		results = odesolve(model, t)
	print 'Solving time-dependent sensitivity analysis for %d parameters' % len(parameters)

	for param in parameters:
		if param.name.startswith('__'):
			continue
		print 'Solving for', param.name
		
		oldval = param.value

		newval_left = param.value * (1. + 0.5 * delta)
		model.parameters[param.name].value = newval_left
		new_results_left = odesolve(model, t)

		newval_right = param.value * (1. - 0.5 * delta)
		model.parameters[param.name].value = newval_right
		new_results_right = odesolve(model, t)

		model.parameters[param.name].value = oldval

		all_param_results[param.name] = sens_res(param.name,
												 new_results_left,
												 new_results_right,
												 delta * param.value)
	
	return all_param_results

def plot_effect_of_pore_reverse_rate():
    ci = color_iter()

    pore_reverse_rates = [1e-2, 1e-4, 1e-6]
    t = np.linspace(0, 5000, 500)
    plt.figure()

    for pore_reverse_rate in pore_reverse_rates:
        params_dict = {'Bax_0': 50., 'Vesicles_0': 50.,
                        'pore_reverse_rate_k': pore_reverse_rate}
        b = one_cpt.Builder(params_dict=params_dict)
        b.build_model_bax_schwarz_reversible()
        x = odesolve(b.model, t)
        avg_pores = x['pores']/b.model.parameters['Vesicles_0'].value
        col = ci.next()
        plt.plot(t, avg_pores, color=col, linestyle='-',
                 label="Pores, reverse rate %g" % pore_reverse_rate)
        plt.plot(t, 1 - np.exp(-avg_pores), color=col, linestyle='--',
                 label="Dye release, reverse rate %g" % pore_reverse_rate)

    plt.legend(loc='upper left')
    plt.xlabel('Time (sec)')
    plt.ylabel('Dye release/Avg. pores')
    plt.title('Dye release/pore formation with varying reverse rates')
    plt.show()

def plot_fraction_bax_bound(model, figure_id=None):
    bax_concs = np.logspace(0, 4, 40)
    t = np.linspace(0, 3000, 1000)

    #plt.figure()
    ss_mBax_values = []
    for bax_conc in bax_concs:
        model.parameters['Bax_0'].value = bax_conc
        x = odesolve(model, t)
        mBax_frac = x['mBax'] / model.parameters['Bax_0'].value
        ss_mBax_value = mBax_frac[-1]
        ss_mBax_values.append(ss_mBax_value)
        #plt.plot(t, mBax_frac)
    #plt.show()

    ss_mBax_values = np.array(ss_mBax_values)

    if figure_id is not None:
        plt.figure(figure_id)
    else:
        plt.figure()

    plt.semilogx(bax_concs, ss_mBax_values, linewidth=2)
    plt.xlabel('Bax concentration (nM)')
    plt.ylabel('Pct. Bax at liposomes')
    plt.ylim([0, 1])
    plt.title('Frac Bax bound vs. Bax conc')
    plt.show()

    def sensitivity(self,tspan=None, parameters_ref=None, sp_SA=None, sampling_method = 'saltelli', analyze_method='sobol',N=1 , verbose=True):
      
        ana_meth = importlib.import_module('SALib.analyze.%s'%analyze_method)
      
        if tspan is not None:
            self.tspan = tspan
        elif self.tspan is None:
            raise Exception("'time t' must be defined.")
        
        if sp_SA is None:
            raise Exception("A species to do the sensitivity analysis on must be defined")
        if sp_SA not in [str(sp) for sp in self.model.species] and sp_SA  not in [str(obs.name) for obs in self.model.observables]:
            raise Exception("Species is not in model species")

        ref = odesolve(self.model, self.tspan, param_values=parameters_ref)

        if verbose: print "Getting parameters information"
        self.pars_info(parameters_ref, sampling_method=sampling_method, N=N)
        
        if verbose: print "Simulating with parameters from sampling"
        p=multiprocessing.Pool(2)
        func = partial(self.sim_model,sp_SA=sp_SA)
        self.sim = p.map(func, self.param_sets)
        self.output = [sum((r - ref[sp_SA])**2) for r in self.sim]
          
        if verbose: print "Sensitivity analysis"
        self.sa_result = ana_meth.analyze(self.problem, numpy.array(self.output), print_to_console=True)

def test_integrate_with_expression():

    Monomer('s1')
    Monomer('s9')
    Monomer('s16')
    Monomer('s20')

    Parameter('ka',2e-5)
    Parameter('ka20', 1e5)

    Initial(s9(), Parameter('s9_0', 10000))

    Observable('s1_obs', s1())
    Observable('s9_obs', s9())
    Observable('s16_obs', s16())
    Observable('s20_obs', s20())

    Expression('keff', (ka*ka20)/(ka20+s9_obs))

    Rule('R1', None >> s16(), ka)
    Rule('R2', None >> s20(), ka)
    Rule('R3', s16() + s20() >> s16() + s1(), keff)

    time = linspace(0, 40, 100)
    x = odesolve(model, time)

def test_simulation():
    time = np.linspace(0, 18000, 1801)
    x = odesolve(model, time)
    expectedx_final = {'IkBan_obs': 3145.8443053128276,
                       'NFkBn_obs': 178.06326194655537,
                       'IkBap_NFkBc_obs': 1.3899039187480429e-07,
                       'IkBa_on_obs': 0.045201392272254615,
                       'IKKKa_obs': 0.00054704675357157236,
                       'A20_off_obs': 1.9547986077277453,
                       'A20t_obs': 6.2762670477040663,
                       'IkBan_NFkBn_obs': 28.012691613691331,
                       'IKKn_obs': 199999.9951281352,
                       'IKKa_obs': 2.3395704210536575e-08,
                       'TNFR1a_obs': 5.3362282855621739e-06,
                       'IkBa_off_obs': 1.9547986077277459,
                       'TNF_ext_obs': 9.5649122910157076e-08,
                       'A20_on_obs': 0.045201392272254573,
                       'IkBap_obs': 2.3942600622297971e-09,
                       'IKKii_obs': 0.0045466692017435365,
                       'IKKKn_obs': 99999.999452953241,
                       'IKKi_obs': 0.00032517206132429724,
                       'TNFR1i_obs': 1999.9999946637713,
                       'IkBat_obs': 6.2762670477040681,
                       'NFkBc_obs': 139.40965738805579,
                       'IkBac_obs': 8607.2535257899253,
                       'A20_obs': 5714.7212868294837,
                       'IkBa_NFkB_obs': 99655.514388905125}
    for i in expectedx_final.keys():
        final_obs = x[i][-1]
        if abs(final_obs - expectedx_final[i]) > 1e-3:
            raise ValueError("%s: Expected %f, got %f, diff %e" %
                            (i, expectedx_final[i], final_obs, final_obs -
                             expectedx_final[i]))

    def tropicalize(self, tspan=None, param_values=None, diff_par=1, ignore=1, epsilon=1, find_passengers_by='imp_nodes',
                    plot_imposed_trace=False, verbose=False):
        """
        tropicalization of driver species
        :param diff_par:
        :param find_passengers_by: Option to find passenger species. 'imp_nodes' finds the nodes that only have one edge.
        'qssa' finds passenger species using the quasi steady state approach
        :param plot_imposed_trace: Option to plot imposed trace
        :param tspan: Time span
        :param param_values: PySB model parameter values
        :param ignore: Initial time points to ignore
        :param epsilon: Order of magnitude difference between solution of ODE and imposed trace to consider species as
        passenger
        :param verbose: Verbose
        :return:
        """

        if verbose:
            print("Solving Simulation")

        if tspan is not None:
            self.tspan = tspan
        else:
            raise SimulatorException("'tspan' must be defined.")

        if param_values is not None:
            # accept vector of parameter values as an argument
            if len(param_values) != len(self.model.parameters):
                raise Exception("param_values must be the same length as model.parameters")
            if not isinstance(param_values, numpy.ndarray):
                param_values = numpy.array(param_values)
        else:
            # create parameter vector from the values in the model
            param_values = numpy.array([p.value for p in self.model.parameters])

        # convert model parameters into dictionary
        self.param_values = dict((p.name, param_values[i]) for i, p in enumerate(self.model.parameters))

        self.y = odesolve(self.model, self.tspan, self.param_values)

        if verbose:
            print("Getting Passenger species")
        if find_passengers_by == 'qssa':
            if plot_imposed_trace:
                self.find_passengers(self.y[ignore:], epsilon, plot=plot_imposed_trace)
            else:
                self.find_passengers(self.y[ignore:], epsilon)
        elif find_passengers_by == 'imp_nodes':
            self.find_nonimportant_nodes()
        else:
            raise Exception("equations to tropicalize must be chosen")

        if verbose:
            print("equation to tropicalize")
        self.equations_to_tropicalize()

        if verbose:
            print("Getting signatures")
        self.signal_signature(self.y[ignore:], diff_par=diff_par)
        return

def plot_liposome_titration(builder=one_cpt.Builder()):
    b = builder
    b.translocate_Bax()
    b.model.parameters['Bax_0'].value = 100
    t = np.linspace(0, 3000, 1000)
    lipo_concs = np.logspace(-2, 4, 40)
    ss_mBax_values = []

    # Plot timecourses
    plt.ion()
    #plt.figure()

    for lipo_conc in lipo_concs:
        b.model.parameters['Vesicles_0'].value = lipo_conc
        x = odesolve(b.model, t)
        mBax_frac = x['mBax'] / b.model.parameters['Bax_0'].value
        plt.plot(t, mBax_frac)
        ss_mBax_value = mBax_frac[-1]
        ss_mBax_values.append(ss_mBax_value)
    plt.show()
    ss_mBax_values = np.array(ss_mBax_values)

    # Plot liposome titration
    plt.figure()
    plt.plot(lipo_concs, ss_mBax_values)
    plt.xlabel('Liposome concentration (nM)')
    plt.ylabel('Pct. Bax at liposomes (Bax_0 = %d nM)' %
               b.model.parameters['Bax_0'].value)
    plt.title('Bax/Liposome binding curve')
    plt.show()

def run_model(tmax=12000, fittype='explin'):
    t = linspace(0, tmax, 1000)
    x = odesolve(model, t)
    figure()
    plot(t, (x['ctBid'])/tBid_0.value, label='ctBid')
    plot(t, (x['mtBid'])/tBid_0.value, label='mtBid')
    plot(t, (x['cBax'])/Bax_0.value, label='cBax')
    plot(t, (x['mBax'])/Bax_0.value, label='mBax')
    #plot(t, x['metBid']/tBid_0.value, label='metBid')
    #plot(t, x['mtBid']/tBid_0.value, label='mtBid')
    #plot(t, x['mftBid']/tBid_0.value, label='mftBid')
    #plot(t, x['mfBax']/Bax_0.value, label='mfBax')
    #plot(t, x['meBax']/Bax_0.value, label='meBax')
    #plot(t, x['tBidBax']/Bax_0.value, label='tBidBax')
    #plot(t, x['tBidiBax']/Bax_0.value, label='tBidiBax')
    #plot(t, (x['iBax']+x['pBax'])/Bax_0.value, label='ipBax')
    #plot(t, (x['iBax'])/Bax_0.value, label='iBax')
    #plot(t, (x['eiBax'])/Bax_0.value, label='eiBax')
    #plot(t, (2*x['ePore'])/Bax_0.value, label='ePore')
    #plot(t, (2*x['fPore'])/Bax_0.value, label='fPore')
    #plot(t, (x['pBax'])/Bax_0.value, label='pBax')
    #plot(t, (2*(x['Bax2']+(2*x['Bax4'])))/Bax_0.value, label='Bax2')
    #plot(t, (2*x['Bax2'])/Bax_0.value, label='Bax2')
    #plot(t, (4*x['Bax4'])/Bax_0.value, label='Bax4')
    #plot(t, x['pBax']/Bax_0.value, label='pBax')
    #plot(t, x['eVes']/Vesicles_0.value, label='eVes')
    legend()
    xlabel("Time (seconds)")
    ylabel("Normalized Concentration")

def test_build_models():
    model_names = [#'t', 'ta', 'tai', 'taid', 'taidt', 'tair', 'taird',
                   #'tairdt', 'tad', 'tadt', 'tar', 'tard', 'tardt',
                    #---
                   #'bax_heat',
                   #'bax_heat_reversible',
                   #'bax_heat_dimer',
                   #'bax_heat_dimer_reversible',
                   #'bax_heat_auto1',
                   #'bax_heat_auto2',
                    #---
                   #'bax_heat_bh3_auto2',
                   #'bax_heat_auto1_reversible_activation',
                   #'bax_heat_auto2_reversible_activation',
                   #'bax_heat_auto1_reversible',
                   #'bax_heat_auto2_reversible',
                   #'bax_heat_auto1_dimer',
                   #'bax_heat_auto2_dimer',
                   #'bax_heat_auto1_dimer_reversible',
                   #'bax_heat_auto2_dimer_reversible',
                   #'bax_heat_bh3_exposure_auto2',
                   #'bax_schwarz',
                   #'bax_schwarz_reversible',
                   #'bax_schwarz_dimer',
                   #'bax_schwarz_dimer_reversible',
                   #'bax_schwarz_tetramer_reversible',
    ]

    for model_name in model_names:
        print "Running model %s" % model_name
        b = Builder()
        eval("b.build_model_%s()" % model_name)
        x = odesolve(b.model, t)
        plot_model(t, x, 'build_model_%s' % model_name)

def test_integrate_with_expression():
    """Ensure a model with Expressions simulates."""

    Monomer('s1')
    Monomer('s9')
    Monomer('s16')
    Monomer('s20')

    # Parameters should be able to contain s(\d+) without error
    Parameter('ks0',2e-5)
    Parameter('ka20', 1e5)

    Initial(s9(), Parameter('s9_0', 10000))

    Observable('s1_obs', s1())
    Observable('s9_obs', s9())
    Observable('s16_obs', s16())
    Observable('s20_obs', s20())

    Expression('keff', (ks0*ka20)/(ka20+s9_obs))

    Rule('R1', None >> s16(), ks0)
    Rule('R2', None >> s20(), ks0)
    Rule('R3', s16() + s20() >> s16() + s1(), keff)

    time = np.linspace(0, 40)
    x = odesolve(model, time)

def plot_model_data():
    """Plots a simulation of the model along with the data."""
    plt.ion()
    x = odesolve(model, tspan)
    plt.plot(tspan, x['A_'])
    plt.plot(tspan, synthetic_data['A_'])
    plt.show()

def run_figure_sim(model):
    """Run the C8 dose-response series shown in Fig. 11 of [Albeck2008]_.

    Returns
    -------
    [t, outputs] : list containing two numpy.array objects
        t: The time coordinates of each timepoint, in seconds, from 0 to
            60,000 seconds (15 hours), at 60-second intervals.
        outputs: A 901 x 7 array. Each row corresponds to the fraction of Smac
            released into the cytosol at each timepoint. Each column correspond
            to a distinct caspase 8 dose, from lowest to highest:
            [1, 5, 10, 50, 100, 500, 1000], in molecules per cell.
    """
    # Set timepoints and c8 doses
    tf = 15 * 3600 # 15 hours
    t = np.linspace(0, tf, (tf / 60) + 1)
    c8_doses = [0.01e2, 0.05e2, 0.1e2, 0.5e2, 1e2, 5e2, 10e2]

    outputs = np.empty((len(t), len(c8_doses)))
    for i, c8_dose in enumerate(c8_doses):
        model.parameters['C8_0'].value = c8_dose
        x = odesolve(model, t, rtol=rtol)
        frac_Smac_release = x['cSmac'] / model.parameters['Smac_0'].value
        outputs[:,i] = frac_Smac_release

    return [t, outputs]

def main():
    t = linspace(0, 60)
    y = pint.odesolve(model, t)
    print y['A_total'][-1], y['B_total'][-1]
    plot(t, y['A_total'])
    plot(t, y['B_total'])
    legend([o.name for o in model.observables])
    show()

def synthetic_data(model, tspan, obs_list=None, sigma=0.1, seed=None):
    from pysb.integrate import odesolve

    random = numpy.random.RandomState(seed)
    ysim = odesolve(model, tspan)
    ysim_array = ysim.view().reshape(len(tspan), len(ysim.dtype))
    ysim_array *= random.randn(*ysim_array.shape) * sigma + 1
    return ysim

    def visualize(self, fig_path='', tspan=None, parameters=None, verbose=False):
        if verbose:
            print("Solving Simulation")

        if tspan is not None:
            self.tspan = tspan
        elif self.tspan is None:
            raise Exception("'tspan' must be defined.")

        if parameters is not None:
            # accept vector of parameter values as an argument
            if len(parameters) != len(self.model.parameters):
                raise Exception("parameters must be the same length as model.parameters")
            if not isinstance(parameters, numpy.ndarray):
                parameters = numpy.array(parameters)
        else:
            # create parameter vector from the values in the model
            parameters = numpy.array([p.value for p in self.model.parameters])

        new_pars = dict((p.name, parameters[i]) for i, p in enumerate(self.model.parameters))
        self.parameters = new_pars

        self.y = odesolve(self.model, self.tspan, self.parameters)

        if verbose:
            print("Creating graph")
        self.species_graph()
        self.sp_graph.layout(prog='dot',
                             args="-Gstart=50 -Gesep=1.5 -Gnodesep=0.8 -Gsep=0.5  -Gsplines=true -Gsize=30.75,10.75\! "
                                  "-Gratio=fill -Grankdir=LR -Gdpi=200! -Gordering=in")
        self.sp_graph.add_node('t',
                               label='time',
                               shape='oval',
                               fillcolor='white', style="filled", color="transparent",
                               fontsize="50",
                               margin="0,0",
                               pos="20,20!")

        self.edges_colors(self.y)
        self.nodes_colors(self.y)

        if os.path.exists(fig_path):
            directory = fig_path
        else:
            directory = os.getcwd() + '/visualizations'
            os.makedirs(directory)

        if verbose:
            "Generating images"
        for kx, time in enumerate(self.tspan):
            self.sp_graph.get_node('t').attr['label'] = 'time:' + ' ' + '%d' % time + ' ' + 'sec'
            map(functools.partial(change_edge_colors, graph=self.sp_graph), list(self.colors_time_edges.index),
                list(self.colors_time_edges.iloc[:, kx]))
            map(functools.partial(change_node_colors, graph=self.sp_graph), list(self.colors_time_nodes.index),
                list(self.colors_time_nodes.iloc[:, kx]))
            self.sp_graph.draw(directory + '/file' + '%03d' % kx + '.png')

def plot_pores_and_efflux(model):
    t = np.linspace(0, 5000, 500)
    x = odesolve(model, t)
    plt.figure()
    avg_pores = x['pores']/model.parameters['Vesicles_0'].value
    plt.plot(t, avg_pores, 'r', label="Pores/Ves")
    plt.plot(t, 1 - np.exp(-avg_pores), 'b', label="Dye release")
    plt.legend(loc='lower right')
    plt.xlabel('Time (sec)')
    plt.show()

    def tropicalize(self, tspan=None, param_values=None, ignore=1, epsilon=1, plot_imposed_trace=False, verbose=False):
        """
        tropicalization of driver species
        :param plot_imposed_trace: Option to plot imposed trace
        :param tspan: Time span
        :param param_values: PySB model parameter values
        :param ignore: Initial time points to ignore
        :param epsilon: Order of magnitude difference between solution of ODE and imposed trace to consider species as
        passenger
        :param verbose: Verbose
        :return:
        """

        if verbose:
            print("Solving Simulation")

        if tspan is not None:
            self.tspan = tspan
        else:
            raise SimulatorException("'tspan' must be defined.")

        if param_values is not None:
            # accept vector of parameter values as an argument
            if len(param_values) != len(self.model.parameters):
                raise Exception("param_values must be the same length as model.parameters")
            if not isinstance(param_values, numpy.ndarray):
                param_values = numpy.array(param_values)
        else:
            # create parameter vector from the values in the model
            param_values = numpy.array([p.value for p in self.model.parameters])

        # convert model parameters into dictionary
        self.param_values = dict((p.name, param_values[i]) for i, p in enumerate(self.model.parameters))

        self.y = odesolve(self.model, self.tspan, self.param_values)

        if verbose:
            print("Getting Passenger species")
        if plot_imposed_trace:
            self.find_passengers(self.y[ignore:], epsilon, plot=plot_imposed_trace)
        else:
            self.find_passengers(self.y[ignore:], epsilon)

        if verbose:
            print("equation to tropicalize")
        self.equations_to_tropicalize()

        if verbose:
            print("Getting tropicalized equations")
        self.final_tropicalization()
        self.data_drivers(self.y[ignore:])
        return

def run_model_ode(tmax=12000):
    from pysb.integrate import odesolve
    t = np.linspace(0, tmax, 1000)
    x = odesolve(model, t)

    ci = color_iter()
    plt.figure(1)
    # Translocation
    plt.plot(t, (x['ctBid'])/tBid_0.value, label='ctBid', color=ci.next())
    plt.plot(t, (x['mtBid'])/tBid_0.value, label='mtBid', color=ci.next())
    plt.plot(t, (x['cBax'])/Bax_0.value, label='cBax', color=ci.next())
    plt.plot(t, (x['mBax'])/Bax_0.value, label='mBax', color=ci.next())
    plt.legend()