Commit bd2b0518 authored by Marcella Wijngaarden's avatar Marcella Wijngaarden
Browse files

Merge branch 'extrinsic-burnin' into BW+QuickCBC

parents a878d5bd 86c396c3
......@@ -1378,152 +1378,14 @@ void CBC_start(struct Net *net, int *mxc, FILE *chainS, double **paramx, double
// Initial search to find CBC intrinsic model signal
MCMC_intrinsic(net, 1, mxc, 10000, chainS, paramx, skyx, pallx, who, heat, history, global, freq, D, SN, N, Tobs, r, rundata);
// Initial (optional) search to find extrinsic model parameters
if (rundata->intrinsic_only == 0)
{
int mty, mtid, mti;
double *logL, *logLsky;
double **data, ***wave;
// double *Larray, **Tarray;
double *DD;
double **WW;
double ***DHc, ***DHs, ***HH;
logLsky = double_vector(NC);
// whitened signal in each detector
wave = double_tensor(NC,net->Nifo,N);
// whitened data in each detector
data = double_matrix(net->Nifo,N);
for (id = 0; id < net->Nifo; ++id)
{
data[id][0] = 0.0;
data[id][N/2] = 0.0;
}
for (i = 1; i < N/2; ++i)
{
for (id = 0; id < net->Nifo; ++id)
{
x = 1.0/sqrt(SN[id][i]);
data[id][i] = D[id][2*i]*x;
data[id][N-i] = D[id][2*i+1]*x;
}
}
DD = double_vector(net->Nifo);
WW = double_matrix(NC,net->Nifo);
DHc = double_tensor(NC,net->Nifo,nt+1);
DHs = double_tensor(NC,net->Nifo,nt+1);
HH = double_tensor(NC,net->Nifo,3);
Lmax = 0.0;
printf("\nFinding sky location\n");
// Find sky location and set up whitend signal arrays
#pragma omp parallel for
for(j = 0; j < NC; j++)
{
int mtid, mti;
double mty, Scale;
double *SNRsq;
double **hwave;
SNRsq = double_vector(net->Nifo);
hwave = double_matrix(net->Nifo,N);
// hwave in [0,N/2] real and [N/2,N] imag
templates(net, hwave, freq, paramx[j], N, rundata);
// might want to put a catch here for any chains that didn't reach a decent logL and re-run the rearch phase
// logLstart[j] = log_likelihood(net, D, paramx[j], freq, SN, N, Tobs, rundata);
mty = 0.0;
for (mtid = 0; mtid < net->Nifo; ++mtid)
{
SNRsq[mtid] = 0.0;
for (mti = 1; mti < N/2; ++mti)
{
// The sky mapping code puts in the amplitude and phase shifts. The signals only differ due to the whitening
SNRsq[mtid] += 4.0*(hwave[0][mti]*hwave[0][mti]+hwave[0][N-mti]*hwave[0][N-mti])/SN[mtid][mti];
}
Scale = 1.0;
if(mtid > 0) Scale = paramx[j][(mtid-1)*3+NX+2]*paramx[j][(mtid-1)*3+NX+2];
SNRsq[mtid] *= Scale;
mty += SNRsq[mtid];
}
printf("%f %f\n", SNRsq[0], SNRsq[1]);
// find a sky location roughly consistent with the time delays
skyring(net, paramx[j], skyx[j], pallx[j], SNRsq, freq, D, SN, N, Tobs, rundata->rvec[j], NX, NS, rundata->gmst);
// // find a sky location roughly consistent with the time delay, amplitude ratio and phase difference
// // skystart(net, paramx[j], skyx[j], pallx[j], SNRsq, freq, D, SN, N, Tobs, rvec[j]);
dshifts(net, skyx[j], paramx[j], NX, rundata);
// map to the geocenter params
pmap_all(net, pallx[j], paramx[j], skyx[j], NX, rundata->gmst);
// get the geocenter reference template. This does depend on the assumed sky location, hence follows skystart
geotemplate(hwave[0], freq, pallx[j], N, rundata);
for (mtid = 0; mtid < net->Nifo; ++mtid)
{
wave[j][mtid][0] = 0.0;
wave[j][mtid][N/2] = 0.0;
for (mti = 1; mti < N/2; ++mti)
{
// The sky mapping code puts in the amplitude and phase shifts. The signals only differ due to the whitening
mty = 1.0/sqrt(SN[mtid][mti]);
wave[j][mtid][mti] = hwave[0][mti]*mty;
wave[j][mtid][N-mti] = hwave[0][N-mti]*mty;
}
}
skylikesetup(net, data, wave[j], DD, WW[j], DHc[j], DHs[j], Tobs, N, bn, nt, imin, imax);
fisherskysetup(net, wave[j], HH[j], Tobs, N);
logLsky[j] = skylike(net, skyx[j], DD, WW[j], DHc[j], DHs[j], dtx, nt, 0, rundata);
logLfull[j] = log_likelihood_full(net, D, pallx[j], freq, SN, rho[j], N, Tobs, rundata);
printf("%d logL initial %f full %f sky %f\n", j, logLstart[j], logLfull[j], logLsky[j]);
free_double_vector(SNRsq);
free_double_matrix(hwave,net->Nifo);
}
chain = fopen("sky.dat","w");
skymcmc(net, 1000000, mxc, chain, paramx, skyx, pallx, rundata->who, heat, dtx, nt, DD, WW, DHc, DHs, HH, Tobs, r, rundata);
fclose(chain);
free_double_vector(logLsky);
free_double_tensor(wave,NC,net->Nifo);
free_double_vector(DD);
free_double_matrix(WW,NC);
free_double_tensor(DHc,NC,net->Nifo);
free_double_tensor(DHs,NC,net->Nifo);
free_double_tensor(HH,NC,net->Nifo);
free_double_matrix(data,net->Nifo);
}
for(j = 0; j < NC; j++)
{
// If we are only doing intrinsic burnin, use sky parameters given by BW
if (rundata->intrinsic_only == 1)
for(k = 0; k < 4; k++)
{
for(k = 0; k < 4; k++)
{
skyx[j][k] = pallx[j][NX+k];
}
skyx[j][k] = pallx[j][NX+k];
}
pmap_all(net, pallx[j], paramx[j], skyx[j], NX, rundata->gmst);
logLfull[j] = log_likelihood_full(net, D, pallx[j], freq, SN, rho[j], N, Tobs, rundata);
......@@ -4556,7 +4418,7 @@ void AmpPhase(double *Af, double *Pf, RealVector *freq, double *params, int N, s
distance = exp(params[6]);
/* ==== TIDAL PARAMETERS ==== */
fRef_in = fref;
fRef_in = rundata->fref;
NRTidal_version = rundata->NRTidal_version;
lambda_type_version = rundata->lambda_type_version;
......@@ -4646,7 +4508,7 @@ void PDwave(double *wavef, RealVector *freq, double *params, int N, struct bayes
distance = exp(params[6]);
/* ==== TIDAL PARAMETERS ==== */
fRef_in = fref;
fRef_in = rundata->fref;
NRTidal_version = rundata->NRTidal_version;
lambda_type_version = rundata->lambda_type_version;
......@@ -4696,130 +4558,6 @@ void PDwave(double *wavef, RealVector *freq, double *params, int N, struct bayes
}
// Note, returns hwave in array structure real[<N/2], imag[>N/2]
// Only used in burnin (extrinsic updates)
void templates(struct Net *net, double **hwave, RealVector *freq, double *params, int N, struct bayesCBC *rundata)
{
AmpPhaseFDWaveform *ap = NULL;
NRTidal_version_type NRTidal_version;
Lambda_type lambda_type_version;
double phi0, fRef_in, mc, q, m1_SI, m2_SI, chi1, chi2, f_min, f_max, distance;
int ret, flag, i, j;
double p, cp, sp;
double f, x, y, deltaF, ts, Tobs, sqT;
double pd, Ar, td, A, fs;
double *extraParams;
double mt, eta, dm;
double lambda1, lambda2;
double lambdaT, dLambdaT, sym_mass_ratio_eta;
int NX = rundata->NX;
int NP = rundata->NP;
int NS = rundata->NS;
int NC = rundata->NC; // number of chains
int NCC = rundata->NCC; // number of cold chains
int NH = rundata->NH; // length of history
int NQ = rundata->NQ; // number of mass ratios in global proposal
int NM = rundata->NM; // number of chirp masses in global proposal
extraParams = (double*)malloc(sizeof(double)* 2);
Tobs = 1.0/freq->data[1];
sqT = 1.0;//sqrt(Tobs);
fs = fbegin(params);
mc = exp(params[0]);
mt = exp(params[1]);
eta = pow((mc/mt), (5.0/3.0));
if(eta > 0.25)
{
dm = 0.0;
}
else
{
dm = sqrt(1.0-4.0*eta);
}
m1_SI = mt*(1.0+dm)/2.0;
m2_SI = mt*(1.0-dm)/2.0;
chi1 = params[2];
chi2 = params[3];
phi0 = params[4];
ts = Tobs-params[5];
distance = exp(params[6]);
/* ==== TIDAL PARAMETERS ==== */
fRef_in = fref;
NRTidal_version = rundata->NRTidal_version;
lambda_type_version = rundata->lambda_type_version;
if (NRTidal_version != NoNRT_V) {
// Compute lambda1 & lambda2
if (lambda_type_version == lambdaTilde) {
lambdaT = params[7];
dLambdaT = params[8];
sym_mass_ratio_eta = m1_SI*m2_SI/((m1_SI+m2_SI)*(m1_SI+m2_SI));
LambdaTsEta2Lambdas(lambdaT, dLambdaT, sym_mass_ratio_eta, &lambda1, &lambda2);
extraParams[0] = lambda1; // lambda1
extraParams[1] = lambda2; // lambda2
} else {
extraParams[0] = params[7]; // lambda1
extraParams[1] = params[8]; // lambda2
}
}
ret = IMRPhenomDGenerateh22FDAmpPhase(&ap, freq, phi0, fRef_in, m1_SI, m2_SI, chi1, chi2, distance, extraParams, NRTidal_version);
for (j=0; j< net->Nifo; j++)
{
pd = 0.0; // phase offset
td = 0.0; // time offset
Ar = 1.0; // amplitude ratio
if(j > 0)
{
pd = params[(j-1)*3+NX]; // phase offset
td = params[(j-1)*3+NX+1]; // time offset
Ar = params[(j-1)*3+NX+2]; // amplitude
}
hwave[j][0] = 0.0;
hwave[j][N/2] = 0.0;
for (i=1; i< N/2; i++)
{
f = freq->data[i];
if(f > fs)
{
A = Ar*h22fac*ap->amp[i]/sqT;
p = ap->phase[i];
x = TPI*f*(ts-td)+pd-p;
hwave[j][i] = A*cos(x);
hwave[j][N-i] = A*sin(x);
}
else
{
hwave[j][i] = 0.0;
hwave[j][N-i] = 0.0;
}
}
}
DestroyAmpPhaseFDWaveform(ap);
free(extraParams);
}
// TODO: Rewrite!
void cbcGeoTemplate(struct bayesCBC *rundata, int ic, int N, double **template)
{
......@@ -4940,7 +4678,7 @@ void geowave(double *gwave, RealVector *freq, double *params, int N, struct baye
distance = exp(params[6]);
/* ==== TIDAL PARAMETERS ==== */
fRef_in = fref;
fRef_in = rundata->fref;
NRTidal_version = rundata->NRTidal_version;
......@@ -4993,108 +4731,6 @@ void geowave(double *gwave, RealVector *freq, double *params, int N, struct baye
}
// Returns hwave in array structure real[<N/2], imag[>N/2]
// Only used in burnin (extrinsic updates)
void geotemplate(double *gwave, RealVector *freq, double *params, int N, struct bayesCBC *rundata)
{
AmpPhaseFDWaveform *ap = NULL;
NRTidal_version_type NRTidal_version;
double phi0, fRef_in, mc, q, m1_SI, m2_SI, chi1, chi2, f_min, f_max, distance;
int ret, flag, i, j, id;
double p, cp, sp;
double f, x, y, deltaF, ts, Tobs, sqT;
double pd, Ar, td, A;
double *extraParams;
double mt, eta, dm, fs;
double lambda1, lambda2;
extraParams = (double*)malloc(sizeof(double)* 2);
fs = fbegin(params);
Tobs = 1.0/freq->data[1];
sqT = 1.0;//sqrt(Tobs);
mc = exp(params[0]);
mt = exp(params[1]);
eta = pow((mc/mt), (5.0/3.0));
if(eta > 0.25)
{
dm = 0.0;
}
else
{
dm = sqrt(1.0-4.0*eta);
}
m1_SI = mt*(1.0+dm)/2.0;
m2_SI = mt*(1.0-dm)/2.0;
chi1 = params[2];
chi2 = params[3];
phi0 = 0.5*params[4]; // I'm holding the GW phase in [4], while PhenomD wants orbital
ts = Tobs-params[5];
distance = exp(params[6]);
/* ==== TIDAL PARAMETERS ==== */
fRef_in = fref;
NRTidal_version = rundata->NRTidal_version;
double lambdaT, dLambdaT, sym_mass_ratio_eta;
Lambda_type lambda_type;
lambda_type = rundata->lambda_type_version;
if (NRTidal_version != NoNRT_V) {
// Compute lambda1 & lambda2
if (lambda_type == lambdaTilde) {
lambdaT = params[7];
dLambdaT = params[8];
sym_mass_ratio_eta = m1_SI*m2_SI/((m1_SI+m2_SI)*(m1_SI+m2_SI));
LambdaTsEta2Lambdas(lambdaT, dLambdaT, sym_mass_ratio_eta, &lambda1, &lambda2);;
extraParams[0] = lambda1; // lambda1
extraParams[1] = lambda2; // lambda2
} else {
extraParams[0] = params[7]; // lambda1
extraParams[1] = params[8]; // lambda2
}
}
ret = IMRPhenomDGenerateh22FDAmpPhase(&ap, freq, phi0, fRef_in, m1_SI, m2_SI, chi1, chi2, distance, extraParams, NRTidal_version);
gwave[0] = 0.0;
gwave[N/2] = 0.0;
for (i=1; i< N/2; i++)
{
f = freq->data[i];
if(f > fs)
{
A = h22fac*ap->amp[i]/sqT;
p = ap->phase[i];
x = TPI*f*ts-p;
gwave[i] = A*cos(x);
gwave[N-i] = A*sin(x);
}
else
{
gwave[i] = 0.0;
gwave[N-i] = 0.0;
}
}
DestroyAmpPhaseFDWaveform(ap);
free(extraParams);
}
void detector_shifts(struct Net *net, double *params, int NX, double gmst)
{
......@@ -5213,7 +4849,7 @@ void fulltemplates(struct Net *net, double **hwave, RealVector *freq, double *pa
distance = exp(params[6]);
/* ==== TIDAL PARAMETERS ==== */
fRef_in = fref;
fRef_in = rundata->fref;
NRTidal_version = rundata->NRTidal_version;
double lambdaT, dLambdaT, sym_mass_ratio_eta;
......@@ -5344,7 +4980,7 @@ void fullphaseamp(struct Net *net, double **amp, double **phase, RealVector *fre
distance = exp(params[6]);
/* ==== TIDAL PARAMETERS ==== */
fRef_in = fref;
fRef_in = rundata->fref;
NRTidal_version = rundata->NRTidal_version;
double lambdaT, dLambdaT, sym_mass_ratio_eta;
......@@ -5603,7 +5239,7 @@ void PDtemplates(double *waveH, double *waveL, RealVector *freq, double *params,
Ar = params[NX+2]; // amplitude ratio L/H
/* ==== TIDAL PARAMETERS ==== */
fRef_in = fref;
fRef_in = rundata->fref;
NRTidal_version = rundata->NRTidal_version;
double lambda1, lambda2, lambdaT, dLambdaT, sym_mass_ratio_eta;
......@@ -5694,59 +5330,7 @@ double Fmag(double *sky, int id, double gmst)
}
void skyring(struct Net *net, double *params, double *sky, double *pall, double *SNRsq, RealVector *freq, double **D, double **SN, int N, double Tobs, gsl_rng * r, int NX, int NS, double gmst)
{
int i, j, k, id;
int iref;
double *delt;
double Fp, Fc, Fs;
double logL, logLmax;
double *stry;
iref = net->labels[0];
stry = double_vector(NS);
delt = double_vector(net->Nifo);
for(id = 1; id< net->Nifo; id++)
{
delt[id] = -params[(id-1)*3+NX+1];
}
// sky parameter order
//[0] alpha, [1] sin(delta) [2] psi [3] cosi [4] scale [5] dt [6] dphi
// Scale changes the overall distance reference at geocenter
// dt shifts the geocenter time
// dphi shifts the geocenter phase
stry[4] = 1.0;
stry[5] = 0.0;
stry[6] = 0.0;
// find sky locations consistent with time delays
ringfind(net, delt, stry, SNRsq, r, gmst);
stry[2] = PI*gsl_rng_uniform(r);
stry[3] = -1.0+2.0*gsl_rng_uniform(r);
pmap_all(net, pall, params, stry, NX, gmst);
for(j = 0; j< NS; j++) sky[j] = stry[j];
free_double_vector(delt);
free_double_vector(stry);
}
void skystart(struct Net *net, double *params, double *sky, double *pall, double *SNRsq, RealVector *freq, double **D, double **SN, int N, double Tobs, gsl_rng * r, struct bayesCBC *rundata)
void skystart(struct Net *net, double *params, double *sky, double *pall, double *SNRsq, RealVector *freq, double **D, double **SN, int N, double Tobs, gsl_rng * r, struct bayesCBC *rundata)
{
int i, j, k, id;
int iref;
......@@ -5816,6 +5400,7 @@ void skystart(struct Net *net, double *params, double *sky, double *pall, double
}
}
}
free_double_vector(rho);
......@@ -5825,79 +5410,6 @@ void skystart(struct Net *net, double *params, double *sky, double *pall, double
}
void dshifts(struct Net *net, double *sky, double *params, int NX, struct bayesCBC *rundata)
{
int i, j, k, id;
double dtimes[4];
double Fplus[4], Fcross[4], F[4];
double lambda[4];
double logL, logLmax, ecc;
// sky parameter order
//[0] alpha, [1] sin(delta) [2] psi [3] cos(iota) [4] scale [5] dt [6] phi0
Times(sky[0], sky[1], dtimes, rundata->gmst);
for(id = 0; id<net->Nifo; id++)
{
ecc = 2.0*sky[3]/(1.0+sky[3]*sky[3]);
F_ant(sky[2], sky[0], sky[1], &Fplus[id], &Fcross[id], net->labels[id], rundata->gmst);
F[id] = sqrt(Fplus[id]*Fplus[id]+ecc*ecc*Fcross[id]*Fcross[id]);
lambda[id] = atan2(ecc*Fcross[id],Fplus[id]);
if(lambda[id] < 0.0) lambda[id] += TPI;
}
// Note that F[] and lambda[] have already used the label array to get the correct detector ordering while dtimes uses a fixed detector ordering so needs the label array to get the correct delays
for(id = 1; id<net->Nifo; id++)
{
params[(id-1)*3+NX] = lambda[id]-lambda[0]; // GW phase offset
if(params[(id-1)*3+NX] < 0.0) params[(id-1)*3+NX] += TPI;
if(params[(id-1)*3+NX] > TPI) params[(id-1)*3+NX] -= TPI;
params[(id-1)*3+NX+1] = dtimes[net->labels[id]]-dtimes[net->labels[0]]; // time offset
params[(id-1)*3+NX+2] = F[id]/F[0]; // amplitude ratio
}
}
void fisherskysetup(struct Net *net, double **wave, double **HH, double Tobs, int n)
{
double f;
int i, j, k, l, bn, ii, id;
double **wavef;
wavef = double_matrix(net->Nifo,n);
for(id = 0; id<net->Nifo; id++)
{
wavef[id][0] = 0.0;
wavef[id][n/2] = 0.0;
for(i = 1; i<n/2; i++)
{
f = (double)(i)/Tobs;
j = i;
k = n-i;
wavef[id][j] = f*wave[id][j];
wavef[id][k] = f*wave[id][k];
}
}
for(id = 0; id<net->Nifo; id++)
{
HH[id][0] = 4.0*f_nwip(wave[id], wave[id], n);
HH[id][1] = 4.0*f_nwip(wave[id], wavef[id], n);
HH[id][2] = 4.0*f_nwip(wavef[id], wavef[id], n);
}
free_double_matrix(wavef,net->Nifo);
}
void upsample(int n, double Tobs, int *nt, int *bn)
{
double dt, dtres;
......@@ -5913,801 +5425,20 @@ void upsample(int n, double Tobs, int *nt, int *bn)
// k = 1;
// upsampled
*bn = n*k;
dt = Tobs/(double)(n*k);
// allow for geocenter - surface light travel time
// and for attempts to slide the waveform by the max delay time
*nt = 2*(int)(dtmax/dt)+1;
}
void skylikesetup(struct Net *net, double **data, double **wave, double *D, double *H, double **DHc, double **DHs, double Tobs, int n, int bn, int nt, int imin, int imax)
{
double *corr, *corrf;
double dt;
int i, j, k, kk, l, ii, id;
dt = Tobs/(double)(bn);
corr = double_vector(bn);
corrf = double_vector(bn);
for(id = 0; id<net->Nifo; id++)
{
D[id] = 4.0*f_nwip(data[id], data[id], n);
H[id] = 4.0*f_nwip(wave[id], wave[id], n);
// The arrays are zero padded to increase the sample cadence
for (i=0; i < bn; i++)
{
corr[i] = 0.0;
corrf[i] = 0.0;
}
// the non-zero part
for (i=1; i < n/2; i++)
{
l=i;
k=bn-i;
kk = n-i;
corr[l] = ( data[id][l]*wave[id][l] + data[id][kk]*wave[id][kk]);
corr[k] = ( data[id][kk]*wave[id][l] - data[id][l]*wave[id][kk]);
corrf[l] = corr[k];
corrf[k] = -corr[l];
}
gsl_fft_halfcomplex_radix2_inverse(corr, 1, bn);
gsl_fft_halfcomplex_radix2_inverse(corrf, 1, bn);
k = (nt-1)/2;
for(i=-k; i<k; i++)