def triangleMAPs(savefilename,basename):
    with open(savefilename,'rb') as savefile:
        bf= numpy.array(pickle.load(savefile))
        samples= numpy.array(pickle.load(savefile))
        bf_g15= numpy.array(pickle.load(savefile))
        samples_g15= numpy.array(pickle.load(savefile))
        bf_zero= numpy.array(pickle.load(savefile))
        samples_zero= numpy.array(pickle.load(savefile))
    labels= []
    for jj in range(samples.shape[2]):
        labels.append(r"$\mathrm{param}\ %i$" % jj)
    maps= define_rcsample.MAPs()
    for ii, map in enumerate(maps.map()):
        if ii >= len(bf): break
        tfeh= numpy.nanmedian(map['FE_H'])
        tafe= numpy.nanmedian(map[define_rcsample._AFETAG])
        for tbf,tsamples,ext in zip([bf,bf_g15,bf_zero],
                                    [samples,samples_g15,samples_zero],
                                    ['fid','g15','zero']):
            try:
                triangle.corner(tsamples[ii,].T,quantiles=[0.16, 0.5, 0.84],
                                labels=labels,
                                show_titles=True,title_args={"fontsize": 12},
                                bins=21)
            except ValueError: pass
            else:
                bovy_plot.bovy_text(r'$[\mathrm{{Fe/H}}] = {feh:.1f},$'\
                                        .format(feh=tfeh)+'\n'
                                    +r'$[\alpha/\mathrm{{Fe}}] = {afe:.2f}$'\
                                        .format(afe=tafe),
                                    top_left=True,size=16.)
                bovy_plot.bovy_end_print(basename+"_%i_%s.png" % (ii,ext))
    return None

def plot_corner(logpost, chain, outdir=None):
    flatchain = chain.reshape((-1, chain.shape[2]))

    tri.corner(flatchain, labels=logpost.pnames, quantiles=[0.05, 0.95])

    if outdir is not None:
        pp.savefig(op.join(outdir, 'corner.pdf'))

 def plot_triangle(self, **kwargs):
     try:
         from triangle import corner
     except ImportError:
         raise ImportError('Plotting requires trianglepy')
     corner(self.posterior_data[self.parameter_names],
            labels=self.parameter_names, **kwargs)

def MCParTemp(nwalkers,steps,ntemps=40,burnin=100,threads=1):
    t1 = time.time()
    # use emcee parallel tempering sampler to find prob. surface
    ndim = 4
    startpos = [[[0,0,0,0] + 5e-3*np.random.randn(ndim) for i in range(nwalkers)] for i in range(ntemps)]

    # time how long it should take
    t0 = time.time()
    testpos = [[[0,0,0,0] + 5e-3*np.random.randn(4) for i in range(8)] for i in range(8)]
    sampler  = emcee.PTSampler(8,8, 4, lnlike, xs.lnprior, threads=threads)
    for ptest,logprobtest, logliketest in sampler.sample(testpos,iterations=10):
        pass
    sampler.reset()
    tdiff = time.time() - t0
    esttime = ntemps/8.*nwalkers/8.*(burnin+steps)/10.*tdiff/60.
    print 'Estimated time = {} minutes'.format(esttime)

    # burn-in for defined number of steps (default=100)
    sampler  = emcee.PTSampler(ntemps,nwalkers, ndim, lnlike, xs.lnprior, threads=threads)
    for p,logprob, loglike in sampler.sample(startpos,iterations=burnin):
        pass
    sampler.reset()

    # run mcmc
    sampler = emcee.PTSampler(ntemps, nwalkers, ndim, lnlike, xs.lnprior, threads=threads)
    for p, logprob, loglike in sampler.sample(p, lnprob0=logprob, lnlike0=loglike,iterations=steps):
        pass

    print 'Time: {} minutes'.format(time.time() - t0)
    samples = sampler.chain[:,:,:,:].reshape((-1, ndim))
    triangle.corner(samples,labels=[r'$\Delta s$',r'$M_A$',r'$F_1^S$',r'$\mu_S$'],quantiles=[0.05,0.5,0.95])
    print sampler.thermodynamic_integration_log_evidence(fburnin=0)
    return samples

def plot_mcmc_results(chain):
    # Pull m and b arrays out of the Markov chain.
    mm = [m for b,m in chain]
    bb = [b for b,m in chain]
    # Scatterplot of m,b posterior samples
    plt.clf()
    plt.contour(bgrid, mgrid, posterior, pdf_contour_levels(posterior))
    plt.plot(bb, mm, 'b.', alpha=0.1)
    plot_mb_setup()
    plt.show()
    
    # Histograms
    import triangle
    triangle.corner(chain, labels=['b','m'], extents=[0.99]*2)
    plt.show()
    
    # Traces
    plt.clf()
    plt.subplot(2,1,1)
    plt.plot(mm, 'k-')
    plt.ylim(mlo,mhi)
    plt.ylabel('m')
    plt.subplot(2,1,2)
    plt.plot(bb, 'k-')
    plt.ylabel('b')
    plt.ylim(blo,bhi)
    plt.show()

    def plotCorner(self):
        model = self.model

        mcmcVParams = self.mcmcRes.vparam_names
        nestVParams = self.nestRes.vparam_names

        mcmcSamples = self.mcmcRes.samples
        nestSamples = self.nestRes.samples

        mcmc_ndim, mcmc_nsamples = len(mcmcVParams), len(mcmcSamples)
        nest_ndim, nest_nsamples = len(nestVParams), len(nestSamples)

        # make figure

        figure_mcmc = triangle.corner(mcmcSamples, labels=[mcmcVParams[0], mcmcVParams[1], mcmcVParams[2], mcmcVParams[3]],
                     truths=[model.get(mcmcVParams[0]), model.get(mcmcVParams[1]),
                             model.get(mcmcVParams[2]), model.get(mcmcVParams[3])],
                     range=mcmc_ndim*[0.9999],
                     show_titles=True, title_args={"fontsize": 12})

        figure_mcmc.gca().annotate("mcmc sampling", xy=(0.5, 1.0), xycoords="figure fraction",
                  xytext=(0, -5), textcoords="offset points",
                  ha="center", va="top")

        figure_nest = triangle.corner(nestSamples, labels=[nestVParams[0], nestVParams[1], nestVParams[2], nestVParams[3]],
                     truths=[model.get(nestVParams[0]), model.get(nestVParams[1]),
                             model.get(nestVParams[2]), model.get(nestVParams[3])],
                     weights=self.nestRes.weights, range=nest_ndim*[0.9999],
                     show_titles=True, title_args={"fontsize": 12})

        figure_nest.gca().annotate("nest sampling", xy=(0.5, 1.0), xycoords="figure fraction",
                  xytext=(0, -5), textcoords="offset points",
                  ha="center", va="top")

        return figure_mcmc, figure_nest

	def plot_PDFtriangle(self,parameterset, labels):		

		if parameterset=='10pars':
			figure = triangle.corner(self.chain.flatchain, labels= labels, plot_contours=True, plot_datapoints = False, show_titles=True, quantiles=[0.16, 0.50, 0.84])
		elif parameterset == 'int_lums':
			figure = triangle.corner(self.int_lums.T, labels= labels,   plot_contours=True, plot_datapoints = False, show_titles=True, quantiles=[0.16, 0.50, 0.84])
		return figure

 def plot_triangle(self, **kwargs):
     try:
         from triangle import corner
     except ImportError:
         raise ImportError('Plotting requires trianglepy')
     data_columns = self.posterior_data.columns[2:]
     corner(self.posterior_data[data_columns],
            weights=self.posterior_data.posterior, **kwargs)

def mcmctriangle(ID):
	filename=('tables/ndim_' + str(ndim) + '_walkers_' + str(nwalkers) + '.fits')
	data = Table.read(filename)
	data_t = np.array([data[key] for key in param_keys]).transpose()
	truths = [np.median(data[key]) for key in param_keys]
	triangle.corner(data_t,labels=param_labels,quantiles=[0.1587,0.5000,0.8413])
	plt.savefig('plots/triangle_ndim_' + str(ndim) + '_walkers_' + str(nwalkers) + '.pdf')
	plt.close()

def convergence_plots(sampler):
    plt.figure()
    plt.plot(sampler.lnprobability.T)
    plt.figure()
    pu.plot_emcee_chains(sampler.chain)
    triangle.corner(sampler.flatchain)

    print 'Autocorrelation lengths: ', ac.emcee_chain_autocorrelation_lengths(sampler.chain)
    print 'Gelman-Rubin R: ', ac.emcee_gelman_rubin_r(sampler.chain)

def PlotTriangle(fileroot,usetruths=True):
	data = np.loadtxt(fileroot+'post_equal_weights.dat', usecols=(0,1,2,3,4,5))
	if (usetruths):
		truths = np.loadtxt(fileroot+'injected_values.txt')
		figure = triangle.corner(data, labels=names, truths=truths, bins=30, quantiles=[0.05, 0.5, 0.95],
								 show_titles=True, title_args={"fontsize": 12})
	else:
		figure = triangle.corner(data, labels=names, bins=30, quantiles=[0.05, 0.5, 0.95],
								 show_titles=True, title_args={"fontsize": 12})
	figure.savefig(fileroot+'posterior_plot.png')

def plot_triangle(samples):
    import triangle
    if samples.shape[1] == 2:
        fig = triangle.corner(samples, labels=["$t_0$", "$P$"])
    if samples.shape[1] == 6:
        fig = triangle.corner(samples, labels=["$t_0$", r"depth", r"duration",
                                               r"$b$", "$q_1$", "$q_2$"])
    elif samples.shape[1] == 8:
        fig = triangle.corner(samples, labels=["$t_0$", r"depth", r"duration",
                                               r"$b$", "$q_1$", "$q_2$", "$q_3$", "$q_4$"])
    plt.show()

def triangle_plot_chain(chain, lnprob, prefix, truth=None, truth_lnprob=0.):
    """
    Make a 7x7 triangle.
    """
    nx, nq = chain.shape
    maxlnp = np.max(lnprob)
    lnpextent = [(maxlnp-14.5, maxlnp+0.5)]
    bar = SixPosition(chain[0]) # temporary variable to get names
    foo = np.concatenate((chain, lnprob.reshape((nx, 1))), axis=1)
    if truth is None:
        truths = None
    else:
        truths = np.append(truth, [truth_lnprob])
    labels = np.append(bar.get_sixpos_names(), [r"$\ln p$"])
    extents = bar.get_sixpos_extents() + lnpextent
    fig = tri.corner(foo, labels=labels, extents=extents,
                     truths=truths, plot_contours=False)
    fn = prefix + "a.png"
    print "triangle_plot_chain(): writing " + fn
    fig.savefig(fn)
    obsfoo = 1. * foo # copy
    for i in range(nx):
        obsfoo[i,:6] = SixPosition(foo[i,:6]).get_observables_array()
    if truth is None:
        trueobs = None
    else:
        trueobs = 1. * truths # copy
        trueobs[:6] = SixPosition(truth).get_observables_array()
    labels = np.append(bar.get_observables_names(), [r"$\ln p$"])
    extents = bar.get_observables_extents() + lnpextent
    fig = tri.corner(obsfoo, labels=labels, extents=extents,
                     truths=trueobs, plot_contours=False)
    fn = prefix + "b.png"
    print "triangle_plot_chain(): writing " + fn
    fig.savefig(fn)
    intfoo = 1. * foo[:,2:] # copy
    for i in range(nx):
        intfoo[i,:4] = SixPosition(foo[i,:6]).get_integrals_of_motion()
    if truth is None:
        trueint = None
    else:
        trueint = 1. * truths[2:] # copy
        trueint[:4] = SixPosition(truth).get_integrals_of_motion()
    labels = np.append(bar.get_integrals_of_motion_names(), [r"$\ln p$"])
    extents = bar.get_integrals_of_motion_extents() + lnpextent
    fig = tri.corner(intfoo, labels=labels, extents=extents,
                     truths=trueint, plot_contours=False)
    fn = prefix + "c.png"
    print "triangle_plot_chain(): writing " + fn
    fig.savefig(fn)
    return None

    def visModel(self, labels=None, fname=None):

        nSamples = 1e6
        samples = self.reg.sample(n_samples=nSamples)
        #quick fix for better plotting
        samples = np.log10(samples)

        if self.pred == None:
            if hasTriangle:
                figure = triangle.corner(samples, labels=labels,
                                         quantiles=[0.16, 0.5, 0.84],
                                         show_titles=True, title_args={"fontsize": 12})
                if fname!=None:
                    plt.savefig('predicted_'+fname)

            elif samples.size[1]==2:
                f, ax = plt.subplots(1)
                ax.hist2d(samples[:,0], samples[:,1])
                if fname!=None:
                    plt.savefig('predicted_'+fname)

            else:
                raise NotImplementedError("Plotting datasets w/ dim > 2 without Triangle not implemented")

        else:
            if hasTriangle:
                figure = triangle.corner(samples, labels=labels,
                                         quantiles=[0.16, 0.5, 0.84],
                                         show_titles=True, title_args={"fontsize": 12})
                if fname!=None:
                    plt.savefig('predicted_'+fname)
                
                figure = triangle.corner(self.X, labels=labels,
                                         quantiles=[0.16, 0.5, 0.84],
                                         show_titles=True, title_args={"fontsize": 12})
                if fname!=None:
                    plt.savefig('original_'+fname)

            elif samples.size[1]==2:
                f, ax = plt.subplots(2)
                ax[0].hist2d(samples[:,0], samples[:,1])
                ax[1].hist2d(self.X[:,0], self.X[:,1])
                ax[1].set_xlabel(labels[0])
                ax[1].set_ylabel(labels[1])
                plt.tight_layout()
                if fname!=None:
                    plt.savefig(fname)

            else:
                raise NotImplementedError("Plotting datasets w/ dim > 2 without Triangle not implemented")

def plot_grid_points(db='phoenix2016.db3',saveplot=False):
    # plot the grid points
    connection = sqlite.connect(db)
    sql_query = 'SELECT teff,logg,mh,alpha FROM parameter_sets'
    tab = pd.read_sql_query(sql_query,connection)

    connection.close()

    arr = np.array([np.array(tab['teff']),np.array(tab['logg']),np.array(tab['mh']),tab['alpha']])
    corner(tab,plot_contours=False,labels=['Teff','log g', '[M/H]', '[alpha/Fe]'])
    if saveplot:
        fname = os.path.basename(db)
        fname = os.path.splitext(fname)
        plt.savefig(fname[0]+'.pdf')

    def triangle_plot(self,*args,**kwargs):
        
        if not args:
            args=self.keys
            
        assert len(args)>1
        
        labels=[]
        idx=[]
        for key in args:
            if key.startswith('_sigma_'):
                name=key.split('_sigma_')[1]
                label=r'$\sigma_{%s}$' % name
            else:
                namestr=key
                for g in greek:
                    if key.startswith(g):
                        namestr=r'\%s' % key

                label='$%s$' % namestr

            labels.append(label)
            idx.append(self.index[key])
        
        fig = triangle.corner(self.samples[:,idx], labels=labels)

def corner_plot(results, showpars=None, start=0, thin=1):
    #just wrap subtriangle
    """
    Make a triangle plot of the (thinned, latter) samples of the posterior
    parameter space.  Optionally make the plot only for a supplied subset
    of the parameters.
    """

    # pull out the parameter names and flatten the thinned chains
    parnames = np.array(results['model'].theta_labels())
    flatchain = results['chain'][:,start::thin,:]
    flatchain = flatchain.reshape(flatchain.shape[0] * flatchain.shape[1],
                                  flatchain.shape[2])
    truths = results['initial_center']

    # restrict to parameters you want to show
    if showpars is not None:
        ind_show = np.array([p in showpars for p in parnames], dtype= bool)
        flatchain = flatchain[:,ind_show]
        truths = truths[ind_show]
        parnames= parnames[ind_show]
        parlabels = [pardict[p] for p in parnames]
    fig = triangle.corner(flatchain, labels = parlabels,
                          quantiles=[0.16, 0.5, 0.84], verbose=False,
                          truths = truths)

    return fig

def sample(ndim, nwalkers, nsteps, start, w, ur, sigma_ur, nuvu, sigma_nuvu, age, pd, ps):
    if len(age) != len(ur):
        raise SystemExit("Number of ages does not coincide with number of galaxies...")
    p0 = [start + 1e-4 * N.random.randn(ndim) for i in range(nwalkers)]
    sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(w, ur, sigma_ur, nuvu, sigma_nuvu, age, pd, ps))
    # burn in
    pos, prob, state = sampler.run_mcmc(p0, 50)
    # reset and run with last positions of burn in run
    sampler.reset()
    print "RESET", pos
    sampler.run_mcmc(pos, nsteps)
    samples = sampler.chain[:, :, :].reshape((-1, ndim))
    samples_save = (
        "/Users/becky/Projects/Green-Valley-Project/bayesian/find_t_tau/gv/not_clean/samples_gv_not_clean_"
        + str(len(samples))
        + "_"
        + str(len(age))
        + "_"
        + str(time.strftime("%H_%M_%d_%m_%y"))
        + ".npy"
    )
    N.save(samples_save, samples)
    fig = triangle.corner(samples, labels=[r"$ t_{smooth} $", r"$ \tau_{smooth} $", r"$ t_{disc} $", r"$ \tau_{disc}$"])
    fig.savefig(
        "triangle_t_tau_gv_not_clean_"
        + str(len(samples))
        + "_"
        + str(len(age))
        + "_"
        + str(time.strftime("%H_%M_%d_%m_%y"))
        + ".pdf"
    )
    return samples, fig

def MCMC(theta, x, y, yerr, M1, ecc, fname, n, sub, plot=False):

    nwalkers, ndim = 32, len(theta)
    p0 = [theta+1e-4*np.random.rand(ndim) for i in range(nwalkers)]
    sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob,
                                    args=(x, y, yerr, M1, ecc))
    p0, lp, state = sampler.run_mcmc(p0, 500)
    sampler.reset()
    p0, lp, state = sampler.run_mcmc(p0, 5000)

    if plot == True:
        fig_labels = ['P', 'M2', 'T0', 'V0', 'omega']
        flatchain = sampler.chain[:, 50:, :].reshape((-1, ndim))
        fig = triangle.corner(flatchain, truths=theta, labels=fig_labels)
        plt.savefig("%s_triangle" % fname)

    print "saving samples"
    f = h5py.File("%s/results/%s_%s_%s_samples" % (DIR, n, fname, sub), "w")
    data = f.create_dataset("samples", np.shape(sampler.chain))
    data[:, :] = np.array(sampler.chain)
    f.close()

    flat = sampler.chain[:, 50:, :].reshape((-1, ndim))
    mcmc_result = map(lambda v: (v[1], v[2]-v[1], v[1]-v[0]),
                      zip(*np.percentile(flat, [16, 50, 84], axis=0)))
    np.savetxt("%s/results/%s_%s_%s_results.txt" % (DIR, n, fname, sub),
               mcmc_result)

    return mcmc_result

def triangle_gen():
    plt.clf()
    plot_chi()
    burnIn = input('Burn off how many steps? ')
    samples = chainDat[:,burnIn:,:].reshape((-1,ndim))
    triangle.corner(samples, labels=['$R_{In}$ [AU]','$\Delta$R [AU]','log($M_{D}$ [$M_{Earth}$])','i [degrees]','PA [degrees]'], truths=[best_disk[0],best_disk[1],best_disk[2],best_disk[3],best_disk[4]], truth_color='g', quantiles=[0.16,.5,.84], show_titles=True, plot_contours=True)
    plt.savefig('MCMCRUNS/'+append+'/'+whatbywhat+'/'+disk_name+'_'+whatbywhat+'_'+append+'_Triangle.png')
    plt.show()
    
    values = np.zeros((ndim,3))
    for i in range(ndim):
            chainSqueeze=chainDat[:,burnIn:,i].reshape(-1)
        quantiles = mstats.mquantiles(chainSqueeze,prob=[.16,.5,.84],axis=None)
        values[i,0] = quantiles[1] #median
        values[i,1] = quantiles[2]-quantiles[1] #plus one sigma
        values[i,2] = quantiles[0]-quantiles[1] #minus one sigma