26 #include <lal/LALInspiral.h>
27 #include <lal/DetResponse.h>
28 #include <lal/SeqFactories.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>
47 #include <lal/LALStdlib.h>
48 #include <lal/LALInferenceClusteredKDE.h>
49 #include <lal/LALInferenceNestedSampler.h>
52 #define UNUSED __attribute__ ((unused))
97 static 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};
100 static const char *
extrinsicNames[] = {
"rightascension",
"declination",
"cosalpha",
"azimuth",
"polarisation",
"distance",
101 "logdistance",
"time",
"costheta_jn",
"t0",
"theta",
"hrss",
"loghrss", NULL};
106 for (
i = 0;
i < 3;
i++) {
119 for (currentIFO =
data; currentIFO; currentIFO = currentIFO->
next) {
122 for (subsequentIFO = currentIFO->
next; subsequentIFO; subsequentIFO = subsequentIFO->
next) {
130 return nIFO - nCollision;
148 sprintf(offopt,
"--proposal-no-%s",
name);
149 sprintf(onopt,
"--proposal-%s",
name);
160 const char *fname =
"LALInferenceAddProposalToCycle";
168 if (cycle->
order == NULL) {
173 for (
i = cycle->
length; i < cycle->length + weight;
i++) {
193 for (
i = cycle->
length - 1;
i > 0;
i--) {
195 j = gsl_rng_uniform_int(rng,
i+1);
211 cycle = thread->
cycle;
225 REAL8 logPropRatio=-INFINITY;
234 while (proposedParams->
head == NULL);
262 INT4 cyclic_reflective_kde = 0;
266 INT4 singleadapt = 1;
332 INT4 marg_timephi = 0;
346 INT4 analytic_test = 0;
363 INT4 adapting = !noAdapt;
373 INT4 sampling_prior = 0;
399 cyclic_reflective_kde = 1;
429 if (nUniqueDet < 2) {
433 if (nUniqueDet != 3) {
437 if (nUniqueDet >= 3) {
485 const INT4 BIGWEIGHT = 20;
486 const INT4 SMALLWEIGHT = 5;
487 const INT4 TINYWEIGHT = 1;
601 REAL8 logPropRatio, sqrttemp,
mu, sigma;
617 varNr = 1 + gsl_rng_uniform_int(rng, dim);
622 fprintf(stderr,
"Attempting to set non-REAL8 parameter with numerical sigma (in %s, %d)\n",
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__);
640 if (sqrttemp == INFINITY) {
644 *((
REAL8 *)param->
value) += gsl_ran_ugaussian(rng) * sigma;
649 *((
REAL8 *)param->
value) += gsl_ran_ugaussian(rng) * (
max - min);
652 *((
REAL8 *)param->
value) += gsl_ran_ugaussian(rng) * sigma * sqrttemp;
678 REAL8 sigma, big_sigma;
686 if (gsl_ran_ugaussian(GSLrandom) < 1.0e-3)
688 if (gsl_ran_ugaussian(GSLrandom) < 1.0e-4)
694 varNr = 1 + gsl_rng_uniform_int(GSLrandom, dim);
700 if (!strcmp(param->
name,
"eta")) {
702 }
else if (!strcmp(param->
name,
"q")) {
704 }
else if (!strcmp(param->
name,
"chirpmass")) {
706 }
else if (!strcmp(param->
name,
"time")) {
708 }
else if (!strcmp(param->
name,
"t0")) {
710 }
else if (!strcmp(param->
name,
"phase")) {
712 }
else if (!strcmp(param->
name,
"distance")) {
714 }
else if (!strcmp(param->
name,
"declination")) {
716 }
else if (!strcmp(param->
name,
"azimuth")) {
718 }
else if (!strcmp(param->
name,
"cosalpha")) {
720 }
else if (!strcmp(param->
name,
"rightascension")) {
722 }
else if (!strcmp(param->
name,
"polarisation")) {
724 }
else if (!strcmp(param->
name,
"costheta_jn")) {
726 }
else if (!strcmp(param->
name,
"a_spin1")) {
728 }
else if (!strcmp(param->
name,
"a_spin2")) {
731 fprintf(stderr,
"Could not find parameter %s!", param->
name);
735 *(
REAL8 *)param->
value += gsl_ran_ugaussian(GSLrandom)*sigma;
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;
742 *(
REAL8 *)param->
value += gsl_ran_ugaussian(GSLrandom)*big_sigma*sigma*0.01;
749 REAL8 logPropRatio = 0.0;
760 gsl_matrix *eigenvectors;
761 REAL8 jumpSize, tmp, inc;
762 REAL8 logPropRatio = 0.0;
775 i = gsl_rng_uniform_int(rng,
N);
776 jumpSize = sqrt(thread->
temperature * eigenvalues->
data[
i]) * gsl_ran_ugaussian(rng);
779 proposeIterator = proposedParams->
head;
780 if (proposeIterator == NULL) {
781 fprintf(stderr,
"Bad proposed params in %s, line %d\n",
790 inc = jumpSize * gsl_matrix_get(eigenvectors,
j,
i);
798 }
while ((proposeIterator = proposeIterator->
next) != NULL &&
j <
N);
810 REAL8 one_deg = 1.0 / (2.0*M_PI);
811 REAL8 logPropRatio = 0.0;
818 jumpX = sigma * gsl_ran_ugaussian(rng);
819 jumpY = sigma * gsl_ran_ugaussian(rng);
825 newDEC =
DEC + jumpY;
871 const char **names) {
872 size_t i,
N, Ndim, nPts;
874 REAL8 maxScale, Y, logmax, X, scale;
881 const char* local_names[
N];
887 while (item != NULL) {
898 for(Ndim=0,
i=0;
names[
i] != NULL;
i++ ) {
906 for(Ndim=0,
i=0;
names[
i] != NULL;
i++ ) {
914 if (dePts == NULL || nPts <= 1) {
921 i = gsl_rng_uniform_int(thread->
GSLrandom, nPts);
934 logmax = log(maxScale);
935 X = 2.0*logmax*Y - logmax;
938 for (
i = 0;
names[
i] != NULL;
i++) {
944 x = other + scale*(cur-other);
950 if (scale < maxScale && scale > (1.0/maxScale))
951 logPropRatio = log(scale)*((
REAL8)Ndim);
953 logPropRatio = -INFINITY;
984 const char **names) {
987 REAL8 logPropRatio = 0.0;
990 const char* local_names[
N];
996 while (item != NULL) {
1018 size_t sample_size=3;
1023 if (dePts == NULL || nPts <= 1) {
1025 return logPropRatio;
1030 UINT4 indices[sample_size];
1031 UINT4 all_indices[nPts];
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));
1037 double univariate_normals[sample_size];
1038 for(
i=0;
i<sample_size;
i++) univariate_normals[
i] = gsl_ran_ugaussian(thread->
GSLrandom);
1044 REAL8 centre_of_mass=0.0;
1046 for(
i=0;
i<sample_size;
i++)
1051 for(
i=0,
w=0.0;
i<sample_size;
i++)
1061 return logPropRatio;
1067 const char **names) {
1068 size_t i,
j,
N, Ndim, nPts;
1072 REAL8 logPropRatio = 0.0;
1075 const char *localnames[
N];
1079 if (
names == NULL) {
1082 while (item != NULL) {
1084 localnames[
i] = item->
name;
1095 for (Ndim=0,
i=0;
names[
i] != NULL;
i++ ) {
1103 if (dePts == NULL || nPts <= 1)
1104 return logPropRatio;
1109 i = gsl_rng_uniform_int(rng, nPts);
1111 j = gsl_rng_uniform_int(rng, nPts);
1117 const REAL8 modeHoppingFrac = 0.5;
1120 if (gsl_rng_uniform(rng) < modeHoppingFrac) {
1125 scale = 2.38/sqrt(Ndim) * exp(log(0.1) + log(100.0) * gsl_rng_uniform(rng));
1128 for (
i = 0;
names[
i] != NULL;
i++) {
1141 return logPropRatio;
1164 return cbrt(
x*(dmax*dmax*dmax - dmin*dmin*dmin) + dmin*dmin*dmin);
1168 REAL8 logdmin, logdmax;
1172 REAL8 dmin=exp(logdmin);
1173 REAL8 dmax=exp(logdmax);
1177 return log(cbrt(
x*(dmax*dmax*dmax - dmin*dmin*dmin) + dmin*dmin*dmin));
1187 return acos(cos(min) -
x*(cos(min) - cos(
max)));
1197 return asin(
x*(sin(
max) - sin(min)) + sin(min));
1207 return min +
x*(
max - min);
1216 mMin56 = pow(min, 5.0/6.0);
1217 mMax56 = pow(
max, 5.0/6.0);
1219 delta = 1.0/mMin56 - 1.0/mMax56;
1223 return pow(1.0/(1.0/mMin56 -
u), 6.0/5.0);
1231 logP += -11.0/6.0*log(Mc);
1265 INT4 analytic_test,
i;
1266 REAL8 logBackwardJump;
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};
1281 if (analytic_test) {
1293 for (
i = 0; flat_params[
i] != NULL;
i++) {
1340 gsl_matrix_set(
eta,
i,
j,
x);
1346 if (analytic_test) {
1353 return logPropRatio;
1357 REAL8 logPropRatio = 0., tmp = 0.;
1370 return logPropRatio;
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];
1380 return sqrt(
x[0]*
x[0] +
x[1]*
x[1] +
x[2]*
x[2]);
1397 return v[0]*
w[0] + v[1]*
w[1] + v[2]*
w[2];
1407 vproj[0] = what[0]*vdotw;
1408 vproj[1] = what[1]*vdotw;
1409 vproj[2] = what[2]*vdotw;
1413 diff[0] =
w[0] - v[0];
1414 diff[1] =
w[1] - v[1];
1415 diff[2] =
w[2] - v[2];
1420 REAL8 n[3], nhat[3], xy[3], xz[3],
pn[3], pnperp[3];
1447 cart[0] = cos(longi)*cos(lat);
1448 cart[1] = sin(longi)*cos(lat);
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]));
1461 SkyPosition currentEqu, currentGeo, newEqu, newGeo;
1465 REAL8 currentLoc[3], newLoc[3];
1466 REAL8 newGeoLat, newGeoLongi;
1520 *newTime = oldTime + oldDt - newDt;
1529 REAL8 component_min,component_max;
1531 gsl_matrix *glitch_f, *glitch_Q, *glitch_A;
1535 prior -= log(component_max - component_min);
1539 prior -= log(component_max - component_min);
1543 prior -= log(component_max - component_min);
1547 prior -= log(component_max - component_min);
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);
1569 REAL8 SNRPEAK = 5.0;
1580 SNR = 20.0 * SNRPEAK * gsl_rng_uniform(
r);
1582 den = SNR/(SNRPEAK*SNRPEAK) * exp(-SNR/SNRPEAK);
1586 alpha = gsl_rng_uniform(
r);
1587 }
while (
alpha > den);
1589 return SNR/sqrt((PIterm*
Q/f));
1596 INT4 ifo, nifo, timeflag=0;
1597 REAL8 logPropRatio = 0.0;
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];
1606 REAL8 pForward, pReverse;
1636 intpart = (
INT4)gmst;
1637 decpart = gmst - (
REAL8)intpart;
1639 gmst = gmst < 0. ? gmst +
LAL_TWOPI : gmst;
1644 k[0] = cos(gmst-ra) * cos(dec);
1645 k[1] =-sin(gmst-ra) * cos(dec);
1653 IFO = gsl_matrix_alloc(nifo, 3);
1655 for(ifo=0; ifo<nifo; ifo++) {
1658 gsl_matrix_set(IFO, ifo,
i, IFOX[
i]);
1667 i=gsl_rng_uniform_int(rng, nifo);
1668 j=gsl_rng_uniform_int(rng, nifo);
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);
1680 for(
i=0;
i<3;
i++) {
1681 n[
i] = IFO1[
i] - IFO2[
i];
1682 normalize +=
n[
i] *
n[
i];
1684 normalize = 1./sqrt(normalize);
1691 omega =
LAL_TWOPI * gsl_rng_uniform(rng);
1692 cosomega = cos(omega);
1693 sinomega = sin(omega);
1694 c1momega = 1.0 - cosomega;
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];
1712 newDec = asin(kp[2]);
1713 newRA = atan2(kp[1], kp[0]) + gmst;
1725 for(
i=0;
i<3;
i++) {
1729 newTime = tx + baryTime - ty;
1738 newPsi =
LAL_PI * gsl_rng_uniform(rng);
1746 pForward = cos(newDec);
1747 pReverse = cos(dec);
1749 gsl_matrix_free(IFO);
1751 logPropRatio = log(pReverse/pForward);
1755 return logPropRatio;
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;
1776 static INT4 warningDelivered = 0;
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;
1793 return logPropRatio;
1810 const REAL8 epsTime = 6
e-6;
1811 const REAL8 epsAngle = 3
e-4;
1813 nRA = gsl_ran_ugaussian(rng);
1814 nDec = gsl_ran_ugaussian(rng);
1815 nTime = gsl_ran_ugaussian(rng);
1817 newRA += epsAngle*nRA;
1818 newDec += epsAngle*nDec;
1819 newTime += epsTime*nTime;
1827 nRefRA = (ra - refRA)/epsAngle;
1828 nRefDec = (dec - refDec)/epsAngle;
1829 nRefTime = (baryTime - refTime)/epsTime;
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);
1840 logPropRatio = log(pReverse/pForward);
1842 return logPropRatio;
1851 REAL8 logPropRatio = 0.0;
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);
1874 return logPropRatio;
1882 INT4 glitchLower, glitchUpper;
1887 REAL8 amparg, phiarg, Ai;
1890 gsl_matrix *glitch_f, *glitch_Q, *glitch_A;
1891 gsl_matrix *glitch_t, *glitch_p;
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);
1929 glitchLower = (
INT4)floor((f0 - 1./tau)/
deltaF);
1930 glitchUpper = (
INT4)floor((f0 + 1./tau)/
deltaF);
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);
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);
1946 Ai = Amp*tau*0.5*sqrt(
LAL_PI)*exp(-amparg*amparg)*asd->
data->
data[
i]/sqrt(Tobs);
1951 gRe -= Ai*cos(phiarg);
1952 gIm -= Ai*sin(phiarg);
1956 gRe += Ai*cos(phiarg);
1957 gIm += Ai*sin(phiarg);
1961 gRe = Ai*cos(phiarg);
1962 gIm = Ai*sin(phiarg);
1970 gsl_matrix_set(glitchFD, ifo, 2*
i, gRe);
1971 gsl_matrix_set(glitchFD, ifo, 2*
i+1, gIm);
1989 REAL8 sqrt3 = 1.7320508;
1998 if(SNR < 5.0) SNR = 5.0;
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;
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;
2020 REAL8 logPropRatio = 0.0;
2035 gsl_matrix *glitchFD, *glitch_f, *glitch_Q, *glitch_A, *glitch_t, *glitch_p;
2069 if (gsize->
data[ifo]==0) {
2075 return logPropRatio;
2079 n = (
INT4)floor(gsl_rng_uniform(rng) * (
REAL8)(gsize->
data[ifo]));
2085 t0 = gsl_matrix_get(glitch_t, ifo,
n);
2086 f0 = gsl_matrix_get(glitch_f, ifo,
n);
2087 Q = gsl_matrix_get(glitch_Q, ifo,
n);
2088 Amp = gsl_matrix_get(glitch_A, ifo,
n);
2089 phi0 = gsl_matrix_get(glitch_p, ifo,
n);
2094 params_x->
data[1] = f0;
2095 params_x->
data[2] =
Q;
2096 params_x->
data[3] = Amp * (0.25*Anorm);
2097 params_x->
data[4] = phi0;
2105 params_y->
data[
i] = params_x->
data[
i] + gsl_ran_ugaussian(rng)*sigmas_x->
data[
i]*scale;
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];
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);
2135 for(
i=0;
i<5;
i++) {
2149 qyx = eyx - log(nyx);
2150 qxy = exy - log(nxy);
2152 logPropRatio = qxy-qyx;
2160 return logPropRatio;
2171 REAL8 t=0, f=0,
Q=0, A=0;
2176 REAL8 pForward = 0.0;
2177 REAL8 pReverse = 0.0;
2178 REAL8 logPropRatio = 0.0;
2180 gsl_matrix *
params = NULL;
2201 draw = gsl_rng_uniform(rng);
2211 if(ny<nmin || ny>=nmax) {
2212 logPropRatio = -INFINITY;
2213 return logPropRatio;
2220 t =
draw_flat(thread,
"morlet_t0_prior");
2221 f =
draw_flat(thread,
"morlet_f0_prior");
2225 gsl_matrix_set(
params, ifo,
nx, t);
2229 gsl_matrix_set(
params, ifo,
nx, f);
2233 val =
draw_flat(thread,
"morlet_Q_prior");
2234 gsl_matrix_set(
params, ifo,
nx, val);
2243 gsl_matrix_set(
params, ifo,
nx, A);
2247 val =
draw_flat(thread,
"morlet_phi_prior");
2248 gsl_matrix_set(
params, ifo,
nx, val);
2267 draw = gsl_rng_uniform(rng);
2275 f = gsl_matrix_get(
params, ifo,
n);
2277 t = gsl_matrix_get(
params, ifo,
n);
2279 Q = gsl_matrix_get(
params, ifo,
n);
2281 A = gsl_matrix_get(
params, ifo,
n);
2286 gsl_matrix_set(
params, ifo,
i, gsl_matrix_get(
params, ifo,
i+1));
2288 gsl_matrix_set(
params, ifo,
i, gsl_matrix_get(
params, ifo,
i+1));
2290 gsl_matrix_set(
params, ifo,
i, gsl_matrix_get(
params, ifo,
i+1));
2292 gsl_matrix_set(
params, ifo,
i, gsl_matrix_get(
params, ifo,
i+1));
2294 gsl_matrix_set(
params, ifo,
i, gsl_matrix_get(
params, ifo,
i+1));
2320 pForward = qxy + qx;
2321 pReverse = qyx + qy;
2323 logPropRatio = pForward-pReverse;
2325 return logPropRatio;
2331 REAL8 logPropRatio = 0.0;
2341 phi = fmod(phi, 2.0*M_PI);
2342 psi = fmod(psi, M_PI);
2347 return logPropRatio;
2356 REAL8 logPropRatio = 0.0;
2372 draw = gsl_rng_uniform(rng);
2379 psi = (
alpha + beta)*0.5;
2380 phi = (
alpha - beta)*0.5;
2388 return logPropRatio;
2396 REAL8 logPropRatio = 0.0;
2403 plusminus = gsl_rng_uniform(thread->
GSLrandom);
2404 if ( plusminus < 0.5 )
2411 return logPropRatio;
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];
2427 REAL8 cosnewIota, cosnewIota2;
2429 INT4 nUniqueDet, det;
2445 cosIota = cos(iota);
2446 cosIota2 = cosIota*cosIota;
2450 for (det=0; det < nUniqueDet; det++) {
2465 y2[
i]=Fcross*Fcross;
2468 R2[
i] = (((1.0+cosIota2)*(1.0+cosIota2))/(4.0*dist2))*Fplus*Fplus
2469 + ((cosIota2)/(dist2))*Fcross*Fcross;
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] -
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]);
2481 (*newPsi) = (2.*atan((b -
a*sqrt((
a2 + b*b)/(
a2)))/
a))/4.;
2483 while ((*newPsi)<0){
2484 (*newPsi)=(*newPsi)+
LAL_PI/4.0;
2487 while ((*newPsi)>
LAL_PI/4.0){
2488 (*newPsi)=(*newPsi)-
LAL_PI/4.0;
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];
2495 newFcross[
i] =
y[
i]*cos(2.0*(*newPsi)) -
x[
i]*sin(2.0*(*newPsi));
2496 newFcross2[
i] = newFcross[
i] * newFcross[
i];
2499 c12 = -2.0*((R2[1]*(newFcross2[2])-R2[2]*(newFcross2[1]))
2500 /(R2[1]*(newFplus2[2])-R2[2]*(newFplus2[1])))-1.0;
2503 c12 = (3.0-c12)/(1.0+c12);
2504 (*newPsi) = (*newPsi)+
LAL_PI/4.0;
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];
2510 newFcross[
i] =
y[
i]*cos(2.0*(*newPsi)) -
x[
i]*sin(2.0*(*newPsi));
2511 newFcross2[
i] = newFcross[
i] * newFcross[
i];
2521 cosnewIota2 = c12-sqrt(c12*c12-1.0);
2522 cosnewIota = sqrt(cosnewIota2);
2523 *newIota = acos(cosnewIota);
2526 ((((1.0+cosnewIota2)*(1.0+cosnewIota2))/(4.0))*newFplus2[1]
2527 + (cosnewIota2)*newFcross2[1])
2530 if (Fplus*newFplus[3]<0){
2531 (*newPsi)=(*newPsi)+
LAL_PI/2.;
2532 newFcross[3]=-newFcross[3];
2535 if (Fcross*cosIota*cosnewIota*newFcross[3]<0){
2536 (*newIota)=
LAL_PI-(*newIota);
2568 REAL8 d_inner_h = MatchedFilterSNR * OptimalSNR;
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;
2578 new_x = xmean + gsl_ran_gaussian(thread->
GSLrandom,xsigma);
2579 REAL8 new_d = 1.0/new_x;
2583 OptimalSNR *= new_x / old_x;
2587 const char *ifonames[5]={
"H1",
"L1",
"V1",
"I1",
"J1"};
2591 sprintf(pname,
"%s_optimal_snr",ifonames[
i]);
2608 logxdjac = 2.0*log(new_d) - 2.0*log(old_d);
2619 REAL8 new_logd = log(new_d);
2621 logxdjac = log(new_d) - log(old_d);
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);
2628 return(log_p_reverse - log_p_forward + logxdjac);
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;
2647 REAL8 logPropRatio = 0.0;
2651 static INT4 warningDelivered = 0;
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;
2669 return logPropRatio;
2687 fprintf(stderr,
"LALInferenceExtrinsicParamProposal: No theta_jn parameter!\n");
2693 reflected_extrinsic_parameters(thread, ra, dec, baryTime,
dist, iota, psi, &newRA, &newDec, &newTime, &newDist, &newIota, &newPsi);
2696 const REAL8 epsDist = 1
e-8;
2697 const REAL8 epsTime = 1
e-8;
2698 const REAL8 epsAngle = 1
e-8;
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);
2707 newRA += epsAngle*nRA;
2708 newDec += epsAngle*nDec;
2709 newTime += epsTime*nTime;
2710 newDist += epsDist*nDist;
2711 newIota += epsAngle*nIota;
2712 newPsi += epsAngle*nPsi;
2716 reflected_extrinsic_parameters(thread, newRA, newDec, newTime, newDist, newIota, newPsi, &refRA, &refDec, &refTime, &refDist, &refIota, &refPsi);
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;
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);
2736 REAL8 logNewDist = log(newDist);
2739 REAL8 newcosIota = cos(newIota);
2743 logPropRatio = pReverse - pForward;
2745 return logPropRatio;
2755 REAL8FFTPlan **plans;
2765 plans =
XLALCalloc(nDet,
sizeof(REAL8FFTPlan *));
2767 for (
i=0;
i<nDet;
i++) {
2771 asds[
i] =
data->noiseASD;
2772 psds[
i] =
data->oneSidedNoisePowerSpectrum;
2774 td_data[
i] =
data->timeData;
2775 fd_data[
i] =
data->freqData;
2777 plans[
i] =
data->freqToTimeFFTPlan;
2794 INT4 no_adapt, adapting;
2795 INT4 adaptTau, adaptableStep, adaptLength, adaptResetBuffer;
2796 REAL8 sigma, s_gamma;
2797 REAL8 logPAtAdaptStart = -INFINITY;
2801 for(
this=
params->head;
this;
this=this->next) {
2803 char *
name = this->name;
2805 if (!strcmp(
name,
"eta") || !strcmp(
name,
"q") || !strcmp(
name,
"time") || !strcmp(
name,
"a_spin2") || !strcmp(
name,
"a_spin1") || !strcmp(
name,
"t0")){
2807 }
else if (!strcmp(
name,
"polarisation") || !strcmp(
name,
"phase") || !strcmp(
name,
"costheta_jn")){
2828 adapting = !no_adapt;
2836 adaptLength = pow(10, adaptTau);
2837 adaptResetBuffer = 100;
2887 INT4 adaptableStep = 0;
2889 REAL8 priorMin, priorMax, dprior, s_gamma;
2890 REAL8 accept, propose, sigma;
2903 if (adaptableStep && adapting) {
2922 if (adaptableStep) {
2935 dprior = priorMax - priorMin;
2949 sigma += s_gamma * (dprior/100.0) * (1.0-targetAcceptance);
2951 sigma -= s_gamma * (dprior/100.0) * (targetAcceptance);
2954 sigma = (sigma > dprior ? dprior : sigma);
2955 sigma = (sigma < DBL_MIN ? DBL_MIN : sigma);
3002 for (
j=0;
j<nCols;
j++)
3006 for (
j=0;
j<nCols;
j++) {
3007 if (!strcmp(
"logl",
params[
j])) {
3012 char* internal_param_name =
XLALCalloc(512,
sizeof(
char));
3016 if (!strcmp(item->
name, internal_param_name) &&
3029 for (item = backwardClusterParams->
head; item; item = item->
next)
3044 if (acl_real8 == INFINITY)
3046 else if (acl_real8 < 1.0)
3049 acl = (
INT4)acl_real8;
3051 INT4 downsampled_size = ceil((
REAL8)nInSamps/acl);
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];
3059 sampleArray = downsampled_array;
3060 nInSamps = downsampled_size;
3069 fprintf(stderr,
"\nERROR: Couldn't build kmeans clustering from the file specified.\n");
3106 INT4 cyclic_reflective,
3110 gsl_matrix_view mview;
3111 char outp_name[256];
3112 char outp_draws_name[256];
3119 mview = gsl_matrix_view_array(array, nSamps, dim);
3120 kde->
kmeans = (*cluster_method)(&mview.matrix, ntrials, thread->
GSLrandom);
3138 printf(
"Thread %i found %i clusters.\n", thread->
id, kde->
kmeans->
k);
3140 sprintf(outp_name,
"clustered_samples.%2.2d", thread->
id);
3141 sprintf(outp_draws_name,
"clustered_draws.%2.2d", thread->
id);
3163 outp = fopen(outp_name,
"w");
3191 outp = fopen(outp_name,
"w");
3195 for (
i=0;
i<nSamps;
i++) {
3226 if (!strcmp(existing_kde->
name, kde->
name)) {
3227 old_kde = existing_kde;
3231 while (existing_kde->
next != NULL) {
3233 if (!strcmp(existing_kde->
next->name, kde->
name)) {
3234 old_kde = existing_kde->
next;
3236 existing_kde->
next = kde;
3239 existing_kde = existing_kde->
next;
3243 existing_kde->
next=kde;
3262 if (proposal != NULL) {
3288 if (effSampleSize > 0)
3289 step = (
INT4) floor(bufferSize/effSampleSize);
3299 for (
i=0;
i < nPoints;
i++)
3300 DEsamples[
i] =
temp + (
i*nPar);
3337 for (item = backwardClusterParams->
head; item; item = item->
next)
3375 return logPropRatio;
3393 REAL8 cumulativeWeight, totalWeight;
3394 REAL8 logPropRatio = 0.0;
3401 return logPropRatio;
3412 totalWeight += kde->
weight;
3420 cumulativeWeight = kde->
weight;
3421 while(cumulativeWeight/totalWeight < randomDraw) {
3423 cumulativeWeight += kde->
weight;
3441 if (propDensity == NULL || *propDensity == -INFINITY)
3444 logCurrentP = *propDensity;
3448 logPropRatio = logCurrentP - logProposedP;
3450 if (propDensity != NULL)
3451 *propDensity = logProposedP;
3456 return logPropRatio;
3475 INT4 nPar, nPoints, nSkip;
3490 for (
i=0;
i < nPoints;
i++)
3491 DEarray[
i] =
temp + (
i*nPar);
3496 if (max_acl == INFINITY)
3498 else if (max_acl < 1.0)
3501 *maxACL = (
INT4)max_acl;
3534 INT4 par=0, lag=0,
i=0, imax;
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;
3550 while (cumACF >=
s) {
3551 ACF = gsl_stats_correlation(array + par, nPar, array + lag*nPar + par, nPar, nPoints-lag);
3552 cumACF += 2.0 * ACF;
3563 else if (gsl_isnan(ACL))
3566 for (
i=0;
i<nPoints;
i++)
3567 array[
i*nPar + par] +=
mean;
3607 INT4 iEff = nPoints/acl;
3625 fprintf(
fp,
"%9.5f\t", exp(logPropRatio));
3636 for(det=thread->
parent->data;det;det=det->
next)
3647 for (
i = 0;
i < nspl;
i++) {
3650 fprintf(stderr,
"variable name too long\n"); exit(1);
3654 fprintf(stderr,
"variable name too long\n"); exit(1);
3659 amp += ampWidth*gsl_ran_ugaussian(thread->
GSLrandom)/sqrt(nifo*(
REAL8)nspl);
3664 ph += phaseWidth*gsl_ran_ugaussian(thread->
GSLrandom)/sqrt(nifo*(
REAL8)nspl);
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.
REAL8 * LALInferenceKmeansDraw(LALInferenceKmeans *kmeans)
Draw a sample from a kmeans-KDE estimate of a distribution.
void LALInferenceKmeansDestroy(LALInferenceKmeans *kmeans)
Destroy a kmeans instance.
LALInferenceKmeans * LALInferenceOptimizedKmeans(gsl_matrix *data, INT4 ntrials, gsl_rng *rng)
Kmeans cluster data, maximizing BIC assuming BIC(k) is concave-down.
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
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)
LALInferenceVariables currentParams
static double double delta
void XLALComputeDetAMResponse(double *fplus, double *fcross, const REAL4 D[3][3], const double ra, const double dec, const double psi, const double gmst)
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.
void LALInferenceReadAsciiHeader(FILE *input, char params[][VARNAME_MAX], INT4 *nCols)
Read column names from an ASCII file.
REAL8 * LALInferenceParseDelimitedAscii(FILE *input, INT4 nCols, INT4 *wantedCols, INT4 *nLines)
Utility for readling in delimited ASCII files.
LALInferenceVariableType LALInferenceGetVariableType(const LALInferenceVariables *vars, const char *name)
Get the LALInferenceVariableType of the parameter named name in vars.
void LALInferenceDiscardPTMCMCHeader(FILE *filestream)
Discard the standard header of a PTMCMC chain file.
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.
gsl_matrix * LALInferenceGetgslMatrixVariable(LALInferenceVariables *vars, const char *name)
REAL8Vector * LALInferenceGetREAL8VectorVariable(LALInferenceVariables *vars, const char *name)
void LALInferenceAddREAL8VectorVariable(LALInferenceVariables *vars, const char *name, REAL8Vector *value, LALInferenceParamVaryType vary)
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.
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)
void LALInferenceCopyVariables(LALInferenceVariables *origin, LALInferenceVariables *target)
Deep copy the variables from one to another LALInferenceVariables structure.
void LALInferenceBurninStream(FILE *filestream, INT4 burnin)
Burn-in a generic ASCII stream.
void LALInferenceAddREAL8Variable(LALInferenceVariables *vars, const char *name, REAL8 value, LALInferenceParamVaryType vary)
INT4 LALInferenceGetINT4Variable(LALInferenceVariables *vars, const char *name)
void LALInferenceAddINT4Variable(LALInferenceVariables *vars, const char *name, INT4 value, LALInferenceParamVaryType vary)
UINT4Vector * LALInferenceGetUINT4VectorVariable(LALInferenceVariables *vars, const char *name)
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.
REAL8(* LALInferenceProposalFunction)(struct tagLALInferenceThreadState *thread, LALInferenceVariables *currentParams, LALInferenceVariables *proposedParams)
Jump proposal distribution Computes proposedParams based on currentParams and additional variables st...
void LALInferenceClearVariables(LALInferenceVariables *vars)
Delete the variables in this structure.
ProcessParamsTable * LALInferenceGetProcParamVal(ProcessParamsTable *procparams, const char *name)
Returns the element of the process params table with "name".
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...
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...
const CHAR * LALInferenceGetstringVariable(LALInferenceVariables *vars, const char *name)
int LALInferenceCheckVariable(LALInferenceVariables *vars, const char *name)
Checks for name being present in vars returns 1(==true) or 0.
LALInferenceVariableItem * LALInferenceGetItemNr(LALInferenceVariables *vars, int idx)
Return the list node for the idx-th item - do not rely on this Indexing starts at 1.
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.
int LALInferenceCheckVariableNonFixed(LALInferenceVariables *vars, const char *name)
Checks for name being present in vars and having type LINEAR or CIRCULAR.
@ LALINFERENCE_PARAM_OUTPUT
A parameter that never changes, functions should respect this.
@ LALINFERENCE_PARAM_FIXED
A parameter that is cyclic, such as an angle between 0 and 2pi.
@ LALINFERENCE_PARAM_LINEAR
@ LALINFERENCE_void_ptr_t
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.
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.
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.
LALInferenceProposalCycle * LALInferenceInitProposalCycle(void)
Create a new proposal cycle.
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.
LALInferenceVariables * LALInferenceParseProposalArgs(LALInferenceRunState *runState)
Go through all logic for deciding proposals to use.
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.
LALInferenceProposalCycle * LALInferenceSetupDefaultInspiralProposalCycle(LALInferenceVariables *propArgs)
A reasonable default 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.
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.
LALInferenceProposal * LALInferenceInitProposal(LALInferenceProposalFunction func, const char *name)
Creates a new proposal object from the given func pointer and name.
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)
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
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
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.
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)
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)
REAL8 XLALGreenwichMeanSiderealTime(const LIGOTimeGPS *gpstime)
LIGOTimeGPS * XLALGPSSetREAL8(LIGOTimeGPS *epoch, REAL8 t)
REAL8 XLALGPSGetREAL8(const LIGOTimeGPS *epoch)
LALDetector detectors[LAL_NUM_DETECTORS]
Struct to hold clustered-KDE estimates.
LALInferenceKmeans * kmeans
struct tagLALInferenceClusteredKDEProposal * next
LALInferenceVariables * params
Structure to contain IFO data.
LALDetector * detector
integration limits for overlap integral in F-domain
struct tagLALInferenceIFOData * next
ROQ data.
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.
LALInferenceVariables * params
Structure for holding a proposal cycle.
LALInferenceProposal ** proposals
char last_proposal_name[VARNAME_MAX]
Counter for cycling through proposals.
INT4 * order
Array of proposals (one per proposal function)
INT4 nProposals
Length of cycle.
Structure for holding a LALInference proposal, along with name and stats.
LALInferenceProposalFunction func
Structure containing inference run state This includes pointers to the function types required to run...
ProcessParamsTable * commandLine
LALInferenceVariables * proposalArgs
The data from the interferometers.
LALInferenceVariables * algorithmParams
Any special arguments for the prior function.
struct tagLALInferenceIFOData * data
Log sample, i.e.
LALInferenceAlgorithm algorithm
A function that returns a new set of variables for the model.
Structure containing chain-specific variables.
LALInferenceVariables * currentParams
Prior boundaries, etc.
struct tagLALInferenceRunState * parent
LALInferenceVariables * priorArgs
Stope things such as output arrays.
LALInferenceModel * model
Cycle of proposals to call.
INT4 accepted
This should be removed, can be given as an algorithmParams entry.
size_t differentialPointsLength
Array of points for differential evolution.
LALInferenceProposalCycle * cycle
The proposal cycle.
LALInferenceVariables * proposalArgs
REAL8 temperature
Array containing multiple proposal densities.
size_t differentialPointsSkip
Size of the differentialPoints memory block (must be >= length of differential points).
LALInferenceVariables ** differentialPoints
Parameters proposed.
The LALInferenceVariableItem list node structure This should only be accessed using the accessor func...
LALInferenceVariableType type
struct tagVariableItem * next
LALInferenceParamVaryType vary
The LALInferenceVariables structure to contain a set of parameters Implemented as a linked list of LA...
LALInferenceVariableItem * head
CHAR value[LIGOMETA_VALUE_MAX]
struct tagProcessParamsTable * next