First commit of lal_skymap

configure.ac

@@ -17,6 +17,7 @@ AC_CONFIG_FILES([ \
 	src/bin/Makefile \
 	src/lib/gstlal.pc \
 	src/lib/Makefile \
+	src/skymap/Makefile \
 	src/plugins/Makefile \
 	src/plugins/python/Makefile \
 	src/python/Makefile \

src/Makefile.am

-SUBDIRS = lib plugins python bin utilities
+SUBDIRS = lib skymap plugins python bin utilities

src/plugins/Makefile.am

 SUBDIRS = python
-INCLUDES = -I$(top_srcdir)/src/lib
+INCLUDES = -I$(top_srcdir)/src/lib -I$(top_srcdir)/src/skymap
 
 plugin_LTLIBRARIES = libgstlal.la
 
@@ -20,6 +20,7 @@ libgstlal_la_SOURCES = \
 	gstlalcollectpads.h gstlalcollectpads.c \
 	gstadder.h gstadder.c \
 	gstlal_coinc.h gstlal_coinc.c \
+	gstlal_skymap.h gstlal_skymap.c \
 	gstlal_triggergen.h gstlal_triggergen.c \
 	gstlal_triggerxmlwriter.h gstlal_triggerxmlwriter.c \
 	gstlal_gate.h gstlal_gate.c \
@@ -34,7 +35,7 @@ libgstlal_la_SOURCES = \
 	gstlal_iirbank.h gstlal_iirbank.c \
 	$(NDSSRC_SOURCES)
 libgstlal_la_CFLAGS = -O2 $(AM_CFLAGS) $(CFLAGS) ${gstreamer_CFLAGS} ${LAL_CFLAGS} ${GSL_CFLAGS} ${FFTW_CFLAGS} ${NDS_CFLAGS}
-libgstlal_la_LIBADD = ${gstreamer_LIBS} ${LAL_LIBS} ${GSL_LIBS} ${FFTW_LIBS} ${NDS_LIBS} $(top_srcdir)/src/lib/libgstlal.la
+libgstlal_la_LIBADD = ${gstreamer_LIBS} ${LAL_LIBS} ${GSL_LIBS} ${FFTW_LIBS} ${NDS_LIBS} $(top_srcdir)/src/skymap/libwanalysis.la $(top_srcdir)/src/lib/libgstlal.la
 libgstlal_la_LDFLAGS = ${GSTREAMER_PLUGIN_LDFLAGS}
 
 noinst_HEADERS = gstlal_plugins.h

src/plugins/gstlal.c

@@ -82,6 +82,7 @@
 #include <gstlal_nxydump.h>
 #include <gstadder.h>
 #include <gstlal_coinc.h>
+#include <gstlal_skymap.h>
 #include <gstlal_triggergen.h>
 #include <gstlal_triggerxmlwriter.h>
 #include <gstlal_gate.h>
@@ -133,6 +134,7 @@ static gboolean plugin_init(GstPlugin *plugin)
 		{"lal_nxydump", GSTLAL_NXYDUMP_TYPE},
 		{"lal_adder", GST_TYPE_ADDER},
 		{"lal_coinc", GSTLAL_COINC_TYPE},
+		{"lal_skymap", GSTLAL_SKYMAP_TYPE},
 		{"lal_triggergen", GSTLAL_TRIGGERGEN_TYPE},
 		{"lal_triggerxmlwriter", GSTLAL_TRIGGERXMLWRITER_TYPE},
 		{"lal_gate", GSTLAL_GATE_TYPE},

src/plugins/gstlal_skymap.h

+/*
+ * Copyright (C) 2009 Leo Singer <leo.singer@ligo.org>
+ * 
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject to
+ * the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+ * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+ * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+ * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ *
+ * You should have received a copy of the GNU Library General Public
+ * License along with this library; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
+ * USA.
+ */
+
+
+/*
+ * ============================================================================
+ *
+ *                                  Preamble
+ *
+ * ============================================================================
+ */
+
+
+#ifndef __GSTLAL_SKYMAP_H__
+#define __GSTLAL_SKYMAP_H__
+
+
+#include <gst/gst.h>
+#include <gst/base/gstcollectpads.h>
+#include <lal/LIGOMetadataTables.h>
+#include "wanalysis.h"
+
+
+G_BEGIN_DECLS
+
+
+/*
+ * ============================================================================
+ *
+ *                            Coincidence Generator
+ *
+ * ============================================================================
+ */
+
+
+#define GSTLAL_SKYMAP_TYPE \
+	(gstlal_skymap_get_type())
+#define GSTLAL_SKYMAP(obj) \
+	(G_TYPE_CHECK_INSTANCE_CAST((obj), GSTLAL_SKYMAP_TYPE, GSTLALSkymap))
+#define GSTLAL_SKYMAP_CLASS(klass) \
+	(G_TYPE_CHECK_CLASS_CAST((klass), GSTLAL_SKYMAP_TYPE, GSTLALSkymapClass))
+#define GST_IS_GSTLAL_SKYMAP(obj) \
+	(G_TYPE_CHECK_INSTANCE_TYPE((obj), GSTLAL_SKYMAP_TYPE))
+#define GST_IS_GSTLAL_SKYMAP_CLASS(klass) \
+	(G_TYPE_CHECK_CLASS_TYPE((klass), GSTLAL_SKYMAP_TYPE))
+
+
+typedef struct {
+	GstElementClass parent_class;
+} GSTLALSkymapClass;
+
+
+typedef struct {
+	GstElement element;
+
+	/* collectpads */
+	GstCollectPads *collect;
+	GstPadEventFunction collect_event;
+
+	/* pads & collectdatas */
+	GstPad* srcpad;
+	GstCollectData* coinc_collectdata;
+	GSList* snr_collectdatas;
+
+	/* timing */
+	guint64 dt;
+	guint64 trigger_present_padding;
+	guint64 trigger_absent_padding;
+
+	/* template bank */
+	GMutex *bank_lock;
+	char *bank_filename;
+	SnglInspiralTable *bank;
+	int ntemplates;
+
+	/* analysis */
+	analysis wanalysis;
+} GSTLALSkymap;
+
+
+GType gstlal_skymap_get_type(void);
+
+
+G_END_DECLS
+
+
+#endif	/* __GSTLAL_SKYMAP_H__ */

src/skymap/Makefile.am

+lib_LTLIBRARIES = libwanalysis.la
+
+libwanalysis_la_SOURCES = wanalysis.c wanalysis.h
+libwanalysis_la_CPPFLAGS = $(LAL_CFLAGS)
+libwanalysis_la_LIBADD = $(LAL_LIBS)
+libwanalysis_la_LDFLAGS = $(LAL_LDFLAGS)

src/skymap/README

+This directory contains files taken from
+
+omega/branches/gstlal/src/bayesian
+
+that specify a common interface through which both Omega and GSTLAL can
+produce a full skymap.  This functionality will eventually be migrated to the
+XLALSkymap package, and this directory deleted.

src/skymap/wanalysis.c

+// C headers
+
+#include <math.h>
+#include <string.h>
+
+// LAL headers
+
+#include <lal/LALConstants.h>
+#include <lal/Skymap.h>
+#include <lal/Random.h>
+
+// Omega headers
+
+#include "wanalysis.h"
+
+void analysis_default_construct(analysis* a)
+{
+    a->n_detectors = 0;
+    a->rate = 0;
+    a->wSw = 0;
+    a->n_directions = 0;
+    a->directions = 0;
+    a->min_t = log(0);
+    a->max_t = -log(0);
+    a->delta_t = 1;
+    a->log_skymap = 0;
+    a->p_realloc = realloc;
+    a->p_free = free;
+    for (size_t i = 0; i != XLALSKYMAP_N; ++i)
+    {
+        a->xSw_real[i] = 0;
+	a->xSw_imag[i] = 0;
+        a->min_ts[i] = 0;
+	a->max_ts[i] = 0;
+        a->calibration_error[i] = 0.1;
+    }
+}
+
+void analysis_default_directions(analysis* a)
+{
+    size_t u = 360 * 10 / 4;
+    size_t v = 180 * 10 / 4;
+
+    a->n_directions = u * v;
+    a->directions = a->p_realloc(0, a->n_directions * 2 * sizeof(double));
+    
+    double* p = a->directions;
+    for (size_t j = 0; j != v; ++j)
+    {
+        double theta = (0.4 * (j + .5)) * LAL_PI / 180; 
+        for (size_t i = 0; i != u; ++i)
+	{
+	    double phi = (0.4 * (i + .5)) * LAL_PI / 180;
+	    *(p++) = theta;
+	    *(p++) = phi;
+	}
+    }
+}
+
+int analysis_identify_detector(const char* s)
+{
+    // truncate the string to first two letters
+    char c[3];
+    c[0] = s[0];
+    c[1] = s[1];
+    c[2] =   0;
+    // compare with detector codes
+    // TODO: we can probably replace this with iteration through
+    //       LALCachedDetectors via LALDetector::frDetector.prefix
+    if (strcmp(c, "T1") == 0) 
+    {
+    	return LAL_TAMA_300_DETECTOR;
+    }
+    else if (strcmp(c, "V1") == 0) 
+    {
+        return LAL_VIRGO_DETECTOR;
+    }
+    else if (strcmp(c, "G1") == 0) 
+    {
+        return LAL_GEO_600_DETECTOR;
+    }
+    else if (strcmp(c, "H2") == 0) 
+    {
+        return LAL_LHO_2K_DETECTOR;
+    }
+    else if (strcmp(c, "H1") == 0) 
+    {
+        return LAL_LHO_4K_DETECTOR;
+    }
+    else if (strcmp(c, "L1") == 0) 
+    {
+         return LAL_LLO_4K_DETECTOR;
+    }
+    else 
+    {
+        return -1;
+    }
+}
+
+analysis* analysis_example(void)
+{
+    analysis* a = (analysis*) malloc(sizeof(analysis));
+    
+    analysis_default_construct(a);
+    analysis_default_directions(a);
+
+    double duration = 1; // seconds of data available
+    RandomParams* rng = XLALCreateRandomParams(0);
+
+    a->n_detectors = 3;
+
+    a->detectors[0] = LAL_LHO_4K_DETECTOR;
+    a->detectors[1] = LAL_LLO_4K_DETECTOR;
+    a->detectors[2] = LAL_VIRGO_DETECTOR;
+    
+    a->rate = 512; // Hz
+    
+    a->wSw = (double*) a->p_realloc(0, sizeof(double) * a->n_detectors);
+
+    a->min_t = 0.4; // seconds at earth barycenter
+    a->max_t = 0.6;
+
+    a->log_skymap = (double*) a->p_realloc(0, sizeof(double) * a->n_directions);
+
+    for (size_t i = 0; i != a->n_detectors; ++i)
+    {
+        a->wSw[i] = 1.0; // white noise and unit template
+	
+	a->xSw_real[i] = (double*) a->p_realloc(0, sizeof(double) * a->rate * duration);
+        a->xSw_imag[i] = (double*) a->p_realloc(0, sizeof(double) * a->rate * duration);
+        
+	for (size_t j = 0; j != (a->rate * duration); ++j)
+	{
+	    // make gaussian noise
+	    a->xSw_real[i][j] = XLALUniformDeviate(rng) * sqrt(a->wSw[i]);
+	    a->xSw_imag[i][j] = XLALUniformDeviate(rng) * sqrt(a->wSw[i]);
+	}
+
+        // allowed time ranges at each detector
+	a->min_ts[i] = 0.49; // seconds from reference time
+	a->max_ts[i] = 0.51;
+
+    }
+
+    analyze(a);
+
+    // struct now contains pointers to all the input data used and the results
+
+    return a;
+}
+
+void analyze(analysis* s)
+{
+
+    XLALSkymapPlanType plan;
+
+    // Construct the network
+    XLALSkymapPlanConstruct(s->rate, s->n_detectors, s->detectors, &plan);
+
+    // Working array
+    double* p_t = 0;
+
+    for (size_t i = 0; i != s->n_directions; ++i) {
+        XLALSkymapDirectionPropertiesType properties;
+        XLALSkymapDirectionPropertiesConstruct(&plan, s->directions + (i * 2), &properties);
+
+        double t_begin = s->min_t;
+        double t_end = s->max_t;
+
+        for (size_t j = 0; j != s->n_detectors; ++j) {
+            t_begin = fmax(t_begin, s->min_ts[j] - properties.delay[j]);
+            t_end = fmin(t_end, s->max_ts[j] - properties.delay[j]);
+        }
+
+        if (t_begin < t_end) {
+
+            t_begin -= 0.0001;
+            t_end   += 0.0001;
+
+            XLALSkymapKernelType kernel;
+            // strict version assumes no calibration error (compatible with LAL release)
+            //XLALSkymapKernelConstruct(&plan, &properties, wSw, &kernel);
+            // loose version allows for amplitude calibration error (requires LAL built from repository head as of 4/18/10)
+            XLALSkymapUncertainKernelConstruct(&plan, &properties, s->wSw, s->calibration_error, &kernel);
+
+            double old_p;
+            double new_p = log(0);
+
+            // take at least 10 samples
+            double dt = fmin(s->delta_t, (t_end - t_begin) * 0.1);
+
+            do {
+                size_t p_t_size = (size_t) (ceil((t_end - t_begin) / dt));
+                p_t = (double*) s->p_realloc(p_t, p_t_size * sizeof(double));
+                double t = t_begin;
+                for (size_t j = 0; j != p_t_size; ++j)
+                {
+                      double log_p_real, log_p_imag;
+                      XLALSkymapApply(&plan, &properties, &kernel, s->xSw_real, t, &log_p_real);
+                      XLALSkymapApply(&plan, &properties, &kernel, s->xSw_imag, t, &log_p_imag);
+                      p_t[j] = log_p_real + log_p_imag;
+                      t += dt;
+                }
+                old_p = new_p;
+                new_p = XLALSkymapLogTotalExp(&p_t[0], &p_t[0] + p_t_size) - log(p_t_size);
+                dt = dt / 2;
+            } while (fabs(new_p - old_p) > .1);
+            s->log_skymap[i] = new_p + log(t_end - t_begin);
+        }
+        else {
+            // the times of interest in each detector exclude this
+            // direction
+            s->log_skymap[i] = log(0);
+        }
+    }
+
+    s->p_free(p_t);
+    
+}
+

src/skymap/wanalysis.h

+// C headers
+
+#include <stdlib.h>
+
+// LAL headers
+
+#include <lal/Skymap.h>
+
+// Structure to hold analysis state
+
+typedef struct
+{
+
+    // Number of detectors <= XLALSKYMAP_N
+
+    size_t n_detectors;
+    
+    // Array holding LAL_DETECTOR codes to identify the instruments
+    
+    int detectors[XLALSKYMAP_N];
+
+    // Sample rate (integer Hz) of timeseries
+    
+    int rate;
+    
+    // Matched filter normalization
+    // Waveform * NoiseCovariance^{-1} * Waveform
+
+    double* wSw;
+    
+    // Time series unnormalized matched filter
+    // Data * NoiseCovariance^{-1} * Waveform
+    
+    double* xSw_real[XLALSKYMAP_N]; 
+    double* xSw_imag[XLALSKYMAP_N];
+    
+    // The waveform is assumed to be of the form A(t) e^{i phi(t)} with
+    // a cosine-like (real) and sine-like (imaginary) component that each
+    // have the same normalization term
+    
+    // Number of directions to compute for
+    
+    size_t n_directions;
+    
+    // Array of packed theta and phi values 
+    // {theta[0], phi[0], theta[1], phi[1], ... }
+    
+    double* directions; // Directions to sample at
+        
+    // Timing information
+    
+    double min_t; // Earth Barycenter time limits
+    double max_t;
+    
+    double delta_t; // Integration step size hint
+
+    double min_ts[XLALSKYMAP_N]; // Individual instrument time limits
+    double max_ts[XLALSKYMAP_N];
+
+    // Array to store the allowed amplitude calibration error
+    
+    double calibration_error[XLALSKYMAP_N];
+        
+    // Output array (with n_directions elements)
+
+    double* log_skymap;
+    
+    // Memory management hooks
+
+    void* (*p_realloc)(void*, size_t);
+    void (*p_free)(void*);
+    
+    
+} analysis;
+
+void analysis_default_construct(analysis* a);
+void analysis_default_directions(analysis* a);
+
+// Perform the analysis
+
+void analyze(analysis* s);
+
+// Convert between detector strings and LALDetectector identifiers
+
+int analysis_identify_detector(const char* c);
+

tests/Makefile.am

@@ -2,6 +2,7 @@ EXTRA_DIST = firbank_test_01.py fixtures.py plots_test_01.py resample_test_01.py
 
 check :
 	python test_coinc.py
+	python test_skymap.py
 	python firbank_test_01.py
 	python plots_test_01.py
 	python resample_test_01.py

tests/test_skymap.py

+#!/usr/bin/env python
+"""
+
+Unit tests for lal_skymap.
+
+"""
+__author__       = "Leo Singer <leo.singer@ligo.org>"
+__copyright__    = "Copyright 2010, Leo Singer"
+
+
+import unittest
+from fixtures import *
+from gstlal.pipeutil import *
+
+
+class TestSkymap(PipelineTestFixture):
+
+	def runTest(self):
+		pass # TODO
+
+if __name__ == '__main__':
+	suite = unittest.main()