Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALInference 4.1.9.1-5e288d3
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LALInferenceProposal.c
Go to the documentation of this file.
1/*
2 * LALInferenceProposal.c: Bayesian Followup, jump proposals.
3 *
4 * Copyright (C) 2011 Ilya Mandel, Vivien Raymond, Christian Roever,
5 * Marc van der Sluys, John Veitch, Will M. Farr, Ben Farr
6 *
7 *
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with with program; see the file COPYING. If not, write to the
20 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
21 * MA 02110-1301 USA
22 */
23
24#include <stdio.h>
25#include <stdlib.h>
26#include <lal/LALInspiral.h>
27#include <lal/DetResponse.h>
28#include <lal/SeqFactories.h>
29#include <lal/Date.h>
30#include <lal/VectorOps.h>
31#include <lal/TimeFreqFFT.h>
32#include <lal/GenerateInspiral.h>
33#include <lal/TimeDelay.h>
34#include <lal/SkyCoordinates.h>
35#include <lal/LALInference.h>
36#include <lal/LALInferenceInit.h>
37#include <lal/LALInferencePrior.h>
38#include <lal/LALInferenceLikelihood.h>
39#include <lal/LALInferenceTemplate.h>
40#include <lal/LALInferenceProposal.h>
41#include <lal/LALDatatypes.h>
42#include <lal/FrequencySeries.h>
43#include <lal/LALSimInspiral.h>
44#include <lal/LALSimNoise.h>
45#include <lal/XLALError.h>
46
47#include <lal/LALStdlib.h>
48#include <lal/LALInferenceClusteredKDE.h>
49#include <lal/LALInferenceNestedSampler.h>
50
51#ifdef __GNUC__
52#define UNUSED __attribute__ ((unused))
53#else
54#define UNUSED
55#endif
56
57typedef enum {
61
62const char *const cycleArrayName = "Proposal Cycle";
63const char *const cycleArrayLengthName = "Proposal Cycle Length";
64const char *const cycleArrayCounterName = "Proposal Cycle Counter";
65
66/* Proposal Names */
67const char *const nullProposalName = "NULL";
68const char *const singleAdaptProposalName = "Single";
69const char *const singleProposalName = "Single";
70const char *const orbitalPhaseJumpName = "OrbitalPhase";
71const char *const covarianceEigenvectorJumpName = "CovarianceEigenvector";
72const char *const skyLocWanderJumpName = "SkyLocWander";
73const char *const differentialEvolutionFullName = "DifferentialEvolutionFull";
74const char *const differentialEvolutionIntrinsicName = "DifferentialEvolutionIntrinsic";
75const char *const differentialEvolutionExtrinsicName = "DifferentialEvolutionExtrinsic";
76const char *const ensembleStretchFullName = "EnsembleStretchFull";
77const char *const ensembleStretchIntrinsicName = "EnsembleStretchIntrinsic";
78const char *const ensembleStretchExtrinsicName = "EnsembleStretchExtrinsic";
79const char *const drawApproxPriorName = "DrawApproxPrior";
80const char *const drawFlatPriorName = "DrawFlatPrior";
81const char *const skyReflectDetPlaneName = "SkyReflectDetPlane";
82const char *const skyRingProposalName = "SkyRingProposal";
83const char *const PSDFitJumpName = "PSDFitJump";
84const char *const polarizationPhaseJumpName = "PolarizationPhase";
85const char *const polarizationCorrPhaseJumpName = "CorrPolarizationPhase";
86const char *const extrinsicParamProposalName = "ExtrinsicParamProposal";
87const char *const frequencyBinJumpName = "FrequencyBin";
88const char *const GlitchMorletJumpName = "glitchMorletJump";
89const char *const GlitchMorletReverseJumpName = "glitchMorletReverseJump";
90const char *const ensembleWalkFullName = "EnsembleWalkFull";
91const char *const ensembleWalkIntrinsicName = "EnsembleWalkIntrinsic";
92const char *const ensembleWalkExtrinsicName = "EnsembleWalkExtrinsic";
93const char *const clusteredKDEProposalName = "ClusteredKDEProposal";
94const char *const splineCalibrationProposalName = "SplineCalibration";
95const char *const distanceLikelihoodProposalName = "DistanceLikelihood";
96
97static const char *intrinsicNames[] = {"chirpmass", "q", "eta", "mass1", "mass2", "a_spin1", "a_spin2","tilt_spin1", "tilt_spin2", "phi12", "phi_jl", "frequency", "quality", "duration","polar_angle", "phase", "polar_eccentricity","dchiMinus2","dchiMinus1","dchi0","dchi1","dchi2","dchi3","dchi3S","dchi3NS","dchi4","dchi4S","dchi4NS","dchi5","dchi5S","dchi5NS","dchi5l","dchi5lS","dchi5lNS","dchi6","dchi6S","dchi6NS","dchi6l","dchi7","dchi7S","dchi7NS","aPPE","alphaPPE","bPPE","betaPPE","betaStep","fStep","dxi1","dxi2","dxi3","dxi4","dxi5","dxi6","dalpha1","dalpha2","dalpha3","dalpha4","dalpha5","dbeta1","dbeta2","dbeta3","dsigma1","dsigma2","dsigma3","dsigma4","lambda1","lambda2","lambdaT","dlambdaT","logp1", "gamma1", "gamma2", "gamma3", "SDgamma0","SDgamma1","SDgamma2","SDgamma3","log10lambda_eff","lambda_eff","nonGR_alpha","LIV_A_sign","dQuadMon1","dQuadMon2","dQuadMonA","dQuadMonS","dchikappaS","dchikappaA", "domega220", "dtau220", "domega210", "dtau210", "domega330", "dtau330", "domega440", "dtau440", "domega550", "dtau550", "db1", "db2", "db3", "db4", "dc1", "dc2", "dc4", "dcl", NULL};
98
99
100static const char *extrinsicNames[] = {"rightascension", "declination", "cosalpha", "azimuth", "polarisation", "distance",
101 "logdistance", "time", "costheta_jn", "t0", "theta","hrss", "loghrss", NULL};
102
104 INT4 i;
105
106 for (i = 0; i < 3; i++) {
107 if (d1->location[i] != d2->location[i])
108 return 0;
109 }
110
111 return 1;
112}
113
115 INT4 nIFO = 0;
116 INT4 nCollision = 0;
117 LALInferenceIFOData *currentIFO = NULL;
118
119 for (currentIFO = data; currentIFO; currentIFO = currentIFO->next) {
120 LALInferenceIFOData *subsequentIFO = NULL;
121 nIFO++;
122 for (subsequentIFO = currentIFO->next; subsequentIFO; subsequentIFO = subsequentIFO->next) {
123 if (same_detector_location(subsequentIFO->detector, currentIFO->detector)) {
124 nCollision++;
125 break;
126 }
127 }
128 }
129
130 return nIFO - nCollision;
131}
132
134{
136 proposal->func = func;
137 proposal->proposed = 0;
138 proposal->accepted = 0;
139 strcpy(proposal->name, name);
140 return proposal;
141}
142
143
144void LALInferenceRegisterProposal(LALInferenceVariables *propArgs, const char *name, INT4 *flag, ProcessParamsTable *command_line) {
145 char offopt[VARNAME_MAX+15];
146 char onopt[VARNAME_MAX+12];
147
148 sprintf(offopt, "--proposal-no-%s", name);
149 sprintf(onopt, "--proposal-%s", name);
150
151 if (LALInferenceGetProcParamVal(command_line, offopt))
152 *flag = 0;
153 else if (LALInferenceGetProcParamVal(command_line, onopt))
154 *flag = 1;
155
157}
158
160 const char *fname = "LALInferenceAddProposalToCycle";
161 INT4 i;
162
163 /* Quit without doing anything if weight = 0. */
164 if (weight == 0)
165 return;
166
167 cycle->order = XLALRealloc(cycle->order, (cycle->length + weight)*sizeof(INT4));
168 if (cycle->order == NULL) {
169 XLALError(fname, __FILE__, __LINE__, XLAL_ENOMEM);
170 exit(1);
171 }
172
173 for (i = cycle->length; i < cycle->length + weight; i++) {
174 cycle->order[i] = cycle->nProposals;
175 }
176
177 cycle->nProposals += 1;
178 cycle->proposals = XLALRealloc(cycle->proposals, (cycle->nProposals)*sizeof(LALInferenceProposal));
179 if (cycle->proposals == NULL) {
180 XLALError(fname, __FILE__, __LINE__, XLAL_ENOMEM);
181 exit(1);
182 }
183 cycle->proposals[cycle->nProposals-1] = prop;
184
185 cycle->length += weight;
186}
187
188
189
191 INT4 i, j, temp;
192
193 for (i = cycle->length - 1; i > 0; i--) {
194 /* Fill in array from right to left, chosen randomly from remaining proposals. */
195 j = gsl_rng_uniform_int(rng, i+1);
196
197 temp = cycle->order[j];
198 cycle->order[j] = cycle->order[i];
199 cycle->order[i] = temp;
200 }
201}
202
203
206 LALInferenceVariables *proposedParams) {
207 INT4 i = 0;
208 LALInferenceProposalCycle *cycle=NULL;
209
210 /* Must have cycle array and cycle array length in propArgs. */
211 cycle = thread->cycle;
212 if (cycle == NULL) {
213 XLALError("LALInferenceCyclicProposal()",__FILE__,__LINE__,XLAL_FAILURE);
214 exit(1);
215 }
216
217 if (cycle->counter >= cycle->length) {
218 XLALError("LALInferenceCyclicProposal()",__FILE__,__LINE__,XLAL_FAILURE);
219 exit(1);
220 }
221
222 /* One instance of each proposal object is stored in cycle->proposals.
223 cycle->order is a list of elements to call from the proposals */
224
225 REAL8 logPropRatio=-INFINITY;
226 do
227 {
228 i = cycle->order[cycle->counter];
229 logPropRatio = cycle->proposals[i]->func(thread, currentParams, proposedParams);
230 strcpy(cycle->last_proposal_name, cycle->proposals[i]->name);
231 cycle->counter = (cycle->counter + 1) % cycle->length;
232 }
233 /* Call proposals until one succeeds */
234 while (proposedParams->head == NULL);
235
236 return logPropRatio;
237}
238
241 strcpy(cycle->last_proposal_name, nullProposalName);
242
243 return cycle;
244}
245
247 XLALFree(cycle->proposals);
248 XLALFree(cycle->order);
249}
250
253 LALInferenceIFOData *ifo = runState->data;
254
255 /* This will copy any existing arguments over from runState. I (John) don't think this should be necessary
256 * as this function is used to initialise these arguments in the first place. */
258 if(runState->proposalArgs && runState->proposalArgs->dimension>0) LALInferenceCopyVariables(runState->proposalArgs, propArgs);
259
260 INT4 Nskip = 1;
261 INT4 noise_only = 0;
262 INT4 cyclic_reflective_kde = 0;
263
264 /* Flags for proposals, initialized with the MCMC defaults */
265
266 INT4 singleadapt = 1; /* Disabled for bug checking */
267 INT4 psiphi = 1;
268 INT4 ext_param = 1;
269 INT4 skywander = 1;
270 INT4 skyreflect = 1;
271 INT4 drawprior = 1;
272 INT4 covjump = 0;
273 INT4 diffevo = 1;
274 INT4 stretch = 0;
275 INT4 walk = 0;
276 INT4 skyring = 1;
277 INT4 distance = 1;
278 INT4 kde = 0;
279 INT4 spline_cal = 0;
280 INT4 psdfit = 0;
281 INT4 glitchfit = 0;
282
284 singleadapt = 0;
285 psiphi = 0;
286 ext_param = 0;
287 skywander = 0;
288 skyreflect = 0;
289 drawprior = 0;
290 covjump = 1;
291 diffevo = 1;
292 stretch = 1;
293 walk = 1;
294 skyring = 0;
295 distance = 1;
296 kde = 0;
297 spline_cal = 0;
298 psdfit = 0;
299 glitchfit = 0;
300 }
301 if (LALInferenceCheckVariable(runState->algorithmParams, "LIB"))
302 distance=0;
303
304 ProcessParamsTable *command_line = runState->commandLine;
305
306 if(LALInferenceGetProcParamVal(command_line, "--margdist")
307 ||LALInferenceGetProcParamVal(command_line,"--margdist-comoving"))
308 distance=0;
309
310 INT4 verbose = 0;
311 if (LALInferenceGetProcParamVal(command_line, "--verbose"))
312 verbose = 1;
314
315 /* Determine the number of iterations between each entry in the DE buffer */
316 if (LALInferenceCheckVariable(runState->algorithmParams, "Nskip"))
317 Nskip = LALInferenceGetINT4Variable(runState->algorithmParams, "Nskip");
318 LALInferenceAddINT4Variable(propArgs, "Nskip", Nskip, LALINFERENCE_PARAM_FIXED);
319
320 /* Count the number of IFOs and uniquely-located IFOs to decide which sky-related proposals to use */
321 INT4 nDet = 0;
322 ifo=runState->data;
323 while (ifo) {
324 nDet++;
325 ifo = ifo->next;
326 }
328
329 INT4 nUniqueDet = numDetectorsUniquePositions(runState->data);
330 LALInferenceAddINT4Variable(propArgs, "nUniqueDet", nUniqueDet, LALINFERENCE_PARAM_FIXED);
331
332 INT4 marg_timephi = 0;
333 if (LALInferenceGetProcParamVal(command_line, "--margtimephi"))
334 marg_timephi = 1;
335
336 INT4 marg_time = 0;
337 if (marg_timephi || LALInferenceGetProcParamVal(command_line, "--margtime"))
338 marg_time = 1;
339 LALInferenceAddINT4Variable(propArgs, "marg_time", marg_time, LALINFERENCE_PARAM_FIXED);
340
341 INT4 marg_phi = 0;
342 if (marg_timephi || LALInferenceGetProcParamVal(command_line, "--margphi"))
343 marg_phi = 1;
344 LALInferenceAddINT4Variable(propArgs, "marg_phi", marg_phi, LALINFERENCE_PARAM_FIXED);
345
346 INT4 analytic_test = 0;
347 if (LALInferenceGetProcParamVal(command_line, "--correlatedGaussianLikelihood") ||
348 LALInferenceGetProcParamVal(command_line, "--bimodalGaussianLikelihood") ||
349 LALInferenceGetProcParamVal(command_line, "--rosenbrockLikelihood")) {
350 analytic_test = 1;
351 distance = 0;
352 }
353 LALInferenceAddINT4Variable(propArgs, "analytical_test", analytic_test, LALINFERENCE_PARAM_FIXED);
354
355 INT4 skyframe = 1;
356 if (LALInferenceGetProcParamVal(command_line, "--no-sky-frame"))
357 skyframe = 0;
358
359 INT4 noAdapt = 0;
360 if (LALInferenceGetProcParamVal(command_line, "--no-adapt") ||
361 LALInferenceGetProcParamVal(command_line, "--noiseonly"))
362 noAdapt = 1;
363 INT4 adapting = !noAdapt;
364 LALInferenceAddINT4Variable(propArgs, "no_adapt", noAdapt, LALINFERENCE_PARAM_LINEAR);
365 LALInferenceAddINT4Variable(propArgs, "adapting", adapting, LALINFERENCE_PARAM_LINEAR);
366
367 INT4 tau = 5;
368 ppt = LALInferenceGetProcParamVal(command_line, "--adapt-tau");
369 if (ppt)
370 tau = atof(ppt->value);
371 LALInferenceAddINT4Variable(propArgs, "adaptTau", tau, LALINFERENCE_PARAM_FIXED);
372
373 INT4 sampling_prior = 0;
374 ppt = LALInferenceGetProcParamVal(command_line, "--zerologlike");
375 if (ppt)
376 sampling_prior = 1;
377 LALInferenceAddINT4Variable(propArgs, "sampling_prior", sampling_prior, LALINFERENCE_PARAM_FIXED);
378
379 if (LALInferenceGetProcParamVal(command_line, "--enable-spline-calibration"))
380 spline_cal = 1;
381
382 ppt = LALInferenceGetProcParamVal(command_line, "--psd-fit");
383 if (!ppt)
384 ppt = LALInferenceGetProcParamVal(command_line, "--psdFit");
385 if (ppt)
386 psdfit = 1;
387
388 if (LALInferenceGetProcParamVal(command_line, "--glitch-fit"))
389 glitchfit = 1;
390
391 /* Convenience option for turning off ensemble moves */
392 if (LALInferenceGetProcParamVal(command_line, "--proposal-no-ensemble")) {
393 stretch = 0;
394 walk = 0;
395 }
396
397 /* Check if imposing cyclic reflective bounds */
398 if (LALInferenceGetProcParamVal(runState->commandLine, "--cyclic-reflective-kde"))
399 cyclic_reflective_kde = 1;
400 LALInferenceAddINT4Variable(propArgs, "cyclic_reflective_kde", cyclic_reflective_kde, LALINFERENCE_PARAM_FIXED);
401
402 if (LALInferenceGetProcParamVal(command_line, "--noiseonly"))
403 noise_only = 1;
404 LALInferenceAddINT4Variable(propArgs, "noiseonly", noise_only, LALINFERENCE_PARAM_FIXED);
405
406 /* Turn off signal proposals if no signal is in the model */
407 if (noise_only) {
408 singleadapt = 0;
409 psiphi = 0;
410 ext_param = 0;
411 skywander = 0;
412 skyreflect = 0;
413 drawprior = 0;
414 covjump = 0;
415 diffevo = 0;
416 stretch = 0;
417 walk = 0;
418 distance = 0;
419 skyring = 0;
420 spline_cal = 0;
421 }
422
423 /* Turn off phi-related proposals if marginalizing over phi in likelihood */
424 if (marg_phi) {
425 psiphi = 0;
426 }
427
428 /* Disable proposals that won't work with the current number of unique detectors */
429 if (nUniqueDet < 2) {
430 skyring = 0;
431 }
432
433 if (nUniqueDet != 3) {
434 skyreflect = 0;
435 }
436
437 if (nUniqueDet >= 3) {
438 ext_param = 0;
439 }
440
441 /* Turn off ra-dec related proposals when using the sky-frame coordinate system */
442 if (skyframe) {
443 ext_param = 0;
444 skywander = 0;
445 skyreflect = 0;
446 skyring = 0;
447 }
448
449 /* Register all proposal functions, check for explicit command-line requests */
450 LALInferenceRegisterProposal(propArgs, "single-adapt", &singleadapt, command_line);
451 LALInferenceRegisterProposal(propArgs, "psiphi", &psiphi, command_line);
452 LALInferenceRegisterProposal(propArgs, "extrinsic-param", &ext_param, command_line);
453 LALInferenceRegisterProposal(propArgs, "sky-wander", &skywander, command_line);
454 LALInferenceRegisterProposal(propArgs, "sky-reflect", &skyreflect, command_line);
455 LALInferenceRegisterProposal(propArgs, "draw-prior", &drawprior, command_line);
456 LALInferenceRegisterProposal(propArgs, "eigenvectors", &covjump, command_line);
457 LALInferenceRegisterProposal(propArgs, "differential-evolution", &diffevo, command_line);
458 LALInferenceRegisterProposal(propArgs, "ensemble-stretch", &stretch, command_line);
459 LALInferenceRegisterProposal(propArgs, "ensemble-walk", &walk, command_line);
460 LALInferenceRegisterProposal(propArgs, "sky-ring", &skyring, command_line);
461 LALInferenceRegisterProposal(propArgs, "distance", &distance, command_line);
462 LALInferenceRegisterProposal(propArgs, "kde", &kde, command_line);
463 LALInferenceRegisterProposal(propArgs, "spline_cal", &spline_cal, command_line);
464 LALInferenceRegisterProposal(propArgs, "psdfit", &psdfit, command_line);
465 LALInferenceRegisterProposal(propArgs, "glitchfit", &glitchfit, command_line);
466
467 /* Setup adaptive proposals */
468 if (singleadapt){
471 XLALFree(model);
472 }
473
474 /* Setup now since we need access to the data */
475 if (glitchfit)
476 LALInferenceSetupGlitchProposal(runState->data, propArgs);
477
478 return propArgs;
479}
480
481
484
485 const INT4 BIGWEIGHT = 20;
486 const INT4 SMALLWEIGHT = 5;
487 const INT4 TINYWEIGHT = 1;
488
490
491 if (LALInferenceGetINT4Variable(propArgs, "single-adapt")) {
493 LALInferenceAddProposalToCycle(cycle, prop, BIGWEIGHT);
494 }
495
496 if (LALInferenceGetINT4Variable(propArgs, "psiphi")) {
498 LALInferenceAddProposalToCycle(cycle, prop, TINYWEIGHT);
499 }
500
501 if (LALInferenceGetINT4Variable(propArgs, "extrinsic-param")) {
503 LALInferenceAddProposalToCycle(cycle, prop, SMALLWEIGHT);
504 }
505
506 if (LALInferenceGetINT4Variable(propArgs, "sky-wander")) {
508 LALInferenceAddProposalToCycle(cycle, prop, SMALLWEIGHT);
509 }
510
511 if (LALInferenceGetINT4Variable(propArgs, "sky-reflect")) {
513 LALInferenceAddProposalToCycle(cycle, prop, TINYWEIGHT);
514 }
515
516 if (LALInferenceGetINT4Variable(propArgs, "draw-prior")) {
518 LALInferenceAddProposalToCycle(cycle, prop, TINYWEIGHT);
519 }
520
521 if (LALInferenceGetINT4Variable(propArgs, "eigenvectors")) {
523 LALInferenceAddProposalToCycle(cycle, prop, BIGWEIGHT);
524 }
525
526 if (LALInferenceGetINT4Variable(propArgs, "differential-evolution")) {
528 LALInferenceAddProposalToCycle(cycle, prop, BIGWEIGHT);
529
531 LALInferenceAddProposalToCycle(cycle, prop, SMALLWEIGHT);
532
534 LALInferenceAddProposalToCycle(cycle, prop, SMALLWEIGHT);
535 }
536
537 if (LALInferenceGetINT4Variable(propArgs, "ensemble-stretch")) {
539 LALInferenceAddProposalToCycle(cycle, prop, BIGWEIGHT);
540
542 LALInferenceAddProposalToCycle(cycle, prop, SMALLWEIGHT);
543
545 LALInferenceAddProposalToCycle(cycle, prop, SMALLWEIGHT);
546 }
547
548 if (LALInferenceGetINT4Variable(propArgs, "ensemble-walk")) {
550 LALInferenceAddProposalToCycle(cycle, prop, BIGWEIGHT);
551
553 LALInferenceAddProposalToCycle(cycle, prop, SMALLWEIGHT);
554
556 LALInferenceAddProposalToCycle(cycle, prop, SMALLWEIGHT);
557 }
558
559 if (LALInferenceGetINT4Variable(propArgs, "sky-ring")) {
561 LALInferenceAddProposalToCycle(cycle, prop, SMALLWEIGHT);
562 }
563
564 /* Distance proposal */
565 if (LALInferenceGetINT4Variable(propArgs, "distance")) {
567 LALInferenceAddProposalToCycle(cycle, prop, SMALLWEIGHT);
568 }
569
570 if (LALInferenceGetINT4Variable(propArgs, "kde")) {
572 LALInferenceAddProposalToCycle(cycle, prop, BIGWEIGHT);
573 }
574
575 if (LALInferenceGetINT4Variable(propArgs, "spline_cal")) {
577 LALInferenceAddProposalToCycle(cycle, prop, SMALLWEIGHT);
578 }
579
580 if (LALInferenceGetINT4Variable(propArgs, "psdfit")) {
582 LALInferenceAddProposalToCycle(cycle, prop, SMALLWEIGHT);
583 }
584
585 if (LALInferenceGetINT4Variable(propArgs, "glitchfit")) {
587 LALInferenceAddProposalToCycle(cycle, prop, SMALLWEIGHT);
588
590 LALInferenceAddProposalToCycle(cycle, prop, SMALLWEIGHT);
591 }
592
593 return cycle;
594}
595
596
599 LALInferenceVariables *proposedParams) {
600 INT4 dim, varNr;
601 REAL8 logPropRatio, sqrttemp, mu, sigma;
602 char tmpname[MAX_STRLEN] = "";
603 LALInferenceVariableItem *param = NULL;
604
607 gsl_rng *rng = thread->GSLrandom;
608
609 if (!LALInferenceGetINT4Variable(args, "no_adapt")) {
610 if (!LALInferenceCheckVariable(args, "adapting"))
612
613 sqrttemp = sqrt(thread->temperature);
614 dim = proposedParams->dimension;
615
616 do {
617 varNr = 1 + gsl_rng_uniform_int(rng, dim);
618 param = LALInferenceGetItemNr(proposedParams, varNr);
619 } while (!LALInferenceCheckVariableNonFixed(proposedParams, param->name) || param->type != LALINFERENCE_REAL8_t);
620
621 if (param->type != LALINFERENCE_REAL8_t) {
622 fprintf(stderr, "Attempting to set non-REAL8 parameter with numerical sigma (in %s, %d)\n",
623 __FILE__, __LINE__);
624 exit(1);
625 }
626
627 sprintf(tmpname,"%s_%s",param->name,ADAPTSUFFIX);
628 if (!LALInferenceCheckVariable(thread->proposalArgs, tmpname)) {
629 fprintf(stderr, "Attempting to draw single-parameter jump for %s but cannot find sigma!\nError in %s, line %d.\n",
630 param->name,__FILE__, __LINE__);
631 exit(1);
632 }
633
634 sigma = LALInferenceGetREAL8Variable(thread->proposalArgs, tmpname);
635
636 /* Save the name of the proposed variable */
637 LALInferenceAddstringVariable(args, "proposedVariableName", param->name, LALINFERENCE_PARAM_OUTPUT);
638
639 /* If temperature is infinite, scale sigma to the prior width */
640 if (sqrttemp == INFINITY) {
641 if (LALInferenceCheckGaussianPrior(thread->priorArgs, param->name)) {
642 LALInferenceGetGaussianPrior(thread->priorArgs, param->name, &mu, &sigma);
643
644 *((REAL8 *)param->value) += gsl_ran_ugaussian(rng) * sigma;
645 } else {
646 REAL8 min, max;
647 LALInferenceGetMinMaxPrior(thread->priorArgs, param->name, &min, &max);
648
649 *((REAL8 *)param->value) += gsl_ran_ugaussian(rng) * (max - min);
650 }
651 } else
652 *((REAL8 *)param->value) += gsl_ran_ugaussian(rng) * sigma * sqrttemp;
653
654 LALInferenceCyclicReflectiveBound(proposedParams, thread->priorArgs);
655
656 /* Set the log of the proposal ratio to zero, since this is a
657 symmetric proposal. */
658 logPropRatio = 0.0;
659
660 INT4 as = 1;
661 LALInferenceSetVariable(args, "adaptableStep", &as);
662
663 } else {
664 /* We are not adaptive, or for some reason don't have a sigma
665 vector---fall back on old proposal. */
666 logPropRatio = LALInferenceSingleProposal(thread, currentParams, proposedParams);
667 }
668
669 return logPropRatio;
670}
671
674 LALInferenceVariables *proposedParams) {
675 LALInferenceVariableItem *param=NULL;
677 gsl_rng * GSLrandom = thread->GSLrandom;
678 REAL8 sigma, big_sigma;
679 INT4 dim, varNr;
680
682
683 sigma = 0.1 * sqrt(thread->temperature); /* Adapt step to temperature. */
684 big_sigma = 1.0;
685
686 if (gsl_ran_ugaussian(GSLrandom) < 1.0e-3)
687 big_sigma = 1.0e1; //Every 1e3 iterations, take a 10x larger jump in a parameter
688 if (gsl_ran_ugaussian(GSLrandom) < 1.0e-4)
689 big_sigma = 1.0e2; //Every 1e4 iterations, take a 100x larger jump in a parameter
690
691 dim = proposedParams->dimension;
692
693 do {
694 varNr = 1 + gsl_rng_uniform_int(GSLrandom, dim);
695 param = LALInferenceGetItemNr(proposedParams, varNr);
696 } while (!LALInferenceCheckVariableNonFixed(proposedParams, param->name) || param->type != LALINFERENCE_REAL8_t);
697
698 /* Scale jumps proposal appropriately for prior sampling */
699 if (LALInferenceGetINT4Variable(args, "sampling_prior")) {
700 if (!strcmp(param->name, "eta")) {
701 sigma = 0.02;
702 } else if (!strcmp(param->name, "q")) {
703 sigma = 0.08;
704 } else if (!strcmp(param->name, "chirpmass")) {
705 sigma = 1.0;
706 } else if (!strcmp(param->name, "time")) {
707 sigma = 0.02;
708 } else if (!strcmp(param->name, "t0")) {
709 sigma = 0.02;
710 } else if (!strcmp(param->name, "phase")) {
711 sigma = 0.6;
712 } else if (!strcmp(param->name, "distance")) {
713 sigma = 10.0;
714 } else if (!strcmp(param->name, "declination")) {
715 sigma = 0.3;
716 } else if (!strcmp(param->name, "azimuth")) {
717 sigma = 0.6;
718 } else if (!strcmp(param->name, "cosalpha")) {
719 sigma = 0.1;
720 } else if (!strcmp(param->name, "rightascension")) {
721 sigma = 0.6;
722 } else if (!strcmp(param->name, "polarisation")) {
723 sigma = 0.6;
724 } else if (!strcmp(param->name,"costheta_jn")) {
725 sigma = 0.3;
726 } else if (!strcmp(param->name, "a_spin1")) {
727 sigma = 0.1;
728 } else if (!strcmp(param->name, "a_spin2")) {
729 sigma = 0.1;
730 } else {
731 fprintf(stderr, "Could not find parameter %s!", param->name);
732 exit(1);
733 }
734
735 *(REAL8 *)param->value += gsl_ran_ugaussian(GSLrandom)*sigma;
736 } else {
737 if (!strcmp(param->name,"eta") || !strcmp(param->name,"q") || !strcmp(param->name,"time") || !strcmp(param->name,"t0") || !strcmp(param->name,"a_spin2") || !strcmp(param->name,"a_spin1")){
738 *(REAL8 *)param->value += gsl_ran_ugaussian(GSLrandom)*big_sigma*sigma*0.001;
739 } else if (!strcmp(param->name,"polarisation") || !strcmp(param->name,"phase") || !strcmp(param->name,"costheta_jn")){
740 *(REAL8 *)param->value += gsl_ran_ugaussian(GSLrandom)*big_sigma*sigma*0.1;
741 } else {
742 *(REAL8 *)param->value += gsl_ran_ugaussian(GSLrandom)*big_sigma*sigma*0.01;
743 }
744 }
745
746 LALInferenceCyclicReflectiveBound(proposedParams, thread->priorArgs);
747
748 /* Symmetric Proposal. */
749 REAL8 logPropRatio = 0.0;
750
751 return logPropRatio;
752}
753
754
757 LALInferenceVariables *proposedParams) {
758 LALInferenceVariableItem *proposeIterator;
759 REAL8Vector *eigenvalues;
760 gsl_matrix *eigenvectors;
761 REAL8 jumpSize, tmp, inc;
762 REAL8 logPropRatio = 0.0;
763 INT4 N, i, j;
764
766
768 gsl_rng *rng = thread->GSLrandom;
769
770 eigenvalues = LALInferenceGetREAL8VectorVariable(args, "covarianceEigenvalues");
771 eigenvectors = LALInferenceGetgslMatrixVariable(args, "covarianceEigenvectors");
772
773 N = eigenvalues->length;
774
775 i = gsl_rng_uniform_int(rng, N);
776 jumpSize = sqrt(thread->temperature * eigenvalues->data[i]) * gsl_ran_ugaussian(rng);
777
778 j = 0;
779 proposeIterator = proposedParams->head;
780 if (proposeIterator == NULL) {
781 fprintf(stderr, "Bad proposed params in %s, line %d\n",
782 __FILE__, __LINE__);
783 exit(1);
784 }
785
786 do {
787 if (LALInferenceCheckVariableNonFixed(proposedParams, proposeIterator->name) &&
788 proposeIterator->type==LALINFERENCE_REAL8_t) {
789 tmp = LALInferenceGetREAL8Variable(proposedParams, proposeIterator->name);
790 inc = jumpSize * gsl_matrix_get(eigenvectors, j, i);
791
792 tmp += inc;
793
794 LALInferenceSetVariable(proposedParams, proposeIterator->name, &tmp);
795
796 j++;
797 }
798 } while ((proposeIterator = proposeIterator->next) != NULL && j < N);
799
800 return logPropRatio;
801}
802
805 LALInferenceVariables *proposedParams) {
806 REAL8 sigma;
807 REAL8 jumpX, jumpY;
808 REAL8 RA, DEC;
809 REAL8 newRA, newDEC;
810 REAL8 one_deg = 1.0 / (2.0*M_PI);
811 REAL8 logPropRatio = 0.0;
812
814
815 gsl_rng *rng = thread->GSLrandom;
816
817 sigma = sqrt(thread->temperature) * one_deg;
818 jumpX = sigma * gsl_ran_ugaussian(rng);
819 jumpY = sigma * gsl_ran_ugaussian(rng);
820
821 RA = LALInferenceGetREAL8Variable(proposedParams, "rightascension");
822 DEC = LALInferenceGetREAL8Variable(proposedParams, "declination");
823
824 newRA = RA + jumpX;
825 newDEC = DEC + jumpY;
826
827 LALInferenceSetVariable(proposedParams, "rightascension", &newRA);
828 LALInferenceSetVariable(proposedParams, "declination", &newDEC);
829
830
831 return logPropRatio;
832}
833
836 LALInferenceVariables *proposedParams) {
837 return(LALInferenceDifferentialEvolutionNames(thread, currentParams, proposedParams, NULL));
838}
839
842 LALInferenceVariables *proposedParams) {
843 return(LALInferenceEnsembleStretchNames(thread, currentParams, proposedParams, NULL));
844}
845
846
849 LALInferenceVariables *proposedParams) {
850
851 return(LALInferenceEnsembleStretchNames(thread, currentParams, proposedParams, intrinsicNames));
852
853}
854
857 LALInferenceVariables *proposedParams) {
858 REAL8 logPropRatio;
859
860 logPropRatio = LALInferenceEnsembleStretchNames(thread, currentParams, proposedParams, extrinsicNames);
861
862 return logPropRatio;
863}
864
865/* This jump uses the current sample 'A' and another randomly
866 * drawn 'B' from the ensemble of live points, and proposes
867 * C = B+Z(A-B) where Z is a scale factor */
870 LALInferenceVariables *proposedParams,
871 const char **names) {
872 size_t i, N, Ndim, nPts;
873 REAL8 logPropRatio;
874 REAL8 maxScale, Y, logmax, X, scale;
875 REAL8 cur, other, x;
877 LALInferenceVariables **dePts;
879
880 N = LALInferenceGetVariableDimension(currentParams) + 1; /* More names than we need. */
881 const char* local_names[N];
882 if (names == NULL) {
883 names = local_names;
884
885 item = currentParams->head;
886 i = 0;
887 while (item != NULL) {
889 names[i] = item->name;
890 i++;
891 }
892 item = item->next;
893 }
894 names[i]=NULL; /* Terminate */
895 }
896
897 Ndim = 0;
898 for(Ndim=0,i=0; names[i] != NULL; i++ ) {
900 Ndim++;
901 }
902
904
905 Ndim = 0;
906 for(Ndim=0, i=0; names[i] != NULL; i++ ) {
907 if (LALInferenceCheckVariableNonFixed(proposedParams, names[i]))
908 Ndim++;
909 }
910
911 dePts = thread->differentialPoints;
912 nPts = thread->differentialPointsLength;
913
914 if (dePts == NULL || nPts <= 1) {
915 logPropRatio = 0.0;
916 return logPropRatio; /* Quit now, since we don't have any points to use. */
917 }
918
919 /* Choose a different sample */
920 do {
921 i = gsl_rng_uniform_int(thread->GSLrandom, nPts);
922 } while (!LALInferenceCompareVariables(currentParams, dePts[i]));
923
924 ptI = dePts[i];
925
926 /* Scale z is chosen according to be symmetric under z -> 1/z */
927 /* so p(x) \propto 1/z between 1/a and a */
928
929 /* TUNABLE PARAMETER (a), must be >1. Larger value -> smaller acceptance */
930 maxScale=3.0;
931
932 /* Draw sample between 1/max and max */
933 Y = gsl_rng_uniform(thread->GSLrandom);
934 logmax = log(maxScale);
935 X = 2.0*logmax*Y - logmax;
936 scale = exp(X);
937
938 for (i = 0; names[i] != NULL; i++) {
939 /* Ignore variable if it's not in each of the params. */
940 if (LALInferenceCheckVariableNonFixed(proposedParams, names[i]) &&
942 cur = LALInferenceGetREAL8Variable(proposedParams, names[i]);
944 x = other + scale*(cur-other);
945
946 LALInferenceSetVariable(proposedParams, names[i], &x);
947 }
948 }
949
950 if (scale < maxScale && scale > (1.0/maxScale))
951 logPropRatio = log(scale)*((REAL8)Ndim);
952 else
953 logPropRatio = -INFINITY;
954
955 return logPropRatio;
956}
957
958
961 LALInferenceVariables *proposedParams) {
962 return(LALInferenceEnsembleWalkNames(thread, currentParams, proposedParams, NULL));
963}
964
965
968 LALInferenceVariables *proposedParams) {
969
970 return(LALInferenceEnsembleWalkNames(thread, currentParams, proposedParams, intrinsicNames));
971
972}
973
976 LALInferenceVariables *proposedParams) {
977 return(LALInferenceEnsembleWalkNames(thread, currentParams, proposedParams, extrinsicNames));
978}
979
980
983 LALInferenceVariables *proposedParams,
984 const char **names) {
985 size_t i;
987 REAL8 logPropRatio = 0.0;
988
989 size_t N = LALInferenceGetVariableDimension(currentParams) + 1; /* More names than we need. */
990 const char* local_names[N];
991 if (names == NULL) {
992 names = local_names;
993
994 item = currentParams->head;
995 i = 0;
996 while (item != NULL) {
998 names[i] = item->name;
999 i++;
1000 }
1001
1002 item = item->next;
1003 }
1004 names[i]=NULL; /* Terminate */
1005 }
1006
1007
1008 /*
1009 size_t Ndim = 0;
1010 for(Ndim=0,i=0; names[i] != NULL; i++ ) {
1011 if(LALInferenceCheckVariableNonFixed(currentParams,names[i]))
1012 Ndim++;
1013 }
1014 */
1015
1016 LALInferenceVariables **pointsPool = thread->differentialPoints;
1017 size_t k=0;
1018 size_t sample_size=3;
1019
1020 LALInferenceVariables **dePts = thread->differentialPoints;
1021 size_t nPts = thread->differentialPointsLength;
1022
1023 if (dePts == NULL || nPts <= 1) {
1024 logPropRatio = 0.0;
1025 return logPropRatio; /* Quit now, since we don't have any points to use. */
1026 }
1027
1029
1030 UINT4 indices[sample_size];
1031 UINT4 all_indices[nPts];
1032
1033 for (i=0;i<nPts;i++) all_indices[i]=i;
1034 gsl_ran_choose(thread->GSLrandom,indices, sample_size, all_indices, nPts, sizeof(UINT4));
1035
1036 double w=0.0;
1037 double univariate_normals[sample_size];
1038 for(i=0;i<sample_size;i++) univariate_normals[i] = gsl_ran_ugaussian(thread->GSLrandom);
1039
1040 /* Note: Simplified this loop on master 2015-08-12, take this version when rebasing */
1041 for(k=0;names[k]!=NULL;k++)
1042 {
1043 if(!LALInferenceCheckVariableNonFixed(proposedParams,names[k]) || LALInferenceGetVariableType(proposedParams,names[k])!=LALINFERENCE_REAL8_t) continue;
1044 REAL8 centre_of_mass=0.0;
1045 /* Compute centre of mass */
1046 for(i=0;i<sample_size;i++)
1047 {
1048 centre_of_mass+=LALInferenceGetREAL8Variable(pointsPool[indices[i]],names[k])/((REAL8)sample_size);
1049 }
1050 /* Compute offset */
1051 for(i=0,w=0.0;i<sample_size;i++)
1052 {
1053 w+= univariate_normals[i] * (LALInferenceGetREAL8Variable(pointsPool[indices[i]],names[k]) - centre_of_mass);
1054 }
1055 REAL8 tmp = LALInferenceGetREAL8Variable(proposedParams,names[k]) + w;
1056 LALInferenceSetVariable(proposedParams,names[k],&tmp);
1057 }
1058
1059 logPropRatio = 0.0;
1060
1061 return logPropRatio;
1062}
1063
1066 LALInferenceVariables *proposedParams,
1067 const char **names) {
1068 size_t i, j, N, Ndim, nPts;
1070 LALInferenceVariables **dePts;
1071 LALInferenceVariables *ptI, *ptJ;
1072 REAL8 logPropRatio = 0.0;
1073 REAL8 scale, x;
1074 N = LALInferenceGetVariableDimension(currentParams) + 1; /* More names than we need. */
1075 const char *localnames[N];
1076
1077 gsl_rng *rng = thread->GSLrandom;
1078
1079 if (names == NULL) {
1080 item = currentParams->head;
1081 i = 0;
1082 while (item != NULL) {
1084 localnames[i] = item->name;
1085 i++;
1086 }
1087
1088 item = item->next;
1089 }
1090 localnames[i]=NULL; /* Terminate */
1091 names = localnames;
1092 }
1093
1094 Ndim = 0;
1095 for (Ndim=0, i=0; names[i] != NULL; i++ ) {
1097 Ndim++;
1098 }
1099
1100 dePts = thread->differentialPoints;
1101 nPts = thread->differentialPointsLength;
1102
1103 if (dePts == NULL || nPts <= 1)
1104 return logPropRatio; /* Quit now, since we don't have any points to use. */
1105
1107
1108
1109 i = gsl_rng_uniform_int(rng, nPts);
1110 do {
1111 j = gsl_rng_uniform_int(rng, nPts);
1112 } while (j == i);
1113
1114 ptI = dePts[i];
1115 ptJ = dePts[j];
1116
1117 const REAL8 modeHoppingFrac = 0.5;
1118 /* Some fraction of the time, we do a "mode hopping" jump,
1119 where we jump exactly along the difference vector. */
1120 if (gsl_rng_uniform(rng) < modeHoppingFrac) {
1121 scale = 1.0;
1122 } else {
1123 /* Otherwise scale is chosen uniform in log between 0.1 and 10 times the
1124 desired jump size. */
1125 scale = 2.38/sqrt(Ndim) * exp(log(0.1) + log(100.0) * gsl_rng_uniform(rng));
1126 }
1127
1128 for (i = 0; names[i] != NULL; i++) {
1132 /* Ignore variable if it's not in each of the params. */
1133 } else {
1135 x += scale * LALInferenceGetREAL8Variable(ptJ, names[i]);
1136 x -= scale * LALInferenceGetREAL8Variable(ptI, names[i]);
1137 LALInferenceSetVariable(proposedParams, names[i], &x);
1138 }
1139 }
1140
1141 return logPropRatio;
1142}
1143
1146 LALInferenceVariables *proposedParams) {
1147
1149}
1150
1153 LALInferenceVariables *proposedParams) {
1155}
1156
1158 REAL8 dmin, dmax, x;
1159
1160 LALInferenceGetMinMaxPrior(thread->priorArgs, "distance", &dmin, &dmax);
1161
1162 x = gsl_rng_uniform(thread->GSLrandom);
1163
1164 return cbrt(x*(dmax*dmax*dmax - dmin*dmin*dmin) + dmin*dmin*dmin);
1165}
1166
1168 REAL8 logdmin, logdmax;
1169
1170 LALInferenceGetMinMaxPrior(thread->priorArgs, "logdistance", &logdmin, &logdmax);
1171
1172 REAL8 dmin=exp(logdmin);
1173 REAL8 dmax=exp(logdmax);
1174
1175 REAL8 x = gsl_rng_uniform(thread->GSLrandom);
1176
1177 return log(cbrt(x*(dmax*dmax*dmax - dmin*dmin*dmin) + dmin*dmin*dmin));
1178}
1179
1180static REAL8 draw_colatitude(LALInferenceThreadState *thread, const char *name) {
1181 REAL8 min, max, x;
1182
1183 LALInferenceGetMinMaxPrior(thread->priorArgs, name, &min, &max);
1184
1185 x = gsl_rng_uniform(thread->GSLrandom);
1186
1187 return acos(cos(min) - x*(cos(min) - cos(max)));
1188}
1189
1191 REAL8 min, max, x;
1192
1193 LALInferenceGetMinMaxPrior(thread->priorArgs, "declination", &min, &max);
1194
1195 x = gsl_rng_uniform(thread->GSLrandom);
1196
1197 return asin(x*(sin(max) - sin(min)) + sin(min));
1198}
1199
1200static REAL8 draw_flat(LALInferenceThreadState *thread, const char *name) {
1201 REAL8 min, max, x;
1202
1203 LALInferenceGetMinMaxPrior(thread->priorArgs, name, &min, &max);
1204
1205 x = gsl_rng_uniform(thread->GSLrandom);
1206
1207 return min + x*(max - min);
1208}
1209
1211 REAL8 min, max, delta;
1212 REAL8 mMin56, mMax56, u;
1213
1214 LALInferenceGetMinMaxPrior(thread->priorArgs, "chirpmass", &min, &max);
1215
1216 mMin56 = pow(min, 5.0/6.0);
1217 mMax56 = pow(max, 5.0/6.0);
1218
1219 delta = 1.0/mMin56 - 1.0/mMax56;
1220
1221 u = delta*gsl_rng_uniform(thread->GSLrandom);
1222
1223 return pow(1.0/(1.0/mMin56 - u), 6.0/5.0);
1224}
1225
1227 REAL8 logP = 0.0;
1228
1229 if(LALInferenceCheckVariable(params, "chirpmass")) {
1230 REAL8 Mc = *(REAL8 *)LALInferenceGetVariable(params, "chirpmass");
1231 logP += -11.0/6.0*log(Mc);
1232 }
1233
1234 /* Flat in time, ra, psi, phi. */
1235
1236 if(LALInferenceCheckVariable(params,"logdistance"))
1237 logP += 3.0* *(REAL8 *)LALInferenceGetVariable(params,"logdistance");
1238 else if(LALInferenceCheckVariable(params,"distance"))
1239 logP += 2.0*log(*(REAL8 *)LALInferenceGetVariable(params,"distance"));
1240
1241 if(LALInferenceCheckVariable(params,"declination"))
1242 logP += log(cos(*(REAL8 *)LALInferenceGetVariable(params, "declination")));
1243
1244 if (LALInferenceCheckVariable(params, "tilt_spin1")) {
1245 logP += log(fabs(sin(*(REAL8 *)LALInferenceGetVariable(params, "tilt_spin1"))));
1246 }
1247
1248 if (LALInferenceCheckVariable(params, "tilt_spin2")) {
1249 logP += log(fabs(sin(*(REAL8 *)LALInferenceGetVariable(params, "tilt_spin2"))));
1250 }
1251
1252 return logP;
1253}
1254
1255/* WARNING: If you add any non-flat draws to this proposal, you MUST
1256 update the above approxLogPrior function with the corresponding
1257 density. The proposal ratio is calculated using approxLogPrior, so
1258 non-flat proposals that do not have a corresponding density term in
1259 approxLogPrior will result in a failure to sample from the
1260 posterior density! */
1263 LALInferenceVariables *proposedParams) {
1264 REAL8 tmp = 0.0;
1265 INT4 analytic_test, i;
1266 REAL8 logBackwardJump;
1267 REAL8 logPropRatio;
1269
1271
1272 const char *flat_params[] = {"q", "eta", "t0", "azimuth", "cosalpha", "time", "phase", "polarisation",
1273 "rightascension", "costheta_jn", "phi_jl",
1274 "phi12", "a_spin1", "a_spin2", "logp1", "gamma1", "gamma2", "gamma3",
1275 "SDgamma0","SDgamma1","SDgamma2","SDgamma3", NULL};
1276
1278
1279 analytic_test = LALInferenceGetINT4Variable(args, "analytical_test");
1280
1281 if (analytic_test) {
1283 while (ptr!=NULL) {
1285 tmp = draw_flat(thread, ptr->name);
1286 LALInferenceSetVariable(proposedParams, ptr->name, &tmp);
1287 }
1288 ptr=ptr->next;
1289 }
1290 } else {
1291 logBackwardJump = approxLogPrior(currentParams);
1292
1293 for (i = 0; flat_params[i] != NULL; i++) {
1294 if (LALInferenceCheckVariableNonFixed(proposedParams, flat_params[i])) {
1295 REAL8 val = draw_flat(thread, flat_params[i]);
1296 LALInferenceSetVariable(proposedParams, flat_params[i], &val);
1297 }
1298 }
1299
1300 if (LALInferenceCheckVariableNonFixed(proposedParams, "chirpmass")) {
1301 REAL8 Mc = draw_chirp(thread);
1302 LALInferenceSetVariable(proposedParams, "chirpmass", &Mc);
1303 }
1304
1305 if (LALInferenceCheckVariableNonFixed(proposedParams, "logdistance")) {
1306 REAL8 logdist = draw_logdistance(thread);
1307 LALInferenceSetVariable(proposedParams, "logdistance", &logdist);
1308 } else if (LALInferenceCheckVariableNonFixed(proposedParams, "distance")) {
1309 REAL8 dist = draw_distance(thread);
1310 LALInferenceSetVariable(proposedParams, "distance", &dist);
1311 }
1312
1313 if (LALInferenceCheckVariableNonFixed(proposedParams, "declination")) {
1314 REAL8 dec = draw_dec(thread);
1315 LALInferenceSetVariable(proposedParams, "declination", &dec);
1316 }
1317
1318 if (LALInferenceCheckVariableNonFixed(proposedParams, "tilt_spin1")) {
1319 REAL8 tilt1 = draw_colatitude(thread, "tilt_spin1");
1320 LALInferenceSetVariable(proposedParams, "tilt_spin1", &tilt1);
1321 }
1322
1323 if (LALInferenceCheckVariableNonFixed(proposedParams, "tilt_spin2")) {
1324 REAL8 tilt2 = draw_colatitude(thread, "tilt_spin2");
1325 LALInferenceSetVariable(proposedParams, "tilt_spin2", &tilt2);
1326 }
1327
1328 if (LALInferenceCheckVariableNonFixed(proposedParams, "psdscale")) {
1329 REAL8 x, min, max;
1330 INT4 j;
1331
1332 min=0.10;
1333 max=10.0;
1334
1335 gsl_matrix *eta = LALInferenceGetgslMatrixVariable(proposedParams, "psdscale");
1336
1337 for(i=0; i<(INT8)eta->size1; i++) {
1338 for(j=0; j<(INT8)eta->size2; j++) {
1339 x = min + gsl_rng_uniform(thread->GSLrandom) * (max - min);
1340 gsl_matrix_set(eta, i, j, x);
1341 }
1342 }
1343 }//end if(psdscale)
1344 }
1345
1346 if (analytic_test) {
1347 /* Flat in every variable means uniform jump probability. */
1348 logPropRatio = 0.0;
1349 } else {
1350 logPropRatio = logBackwardJump - approxLogPrior(proposedParams);
1351 }
1352
1353 return logPropRatio;
1354}
1355
1357 REAL8 logPropRatio = 0., tmp = 0.;
1358
1361
1362 while(ptr!=NULL) {
1363 if(ptr->vary != LALINFERENCE_PARAM_FIXED && LALInferenceCheckMinMaxPrior(threadState->priorArgs, ptr->name ) ) {
1364 tmp = draw_flat(threadState, ptr->name);
1365 LALInferenceSetVariable(proposedParams, ptr->name, &tmp);
1366 }
1367 ptr=ptr->next;
1368 }
1369
1370 return logPropRatio;
1371}
1372
1373static void cross_product(REAL8 x[3], const REAL8 y[3], const REAL8 z[3]) {
1374 x[0] = y[1]*z[2]-y[2]*z[1];
1375 x[1] = y[2]*z[0]-y[0]*z[2];
1376 x[2] = y[0]*z[1]-y[1]*z[0];
1377}
1378
1379static REAL8 norm(const REAL8 x[3]) {
1380 return sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
1381}
1382
1383static void unit_vector(REAL8 v[3], const REAL8 w[3]) {
1384 REAL8 n = norm(w);
1385
1386 if (n == 0.0) {
1387 XLALError("unit_vector", __FILE__, __LINE__, XLAL_FAILURE);
1388 exit(1);
1389 } else {
1390 v[0] = w[0] / n;
1391 v[1] = w[1] / n;
1392 v[2] = w[2] / n;
1393 }
1394}
1395
1396static REAL8 dot(const REAL8 v[3], const REAL8 w[3]) {
1397 return v[0]*w[0] + v[1]*w[1] + v[2]*w[2];
1398}
1399
1400static void project_along(REAL8 vproj[3], const REAL8 v[3], const REAL8 w[3]) {
1401 REAL8 what[3];
1402 REAL8 vdotw;
1403
1404 unit_vector(what, w);
1405 vdotw = dot(v, w);
1406
1407 vproj[0] = what[0]*vdotw;
1408 vproj[1] = what[1]*vdotw;
1409 vproj[2] = what[2]*vdotw;
1410}
1411
1412static void vsub(REAL8 diff[3], const REAL8 w[3], const REAL8 v[3]) {
1413 diff[0] = w[0] - v[0];
1414 diff[1] = w[1] - v[1];
1415 diff[2] = w[2] - v[2];
1416}
1417
1418static void reflect_plane(REAL8 pref[3], const REAL8 p[3],
1419 const REAL8 x[3], const REAL8 y[3], const REAL8 z[3]) {
1420 REAL8 n[3], nhat[3], xy[3], xz[3], pn[3], pnperp[3];
1421
1422 vsub(xy, y, x);
1423 vsub(xz, z, x);
1424
1425 cross_product(n, xy, xz);
1426 unit_vector(nhat, n);
1427
1428 project_along(pn, p, nhat);
1429 vsub(pnperp, p, pn);
1430
1431 vsub(pref, pnperp, pn);
1432}
1433
1434static void get_detectors(LALInferenceIFOData *data, LALDetector **detectors) {
1435 INT4 i;
1436 INT4 nDet = 0;
1438 for (i=0,ifo=data; ifo; i++,ifo=ifo->next)
1439 nDet++;
1440
1441 *detectors = XLALCalloc(nDet, sizeof(LALDetector));
1442 for (i=0,ifo=data; ifo; i++,ifo=ifo->next)
1443 *detectors[i] = *(ifo->detector);
1444}
1445
1446static void sph_to_cart(REAL8 cart[3], const REAL8 lat, const REAL8 longi) {
1447 cart[0] = cos(longi)*cos(lat);
1448 cart[1] = sin(longi)*cos(lat);
1449 cart[2] = sin(lat);
1450}
1451
1452static void cart_to_sph(const REAL8 cart[3], REAL8 *lat, REAL8 *longi) {
1453 *longi = atan2(cart[1], cart[0]);
1454 *lat = asin(cart[2] / sqrt(cart[0]*cart[0] + cart[1]*cart[1] + cart[2]*cart[2]));
1455}
1456
1457static void reflected_position_and_time(LALInferenceThreadState *thread, const REAL8 ra, const REAL8 dec,
1458 const REAL8 oldTime, REAL8 *newRA, REAL8 *newDec, REAL8 *newTime) {
1460 memset(&status, 0, sizeof(status));
1461 SkyPosition currentEqu, currentGeo, newEqu, newGeo;
1464 REAL8 x[3], y[3], z[3];
1465 REAL8 currentLoc[3], newLoc[3];
1466 REAL8 newGeoLat, newGeoLongi;
1467 REAL8 oldDt, newDt;
1468 LALDetector xD, yD, zD;
1469
1470 currentEqu.latitude = dec;
1471 currentEqu.longitude = ra;
1474
1475 epoch = thread->parent->data->epoch;
1476 get_detectors(thread->parent->data, &detectors);
1477
1478 LALEquatorialToGeographic(&status, &currentGeo, &currentEqu, &epoch);
1479
1480 /* This function should only be called when we know that we have
1481 three detectors, or the following will crash. */
1482 xD = detectors[0];
1483 memcpy(x, xD.location, 3*sizeof(REAL8));
1484
1485 INT4 det = 1;
1486 yD = detectors[det];
1487 while (same_detector_location(&yD, &xD)) {
1488 det++;
1489 yD = detectors[det];
1490 }
1491 memcpy(y, yD.location, 3*sizeof(REAL8));
1492 det++;
1493
1494 zD = detectors[det];
1495 while (same_detector_location(&zD, &yD) || same_detector_location(&zD, &xD)) {
1496 det++;
1497 zD = detectors[det];
1498 }
1499 memcpy(z, zD.location, 3*sizeof(REAL8));
1500
1501 sph_to_cart(currentLoc, currentGeo.latitude, currentGeo.longitude);
1502
1503 reflect_plane(newLoc, currentLoc, x, y, z);
1504
1505 cart_to_sph(newLoc, &newGeoLat, &newGeoLongi);
1506
1507 newGeo.latitude = newGeoLat;
1508 newGeo.longitude = newGeoLongi;
1511 LALGeographicToEquatorial(&status, &newEqu, &newGeo, &epoch);
1512
1513 oldDt = XLALTimeDelayFromEarthCenter(detectors[0].location, currentEqu.longitude,
1514 currentEqu.latitude, &epoch);
1515 newDt = XLALTimeDelayFromEarthCenter(detectors[0].location, newEqu.longitude,
1516 newEqu.latitude, &epoch);
1517
1518 *newRA = newEqu.longitude;
1519 *newDec = newEqu.latitude;
1520 *newTime = oldTime + oldDt - newDt;
1521
1523}
1524
1526 LALInferenceVariables *proposedParams,
1527 INT4 ifo, INT4 k) {
1528 REAL8 prior = 0.0;
1529 REAL8 component_min,component_max;
1530 REAL8 A, f, Q, Anorm;
1531 gsl_matrix *glitch_f, *glitch_Q, *glitch_A;
1532
1533 component_min = LALInferenceGetREAL8Variable(thread->priorArgs,"morlet_f0_prior_min");
1534 component_max = LALInferenceGetREAL8Variable(thread->priorArgs,"morlet_f0_prior_max");
1535 prior -= log(component_max - component_min);
1536
1537 component_min = LALInferenceGetREAL8Variable(thread->priorArgs, "morlet_Q_prior_min");
1538 component_max = LALInferenceGetREAL8Variable(thread->priorArgs, "morlet_Q_prior_max");
1539 prior -= log(component_max - component_min);
1540
1541 component_min = LALInferenceGetREAL8Variable(thread->priorArgs, "morlet_t0_prior_min");
1542 component_max = LALInferenceGetREAL8Variable(thread->priorArgs, "morlet_t0_prior_max");
1543 prior -= log(component_max - component_min);
1544
1545 component_min = LALInferenceGetREAL8Variable(thread->priorArgs, "morlet_phi_prior_min");
1546 component_max = LALInferenceGetREAL8Variable(thread->priorArgs, "morlet_phi_prior_max");
1547 prior -= log(component_max - component_min);
1548
1549 //"Malmquist" prior on A
1550 glitch_f = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_f0");
1551 glitch_Q = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_Q");
1552 glitch_A = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_Amp");
1553
1554 A = gsl_matrix_get(glitch_A, ifo, k);
1555 Q = gsl_matrix_get(glitch_Q, ifo, k);
1556 f = gsl_matrix_get(glitch_f, ifo, k);
1557
1558 Anorm = LALInferenceGetREAL8Variable(thread->priorArgs, "glitch_norm");
1559
1560 prior += logGlitchAmplitudeDensity(A*Anorm, Q, f);
1561
1562 return prior;
1563}
1564
1565
1566static REAL8 glitchAmplitudeDraw(REAL8 Q, REAL8 f, gsl_rng *r) {
1567 REAL8 SNR;
1568 REAL8 PIterm = 0.5*LAL_2_SQRTPI*LAL_SQRT1_2;
1569 REAL8 SNRPEAK = 5.0;
1570
1571 REAL8 den=0.0, alpha=1.0;
1572 REAL8 max= 1.0/(SNRPEAK*LAL_E);;
1573
1574 // x/a^2 exp(-x/a) prior on SNR. Peaks at x = a. Good choice is a=5
1575
1576 // rejection sample. Envelope function on runs out to ten times the peak
1577 // don't even bother putting in this minute correction to the normalization
1578 // (it is a 5e-4 correction).
1579 do {
1580 SNR = 20.0 * SNRPEAK * gsl_rng_uniform(r);
1581
1582 den = SNR/(SNRPEAK*SNRPEAK) * exp(-SNR/SNRPEAK);
1583
1584 den /= max;
1585
1586 alpha = gsl_rng_uniform(r);
1587 } while (alpha > den);
1588
1589 return SNR/sqrt((PIterm*Q/f));
1590}
1591
1594 LALInferenceVariables *proposedParams) {
1595 INT4 i, j, l;
1596 INT4 ifo, nifo, timeflag=0;
1597 REAL8 logPropRatio = 0.0;
1598 REAL8 ra, dec;
1599 REAL8 baryTime, gmst;
1600 REAL8 newRA, newDec, newTime, newPsi;
1601 REAL8 intpart, decpart;
1602 REAL8 omega, cosomega, sinomega, c1momega;
1603 REAL8 IFO1[3], IFO2[3];
1604 REAL8 IFOX[3], k[3];
1605 REAL8 normalize;
1606 REAL8 pForward, pReverse;
1607 REAL8 n[3];
1608 REAL8 kp[3];
1609 LIGOTimeGPS GPSlal, epoch;
1611 gsl_matrix *IFO;
1612
1614
1616 gsl_rng *rng = thread->GSLrandom;
1617
1618 epoch = thread->parent->data->epoch;
1619 get_detectors(thread->parent->data, &detectors);
1620
1621 ra = LALInferenceGetREAL8Variable(proposedParams, "rightascension");
1622 dec = LALInferenceGetREAL8Variable(proposedParams, "declination");
1623
1624 if (LALInferenceCheckVariable(proposedParams, "time")){
1625 baryTime = LALInferenceGetREAL8Variable(proposedParams, "time");
1626 timeflag = 1;
1627 } else {
1628 baryTime = XLALGPSGetREAL8(&epoch);
1629 }
1630
1631 XLALGPSSetREAL8(&GPSlal, baryTime);
1632 gmst = XLALGreenwichMeanSiderealTime(&GPSlal);
1633
1634 //remap gmst back to [0:2pi]
1635 gmst /= LAL_TWOPI;
1636 intpart = (INT4)gmst;
1637 decpart = gmst - (REAL8)intpart;
1638 gmst = decpart*LAL_TWOPI;
1639 gmst = gmst < 0. ? gmst + LAL_TWOPI : gmst;
1640
1641 /*
1642 line-of-sight vector
1643 */
1644 k[0] = cos(gmst-ra) * cos(dec);
1645 k[1] =-sin(gmst-ra) * cos(dec);
1646 k[2] = sin(dec);
1647
1648 /*
1649 Store location for each detector
1650 */
1651 nifo = LALInferenceGetINT4Variable(args, "nDet");
1652
1653 IFO = gsl_matrix_alloc(nifo, 3);
1654
1655 for(ifo=0; ifo<nifo; ifo++) {
1656 memcpy(IFOX, detectors[ifo].location, 3*sizeof(REAL8));
1657 for (i=0; i<3; i++)
1658 gsl_matrix_set(IFO, ifo, i, IFOX[i]);
1659 }
1660
1661 /*
1662 Randomly select two detectors from the network
1663 -this assumes there are no co-located detectors
1664 */
1665 i = j = 0;
1666 while (i==j) {
1667 i=gsl_rng_uniform_int(rng, nifo);
1668 j=gsl_rng_uniform_int(rng, nifo);
1669 }
1670
1671 for(l=0; l<3; l++) {
1672 IFO1[l] = gsl_matrix_get(IFO, i, l);
1673 IFO2[l] = gsl_matrix_get(IFO, j, l);
1674 }
1675
1676 /*
1677 detector axis
1678 */
1679 normalize=0.0;
1680 for(i=0; i<3; i++) {
1681 n[i] = IFO1[i] - IFO2[i];
1682 normalize += n[i] * n[i];
1683 }
1684 normalize = 1./sqrt(normalize);
1685 for(i=0; i<3; i++)
1686 n[i] *= normalize;
1687
1688 /*
1689 rotation angle
1690 */
1691 omega = LAL_TWOPI * gsl_rng_uniform(rng);
1692 cosomega = cos(omega);
1693 sinomega = sin(omega);
1694 c1momega = 1.0 - cosomega;
1695
1696 /*
1697 rotate k' = Rk
1698 */
1699 kp[0] = (c1momega*n[0]*n[0] + cosomega) * k[0]
1700 + (c1momega*n[0]*n[1] - sinomega*n[2]) * k[1]
1701 + (c1momega*n[0]*n[2] + sinomega*n[1]) * k[2];
1702 kp[1] = (c1momega*n[0]*n[1] + sinomega*n[2]) * k[0]
1703 + (c1momega*n[1]*n[1] + cosomega) * k[1]
1704 + (c1momega*n[1]*n[2] - sinomega*n[0]) * k[2];
1705 kp[2] = (c1momega*n[0]*n[2] - sinomega*n[1]) * k[0]
1706 + (c1momega*n[1]*n[2] + sinomega*n[0]) * k[1]
1707 + (c1momega*n[2]*n[2] + cosomega) * k[2];
1708
1709 /*
1710 convert k' back to ra' and dec'
1711 */
1712 newDec = asin(kp[2]);
1713 newRA = atan2(kp[1], kp[0]) + gmst;
1714 if (newRA < 0.0)
1715 newRA += LAL_TWOPI;
1716 else if (newRA >= LAL_TWOPI)
1717 newRA -= LAL_TWOPI;
1718 /*
1719 compute new geocenter time using
1720 fixed arrival time at IFO1 (arbitrary)
1721 */
1722 REAL8 tx; //old time shift = k * n
1723 REAL8 ty; //new time shift = k'* n
1724 tx=ty=0;
1725 for(i=0; i<3; i++) {
1726 tx += -IFO1[i]*k[i] /LAL_C_SI;
1727 ty += -IFO1[i]*kp[i]/LAL_C_SI;
1728 }
1729 newTime = tx + baryTime - ty;
1730
1731 XLALGPSSetREAL8(&GPSlal, newTime);
1732
1733 /*
1734 draw new polarisation angle uniformally
1735 for now
1736 MARK: Need to be smarter about psi in sky-ring jump
1737 */
1738 newPsi = LAL_PI * gsl_rng_uniform(rng);
1739
1740 LALInferenceSetVariable(proposedParams, "polarisation", &newPsi);
1741 LALInferenceSetVariable(proposedParams, "rightascension", &newRA);
1742 LALInferenceSetVariable(proposedParams, "declination", &newDec);
1743 if (timeflag)
1744 LALInferenceSetVariable(proposedParams, "time", &newTime);
1745
1746 pForward = cos(newDec);
1747 pReverse = cos(dec);
1748
1749 gsl_matrix_free(IFO);
1750
1751 logPropRatio = log(pReverse/pForward);
1752
1754
1755 return logPropRatio;
1756}
1757
1760 LALInferenceVariables *proposedParams) {
1761 INT4 timeflag=0;
1762 REAL8 ra, dec, baryTime;
1763 REAL8 newRA, newDec, newTime;
1764 REAL8 nRA, nDec, nTime;
1765 REAL8 refRA, refDec, refTime;
1766 REAL8 nRefRA, nRefDec, nRefTime;
1767 REAL8 pForward, pReverse;
1768 REAL8 logPropRatio = 0.0;
1769 INT4 nUniqueDet;
1771
1772
1773 /* Find the number of distinct-position detectors. */
1774 /* Exit with same parameters (with a warning the first time) if
1775 there are not three detectors. */
1776 static INT4 warningDelivered = 0;
1777
1779 gsl_rng *rng = thread->GSLrandom;
1780
1781 epoch = thread->parent->data->epoch;
1782 nUniqueDet = LALInferenceGetINT4Variable(args, "nUniqueDet");
1783
1784 if (nUniqueDet != 3) {
1785 if (!warningDelivered) {
1786 fprintf(stderr, "WARNING: trying to reflect through the decector plane with %d\n", nUniqueDet);
1787 fprintf(stderr, "WARNING: geometrically independent locations,\n");
1788 fprintf(stderr, "WARNING: but this proposal should only be used with exactly 3 independent detectors.\n");
1789 fprintf(stderr, "WARNING: %s, line %d\n", __FILE__, __LINE__);
1790 warningDelivered = 1;
1791 }
1792
1793 return logPropRatio;
1794 }
1796
1797 ra = LALInferenceGetREAL8Variable(currentParams, "rightascension");
1798 dec = LALInferenceGetREAL8Variable(currentParams, "declination");
1799
1801 baryTime = LALInferenceGetREAL8Variable(currentParams, "time");
1802 timeflag=1;
1803 } else {
1804 baryTime = XLALGPSGetREAL8(&epoch);
1805 }
1806
1807 reflected_position_and_time(thread, ra, dec, baryTime, &newRA, &newDec, &newTime);
1808
1809 /* Unit normal deviates, used to "fuzz" the state. */
1810 const REAL8 epsTime = 6e-6; /* 1e-1 / (16 kHz) */
1811 const REAL8 epsAngle = 3e-4; /* epsTime*c/R_Earth */
1812
1813 nRA = gsl_ran_ugaussian(rng);
1814 nDec = gsl_ran_ugaussian(rng);
1815 nTime = gsl_ran_ugaussian(rng);
1816
1817 newRA += epsAngle*nRA;
1818 newDec += epsAngle*nDec;
1819 newTime += epsTime*nTime;
1820
1821 /* And the doubly-reflected position (near the original, but not
1822 exactly due to the fuzzing). */
1823 reflected_position_and_time(thread, newRA, newDec, newTime, &refRA, &refDec, &refTime);
1824
1825 /* The Gaussian increments required to shift us back to the original
1826 position from the doubly-reflected position. */
1827 nRefRA = (ra - refRA)/epsAngle;
1828 nRefDec = (dec - refDec)/epsAngle;
1829 nRefTime = (baryTime - refTime)/epsTime;
1830
1831 pForward = gsl_ran_ugaussian_pdf(nRA) * gsl_ran_ugaussian_pdf(nDec) * gsl_ran_ugaussian_pdf(nTime);
1832 pReverse = gsl_ran_ugaussian_pdf(nRefRA) * gsl_ran_ugaussian_pdf(nRefDec) * gsl_ran_ugaussian_pdf(nRefTime);
1833
1834 LALInferenceSetVariable(proposedParams, "rightascension", &newRA);
1835 LALInferenceSetVariable(proposedParams, "declination", &newDec);
1836
1837 if (timeflag)
1838 LALInferenceSetVariable(proposedParams, "time", &newTime);
1839
1840 logPropRatio = log(pReverse/pForward);
1841
1842 return logPropRatio;
1843}
1844
1847 LALInferenceVariables *proposedParams) {
1848 INT4 i,j;
1849 INT4 N, nifo;
1850 REAL8 draw=0.0;
1851 REAL8 logPropRatio = 0.0;
1852 REAL8Vector *var;
1853 gsl_matrix *ny;
1854
1856
1857 var = LALInferenceGetREAL8VectorVariable(thread->proposalArgs, "psdsigma");
1858
1859 //Get current state of chain into workable form
1860 ny = LALInferenceGetgslMatrixVariable(proposedParams, "psdscale");
1861
1862 //Get size of noise parameter array
1863 nifo = (INT4)ny->size1;
1864 N = (INT4)ny->size2;
1865
1866 //perturb noise parameter
1867 for(i=0; i<nifo; i++) {
1868 for(j=0; j<N; j++) {
1869 draw = gsl_matrix_get(ny, i, j) + gsl_ran_ugaussian(thread->GSLrandom) * var->data[j];
1870 gsl_matrix_set(ny, i, j, draw);
1871 }
1872 }
1873
1874 return logPropRatio;
1875}
1876
1878 LALInferenceVariables *proposedParams,
1879 gsl_matrix *glitchFD, INT4 ifo, INT4 n, INT4 flag) {
1880 INT4 i=0;
1881 INT4 lower, upper;
1882 INT4 glitchLower, glitchUpper;
1883 REAL8FrequencySeries **asds, *asd = NULL;
1884 REAL8Vector *flows;
1885 REAL8 deltaT, Tobs, deltaF;
1886 REAL8 Q, Amp, t0, ph0, f0; //sine-Gaussian parameters
1887 REAL8 amparg, phiarg, Ai;//helpers for computing sineGaussian
1888 REAL8 gRe, gIm; //real and imaginary parts of current glitch model
1889 REAL8 tau;
1890 gsl_matrix *glitch_f, *glitch_Q, *glitch_A;
1891 gsl_matrix *glitch_t, *glitch_p;
1892 REAL8TimeSeries **td_data;
1893
1895
1896 /* FIXME: can't store arrays of REAL8FrequencySeries in LI
1897 Variables (any more?)
1898
1899 Need instead Nifo variables called "asdH1", "asdL1", ...
1900 */
1902 flows = LALInferenceGetREAL8VectorVariable(args, "flows");
1903 td_data = *(REAL8TimeSeries ***)LALInferenceGetVariable(args, "td_data");
1904
1905 /* get dataPtr pointing to correct IFO */
1906 asd = asds[ifo];
1907
1908 deltaT = td_data[ifo]->deltaT;
1909 Tobs = (((REAL8)td_data[ifo]->data->length) * deltaT);
1910 deltaF = 1.0 / Tobs;
1911
1912 lower = (INT4)ceil(flows->data[ifo] / deltaF);
1913 upper = (INT4)floor(flows->data[ifo] / deltaF);
1914
1915 glitch_f = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_f0");
1916 glitch_Q = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_Q");
1917 glitch_A = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_Amp");
1918 glitch_t = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_t0");
1919 glitch_p = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_phi");
1920
1921 Q = gsl_matrix_get(glitch_Q, ifo, n);
1922 Amp = gsl_matrix_get(glitch_A, ifo, n);
1923 t0 = gsl_matrix_get(glitch_t, ifo, n);
1924 ph0 = gsl_matrix_get(glitch_p, ifo, n);
1925 f0 = gsl_matrix_get(glitch_f, ifo, n);
1926
1927 //6 x decay time of sine Gaussian (truncate how much glitch we compute)
1928 tau = Q/LAL_TWOPI/f0;
1929 glitchLower = (INT4)floor((f0 - 1./tau)/deltaF);
1930 glitchUpper = (INT4)floor((f0 + 1./tau)/deltaF);
1931
1932 //set glitch model to zero
1933 if (flag==0) {
1934 for(i=lower; i<=upper; i++) {
1935 gsl_matrix_set(glitchFD, ifo, 2*i, 0.0);
1936 gsl_matrix_set(glitchFD, ifo, 2*i+1, 0.0);
1937 }
1938 }
1939
1940 for (i=glitchLower; i<glitchUpper; i++) {
1941 if (i>=lower && i<=upper) {
1942 gRe = gsl_matrix_get(glitchFD, ifo, 2*i);
1943 gIm = gsl_matrix_get(glitchFD, ifo, 2*i+1);
1944 amparg = ((REAL8)i*deltaF - f0)*LAL_PI*tau;
1945 phiarg = LAL_PI*(REAL8)i + ph0 - LAL_TWOPI*(REAL8)i*deltaF*(t0-Tobs/2.);//TODO: SIMPLIFY PHASE FOR SINEGAUSSIAN
1946 Ai = Amp*tau*0.5*sqrt(LAL_PI)*exp(-amparg*amparg)*asd->data->data[i]/sqrt(Tobs);
1947
1948 switch(flag) {
1949 // Remove wavelet from model
1950 case -1:
1951 gRe -= Ai*cos(phiarg);
1952 gIm -= Ai*sin(phiarg);
1953 break;
1954 // Add wavelet to model
1955 case 1:
1956 gRe += Ai*cos(phiarg);
1957 gIm += Ai*sin(phiarg);
1958 break;
1959 // Replace model with wavelet
1960 case 0:
1961 gRe = Ai*cos(phiarg);
1962 gIm = Ai*sin(phiarg);
1963 break;
1964 //Do nothing
1965 default:
1966 break;
1967 }//end switch
1968
1969 //update glitch model
1970 gsl_matrix_set(glitchFD, ifo, 2*i, gRe);
1971 gsl_matrix_set(glitchFD, ifo, 2*i+1, gIm);
1972
1973 }//end upper/lower check
1974 }//end loop over glitch samples
1975}
1976
1978 REAL8 f0;
1979 REAL8 Q;
1980 REAL8 Amp;
1981
1982 REAL8 sigma_t0;
1983 REAL8 sigma_f0;
1984 REAL8 sigma_Q;
1985 REAL8 sigma_Amp;
1986 REAL8 sigma_phi0;
1987
1988 REAL8 SNR = 0.0;
1989 REAL8 sqrt3 = 1.7320508;
1990
1991 f0 = params->data[1];
1992 Q = params->data[2];
1993 Amp = params->data[3];
1994
1995 SNR = Amp*sqrt(Q/(2.0*sqrt(LAL_TWOPI)*f0));
1996
1997 // this caps the size of the proposed jumps
1998 if(SNR < 5.0) SNR = 5.0;
1999
2000 sigma_t0 = 1.0/(LAL_TWOPI*f0*SNR);
2001 sigma_f0 = 2.0*f0/(Q*SNR);
2002 sigma_Q = 2.0*Q/(sqrt3*SNR);
2003 sigma_Amp = Amp/SNR;
2004 sigma_phi0 = 1.0/SNR;
2005
2006 // Map diagonal Fisher elements to sigmas vector
2007 sigmas->data[0] = sigma_t0;
2008 sigmas->data[1] = sigma_f0;
2009 sigmas->data[2] = sigma_Q;
2010 sigmas->data[3] = sigma_Amp;
2011 sigmas->data[4] = sigma_phi0;
2012}
2013
2016 LALInferenceVariables *proposedParams) {
2017 INT4 i, ifo;
2018 INT4 n;
2019
2020 REAL8 logPropRatio = 0.0;
2021 REAL8 t0;
2022 REAL8 f0;
2023 REAL8 Q;
2024 REAL8 Amp;
2025 REAL8 phi0;
2026
2027 REAL8 scale;
2028
2029 REAL8 qyx;
2030 REAL8 qxy;
2031 REAL8 Anorm;
2032
2034
2035 gsl_matrix *glitchFD, *glitch_f, *glitch_Q, *glitch_A, *glitch_t, *glitch_p;
2036
2037 gsl_rng *rng = thread->GSLrandom;
2038 /*
2039 Vectors to store wavelet parameters.
2040 Order:
2041 [0] t0
2042 [1] f0
2043 [2] Q
2044 [3] Amp
2045 [4] phi0
2046 */
2047 REAL8Vector *params_x = XLALCreateREAL8Vector(5);
2048 REAL8Vector *params_y = XLALCreateREAL8Vector(5);
2049
2050 REAL8Vector *sigmas_x = XLALCreateREAL8Vector(5);
2051 REAL8Vector *sigmas_y = XLALCreateREAL8Vector(5);
2052
2053 /* Get glitch meta paramters (dimnsion, proposal) */
2054 UINT4Vector *gsize = LALInferenceGetUINT4VectorVariable(proposedParams, "glitch_size");
2055
2056 glitchFD = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_FD");
2057 glitch_f = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_f0");
2058 glitch_Q = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_Q");
2059 glitch_A = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_Amp");
2060 glitch_t = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_t0");
2061 glitch_p = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_phi");
2062
2063 Anorm = LALInferenceGetREAL8Variable(thread->priorArgs, "glitch_norm");
2064
2065 /* Choose which IFO */
2066 ifo = (INT4)floor(gsl_rng_uniform(rng) * (REAL8)(gsize->length));
2067
2068 /* Bail out of proposal if no wavelets */
2069 if (gsize->data[ifo]==0) {
2070 XLALDestroyREAL8Vector(params_x);
2071 XLALDestroyREAL8Vector(params_y);
2072 XLALDestroyREAL8Vector(sigmas_x);
2073 XLALDestroyREAL8Vector(sigmas_y);
2074
2075 return logPropRatio;
2076 }
2077
2078 /* Choose which glitch */
2079 n = (INT4)floor(gsl_rng_uniform(rng) * (REAL8)(gsize->data[ifo]));
2080
2081 /* Remove wavlet form linear combination */
2082 UpdateWaveletSum(thread, proposedParams, glitchFD, ifo, n, -1);
2083
2084 /* Get parameters of n'th glitch int params vector */
2085 t0 = gsl_matrix_get(glitch_t, ifo, n); //Centroid time
2086 f0 = gsl_matrix_get(glitch_f, ifo, n); //Frequency
2087 Q = gsl_matrix_get(glitch_Q, ifo, n); //Quality
2088 Amp = gsl_matrix_get(glitch_A, ifo, n); //Amplitude
2089 phi0 = gsl_matrix_get(glitch_p, ifo, n); //Centroid phase
2090
2091
2092 /* Map to params Vector and compute Fisher */
2093 params_x->data[0] = t0;
2094 params_x->data[1] = f0;
2095 params_x->data[2] = Q;
2096 params_x->data[3] = Amp * (0.25*Anorm);//TODO: What is the 0.25*Anorm about?
2097 params_x->data[4] = phi0;
2098
2099 MorletDiagonalFisherMatrix(params_x, sigmas_x);
2100
2101 /* Jump from x -> y: y = x + N[0,sigmas_x]*scale */
2102 scale = 0.4082482; // 1/sqrt(6)
2103
2104 for(i=0; i<5; i++)
2105 params_y->data[i] = params_x->data[i] + gsl_ran_ugaussian(rng)*sigmas_x->data[i]*scale;
2106
2107 /* Set parameters of n'th glitch int params vector */
2108 /* Map to params Vector and compute Fisher */
2109 t0 = params_y->data[0];
2110 f0 = params_y->data[1];
2111 Q = params_y->data[2];
2112 Amp = params_y->data[3]/(0.25*Anorm);
2113 phi0 = params_y->data[4];
2114
2115 gsl_matrix_set(glitch_t, ifo, n, t0);
2116 gsl_matrix_set(glitch_f, ifo, n, f0);
2117 gsl_matrix_set(glitch_Q, ifo, n, Q);
2118 gsl_matrix_set(glitch_A, ifo, n, Amp);
2119 gsl_matrix_set(glitch_p, ifo, n, phi0);
2120
2121 /* Add wavlet to linear combination */
2122 UpdateWaveletSum(thread, proposedParams, glitchFD, ifo, n, 1);
2123
2124 /* Now compute proposal ratio using Fisher at y */
2125 MorletDiagonalFisherMatrix(params_y, sigmas_y);
2126
2127 REAL8 sx = 1.0; // sigma
2128 REAL8 sy = 1.0;
2129 REAL8 dx = 1.0; // (params_x - params_y)/sigma
2130 REAL8 dy = 1.0;
2131 REAL8 exy = 0.0; // argument of exponential part of q
2132 REAL8 eyx = 0.0;
2133 REAL8 nxy = 1.0; // argument of normalization part of q
2134 REAL8 nyx = 1.0; // (computed as product to avoid too many log()s )
2135 for(i=0; i<5; i++) {
2136 sx = scale*sigmas_x->data[i];
2137 sy = scale*sigmas_y->data[i];
2138
2139 dx = (params_x->data[i] - params_y->data[i])/sx;
2140 dy = (params_x->data[i] - params_y->data[i])/sy;
2141
2142 nxy *= sy;
2143 nyx *= sx;
2144
2145 exy += -dy*dy/2.0;
2146 eyx += -dx*dx/2.0;
2147 }
2148
2149 qyx = eyx - log(nyx); //probabiltiy of proposing y given x
2150 qxy = exy - log(nxy); //probability of proposing x given y
2151
2152 logPropRatio = qxy-qyx;
2153
2154 XLALDestroyREAL8Vector(params_x);
2155 XLALDestroyREAL8Vector(params_y);
2156
2157 XLALDestroyREAL8Vector(sigmas_x);
2158 XLALDestroyREAL8Vector(sigmas_y);
2159
2160 return logPropRatio;
2161}
2162
2165 LALInferenceVariables *proposedParams) {
2166 INT4 i,n;
2167 INT4 ifo;
2168 INT4 rj, nx, ny;
2169 REAL8 draw;
2170 REAL8 val, Anorm;
2171 REAL8 t=0, f=0, Q=0, A=0;
2172 REAL8 qx = 0.0; //log amp proposals
2173 REAL8 qy = 0.0;
2174 REAL8 qyx = 0.0; //log pixel proposals
2175 REAL8 qxy = 0.0;
2176 REAL8 pForward = 0.0; //combined p() & q() probabilities for ...
2177 REAL8 pReverse = 0.0; //...RJMCMC hastings ratio
2178 REAL8 logPropRatio = 0.0;
2179 INT4 adapting=1;
2180 gsl_matrix *params = NULL;
2181
2182 gsl_rng *rng = thread->GSLrandom;
2183 LALInferenceVariables *propArgs = thread->proposalArgs;
2184
2186
2187 UINT4Vector *gsize = LALInferenceGetUINT4VectorVariable(proposedParams, "glitch_size");
2188 gsl_matrix *glitchFD = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_FD");
2189
2190 INT4 nmin = (INT4)(LALInferenceGetREAL8Variable(thread->priorArgs,"glitch_dim_min"));
2191 INT4 nmax = (INT4)(LALInferenceGetREAL8Variable(thread->priorArgs,"glitch_dim_max"));
2192
2193 if (LALInferenceCheckVariable(propArgs, "adapting"))
2194 adapting = LALInferenceGetINT4Variable(propArgs, "adapting");
2195
2196 /* Choose which IFO */
2197 ifo = (INT4)floor(gsl_rng_uniform(rng) * (REAL8)(gsize->length) );
2198 nx = gsize->data[ifo];
2199
2200 /* Choose birth or death move */
2201 draw = gsl_rng_uniform(rng);
2202 if (draw < 0.5)
2203 rj = 1;
2204 else
2205 rj = -1;
2206
2207 /* find dimension of proposed model */
2208 ny = nx + rj;
2209
2210 /* Check that new dimension is allowed */
2211 if(ny<nmin || ny>=nmax) {
2212 logPropRatio = -INFINITY;
2213 return logPropRatio;
2214 }
2215
2216 switch(rj) {
2217 /* Birth */
2218 case 1:
2219 //Add new wavelet to glitch model
2220 t = draw_flat(thread, "morlet_t0_prior");
2221 f = draw_flat(thread, "morlet_f0_prior");
2222
2223 //Centroid time
2224 params = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_t0");
2225 gsl_matrix_set(params, ifo, nx, t);
2226
2227 //Frequency
2228 params = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_f0");
2229 gsl_matrix_set(params, ifo, nx, f);
2230
2231 //Quality
2232 params = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_Q");
2233 val = draw_flat(thread, "morlet_Q_prior");
2234 gsl_matrix_set(params, ifo, nx, val);
2235 Q = val;
2236
2237 //Amplitude
2238 params = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_Amp");
2239 val = glitchAmplitudeDraw(Q, f, rng);
2240 Anorm = LALInferenceGetREAL8Variable(thread->priorArgs, "glitch_norm");
2241 A = val/Anorm;
2242
2243 gsl_matrix_set(params, ifo, nx, A);
2244
2245 //Centroid phase
2246 params = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_phi");
2247 val = draw_flat(thread, "morlet_phi_prior");
2248 gsl_matrix_set(params, ifo, nx, val);
2249
2250 //Add wavlet to linear combination
2251 UpdateWaveletSum(thread, proposedParams, glitchFD, ifo, nx, 1);
2252
2253 //Compute probability of drawing parameters
2254 qy = evaluate_morlet_proposal(thread, proposedParams, ifo, nx);// + log(gsl_matrix_get(power,ifo,k));
2255
2256 //Compute reverse probability of dismissing k
2257 qxy = 0.0;//-log((double)runState->data->freqData->data->length);//-log( (REAL8)ny );
2258
2259 if (adapting)
2260 qy += 10.0;
2261
2262 break;
2263
2264 /* Death */
2265 case -1:
2266 //Choose wavelet to remove from glitch model
2267 draw = gsl_rng_uniform(rng);
2268 n = (INT4)(floor(draw * (REAL8)nx)); //choose which hot pixel
2269
2270 // Remove wavlet from linear combination
2271 UpdateWaveletSum(thread, proposedParams, glitchFD, ifo, n, -1);
2272
2273 //Get t and f of removed wavelet
2274 params = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_f0");
2275 f = gsl_matrix_get(params, ifo, n);
2276 params = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_t0");
2277 t = gsl_matrix_get(params, ifo, n);
2278 params = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_Q");
2279 Q = gsl_matrix_get(params, ifo, n);
2280 params = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_Amp");
2281 A = gsl_matrix_get(params, ifo, n);
2282
2283 //Shift morlet parameters to fill in array
2284 for(i=n; i<ny; i++) {
2285 params = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_f0");
2286 gsl_matrix_set(params, ifo, i, gsl_matrix_get(params, ifo, i+1));
2287 params = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_Q");
2288 gsl_matrix_set(params, ifo, i, gsl_matrix_get(params, ifo, i+1));
2289 params = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_Amp");
2290 gsl_matrix_set(params, ifo, i, gsl_matrix_get(params, ifo, i+1));
2291 params = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_t0");
2292 gsl_matrix_set(params, ifo, i, gsl_matrix_get(params, ifo, i+1));
2293 params = LALInferenceGetgslMatrixVariable(proposedParams, "morlet_phi");
2294 gsl_matrix_set(params, ifo, i, gsl_matrix_get(params, ifo, i+1));
2295 }
2296
2297 //Compute reverse probability of drawing parameters
2298 //find TF pixel
2299
2300 qx = evaluate_morlet_proposal(thread, currentParams, ifo, n);// + log(gsl_matrix_get(power,ifo,k));
2301
2302 //Compute forward probability of dismissing k
2303 qyx = 0.0;//-log((double)runState->data->freqData->data->length);//0.0;//-log( (REAL8)nx );
2304
2305 if(adapting)
2306 qx += 10.0;
2307
2308 break;
2309
2310 default:
2311 break;
2312 }
2313
2314 /* Update proposal structure for return to MCMC */
2315
2316 //Update model meta-date
2317 gsize->data[ifo] = ny;
2318
2319 //Re-package prior and proposal ratios into runState
2320 pForward = qxy + qx;
2321 pReverse = qyx + qy;
2322
2323 logPropRatio = pForward-pReverse;
2324
2325 return logPropRatio;
2326}
2327
2330 LALInferenceVariables *proposedParams) {
2331 REAL8 logPropRatio = 0.0;
2332
2334
2335 REAL8 psi = LALInferenceGetREAL8Variable(proposedParams, "polarisation");
2336 REAL8 phi = LALInferenceGetREAL8Variable(proposedParams, "phase");
2337
2338 phi += M_PI;
2339 psi += M_PI/2;
2340
2341 phi = fmod(phi, 2.0*M_PI);
2342 psi = fmod(psi, M_PI);
2343
2344 LALInferenceSetVariable(proposedParams, "polarisation", &psi);
2345 LALInferenceSetVariable(proposedParams, "phase", &phi);
2346
2347 return logPropRatio;
2348}
2349
2352 LALInferenceVariables *proposedParams) {
2353 REAL8 alpha,beta;
2354 REAL8 draw;
2355 REAL8 psi, phi;
2356 REAL8 logPropRatio = 0.0;
2357
2359
2360 gsl_rng *rng = thread->GSLrandom;
2361
2362 psi = LALInferenceGetREAL8Variable(proposedParams, "polarisation");
2363 phi = LALInferenceGetREAL8Variable(proposedParams, "phase");
2364
2365 alpha = psi + phi;
2366 beta = psi - phi;
2367
2368 //alpha => 0:3pi
2369 //beta => -2pi:pi
2370
2371 //big jump in either alpha (beta) or beta (alpha)
2372 draw = gsl_rng_uniform(rng);
2373 if (draw < 0.5)
2374 alpha = gsl_rng_uniform(rng)*3.0*LAL_PI;
2375 else
2376 beta = -LAL_TWOPI + gsl_rng_uniform(rng)*3.0*LAL_PI;
2377
2378 //transform back to psi,phi space
2379 psi = (alpha + beta)*0.5;
2380 phi = (alpha - beta)*0.5;
2381
2382 //map back in range
2383 LALInferenceCyclicReflectiveBound(proposedParams, thread->priorArgs);
2384
2385 LALInferenceSetVariable(proposedParams, "polarisation", &psi);
2386 LALInferenceSetVariable(proposedParams, "phase", &phi);
2387
2388 return logPropRatio;
2389}
2390
2393 LALInferenceVariables *proposedParams) {
2394 REAL8 f0, df;
2395 REAL8 plusminus;
2396 REAL8 logPropRatio = 0.0;
2397
2399
2400 f0 = LALInferenceGetREAL8Variable(proposedParams, "f0");
2401 df = LALInferenceGetREAL8Variable(proposedParams, "df");
2402
2403 plusminus = gsl_rng_uniform(thread->GSLrandom);
2404 if ( plusminus < 0.5 )
2405 f0 -= df;
2406 else
2407 f0 += df;
2408
2409 LALInferenceSetVariable(proposedParams, "f0", &f0);
2410
2411 return logPropRatio;
2412}
2413
2414//This proposal needs to be called with exactly 3 independent detector locations.
2416 const REAL8 baryTime, const REAL8 dist, const REAL8 iota, const REAL8 psi,
2417 REAL8 *newRA, REAL8 *newDec, REAL8 *newTime,
2418 REAL8 *newDist, REAL8 *newIota, REAL8 *newPsi) {
2419 REAL8 R2[4];
2420 REAL8 dist2;
2421 REAL8 gmst, newGmst;
2422 REAL8 cosIota, cosIota2;
2423 REAL8 Fplus, Fcross, psi_temp;
2424 REAL8 x[4], y[4], x2[4], y2[4];
2425 REAL8 newFplus[4], newFplus2[4], newFcross[4], newFcross2[4];
2426 REAL8 a, a2, b, c12;
2427 REAL8 cosnewIota, cosnewIota2;
2428 LIGOTimeGPS GPSlal;
2429 INT4 nUniqueDet, det;
2431
2432 get_detectors(thread->parent->data, &detectors);
2433 nUniqueDet = LALInferenceGetINT4Variable(thread->proposalArgs, "nUniqueDet");
2434
2435 XLALGPSSetREAL8(&GPSlal, baryTime);
2436 gmst = XLALGreenwichMeanSiderealTime(&GPSlal);
2437
2438 reflected_position_and_time(thread, ra, dec, baryTime, newRA, newDec, newTime);
2439
2440 XLALGPSSetREAL8(&GPSlal, *newTime);
2441 newGmst = XLALGreenwichMeanSiderealTime(&GPSlal);
2442
2443 dist2 = dist*dist;
2444
2445 cosIota = cos(iota);
2446 cosIota2 = cosIota*cosIota;
2447
2448 /* Loop over interferometers */
2449 INT4 i=1, j=0;
2450 for (det=0; det < nUniqueDet; det++) {
2451 psi_temp = 0.0;
2452
2453 XLALComputeDetAMResponse(&Fplus, &Fcross, (const REAL4(*)[3])detectors[det].response, *newRA, *newDec, psi_temp, newGmst);
2454 j=i-1;
2455 while (j>0){
2456 if (Fplus==x[j]){
2457 det++;
2458 XLALComputeDetAMResponse(&Fplus, &Fcross, (const REAL4(*)[3])detectors[det].response, *newRA, *newDec, psi_temp, newGmst);
2459 }
2460 j--;
2461 }
2462 x[i]=Fplus;
2463 x2[i]=Fplus*Fplus;
2464 y[i]=Fcross;
2465 y2[i]=Fcross*Fcross;
2466
2467 XLALComputeDetAMResponse(&Fplus, &Fcross, (const REAL4(*)[3])detectors[det].response, ra, dec, psi, gmst);
2468 R2[i] = (((1.0+cosIota2)*(1.0+cosIota2))/(4.0*dist2))*Fplus*Fplus
2469 + ((cosIota2)/(dist2))*Fcross*Fcross;
2470
2471 i++;
2472 }
2473
2474 a = (R2[3]*x2[2]*y2[1] - R2[2]*x2[3]*y2[1] - R2[3]*x2[1]*y2[2] + R2[1]*x2[3]*y2[2] + R2[2]*x2[1]*y2[3] -
2475 R2[1]*x2[2]*y2[3]);
2476 a2 = a*a;
2477 b = (-(R2[3]*x[1]*x2[2]*y[1]) + R2[2]*x[1]*x2[3]*y[1] + R2[3]*x2[1]*x[2]*y[2] - R2[1]*x[2]*x2[3]*y[2] +
2478 R2[3]*x[2]*y2[1]*y[2] - R2[3]*x[1]*y[1]*y2[2] - R2[2]*x2[1]*x[3]*y[3] + R2[1]*x2[2]*x[3]*y[3] - R2[2]*x[3]*y2[1]*y[3] + R2[1]*x[3]*y2[2]*y[3] +
2479 R2[2]*x[1]*y[1]*y2[3] - R2[1]*x[2]*y[2]*y2[3]);
2480
2481 (*newPsi) = (2.*atan((b - a*sqrt((a2 + b*b)/(a2)))/a))/4.;
2482
2483 while ((*newPsi)<0){
2484 (*newPsi)=(*newPsi)+LAL_PI/4.0;
2485 }
2486
2487 while ((*newPsi)>LAL_PI/4.0){
2488 (*newPsi)=(*newPsi)-LAL_PI/4.0;
2489 }
2490
2491 for (i = 1; i < 4; i++){
2492 newFplus[i] = x[i]*cos(2.0*(*newPsi)) + y[i]*sin(2.0*(*newPsi));
2493 newFplus2[i] = newFplus[i] * newFplus[i];
2494
2495 newFcross[i] = y[i]*cos(2.0*(*newPsi)) - x[i]*sin(2.0*(*newPsi));
2496 newFcross2[i] = newFcross[i] * newFcross[i];
2497 }
2498
2499 c12 = -2.0*((R2[1]*(newFcross2[2])-R2[2]*(newFcross2[1]))
2500 /(R2[1]*(newFplus2[2])-R2[2]*(newFplus2[1])))-1.0;
2501
2502 if (c12<1.0){
2503 c12 = (3.0-c12)/(1.0+c12);
2504 (*newPsi) = (*newPsi)+LAL_PI/4.0;
2505
2506 for (i = 1; i < 4; i++){
2507 newFplus[i] = x[i]*cos(2.0*(*newPsi)) + y[i]*sin(2.0*(*newPsi));
2508 newFplus2[i] = newFplus[i] * newFplus[i];
2509
2510 newFcross[i] = y[i]*cos(2.0*(*newPsi)) - x[i]*sin(2.0*(*newPsi));
2511 newFcross2[i] = newFcross[i] * newFcross[i];
2512 }
2513 }
2514
2515 if (c12<1){
2516 *newIota = iota;
2517 *newDist = dist;
2518 return;
2519 }
2520
2521 cosnewIota2 = c12-sqrt(c12*c12-1.0);
2522 cosnewIota = sqrt(cosnewIota2);
2523 *newIota = acos(cosnewIota);
2524
2525 *newDist = sqrt((
2526 ((((1.0+cosnewIota2)*(1.0+cosnewIota2))/(4.0))*newFplus2[1]
2527 + (cosnewIota2)*newFcross2[1])
2528 )/ R2[1]);
2529
2530 if (Fplus*newFplus[3]<0){
2531 (*newPsi)=(*newPsi)+LAL_PI/2.;
2532 newFcross[3]=-newFcross[3];
2533 }
2534
2535 if (Fcross*cosIota*cosnewIota*newFcross[3]<0){
2536 (*newIota)=LAL_PI-(*newIota);
2537 }
2538
2540}
2541
2544 LALInferenceVariables *proposedParams) {
2545 REAL8 old_d=0.0;
2546 DistanceParam distParam;
2547
2548 if (LALInferenceCheckVariable(currentParams, "distance")) {
2549 distParam = USES_DISTANCE_VARIABLE;
2550 old_d = LALInferenceGetREAL8Variable(currentParams,"distance");
2551 } else if (LALInferenceCheckVariable(currentParams, "logdistance")) {
2552 distParam = USES_LOG_DISTANCE_VARIABLE;
2553 old_d = exp(LALInferenceGetREAL8Variable(currentParams,"logdistance"));
2554 } else {
2555 XLAL_ERROR_REAL8(XLAL_FAILURE, "could not find 'distance' or 'logdistance' in current params");
2556 }
2557 if(!LALInferenceCheckVariable(currentParams,"optimal_snr") || !LALInferenceCheckVariable(currentParams,"matched_filter_snr"))
2558 /* Force a likelihood calculation to generate the parameters */
2559 thread->parent->likelihood(currentParams,thread->parent->data,thread->model);
2560 if(!LALInferenceCheckVariable(currentParams,"optimal_snr") || !LALInferenceCheckVariable(currentParams,"matched_filter_snr"))
2561 {
2562 /* If still can't find the required parameters, reject this jump */
2564 return(-INFINITY);
2565 }
2566 REAL8 OptimalSNR = LALInferenceGetREAL8Variable(currentParams,"optimal_snr");
2567 REAL8 MatchedFilterSNR = LALInferenceGetREAL8Variable(currentParams,"matched_filter_snr");
2568 REAL8 d_inner_h = MatchedFilterSNR * OptimalSNR;
2569
2570 /* Get params of sampling distribution */
2571 REAL8 xmean = d_inner_h/(OptimalSNR*OptimalSNR*old_d);
2572 REAL8 xsigma = 1.0/(OptimalSNR*old_d);
2573 REAL8 old_x = 1.0/old_d;
2574
2575 /* Sample new x. Do not check for x<0 since that can throw the code
2576 * into a difficult-to-terminate loop */
2577 REAL8 new_x;
2578 new_x = xmean + gsl_ran_gaussian(thread->GSLrandom,xsigma);
2579 REAL8 new_d = 1.0/new_x;
2580
2582 /* Adjust SNRs */
2583 OptimalSNR *= new_x / old_x;
2584 LALInferenceSetREAL8Variable(proposedParams,"optimal_snr",OptimalSNR);
2585
2586 /* Update individual detector information */
2587 const char *ifonames[5]={"H1","L1","V1","I1","J1"};
2588 char pname[64]="";
2589 for(UINT4 i=0;i<5;i++)
2590 {
2591 sprintf(pname,"%s_optimal_snr",ifonames[i]);
2593 LALInferenceSetREAL8Variable(proposedParams,pname,LALInferenceGetREAL8Variable(currentParams,pname) * (new_x/old_x));
2594 }
2595
2596 REAL8 logxdjac;
2597 /* Set distance */
2598 if(distParam == USES_DISTANCE_VARIABLE)
2599 {
2600 /* The natural proposal density is in x, but we need to compute
2601 the proposal ratio density in d. So, we need a factor of
2602
2603 |dx_old/dd_old| / |dx_new/dd_new| = d_new^2 / d_old^2
2604
2605 to correct the x proposal density.
2606 */
2607 LALInferenceSetVariable(proposedParams,"distance",&new_d);
2608 logxdjac = 2.0*log(new_d) - 2.0*log(old_d);
2609 }
2610 else
2611 {
2612 /* The natural proposal density is in x, but we need to compute
2613 the proposal ratio in log(d). So, we need a factor of
2614
2615 |dx_old/dlog(d_old)| / |dx_new/dlog(d_new)| = d_new / d_old
2616
2617 to correct the x proposal density.
2618 */
2619 REAL8 new_logd = log(new_d);
2620 LALInferenceSetVariable(proposedParams,"logdistance",&new_logd);
2621 logxdjac = log(new_d) - log(old_d);
2622 }
2623 /* Proposal ratio is not symmetric */
2624 /* P.R. = p(xold|xmean,xsigma) / p(xnew|xmean,xsigma) * jacobian */
2625 REAL8 log_p_reverse = -0.5*(old_x-xmean)*(old_x-xmean)/(xsigma*xsigma);
2626 REAL8 log_p_forward = -0.5*(new_x-xmean)*(new_x-xmean)/(xsigma*xsigma);
2627
2628 return(log_p_reverse - log_p_forward + logxdjac);
2629}
2630
2631
2634 LALInferenceVariables *proposedParams) {
2635 INT4 nUniqueDet;
2636 INT4 timeflag=0;
2637 REAL8 baryTime;
2638 REAL8 ra, dec;
2639 REAL8 psi, dist;
2640 REAL8 newRA, newDec, newTime, newDist, newIota, newPsi;
2641 REAL8 nRA, nDec, nTime, nDist, nIota, nPsi;
2642 REAL8 refRA, refDec, refTime, refDist, refIota, refPsi;
2643 REAL8 nRefRA, nRefDec, nRefTime, nRefDist, nRefIota, nRefPsi;
2644 REAL8 pForward, pReverse;
2645 REAL8 cst;
2646 REAL8 iota=0.0;
2647 REAL8 logPropRatio = 0.0;
2648 /* Find the number of distinct-position detectors. */
2649 /* Exit with same parameters (with a warning the first time) if
2650 there are not EXACTLY three unique detector locations. */
2651 static INT4 warningDelivered = 0;
2653
2654
2656 gsl_rng *rng = thread->GSLrandom;
2657 epoch = thread->parent->data->epoch;
2658
2659 nUniqueDet = LALInferenceGetINT4Variable(args, "nUniqueDet");
2660 if (nUniqueDet != 3) {
2661 if (!warningDelivered) {
2662 fprintf(stderr, "WARNING: trying to reflect through the decector plane with %d\n", nUniqueDet);
2663 fprintf(stderr, "WARNING: geometrically independent locations,\n");
2664 fprintf(stderr, "WARNING: but this proposal should only be used with exactly 3 independent detectors.\n");
2665 fprintf(stderr, "WARNING: %s, line %d\n", __FILE__, __LINE__);
2666 warningDelivered = 1;
2667 }
2668
2669 return logPropRatio;
2670 }
2672
2673
2674 ra = LALInferenceGetREAL8Variable(proposedParams, "rightascension");
2675 dec = LALInferenceGetREAL8Variable(proposedParams, "declination");
2676
2677 if (LALInferenceCheckVariable(proposedParams,"time")){
2678 baryTime = LALInferenceGetREAL8Variable(proposedParams, "time");
2679 timeflag = 1;
2680 } else {
2681 baryTime = XLALGPSGetREAL8(&epoch);
2682 }
2683
2684 if (LALInferenceCheckVariable(proposedParams,"costheta_jn"))
2685 iota = acos(LALInferenceGetREAL8Variable(proposedParams, "costheta_jn"));
2686 else
2687 fprintf(stderr, "LALInferenceExtrinsicParamProposal: No theta_jn parameter!\n");
2688
2689 psi = LALInferenceGetREAL8Variable(proposedParams, "polarisation");
2690
2691 dist = exp(LALInferenceGetREAL8Variable(proposedParams, "logdistance"));
2692
2693 reflected_extrinsic_parameters(thread, ra, dec, baryTime, dist, iota, psi, &newRA, &newDec, &newTime, &newDist, &newIota, &newPsi);
2694
2695 /* Unit normal deviates, used to "fuzz" the state. */
2696 const REAL8 epsDist = 1e-8;
2697 const REAL8 epsTime = 1e-8;
2698 const REAL8 epsAngle = 1e-8;
2699
2700 nRA = gsl_ran_ugaussian(rng);
2701 nDec = gsl_ran_ugaussian(rng);
2702 nTime = gsl_ran_ugaussian(rng);
2703 nDist = gsl_ran_ugaussian(rng);
2704 nIota = gsl_ran_ugaussian(rng);
2705 nPsi = gsl_ran_ugaussian(rng);
2706
2707 newRA += epsAngle*nRA;
2708 newDec += epsAngle*nDec;
2709 newTime += epsTime*nTime;
2710 newDist += epsDist*nDist;
2711 newIota += epsAngle*nIota;
2712 newPsi += epsAngle*nPsi;
2713
2714 /* And the doubly-reflected position (near the original, but not
2715 exactly due to the fuzzing). */
2716 reflected_extrinsic_parameters(thread, newRA, newDec, newTime, newDist, newIota, newPsi, &refRA, &refDec, &refTime, &refDist, &refIota, &refPsi);
2717
2718 /* The Gaussian increments required to shift us back to the original
2719 position from the doubly-reflected position. */
2720 nRefRA = (ra - refRA)/epsAngle;
2721 nRefDec = (dec - refDec)/epsAngle;
2722 nRefTime = (baryTime - refTime)/epsTime;
2723 nRefDist = (dist - refDist)/epsDist;
2724 nRefIota = (iota - refIota)/epsAngle;
2725 nRefPsi = (psi - refPsi)/epsAngle;
2726
2727 cst = log(1./(sqrt(2.*LAL_PI)));
2728 pReverse = 6*cst-0.5*(nRefRA*nRefRA+nRefDec*nRefDec+nRefTime*nRefTime+nRefDist*nRefDist+nRefIota*nRefIota+nRefPsi*nRefPsi);
2729 pForward = 6*cst-0.5*(nRA*nRA+nDec*nDec+nTime*nTime+nDist*nDist+nIota*nIota+nPsi*nPsi);
2730
2731 LALInferenceSetVariable(proposedParams, "rightascension", &newRA);
2732 LALInferenceSetVariable(proposedParams, "declination", &newDec);
2733 if (timeflag)
2734 LALInferenceSetVariable(proposedParams, "time", &newTime);
2735
2736 REAL8 logNewDist = log(newDist);
2737 LALInferenceSetVariable(proposedParams, "logdistance", &logNewDist);
2738
2739 REAL8 newcosIota = cos(newIota);
2740 LALInferenceSetVariable(proposedParams, "costheta_jn", &newcosIota);
2741 LALInferenceSetVariable(proposedParams, "polarisation", &newPsi);
2742
2743 logPropRatio = pReverse - pForward;
2744
2745 return logPropRatio;
2746}
2747
2748
2750 INT4 i, nDet;
2751 REAL8Vector *flows, *fhighs;
2752 REAL8FrequencySeries **asds, **psds;
2753 REAL8TimeSeries **td_data;
2754 COMPLEX16FrequencySeries **fd_data;
2755 REAL8FFTPlan **plans;
2756
2757 nDet = LALInferenceGetINT4Variable(propArgs, "nDet");
2758
2759 flows = XLALCreateREAL8Vector(nDet);
2760 fhighs = XLALCreateREAL8Vector(nDet);
2761 asds = XLALCalloc(nDet, sizeof(REAL8FrequencySeries *));
2762 psds = XLALCalloc(nDet, sizeof(REAL8FrequencySeries *));
2763 td_data = XLALCalloc(nDet, sizeof(REAL8TimeSeries *));
2764 fd_data = XLALCalloc(nDet, sizeof(COMPLEX16FrequencySeries *));
2765 plans = XLALCalloc(nDet, sizeof(REAL8FFTPlan *));
2766
2767 for (i=0; i<nDet; i++) {
2768 flows->data[i] = data->fLow;
2769 fhighs->data[i] = data->fHigh;
2770
2771 asds[i] = data->noiseASD;
2772 psds[i] = data->oneSidedNoisePowerSpectrum;
2773
2774 td_data[i] = data->timeData;
2775 fd_data[i] = data->freqData;
2776
2777 plans[i] = data->freqToTimeFFTPlan;
2778 data = data->next;
2779 }
2780
2788
2789}
2790
2791
2792/** Setup adaptive proposals. Should be called when state->currentParams is already filled with an initial sample */
2794 INT4 no_adapt, adapting;
2795 INT4 adaptTau, adaptableStep, adaptLength, adaptResetBuffer;
2796 REAL8 sigma, s_gamma;
2797 REAL8 logPAtAdaptStart = -INFINITY;
2798
2800
2801 for(this=params->head; this; this=this->next) {
2803 char *name = this->name;
2804
2805 if (!strcmp(name, "eta") || !strcmp(name, "q") || !strcmp(name, "time") || !strcmp(name, "a_spin2") || !strcmp(name, "a_spin1") || !strcmp(name,"t0")){
2806 sigma = 0.001;
2807 } else if (!strcmp(name, "polarisation") || !strcmp(name, "phase") || !strcmp(name, "costheta_jn")){
2808 sigma = 0.1;
2809 } else {
2810 sigma = 0.01;
2811 }
2812
2813 /* Set up variables to store current sigma, proposed and accepted */
2814 char varname[MAX_STRLEN] = "";
2815 sprintf(varname, "%s_%s", name, ADAPTSUFFIX);
2817
2818 sigma = 0.0;
2819 sprintf(varname, "%s_%s", name, ACCEPTSUFFIX);
2821
2822 sprintf(varname, "%s_%s", name, PROPOSEDSUFFIX);
2824 }
2825 }
2826
2827 no_adapt = LALInferenceGetINT4Variable(propArgs, "no_adapt");
2828 adapting = !no_adapt; // Indicates if current iteration is being adapted
2829 LALInferenceAddINT4Variable(propArgs, "adapting", adapting, LALINFERENCE_PARAM_LINEAR);
2830
2831 char name[MAX_STRLEN] = "none";
2832 LALInferenceAddstringVariable(propArgs, "proposedVariableName", name, LALINFERENCE_PARAM_OUTPUT);
2833
2834 adaptableStep = 0;
2835 adaptTau = LALInferenceGetINT4Variable(propArgs, "adaptTau"); // Sets decay of adaption function
2836 adaptLength = pow(10, adaptTau); // Number of iterations to adapt before turning off
2837 adaptResetBuffer = 100; // Number of iterations before adapting after a restart
2838 s_gamma = 1.0; // Sets the size of changes to jump size during adaptation
2839
2840 LALInferenceAddINT4Variable(propArgs, "adaptableStep", adaptableStep, LALINFERENCE_PARAM_LINEAR);
2841 LALInferenceAddINT4Variable(propArgs, "adaptLength", adaptLength, LALINFERENCE_PARAM_LINEAR);
2842 LALInferenceAddINT4Variable(propArgs, "adaptResetBuffer", adaptResetBuffer, LALINFERENCE_PARAM_LINEAR);
2843 LALInferenceAddREAL8Variable(propArgs, "s_gamma", s_gamma, LALINFERENCE_PARAM_LINEAR);
2844 LALInferenceAddREAL8Variable(propArgs, "logPAtAdaptStart", logPAtAdaptStart, LALINFERENCE_PARAM_LINEAR);
2845
2846 return;
2847}
2848
2849
2850/** Update proposal statistics if tracking */
2852 INT4 i = 0;
2853
2854 LALInferenceProposal *prop = thread->cycle->proposals[i];
2855
2856 /* Find the proposal that was last called (by name) */
2857 while (strcmp(prop->name, thread->cycle->last_proposal_name)) {
2858 i++;
2859 prop = thread->cycle->proposals[i];
2860 }
2861
2862 /* Update proposal statistics */
2863 prop->proposed++;
2864 if (thread->accepted == 1){
2865 prop->accepted++;
2866 }
2867
2868 return;
2869}
2870
2871/* Zero out proposal statistics counters */
2873 INT4 i=0;
2874
2875 for (i=0; i<cycle->nProposals; i++) {
2876 LALInferenceProposal *prop = cycle->proposals[i];
2877
2878 prop->proposed = 0;
2879 prop->accepted = 0;
2880 }
2881
2882 return;
2883}
2884
2885/** Update the adaptive proposal. Whether or not a jump was accepted is passed with accepted */
2887 INT4 adaptableStep = 0;
2888 INT4 adapting = 0;
2889 REAL8 priorMin, priorMax, dprior, s_gamma;
2890 REAL8 accept, propose, sigma;
2891 const char *name;
2892
2894
2895 if (LALInferenceCheckVariable(args, "adaptableStep" ) &&
2896 LALInferenceCheckVariable(args, "adapting" )){
2897 adaptableStep = LALInferenceGetINT4Variable(args, "adaptableStep");
2898 adapting = LALInferenceGetINT4Variable(args, "adapting");
2899 }
2900 /* Don't do anything if these are not found */
2901 else return;
2902
2903 if (adaptableStep && adapting) {
2904 char tmpname[MAX_STRLEN] = "";
2905 name = LALInferenceGetstringVariable(args, "proposedVariableName");
2906
2907 sprintf(tmpname, "%s_%s", name, PROPOSEDSUFFIX);
2908 propose = LALInferenceGetREAL8Variable(args, tmpname) + 1;
2909 LALInferenceSetVariable(args, tmpname, &propose);
2910
2911 sprintf(tmpname, "%s_%s", name, ACCEPTSUFFIX);
2912 accept = LALInferenceGetREAL8Variable(args, tmpname) + thread->accepted;
2913 LALInferenceSetVariable(args, tmpname, &accept);
2914 }
2915
2916 /* Adapt if desired. */
2917 if (LALInferenceCheckVariable(args, "proposedVariableName") &&
2918 LALInferenceCheckVariable(args, "s_gamma") &&
2919 LALInferenceCheckVariable(args, "adapting") &&
2920 LALInferenceCheckVariable(args, "adaptableStep")) {
2921
2922 if (adaptableStep) {
2923 name = LALInferenceGetstringVariable(args, "proposedVariableName");
2924 char tmpname[MAX_STRLEN]="";
2925
2926 sprintf(tmpname, "%s_%s", name, ADAPTSUFFIX);
2927
2928 s_gamma = LALInferenceGetREAL8Variable(args, "s_gamma");
2929
2930 sigma = LALInferenceGetREAL8Variable(args, tmpname);
2931
2933 {
2934 LALInferenceGetMinMaxPrior(thread->priorArgs, name, &priorMin, &priorMax);
2935 dprior = priorMax - priorMin;
2936 }
2938 {
2939 REAL8 mu;
2940 /* Get std. dev. as dprior */
2942 }
2943 else
2944 {
2945 return;
2946 }
2947
2948 if (thread->accepted == 1){
2949 sigma += s_gamma * (dprior/100.0) * (1.0-targetAcceptance);
2950 } else {
2951 sigma -= s_gamma * (dprior/100.0) * (targetAcceptance);
2952 }
2953
2954 sigma = (sigma > dprior ? dprior : sigma);
2955 sigma = (sigma < DBL_MIN ? DBL_MIN : sigma);
2956 LALInferenceSetVariable(args, tmpname, &sigma);
2957 }
2958 }
2959
2960 /* Make sure we don't do this again until we take another adaptable step.*/
2961 adaptableStep = 0;
2962 LALInferenceSetVariable(args, "adaptableStep", &adaptableStep);
2963}
2964
2965
2966
2967
2968/**
2969 * Setup all clustered-KDE proposals with samples read from file.
2970 *
2971 * Constructed clustered-KDE proposals from all sample lists provided in
2972 * files given on the command line.
2973 * @param thread The LALInferenceThreadState to get command line options from and to the proposal cycle of.
2974 * @param input The input file pointer
2975 * @param burnin The number of burn-in points
2976 * @param weight The weight?
2977 * @param ptmcmc Flag to determine if using parallel tempering MCMC
2978 */
2981 INT4 j=0, k=0;
2982
2983 INT4 cyclic_reflective = LALInferenceGetINT4Variable(thread->proposalArgs, "cyclic_reflective_kde");
2984
2986
2987 INT4 nInSamps;
2988 INT4 nCols;
2989 REAL8 *sampleArray;
2990
2991 if (ptmcmc)
2993
2994 char params[128][VARNAME_MAX];
2995 LALInferenceReadAsciiHeader(input, params, &nCols);
2996
2997 LALInferenceVariables *backwardClusterParams = XLALCalloc(1, sizeof(LALInferenceVariables));
2998
2999 /* Only cluster parameters that are being sampled */
3000 INT4 nValidCols=0;
3001 INT4 *validCols = XLALCalloc(nCols, sizeof(INT4));
3002 for (j=0; j<nCols; j++)
3003 validCols[j] = 0;
3004
3005 INT4 logl_idx = 0;
3006 for (j=0; j<nCols; j++) {
3007 if (!strcmp("logl", params[j])) {
3008 logl_idx = j;
3009 continue;
3010 }
3011
3012 char* internal_param_name = XLALCalloc(512, sizeof(char));
3014
3015 for (item = thread->currentParams->head; item; item = item->next) {
3016 if (!strcmp(item->name, internal_param_name) &&
3018 nValidCols++;
3019 validCols[j] = 1;
3020 LALInferenceAddVariable(backwardClusterParams, item->name, item->value, item->type, item->vary);
3021 break;
3022 }
3023 }
3024 }
3025
3026 /* LALInferenceAddVariable() builds the array backwards, so reverse it. */
3027 LALInferenceVariables *clusterParams = XLALCalloc(1, sizeof(LALInferenceVariables));
3028
3029 for (item = backwardClusterParams->head; item; item = item->next)
3030 LALInferenceAddVariable(clusterParams, item->name, item->value, item->type, item->vary);
3031
3032 /* Burn in samples and parse the remainder */
3033 if (ptmcmc)
3034 LALInferenceBurninPTMCMC(input, logl_idx, nValidCols);
3035 else
3036 LALInferenceBurninStream(input, burnin);
3037
3038 sampleArray = LALInferenceParseDelimitedAscii(input, nCols, validCols, &nInSamps);
3039
3040 /* Downsample PTMCMC file to have independent samples */
3041 if (ptmcmc) {
3042 REAL8 acl_real8 = LALInferenceComputeMaxAutoCorrLen(sampleArray, nInSamps, nValidCols);
3043 INT4 acl;
3044 if (acl_real8 == INFINITY)
3045 acl = INT_MAX;
3046 else if (acl_real8 < 1.0)
3047 acl = 1;
3048 else
3049 acl = (INT4)acl_real8;
3050
3051 INT4 downsampled_size = ceil((REAL8)nInSamps/acl);
3052 REAL8 *downsampled_array = (REAL8 *)XLALCalloc(downsampled_size * nValidCols, sizeof(REAL8));
3053 printf("Downsampling to achieve %i samples.\n", downsampled_size);
3054 for (k=0; k < downsampled_size; k++) {
3055 for (j=0; j < nValidCols; j++)
3056 downsampled_array[k*nValidCols + j] = sampleArray[k*nValidCols*acl + j];
3057 }
3058 XLALFree(sampleArray);
3059 sampleArray = downsampled_array;
3060 nInSamps = downsampled_size;
3061 }
3062
3063 /* Build the KDE estimate and add to the KDE proposal set */
3064 INT4 ntrials = 50; // Number of trials at fixed-k to find optimal BIC
3065 LALInferenceInitClusteredKDEProposal(thread, kde, sampleArray, nInSamps, clusterParams, clusteredKDEProposalName, weight, LALInferenceOptimizedKmeans, cyclic_reflective, ntrials);
3066
3067 /* If kmeans construction failed, halt the run */
3068 if (!kde->kmeans) {
3069 fprintf(stderr, "\nERROR: Couldn't build kmeans clustering from the file specified.\n");
3070 XLALFree(kde);
3071 exit(-1);
3072 }
3073
3075
3076 LALInferenceClearVariables(backwardClusterParams);
3077 XLALFree(backwardClusterParams);
3078 XLALFree(sampleArray);
3079}
3080
3081
3082/**
3083 * Initialize a clustered-KDE proposal.
3084 *
3085 * Estimates the underlying distribution of a set of points with a clustered kernel density estimate
3086 * and constructs a jump proposal from the estimated distribution.
3087 * @param thread The current LALInferenceThreadState.
3088 * @param[out] kde An empty proposal structure to populate with the clustered-KDE estimate.
3089 * @param[in] array The data to estimate the underlying distribution of.
3090 * @param[in] nSamps Number of samples contained in \a array.
3091 * @param[in] params The parameters contained in \a array.
3092 * @param[in] name The name of the proposal being constructed.
3093 * @param[in] weight The relative weight this proposal is to have against other KDE proposals.
3094 * @param[in] cluster_method A pointer to the clustering function to be used.
3095 * @param[in] cyclic_reflective Flag to check for cyclic/reflective bounds.
3096 * @param[in] ntrials Number of kmeans attempts at fixed k to find optimal BIC.
3097 */
3100 REAL8 *array,
3101 INT4 nSamps,
3103 const char *name,
3104 REAL8 weight,
3105 LALInferenceKmeans* (*cluster_method)(gsl_matrix*, INT4, gsl_rng*),
3106 INT4 cyclic_reflective,
3107 INT4 ntrials) {
3108 INT4 dim;
3109 INT4 ndraws = 1000;
3110 gsl_matrix_view mview;
3111 char outp_name[256];
3112 char outp_draws_name[256];
3113
3114 strcpy(kde->name, name);
3116
3117 /* If kmeans is already assigned, assume it was calculated elsewhere */
3118 if (!kde->kmeans) {
3119 mview = gsl_matrix_view_array(array, nSamps, dim);
3120 kde->kmeans = (*cluster_method)(&mview.matrix, ntrials, thread->GSLrandom);
3121 }
3122
3123 /* Return if kmeans setup failed */
3124 if (!kde->kmeans)
3125 return;
3126
3127 kde->dimension = kde->kmeans->dim;
3128 kde->params = params;
3129
3130 kde->weight = weight;
3131 kde->next = NULL;
3132
3133 /* Selectivey impose bounds on KDEs */
3134 LALInferenceKmeansImposeBounds(kde->kmeans, params, thread->priorArgs, cyclic_reflective);
3135
3136 /* Print out clustered samples, assignments, and PDF values if requested */
3137 if (LALInferenceGetINT4Variable(thread->proposalArgs, "verbose")) {
3138 printf("Thread %i found %i clusters.\n", thread->id, kde->kmeans->k);
3139
3140 sprintf(outp_name, "clustered_samples.%2.2d", thread->id);
3141 sprintf(outp_draws_name, "clustered_draws.%2.2d", thread->id);
3142
3143 LALInferenceDumpClusteredKDE(kde, outp_name, array);
3144 LALInferenceDumpClusteredKDEDraws(kde, outp_draws_name, ndraws);
3145 }
3146}
3147
3148
3149/**
3150 * Dump draws from a KDE to file.
3151 *
3152 * Print out the samples used to estimate the distribution, along with their
3153 * cluster assignments, and the PDF evaluated at each sample.
3154 * @param[in] kde The clustered KDE to dump the info of.
3155 * @param[in] outp_name The name of the output file.
3156 * @param[in] array The array of samples used for the KDE (it only stores a whitened version).
3157 */
3159 FILE *outp;
3160 REAL8 PDF;
3161 INT4 i, j;
3162
3163 outp = fopen(outp_name, "w");
3165 fprintf(outp, "cluster\tweight\tPDF\n");
3166
3167 for (i=0; i<kde->kmeans->npts; i++) {
3168 PDF = LALInferenceKmeansPDF(kde->kmeans, array + i*kde->dimension);
3169 for (j=0; j<kde->dimension; j++)
3170 fprintf(outp, "%g\t", array[i*kde->dimension + j]);
3171 fprintf(outp, "%i\t%f\t%g\n", kde->kmeans->assignments[i], kde->kmeans->weights[kde->kmeans->assignments[i]], PDF);
3172 }
3173 fclose(outp);
3174}
3175
3176
3177/**
3178 * Dump clustered KDE information to file.
3179 *
3180 * Dump a requested number of draws from a clustered-KDE to file,
3181 * along with the value of the PDF at each point.
3182 * @param[in] kde The clustered-KDE proposal to draw from.
3183 * @param[in] outp_name The name of the file to write to.
3184 * @param[in] nSamps The number of draws to write.
3185 */
3187 FILE *outp;
3188 INT4 i, j;
3189 REAL8 *draw, PDF;
3190
3191 outp = fopen(outp_name, "w");
3193 fprintf(outp, "PDF\n");
3194
3195 for (i=0; i<nSamps; i++) {
3196 draw = LALInferenceKmeansDraw(kde->kmeans);
3197 PDF = LALInferenceKmeansPDF(kde->kmeans, draw);
3198 for (j=0; j<kde->dimension; j++)
3199 fprintf(outp, "%g\t", draw[j]);
3200 fprintf(outp, "%g\n", PDF);
3201 XLALFree(draw);
3202 }
3203 fclose(outp);
3204}
3205
3206
3207/**
3208 * Add a KDE proposal to the KDE proposal set.
3209 *
3210 * If other KDE proposals already exist, the provided KDE is appended to the list, otherwise it is added
3211 * as the first of such proposals.
3212 * @param propArgs The proposal arguments to be added to.
3213 * @param[in] kde The proposal to be added to \a thread->cycle.
3214 */
3216 /* If proposal doesn't already exist, add to proposal args */
3219
3220 /* If proposals already exist, add to the end */
3221 } else {
3223 LALInferenceClusteredKDE *old_kde = NULL;
3224
3225 /* If the first proposal has the same name, replace it */
3226 if (!strcmp(existing_kde->name, kde->name)) {
3227 old_kde = existing_kde;
3228 kde->next = existing_kde->next;
3229 LALInferenceSetVariable(propArgs, clusteredKDEProposalName, (void *)&kde);
3230 } else {
3231 while (existing_kde->next != NULL) {
3232 /* Replace proposal with the same name if found */
3233 if (!strcmp(existing_kde->next->name, kde->name)) {
3234 old_kde = existing_kde->next;
3235 kde->next = old_kde->next;
3236 existing_kde->next = kde;
3237 break;
3238 }
3239 existing_kde = existing_kde->next;
3240 }
3241
3242 /* If a proposal was not replaced, add the proposal to the end of the list */
3243 existing_kde->next=kde;
3244 }
3245
3247 }
3248
3249 return;
3250}
3251
3252
3253/**
3254 * Destroy an existing clustered-KDE proposal.
3255 *
3256 * Convenience function for freeing a clustered KDE proposal that
3257 * already exists. This is particularly useful for a proposal that
3258 * is updated during a run.
3259 * @param proposal The proposal to be destroyed.
3260 */
3262 if (proposal != NULL) {
3265 XLALFree(proposal->params);
3266 }
3267 return;
3268}
3269
3270
3271/**
3272 * Setup a clustered-KDE proposal from the differential evolution buffer.
3273 *
3274 * Reads the samples currently in the differential evolution buffer and construct a
3275 * jump proposal from its clustered kernel density estimate.
3276 * @param thread The LALInferenceThreadState to get the buffer from and add the proposal to.
3277 */
3279 INT4 i;
3280
3281 /* If ACL can be estimated, thin DE buffer to only have independent samples */
3282 REAL8 bufferSize = (REAL8) thread->differentialPointsLength;
3283 REAL8 effSampleSize = (REAL8) LALInferenceComputeEffectiveSampleSize(thread);
3284
3285 /* Correlations wont effect the proposal much, so floor is taken instead of ceil
3286 * when determining the step size */
3287 INT4 step = 1;
3288 if (effSampleSize > 0)
3289 step = (INT4) floor(bufferSize/effSampleSize);
3290
3291 if (step == 0)
3292 step = 1;
3293 INT4 nPoints = (INT4) ceil(bufferSize/(REAL8)step);
3294
3295 /* Get points to be clustered from the differential evolution buffer. */
3297 REAL8** DEsamples = (REAL8**) XLALCalloc(nPoints, sizeof(REAL8*));
3298 REAL8* temp = (REAL8*) XLALCalloc(nPoints * nPar, sizeof(REAL8));
3299 for (i=0; i < nPoints; i++)
3300 DEsamples[i] = temp + (i*nPar);
3301
3302 LALInferenceThinnedBufferToArray(thread, DEsamples, step);
3303
3304 /* Check if imposing cyclic reflective bounds */
3305 INT4 cyclic_reflective = LALInferenceGetINT4Variable(thread->proposalArgs, "cyclic_reflective_kde");
3306
3307 INT4 ntrials = 5;
3308 LALInferenceSetupClusteredKDEProposalFromRun(thread, DEsamples[0], nPoints, cyclic_reflective, ntrials);
3309
3310 /* The proposal copies the data, so the local array can be freed */
3311 XLALFree(temp);
3312 XLALFree(DEsamples);
3313}
3314
3315/**
3316 * Setup a clustered-KDE proposal from the parameters in a run.
3317 *
3318 * Reads the samples currently in the differential evolution buffer and construct a
3319 * jump proposal from its clustered kernel density estimate.
3320 * @param thread The LALInferenceThreadState to get the buffer from and add the proposal to.
3321 * @param samples The samples to estimate the distribution of. Column order expected to match
3322 * the order in \a thread->currentParams.
3323 * @param size Number of samples in \a samples.
3324 * @param cyclic_reflective Flag to check for cyclic/reflective bounds.
3325 * @param ntrials Number of tirals at fixed-k to find optimal BIC
3326 */
3328 REAL8 weight=2.;
3329
3330 /* Keep track of clustered parameter names */
3331 LALInferenceVariables *backwardClusterParams = XLALCalloc(1, sizeof(LALInferenceVariables));
3332 LALInferenceVariables *clusterParams = XLALCalloc(1, sizeof(LALInferenceVariables));
3334 for (item = thread->currentParams->head; item; item = item->next)
3336 LALInferenceAddVariable(backwardClusterParams, item->name, item->value, item->type, item->vary);
3337 for (item = backwardClusterParams->head; item; item = item->next)
3338 LALInferenceAddVariable(clusterParams, item->name, item->value, item->type, item->vary);
3339
3340 /* Build the proposal */
3342 LALInferenceInitClusteredKDEProposal(thread, proposal, samples, size, clusterParams, clusteredKDEProposalName, weight, LALInferenceOptimizedKmeans, cyclic_reflective, ntrials);
3343
3344 /* Only add the kmeans was successfully setup */
3345 if (proposal->kmeans)
3347 else {
3348 LALInferenceClearVariables(clusterParams);
3349 XLALFree(clusterParams);
3350 XLALFree(proposal);
3351 }
3352
3353 LALInferenceClearVariables(backwardClusterParams);
3354 XLALFree(backwardClusterParams);
3355}
3356
3357
3358/**
3359 * A proposal based on the clustered kernal density estimate of a set of samples.
3360 *
3361 * Proposes samples from the estimated distribution of a collection of points.
3362 * The distribution is estimated with a clustered kernel density estimator. This
3363 * proposal is added to the proposal cycle with a specified weight, and in turn
3364 * chooses at random a KDE-estimate from a linked list.
3365 * @param thread The current LALInferenceThreadState.
3366 * @param currentParams The current parameters.
3367 * @param[out] proposedParams The proposed parameters.
3368 * @return proposal_ratio The (log) proposal ratio for maintaining detailed balance
3369 */
3371 REAL8 logPropRatio;
3372
3373 logPropRatio = LALInferenceStoredClusteredKDEProposal(thread, currentParams, proposedParams, NULL);
3374
3375 return logPropRatio;
3376}
3377
3378/**
3379 * An interface to the KDE proposal that avoids a KDE evaluation if possible.
3380 *
3381 * If the value of the KDE at the current location is known, use it. Otherwise
3382 * calculate and return.
3383 * @param thread The current LALInferenceThreadState.
3384 * @param currentParams The current parameters.
3385 * @param[out] proposedParams The proposed parameters.
3386 * @param propDensity If input is not NULL or >-INFINITY, assume this is the
3387 * proposal density at \a currentParams, otherwise
3388 * calculate. It is then replaced with the proposal
3389 * density at \a proposedParams.
3390 * @return proposal_ratio The (log) proposal ratio for maintaining detailed balance
3391 */
3393 REAL8 cumulativeWeight, totalWeight;
3394 REAL8 logPropRatio = 0.0;
3395
3397 LALInferenceVariables *propArgs = thread->proposalArgs;
3398
3400 LALInferenceClearVariables(proposedParams);
3401 return logPropRatio; /* Quit now, since there is no proposal to call */
3402 }
3403
3405
3406 /* Clustered KDE estimates are stored in a linked list, with possibly different weights */
3408
3409 totalWeight = 0.;
3410 LALInferenceClusteredKDE *kde = kdes;
3411 while (kde!=NULL) {
3412 totalWeight += kde->weight;
3413 kde = kde->next;
3414 }
3415
3416 /* If multiple KDE estimates exists, draw one at random */
3417 REAL8 randomDraw = gsl_rng_uniform(thread->GSLrandom);
3418
3419 kde = kdes;
3420 cumulativeWeight = kde->weight;
3421 while(cumulativeWeight/totalWeight < randomDraw) {
3422 kde = kde->next;
3423 cumulativeWeight += kde->weight;
3424 }
3425
3426 /* Draw a sample and fill the proposedParams variable with the parameters described by the KDE */
3427 REAL8 *current = XLALCalloc(kde->dimension, sizeof(REAL8));
3428 REAL8 *proposed = LALInferenceKmeansDraw(kde->kmeans);
3429
3430 INT4 i=0;
3431 for (item = kde->params->head; item; item = item->next) {
3433 current[i] = *(REAL8 *) LALInferenceGetVariable(currentParams, item->name);
3434 LALInferenceSetVariable(proposedParams, item->name, &(proposed[i]));
3435 i++;
3436 }
3437 }
3438
3439 /* Calculate the proposal ratio */
3440 REAL8 logCurrentP;
3441 if (propDensity == NULL || *propDensity == -INFINITY)
3442 logCurrentP = LALInferenceKmeansPDF(kde->kmeans, current);
3443 else
3444 logCurrentP = *propDensity;
3445
3446 REAL8 logProposedP = LALInferenceKmeansPDF(kde->kmeans, proposed);
3447
3448 logPropRatio = logCurrentP - logProposedP;
3449
3450 if (propDensity != NULL)
3451 *propDensity = logProposedP;
3452
3453 XLALFree(current);
3454 XLALFree(proposed);
3455
3456 return logPropRatio;
3457}
3458
3459
3460/**
3461 * A wrapper for the KDE proposal that doesn't store KDE evaluations.
3462 *
3463 */
3464
3465
3466/**
3467 * Compute the maximum ACL from the differential evolution buffer.
3468 *
3469 * Given the current differential evolution buffer, the maximum
3470 * one-dimensional autocorrelation length is found.
3471 * @param thread The thread state containing the differential evolution buffer.
3472 * @param maxACL UNDOCUMENTED
3473*/
3475 INT4 nPar, nPoints, nSkip;
3476 INT4 i;
3477 REAL8** DEarray;
3478 REAL8* temp;
3479 REAL8 max_acl;
3480
3482 nPoints = thread->differentialPointsLength;
3483
3484 /* Determine the number of iterations between each entry in the DE buffer */
3485 nSkip = thread->differentialPointsSkip;
3486
3487 /* Prepare 2D array for DE points */
3488 DEarray = (REAL8**) XLALCalloc(nPoints, sizeof(REAL8*));
3489 temp = (REAL8*) XLALCalloc(nPoints * nPar, sizeof(REAL8));
3490 for (i=0; i < nPoints; i++)
3491 DEarray[i] = temp + (i*nPar);
3492
3493 LALInferenceBufferToArray(thread, DEarray);
3494 max_acl = nSkip * LALInferenceComputeMaxAutoCorrLen(DEarray[nPoints/2], nPoints-nPoints/2, nPar);
3495
3496 if (max_acl == INFINITY)
3497 max_acl = INT_MAX;
3498 else if (max_acl < 1.0)
3499 max_acl = 1.0;
3500
3501 *maxACL = (INT4)max_acl;
3502 XLALFree(temp);
3503 XLALFree(DEarray);
3504}
3505
3506/**
3507 * Compute the maximum single-parameter autocorrelation length. Each
3508 * parameter's ACL is the smallest s such that
3509 *
3510 * 1 + 2*ACF(1) + 2*ACF(2) + ... + 2*ACF(M*s) < s,
3511 *
3512 * The parameter M controls the length of the window over which we sum
3513 * the ACF to estimate the ACL; there must be at least M*ACL samples
3514 * summed. This ensures that we obtain a reliable estimate of the ACL
3515 * by incorporating lags that are much longer that the estimated ACL.
3516 *
3517 * The maximum window length is also restricted to be N/K as a safety
3518 * precaution against relying on data near the extreme of the lags in
3519 * the ACF, where there is a lot of noise.
3520 *
3521 * By default, safe parameters are M = 5, K = 2.
3522 *
3523 * If no estimate can be obtained, then return Infinity.
3524 *
3525 * @param array Array with rows containing samples.
3526 * @param nPoints UNDOCUMENTED
3527 * @param nPar UNDOCUMENTED
3528 * @return The maximum one-dimensional autocorrelation length
3529*/
3531 INT4 M=5, K=2;
3532
3533 REAL8 mean, ACL, ACF, maxACL=0;
3534 INT4 par=0, lag=0, i=0, imax;
3535 REAL8 cumACF, s;
3536
3537 if (nPoints > 1) {
3538 imax = nPoints/K;
3539
3540 for (par=0; par<nPar; par++) {
3541 mean = gsl_stats_mean(array+par, nPar, nPoints);
3542 for (i=0; i<nPoints; i++)
3543 array[i*nPar + par] -= mean;
3544
3545 lag=1;
3546 ACL=1.0;
3547 ACF=1.0;
3548 s=1.0/(REAL8)M;
3549 cumACF=1.0;
3550 while (cumACF >= s) {
3551 ACF = gsl_stats_correlation(array + par, nPar, array + lag*nPar + par, nPar, nPoints-lag);
3552 cumACF += 2.0 * ACF;
3553 lag++;
3554 s = (REAL8)lag/(REAL8)M;
3555 if (lag > imax) {
3556 return INFINITY; /* Short circuit: this parameter has indeterminate ACL */
3557 }
3558 }
3559 ACL = cumACF;
3560
3561 if (ACL > maxACL)
3562 maxACL = ACL;
3563 else if (gsl_isnan(ACL))
3564 return INFINITY; /* Short circuit: this parameter has indeterminate ACL */
3565
3566 for (i=0; i<nPoints; i++)
3567 array[i*nPar + par] += mean;
3568 }
3569 } else {
3570 maxACL = INFINITY;
3571 }
3572
3573 return maxACL;
3574}
3575
3576/**
3577 * Update the estimatate of the autocorrelation length.
3578 *
3579 * @param thread The current LALInferenceThreadState.
3580*/
3582 // Calculate ACL with latter half of data to avoid ACL overestimation from chain equilibrating after adaptation
3583 INT4 acl;
3584
3586
3587 LALInferenceSetVariable(thread->proposalArgs, "acl", &acl);
3588}
3589
3590/**
3591 * Determine the effective sample size based on the DE buffer.
3592 *
3593 * Compute the number of independent samples in the differential evolution
3594 * buffer.
3595 * @param thread The current LALInferenceThreadState.
3596 */
3598 /* Update the ACL estimate, assuming a thinned DE buffer if ACL isn't available */
3599 INT4 acl = 1;
3600 if (LALInferenceCheckVariable(thread->proposalArgs, "acl")) {
3602 acl = LALInferenceGetINT4Variable(thread->proposalArgs, "acl");
3603 }
3604
3605 /* Estimate the total number of samples post-burnin based on samples in DE buffer */
3606 INT4 nPoints = thread->differentialPointsLength * thread->differentialPointsSkip;
3607 INT4 iEff = nPoints/acl;
3608 return iEff;
3609}
3610
3611
3613 fprintf(fp, "proposal\t");
3616 fprintf(fp, "prop_ratio\taccepted\t");
3617 fprintf(fp, "\n");
3618 return 0;
3619}
3620
3622 fprintf(fp, "%s\t", cycle->proposals[cycle->counter]->name);
3624 LALInferencePrintSampleNonFixed(fp, theta_prime);
3625 fprintf(fp, "%9.5f\t", exp(logPropRatio));
3626 fprintf(fp, "%d\t", accepted);
3627 fprintf(fp, "\n");
3628 return;
3629}
3630
3632 INT4 nifo = LALInferenceGetINT4Variable(thread->proposalArgs, "nDet");
3633 LALInferenceIFOData *det=NULL;
3635
3636 for(det=thread->parent->data;det;det=det->next)
3637 {
3638 UINT4 i;
3639
3640 char ampName[VARNAME_MAX];
3641 char phaseName[VARNAME_MAX];
3642 REAL8 dummy;
3643
3644 REAL8 ampWidth;
3645 REAL8 phaseWidth;
3646 UINT4 nspl = LALInferenceGetUINT4Variable(proposedParams, "spcal_npts");
3647 for (i = 0; i < nspl; i++) {
3648 if((VARNAME_MAX <= snprintf(ampName, VARNAME_MAX, "%s_spcal_amp_%i", det->name, i)))
3649 {
3650 fprintf(stderr,"variable name too long\n"); exit(1);
3651 }
3652 if((VARNAME_MAX <= snprintf(phaseName, VARNAME_MAX, "%s_spcal_phase_%i", det->name, i)))
3653 {
3654 fprintf(stderr,"variable name too long\n"); exit(1);
3655 }
3656
3657 LALInferenceGetGaussianPrior(thread->priorArgs, ampName, &dummy, &ampWidth);
3658 REAL8 amp = LALInferenceGetREAL8Variable(proposedParams, ampName);
3659 amp += ampWidth*gsl_ran_ugaussian(thread->GSLrandom)/sqrt(nifo*(REAL8)nspl);
3660 LALInferenceSetREAL8Variable(proposedParams, ampName, amp);
3661
3662 LALInferenceGetGaussianPrior(thread->priorArgs, phaseName, &dummy, &phaseWidth);
3663 REAL8 ph = LALInferenceGetREAL8Variable(proposedParams, phaseName);
3664 ph += phaseWidth*gsl_ran_ugaussian(thread->GSLrandom)/sqrt(nifo*(REAL8)nspl);
3665 LALInferenceSetREAL8Variable(proposedParams, phaseName, ph);
3666 }
3667 };
3668
3669 return 0.0;
3670}
void LALInferenceKmeansImposeBounds(LALInferenceKmeans *kmeans, LALInferenceVariables *params, LALInferenceVariables *priorArgs, INT4 cyclic_reflective)
Impose boundaries on individual KDEs.
REAL8 LALInferenceKmeansPDF(LALInferenceKmeans *kmeans, REAL8 *pt)
Evaluate the estimated (log) PDF from a clustered-KDE at a point.
LALInferenceKmeans * LALInferenceOptimizedKmeans(gsl_matrix *data, INT4 ntrials, gsl_rng *rng)
Kmeans cluster data, maximizing BIC assuming BIC(k) is concave-down.
void LALInferenceKmeansDestroy(LALInferenceKmeans *kmeans)
Destroy a kmeans instance.
REAL8 * LALInferenceKmeansDraw(LALInferenceKmeans *kmeans)
Draw a sample from a kmeans-KDE estimate of a distribution.
LALInferenceModel * LALInferenceInitCBCModel(LALInferenceRunState *state)
Initialise state variables needed for LALInferenceNest or LALInferenceMCMC to run on a CBC signal.
static REAL8 mean(REAL8 *array, int N)
static INT4 numDetectorsUniquePositions(LALInferenceIFOData *data)
static void UpdateWaveletSum(LALInferenceThreadState *thread, LALInferenceVariables *proposedParams, gsl_matrix *glitchFD, INT4 ifo, INT4 n, INT4 flag)
static REAL8 draw_colatitude(LALInferenceThreadState *thread, const char *name)
REAL8 LALInferenceDifferentialEvolutionNames(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams, const char **names)
static void MorletDiagonalFisherMatrix(REAL8Vector *params, REAL8Vector *sigmas)
static void cart_to_sph(const REAL8 cart[3], REAL8 *lat, REAL8 *longi)
static void reflected_position_and_time(LALInferenceThreadState *thread, const REAL8 ra, const REAL8 dec, const REAL8 oldTime, REAL8 *newRA, REAL8 *newDec, REAL8 *newTime)
static void reflected_extrinsic_parameters(LALInferenceThreadState *thread, const REAL8 ra, const REAL8 dec, const REAL8 baryTime, const REAL8 dist, const REAL8 iota, const REAL8 psi, REAL8 *newRA, REAL8 *newDec, REAL8 *newTime, REAL8 *newDist, REAL8 *newIota, REAL8 *newPsi)
static REAL8 draw_logdistance(LALInferenceThreadState *thread)
const char *const distanceLikelihoodProposalName
static const char * intrinsicNames[]
static void unit_vector(REAL8 v[3], const REAL8 w[3])
static REAL8 norm(const REAL8 x[3])
static const char * extrinsicNames[]
const char *const splineCalibrationProposalName
@ USES_LOG_DISTANCE_VARIABLE
@ USES_DISTANCE_VARIABLE
static void cross_product(REAL8 x[3], const REAL8 y[3], const REAL8 z[3])
static REAL8 draw_flat(LALInferenceThreadState *thread, const char *name)
static void get_detectors(LALInferenceIFOData *data, LALDetector **detectors)
static REAL8 evaluate_morlet_proposal(LALInferenceThreadState *thread, LALInferenceVariables *proposedParams, INT4 ifo, INT4 k)
static REAL8 draw_chirp(LALInferenceThreadState *thread)
static REAL8 draw_distance(LALInferenceThreadState *thread)
static void project_along(REAL8 vproj[3], const REAL8 v[3], const REAL8 w[3])
static void vsub(REAL8 diff[3], const REAL8 w[3], const REAL8 v[3])
static void sph_to_cart(REAL8 cart[3], const REAL8 lat, const REAL8 longi)
static REAL8 draw_dec(LALInferenceThreadState *thread)
static INT4 same_detector_location(LALDetector *d1, LALDetector *d2)
static REAL8 approxLogPrior(LALInferenceVariables *params)
REAL8 LALInferencePolarizationPhaseJump(UNUSED LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
static REAL8 dot(const REAL8 v[3], const REAL8 w[3])
static void reflect_plane(REAL8 pref[3], const REAL8 p[3], const REAL8 x[3], const REAL8 y[3], const REAL8 z[3])
static REAL8 glitchAmplitudeDraw(REAL8 Q, REAL8 f, gsl_rng *r)
#define max(a, b)
LALInferenceVariables currentParams
ProcessParamsTable * ppt
int j
ProcessParamsTable * ptr
int k
static double double delta
#define fprintf
const char *const name
const char * names[]
int s
int l
double e
double theta
const double sx
const double Q
const double pn
const double sy
const double u
const double ny
const double a2
const double w
sigmaKerr data[0]
const double nx
void XLALComputeDetAMResponse(double *fplus, double *fcross, const REAL4 D[3][3], const double ra, const double dec, const double psi, const double gmst)
#define LAL_2_SQRTPI
#define LAL_C_SI
#define LAL_E
#define LAL_PI
#define LAL_TWOPI
#define LAL_SQRT1_2
double REAL8
int64_t INT8
uint32_t UINT4
int32_t INT4
float REAL4
INT4 LALInferenceBufferToArray(LALInferenceThreadState *thread, REAL8 **DEarray)
INT4 LALInferenceThinnedBufferToArray(LALInferenceThreadState *thread, REAL8 **DEarray, INT4 step)
LALInferenceVariable buffer to array and vica versa.
void LALInferencePrintSampleNonFixed(FILE *fp, LALInferenceVariables *sample)
Output only non-fixed parameters.
Definition: LALInference.c:948
void LALInferenceReadAsciiHeader(FILE *input, char params[][VARNAME_MAX], INT4 *nCols)
Read column names from an ASCII file.
LALInferenceVariableType LALInferenceGetVariableType(const LALInferenceVariables *vars, const char *name)
Get the LALInferenceVariableType of the parameter named name in vars.
Definition: LALInference.c:321
void LALInferenceDiscardPTMCMCHeader(FILE *filestream)
Discard the standard header of a PTMCMC chain file.
#define VARNAME_MAX
Definition: LALInference.h:50
int LALInferenceFprintParameterNonFixedHeaders(FILE *out, LALInferenceVariables *params)
Print the parameters which do not vary to a file as a tab-separated ASCII line.
UINT4 LALInferenceGetUINT4Variable(LALInferenceVariables *vars, const char *name)
INT4 LALInferenceGetVariableDimensionNonFixed(LALInferenceVariables *vars)
Get number of dimensions in vars which are not fixed to a certain value.
Definition: LALInference.c:255
void LALInferenceAddREAL8VectorVariable(LALInferenceVariables *vars, const char *name, REAL8Vector *value, LALInferenceParamVaryType vary)
gsl_matrix * LALInferenceGetgslMatrixVariable(LALInferenceVariables *vars, const char *name)
void LALInferenceAddVariable(LALInferenceVariables *vars, const char *name, const void *value, LALInferenceVariableType type, LALInferenceParamVaryType vary)
Add a variable named name to vars with initial value referenced by value.
Definition: LALInference.c:395
REAL8 LALInferenceGetREAL8Variable(LALInferenceVariables *vars, const char *name)
void LALInferenceSetREAL8Variable(LALInferenceVariables *vars, const char *name, REAL8 value)
int LALInferenceCompareVariables(LALInferenceVariables *var1, LALInferenceVariables *var2)
Check for equality in two variables.
void LALInferenceAddstringVariable(LALInferenceVariables *vars, const char *name, const CHAR *value, LALInferenceParamVaryType vary)
LALInferenceVariableItem * LALInferenceGetItemNr(LALInferenceVariables *vars, int idx)
Return the list node for the idx-th item - do not rely on this Indexing starts at 1.
Definition: LALInference.c:207
void LALInferenceCopyVariables(LALInferenceVariables *origin, LALInferenceVariables *target)
Deep copy the variables from one to another LALInferenceVariables structure.
Definition: LALInference.c:558
void LALInferenceBurninStream(FILE *filestream, INT4 burnin)
Burn-in a generic ASCII stream.
void * LALInferenceGetVariable(const LALInferenceVariables *vars, const char *name)
Return a pointer to the memory the variable vars is stored in specified by name User must cast this p...
Definition: LALInference.c:238
REAL8Vector * LALInferenceGetREAL8VectorVariable(LALInferenceVariables *vars, const char *name)
void LALInferenceAddREAL8Variable(LALInferenceVariables *vars, const char *name, REAL8 value, LALInferenceParamVaryType vary)
INT4 LALInferenceGetINT4Variable(LALInferenceVariables *vars, const char *name)
ProcessParamsTable * LALInferenceGetProcParamVal(ProcessParamsTable *procparams, const char *name)
Returns the element of the process params table with "name".
void LALInferenceAddINT4Variable(LALInferenceVariables *vars, const char *name, INT4 value, LALInferenceParamVaryType vary)
REAL8 * LALInferenceParseDelimitedAscii(FILE *input, INT4 nCols, INT4 *wantedCols, INT4 *nLines)
Utility for readling in delimited ASCII files.
void LALInferenceBurninPTMCMC(FILE *filestream, INT4 logl_idx, INT4 nPar)
Determine burnin cycle from delta-logl criteria.
INT4 LALInferenceGetVariableDimension(LALInferenceVariables *vars)
Get number of dimensions in variable vars.
Definition: LALInference.c:250
REAL8(* LALInferenceProposalFunction)(struct tagLALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Jump proposal distribution Computes proposedParams based on currentParams and additional variables st...
Definition: LALInference.h:391
const CHAR * LALInferenceGetstringVariable(LALInferenceVariables *vars, const char *name)
UINT4Vector * LALInferenceGetUINT4VectorVariable(LALInferenceVariables *vars, const char *name)
void LALInferenceClearVariables(LALInferenceVariables *vars)
Delete the variables in this structure.
Definition: LALInference.c:532
INT4 LALInferenceFprintParameterNonFixedHeadersWithSuffix(FILE *out, LALInferenceVariables *params, const char *suffix)
Print the parameters which do not vary to a file as a tab-separated ASCII line, adding the given suff...
int LALInferenceCheckVariable(LALInferenceVariables *vars, const char *name)
Checks for name being present in vars returns 1(==true) or 0.
Definition: LALInference.c:525
void LALInferenceTranslateExternalToInternalParamName(char *outName, const char *inName)
Converts between externally used parameter names and those internal.
void LALInferenceSetVariable(LALInferenceVariables *vars, const char *name, const void *value)
Set a variable named name in vars with a value.
Definition: LALInference.c:352
int LALInferenceCheckVariableNonFixed(LALInferenceVariables *vars, const char *name)
Checks for name being present in vars and having type LINEAR or CIRCULAR.
Definition: LALInference.c:503
@ LALINFERENCE_PARAM_OUTPUT
A parameter that never changes, functions should respect this.
Definition: LALInference.h:131
@ LALINFERENCE_PARAM_FIXED
A parameter that is cyclic, such as an angle between 0 and 2pi.
Definition: LALInference.h:130
@ LALINFERENCE_PARAM_LINEAR
Definition: LALInference.h:128
@ LALINFERENCE_REAL8_t
Definition: LALInference.h:109
@ LALINFERENCE_void_ptr_t
Definition: LALInference.h:119
void LALInferenceNestedSamplingAlgorithm(LALInferenceRunState *runState)
NestedSamplingAlgorithm implements the nested sampling algorithm, see e.g.
int LALInferenceCheckMinMaxPrior(LALInferenceVariables *priorArgs, const char *name)
Check for types of standard prior.
void LALInferenceCyclicReflectiveBound(LALInferenceVariables *parameter, LALInferenceVariables *priorArgs)
Apply cyclic and reflective boundaries to parameter to bring it back within the allowed prior ranges ...
void LALInferenceGetMinMaxPrior(LALInferenceVariables *priorArgs, const char *name, REAL8 *min, REAL8 *max)
Get the minimum and maximum values of the uniform prior from the priorArgs list, given a name.
int LALInferenceCheckGaussianPrior(LALInferenceVariables *priorArgs, const char *name)
Check for a Gaussian prior (with a mean and variance)
void LALInferenceGetGaussianPrior(LALInferenceVariables *priorArgs, const char *name, REAL8 *mu, REAL8 *sigma)
Get the mu and sigma values of the Gaussian prior from the priorArgs list, given a name.
REAL8 logGlitchAmplitudeDensity(REAL8 A, REAL8 Q, REAL8 f)
Return the log Prior for the glitch amplitude.
REAL8 LALInferenceEnsembleWalkNames(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams, const char **names)
REAL8 LALInferenceEnsembleStretchFull(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Ensemble stretch moves - see http://dx.doi.org/10.2140/camcos.2010.5.65.
const char *const singleAdaptProposalName
void LALInferenceDumpClusteredKDEDraws(LALInferenceClusteredKDE *kde, char *outp_name, INT4 nSamps)
Dump clustered KDE information to file.
const char *const polarizationPhaseJumpName
REAL8 LALInferenceSkyRingProposal(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
const char *const GlitchMorletJumpName
void LALInferenceDestroyClusteredKDEProposal(LALInferenceClusteredKDE *proposal)
Destroy an existing clustered-KDE proposal.
REAL8 LALInferenceSkyLocWanderJump(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Jump around by 0.01 radians in angle on the sky.
#define ADAPTSUFFIX
const char *const frequencyBinJumpName
REAL8 LALInferenceEnsembleStretchIntrinsic(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
const char *const ensembleWalkIntrinsicName
const char *const singleProposalName
void LALInferenceSetupAdaptiveProposals(LALInferenceVariables *propArgs, LALInferenceVariables *params)
Setup adaptive proposals.
const char *const extrinsicParamProposalName
const char *const covarianceEigenvectorJumpName
const char *const ensembleStretchExtrinsicName
REAL8 LALInferenceEnsembleWalkFull(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
REAL8 LALInferenceCovarianceEigenvectorJump(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Choose a random covariance matrix eigenvector to jump along.
const char *const cycleArrayName
void LALInferenceUpdateMaxAutoCorrLen(LALInferenceThreadState *thread)
Update the estimatate of the autocorrelation length.
REAL8 LALInferenceGlitchMorletProposal(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
REAL8 LALInferenceExtrinsicParamProposal(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Proposal for the extrinsic parameters.
const char *const drawFlatPriorName
REAL8 LALInferenceGlitchMorletReverseJump(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
const char *const differentialEvolutionExtrinsicName
const char *const GlitchMorletReverseJumpName
REAL8 LALInferenceSplineCalibrationProposal(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Proposes jumps in the spline calibration parameters, if present.
LALInferenceProposalCycle * LALInferenceSetupDefaultInspiralProposalCycle(LALInferenceVariables *propArgs)
A reasonable default proposal.
REAL8 LALInferenceSingleAdaptProposal(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Like LALInferenceSingleProposal() but will use adaptation if the –adapt command-line flag given.
const char *const cycleArrayCounterName
const char *const ensembleStretchIntrinsicName
REAL8 LALInferenceComputeMaxAutoCorrLen(REAL8 *array, INT4 nPoints, INT4 nPar)
Compute the maximum single-parameter autocorrelation length.
const char *const PSDFitJumpName
REAL8 LALInferenceEnsembleWalkIntrinsic(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
REAL8 LALInferenceDifferentialEvolutionExtrinsic(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Perform a differential evolution step on only the extrinsic parameters.
REAL8 LALInferenceStoredClusteredKDEProposal(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams, REAL8 *propDensity)
An interface to the KDE proposal that avoids a KDE evaluation if possible.
const char *const skyRingProposalName
const char *const ensembleWalkExtrinsicName
const char *const clusteredKDEProposalName
REAL8 LALInferenceSingleProposal(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Non-adaptive, sigle-variable update proposal with reasonable widths in each dimension.
void LALInferenceInitClusteredKDEProposal(LALInferenceThreadState *thread, LALInferenceClusteredKDE *kde, REAL8 *array, INT4 nSamps, LALInferenceVariables *params, const char *name, REAL8 weight, LALInferenceKmeans *(*cluster_method)(gsl_matrix *, INT4, gsl_rng *), INT4 cyclic_reflective, INT4 ntrials)
Initialize a clustered-KDE proposal.
REAL8 LALInferenceEnsembleStretchExtrinsic(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
REAL8 LALInferenceDistanceLikelihoodProposal(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Proposal which draws a sample from the distance likelihood function Requires the currentParams to hav...
void LALInferenceUpdateAdaptiveJumps(LALInferenceThreadState *thread, REAL8 targetAcceptance)
Update the adaptive proposal.
void LALInferenceAddProposalToCycle(LALInferenceProposalCycle *cycle, LALInferenceProposal *prop, INT4 weight)
Adds weight copies of the proposal prop to the end of the proposal cycle.
const char *const drawApproxPriorName
INT4 LALInferencePrintProposalTrackingHeader(FILE *fp, LALInferenceVariables *params)
Output proposal tracking header to file *fp.
LALInferenceVariables * LALInferenceParseProposalArgs(LALInferenceRunState *runState)
Go through all logic for deciding proposals to use.
void LALInferenceDeleteProposalCycle(LALInferenceProposalCycle *cycle)
Completely remove the current proposal cycle, freeing the associated memory.
const char *const polarizationCorrPhaseJumpName
void LALInferenceComputeMaxAutoCorrLenFromDE(LALInferenceThreadState *thread, INT4 *maxACL)
A wrapper for the KDE proposal that doesn't store KDE evaluations.
void LALInferenceTrackProposalAcceptance(LALInferenceThreadState *thread)
Update proposal statistics if tracking.
void LALInferenceSetupClusteredKDEProposalFromDEBuffer(LALInferenceThreadState *thread)
Setup a clustered-KDE proposal from the differential evolution buffer.
void LALInferencePrintProposalTracking(FILE *fp, LALInferenceProposalCycle *cycle, LALInferenceVariables *theta, LALInferenceVariables *theta_prime, REAL8 logPropRatio, INT4 accepted)
Output proposal tracking information to file *fp.
void LALInferenceDumpClusteredKDE(LALInferenceClusteredKDE *kde, char *outp_name, REAL8 *array)
Dump draws from a KDE to file.
void LALInferenceZeroProposalStats(LALInferenceProposalCycle *cycle)
INT4 LALInferenceComputeEffectiveSampleSize(LALInferenceThreadState *thread)
Determine the effective sample size based on the DE buffer.
void LALInferenceRegisterProposal(LALInferenceVariables *propArgs, const char *name, INT4 *flag, ProcessParamsTable *command_line)
Use a default flag and a check of the command-line to set proposal flags in proposal args.
void LALInferenceSetupClusteredKDEProposalsFromASCII(LALInferenceThreadState *thread, FILE *input, INT4 burnin, REAL8 weight, INT4 ptmcmc)
Setup all clustered-KDE proposals with samples read from file.
const char *const orbitalPhaseJumpName
REAL8 LALInferenceCorrPolarizationPhaseJump(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Polarization-phase correlation jump.
const char *const skyReflectDetPlaneName
const char *const skyLocWanderJumpName
const char *const ensembleWalkFullName
void LALInferenceAddClusteredKDEProposalToSet(LALInferenceVariables *propArgs, LALInferenceClusteredKDE *kde)
Add a KDE proposal to the KDE proposal set.
void LALInferenceRandomizeProposalCycle(LALInferenceProposalCycle *cycle, gsl_rng *rng)
Randomizes the order of the proposals in the proposal cycle.
REAL8 LALInferenceEnsembleStretchNames(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams, const char **names)
#define PROPOSEDSUFFIX
LALInferenceProposal * LALInferenceInitProposal(LALInferenceProposalFunction func, const char *name)
Creates a new proposal object from the given func pointer and name.
REAL8 LALInferenceDifferentialEvolutionFull(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Differential evolution, on all non-fixed, non-output parameters.
REAL8 LALInferenceDifferentialEvolutionIntrinsic(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Perform differential evolution on only the intrinsic parameters.
REAL8 LALInferenceDrawFlatPrior(LALInferenceThreadState *threadState, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Draws from a flat prior for all variables where are flat prior is specified.
const char *const differentialEvolutionIntrinsicName
LALInferenceProposalCycle * LALInferenceInitProposalCycle(void)
Create a new proposal cycle.
REAL8 LALInferenceCyclicProposal(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Proposes a jump from the next proposal in the proposal cycle.
void LALInferenceSetupGlitchProposal(LALInferenceIFOData *data, LALInferenceVariables *propArgs)
Setup glitch-related stuff.
const char *const cycleArrayLengthName
const char *const ensembleStretchFullName
#define MAX_STRLEN
void LALInferenceSetupClusteredKDEProposalFromRun(LALInferenceThreadState *thread, REAL8 *samples, INT4 size, INT4 cyclic_reflective, INT4 ntrials)
Setup a clustered-KDE proposal from the parameters in a run.
REAL8 LALInferenceDrawApproxPrior(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Draws from an approximation to the true prior.
REAL8 LALInferencePSDFitJump(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
REAL8 LALInferenceSkyReflectDetPlane(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Reflects the sky location through the plane formed by three detectors.
REAL8 LALInferenceEnsembleWalkExtrinsic(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
REAL8 LALInferenceClusteredKDEProposal(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
A proposal based on the clustered kernal density estimate of a set of samples.
#define ACCEPTSUFFIX
REAL8 LALInferenceFrequencyBinJump(LALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Proposal to jump in frequency by one frequency bin.
const char *const differentialEvolutionFullName
const char *const nullProposalName
void * XLALCalloc(size_t m, size_t n)
void * XLALRealloc(void *p, size_t n)
void XLALFree(void *p)
static const INT4 r
static const INT4 a
COORDINATESYSTEM_GEOGRAPHIC
COORDINATESYSTEM_EQUATORIAL
void LALEquatorialToGeographic(LALStatus *stat, SkyPosition *output, SkyPosition *input, LIGOTimeGPS *gpsTime)
void LALGeographicToEquatorial(LALStatus *stat, SkyPosition *output, SkyPosition *input, LIGOTimeGPS *gpsTime)
REAL8 XLALTimeDelayFromEarthCenter(const double detector_earthfixed_xyz_metres[3], double source_right_ascension_radians, double source_declination_radians, const LIGOTimeGPS *gpstime)
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
void XLALDestroyREAL8Vector(REAL8Vector *vector)
#define XLAL_ERROR_REAL8(...)
void XLALError(const char *func, const char *file, int line, int errnum)
XLAL_ENOMEM
XLAL_FAILURE
REAL8 XLALGreenwichMeanSiderealTime(const LIGOTimeGPS *gpstime)
LIGOTimeGPS * XLALGPSSetREAL8(LIGOTimeGPS *epoch, REAL8 t)
REAL8 XLALGPSGetREAL8(const LIGOTimeGPS *epoch)
M
list y
list mu
args
type
DEC
RA
LALDetector detectors[LAL_NUM_DETECTORS]
double alpha
double t0
REAL8 location[3]
Struct to hold clustered-KDE estimates.
LALInferenceKmeans * kmeans
struct tagLALInferenceClusteredKDEProposal * next
LALInferenceVariables * params
Structure to contain IFO data.
Definition: LALInference.h:625
LALDetector * detector
integration limits for overlap integral in F-domain
Definition: LALInference.h:651
char name[DETNAMELEN]
Definition: LALInference.h:626
struct tagLALInferenceIFOData * next
ROQ data.
Definition: LALInference.h:659
Structure for performing the kmeans clustering algorithm on a set of samples.
INT4 npts
Number of points being clustered.
INT4 * assignments
Cluster assignments.
INT4 k
Number of clusters.
REAL8 * weights
Fraction of data points in each cluster.
INT4 dim
Dimension of data.
Structure to constain a model and its parameters.
Definition: LALInference.h:436
LALInferenceVariables * params
Definition: LALInference.h:437
Structure for holding a proposal cycle.
Definition: LALInference.h:526
LALInferenceProposal ** proposals
Definition: LALInference.h:527
char last_proposal_name[VARNAME_MAX]
Counter for cycling through proposals.
Definition: LALInference.h:532
INT4 * order
Array of proposals (one per proposal function)
Definition: LALInference.h:528
INT4 nProposals
Length of cycle.
Definition: LALInference.h:530
Structure for holding a LALInference proposal, along with name and stats.
Definition: LALInference.h:512
char name[VARNAME_MAX]
Definition: LALInference.h:514
LALInferenceProposalFunction func
Definition: LALInference.h:513
Structure containing inference run state This includes pointers to the function types required to run...
Definition: LALInference.h:592
ProcessParamsTable * commandLine
Definition: LALInference.h:593
LALInferenceVariables * proposalArgs
The data from the interferometers.
Definition: LALInference.h:602
LALInferenceVariables * algorithmParams
Any special arguments for the prior function.
Definition: LALInference.h:604
struct tagLALInferenceIFOData * data
Log sample, i.e.
Definition: LALInference.h:601
LALInferenceAlgorithm algorithm
A function that returns a new set of variables for the model.
Definition: LALInference.h:595
Structure containing chain-specific variables.
Definition: LALInference.h:541
LALInferenceVariables * currentParams
Prior boundaries, etc.
Definition: LALInference.h:556
struct tagLALInferenceRunState * parent
Definition: LALInference.h:578
LALInferenceVariables * priorArgs
Stope things such as output arrays.
Definition: LALInference.h:553
LALInferenceModel * model
Cycle of proposals to call.
Definition: LALInference.h:548
INT4 accepted
This should be removed, can be given as an algorithmParams entry.
Definition: LALInference.h:574
size_t differentialPointsLength
Array of points for differential evolution.
Definition: LALInference.h:560
LALInferenceProposalCycle * cycle
The proposal cycle.
Definition: LALInference.h:547
LALInferenceVariables * proposalArgs
Definition: LALInference.h:551
REAL8 temperature
Array containing multiple proposal densities.
Definition: LALInference.h:550
size_t differentialPointsSkip
Size of the differentialPoints memory block (must be >= length of differential points).
Definition: LALInference.h:566
LALInferenceVariables ** differentialPoints
Parameters proposed.
Definition: LALInference.h:559
The LALInferenceVariableItem list node structure This should only be accessed using the accessor func...
Definition: LALInference.h:154
char name[VARNAME_MAX]
Definition: LALInference.h:155
LALInferenceVariableType type
Definition: LALInference.h:157
struct tagVariableItem * next
Definition: LALInference.h:159
LALInferenceParamVaryType vary
Definition: LALInference.h:158
The LALInferenceVariables structure to contain a set of parameters Implemented as a linked list of LA...
Definition: LALInference.h:170
LALInferenceVariableItem * head
Definition: LALInference.h:171
CHAR value[LIGOMETA_VALUE_MAX]
struct tagProcessParamsTable * next
REAL8Sequence * data
REAL8 * data
REAL8 longitude
REAL8 latitude
CoordinateSystem system
UINT4 * data
FILE * fp
enum @1 epoch
double df