Skip to content
Snippets Groups Projects
Commit 053bf2d0 authored by Rayne Liu's avatar Rayne Liu
Browse files

Finally found the bug with the pool but did not understand it - but PyCBC pool...

Finally found the bug with the pool but did not understand it - but PyCBC pool also worked with ptemcee
parent 0a01c216
Branches
No related tags found
No related merge requests found
......@@ -16,6 +16,8 @@ from matplotlib import rc
plt.rcParams.update({'font.size': 16.5})
import ptemcee
from pycbc.pool import choose_pool
#from multiprocessing import Pool
import math
import h5py
import inspect
......@@ -32,18 +34,19 @@ from scipy.optimize import minimize
#tshift: time shift after the strain peak
#vary_fund: whether you vary the fundamental frequency. Works in the model_dv function.
rootpath="/work/rayne.liu" #"/Users/RayneLiu"
rootpath="/Users/RayneLiu"# "/work/rayne.liu"
nmax=1
tshift=19
vary_fund = True
#sampler parameters
npoints=101
nwalkers = 32
nwalkers = 42
ntemps=12
ndim = int(4*(nmax+1))
burnin = 10 #How many points do you burn before doing the corner plot. You need to watch the convergence of the chain plot a bit.
#This is trivial but often forgotten: this cannot be more than npoints! Usually 1/5~1/4 npoints is what I observe.
numbins = 21 #corner plot parameter - how many bins you want
datacolor = '#105670' #'#4fa3a7'
pkcolor = '#f2c977' #'#ffb45f'
......@@ -160,8 +163,8 @@ def log_likelihood(theta):
# The evidence is just a normalization factor
def log_probability(theta):
lp = log_prior(theta)
print('lp:')
print(lp)
#print('lp:')
#print(lp)
if not np.isfinite(lp):
return -np.inf
return lp + log_likelihood(theta)
......@@ -170,6 +173,9 @@ def log_probability(theta):
#Fit with ptemcee
#Set the number of cores of your processors
pool = choose_pool(4)
pool.size = 4
vary_param = float(vary_fund)
pos = np.array([[random.uniform(-0.1,0.1), random.uniform(-0.1,0.1), 4.28313743e-01, random.uniform(2.5, 2.6) + (1-vary_param) * np.pi]])
for i in range (1,nmax+1):
......@@ -180,7 +186,8 @@ pos = pos.T.flatten()
pos = list(pos)
pos += 1e-5 * np.random.randn(ntemps, nwalkers, ndim)
#print(pos)
sampler = ptemcee.Sampler(nwalkers, ndim, log_likelihood, log_prior, ntemps=ntemps)
sampler = ptemcee.Sampler(nwalkers, ndim, log_likelihood, log_prior, ntemps=ntemps, pool=pool)
sampler.run_mcmc(pos,npoints)
......@@ -291,13 +298,3 @@ plt.xlabel("t")
plt.ylabel("h")
figband.savefig(rootpath+'/git/rdstackingproject/plotsmc/vary'+str(vary_fund)+'nmax='+str(nmax)+'_tshift='+str(tshift)+'_'+str(npoints)+'pt_band.pdf', format = 'pdf')
"""
"""
......@@ -16,6 +16,7 @@ from matplotlib import rc
plt.rcParams.update({'font.size': 16.5})
import ptemcee
from pycbc.pool import choose_pool
import math
import h5py
import inspect
......@@ -170,6 +171,9 @@ def log_probability(theta):
#Fit with ptemcee
#Set the number of cores of your processors
pool = choose_pool(4)
pool.size = 4
vary_param = float(vary_fund)
pos = np.array([[random.uniform(-0.1,0.1), random.uniform(-0.1,0.1), 4.28313743e-01, random.uniform(2.5, 2.6) + (1-vary_param) * np.pi]])
for i in range (1,nmax+1):
......@@ -180,12 +184,11 @@ pos = pos.T.flatten()
pos = list(pos)
pos += 1e-5 * np.random.randn(ntemps, nwalkers, ndim)
#print(pos)
sampler = ptemcee.Sampler(nwalkers, ndim, log_likelihood, log_prior, ntemps=ntemps)
sampler = ptemcee.Sampler(nwalkers, ndim, log_likelihood, log_prior, ntemps=ntemps, pool=pool)
sampler.run_mcmc(pos,npoints)
#Define labels and start plotting
paramlabels_a = [r'$\alpha_'+str(i)+'$' for i in range (nmax+1)]
paramlabels_b = [r'$\beta_'+str(i)+'$' for i in range (nmax+1)]
......
......@@ -16,6 +16,7 @@ from matplotlib import rc
plt.rcParams.update({'font.size': 16.5})
import ptemcee
from pycbc.pool import choose_pool
import math
import h5py
import inspect
......@@ -170,6 +171,9 @@ def log_probability(theta):
#Fit with ptemcee
#Set the number of cores of your processors
pool = choose_pool(4)
pool.size = 4
vary_param = float(vary_fund)
pos = np.array([[random.uniform(-0.1,0.1), random.uniform(-0.1,0.1), 4.28313743e-01, random.uniform(2.5, 2.6) + (1-vary_param) * np.pi]])
for i in range (1,nmax+1):
......@@ -180,7 +184,7 @@ pos = pos.T.flatten()
pos = list(pos)
pos += 1e-5 * np.random.randn(ntemps, nwalkers, ndim)
#print(pos)
sampler = ptemcee.Sampler(nwalkers, ndim, log_likelihood, log_prior, ntemps=ntemps)
sampler = ptemcee.Sampler(nwalkers, ndim, log_likelihood, log_prior, ntemps=ntemps, pool=pool)
sampler.run_mcmc(pos,npoints)
......
......@@ -16,6 +16,7 @@ from matplotlib import rc
plt.rcParams.update({'font.size': 16.5})
import ptemcee
from pycbc.pool import choose_pool
import math
import h5py
import inspect
......@@ -170,6 +171,9 @@ def log_probability(theta):
#Fit with ptemcee
#Set the number of cores of your processors
pool = choose_pool(4)
pool.size = 4
vary_param = float(vary_fund)
pos = np.array([[random.uniform(-0.1,0.1), random.uniform(-0.1,0.1), 4.28313743e-01, random.uniform(2.5, 2.6) + (1-vary_param) * np.pi]])
for i in range (1,nmax+1):
......@@ -180,7 +184,7 @@ pos = pos.T.flatten()
pos = list(pos)
pos += 1e-5 * np.random.randn(ntemps, nwalkers, ndim)
#print(pos)
sampler = ptemcee.Sampler(nwalkers, ndim, log_likelihood, log_prior, ntemps=ntemps)
sampler = ptemcee.Sampler(nwalkers, ndim, log_likelihood, log_prior, ntemps=ntemps, pool=pool)
sampler.run_mcmc(pos,npoints)
......
......@@ -16,6 +16,7 @@ from matplotlib import rc
plt.rcParams.update({'font.size': 16.5})
import ptemcee
from pycbc.pool import choose_pool
import math
import h5py
import inspect
......@@ -170,6 +171,9 @@ def log_probability(theta):
#Fit with ptemcee
#Set the number of cores of your processors
pool = choose_pool(4)
pool.size = 4
vary_param = float(vary_fund)
pos = np.array([[random.uniform(-0.1,0.1), random.uniform(-0.1,0.1), 4.28313743e-01, random.uniform(2.5, 2.6) + (1-vary_param) * np.pi]])
for i in range (1,nmax+1):
......@@ -180,7 +184,7 @@ pos = pos.T.flatten()
pos = list(pos)
pos += 1e-5 * np.random.randn(ntemps, nwalkers, ndim)
#print(pos)
sampler = ptemcee.Sampler(nwalkers, ndim, log_likelihood, log_prior, ntemps=ntemps)
sampler = ptemcee.Sampler(nwalkers, ndim, log_likelihood, log_prior, ntemps=ntemps, pool=pool)
sampler.run_mcmc(pos,npoints)
......
......@@ -12,8 +12,8 @@ initialdir = .
notify_user = rl746@cornell.edu
notification = Complete
arguments = "-processid $(Process)"
request_memory = 16GB
request_cpus = 1
request_memory = 8GB
request_cpus = 4
on_exit_remove = (ExitBySignal == False) || ((ExitBySignal == True) && (ExitSignal != 11))
accounting_group = aei.dev.test_dynesty
queue 1
......
......@@ -12,8 +12,8 @@ initialdir = .
notify_user = rl746@cornell.edu
notification = Complete
arguments = "-processid $(Process)"
request_memory = 16GB
request_cpus = 1
request_memory = 8GB
request_cpus = 4
on_exit_remove = (ExitBySignal == False) || ((ExitBySignal == True) && (ExitSignal != 11))
accounting_group = aei.dev.test_dynesty
queue 1
......
......@@ -12,8 +12,8 @@ initialdir = .
notify_user = rl746@cornell.edu
notification = Complete
arguments = "-processid $(Process)"
request_memory = 16GB
request_cpus = 1
request_memory = 8GB
request_cpus = 4
on_exit_remove = (ExitBySignal == False) || ((ExitBySignal == True) && (ExitSignal != 11))
accounting_group = aei.dev.test_dynesty
queue 1
......
......@@ -12,8 +12,8 @@ initialdir = .
notify_user = rl746@cornell.edu
notification = Complete
arguments = "-processid $(Process)"
request_memory = 16GB
request_cpus = 1
request_memory = 8GB
request_cpus = 4
on_exit_remove = (ExitBySignal == False) || ((ExitBySignal == True) && (ExitSignal != 11))
accounting_group = aei.dev.test_dynesty
queue 1
......
......@@ -12,8 +12,8 @@ initialdir = .
notify_user = rl746@cornell.edu
notification = Complete
arguments = "-processid $(Process)"
request_memory = 16GB
request_cpus = 1
request_memory = 8GB
request_cpus = 4
on_exit_remove = (ExitBySignal == False) || ((ExitBySignal == True) && (ExitSignal != 11))
accounting_group = aei.dev.test_dynesty
queue 1
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment