Skip to content
Snippets Groups Projects
Commit 55d0901a authored by Branson Stephens's avatar Branson Stephens
Browse files

Added some features to gstlalcbc_report (semi-adaptive mass bins, default...

Added some features to gstlalcbc_report (semi-adaptive mass bins, default query, end_time template filter). Also refactored gstlalcbc_report into reports.py. Some other cleaning up.
parent 0f11c0d8
No related branches found
No related tags found
No related merge requests found
from django.http import HttpResponse
from django.http import HttpResponseForbidden
from django.template import RequestContext
from django.shortcuts import render_to_response
from django.conf import settings
......@@ -7,7 +8,26 @@ from django.conf import settings
from gracedb.models import Event
from django.db.models import Q
import os, datetime, json, time
import os, json
from django.core.urlresolvers import reverse
from models import CoincInspiralEvent ,SingleInspiral
from forms import SimpleSearchFormWithSubclasses
from query import parseQuery
from django.db.models import Max, Min
import matplotlib
matplotlib.use('Agg')
import numpy
import matplotlib.pyplot as plot
import StringIO
import base64
import sys
import time
from datetime import datetime, timedelta
from utils import posixToGpsTime
def histo(request):
......@@ -46,8 +66,8 @@ def histo(request):
def rate_data(request):
# XXX there is a better way -- should be using group_by or something.
# WAAY too many queries (~300) going on here.
now = datetime.datetime.now()
day = datetime.timedelta(1)
now = datetime.now()
day = timedelta(1)
ts_min = now - 60 * day
ts_max = now
......@@ -84,3 +104,193 @@ def rate_data(request):
return series
# XXX This should be configurable / moddable or something
MAX_QUERY_RESULTS = 1000
# The following two util routines are for gstlalcbc_report. This is messy.
def cluster(events):
# FIXME N^2 clustering, but event list should always be small anyway...
def quieter(e1, events = events, win = 5):
for e2 in events:
if e2.far == e1.far:
# XXX I need subclass attributes here.
if (e2.gpstime < e1.gpstime + win) and (e2.gpstime > e1.gpstime - win) and (e2.snr > e1.snr):
return True
else:
if (e2.gpstime < e1.gpstime + win) and (e2.gpstime > e1.gpstime - win) and (e2.far < e1.far):
return True
return False
return [e for e in events if not quieter(e)]
def to_png_image(out = sys.stdout):
f = StringIO.StringIO()
plot.savefig(f, format="png")
return base64.b64encode(f.getvalue())
def gstlalcbc_report(request, format=""):
if not request.user or not request.user.is_authenticated():
return HttpResponseForbidden("Forbidden")
if request.method == "GET":
if "query" not in request.GET:
# Use default query. LowMass events from the past week.
t_high = datetime.now()
dt = timedelta(days=7)
t_low = t_high - dt
t_high = posixToGpsTime(time.mktime(t_high.timetuple()))
t_low = posixToGpsTime(time.mktime(t_low.timetuple()))
query = 'CBC LowMass %d .. %d' % (t_low, t_high)
rawquery = query
form = SimpleSearchFormWithSubclasses({'query': query})
else:
form = SimpleSearchFormWithSubclasses(request.GET)
rawquery = request.GET['query']
else:
form = SimpleSearchFormWithSubclasses(request.POST)
rawquery = request.POST['query']
if form.is_valid():
objects = form.cleaned_data['query']
# Check for foreign objects.
for obj in objects:
if not isinstance(obj, CoincInspiralEvent):
errormsg = 'Your query returned items that are not CoincInspiral Events. '
errormsg += 'Please try again.'
form = SimpleSearchFormWithSubclasses()
return render_to_response('gracedb/gstlalcbc_report.html',
{ 'form':form, 'message':errormsg},
context_instance=RequestContext(request))
# Check that we have a well-defined GPS time range.
# Find the gpstime limits of the query by diving into the query object.
# XXX Down the rabbit hole!
qthing = parseQuery(rawquery)
gpsrange = None
for child in qthing.children:
if 'gpstime' in str(child):
for subchild in child.children:
if isinstance(subchild, tuple):
gpsrange = subchild[1]
if not gpsrange:
# Bounce back to the user with an error message
errormsg = 'Your query does not have a gpstime range. Please try again.'
form = SimpleSearchFormWithSubclasses()
return render_to_response('gracedb/gstlalcbc_report.html',
{ 'form':form, 'message':errormsg},
context_instance=RequestContext(request))
lt = int(gpsrange[1]) - int(gpsrange[0])
# Check that there aren't too many objects.
# XXX Hardcoded limit
if objects.count() > 2000:
errormsg = 'Your query returned too many events. Please try again.'
return render_to_response('gracedb/gstlalcbc_report.html',
{ 'form':form, 'message':errormsg},
context_instance=RequestContext(request))
# Zero events will break the min/max thing that comes next.
if objects.count() < 1:
errormsg = 'Your query returned no events. Please try again.'
return render_to_response('gracedb/gstlalcbc_report.html',
{ 'form':form, 'message':errormsg},
context_instance=RequestContext(request))
# Find the min and max on the set of objects.
#gpstime_limits = objects.aggregate(Max('gpstime'), Min('gpstime'))
mchirp_limits = objects.aggregate(Max('coincinspiralevent__mchirp'),
Min('coincinspiralevent__mchirp'))
mass_limits = objects.aggregate(Max('coincinspiralevent__mass'),
Min('coincinspiralevent__mass'))
clustered_events = cluster(objects)
clustered_events = sorted(clustered_events, None, key=lambda x: x.far)
# Make IFAR plot.
ifars = numpy.array(sorted([1.0 / e.far for e in clustered_events])[::-1])
N = numpy.linspace(1, len(ifars), len(ifars))
eN = numpy.linspace(1, 1000 * len(ifars), 1000 * len(ifars)) / 1000.
expected_ifars = lt / eN
up = eN + eN**.5
down = eN - eN**.5
down[down < 0.9] = 0.9
plot.figure(figsize=(6,5))
plot.loglog(ifars[::-1], N[::-1])
plot.fill_between(expected_ifars[::-1], down[::-1], up[::-1], alpha=0.1)
plot.loglog(expected_ifars[::-1], eN[::-1])
plot.ylim([0.9, len(ifars)])
plot.xlabel('IFAR (s)')
plot.ylabel('N')
plot.grid()
ifar_plot = to_png_image()
# Set the color map for loudest event table. Depends on lt.
#FAR_color_map = [ { 'max_far' : 1.e-12, 'color' : 'gold', 'desc' : '< 1/10000 yrs'},
# { 'max_far' : 3.e-10, 'color' : 'silver', 'desc' : '< 1/100 yrs'},
# { 'max_far' : 3.e-8, 'color' : '#A67D3D', 'desc' : '< 1/yr'},
# { 'max_far' : 1.0/lt, 'color' : '#B2C248', 'desc' : 'louder than expected'}, ]
# XXX Okay, this sucks. There is no switch/case structure in the django template
# language. So I couldn't think of any way to do this without checking ranges.
# And the range is zero to infinity. So here I go...
FAR_color_map = [ { 'min_far' : 0.0, 'max_far' : 1.e-12, 'color' : 'gold', 'desc' : '< 1/10000 yrs'},
{ 'min_far' : 1.e-12, 'max_far' : 3.e-10, 'color' : 'silver', 'desc' : '< 1/100 yrs'},
{ 'min_far' : 3.e-10, 'max_far' : 3.e-8, 'color' : '#A67D3D', 'desc' : '< 1/yr'},
{ 'min_far' : 3.e-8, 'max_far' : 1.0/lt, 'color' : '#B2C248', 'desc' : 'louder than expected'},
{ 'min_far' : 1.0/lt, 'max_far' : float('inf'), 'color' : 'white', 'desc' : ''}, ]
# Make mass distribution plots
mchirp = numpy.array([e.mchirp for e in clustered_events])
plot.figure(figsize=(6,4))
lower = int(mchirp_limits['coincinspiralevent__mchirp__min'])
upper = int(mchirp_limits['coincinspiralevent__mchirp__max']) + 1
# How to decide the number of bins?
N_bins = 20
delta = max(float(upper-lower)/N_bins,0.2)
plot.hist(mchirp, bins=numpy.arange(lower,upper,delta))
plot.xlabel('Chirp Mass')
plot.ylabel('count')
plot.grid()
mchirp_dist = to_png_image()
mass = numpy.array([e.mass for e in clustered_events])
plot.figure(figsize=(6,4))
lower = int(mass_limits['coincinspiralevent__mass__min'])
upper = int(mass_limits['coincinspiralevent__mass__max']) + 1
N_bins = 20
delta = max(float(upper-lower)/N_bins,0.2)
plot.hist(mass, bins=numpy.arange(lower,upper,delta))
plot.xlabel('Total Mass')
plot.ylabel('count')
plot.grid()
mass_dist = to_png_image()
if objects.count() == 1:
title = "Query returned %s event." % objects.count()
else:
title = "Query returned %s events." % objects.count()
context = {
'title': title,
'form': form,
'formAction': reverse(gstlalcbc_report),
'count' : objects.count(),
'rawquery' : rawquery,
'FAR_color_map' : FAR_color_map,
'mass_dist' : mass_dist,
'mchirp_dist' : mchirp_dist,
'ifar_plot' : ifar_plot,
'clustered_events' : clustered_events,
}
return render_to_response('gracedb/gstlalcbc_report.html', context,
context_instance=RequestContext(request))
return render_to_response('gracedb/gstlalcbc_report.html',
{ 'form' : form,
},
context_instance=RequestContext(request))
......@@ -168,8 +168,9 @@ def end_time(event,digits=4):
# in the django template, so that the final digits get truncated and padded
# with zeros. Why? The current gpstime is only 10 digits. So I'm going to have to
# make this into a string. Totally sick, I know.
decimal_part = float(event.end_time_ns)/1.e9
decimal_part = round(decimal_part,digits)
return str(event.end_time) + str(decimal_part)[1:]
decimal_part = round(float(event.end_time_ns)/1.e9,digits)
# ugh. must pad with zeros to the right.
decimal_part = str(decimal_part)[1:].ljust(digits+1,'0')
return str(event.end_time) + decimal_part
except Exception, e:
return None
......@@ -9,7 +9,6 @@ urlpatterns = patterns('gracedb.views',
url (r'^$', 'index', name="home"),
url (r'^create/$', 'create', name="create"),
url (r'^search/(?P<format>(json|flex))?$', 'search', name="search"),
url (r'^gstlalcbc_report/(?P<format>(json|flex))?$', 'gstlalcbc_report', name="gstlalcbc_report"),
url (r'^view/(?P<graceid>[GEHT]\d+)', 'view', name="view"),
url (r'^voevent/(?P<graceid>[GEHT]\d+)', 'voevent', name="voevent"),
url (r'^skyalert/(?P<graceid>[GEHT]\d+)', 'skyalert', name="skyalert"),
......
......@@ -38,16 +38,6 @@ from templatetags.scientific import scientific
from buildVOEvent import buildVOEvent, submitToSkyalert
from django.db.models import Max, Min
import matplotlib
matplotlib.use('Agg')
import numpy
import matplotlib.pyplot as plot
import StringIO
import base64
import sys
import logging
# XXX This should be configurable / moddable or something
MAX_QUERY_RESULTS = 1000
......@@ -808,173 +798,6 @@ def search(request, format=""):
},
context_instance=RequestContext(request))
# The following two util routines are for gstlalcbc_report. This is messy.
def cluster(events):
# FIXME N^2 clustering, but event list should always be small anyway...
def quieter(e1, events = events, win = 5):
for e2 in events:
if e2.far == e1.far:
# XXX I need subclass attributes here.
if (e2.gpstime < e1.gpstime + win) and (e2.gpstime > e1.gpstime - win) and (e2.snr > e1.snr):
return True
else:
if (e2.gpstime < e1.gpstime + win) and (e2.gpstime > e1.gpstime - win) and (e2.far < e1.far):
return True
return False
return [e for e in events if not quieter(e)]
def to_png_image(out = sys.stdout):
f = StringIO.StringIO()
plot.savefig(f, format="png")
return base64.b64encode(f.getvalue())
def gstlalcbc_report(request, format=""):
if not request.user or not request.user.is_authenticated():
return HttpResponseForbidden("Forbidden")
logger = logging.getLogger(__name__)
if request.method == "GET" and "query" not in request.GET:
form = SimpleSearchFormWithSubclasses()
else:
if request.method == "GET":
form = SimpleSearchFormWithSubclasses(request.GET)
rawquery = request.GET['query']
else:
form = SimpleSearchFormWithSubclasses(request.POST)
rawquery = request.POST['query']
if form.is_valid():
objects = form.cleaned_data['query']
# Check for foreign objects.
for obj in objects:
logger.debug("object is a %s" % obj.__class__.__name__)
if not isinstance(obj, CoincInspiralEvent):
errormsg = 'Your query returned item(s) that are not CoincInspiral Events. '
errormsg += 'Please try again.'
form = SimpleSearchFormWithSubclasses()
return render_to_response('gracedb/gstlalcbc_report.html',
{ 'form':form, 'message':errormsg},
context_instance=RequestContext(request))
# Check that we have a well-defined GPS time range.
# Find the gpstime limits of the query by diving into the query object.
# XXX Down the rabbit hole!
qthing = parseQuery(rawquery)
gpsrange = None
for child in qthing.children:
if 'gpstime' in str(child):
for subchild in child.children:
if isinstance(subchild, tuple):
gpsrange = subchild[1]
if not gpsrange:
# Bounce back to the user with an error message
errormsg = 'Your query does not have a gpstime range. Please try again.'
form = SimpleSearchFormWithSubclasses()
return render_to_response('gracedb/gstlalcbc_report.html',
{ 'form':form, 'message':errormsg},
context_instance=RequestContext(request))
lt = int(gpsrange[1]) - int(gpsrange[0])
# Check that there aren't too many objects.
# XXX Hardcoded limit
if objects.count() > 2000:
errormsg = 'Your query returned too many events. Please try again.'
return render_to_response('gracedb/gstlalcbc_report.html',
{ 'form':form, 'message':errormsg},
context_instance=RequestContext(request))
# Find the min and max on the set of objects.
#gpstime_limits = objects.aggregate(Max('gpstime'), Min('gpstime'))
mchirp_limits = objects.aggregate(Max('coincinspiralevent__mchirp'),
Min('coincinspiralevent__mchirp'))
mass_limits = objects.aggregate(Max('coincinspiralevent__mass'),
Min('coincinspiralevent__mass'))
clustered_events = cluster(objects)
clustered_events = sorted(clustered_events, None, key=lambda x: x.far)
# Make IFAR plot.
ifars = numpy.array(sorted([1.0 / e.far for e in clustered_events])[::-1])
N = numpy.linspace(1, len(ifars), len(ifars))
eN = numpy.linspace(1, 1000 * len(ifars), 1000 * len(ifars)) / 1000.
expected_ifars = lt / eN
up = eN + eN**.5
down = eN - eN**.5
down[down < 0.9] = 0.9
plot.figure(figsize=(6,5))
plot.loglog(ifars[::-1], N[::-1])
plot.fill_between(expected_ifars[::-1], down[::-1], up[::-1], alpha=0.1)
plot.loglog(expected_ifars[::-1], eN[::-1])
plot.ylim([0.9, len(ifars)])
plot.xlabel('IFAR (s)')
plot.ylabel('N')
plot.grid()
ifar_plot = to_png_image()
# Set the color map for loudest event table. Depends on lt.
#FAR_color_map = [ { 'max_far' : 1.e-12, 'color' : 'gold', 'desc' : '< 1/10000 yrs'},
# { 'max_far' : 3.e-10, 'color' : 'silver', 'desc' : '< 1/100 yrs'},
# { 'max_far' : 3.e-8, 'color' : '#A67D3D', 'desc' : '< 1/yr'},
# { 'max_far' : 1.0/lt, 'color' : '#B2C248', 'desc' : 'louder than expected'}, ]
# XXX Okay, this sucks. There is no switch/case structure in the django template
# language. So I couldn't think of any way to do this without checking ranges.
# And the range is zero to infinity. So here I go...
FAR_color_map = [ { 'min_far' : 0.0, 'max_far' : 1.e-12, 'color' : 'gold', 'desc' : '< 1/10000 yrs'},
{ 'min_far' : 1.e-12, 'max_far' : 3.e-10, 'color' : 'silver', 'desc' : '< 1/100 yrs'},
{ 'min_far' : 3.e-10, 'max_far' : 3.e-8, 'color' : '#A67D3D', 'desc' : '< 1/yr'},
{ 'min_far' : 3.e-8, 'max_far' : 1.0/lt, 'color' : '#B2C248', 'desc' : 'louder than expected'},
{ 'min_far' : 1.0/lt, 'max_far' : float('inf'), 'color' : 'white', 'desc' : ''}, ]
# Make mass distribution plots
mchirp = numpy.array([e.mchirp for e in clustered_events])
plot.figure(figsize=(6,4))
lower = int(mchirp_limits['coincinspiralevent__mchirp__min'])
upper = int(mchirp_limits['coincinspiralevent__mchirp__max']) + 1
plot.hist(mchirp, bins=numpy.arange(lower,upper,0.25))
plot.xlabel('Chirp Mass')
plot.ylabel('count')
plot.grid()
mchirp_dist = to_png_image()
mass = numpy.array([e.mass for e in clustered_events])
plot.figure(figsize=(6,4))
lower = int(mass_limits['coincinspiralevent__mass__min'])
upper = int(mass_limits['coincinspiralevent__mass__max']) + 1
plot.hist(mass, bins=numpy.arange(lower,upper,5/8.0))
plot.xlabel('Total Mass')
plot.ylabel('count')
plot.grid()
mass_dist = to_png_image()
if objects.count() == 1:
title = "Query returned %s event." % objects.count()
else:
title = "Query returned %s events." % objects.count()
context = {
'title': title,
'form': form,
'formAction': reverse(gstlalcbc_report),
'count' : objects.count(),
'rawquery' : rawquery,
'FAR_color_map' : FAR_color_map,
'mass_dist' : mass_dist,
'mchirp_dist' : mchirp_dist,
'ifar_plot' : ifar_plot,
'clustered_events' : clustered_events,
}
return render_to_response('gracedb/gstlalcbc_report.html', context,
context_instance=RequestContext(request))
return render_to_response('gracedb/gstlalcbc_report.html',
{ 'form' : form,
},
context_instance=RequestContext(request))
def oldsearch(request):
assert request.user
if request.method == 'GET':
......
......@@ -5,6 +5,9 @@
{# Analysis-specific attributes for an LM event#}
{% block analysis_specific %}
{# Test whether the object has the end_time property. Older (non CoincInspiral) events don't. #}
{% if object.end_time %}
<div id="container" style="display:table;width:100%">
<div style="display:table-row;width:100%">
......@@ -364,4 +367,6 @@
{% endif %}
{% endif %} <!-- object has end_time property. -->
{% endblock %}
......@@ -83,7 +83,6 @@ onload="document.search_form.query.focus();"
<th>Rank</th>
<th>GraceID</th>
<th>gpstime</th>
<th>end_time_ns</th>
<th>FAR</th>
<th>IFOs</th>
<th>Total Mass</th>
......@@ -98,8 +97,7 @@ onload="document.search_form.query.focus();"
{% endfor %}
<td>{{ forloop.counter }} </td>
<td><a href="{% url view obj.graceid %}">{{ obj.graceid }}</a></td>
<td> {{ obj.gpstime }} </td>
<td> {{ obj.end_time_ns }} </td>
<td> {{ obj|end_time }} </td>
<td> {{ obj.far|scientific }} </td>
<td> {{ obj.instruments }} </td>
<td> {{ obj.mass|floatformat:4 }} </td>
......
......@@ -82,6 +82,10 @@ function toggle(id) {
{% block content %}
<br/>
<a href="{% url gstlalcbc_report %}"><h3>Dynamic CBC Report</h3></a>
<br/>
<br/>
<a name="latency" href="javascript:toggle('latency');"><h3>Latency</h3></a>
<div id="latency" style="display:none;">
......
......@@ -30,6 +30,7 @@ urlpatterns = patterns('',
url (r'^feeds/$', feedview, name="feeds"),
url (r'^reports/$', 'gracedb.reports.histo', name="reports"),
url (r'^reports/gstlalcbc_report/(?P<format>(json|flex))?$', 'gracedb.reports.gstlalcbc_report', name="gstlalcbc_report"),
(r'^reports/(?P<path>.+)$', 'django.views.static.serve',
{'document_root': settings.LATENCY_REPORT_DEST_DIR}),
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment