gstlal_llsnrhistory

#!/usr/bin/env python
#
# Copyright (C) 2011  Chad Hanna
#
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 2 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

import sys
import cgi
import cgitb
import os
os.environ["MPLCONFIGDIR"] = "/tmp"
import matplotlib
matplotlib.use('Agg')
import numpy
import matplotlib.pyplot as plt
import time
import StringIO
import base64
from urlparse import urlparse
cgitb.enable()
form = cgi.FieldStorage()
import pylab
from pylal import rate
import lal


def now():
    return float(lal.UTCToGPS(time.gmtime()))

def to_png_image():
    f = StringIO.StringIO()
    plt.savefig(f, format="png")
    sys.stdout.write( "Content-type: image/png\r\n\r\n" + f.getvalue() )
    f.close()

def read_registry(dir, dataurl, ids):
    nodedict = {}
    for id in ids:
        url = '%s/%s%s' % (dir, id, dataurl)
        try:
            tmp = open(url,"r")
            nodedict[id] = urlparse(tmp.readline()).netloc
            tmp.close()
        except IOError:
            nodedict[id] = ""
    return nodedict
    
def setup_plot():
    fig = plt.figure(figsize=(20,5),)
    fig.patch.set_alpha(0.0)
    h = fig.add_subplot(111, axisbg = 'k')
    plt.subplots_adjust(left = .062, right = 0.98, bottom = 0.3)
    return fig, h

def finish_plot(ids, registry, ylim, title=''):
    plt.grid(color=(0.1,0.4,0.5), linewidth=2)  
    ticks = ["%s : %s " % (id, registry[id]) for id in ids]
    plt.xticks(numpy.arange(len(ids))+.3, ticks, rotation=90, fontsize = 10)
    plt.xlim([0, len(ids)])
    plt.ylim(ylim)
    tickpoints = numpy.linspace(ylim[0], ylim[1], 8)
    ticks = ["%.1e" % (10.**t,) for t in tickpoints]
    plt.yticks(tickpoints, ticks, fontsize = 14)
    plt.title(title, fontsize = 18)
    to_png_image()


def get_ids(form):
    idrange = [int(n) for n in form.getvalue("id").split(",")]
    #FIXME relies on 4 digit ids
    ids = ['%04d' % (job,) for job in range(idrange[0], idrange[1]+1)]
    return ids

if "dir" not in form:
    raise ValueError("must specify dir")
if "id" not in form:
    raise ValueError("must specify id")

ids = get_ids(form)
directory = form.getvalue("dir")
start = form.getvalue("start")
stop = form.getvalue("stop")
reg = read_registry(form.getvalue("dir"), "_registry.txt", ids)

time_now = now()
yaxis = []

for i, id in enumerate(ids):
    try:
        fname = "%s/%s/%s.txt" % (directory, id, "snr_history")
        yaxis.append(float(open("%s/%s/%s.txt" % (directory, id, "bank")).readline().split()[2]))
            
    except IOError as e:
        print>>sys.stderr, "couldn't open ", fname

yaxis = numpy.array(sorted(set(yaxis)))


if start is not None and stop is not None:
    ba = rate.BinnedArray(rate.NDBins((rate.LinearBins(float(start) - time_now, float(stop) - time_now, 900), rate.IrregularBins(yaxis))))
else:
    ba = rate.BinnedArray(rate.NDBins((rate.LinearBins(-4500,-900,900), rate.IrregularBins(yaxis))))

print >>sys.stderr, min(yaxis), max(yaxis)

window = rate.gaussian_window(10, 1, sigma = 3)
print >> sys.stderr, ba.array.shape
ba.array[:] = 1.

for i, id in enumerate(ids):
    try:
        fname = "%s/%s/%s.txt" % (directory, id, "snr_history")
        dur = float(open("%s/%s/%s.txt" % (directory, id, "bank")).readline().split()[2])
        for line in open(fname):
            thist, thissnr = line.split()
            t = float(thist)-time_now
            try:
                if ba[t,dur] < thissnr:
                    ba[t,dur] = thissnr
            except IndexError:
                pass
            
    except IOError as e:
        print>>sys.stderr, "couldn't open ", fname

matplotlib.rcParams.update({"text.usetex": True})

rate.filter_array(ba.array, window)
fig = pylab.figure(figsize=(20,4.0))
pylab.title('max $\\rho_c$ in coincidence')
pylab.xlabel("time since %d (s)" % (int(time_now),))
pylab.ylabel("$\mathcal{M} (M_\odot)$")
pylab.xlim([ba.centres()[0][0], ba.centres()[0][-1]])
pylab.ylim([ba.centres()[1][0], ba.centres()[1][-1]])
pylab.pcolormesh(ba.centres()[0], ba.centres()[1], ba.array.T, norm = matplotlib.colors.LogNorm())
pylab.colorbar()
to_png_image()