bayestar_sky_map.c 55.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
/*                                           >y#
                                            ~'#o+
                                           '~~~md~
                '|+>#!~'':::::....        .~'~'cY#
            .+oy+>|##!~~~''':::......     ~:'':md! .
          #rcmory+>|#~''':::'::...::.::. :..'''Yr:...
        'coRRaamuyb>|!~'''::::':...........  .+n|.::..
       !maMMNMRYmuybb|!~'''':.........::::::: ro'..::..
      .cODDMYouuurub!':::...........:::~'.. |o>::...:..
      >BDNCYYmroyb>|#~:::::::::::::~':.:: :ob::::::::..
      uOCCNAa#'''''||':::.                :oy':::::::::.
    :rRDn!  :~::'y+::':  ... ...:::.     :ob':::::::::::.
   yMYy:   :>yooCY'.':.   .:'':......    ~u+~::::::::::::.
  >>:'. .~>yBDMo!.'': . .:'':.   .      >u|!:::::::::::::.
    ':'~|mYu#:'~'''. :.~:':...         yy>|~:::::::::::::..
    :!ydu>|!rDu::'. +'#~::!#'.~:     |r++>#':::::::::::::..
    mn>>>>>YNo:'': !# >'::::...  ..:cyb++>!:::::::::..:::...
    :ouooyodu:'': .!:.!:::.       yobbbb+>~::::::::....:....
     'cacumo~''' .'~ :~'.::.    :aybbbbbb>':::'~''::::....
      .mamd>'''. :~' :':'.:.   om>bbbyyyb>'.#b>|#~~~'':..
      .yYYo''': .:~' .'::'   .ny>+++byyoao!b+|||#!~~~''''''::.
      .#RUb:''. .:'' .:':   |a#|>>>>yBMdb #yb++b|':::::''':'::::::.
      .'CO!'''  .:'' .'    uu~##|+mMYy>+:|yyo+:::'::.         .::::::
      .:RB~''' ..::'.':   o>~!#uOOu>bby'|yB>.'::  '~!!!!!~':. ..  .::::
       :Rm''': ..:~:!:  'c~~+YNnbyyybb~'mr.':  !+yoy+>||!~'::.       :::.
      ..Oo''': .'' ~:  !+|BDCryuuuuub|#B!::  !rnYaocob|#!~'':.  ..    .::.
      . nB''': :  .'  |dNNduroomnddnuun::.  ydNAMMOary+>#~:.:::...      .:
       .uC~'''    :. yNRmmmadYUROMMBmm.:   bnNDDDMRBoy>|#~':....:.      .:
                 :' ymrmnYUROMAAAAMYn::. .!oYNDDMYmub|!~'::....:..     :
                 !'#booBRMMANDDDNNMO!:. !~#ooRNNAMMOOmuy+#!':::.......    :.
                .!'!#>ynCMNDDDDDNMRu.. '|:!raRMNAMOOdooy+|!~:::........   .:
                 : .'rdbcRMNNNNAMRB!:  |!:~bycmdYYBaoryy+|!~':::.::::.  ..
                 ..~|RMADnnONAMMRdy:. .>#::yyoroccruuybb>#!~'::':...::.
                  :'oMOMOYNMnybyuo!.  :>#::b+youuoyyy+>>|!~':.    :::::
                  ''YMCOYYNMOCCCRdoy##~~~: !b>bb+>>>||#~:..:::     ::::.
                  .:OMRCoRNAMOCROYYUdoy|>~:.~!!~!~~':...:'::::.   :::::.
                  ''oNOYyMNAMMMRYnory+|!!!:.....     ::.  :'::::::::::::
                 .:..uNabOAMMCOdcyb+|!~':::.          !!'.. :~:::::'''':.
                  .   +Y>nOORYauyy>!!'':....           !#~..  .~:''''''':.

****************  ____  _____  ______________________    ____     **************
***************  / __ )/   \ \/ / ____/ ___/_  __/   |  / __ \   ***************
**************  / __  / /| |\  / __/  \__ \ / / / /| | / /_/ /  ****************
*************  / /_/ / ___ |/ / /___ ___/ // / / ___ |/ _, _/  *****************
************  /_____/_/  |_/_/_____//____//_/ /_/  |_/_/ |_|  ******************
*/


/*
50
 * Copyright (C) 2013-2017  Leo Singer
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
 *
 * 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 with program; see the file COPYING. If not, write to the
 * Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
 * MA  02111-1307  USA
 */

68
#include "config.h"
69
#include "bayestar_sky_map.h"
70
#include "bayestar_distance.h"
71
#include "bayestar_moc.h"
72

73
#include <assert.h>
74
75
#include <float.h>
#include <math.h>
76
#include <stdlib.h>
77
78
79
#include <string.h>

#include <lal/DetResponse.h>
80
#include <lal/InspiralInjectionParams.h>
81
82
#include <lal/FrequencySeries.h>
#include <lal/LALSimInspiral.h>
83
#include <lal/LALSimulation.h>
84
85
86

#include <chealpix.h>

87
#include <gsl/gsl_cdf.h>
88
#include <gsl/gsl_integration.h>
89
#include <gsl/gsl_math.h>
90
#include <gsl/gsl_sf_bessel.h>
91
#include <gsl/gsl_sf_exp.h>
92
93
94
#include <gsl/gsl_sf_expint.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_test.h>
95
96
97

#include "logaddexp.h"

98
99
100
101
#ifndef _OPENMP
#define omp ignore
#endif

102

103
104
105
106
107
108
109
/* Compute |z|^2. Hopefully a little faster than gsl_pow_2(cabs(z)), because no
 * square roots are necessary. */
static double cabs2(double complex z) {
    return gsl_pow_2(creal(z)) + gsl_pow_2(cimag(z));
}


110
111
112
113
114
115
116
117
118
/*
 * Catmull-Rom cubic spline interpolant of x(t) for regularly gridded
 * samples x_i(t_i), assuming:
 *
 *     t_0 = -1, x_0 = x[0],
 *     t_1 = 0,  x_1 = x[1],
 *     t_2 = 1,  x_2 = x[2],
 *     t_3 = 2,  x_3 = x[3].
 */
119
static double complex complex_catrom(
120
121
122
123
    float complex x0,
    float complex x1,
    float complex x2,
    float complex x3,
124
125
126
127
128
129
130
131
132
    double t
) {
    return x1
        + t*(-0.5*x0 + 0.5*x2
        + t*(x0 - 2.5*x1 + 2.*x2 - 0.5*x3
        + t*(-0.5*x0 + 1.5*x1 - 1.5*x2 + 0.5*x3)));
}


133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
/*
 * Catmull-Rom cubic spline interpolant of x(t) for regularly gridded
 * samples x_i(t_i), assuming:
 *
 *     t_0 = -1, x_0 = x[0],
 *     t_1 = 0,  x_1 = x[1],
 *     t_2 = 1,  x_2 = x[2],
 *     t_3 = 2,  x_3 = x[3].
 *
 * I am careful to handle certain situations where some of
 * the samples are infinite or not-a-number:
 *
 *     * If t <= 0, then return x1.
 *     * If t >= 1, then return x2.
 *     * Otherwise, if either x0 or x3 are non-real, then fall back to
 *       linear interpolation between x1 and x2.
 *     * Otherwise, if x1 and/or x2 are infinite and both have the same
 *       then return infinity of that sign.
 *     * If x1 and x2 are infinities of different signs or if either are
 *       NaN, then return NaN.
 *     * Otherwise, if x0, x1, x2, and x3 are all finite, then return
 *       the standard Catmull-Rom formula.
 *
 * ***IMPORTANT NOTE*** I use the ISO C99 function isfinite() instead of
 * isinf(), the non-standard finite(), or gsl_finite(), because isfinite()
 * is reliably faster than the alternatives on Mac OS and Scientific Linux.
 *
 */
static double real_catrom(
    double x0,
    double x1,
    double x2,
    double x3,
    double t
) {
    double result;

    if (t <= 0)
        result = x1;
    else if (t >= 1)
        result = x2;
    else if (!(isfinite(x0 + x3)))
        result = x1 + (1 - t) * x2;
176
    else if (isfinite(x1 + x2))
177
178
179
180
181
182
183
184
185
186
187
        result = x1
            + t*(-0.5*x0 + 0.5*x2
            + t*(x0 - 2.5*x1 + 2.*x2 - 0.5*x3
            + t*(-0.5*x0 + 1.5*x1 - 1.5*x2 + 0.5*x3)));
    else
        result = x1 + x2;

    return result;
}


188
189
/* Evaluate a complex time series using cubic spline interpolation, assuming
 * that the vector x gives the samples of the time series at times
190
191
 * 0, 1, ..., nsamples-1. */
static double complex eval_snr(
192
    const float complex *x,
193
194
195
    size_t nsamples,
    double t
) {
196
    ssize_t i;
197
198
199
    double f;
    double complex y;

200
    /* Break |t| into integer and fractional parts. */
201
202
    {
        double dbl_i;
203
        f = modf(t, &dbl_i);
204
205
206
        i = dbl_i;
    }

207
    if (i >= 1 && i < (ssize_t)nsamples - 2)
208
        y = complex_catrom(x[i-1], x[i], x[i+1], x[i+2], f);
209
210
211
212
213
214
215
    else
        y = 0;

    return y;
}


216
217
218
#define INTERP_SIZE 400


219
typedef struct {
220
    double fx, xmin;
221
    double y[INTERP_SIZE];
222
223
224
225
226
227
} cubic_interp;


static double cubic_interp_eval(const cubic_interp *interp, double x)
{
    double i;
228
    const double u = modf((x - interp->xmin) * interp->fx, &i);
229

230
    #define CLAMP(_) ((int) ((_) < 0 ? 0 : ((_) < INTERP_SIZE ? (_) : INTERP_SIZE - 1)))
231
232
233
234
235
236
    const double y0 = interp->y[CLAMP(i - 1)];
    const double y1 = interp->y[CLAMP(i)];
    const double y2 = interp->y[CLAMP(i + 1)];
    const double y3 = interp->y[CLAMP(i + 2)];
    #undef CLAMP

237
    return real_catrom(y0, y1, y2, y3, u);
238
239
240
241
}


typedef struct {
242
    double xmin, ymin, fx, fy;
243
    double a[INTERP_SIZE][INTERP_SIZE][4][4];
244
245
246
} bicubic_interp;


247
static void bicubic_interp_init(bicubic_interp *interp, const double z[INTERP_SIZE][INTERP_SIZE], double xmin, double ymin, double dx, double dy)
248
{
249
250
251
252
    interp->xmin = xmin;
    interp->ymin = ymin;
    interp->fx = 1 / dx;
    interp->fy = 1 / dy;
253

254
255
256
    #define CLAMP(_) ((_) < 0 ? 0 : ((_) < INTERP_SIZE ? (_) : INTERP_SIZE - 1))

    for (int x = 0; x < INTERP_SIZE; x ++)
257
    {
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
        for (int y = 0; y < INTERP_SIZE; y ++)
        {
            double m[4][4], a[4][4];
            for (int x_ = 0; x_ < 4; x_ ++)
                for (int y_ = 0; y_ < 4; y_ ++)
                    m[x_][y_] = z[CLAMP(x + x_ - 1)][CLAMP(y + y_ - 1)];

            a[0][0] = m[1][1];
            a[0][1] = -0.5*m[0][1] + 0.5*m[2][1];
            a[0][2] = m[0][1] - 2.5*m[1][1] + 2.0*m[2][1] - 0.5*m[3][1];
            a[0][3] = -0.5*m[0][1] + 1.5*m[1][1] - 1.5*m[2][1] + 0.5*m[3][1];
            a[1][0] = -0.5*m[1][0] + 0.5*m[1][2];
            a[1][1] = 0.25*m[0][0] - 0.25*m[0][2] - 0.25*m[2][0] + 0.25*m[2][2];
            a[1][2] = -0.5*m[0][0] + 0.5*m[0][2] + 1.25*m[1][0] - 1.25*m[1][2] - 1.0*m[2][0] + 1.0*m[2][2] + 0.25*m[3][0] - 0.25*m[3][2];
            a[1][3] = 0.25*m[0][0] - 0.25*m[0][2] - 0.75*m[1][0] + 0.75*m[1][2] + 0.75*m[2][0] - 0.75*m[2][2] - 0.25*m[3][0] + 0.25*m[3][2];
            a[2][0] = m[1][0] - 2.5*m[1][1] + 2.0*m[1][2] - 0.5*m[1][3];
            a[2][1] = -0.5*m[0][0] + 1.25*m[0][1] - 1.0*m[0][2] + 0.25*m[0][3] + 0.5*m[2][0] - 1.25*m[2][1] + 1.0*m[2][2] - 0.25*m[2][3];
            a[2][2] = m[0][0] - 2.5*m[0][1] + 2.0*m[0][2] - 0.5*m[0][3] - 2.5*m[1][0] + 6.25*m[1][1] - 5.0*m[1][2] + 1.25*m[1][3] + 2.0*m[2][0] - 5.0*m[2][1] + 4.0*m[2][2] - 1.0*m[2][3] - 0.5*m[3][0] + 1.25*m[3][1] - 1.0*m[3][2] + 0.25*m[3][3];
            a[2][3] = -0.5*m[0][0] + 1.25*m[0][1] - 1.0*m[0][2] + 0.25*m[0][3] + 1.5*m[1][0] - 3.75*m[1][1] + 3.0*m[1][2] - 0.75*m[1][3] - 1.5*m[2][0] + 3.75*m[2][1] - 3.0*m[2][2] + 0.75*m[2][3] + 0.5*m[3][0] - 1.25*m[3][1] + 1.0*m[3][2] - 0.25*m[3][3];
            a[3][0] = -0.5*m[1][0] + 1.5*m[1][1] - 1.5*m[1][2] + 0.5*m[1][3];
            a[3][1] = 0.25*m[0][0] - 0.75*m[0][1] + 0.75*m[0][2] - 0.25*m[0][3] - 0.25*m[2][0] + 0.75*m[2][1] - 0.75*m[2][2] + 0.25*m[2][3];
            a[3][2] = -0.5*m[0][0] + 1.5*m[0][1] - 1.5*m[0][2] + 0.5*m[0][3] + 1.25*m[1][0] - 3.75*m[1][1] + 3.75*m[1][2] - 1.25*m[1][3] - 1.0*m[2][0] + 3.0*m[2][1] - 3.0*m[2][2] + 1.0*m[2][3] + 0.25*m[3][0] - 0.75*m[3][1] + 0.75*m[3][2] - 0.25*m[3][3];
            a[3][3] = 0.25*m[0][0] - 0.75*m[0][1] + 0.75*m[0][2] - 0.25*m[0][3] - 0.75*m[1][0] + 2.25*m[1][1] - 2.25*m[1][2] + 0.75*m[1][3] + 0.75*m[2][0] - 2.25*m[2][1] + 2.25*m[2][2] - 0.75*m[2][3] - 0.25*m[3][0] + 0.75*m[3][1] - 0.75*m[3][2] + 0.25*m[3][3];

            for (int x_ = 0; x_ < 4; x_ ++)
                for (int y_ = 0; y_ < 4; y_ ++)
                    interp->a[x][y][x_][y_] = a[x_][y_];
        }
286
287
    }

288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
    #undef CLAMP
}


static double bicubic_interp_eval(const bicubic_interp *interp, double x, double y)
{
    x -= interp->xmin;
    x *= interp->fx;
    y -= interp->ymin;
    y *= interp->fy;

    if (x < 0)
        x = 0;
    else if (x > INTERP_SIZE - 1)
        x = INTERP_SIZE - 1;
    if (y < 0)
        y = 0;
    else if (y > INTERP_SIZE - 1)
        y = INTERP_SIZE - 1;

    double i, j;
    x = modf(x, &i);
    y = modf(y, &j);

    const double (*a)[4] = interp->a[(int) i][(int) j];
    const double u[] = {1, x, x * x, x * x * x};
    const double v[] = {1, y, y * y, y * y * y};

    return v[0] * a[0][0] * u[0] + v[0] * a[0][1] * u[1]
         + v[0] * a[0][2] * u[2] + v[0] * a[0][3] * u[3]
         + v[1] * a[1][0] * u[0] + v[1] * a[1][1] * u[1]
         + v[1] * a[1][2] * u[2] + v[1] * a[1][3] * u[3]
         + v[2] * a[2][0] * u[0] + v[2] * a[2][1] * u[1]
         + v[2] * a[2][2] * u[2] + v[2] * a[2][3] * u[3]
         + v[3] * a[3][0] * u[0] + v[3] * a[3][1] * u[1]
         + v[3] * a[3][2] * u[2] + v[3] * a[3][3] * u[3];
324
325
326
}


327
typedef struct {
328
329
330
    bicubic_interp region0;
    cubic_interp region1;
    cubic_interp region2;
331
332
    double xmax, ymax, vmax, r1, r2;
    int k;
333
} log_radial_integrator;
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373


typedef struct {
    double scale;
    double p;
    double b;
    int k;
} radial_integrand_params;


static double radial_integrand(double r, void *params)
{
    const radial_integrand_params *integrand_params = params;
    const double scale = integrand_params->scale;
    const double p = integrand_params->p;
    const double b = integrand_params->b;
    const int k = integrand_params->k;
    return gsl_sf_exp_mult(
        scale - gsl_pow_2(p / r - 0.5 * b / p),
        gsl_sf_bessel_I0_scaled(b / r) * gsl_pow_int(r, k));
}


static double log_radial_integrand(double r, void *params)
{
    const radial_integrand_params *integrand_params = params;
    const double scale = integrand_params->scale;
    const double p = integrand_params->p;
    const double b = integrand_params->b;
    const int k = integrand_params->k;
    return log(gsl_sf_bessel_I0_scaled(b / r) * gsl_pow_int(r, k))
        + scale - gsl_pow_2(p / r - 0.5 * b / p);
}


static double log_radial_integral(double r1, double r2, double p, double b, int k)
{
    radial_integrand_params params = {0, p, b, k};
    double breakpoints[5];
    unsigned char nbreakpoints = 0;
374
    double result = 0, abserr, log_offset = -INFINITY;
375
376
    int ret;

377
    if (b != 0) {
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
        /* Calculate the approximate distance at which the integrand attains a
         * maximum (middle) and a fraction eta of the maximum (left and right).
         * This neglects the scaled Bessel function factors and the power-law
         * distance prior. It assumes that the likelihood is approximately of
         * the form
         *
         *    -p^2/r^2 + B/r.
         *
         * Then the middle breakpoint occurs at 1/r = -B/2A, and the left and
         * right breakpoints occur when
         *
         *   A/r^2 + B/r = log(eta) - B^2/4A.
         */

        static const double eta = 0.01;
        const double middle = 2 * gsl_pow_2(p) / b;
        const double left = 1 / (1 / middle + sqrt(-log(eta)) / p);
        const double right = 1 / (1 / middle - sqrt(-log(eta)) / p);

        /* Use whichever of the middle, left, and right points lie within the
         * integration limits as initial subdivisions for the adaptive
         * integrator. */

        breakpoints[nbreakpoints++] = r1;
        if(left > breakpoints[nbreakpoints-1] && left < r2)
            breakpoints[nbreakpoints++] = left;
        if(middle > breakpoints[nbreakpoints-1] && middle < r2)
            breakpoints[nbreakpoints++] = middle;
        if(right > breakpoints[nbreakpoints-1] && right < r2)
            breakpoints[nbreakpoints++] = right;
        breakpoints[nbreakpoints++] = r2;
409
410
411
412
    } else {
        /* Inner breakpoints are undefined because b = 0. */
        breakpoints[nbreakpoints++] = r1;
        breakpoints[nbreakpoints++] = r2;
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
    }

    /* Re-scale the integrand so that the maximum value at any of the
     * breakpoints is 1. Note that the initial value of the constant term
     * is overwritten. */

    for (unsigned char i = 0; i < nbreakpoints; i++)
    {
        double new_log_offset = log_radial_integrand(breakpoints[i], &params);
        if (new_log_offset > log_offset)
            log_offset = new_log_offset;
    }

    /* If the largest value of the log integrand was -INFINITY, then the
     * integrand is 0 everywhere. Set log_offset to 0, because subtracting
     * -INFINITY would make the integrand infinite. */
    if (log_offset == -INFINITY)
        log_offset = 0;

    params.scale = -log_offset;

    {
        /* Maximum number of subdivisions for adaptive integration. */
        static const size_t n = 64;

        /* Allocate workspace on stack. Hopefully, a little bit faster than
         * using the heap in multi-threaded code. */

        double alist[n];
        double blist[n];
        double rlist[n];
        double elist[n];
        size_t order[n];
        size_t level[n];
        gsl_integration_workspace workspace = {
            .alist = alist,
            .blist = blist,
            .rlist = rlist,
            .elist = elist,
            .order = order,
            .level = level,
            .limit = n
        };

        /* Set up integrand data structure. */
        const gsl_function func = {radial_integrand, &params};

        /* Perform adaptive Gaussian quadrature. */
        ret = gsl_integration_qagp(&func, breakpoints, nbreakpoints,
            DBL_MIN, 1e-8, n, &workspace, &result, &abserr);

        /* FIXME: do we care to keep the error estimate around? */
    }

    /* FIXME: do something with ret */
    (void)ret;

    /* Done! */
    return log_offset + log(result);
}


475
static void log_radial_integrator_init(log_radial_integrator *integrator, double r1, double r2, int k, double pmax)
476
477
478
479
480
481
482
483
{
    const double alpha = 4;
    const double p0 = 0.5 * (k >= 0 ? r2 : r1);
    const double xmax = log(pmax);
    const double x0 = GSL_MIN_DBL(log(p0), xmax);
    const double xmin = x0 - (1 + M_SQRT2) * alpha;
    const double ymax = x0 + alpha;
    const double ymin = 2 * x0 - M_SQRT2 * alpha - xmax;
484
    const double d = (xmax - xmin) / (INTERP_SIZE - 1); /* dx = dy = du */
485
486
    const double umin = - (1 + M_SQRT1_2) * alpha;
    const double vmax = x0 - M_SQRT1_2 * alpha;
487
    double z[INTERP_SIZE][INTERP_SIZE];
488
489
490
    /* const double umax = xmax - vmax; */ /* unused */

    #pragma omp parallel for
491
    for (size_t i = 0; i < INTERP_SIZE * INTERP_SIZE; i ++)
492
    {
493
494
        const size_t ix = i / INTERP_SIZE;
        const size_t iy = i % INTERP_SIZE;
495
496
497
498
499
500
        const double x = xmin + ix * d;
        const double y = ymin + iy * d;
        const double p = exp(x);
        const double r0 = exp(y);
        const double b = 2 * gsl_pow_2(p) / r0;
        /* Note: using this where p > r0; could reduce evaluations by half */
501
        z[ix][iy] = log_radial_integral(r1, r2, p, b, k);
502
    }
503
    bicubic_interp_init(&integrator->region0, z, xmin, ymin, d, d);
504

505
506
    integrator->region1.xmin = xmin;
    integrator->region1.fx = 1 / d;
507

508
509
    for (size_t i = 0; i < INTERP_SIZE; i ++)
        integrator->region1.y[i] = z[i][INTERP_SIZE - 1];
510

511
512
    integrator->region2.xmin = umin;
    integrator->region2.fx = 1 / d;
513

514
515
    for (size_t i = 0; i < INTERP_SIZE; i ++)
        integrator->region2.y[i] = z[i][INTERP_SIZE - 1 - i];
516
517
518
519

    integrator->xmax = xmax;
    integrator->ymax = ymax;
    integrator->vmax = vmax;
520
521
522
    integrator->r1 = r1;
    integrator->r2 = r2;
    integrator->k = k;
523
524
525
}


526
static double log_radial_integrator_eval(const log_radial_integrator *integrator, double p, double b, double log_p, double log_b)
527
{
528
529
    const double x = log_p;
    const double y = M_LN2 + 2 * log_p - log_b;
530
531
532
    double result;
    assert(x <= integrator->xmax);

533
534
535
536
537
538
    if (p == 0) {
        /* note: p2 == 0 implies b == 0 */
        assert(b == 0);
        int k1 = integrator->k + 1;

        if (k1 == 0)
539
        {
540
            result = log(log(integrator->r2 / integrator->r1));
541
        } else {
542
            result = log((gsl_pow_int(integrator->r2, k1) - gsl_pow_int(integrator->r1, k1)) / k1);
543
        }
544
545
    } else {
        if (y >= integrator->ymax) {
546
            result = cubic_interp_eval(&integrator->region1, x);
547
548
549
550
551
        } else {
            const double v = 0.5 * (x + y);
            if (v <= integrator->vmax)
            {
                const double u = 0.5 * (x - y);
552
                result = cubic_interp_eval(&integrator->region2, u);
553
            } else {
554
                result = bicubic_interp_eval(&integrator->region0, x, y);
555
556
557
            }
        }
        result += gsl_pow_2(0.5 * b / p);
558
559
    }

560
    return result;
561
562
563
}


564
/* Find error in time of arrival. */
565
566
static void toa_errors(
    double *dt,
567
568
569
    double theta,
    double phi,
    double gmst,
570
    int nifos,
571
572
    const double **locs, /* Input: detector position. */
    const double *toas /* Input: time of arrival. */
573
574
575
576
577
) {
    /* Convert to Cartesian coordinates. */
    double n[3];
    ang2vec(theta, phi - gmst, n);

578
    for (int i = 0; i < nifos; i ++)
579
    {
580
581
        double dot = 0;
        for (int j = 0; j < 3; j ++)
582
583
584
585
586
        {
            dot += locs[i][j] * n[j];
        }
        dt[i] = toas[i] + dot / LAL_C_SI;
    }
587
588
589
}


590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
/* Compute antenna factors from the detector response tensor and source
 * sky location, and return as a complex number F_plus + i F_cross. */
static double complex complex_antenna_factor(
    const float response[3][3],
    double ra,
    double dec,
    double gmst
) {
    double complex F;
    XLALComputeDetAMResponse(
        (double *)&F,     /* Type-punned real part */
        1 + (double *)&F, /* Type-punned imag part */
        response, ra, dec, 0, gmst);
    return F;
}


607
608
609
610
611
612
/* Expression for complex amplitude on arrival (without 1/distance factor) */
static double complex signal_amplitude_model(
    double complex F,               /* Antenna factor */
    double complex exp_i_twopsi,    /* e^(i*2*psi), for polarization angle psi */
    double u,                       /* cos(inclination) */
    double u2                       /* cos^2(inclination */
613
) {
614
    const double complex tmp = F * conj(exp_i_twopsi);
615
    return 0.5 * (1 + u2) * creal(tmp) - I * u * cimag(tmp);
616
617
618
}


619
620
621
622
static const unsigned int nglfixed = 10;
static const unsigned int ntwopsi = 10;


623
624
625
626
627
static double complex exp_i(double phi) {
    return cos(phi) + I * sin(phi);
}


628
629
630
631
632
/* Compare two pixels by contained probability. */
static int bayestar_pixel_compare_prob(const void *a, const void *b)
{
    const bayestar_pixel *apix = a;
    const bayestar_pixel *bpix = b;
633

634
635
    const double delta_logp = (apix->value[0] - bpix->value[0])
        - 2 * M_LN2 * (uniq2order64(apix->uniq) - uniq2order64(bpix->uniq));
636

637
638
639
640
641
642
643
    if (delta_logp < 0)
        return -1;
    else if (delta_logp > 0)
        return 1;
    else
        return 0;
}
644
645


646
static void bayestar_pixels_sort_prob(bayestar_pixel *pixels, size_t len)
647
{
648
649
    qsort(pixels, len, sizeof(bayestar_pixel), bayestar_pixel_compare_prob);
}
650
651


652
653
654
655
656
657
658
659
660
661
662
663
/* Compare two pixels by contained probability. */
static int bayestar_pixel_compare_uniq(const void *a, const void *b)
{
    const bayestar_pixel *apix = a;
    const bayestar_pixel *bpix = b;
    const unsigned long long auniq = apix->uniq;
    const unsigned long long buniq = bpix->uniq;

    if (auniq < buniq)
        return -1;
    else if (auniq > buniq)
        return 1;
664
    else
665
        return 0;
666
667
668
}


669
static void bayestar_pixels_sort_uniq(bayestar_pixel *pixels, size_t len)
670
{
671
    qsort(pixels, len, sizeof(bayestar_pixel), bayestar_pixel_compare_uniq);
672
673
674
675
676
677
678
679
680
}


static void *realloc_or_free(void *ptr, size_t size)
{
    void *new_ptr = realloc(ptr, size);
    if (!new_ptr)
    {
        free(ptr);
681
        XLAL_ERROR_NULL(XLAL_ENOMEM, "not enough memory to resize array");
682
683
684
685
686
687
    }
    return new_ptr;
}


/* Subdivide the final last_n pixels of an adaptively refined sky map. */
688
689
static bayestar_pixel *bayestar_pixels_refine(
    bayestar_pixel *pixels, size_t *len, size_t last_n
690
) {
691
    assert(last_n <= *len);
692
693

    /* New length: adding 4*last_n new pixels, removing last_n old pixels. */
694
695
    const size_t new_len = *len + 3 * last_n;
    const size_t new_size = new_len * sizeof(bayestar_pixel);
696

697
698
    pixels = realloc_or_free(pixels, new_size);
    if (pixels)
699
700
701
    {
        for (size_t i = 0; i < last_n; i ++)
        {
702
            const uint64_t uniq = 4 * pixels[*len - i - 1].uniq;
703
            for (unsigned char j = 0; j < 4; j ++)
704
                pixels[new_len - (4 * i + j) - 1].uniq = j + uniq;
705
        }
706
        *len = new_len;
707
    }
708
    return pixels;
709
710
711
}


712
static bayestar_pixel *bayestar_pixels_alloc(size_t *len, unsigned char order)
713
{
714
715
716
    const uint64_t nside = (uint64_t)1 << order;
    const uint64_t npix = nside2npix64(nside);
    const size_t size = npix * sizeof(bayestar_pixel);
717

718
719
    bayestar_pixel *pixels = malloc(size);
    if (!pixels)
720
        XLAL_ERROR_NULL(XLAL_ENOMEM, "not enough memory to allocate sky map");
721

722
723
724
725
    *len = npix;
    for (unsigned long long ipix = 0; ipix < npix; ipix ++)
        pixels[ipix].uniq = nest2uniq64(order, ipix);
    return pixels;
726
727
728
}


729
730
bayestar_pixel *bayestar_sky_map_toa_phoa_snr(
    size_t *out_len,                /* Number of returned pixels */
731
732
733
734
    /* Prior */
    double min_distance,            /* Minimum distance */
    double max_distance,            /* Maximum distance */
    int prior_distance_power,       /* Power of distance in prior */
735
    /* Data */
736
737
    double gmst,                    /* GMST (rad) */
    unsigned int nifos,             /* Number of detectors */
738
    unsigned long nsamples,         /* Length of SNR series */
739
    double sample_rate,             /* Sample rate in seconds */
740
    const double *epochs,           /* Timestamps of SNR time series */
741
    const float complex **snrs,     /* Complex SNR series */
742
743
    const float (**responses)[3],   /* Detector responses */
    const double **locations,       /* Barycentered Cartesian geographic detector positions (m) */
744
    const double *horizons          /* SNR=1 horizon distances for each detector */
745
) {
746
    log_radial_integrator *integrators[] = {NULL, NULL, NULL};
747
748
749
750
751
752
753
    {
        double pmax = 0;
        for (unsigned int iifo = 0; iifo < nifos; iifo ++)
        {
            pmax += gsl_pow_2(horizons[iifo] / max_distance);
        }
        pmax = sqrt(0.5 * pmax);
754
755
        for (unsigned char k = 0; k < 3; k ++)
        {
756
            integrators[k] = malloc(sizeof(log_radial_integrator));
757
758
759
            if (!integrators[k])
            {
                for (unsigned char kk = 0; kk < k; kk ++)
760
                    free(integrators[kk]);
761
762
                return NULL;
            }
763
764
765
            log_radial_integrator_init(
                integrators[k], min_distance, max_distance,
                prior_distance_power + k, pmax);
766
        }
767
768
    }

769
770
    static const unsigned char order0 = 4;
    unsigned char level = order0;
771
772
773
    size_t len;
    bayestar_pixel *pixels = bayestar_pixels_alloc(&len, order0);
    if (!pixels)
774
    {
775
        for (unsigned char k = 0; k < 3; k ++)
776
            free(integrators[k]);
777
        return NULL;
778
    }
779
    const unsigned long npix0 = len;
780

781
782
783
784
785
786
787
788
789
790
791
    /* Look up Gauss-Legendre quadrature rule for integral over cos(i). */
    gsl_integration_glfixed_table *gltable
        = gsl_integration_glfixed_table_alloc(nglfixed);

    /* Don't bother checking the return value. GSL has static, precomputed
     * values for certain orders, and for the order I have picked it will
     * return a pointer to one of these. See:
     *
     * http://git.savannah.gnu.org/cgit/gsl.git/tree/integration/glfixed.c
     */
    assert(gltable);
792
    assert(gltable->precomputed); /* We don't have to free it. */
793

794
    while (1)
795
    {
796
797
        #pragma omp parallel for
        for (unsigned long i = 0; i < npix0; i ++)
798
        {
799
            bayestar_pixel *const pixel = &pixels[len - npix0 + i];
800
801
            double complex F[nifos];
            double dt[nifos];
802
            double accum[3] = {-INFINITY, -INFINITY, -INFINITY};
803

804
805
            {
                double theta, phi;
806
                pix2ang_uniq64(pixel->uniq, &theta, &phi);
807

808
809
810
811
                /* Look up antenna factors */
                for (unsigned int iifo = 0; iifo < nifos; iifo++)
                    F[iifo] = complex_antenna_factor(
                        responses[iifo], phi, M_PI_2-theta, gmst) * horizons[iifo];
812

813
                toa_errors(dt, theta, phi, gmst, nifos, locations, epochs);
814
            }
815

816
817
            /* Integrate over 2*psi */
            for (unsigned int itwopsi = 0; itwopsi < ntwopsi; itwopsi++)
818
            {
819
820
                const double twopsi = (2 * M_PI / ntwopsi) * itwopsi;
                const double complex exp_i_twopsi = exp_i(twopsi);
821
                double accum1[3] = {-INFINITY, -INFINITY, -INFINITY};
822
823
824

                /* Integrate over u from -1 to 1. */
                for (unsigned int iu = 0; iu < nglfixed; iu++)
825
                {
826
                    double u, weight;
827
828
                    double accum2[nsamples][3];

829
830
831
832
                    {
                        /* Look up Gauss-Legendre abscissa and weight. */
                        int ret = gsl_integration_glfixed_point(
                            -1, 1, iu, &u, &weight, gltable);
833

834
835
836
837
838
839
840
                        /* Don't bother checking return value; the only
                         * possible failure is in index bounds checking. */
                        assert(ret == GSL_SUCCESS);
                    }
                    const double u2 = gsl_pow_2(u);
                    double complex z_times_r[nifos];
                    double p2 = 0;
841

842
843
844
845
846
847
848
                    for (unsigned int iifo = 0; iifo < nifos; iifo ++)
                    {
                        p2 += cabs2(
                            z_times_r[iifo] = signal_amplitude_model(
                                F[iifo], exp_i_twopsi, u, u2));
                    }
                    p2 *= 0.5;
849
850
                    const double p = sqrt(p2);
                    const double log_p = log(p);
851

852
                    for (long isample = 0;
853
                        isample < (long)nsamples; isample++)
854
                    {
855
                        double b;
856
                        {
857
858
859
                            double complex I0arg_complex_times_r = 0;
                            for (unsigned int iifo = 0; iifo < nifos; iifo ++)
                            {
860
861
862
                                I0arg_complex_times_r += conj(z_times_r[iifo])
                                    * eval_snr(snrs[iifo], nsamples,
                                        isample - dt[iifo] * sample_rate - 0.5 * (nsamples - 1));
863
864
                            }
                            b = cabs(I0arg_complex_times_r);
865
                        }
866
                        const double log_b = log(b);
867

868
869
                        for (unsigned char k = 0; k < 3; k ++)
                        {
870
                            accum2[isample][k] = log_radial_integrator_eval(
871
                                integrators[k], p, b, log_p, log_b);
872
                        }
873
                    }
874

875
876
877
878
879
880
881
882
883
                    double max_accum2[3] = {-INFINITY, -INFINITY, -INFINITY};
                    for (long isample = 0; isample < 3; isample ++)
                        for (unsigned char k = 0; k < 3; k ++)
                            if (accum2[isample][k] > max_accum2[k])
                                max_accum2[k] = accum2[isample][k];
                    double sum_accum2[3] = {0, 0, 0};
                    for (long isample = 0; isample < 3; isample ++)
                        for (unsigned char k = 0; k < 3; k ++)
                            sum_accum2[k] += exp(accum2[isample][k] - max_accum2[k]);
884
                    for (unsigned char k = 0; k < 3; k ++)
Leo P. Singer's avatar
Leo P. Singer committed
885
                        accum1[k] = logaddexp(accum1[k], log(sum_accum2[k] * weight) + max_accum2[k]);
886
                }
887

888
889
890
891
                for (unsigned char k = 0; k < 3; k ++)
                {
                    accum[k] = logaddexp(accum[k], accum1[k]);
                }
892
            }
893

894
            /* Record logarithm base 4 of posterior. */
895
896
            for (unsigned char k = 0; k < 3; k ++)
            {
897
                pixel->value[k] = accum[k];
898
            }
899
900
        }

Leo Singer's avatar
Leo Singer committed
901
        /* Sort pixels by ascending posterior probability. */
902
        bayestar_pixels_sort_prob(pixels, len);
903

904
905
906
        /* If we have reached order=11 (nside=2048), stop. */
        if (level++ >= 11)
            break;
907

908
        /* Adaptively refine the pixels that contain the most probability. */
909
910
        pixels = bayestar_pixels_refine(pixels, &len, npix0 / 4);
        if (!pixels)
911
            return NULL;
912
913
    }

914
    for (unsigned char k = 0; k < 3; k ++)
915
        free(integrators[k]);
916

917
918
919
920
921
922
923
924
925
926
    /* Rescale so that log(max) = 0. */
    const double max_logp = pixels[len - 1].value[0];
    for (ssize_t i = (ssize_t)len - 1; i >= 0; i --)
        for (unsigned char k = 0; k < 3; k ++)
            pixels[i].value[k] -= max_logp;

    /* Determine normalization of map. */
    double norm = 0;
    for (ssize_t i = (ssize_t)len - 1; i >= 0; i --)
    {
927
        const double dA = uniq2pixarea64(pixels[i].uniq);
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
        const double dP = gsl_sf_exp_mult(pixels[i].value[0], dA);
        if (dP <= 0)
            break; /* We have reached underflow. */
        norm += dP;
    }
    norm = 1 / norm;

    /* Rescale, normalize, and prepare output. */
    for (ssize_t i = (ssize_t)len - 1; i >= 0; i --)
    {
        const double prob = gsl_sf_exp_mult(pixels[i].value[0], norm);
        double rmean = exp(pixels[i].value[1] - pixels[i].value[0]);
        double rstd = exp(pixels[i].value[2] - pixels[i].value[0]) - gsl_pow_2(rmean);
        if (rstd >= 0)
        {
            rstd = sqrt(rstd);
        } else {
            rmean = INFINITY;
            rstd = 1;
        }
        pixels[i].value[0] = prob;
        pixels[i].value[1] = rmean;
        pixels[i].value[2] = rstd;
    }

    /* Sort pixels by ascending NUNIQ index. */
    bayestar_pixels_sort_uniq(pixels, len);
955
956

    /* Done! */
957
958
    *out_len = len;
    return pixels;
959
960
961
}


962
963
964
965
966
967
968
969
double bayestar_log_likelihood_toa_phoa_snr(
    /* Parameters */
    double ra,                      /* Right ascension (rad) */
    double sin_dec,                 /* Sin(declination) */
    double distance,                /* Distance */
    double u,                       /* Cos(inclination) */
    double twopsi,                  /* Twice polarization angle (rad) */
    double t,                       /* Barycentered arrival time (s) */
970
    /* Data */
971
972
    double gmst,                    /* GMST (rad) */
    unsigned int nifos,             /* Number of detectors */
973
    unsigned long nsamples,         /* Lengths of SNR series */
974
    double sample_rate,             /* Sample rate in seconds */
975
    const double *epochs,           /* Timestamps of SNR time series */
976
    const float complex **snrs,     /* Complex SNR series */
977
978
    const float (**responses)[3],   /* Detector responses */
    const double **locations,       /* Barycentered Cartesian geographic detector positions (m) */
979
    const double *horizons          /* SNR=1 horizon distances for each detector */
980
981
982
983
984
985
986
987
) {
    const double dec = asin(sin_dec);
    const double u2 = gsl_pow_2(u);
    const double complex exp_i_twopsi = exp_i(twopsi);
    const double one_by_r = 1 / distance;

    /* Compute time of arrival errors */
    double dt[nifos];
988
    toa_errors(dt, M_PI_2 - dec, ra, gmst, nifos, locations, epochs);
989
990

    double complex i0arg_complex_times_r = 0;
991
992
993
    double A = 0;

    /* Loop over detectors */
994
    for (unsigned int iifo = 0; iifo < nifos; iifo++)
995
    {
996
997
        const double complex F = complex_antenna_factor(
            responses[iifo], ra, dec, gmst) * horizons[iifo];
998

999
1000
        const double complex z_times_r =
             signal_amplitude_model(F, exp_i_twopsi, u, u2);
For faster browsing, not all history is shown. View entire blame