Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

added 2d pdf plotting, flux range options; made percentile/prob excee…

…d more robust
  • Loading branch information...
commit c7be5501a84c25475321a6872c29fe5c339f2a47 1 parent 6f19720
@eric-switzer authored
View
6 data/spt_catalog.ini
@@ -18,7 +18,10 @@ omega_prior = 1.
range_alpha = [-5., 5.]
prior_alpha = [-3., 5.]
# use a dN/dS model identical to that used in SPT
-use_spt_model = True
+use_spt_model = False
+dnds_minflux = 1.e-8
+dnds_maxflux = 1.5
+dnds_numflux = 1e5
[FieldNames]
# this identifies the flux and error columns in the catalog
@@ -37,3 +40,4 @@ comparison_catalog = ../data/source_table_vieira09_3sigma.dat
augmented_catalog = ../data/augmented_spt_catalog.shelve
catalog_in_mjy = True
verbose = True
+make_2dplot = True
View
101 python/plot_2d_pdf.py
@@ -0,0 +1,101 @@
+import numpy as np
+import tempfile
+import subprocess
+import matplotlib
+
+
+def gnuplot_2D(outfilename, region, xaxis, yaxis, xylabels,
+ aspect, fulltitle, cbar_title, eps_outfile=None,
+ draw_objects=[], density='300', size=5., logscale=False):
+ r"""plot a 2D matrix
+ """
+
+ print fulltitle, repr(region.shape)
+
+ if logscale:
+ # TODO: this is a kludge
+ region += np.max(region)/1.e6
+ try:
+ region = np.log10(region)
+ except FloatingPointError:
+ print "plotting error (all zeros, log requested?)"
+ region = np.zeros_like(region)
+
+ input_data_file = tempfile.NamedTemporaryFile()
+ gplfile = tempfile.NamedTemporaryFile(suffix=".gpl")
+ outplot_file = tempfile.NamedTemporaryFile(suffix=".eps")
+
+ deltax = abs(xaxis[1] - xaxis[0])
+ deltay = abs(yaxis[1] - yaxis[0])
+ xminmax = (np.min(xaxis) - 0.5 * deltax, np.max(xaxis) + 0.5 * deltax)
+ yminmax = (np.min(yaxis) - 0.5 * deltay, np.max(yaxis) + 0.5 * deltay)
+ aspect = (yminmax[1] - yminmax[0]) / (xminmax[1] - xminmax[0])
+
+ # TODO: optimize this
+ for xind in range(len(xaxis)):
+ for yind in range(len(yaxis)):
+ outstring = "%g %g %g\n" % \
+ (xaxis[xind], yaxis[yind], region[xind, yind])
+ input_data_file.write(outstring)
+
+ input_data_file.flush()
+
+ gplfile.write("set view map\n")
+ gplfile.write("set xrange [%g:%g]\n" % xminmax)
+ gplfile.write("set yrange [%g:%g]\n" % yminmax)
+ #gplfile.write("set cbrange [%g:%g]\n" % cminmax)
+ gplfile.write("set size ratio -1\n")
+ gplfile.write("set nokey\n")
+ gplfile.write("set tics out\n")
+ gplfile.write('set xlabel "%s"\n' % xylabels[0])
+ gplfile.write('set ylabel "%s"\n' % xylabels[1])
+ gplfile.write('set title "%s"\n' % fulltitle)
+ gplfile.write('set cblabel "%s"\n' % cbar_title)
+ gplfile.write('set rmargin 10\n')
+ gplfile.write("set terminal postscript eps color size %g, %g\n" % \
+ (size, size * aspect * 1.1))
+
+ if eps_outfile is None:
+ gplfile.write('set output "%s"\n' % outplot_file.name)
+ else:
+ gplfile.write('set output "%s"\n' % eps_outfile)
+
+ objnum = 10
+ for drawing in draw_objects:
+ drawing["objnum"] = objnum
+
+ objcmd = "set obj %(objnum)d %(primitive)s " % drawing
+ objcmd += "at graph %(center_x)g, %(center_y)g " % drawing
+
+ if drawing["primitive"] == "circle":
+ objcmd += "size %(radius)g front\n" % drawing
+
+ if drawing["primitive"] == "rect":
+ objcmd += "size %(size_x)g, %(size_y)g front\n" % drawing
+
+ gplfile.write(objcmd)
+
+ objcmd = "set obj %(objnum)d lw %(width)d " % drawing
+ objcmd += "fill empty border rgb \"%(color)s\"\n" % drawing
+ gplfile.write(objcmd)
+
+ objnum += 1
+
+ gplfile.write("set palette rgbformulae 22, 13, -31\n")
+
+ gplfile.write('plot "%s" using 1:2:3 with image\n' % input_data_file.name)
+ gplfile.flush()
+
+ gnuplot = '/cita/h/home-2/eswitzer/local/bin/gnuplot'
+ subprocess.check_call((gnuplot, gplfile.name))
+ gplfile.close()
+ input_data_file.close()
+
+ if eps_outfile is None:
+ subprocess.check_call(('convert', '-density', density, '-trim', '+repage',
+ '-border', '40x40', '-bordercolor', 'white',
+ outplot_file.name, outfilename))
+
+ outplot_file.close()
+
+
View
29 python/process_catalog.py
@@ -3,6 +3,7 @@
import source_count_models as dnds
import utilities as utils
import two_band_posterior_flux as tbpf
+import plot_2d_pdf
def process_ptsrc_catalog_alpha(catalog, gp):
@@ -53,7 +54,9 @@ def process_ptsrc_catalog_alpha(catalog, gp):
# TODO: this is really not very generic now
# frequencies are hard-coded into particular lookup tables
# TODO: why let this extend beyond the log-stepped axis (1.)?
- input_s_linear = np.linspace(1.e-8, 1.5, 1e5, endpoint=False)
+ input_s_linear = np.linspace(gp['dnds_minflux'],
+ gp['dnds_maxflux'],
+ gp['dnds_numflux'], endpoint=False)
if gp['use_spt_model']:
import spt_source_count_models as sptdnds
@@ -142,6 +145,30 @@ def process_ptsrc_catalog_alpha(catalog, gp):
source_entry["alpha_posterior_det"] = posterior_swap[2]
source_entry["prob_exceed_det"] = posterior_swap[3]
+ # plot the 2D posteriors
+ if gp['make_2dplot']:
+ full_package = posterior[4]
+ full_package_swap = posterior_swap[4]
+
+ if ((flux1 / sigma1) > (flux2 / sigma2)):
+ plot_2d_pdf.gnuplot_2D("../plots/%s_fluxflux.png" % srcname,
+ full_package["posterior_fluxflux"],
+ full_package["flux_axis"] * 1000.,
+ full_package["flux_axis"] * 1000.,
+ ["148 GHz Flux (mJy)",
+ "220 GHz Flux (mJy)"],
+ 1., srcname, "", logscale=False)
+
+ else:
+ plot_2d_pdf.gnuplot_2D(
+ "../plots/%s_fluxflux_swap.png" % srcname,
+ np.transpose(full_package_swap["posterior_fluxflux"]),
+ full_package_swap["flux_axis"] * 1000.,
+ full_package_swap["flux_axis"] * 1000.,
+ ["148 GHz Flux (mJy)",
+ "220 GHz Flux (mJy)"],
+ 1., srcname, "", logscale=False)
+
augmented_catalog[srcname] = source_entry
prefix = "The percentile points of the posterior "
View
11 python/two_band_posterior_flux.py
@@ -70,6 +70,8 @@ def two_band_posterior_flux(flux1, flux2, sigma1, sigma2, sigma12, s_in, dnds1,
print "joint code using cal covariance:" + repr(cov_calibration_jy)
print "joint code using covariance:" + repr(total_covariance)
+ full_package = {}
+
if swap_flux:
if gp['verbose']:
print "running band-2 selected source case (swapped)"
@@ -91,6 +93,7 @@ def two_band_posterior_flux(flux1, flux2, sigma1, sigma2, sigma12, s_in, dnds1,
#flux2_dist = np.sum(posterior_fluxindex, axis=1)
flux1_dist = np.sum(posterior_fluxflux, axis=0)
flux2_dist = np.sum(posterior_fluxflux, axis=1)
+
else:
if gp['verbose']:
print "running band-1 selected source case"
@@ -110,6 +113,12 @@ def two_band_posterior_flux(flux1, flux2, sigma1, sigma2, sigma12, s_in, dnds1,
flux1_dist = np.sum(posterior_fluxflux, axis=1)
flux2_dist = np.sum(posterior_fluxflux, axis=0)
+
+ full_package["posterior_fluxindex"] = posterior_fluxindex
+ full_package["posterior_fluxflux"] = posterior_fluxflux
+ full_package["flux_axis"] = flux_axis
+ full_package["alpha_axis"] = alpha_axis
+
# calculate the summaries of the various output PDFs
flux1_percentiles = utils.percentile_points(flux_axis, flux1_dist,
gp['percentiles'])
@@ -124,7 +133,7 @@ def two_band_posterior_flux(flux1, flux2, sigma1, sigma2, sigma12, s_in, dnds1,
gp['spectral_threshold'])
return (flux1_percentiles, flux2_percentiles, \
- alpha_percentiles, probexceed)
+ alpha_percentiles, probexceed, full_package)
# TODO: numpy.einsum may be able to do some operations faster
View
23 python/utilities.py
@@ -139,10 +139,13 @@ def percentile_points(axis, pdf, percentiles):
cdf_points = np.zeros(len(percentiles), dtype=float)
for index, percentile in enumerate(percentiles):
- cdf = np.cumsum(pdf) / float(np.sum(pdf))
- #print percentile
- #print cdf
- cdf_points[index] = np.nanmin(axis[cdf >= percentile])
+ try:
+ cdf = np.cumsum(pdf) / float(np.sum(pdf))
+ cdf_points[index] = np.nanmin(axis[cdf >= percentile])
+ except FloatingPointError:
+ print "percentile_points: invalid PDF"
+ cdf = np.zeros_like(pdf)
+ cdf_points[index] = -1.
return cdf_points
@@ -155,9 +158,15 @@ def fraction_exceeds(vector, threshold):
def prob_exceed(axis, probability, threshold):
"""given x and P(x), find P(x>t)"""
- exceeding = np.where(axis > threshold)
- integral = integrate.simps(probability[exceeding], axis[exceeding])
- return integral / integrate.simps(probability, axis)
+ try:
+ exceeding = np.where(axis > threshold)
+ integral = integrate.simps(probability[exceeding], axis[exceeding])
+ retval = integral / integrate.simps(probability, axis)
+ except FloatingPointError:
+ print "prob_exceed: invalid PDF"
+ retval = -1.
+
+ return retval
# TODO: remove trailing space
Please sign in to comment.
Something went wrong with that request. Please try again.