diff options
| -rwxr-xr-x | fr.py | 30 | ||||
| -rw-r--r-- | submitter/make_dag.py | 18 | ||||
| -rw-r--r-- | submitter/submit.sub | 2 | ||||
| -rw-r--r-- | utils/gf.py | 21 | ||||
| -rw-r--r-- | utils/likelihood.py | 11 | ||||
| -rw-r--r-- | utils/mcmc.py | 16 | ||||
| -rw-r--r-- | utils/misc.py | 15 | ||||
| -rw-r--r-- | utils/plot.py | 72 |
8 files changed, 116 insertions, 69 deletions
@@ -235,35 +235,29 @@ def main(): else: fitter = None - triangle_llh = partial( - llh_utils.triangle_llh, args=args, mcmc_paramset=mcmc_paramset, - asimov_paramset=asimov_paramset, fitter=fitter - ) - lnprior = partial( - llh_utils.lnprior, paramset=mcmc_paramset + ln_prob = partial( + llh_utils.ln_prob, args=args, fitter=fitter, + asimov_paramset=asimov_paramset, mcmc_paramset=mcmc_paramset ) ndim = len(mcmc_paramset) - ntemps = 1 if args.mcmc_seed_type == MCMCSeedType.UNIFORM: p0 = mcmc_utils.flat_seed( - mcmc_paramset, ntemps=ntemps, nwalkers=args.nwalkers + mcmc_paramset, nwalkers=args.nwalkers ) elif args.mcmc_seed_type == MCMCSeedType.GAUSSIAN: p0 = mcmc_utils.gaussian_seed( - mcmc_paramset, ntemps=ntemps, nwalkers=args.nwalkers + mcmc_paramset, nwalkers=args.nwalkers ) samples = mcmc_utils.mcmc( - p0 = p0, - triangle_llh = triangle_llh, - lnprior = lnprior, - ndim = ndim, - nwalkers = args.nwalkers, - burnin = args.burnin, - nsteps = args.nsteps, - ntemps = ntemps, - threads = 1 + p0 = p0, + ln_prob = ln_prob, + ndim = ndim, + nwalkers = args.nwalkers, + burnin = args.burnin, + nsteps = args.nsteps, + threads = 1 # TODO(shivesh): broken because you cannot pickle a GolemFitPy object # threads = misc_utils.thread_factors(args.threads)[0] ) diff --git a/submitter/make_dag.py b/submitter/make_dag.py index f5a06dd..216ca72 100644 --- a/submitter/make_dag.py +++ b/submitter/make_dag.py @@ -16,32 +16,32 @@ full_scan_mfr = [ ] fix_sfr_mfr = [ - # (1, 1, 1, 1, 0, 0), - # (1, 1, 1, 0, 1, 0), + (1, 1, 1, 1, 0, 0), + (1, 1, 1, 0, 1, 0), # (1, 1, 1, 0, 0, 1), (1, 1, 1, 1, 2, 0), # (1, 1, 0, 0, 1, 0), - # (1, 1, 0, 1, 2, 0), + (1, 1, 0, 1, 2, 0), # (1, 1, 0, 1, 0, 0), # (1, 0, 0, 1, 0, 0), - # (0, 1, 0, 0, 1, 0), + (0, 1, 0, 0, 1, 0), # (1, 2, 0, 0, 1, 0), # (1, 2, 0, 1, 2, 0) ] # MCMC run_mcmc = 'True' -burnin = 50 -nsteps = 500 -nwalkers = 200 +burnin = 1000 +nsteps = 4000 +nwalkers = 70 seed = 24 threads = 12 mcmc_seed_type = 'uniform' # FR -dimension = [3] +dimension = [3, 6] energy = [1e6] -likelihood = 'gaussian' +likelihood = 'golemfit' no_bsm = 'False' sigma_ratio = ['0.01'] scale = "1E-20 1E-30" diff --git a/submitter/submit.sub b/submitter/submit.sub index f7c4e7e..b47ac5c 100644 --- a/submitter/submit.sub +++ b/submitter/submit.sub @@ -15,7 +15,7 @@ getenv = True # but do not try to copy outputs back (see: https://htcondor-wiki.cs.wisc.edu/index.cgi/tktview?tn=3081) # +TransferOutput="" -request_memory = 1GB +request_memory = 2GB request_cpus = 12 Universe = vanilla diff --git a/utils/gf.py b/utils/gf.py index 5ad2869..ebc8538 100644 --- a/utils/gf.py +++ b/utils/gf.py @@ -22,21 +22,15 @@ from utils.misc import enum_parse, thread_factors def steering_params(args): steering_categ = args.ast params = gf.SteeringParams() - if 'cobalt0' in socket.gethostname().split('.')[0]: - params.quiet = False - else: - params.quiet = True - # TODO(shivesh): figure out right number for missid - params.track_to_shower_missid = 0.3 + params.quiet = False params.fastmode = True - # params.fastmode = False - # params.readCompact = True - params.readCompact = False params.simToLoad= steering_categ.name.lower() + params.spline_dom_efficiency = False + params.spline_hole_ice = False + params.spline_anisotrophy = False params.evalThreads = args.threads # params.evalThreads = thread_factors(args.threads)[1] - params.spline_hole_ice = True - params.spline_dom_efficiency = True + params.diffuse_fit_type = gf.DiffuseFitType.SinglePowerLaw return params @@ -50,9 +44,12 @@ def set_up_as(fitter, params): def get_llh(fitter, params): fitparams = gf.FitParameters(gf.sampleTag.HESE) + # print params for parm in params: fitparams.__setattr__(parm.name, parm.value) - return fitter.EvalLLH(fitparams) + llh = -fitter.EvalLLH(fitparams) + # print '=== llh = {0}'.format(llh) + return llh def data_distributions(fitter): diff --git a/utils/likelihood.py b/utils/likelihood.py index 9bba79a..f590200 100644 --- a/utils/likelihood.py +++ b/utils/likelihood.py @@ -120,6 +120,8 @@ def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter): for idx, param in enumerate(hypo_paramset.from_tag(ParamTag.BESTFIT)): param.value = fr[idx] + # print 'hypo_paramset', hypo_paramset + if args.likelihood is Likelihood.FLAT: return 1. elif args.likelihood is Likelihood.GAUSSIAN: @@ -127,3 +129,12 @@ def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter): return gaussian_llh(fr, fr_bf, args.sigma_ratio) elif args.likelihood is Likelihood.GOLEMFIT: return gf_utils.get_llh(fitter, hypo_paramset) + +def ln_prob(theta, args, fitter, asimov_paramset, mcmc_paramset): + lp = lnprior(theta, paramset=mcmc_paramset) + if not np.isfinite(lp): + return -np.inf + return lp + triangle_llh( + theta, args=args, asimov_paramset=asimov_paramset, + mcmc_paramset=mcmc_paramset, fitter=fitter + ) diff --git a/utils/mcmc.py b/utils/mcmc.py index f898b83..c712c3a 100644 --- a/utils/mcmc.py +++ b/utils/mcmc.py @@ -21,10 +21,10 @@ from utils.enums import MCMCSeedType from utils.misc import enum_parse, make_dir, parse_bool -def mcmc(p0, triangle_llh, lnprior, ndim, nwalkers, burnin, nsteps, ntemps=1, threads=1): +def mcmc(p0, ln_prob, ndim, nwalkers, burnin, nsteps, threads=1): """Run the MCMC.""" - sampler = emcee.PTSampler( - ntemps, nwalkers, ndim, triangle_llh, lnprior, threads=threads + sampler = emcee.EnsembleSampler( + nwalkers, ndim, ln_prob, threads=threads ) print "Running burn-in" @@ -38,7 +38,7 @@ def mcmc(p0, triangle_llh, lnprior, ndim, nwalkers, burnin, nsteps, ntemps=1, th pass print "Finished" - samples = sampler.chain[0, :, :, :].reshape((-1, ndim)) + samples = sampler.chain.reshape((-1, ndim)) print 'acceptance fraction', sampler.acceptance_fraction print 'sum of acceptance fraction', np.sum(sampler.acceptance_fraction) print 'np.unique(samples[:,0]).shape', np.unique(samples[:,0]).shape @@ -74,22 +74,22 @@ def mcmc_argparse(parser): ) -def flat_seed(paramset, ntemps, nwalkers): +def flat_seed(paramset, nwalkers): """Get gaussian seed values for the MCMC.""" ndim = len(paramset) low = np.array(paramset.seeds).T[0] high = np.array(paramset.seeds).T[1] p0 = np.random.uniform( - low=low, high=high, size=[ntemps, nwalkers, ndim] + low=low, high=high, size=[nwalkers, ndim] ) return p0 -def gaussian_seed(paramset, ntemps, nwalkers): +def gaussian_seed(paramset, nwalkers): """Get gaussian seed values for the MCMC.""" ndim = len(paramset) p0 = np.random.normal( - paramset.values, paramset.stds, size=[ntemps, nwalkers, ndim] + paramset.values, paramset.stds, size=[nwalkers, ndim] ) return p0 diff --git a/utils/misc.py b/utils/misc.py index 6a81bbf..a189f07 100644 --- a/utils/misc.py +++ b/utils/misc.py @@ -88,13 +88,14 @@ class ParamSet(Sequence): except TypeError: param_sequence.append(arg) - # Disallow duplicated params - all_names = [p.name for p in param_sequence] - unique_names = set(all_names) - if len(unique_names) != len(all_names): - duplicates = set([x for x in all_names if all_names.count(x) > 1]) - raise ValueError('Duplicate definitions found for param(s): ' + - ', '.join(str(e) for e in duplicates)) + if len(param_sequence) != 0: + # Disallow duplicated params + all_names = [p.name for p in param_sequence] + unique_names = set(all_names) + if len(unique_names) != len(all_names): + duplicates = set([x for x in all_names if all_names.count(x) > 1]) + raise ValueError('Duplicate definitions found for param(s): ' + + ', '.join(str(e) for e in duplicates)) # Elements of list must be Param type assert all([isinstance(x, Param) for x in param_sequence]), \ diff --git a/utils/plot.py b/utils/plot.py index a659b78..87a6f5c 100644 --- a/utils/plot.py +++ b/utils/plot.py @@ -28,6 +28,43 @@ rc('text', usetex=False) rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) +def calc_nbins(x): + n = (np.max(x) - np.min(x)) / (2 * len(x)**(-1./3) * (np.percentile(x, 75) - np.percentile(x, 25))) + return np.floor(n) + + +def most_likely(arr): + """Return the densest region given a 1D array of data.""" + binning = calc_bins(arr) + harr = np.histogram(arr, binning)[0] + return centers(binning)[np.argmax(harr)] + + +def interval(arr, percentile=68.): + """Returns the *percentile* shortest interval around the mode.""" + center = most_likely(arr) + sarr = sorted(arr) + delta = np.abs(sarr - center) + curr_low = np.argmin(delta) + curr_up = curr_low + npoints = len(sarr) + while curr_up - curr_low < percentile/100.*npoints: + if curr_low == 0: + curr_up += 1 + elif curr_up == npoints-1: + curr_low -= 1 + elif sarr[curr_up]-sarr[curr_low-1] < sarr[curr_up+1]-sarr[curr_low]: + curr_low -= 1 + elif sarr[curr_up]-sarr[curr_low-1] > sarr[curr_up+1]-sarr[curr_low]: + curr_up += 1 + elif (curr_up - curr_low) % 2: + # they are equal so step half of the time up and down + curr_low -= 1 + else: + curr_up += 1 + return sarr[curr_low], center, sarr[curr_up] + + def plot_argparse(parser): """Arguments for plotting.""" parser.add_argument( @@ -136,21 +173,28 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset): g.export(outfile+'_angles.'+of) if args.plot_elements: - nu_index = mcmc_paramset.idx_from_tag(ParamTag.NUISANCE) - fr_index = mcmc_paramset.idx_from_tag(ParamTag.MMANGLES) - sc_index = mcmc_paramset.idx_from_tag(ParamTag.SCALE) - sr_index = mcmc_paramset.idx_from_tag(ParamTag.SRCANGLES) + nu_index = mcmc_paramset.from_tag(ParamTag.NUISANCE, index=True) + fr_index = mcmc_paramset.from_tag(ParamTag.MMANGLES, index=True) + sc_index = mcmc_paramset.from_tag(ParamTag.SCALE, index=True) + if not args.fix_source_ratio: + sr_index = mcmc_paramset.from_tag(ParamTag.SRCANGLES, index=True) nu_elements = raw[:,nu_index] - fr_elements = np.array(map(angles_to_fr, raw[:,fr_index])) + fr_elements = np.array(map(flat_angles_to_u, raw[:,fr_index])) sc_elements = raw[:,sc_index] - sr_elements = np.array(map(flat_angles_to_u, raw[:,sr_index])) - Tchain = np.column_stack( - [nu_elements, fr_elements, sc_elements, sr_elements] - ) + if not args.fix_source_ratio: + sr_elements = np.array(map(angles_to_fr, raw[:,sr_index])) + if args.fix_source_ratio: + Tchain = np.column_stack( + [nu_elements, fr_elements, sc_elements] + ) + else: + Tchain = np.column_stack( + [nu_elements, fr_elements, sc_elements, sr_elements] + ) - trns_ranges = ranges[nu_index,] - trns_axes_labels = axes_labels[nu_index,] + trns_ranges = np.array(ranges)[nu_index,].tolist() + trns_axes_labels = np.array(axes_labels)[nu_index,].tolist() if not args.fix_mixing: trns_axes_labels += \ [r'\mid \tilde{U}_{e1} \mid' , r'\mid \tilde{U}_{e2} \mid' , r'\mid \tilde{U}_{e3} \mid' , \ @@ -158,8 +202,8 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset): r'\mid \tilde{U}_{\tau1} \mid' , r'\mid \tilde{U}_{\tau2} \mid' , r'\mid \tilde{U}_{\tau3} \mid'] trns_ranges += [(0, 1)] * 9 if not args.fix_scale: - trns_axes_labels += [axes_labels[sc_index]] - trns_ranges += [ranges[sc_index]] + trns_axes_labels += [np.array(axes_labels)[sc_index].tolist()] + trns_ranges += [np.array(ranges)[sc_index].tolist()] if not args.fix_source_ratio: trns_axes_labels += [r'\phi_e', r'\phi_\mu', r'\phi_\tau'] trns_ranges += [(0, 1)] * 3 @@ -171,4 +215,4 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset): else: mpl.pyplot.figtext(0.5, 0.7, fig_text, fontsize=15) for of in outformat: - g.export(outfile+'_angles.'+of) + g.export(outfile+'_elements.'+of) |