Python ChainConsumer.plot方法代码示例

# 需要导入模块: from chainconsumer import ChainConsumer [as 别名]
# 或者: from chainconsumer.ChainConsumer import plot [as 别名]
def plot_results(chain, param, chainf, chainf2, paramf, t0, x0, x1, c, temp_dir, seed, interped):
    cc = ChainConsumer()
    cc.add_chain(chain, parameters=param, name="Posterior")
    cc.add_chain(chainf, parameters=paramf, name="Minuit")
    cc.add_chain(chainf2, parameters=paramf, name="Emcee")
    truth = {"$t_0$": t0, "$x_0$": x0, "$x_1$": x1, "$c$": c, r"$\mu$": get_mu(interped, x0, x1, c)}
    cc.plot(filename=temp_dir + "/surfaces_%d.png" % seed, truth=truth)

# 需要导入模块: from chainconsumer import ChainConsumer [as 别名]
# 或者: from chainconsumer.ChainConsumer import plot [as 别名]
def debug_plots(std):
    print(std)

    res = load_stan_from_folder(std, merge=True, cut=False)
    chain, posterior, t, p, f, l, w, ow = res
    # print(w.mean())
    # import matplotlib.pyplot as plt
    # plt.hist(np.log(w), 100)
    # plt.show()
    # exit()
    logw = np.log(w)
    m = np.mean(logw)
    s = np.std(logw)
    print(m, s)
    logw -= (m + 3 * s)
    good = logw < 0
    logw *= good
    w = np.exp(logw)

    c = ChainConsumer()
    c.add_chain(chain, weights=w, name="corrected")
    c.configure(summary=True)
    c.plot(figsize=2.0, filename="output.png", parameters=9)

    c = ChainConsumer()
    c.add_chain(chain, name="uncorrected")
    c.add_chain(chain, weights=w, name="corrected")
    # c.add_chain(chain, name="calib")
    c.plot(filename="output_comparison.png", parameters=9, figsize=1.3)
    c.plot_walks(chains=1, filename="walks.png")

# 需要导入模块: from chainconsumer import ChainConsumer [as 别名]
# 或者: from chainconsumer.ChainConsumer import plot [as 别名]
def plot_all_no_weight(folder, output):
    """ Plot all chains as one, with and without weights applied """
    print("Plotting all as one, with old and new weights")
    chain, posterior, t, p, f, l, w, ow = load_stan_from_folder(folder, merge=True)
    c = ChainConsumer()
    c.add_chain(chain, posterior=posterior, walkers=l)
    c.plot(filename=output, truth=t, figsize=0.75)

# 需要导入模块: from chainconsumer import ChainConsumer [as 别名]
# 或者: from chainconsumer.ChainConsumer import plot [as 别名]
def plot_separate(folder, output):
    """ Plot separate cosmologies """
    print("Plotting all cosmologies separately")
    res = load_stan_from_folder(folder, merge=False)
    c = ChainConsumer()
    for i, (chain, posterior, t, p, f, l, w, ow) in enumerate(res):
        c.add_chain(chain, weights=w, posterior=posterior, walkers=l, name="%d"%i)
    c.plot(filename=output, truth=t, figsize=0.75)

# 需要导入模块: from chainconsumer import ChainConsumer [as 别名]
# 或者: from chainconsumer.ChainConsumer import plot [as 别名]
def plot_single_cosmology_weight(folder, output, i=0):
    print("Plotting cosmology realisation %d" % i)
    res = load_stan_from_folder(folder, merge=False)
    c = ChainConsumer()
    chain, posterior, t, p, f, l, w, ow = res[i]
    c.add_chain(chain, posterior=posterior, walkers=l, name="Uncorrected %d"%i)
    c.add_chain(chain, weights=w, posterior=posterior, walkers=l, name="Corrected %d"%i)
    c.plot(filename=output, truth=t, figsize=0.75)

# 需要导入模块: from chainconsumer import ChainConsumer [as 别名]
# 或者: from chainconsumer.ChainConsumer import plot [as 别名]
def plot_all(folder, output, output_walk=None):
    """ Plot all chains as one """
    print("Plotting all as one")
    chain, posterior, t, p, f, l, w, ow = load_stan_from_folder(folder, merge=True)
    c = ChainConsumer()
    c.add_chain(chain, weights=w, posterior=posterior, walkers=l)
    c.plot(filename=output, truth=t, figsize=0.75)
    if output_walk is not None:
        c.plot_walks(filename=output_walk)

# 需要导入模块: from chainconsumer import ChainConsumer [as 别名]
# 或者: from chainconsumer.ChainConsumer import plot [as 别名]
def plot_separate_weight(folder, output):
    """ Plot separate cosmologies, with and without weights applied """
    print("Plotting all cosmologies separately, with old and new weights")
    res = load_stan_from_folder(folder, merge=False)
    c = ChainConsumer()
    ls = []
    for i, (chain, posterior, t, p, f, l, w, ow) in enumerate(res):
        c.add_chain(chain, posterior=posterior, walkers=l, name="Uncorrected %d"%i)
        c.add_chain(chain, weights=w, posterior=posterior, walkers=l, name="Corrected %d"%i)
        ls += ["-", "--"]
    c.configure_general(linestyles=ls)
    c.plot(filename=output, truth=t, figsize=0.75)

# 需要导入模块: from chainconsumer import ChainConsumer [as 别名]
# 或者: from chainconsumer.ChainConsumer import plot [as 别名]
    
    daval = (alpha/(1+epsilon)) * da #/ rs_fid
    
    hrc = hs / (alpha * (1 + epsilon) * (1 + epsilon)) #* rs_fid / (alpha * (1 + epsilon) * (1 + epsilon)) / c
    res = np.vstack((omch2, daval, hrc)).T
    return res
    
    
    
if False:
    consumer = ChainConsumer()
    consumer.configure_contour(sigmas=[0,1,2])
    consumer.add_chain(load_directory("../bWigMpBin/bWigMpBin_z0"), parameters=[r"$\Omega_c h^2$", r"$\alpha$", r"$\epsilon$"], name="$0.2<z<0.6$")
    consumer.add_chain(load_directory("../bWigMpBin/bWigMpBin_z1"), parameters=[r"$\Omega_c h^2$", r"$\alpha$", r"$\epsilon$"], name="$0.4<z<0.8$")
    consumer.add_chain(load_directory("../bWigMpBin/bWigMpBin_z2"), parameters=[r"$\Omega_c h^2$", r"$\alpha$", r"$\epsilon$"], name="$0.6<z<1.0$")
    consumer.plot(figsize="column", filename="wigglez_multipole_alphaepsilon.pdf")
    #print(consumer.get_latex_table())

if True:
    c = ChainConsumer()
    c.add_chain(convert_directory("../bWigMpBin/bWigMpBin_z0", 0.44), parameters=[r"$\Omega_c h^2$", r"$D_A(z)$", r"$H(z)$"], name="$0.2<z<0.6$")
    c.add_chain(convert_directory("../bWigMpBin/bWigMpBin_z1", 0.60), parameters=[r"$\Omega_c h^2$", r"$D_A(z)$", r"$H(z)$"], name="$0.4<z<0.8$")
    c.add_chain(convert_directory("../bWigMpBin/bWigMpBin_z2", 0.73), parameters=[r"$\Omega_c h^2$", r"$D_A(z)$", r"$H(z)$"], name="$0.6<z<1.0$")
    for n, chain in zip(c.names, c.chains):
        print(n)
        print(chain.mean(axis=0))
        print(np.std(chain, axis=0))
        print("----")
    c.configure_contour(sigmas=[0,1,2])
    c.configure_general(bins=0.7)    
    c.plot(figsize="column", filename="wigglez_multipole_dah.pdf")

# 需要导入模块: from chainconsumer import ChainConsumer [as 别名]
# 或者: from chainconsumer.ChainConsumer import plot [as 别名]
        model_good = EfficiencyModelUncorrected(cs, zs, ts, calibration, zeros, ls, ss, t0s, name="Good%d" % i)
        model_good.fit(sampler, chain_consumer=c)

        model_un = EfficiencyModelUncorrected(cs[mask], zs[mask], ts[mask], calibration,
                                              zeros, ls[mask], ss[mask], t0s[mask], name="Uncorrected%d" % i)
        model_un.fit(sampler, chain_consumer=c)

        biased_chain = c.chains[-1]
        # model_cor.fit(sampler, chain_consumer=c)

        filename = dir_name + "/output/weights.txt"
        if not os.path.exists(filename):
            weights = []
            for i, row in enumerate(biased_chain):
                weights.append(get_weights(row[0], row[1], row[2], row[3], row[4], row[5], threshold))
                print(100.0 * i / biased_chain.shape[0])
            weights = np.array(weights)
            np.savetxt(filename, weights)
        else:
            weights = np.loadtxt(filename)
        weights = (1 / np.power(weights, mask.sum()))
        c.add_chain(biased_chain, name="Importance Sampled", weights=weights)

    c.configure_bar(shade=True)
    c.configure_general(bins=1.0, colours=colours)
    c.configure_contour(sigmas=[0, 0.01, 1, 2], contourf=True, contourf_alpha=0.2)
    c.plot(filename=plot_file, truth=theta, figsize=(7, 7), legend=False, parameters=6)
    for i in range(len(c.chains)):
        c.plot_walks(filename=walk_file % c.names[i], chain=i, truth=theta)

# 需要导入模块: from chainconsumer import ChainConsumer [as 别名]
# 或者: from chainconsumer.ChainConsumer import plot [as 别名]
    t = os.path.abspath(dir_name + "/output/data_%d")
    plot_file = os.path.abspath(dir_name + "/output/surfaces.png")
    walk_file = os.path.abspath(dir_name + "/output/walk_%s.png")

    c = ChainConsumer()
    n = 2
    colours = ["#4CAF50", "#D32F2F", "#1E88E5"] * n  # , "#FFA000"] * n
    for i in range(n):
        mean, sigma, cut, observed, mask = get_data(seed=i)

        model_good = EfficiencyModelUncorrected(observed, name="Good")
        model_un = EfficiencyModelUncorrected(observed[mask])
        model_cor = EfficiencyModelCorrected(observed[mask], cut)

        sampler = EnsembleSampler(num_steps=25000, num_burn=1000, temp_dir=t % i)
        model_good.fit(sampler, chain_consumer=c)
        model_un.fit(sampler, chain_consumer=c)
        biased_chain = c.chains[-1]
        # model_cor.fit(sampler, chain_consumer=c)

        mus = biased_chain[:, 0]
        sigmas = biased_chain[:, 1]
        weights = 1 / get_weights(cut, mus, sigmas, mask.sum())

        c.add_chain(biased_chain, name="Importance Sampled", weights=weights)

    c.configure_bar(shade=True)
    c.configure_general(colours=colours, bins=0.5)
    c.configure_contour(contourf=True, contourf_alpha=0.2)
    c.plot(filename=plot_file, figsize=(5, 5), truth=[mean, sigma], legend=False)

# 需要导入模块: from chainconsumer import ChainConsumer [as 别名]
# 或者: from chainconsumer.ChainConsumer import plot [as 别名]
    print(mask.sum(), n, observed.mean())
    return mean, observed, errors, alpha


if __name__ == "__main__":
    logging.basicConfig(level=logging.DEBUG)
    dir_name = os.path.dirname(__file__)
    t = os.path.abspath(dir_name + "/output/run_%d")
    plot_file = os.path.abspath(dir_name + "/output/surfaces.png")
    walk_file = os.path.abspath(dir_name + "/output/walk_%s.png")

    c = ChainConsumer()
    n = 3
    colours = ["#D32F2F", "#1E88E5"] * n
    for i in range(n):
        mean, observed, errors, alpha = get_data(seed=i)

        model_un = EfficiencyModelUncorrected(observed, errors)
        model_cor = EfficiencyModelCorrected(observed, errors, alpha)

        pgm_file = os.path.abspath(dir_name + "/output/pgm.png")
        fig = model_cor.get_pgm(pgm_file)

        sampler = EnsembleSampler(num_steps=5000, num_burn=1000, temp_dir=t % i, num_walkers=50)
        model_un.fit(sampler, chain_consumer=c)
        model_cor.fit(sampler, chain_consumer=c)

    c.configure_bar(shade=True)
    c.configure_general(colours=colours)
    c.plot(filename=plot_file, figsize=(5, 3), truth=[mean], legend=False)

# 需要导入模块: from chainconsumer import ChainConsumer [as 别名]
# 或者: from chainconsumer.ChainConsumer import plot [as 别名]
    for i in range(n):
        mean, std, observed, errors, alpha, actual, uo, oe, am = get_data(seed=i)
        theta_good = [mean, std] + actual.tolist()
        theta_bias = [mean, std] + am.tolist()
        kwargs = {"num_steps": 70000, "num_burn": 20000, "save_interval": 300,
                  "plot_covariance": True, "unify_latent": True}  # , "callback": v.callback
        sampler = BatchMetropolisHastings(num_walkers=w, kwargs=kwargs, temp_dir=t % i, num_cores=4)

        model_good = EfficiencyModelUncorrected(uo, oe, name="Good%d" % i)
        model_good.fit(sampler, chain_consumer=c)
        print("Good ", model_good.get_log_posterior(theta_good), c.posteriors[-1][-1])

        model_un = EfficiencyModelUncorrected(observed, errors, name="Uncorrected%d" % i)
        model_un.fit(sampler, chain_consumer=c)
        print("Uncorrected ", model_un.get_log_posterior(theta_bias), c.posteriors[-1][-1])

        model_cor = EfficiencyModelCorrected(observed, errors, alpha, name="Corrected%d" % i)
        model_cor.fit(sampler, chain_consumer=c)
        print("Corrected ", model_cor.get_log_posterior(theta_bias), c.posteriors[-1][-1])

    c.configure_bar(shade=True)
    c.configure_general(bins=1.0, colours=colours)
    c.configure_contour(sigmas=[0, 0.01, 1, 2], shade=True, shade_alpha=0.3)
    c.plot(filename=plot_file, truth=theta_bias, figsize=(5, 5), legend=False)
    for i in range(len(c.chains)):
        c.plot_walks(filename=walk_file % c.names[i], chain=i, truth=[mean, std])
        # c.divide_chain(i, w).configure_general(rainbow=True) \
        #     .plot(figsize=(5, 5), filename=plot_file.replace(".png", "_%s.png" % c.names[i]),
        #           truth=theta_bias)

# 需要导入模块: from chainconsumer import ChainConsumer [as 别名]
# 或者: from chainconsumer.ChainConsumer import plot [as 别名]
def random_obs(temp_dir, seed):
    np.random.seed(seed)
    interp = generate_and_return()
    x1 = np.random.normal()
    # colour = np.random.normal(scale=0.1)
    colour = 0
    x0 = 1e-5
    # t0 = np.random.uniform(low=1000, high=2000)
    t0 = 1000
    z = np.random.uniform(low=0.1, high=1.0)

    # deltat = np.random.uniform(low=-20, high=0)
    # num_obs = np.random.randint(low=10, high=40)
    num_obs = 20
    deltat = -35

    filename = temp_dir + "/save_%d.npy" % seed

    if not os.path.exists(filename):
        ts = np.arange(t0 + deltat, (t0 + deltat) + 5 * num_obs, 5)

        times = np.array([[t, t + 0.05, t + 0.1, t + 0.2] for t in ts]).flatten()
        bands = [b for t in ts for b in ["desg", "desr", "desi", "desz"]]
        gains = np.ones(times.shape)
        skynoise = np.random.uniform(low=20, high=800) * np.ones(times.shape)
        zp = 30 * np.ones(times.shape)
        zpsys = ["ab"] * times.size

        obs = Table({"time": times, "band": bands, "gain": gains, "skynoise": skynoise, "zp": zp, "zpsys": zpsys})
        model = sncosmo.Model(source="salt2")
        p = {"z": z, "t0": t0, "x0": x0, "x1": x1, "c": colour}
        model.set(z=z)
        print(seed, " Vals are ", p)
        lc = sncosmo.realize_lcs(obs, model, [p])[0]
        ston = (lc["flux"] / lc["fluxerr"]).max()

        model.set(t0=t0, x1=x1, c=colour, x0=x0)
        try:
            res, fitted_model = sncosmo.fit_lc(
                lc, model, ["t0", "x0", "x1", "c"], guess_amplitude=False, guess_t0=False
            )
        except ValueError:
            return np.nan, np.nan, x1, colour, num_obs, ston, deltat, z, 0

        fig = sncosmo.plot_lc(lc, model=fitted_model, errors=res.errors)
        fig.savefig(temp_dir + os.sep + "lc_%d.png" % seed, bbox_inches="tight", dpi=300)
        my_model = PerfectRedshift([lc], [z], t0, name="posterior%d" % seed)
        sampler = EnsembleSampler(temp_dir=temp_dir, num_burn=400, num_steps=1500)
        c = ChainConsumer()
        my_model.fit(sampler, chain_consumer=c)
        map = {"x0": "$x_0$", "x1": "$x_1$", "c": "$c$", "t0": "$t_0$"}
        parameters = [map[a] for a in res.vparam_names]

        mu1 = get_mu_from_chain(interped, c.chains[-1], c.parameters[-1])
        c.parameteers[-1].append(r"$\mu$")
        c.chains[-1] = np.hstack((c.chains[-1], mu1[:, None]))

        chain2 = np.random.multivariate_normal(res.parameters[1:], res.covariance, size=int(1e5))
        chain2 = np.hstack((chain2, get_mu_from_chain(interp, chain2, parameters)[:, None]))
        c.add_chain(chain2, parameters=parameters, name="Gaussian")
        figfilename = filename.replace(".npy", ".png")
        c.plot(filename=figfilename, truth={"$t_0$": t0, "$x_0$": x0, "$x_1$": x1, "$c$": colour})

        means = []
        stds = []
        isgood = (
            (np.abs(x1 - res.parameters[3]) < 4) & (np.abs(colour - res.parameters[4]) < 2) & (res.parameters[2] > 0.0)
        )
        isgood *= 1.0

        for i in range(len(c.chains)):
            a = c.chains[i][:, -1]
            means.append(a.mean())
            stds.append(np.std(a))
        diffmu = np.diff(means)[0]
        diffstd = np.diff(stds)[0]
        np.save(filename, np.array([diffmu, diffstd, ston, 1.0 * isgood]))

    else:
        vals = np.load(filename)
        diffmu = vals[0]
        diffstd = vals[1]
        ston = vals[2]
        isgood = vals[3]

    return diffmu, diffstd, x1, colour, num_obs, ston, deltat, z, isgood

# 需要导入模块: from chainconsumer import ChainConsumer [as 别名]
# 或者: from chainconsumer.ChainConsumer import plot [as 别名]
    c.add_chain(np.random.multivariate_normal(res.parameters[1:], res.covariance, size=int(1e7)),
                name="Summary Stats", parameters=["$t_0$", "$x_0$", "$x_1$", "$c$"])

    if False:
        if not os.path.exists(mcmc_chain):
            res2, fitted_model2 = sncosmo.mcmc_lc(lcs[0], model, ['t0', 'x0', 'x1', 'c'], nwalkers=20,
                                                  nburn=500, nsamples=4000)
            mcchain = res2.samples
            np.save(mcmc_chain, mcchain)
        else:
            mcchain = np.load(mcmc_chain)
        c.add_chain(mcchain, name="sncosmo mcmc", parameters=["$t_0$", "$x_0$", "$x_1$", "$c$"])
    print("Plot surfaces")
    c.configure_contour(shade=True, shade_alpha=0.2, sigmas=[0.0, 1.0, 2.0, 3.0])
    c.configure_bar(shade=True)
    c.plot(filename=surface, figsize=(7, 7))
    if False:
        fig = sncosmo.plot_lc(lcs[0], model=fitted_model, errors=res.errors)
        fig.savefig(temp_dir + os.sep + "lc_simple.png", bbox_inches="tight", dpi=300)

    alpha = 0.14
    beta = 3.15

    c2 = ChainConsumer()
    means = []
    stds = []
    print("Add chains")
    for i in range(len(c.chains)):
        chain = c.chains[i]
        apparent_interp = generate_and_return()
        x0s = chain[:, c.parameters[i].index("$x_0$")]

# 需要导入模块: from chainconsumer import ChainConsumer [as 别名]
# 或者: from chainconsumer.ChainConsumer import plot [as 别名]
import numpy as np
import pickle
import os
from chainconsumer import ChainConsumer

dir_name = os.path.dirname(__file__)
c = ChainConsumer()

for i in range(10):
    t = os.path.abspath(dir_name + "/output/temp%d.pkl" % i)
    with open(t, 'rb') as output:
        chain = pickle.load(output)
        chain = np.vstack((chain["mu"], chain["sigma"]))

    c.add_chain(chain.T, parameters=[r"$\mu$", r"$\sigma$"])
c.configure_bar(shade=True)
c.configure_contour(sigmas=[0, 0.01, 1, 2], contourf=True, contourf_alpha=0.1)
c.plot(display=True, truth=[100, 20])