EMCEE sampling script

Emcee Snippet

# Initialize the MPI-based pool used for parallelization.
pool = emcee.utils.MPIPool()

if not pool.is_master():
    # Wait for instructions from the master process.
    pool.wait()
    sys.exit(0)

nwalkers = 1500
ndim = len(params)
sampler = emcee.EnsembleSampler(nwalkers, ndim, cost,  pool=pool)
t0 = time.time()
# Run 100 steps as a burn-in.

eps = 1e-3
p0 = (np.random.randn(len(params), nwalkers) * eps + params[:,np.newaxis]).T

p0[p0 >= 12] = 11.99
p0[p0 <= -12] = -11.99
# Avoid infinite log prior at the start


# If run on cluster:

if not "keu" in cwd:
    base_output_dir = os.path.join(cwd, '..', '..', 'Output', 'MCMC')
else:
    base_output_dir = "/panfs/panfem/fem1/keu_vaggi/p53/Output/MCMC"

output_dir = os.path.join(base_output_dir, 'All_Experiments', model_name, settings_name)

if not os.path.exists(output_dir):
    os.makedirs(output_dir)
f = open(os.path.join(output_dir, "chain.dat"), "w")
f.close()

t0 = time.time()
# Sample burn-in
for result in sampler.sample(p0, iterations=3000, storechain=False):
    pass

position = result[0]

t1 = time.time()

print "Partial time: %.4f" % (t1-t0)

iter_n = 0
for result in sampler.sample(position, iterations=1000, storechain=False):
    iter_n += 1
    if np.mod(iter_n, 10) == 0:
        position = result[0]
        f = open(os.path.join(output_dir, "chain.dat"), "a")
        np.savetxt(f, position, fmt = '%.8f', delimiter = ',')
        f.close()

pool.close()
t2 = time.time()

print "Total time: %.4f" % (t2-t0)