diff options
Diffstat (limited to 'utils/plot.py')
| -rw-r--r-- | utils/plot.py | 72 |
1 files changed, 58 insertions, 14 deletions
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) |