diff options
Diffstat (limited to 'utils')
| -rw-r--r-- | utils/fr.py | 2 | ||||
| -rw-r--r-- | utils/gf.py | 5 | ||||
| -rw-r--r-- | utils/likelihood.py | 4 | ||||
| -rw-r--r-- | utils/misc.py | 30 | ||||
| -rw-r--r-- | utils/multinest.py | 10 | ||||
| -rw-r--r-- | utils/param.py | 14 | ||||
| -rw-r--r-- | utils/plot.py | 381 |
7 files changed, 349 insertions, 97 deletions
diff --git a/utils/fr.py b/utils/fr.py index 61b5a61..5acd6f8 100644 --- a/utils/fr.py +++ b/utils/fr.py @@ -212,7 +212,7 @@ def estimate_scale(args): upper_s = (m_eign/args.binning[0]) / (args.binning[0]**(args.dimension-3)) scale = np.power(10, np.average(np.log10([lower_s, upper_s]))) diff = upper_s / lower_s - scale_region = (lower_s/diff, upper_s*diff) + scale_region = (lower_s/np.power(10, args.dimension), upper_s*diff*np.power(10, args.dimension)) scale_region = [np.power(10, np.round(np.log10(x))) for x in scale_region] return scale, scale_region diff --git a/utils/gf.py b/utils/gf.py index ea4f61f..098faa9 100644 --- a/utils/gf.py +++ b/utils/gf.py @@ -57,9 +57,14 @@ def steering_params(args): steering_categ = args.ast params = gf.SteeringParams() params.quiet = False + params.fastmode = True params.simToLoad= steering_categ.name.lower() params.evalThreads = args.threads # params.evalThreads = thread_factors(args.threads)[1] + + # For Tianlu + # params.years = [999] + return params diff --git a/utils/likelihood.py b/utils/likelihood.py index 2e4a22d..c2822a1 100644 --- a/utils/likelihood.py +++ b/utils/likelihood.py @@ -167,9 +167,6 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter): for idx, param in enumerate(hypo_paramset.from_tag(ParamTag.BESTFIT)): param.value = flavour_angles[idx] - print 'llh_paramset', llh_paramset - print 'hypo_paramset', hypo_paramset - print 'fr', fr if args.likelihood is Likelihood.FLAT: llh = 1. elif args.likelihood is Likelihood.GAUSSIAN: @@ -179,7 +176,6 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter): llh = gf_utils.get_llh(fitter, hypo_paramset) elif args.likelihood is Likelihood.GF_FREQ: lhh = gf_utils.get_llh_freq(fitter, hypo_paramset) - print 'llh', llh return llh diff --git a/utils/misc.py b/utils/misc.py index cad03bc..5483675 100644 --- a/utils/misc.py +++ b/utils/misc.py @@ -12,6 +12,7 @@ from __future__ import absolute_import, division import os import errno import multiprocessing +from fractions import gcd import argparse from operator import attrgetter @@ -33,14 +34,18 @@ class SortingHelpFormatter(argparse.HelpFormatter): super(SortingHelpFormatter, self).add_arguments(actions) +def solve_ratio(fr): + denominator = reduce(gcd, fr) + return [int(x/denominator) for x in fr] + + def gen_identifier(args): f = '_DIM{0}'.format(args.dimension) - mr1, mr2, mr3 = args.measured_ratio + mr1, mr2, mr3 = solve_ratio(args.measured_ratio) if args.fix_source_ratio: - sr1, sr2, sr3 = args.source_ratio - f += '_sfr_{0:03d}_{1:03d}_{2:03d}_mfr_{3:03d}_{4:03d}_{5:03d}'.format( - int(sr1*100), int(sr2*100), int(sr3*100), - int(mr1*100), int(mr2*100), int(mr3*100) + sr1, sr2, sr3 = solve_ratio(args.source_ratio) + f += '_sfr_{0:G}_{1:G}_{2:G}_mfr_{3:G}_{4:G}_{5:G}'.format( + sr1, sr2, sr3, mr1, mr2, mr3 ) if args.fix_mixing: f += '_fix_mixing' @@ -128,6 +133,21 @@ def make_dir(outfile): else: raise +def remove_option(parser, arg): + for action in parser._actions: + if (vars(action)['option_strings'] + and vars(action)['option_strings'][0] == arg) \ + or vars(action)['dest'] == arg: + parser._remove_action(action) + + for action in parser._action_groups: + vars_action = vars(action) + var_group_actions = vars_action['_group_actions'] + for x in var_group_actions: + if x.dest == arg: + var_group_actions.remove(x) + return + def seed_parse(s): if s.lower() == 'none': diff --git a/utils/multinest.py b/utils/multinest.py index 9dd0742..f3595f9 100644 --- a/utils/multinest.py +++ b/utils/multinest.py @@ -16,7 +16,7 @@ import numpy as np from pymultinest import analyse, run from utils import likelihood -from utils.misc import gen_outfile_name, make_dir +from utils.misc import gen_identifier, make_dir def CubePrior(cube, ndim, n_params): @@ -37,13 +37,15 @@ def lnProb(cube, ndim, n_params, mn_paramset, llh_paramset, asimov_paramset, llh_paramset[pm].value = mn_paramset[pm].value theta = llh_paramset.values # print 'llh_paramset', llh_paramset - return likelihood.ln_prob( + llh = likelihood.ln_prob( theta=theta, args=args, asimov_paramset=asimov_paramset, llh_paramset=llh_paramset, fitter=fitter ) + # print 'llh', llh + return llh def mn_argparse(parser): @@ -77,8 +79,8 @@ def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter): fitter = fitter ) - prefix = './mnrun/DIM{0}/{1}_{2:>019}_'.format( - args.dimension, gen_outfile_name(args), np.random.randint(0, 2**63) + prefix = './mnrun/DIM{0}/{1}_{2:>010}_'.format( + args.dimension, gen_identifier(args), np.random.randint(0, 2**30) ) make_dir(prefix) print 'Running evidence calculation for {0}'.format(prefix) diff --git a/utils/param.py b/utils/param.py index f861264..25558c0 100644 --- a/utils/param.py +++ b/utils/param.py @@ -252,10 +252,16 @@ def get_paramsets(args, nuisance_paramset): ]) if not args.fix_scale: tag = ParamTag.SCALE - llh_paramset.append( - Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3, - tex=r'{\rm log}_{10}\Lambda^{-1}'+get_units(args.dimension), tag=tag) - ) + if hasattr(args, 'dimension'): + llh_paramset.append( + Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3, + tex=r'{\rm log}_{10}\Lambda^{-1}'+get_units(args.dimension), tag=tag) + ) + elif hasattr(args, 'dimensions'): + llh_paramset.append( + Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3, + tex=r'{\rm log}_{10}\Lambda^{-1} / GeV^{-d+4}', tag=tag) + ) if not args.fix_source_ratio: tag = ParamTag.SRCANGLES llh_paramset.extend([ diff --git a/utils/plot.py b/utils/plot.py index 3d94cc1..8f4eba8 100644 --- a/utils/plot.py +++ b/utils/plot.py @@ -10,6 +10,7 @@ Plotting functions for the BSM flavour ratio analysis from __future__ import absolute_import, division import os +from copy import deepcopy import numpy as np import matplotlib as mpl @@ -112,18 +113,18 @@ def plot_Tchain(Tchain, axes_labels, ranges): def gen_figtext(args): """Generate the figure text.""" t = '' - mr1, mr2, mr3 = args.measured_ratio + mr1, mr2, mr3 = misc_utils.solve_ratio(args.measured_ratio) if args.fix_source_ratio: - sr1, sr2, sr3 = args.source_ratio - t += 'Source flavour ratio = [{0:.2f}, {1:.2f}, {2:.2f}]\nIC ' \ - 'observed flavour ratio = [{3:.2f}, {4:.2f}, ' \ - '{5:.2f}]\nDimension = {6}'.format( + sr1, sr2, sr3 = misc_utils.solve_ratio(args.source_ratio) + t += 'Source flavour ratio = [{0}, {1}, {2}]\nIC ' \ + 'observed flavour ratio = [{3}, {4}, ' \ + '{5}]\nDimension = {6}'.format( sr1, sr2, sr3, mr1, mr2, mr3, args.dimension, int(args.energy) ) else: - t += 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, ' \ - '{2:.2f}]\nDimension = {3}'.format( + t += 'IC observed flavour ratio = [{0}, {1}, ' \ + '{2}]\nDimension = {3}'.format( mr1, mr2, mr3, args.dimension, int(args.energy) ) if args.fix_scale: @@ -229,28 +230,33 @@ def chainer_plot(infile, outfile, outformat, args, llh_paramset): g.export(outfile+'_elements.'+of) +def myround(x, base=5, up=False, down=False): + if up == down and up is True: assert 0 + if up: return int(base * np.round(float(x)/base-0.5)) + elif down: return int(base * np.round(float(x)/base+0.5)) + else: int(base * np.round(float(x)/base)) + + def plot_statistic(data, outfile, outformat, args, scale_param, label=None): """Make MultiNest factor or LLH value plot.""" - print "Making Statistic plot" + print 'Making Statistic plot' fig_text = gen_figtext(args) if label is not None: fig_text += '\n' + label print 'data', data print 'data.shape', data.shape scales, statistic = data.T + min_idx = np.argmin(scales) + null = statistic[min_idx] if args.stat_method is StatCateg.BAYESIAN: - min_idx = np.argmin(scales) - null = statistic[min_idx] reduced_ev = -(statistic - null) elif args.stat_method is StatCateg.FREQUENTIST: - min_idx = np.argmin(scales) - null = statistic[min_idx] reduced_ev = -2*(statistic - null) fig = plt.figure(figsize=(7, 5)) ax = fig.add_subplot(111) - ax.set_xlim(scale_param.ranges) + ax.set_xlim(np.log10(args.scale_region)) ax.set_xlabel('$'+scale_param.tex+'$') if args.stat_method is StatCateg.BAYESIAN: ax.set_ylabel(r'Bayes Factor') @@ -271,85 +277,302 @@ def plot_statistic(data, outfile, outformat, args, scale_param, label=None): ax.add_artist(at) misc_utils.make_dir(outfile) for of in outformat: + print 'Saving as {0}'.format(outfile+'.'+of) fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150) -def myround(x, base=5, up=False, down=False): - if up == down and up is True: assert 0 - if up: return int(base * np.round(float(x)/base-0.5)) - elif down: return int(base * np.round(float(x)/base+0.5)) - else: int(base * np.round(float(x)/base)) +def plot_sens_full(data, outfile, outformat, args): + print 'Making FULL sensitivity plot' - -def plot_BSM_angles_limit(dirname, outfile, outformat, args, bayesian): - """Make BSM angles vs scale limit plot.""" - if not args.plot_angles_limit: return - print "Making BSM angles limit plot.""" - fig_text = gen_figtext(args) - xticks = [r'$\mathcal{O}_{12}$', r'$\mathcal{O}_{13}$', r'$\mathcal{O}_{23}$'] - - raw = [] - for root, dirs, filenames in os.walk(dirname): - for fn in filenames: - if fn[-4:] == '.npy': - raw.append(np.load(os.path.join(root, fn))) - raw = np.vstack(raw) - raw = raw[np.argsort(raw[:,0])] - print 'raw', raw - print 'raw.shape', raw.shape - sc_ranges = ( - myround(np.min(raw[0][:,0]), up=True), - myround(np.max(raw[0][:,0]), down=True) - ) - - proc = [] - if bayesian: - for idx, theta in enumerate(raw): - scale, evidences = theta.T - null = evidences[0] - reduced_ev = -(evidences - null) - al = scale[reduced_ev > np.log(10**(1/2.))] - proc.append((idx+1, al[0])) - else: - for idx, theta in enumerate(raw): - scale, llh = theta.T - delta_llh = -2 * (llh - np.max(llh)) - # 90% CL for 1 dof - al = scale[delta_llh > 2.71] - proc.append((idx+1, al[0])) - - limits = np.array(proc) - print 'limits', limits + colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'} + xticks = ['{0}'.format(x) for x in range(np.min(args.dimensions), + np.max(args.dimensions)+1)] + yranges = [np.inf, -np.inf] + legend_handles = [] fig = plt.figure(figsize=(7, 5)) ax = fig.add_subplot(111) - ax.set_xticklabels(['']+xticks+['']) - ax.set_xlim(0, len(xticks)+1) - ax.set_ylim(sc_ranges[0], sc_ranges[-1]) - ax.set_xlabel(r'BSM angle') - ylabel = r'${\rm log}_{10} \Lambda' + get_units(args.dimension) + r'$' - ax.set_ylabel(ylabel) + ax.set_xlim(args.dimensions[0]-1, args.dimensions[-1]+1) + ax.set_xticklabels([''] + xticks + ['']) + ax.set_xlabel(r'BSM operator dimension ' + r'$d$') + ax.set_ylabel(r'${\rm log}_{10} \Lambda^{-1} / GeV^{-d+4}$') + + argsc = deepcopy(args) + for idim in xrange(len(data)): + dim = args.dimensions[idim] + print 'dim', dim + argsc.dimension = dim + for isrc in xrange(len(data[idim])): + src = args.source_ratios[isrc] + argsc.source_ratio = src + + # fig_text = gen_figtext(argsc) + scales, statistic = data[idim][isrc].T + min_idx = np.argmin(scales) + null = statistic[min_idx] + if args.stat_method is StatCateg.BAYESIAN: + reduced_ev = -(statistic - null) + # al = scales[reduced_ev > np.log(10**(3/2.))] # Strong degree of belief + al = scales[reduced_ev > 0.4] # Testing + elif args.stat_method is StatCateg.FREQUENTIST: + reduced_ev = -2*(statistic - null) + al = scales[reduced_ev > 2.71] # 90% CL for 1 DOF via Wilks + if len(al) == 0: + print 'No points for DIM {0} FRS {1} NULL {2}!'.format( + dim, misc_utils.solve_ratio(src), null + ) + print 'Reduced EV {0}'.format(reduced_ev) + continue + lim = al[0] + label = '[{0}, {1}, {2}]'.format(*misc_utils.solve_ratio(src)) + if lim < yranges[0]: yranges[0] = lim + if lim > yranges[1]: yranges[1] = lim+4 + line = plt.Line2D( + (dim-0.1, dim+0.1), (lim, lim), lw=3, color=colour[isrc], label=label + ) + ax.add_line(line) + if idim == 0: legend_handles.append(line) + x_offset = isrc*0.05 - 0.05 + ax.annotate( + s='', xy=(dim+x_offset, lim), xytext=(dim+x_offset, lim+3), + arrowprops={'arrowstyle': '<-', 'lw': 1.2, 'color':colour[isrc]} + ) - for l in limits: - line = plt.Line2D( - (l[0]-0.1, l[0]+0.1), (l[1], l[1]), lw=3, color='r' - ) - ax.add_line(line) - # ax.arrow( - # l[0], l[1], 0, -1.5, head_width=0.05, head_length=0.2, fc='r', - # ec='r', lw=2 - # ) - ax.annotate( - s='', xy=l, xytext=(l[0], l[1]+1.5), - arrowprops={'arrowstyle': '<-', 'lw': 1.5, 'color':'r'} - ) + try: + yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True)) + ax.set_ylim(yranges) + except: pass + ax.legend(handles=legend_handles, prop=dict(size=8), loc='upper right') for ymaj in ax.yaxis.get_majorticklocs(): - ax.axhline(y=ymaj, ls='-', color='gray', alpha=0.4, linewidth=1) + ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.4, linewidth=1) for xmaj in ax.xaxis.get_majorticklocs(): - ax.axvline(x=xmaj, ls='-', color='gray', alpha=0.4, linewidth=1) + ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1) + misc_utils.make_dir(outfile) for of in outformat: + print 'Saving plot as {0}'.format(outfile+'.'+of) fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150) + +def plot_sens_fixed_angle(data, outfile, outformat, args): + print 'Making FIXED_ANGLE sensitivity plot' + + colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'} + xticks = [r'$\mathcal{O}_{12}$', r'$\mathcal{O}_{13}$', r'$\mathcal{O}_{23}$'] + argsc = deepcopy(args) + for idim in xrange(len(data)): + dim = args.dimensions[idim] + print '= dim', dim + argsc.dimension = dim + + yranges = [np.inf, -np.inf] + legend_handles = [] + + fig = plt.figure(figsize=(7, 5)) + ax = fig.add_subplot(111) + ax.set_xlim(0, len(xticks)+1) + ax.set_xticklabels([''] + xticks + ['']) + ax.set_xlabel(r'BSM operator angle') + ax.set_ylabel(r'${\rm log}_{10} \Lambda^{-1} / GeV^{-d+4}$') + + for isrc in xrange(len(data[idim])): + src = args.source_ratios[isrc] + argsc.source_ratio = src + print '== src', src + + for ian in xrange(len(data[idim][isrc])): + print '=== an', ian + scales, statistic = data[idim][isrc][ian].T + min_idx = np.argmin(scales) + null = statistic[min_idx] + if args.stat_method is StatCateg.BAYESIAN: + reduced_ev = -(statistic - null) + al = scales[reduced_ev > np.log(10**(3/2.))] # Strong degree of belief + elif args.stat_method is StatCateg.FREQUENTIST: + reduced_ev = -2*(statistic - null) + al = scales[reduced_ev > 2.71] # 90% CL for 1 DOF via Wilks + print 'reduced_ev', reduced_ev + if len(al) == 0: + print 'No points for DIM {0} FRS {1}\nreduced_ev {2}!'.format( + dim, src, reduced_ev + ) + continue + lim = al[0] + print 'limit = {0}'.format(lim) + label = '[{0}, {1}, {2}]'.format(*misc_utils.solve_ratio(src)) + if lim < yranges[0]: yranges[0] = lim + if lim > yranges[1]: yranges[1] = lim+4 + line = plt.Line2D( + (ian+1-0.1, ian+1+0.1), (lim, lim), lw=3, color=colour[isrc], label=label + ) + ax.add_line(line) + if len(legend_handles) < isrc+1: + legend_handles.append(line) + x_offset = isrc*0.05 - 0.05 + ax.annotate( + s='', xy=(ian+1+x_offset, lim), xytext=(ian+1+x_offset, lim+3), + arrowprops={'arrowstyle': '<-', 'lw': 1.2, 'color':colour[isrc]} + ) + + try: + yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True)) + ax.set_ylim(yranges) + except: pass + + ax.legend(handles=legend_handles, prop=dict(size=8), loc='upper right') + for ymaj in ax.yaxis.get_majorticklocs(): + ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.4, linewidth=1) + for xmaj in ax.xaxis.get_majorticklocs(): + ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1) + + out = outfile + '_DIM{0}'.format(dim) + misc_utils.make_dir(out) + for of in outformat: + print 'Saving plot as {0}'.format(out+'.'+of) + fig.savefig(out+'.'+of, bbox_inches='tight', dpi=150) + + + for ymaj in ax.yaxis.get_majorticklocs(): + ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.4, linewidth=1) + for xmaj in ax.xaxis.get_majorticklocs(): + ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1) + + +def plot_sens_corr_angle(data, outfile, outformat, args): + print 'Making CORR_ANGLE sensitivity plot' + + labels = [r'$sin^2(2\mathcal{O}_{12})$', + r'$sin^2(2\mathcal{O}_{13})$', + r'$sin^2(2\mathcal{O}_{23})$'] + + argsc = deepcopy(args) + for idim in xrange(len(data)): + dim = args.dimensions[idim] + print '= dim', dim + argsc.dimension = dim + for isrc in xrange(len(data[idim])): + src = args.source_ratios[isrc] + argsc.source_ratio = src + print '== src', src + for ian in xrange(len(data[idim][isrc])): + print '=== an', ian + + d = data[idim][isrc][ian].reshape( + len(data[idim][isrc][ian])**2, 3 + ) + + fig = plt.figure(figsize=(7, 5)) + ax = fig.add_subplot(111) + ax.set_ylim(0, 1) + ax.set_ylabel(labels[ian]) + ax.set_xlabel(r'${\rm log}_{10} \Lambda^{-1}'+get_units(dim)+r'$') + + xranges = [np.inf, -np.inf] + legend_handles = [] + + y = d[:,0] + x = d[:,1] + z = d[:,2] + + print 'x', x + print 'y', y + print 'z', z + null_idx = np.argmin(x) + null = z[null_idx] + print 'null = {0}, for scale = {1}'.format(null, x[null_idx]) + + if args.stat_method is StatCateg.BAYESIAN: + z_r = -(z - null) + elif args.stat_method is StatCateg.FREQUENTIST: + z_r = -2*(z - null) + print 'scale', x + print 'reduced ev', z_r + + pdat = np.array([x, y, z_r, np.ones(x.shape)]).T + print 'pdat', pdat + pdat_clean = [] + for d in pdat: + if not np.any(np.isnan(d)): pdat_clean.append(d) + pdat = np.vstack(pdat_clean) + sort_column = 3 + pdat_sorted = pdat[pdat[:,sort_column].argsort()] + uni, c = np.unique(pdat[:,sort_column], return_counts=True) + print uni, c + print len(uni) + print np.unique(c) + + n = len(uni) + assert len(np.unique(c)) == 1 + c = c[0] + col_array = [] + for col in pdat_sorted.T: + col_array.append(col.reshape(n, c)) + col_array = np.stack(col_array) + sep_arrays = [] + for x_i in xrange(n): + sep_arrays.append(col_array[:,x_i]) + + print len(sep_arrays) + sep_arrays = sep_arrays[0][:3] + print sep_arrays + + if args.stat_method is StatCateg.BAYESIAN: + reduced_pdat_mask = (sep_arrays[2] > log(10**(3/2.))) # Strong degree of belief + elif args.stat_method is StatCateg.FREQUENTIST: + reduced_pdat_mask = (sep_arrays[2] > 4.61) # 90% CL for 2 DOFS via Wilks + reduced_pdat = sep_arrays.T[reduced_pdat_mask].T + print 'reduced_pdat', reduced_pdat + + ax.tick_params(axis='x', labelsize=11) + ax.tick_params(axis='y', labelsize=11) + + mini, maxi = np.min(x), np.max(x) + ax.set_xlim((mini, maxi)) + ax.set_ylim(0, 1) + ax.grid(b=False) + + x_v = reduced_pdat[0].round(decimals=4) + y_v = reduced_pdat[1].round(decimals=4) + uniques = np.unique(x_v) + print 'uniques', uniques + if len(uniques) == 1: continue + bw = np.min(np.diff(uniques)) + print 'bw', bw + print np.diff(uniques) + uni_x_split = np.split(uniques, np.where(np.diff(uniques) > bw*(1e20))[0] + 1) + print 'len(uni_x_split)', len(uni_x_split) + for uni_x in uni_x_split: + p_x_l, p_y_l = [], [] + p_x_u, p_y_u = [], [] + for uni in uni_x: + idxes = np.where(x_v == uni)[0] + ymin, ymax = 1, 0 + for idx in idxes: + if y_v[idx] < ymin: ymin = y_v[idx] + if y_v[idx] > ymax: ymax = y_v[idx] + p_x_l.append(uni) + p_y_l.append(ymin) + p_x_u.append(uni) + p_y_u.append(ymax) + p_x_l, p_y_l = np.array(p_x_l, dtype=np.float64), np.array(p_y_l, dtype=np.float64) + p_x_u, p_y_u = np.array(list(reversed(p_x_u)), dtype=np.float64), np.array(list(reversed(p_y_u)), dtype=np.float64) + p_x = np.hstack([p_x_l, p_x_u]) + p_y = np.hstack([p_y_l, p_y_u]) + p_x = np.r_[p_x, p_x[0]] + p_y = np.r_[p_y, p_y[0]] + print 'p_x', p_x + print 'p_y', p_y + try: + tck, u = interpolate.splprep([p_x, p_y], s=1e-5, per=True) + xi, yi = interpolate.splev(np.linspace(0, 1, 1000), tck) + except: + xi, yi = p_x, p_y + ax.fill(xi, yi, 'r', edgecolor='r', linewidth=1) + + out = outfile + '_DIM{0}_SRC{1}_AN{2}'.format(dim, isrc, ian) + misc_utils.make_dir(out) + for of in outformat: + print 'Saving plot as {0}'.format(out+'.'+of) + fig.savefig(out+'.'+of, bbox_inches='tight', dpi=150) |