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

bugfixes and convert timedelays between BW/bCBC

parent c98c6560
Pipeline #163698 passed with stages
in 1 minute and 47 seconds
......@@ -829,7 +829,7 @@ void printwave(struct Net *net, int N, RealVector *freq, double **paramx, int *w
void pmap(struct Net *net, double *pall, double *param, double *sky, int NX, double gmst)
{
int j;
double ciota, Ap, Ac, alpha, sindelta, psi, Fs;
double ciota, Ap, Ac, alpha, sindelta, psi, Fs, ecc;
double *dtimes;
double Fplus, Fcross, lambda;
......@@ -847,6 +847,7 @@ void pmap(struct Net *net, double *pall, double *param, double *sky, int NX, dou
ciota = sky[3];
Ap = (1.0+ciota*ciota)/2.0;
Ac = ciota;
ecc = Ac/Ap;
alpha = sky[0];
sindelta = sky[1];
psi = sky[2];
......@@ -871,56 +872,6 @@ void pmap(struct Net *net, double *pall, double *param, double *sky, int NX, dou
free(dtimes);
}
/**
* Map all parameters in sky & paramx array to `pall` array (which is relative to geocenter
* and holds sky params too). Version of pmap that updates extrinsic projection (Fcross/Fplus/dtimes).
*/
void pmap_all(struct Net *net, double *pall, double *param, double *sky, int NX, double gmst)
{
int j;
double ciota, Ap, Ac, alpha, sindelta, psi, Fs;
double *dtimes;
double Fplus, Fcross, lambda;
// sky [0] alpha, [1] sin(delta) [2] psi [3] ciota [4] scale [5] dphi [6] dt
// param [0] log(Mc) [1] log(Mt) [2] chi1 [3] chi2 [4] phi0 [5] tp0 [6] log(DL0) then relative amplitudes, time, phases
// pall [0] log(Mc) [1] log(Mt) [2] chi1 [3] chi2 [4] 2phi0 [5] tp [6] log(DL) [7] alpha [8] sindelta [9] psi [10] ciota
dtimes = (double*)malloc(sizeof(double)* 5);
for (j = 0; j < NX; ++j)
{
pall[j] = param[j];
}
ciota = sky[3];
Ap = (1.0+ciota*ciota)/2.0;
Ac = ciota;
alpha = sky[0];
sindelta = sky[1];
psi = sky[2];
Times(alpha, sindelta, net->dtimes, gmst);
F_ant(psi, alpha, sindelta, &net->Fplus[net->labels[0]], &net->Fcross[net->labels[0]], net->labels[0], gmst);
Fs = sqrt(Ap*Ap*net->Fplus[net->labels[0]]*net->Fplus[net->labels[0]]+Ac*Ac*net->Fcross[net->labels[0]]*net->Fcross[net->labels[0]]);
lambda = atan2(Ac*net->Fcross[net->labels[0]],Ap*net->Fplus[net->labels[0]]);
if(lambda < 0.0) lambda += TPI;
// move reference point from ref detector to geocenter
pall[4] = 2.0*param[4] - lambda; // I'm holding the GW phase in pall[4]
pall[5] = param[5] - net->dtimes[net->labels[0]];
pall[6] = param[6] + log(Fs);
pall[NX] = alpha;
pall[NX+1] = sindelta;
pall[NX+2] = psi;
pall[NX+3] = ciota;
free(dtimes);
}
void pmap_back(struct Net *net, double *pall, double *param, double *sky, int NX, double gmst)
{
int j;
......@@ -1901,90 +1852,6 @@ void SNRvsf(struct Net *net, double **D, double *params, RealVector *freq, doubl
}
double log_likelihood_full(struct Net *net, double **D, double *params, RealVector *freq, double **SN, double *rho, int N, double Tobs, struct bayesCBC *rundata)
{
int id;
double logL, logLdet;
double **twave;
double HH, HD, x;
int imin, imax;
double fmin, fmax;
double *ampphase;
double *sky;
double ciota, Ap, Ac, alpha, sindelta, psi;
int extrinsic_update_flag;
extrinsic_update_flag = 0;
fmin = rundata->fmin;
fmax = rundata->fmax;
logL = 0.0;
if(rundata->constantLogLFlag == 0)
{
imin = (int)(Tobs*fmin);
imax = (int)(Tobs*fmax);
if(imax > N/2) imax = N/2;
twave = double_matrix(net->Nifo,N);
ampphase = double_vector(N);
sky = double_vector(rundata->NS);
ciota = params[rundata->NX+3];
Ap = (1.0+ciota*ciota)/2.0;
Ac = ciota;
alpha = params[rundata->NX];
sindelta = params[rundata->NX+1];
psi = params[rundata->NX+2];
sky[0] = alpha;
sky[1] = sindelta;
sky[2] = psi;
sky[3] = ciota;
sky[4] = 1.0;
sky[5] = 0.0;
sky[6] = 0.0;
// Update dtimes, projections (if wanted - currently always done)
if (extrinsic_update_flag == 0)
{
Times(alpha, sindelta, net->dtimes, rundata->gmst);
for(id = 0; id < net->Nifo; id++) F_ant(psi, alpha, sindelta, &net->Fplus[id], &net->Fcross[id], net->labels[id], rundata->gmst);
}
// Compute template with current parameters
geowave(ampphase, freq, params, N, rundata);
projectCBCWaveform(ampphase, N, net->Nifo, fmin, Tobs, sky, twave, net->dtimes, net->Fplus, net->Fcross);
logL = 0.0;
for(id = 0; id < net->Nifo; id++)
{
HH = fourier_nwip_cbc(twave[id], twave[id], SN[id], imin, imax, N);
HD = fourier_nwip_cbc(D[id], twave[id], SN[id], imin, imax, N);
rho[id] = HD/sqrt(HH);
logLdet = HD-0.5*HH;
logL += logLdet;
// Store likelihood for individual detectors for BW
rundata->detLogL[id] = logLdet;
}
free_double_matrix(twave,net->Nifo);
free_double_vector(ampphase);
free_double_vector(sky);
}
return(logL);
}
double log_likelihood_max(struct Net *net, double **D, double *params, RealVector *freq, double **SN, int N, double Tobs, double tmin, double tmax, struct bayesCBC *rundata)
{
int i, j, k, NMAX;
......@@ -2629,26 +2496,26 @@ void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, dou
{
logL = log_likelihood(net, D, paramx[q], freq, SN, N, Tobs, rundata);
if (NX > 7) {
if (lambda_type_version == lambdaTilde) {
lambdaT = paramx[q][7];
dLambdaT = paramx[q][8];
LambdaTsEta2Lambdas(lambdaT, dLambdaT, eta, &lambda1, &lambda2);
} else {
lambda1 = paramx[q][7];
lambda2 = paramx[q][8];
}
fprintf(chain,"%d %e %e %e %e %e %e %e %e %e ", mc/10, logLx[q], logL, Mchirp, Mtot, chieff, m1, m2, lambda1, lambda2);
} else {
fprintf(chain,"%d %e %e %e %e %e %e %e ", mc/10, logLx[q], logL, Mchirp, Mtot, chieff, m1, m2);
}
// if (NX > 7) {
// if (lambda_type_version == lambdaTilde) {
// lambdaT = paramx[q][7];
// dLambdaT = paramx[q][8];
// LambdaTsEta2Lambdas(lambdaT, dLambdaT, eta, &lambda1, &lambda2);
// } else {
// lambda1 = paramx[q][7];
// lambda2 = paramx[q][8];
// }
// fprintf(chain,"%d %e %e %e %e %e %e %e %e %e ", mc/10, logLx[q], logL, Mchirp, Mtot, chieff, m1, m2, lambda1, lambda2);
// } else {
// fprintf(chain,"%d %e %e %e %e %e %e %e ", mc/10, logLx[q], logL, Mchirp, Mtot, chieff, m1, m2);
// }
for(i = NX; i < NX+3*net->Nifo; i++) fprintf(chain,"%e ", paramx[q][i]);
fprintf(chain,"\n");
// for(i = NX; i < NX+3*net->Nifo; i++) fprintf(chain,"%e ", paramx[q][i]);
// fprintf(chain,"\n");
}
else
{
fprintf(chain,"%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n", mxc[0], logLx[q], Mchirp, Mtot, chieff, pallx[q][4], \
// fprintf(chain,"%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n", mxc[0], logLx[q], Mchirp, Mtot, chieff, pallx[q][4], \
pallx[q][5]-(Tobs-2.0), DL, skyx[q][0], skyx[q][1], \
skyx[q][2], ciota, z, Mchirp/(1.0+z), Mtot/(1.0+z), m1/(1.0+z), m2/(1.0+z), m1, m2);
mxc[0] += 1;
......@@ -4594,8 +4461,9 @@ void projectCBCWaveform(double *ampphase, int N, int NI, double fmin, double Tob
double Fs, lambda;
double A, x, td, f;
double epsilon = extParams[3];
double ciota =(1.0-sqrt(1.0-epsilon*epsilon))/epsilon;
double ecc = extParams[3];
double phi0 = extParams[4];
double ciota =(1.0-sqrt(1.0-ecc*ecc))/ecc;
double Ap = (1.0+ciota*ciota)/2.0;
double Ac = ciota;
......@@ -4606,6 +4474,9 @@ void projectCBCWaveform(double *ampphase, int N, int NI, double fmin, double Tob
{
Fs = sqrt(Ap*Ap*Fplus[id]*Fplus[id]+Ac*Ac*Fcross[id]*Fcross[id]); // magnitude of response
lambda = atan2(Ac*Fcross[id],Ap*Fplus[id]);
// Fs = sqrt(Fplus[id]*Fplus[id]+ecc*ecc*Fcross[id]*Fcross[id]);
// lambda = atan2(ecc*Fcross[id],Fplus[id]);
// td = -1*dtimes[id];
td = dtimes[id];
if(lambda < 0.0) lambda += TPI;
......
......@@ -215,7 +215,7 @@ void skymcmc(struct Net *net, int MCX, int *mxc, FILE *chain, double **paramx, d
alpha = gsl_rng_uniform(r);
if(alpha > 0.8 && net->Nifo > 1) // ring
if(alpha > 1.0 && net->Nifo > 1) // ring
{
......@@ -356,41 +356,41 @@ void skymcmc(struct Net *net, int MCX, int *mxc, FILE *chain, double **paramx, d
} // ends choice of update
// if(mc%100 == 0)
// {
// ic = who[1];
// phi = skyx[ic][0];
// if(phi > TPI) phi -= TPI;
// if(phi < 0.0) phi += TPI;
// skyx[ic][0] = phi;
// Mchirp = exp(paramx[ic][0])/MSUN_SI;
// Mtot = exp(paramx[ic][1])/MSUN_SI;
// eta = pow((Mchirp/Mtot), (5.0/3.0));
// dm = sqrt(1.0-4.0*eta);
// m1 = Mtot*(1.0+dm)/2.0;
// m2 = Mtot*(1.0-dm)/2.0;
// chieff = (m1*paramx[ic][2]+m2*paramx[ic][3])/Mtot;
if(mc%100 == 0)
{
ic = who[1];
phi = skyx[ic][0];
if(phi > TPI) phi -= TPI;
if(phi < 0.0) phi += TPI;
skyx[ic][0] = phi;
Mchirp = exp(paramx[ic][0])/MSUN_SI;
Mtot = exp(paramx[ic][1])/MSUN_SI;
eta = pow((Mchirp/Mtot), (5.0/3.0));
dm = sqrt(1.0-4.0*eta);
m1 = Mtot*(1.0+dm)/2.0;
m2 = Mtot*(1.0-dm)/2.0;
chieff = (m1*paramx[ic][2]+m2*paramx[ic][3])/Mtot;
// // counter, log likelihood, chirp mass, total mass, effective spin, phase shift , time shift, distance, RA , sine of DEC,
// // polarization angle, cos inclination
// counter, log likelihood, chirp mass, total mass, effective spin, phase shift , time shift, distance, RA , sine of DEC,
// polarization angle, cos inclination
// DL = exp(pallx[ic][6])/(1.0e6*PC_SI*skyx[ic][4]);
// z = z_DL(DL);
DL = exp(pallx[ic][6])/(1.0e6*PC_SI*skyx[ic][4]);
z = z_DL(DL);
// // Note that skyx[ic][5], skyx[ic][6], hold different quantities than what is printed by the other MCMCs
// Note that skyx[ic][5], skyx[ic][6], hold different quantities than what is printed by the other MCMCs
// // fprintf(rundata->chainS,"%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n", mxc[1], logLx[ic], Mchirp, Mtot, chieff, skyx[ic][5], \
// skyx[ic][6], DL, skyx[ic][0], skyx[ic][1], \
// skyx[ic][2], skyx[ic][3], z, Mchirp/(1.0+z), Mtot/(1.0+z), m1/(1.0+z), m2/(1.0+z), m1, m2);
fprintf(rundata->chainS,"%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n", mxc[1], logLx[ic], Mchirp, Mtot, chieff, skyx[ic][5], \
skyx[ic][6], DL, skyx[ic][0], skyx[ic][1], \
skyx[ic][2], skyx[ic][3], z, Mchirp/(1.0+z), Mtot/(1.0+z), m1/(1.0+z), m2/(1.0+z), m1, m2);
// mxc[1] += 1;
mxc[1] += 1;
// }
}
if(mc%10000 == 0)
{
......@@ -534,7 +534,6 @@ double skylike(struct Net *net, double *params, double *D, double *H, double **D
for(id = 0; id<net->Nifo; id++)
{
ComputeDetFant(net, psi, alpha, sindelta, &Fplus, &Fcross, net->labels[id], rundata->gmst);
// XLALComputeDetAMResponse(&Fplus, &Fcross, (const REAL4(*)[3]) net->response[id], alpha, asin(sindelta), psi, rundata->gmst);
F[id] = sqrt(Ap*Ap*Fplus*Fplus+Ac*Ac*Fcross*Fcross); // magnitude of response
lambda[id] = atan2(Ac*Fcross,Ap*Fplus);
if(lambda[id] < 0.0) lambda[id] += TPI;
......@@ -1245,6 +1244,143 @@ void Ring(struct Net *net, double *skyx, double *skyy, int d1, int d2, gsl_rng *
}
/**
* Map all parameters in sky & paramx array to `pall` array (which is relative to geocenter
* and holds sky params too). Version of pmap that updates extrinsic projection (Fcross/Fplus/dtimes).
*/
void pmap_all(struct Net *net, double *pall, double *param, double *sky, int NX, double gmst)
{
int j;
double ciota, Ap, Ac, alpha, sindelta, psi, Fs;
double *dtimes;
double Fplus, Fcross, lambda;
// sky [0] alpha, [1] sin(delta) [2] psi [3] ciota [4] scale [5] dphi [6] dt
// param [0] log(Mc) [1] log(Mt) [2] chi1 [3] chi2 [4] phi0 [5] tp0 [6] log(DL0) then relative amplitudes, time, phases
// pall [0] log(Mc) [1] log(Mt) [2] chi1 [3] chi2 [4] 2phi0 [5] tp [6] log(DL) [7] alpha [8] sindelta [9] psi [10] ciota
dtimes = (double*)malloc(sizeof(double)* 5);
for (j = 0; j < NX; ++j)
{
pall[j] = param[j];
}
ciota = sky[3];
Ap = (1.0+ciota*ciota)/2.0;
Ac = ciota;
alpha = sky[0];
sindelta = sky[1];
psi = sky[2];
TimeDelays(net, alpha, sindelta, net->dtimes, gmst);
ComputeDetFant(net, psi, alpha, sindelta, &net->Fplus[net->labels[0]], &net->Fcross[net->labels[0]], net->labels[0], gmst);
Fs = sqrt(Ap*Ap*net->Fplus[net->labels[0]]*net->Fplus[net->labels[0]]+Ac*Ac*net->Fcross[net->labels[0]]*net->Fcross[net->labels[0]]);
lambda = atan2(Ac*net->Fcross[net->labels[0]],Ap*net->Fplus[net->labels[0]]);
if(lambda < 0.0) lambda += TPI;
// move reference point from ref detector to geocenter
pall[4] = 2.0*param[4] - lambda; // I'm holding the GW phase in pall[4]
pall[5] = param[5] - net->dtimes[net->labels[0]];
pall[6] = param[6] + log(Fs);
pall[NX] = alpha;
pall[NX+1] = sindelta;
pall[NX+2] = psi;
pall[NX+3] = ciota;
free(dtimes);
}
double log_likelihood_full(struct Net *net, double **D, double *params, RealVector *freq, double **SN, double *rho, int N, double Tobs, struct bayesCBC *rundata)
{
int id;
double logL, logLdet;
double **twave;
double HH, HD, x;
int imin, imax;
double fmin, fmax;
double *ampphase;
double *sky;
double ciota, Ap, Ac, alpha, sindelta, psi;
int extrinsic_update_flag;
extrinsic_update_flag = 0;
fmin = rundata->fmin;
fmax = rundata->fmax;
logL = 0.0;
if(rundata->constantLogLFlag == 0)
{
imin = (int)(Tobs*fmin);
imax = (int)(Tobs*fmax);
if(imax > N/2) imax = N/2;
twave = double_matrix(net->Nifo,N);
ampphase = double_vector(N);
sky = double_vector(rundata->NS);
ciota = params[rundata->NX+3];
Ap = (1.0+ciota*ciota)/2.0;
Ac = ciota;
alpha = params[rundata->NX];
sindelta = params[rundata->NX+1];
psi = params[rundata->NX+2];
sky[0] = alpha;
sky[1] = sindelta;
sky[2] = psi;
sky[3] = ciota;
sky[4] = 1.0;
sky[5] = 0.0;
sky[6] = 0.0;
// Update dtimes, projections (if wanted - currently always done)
if (extrinsic_update_flag == 0)
{
TimeDelays(net, alpha, sindelta, net->dtimes, rundata->gmst);
for(id = 0; id < net->Nifo; id++) ComputeDetFant(net, psi, alpha, sindelta, &net->Fplus[id], &net->Fcross[id], net->labels[id], rundata->gmst);
}
sky[3] = 2*ciota/(1+ciota*ciota);
// Compute template with current parameters
geowave(ampphase, freq, params, N, rundata);
projectCBCWaveform(ampphase, N, net->Nifo, fmin, Tobs, sky, twave, net->dtimes, net->Fplus, net->Fcross);
logL = 0.0;
for(id = 0; id < net->Nifo; id++)
{
HH = fourier_nwip_cbc(twave[id], twave[id], SN[id], imin, imax, N);
HD = fourier_nwip_cbc(D[id], twave[id], SN[id], imin, imax, N);
rho[id] = HD/sqrt(HH);
logLdet = HD-0.5*HH;
logL += logLdet;
// Store likelihood for individual detectors for BW
rundata->detLogL[id] = logLdet;
}
sky[3] = ciota;
free_double_matrix(twave,net->Nifo);
free_double_vector(ampphase);
free_double_vector(sky);
}
return(logL);
}
/* ********************************************************************************** */
/* */
/* Wrappers for LAL */
......@@ -1266,6 +1402,7 @@ void TimeDelays(struct Net *net, double alpha, double sindelta, double *dtimes,
for(ifo=0; ifo<net->Nifo; ifo++)
{
dtimes[ifo] = XLALTimeDelayFromEarthCenter(net->location[ifo], alpha, asin(sindelta), XLALGreenwichMeanSiderealTimeToGPS((REAL8) gmst, &gps));
dtimes[ifo] = -1*dtimes[ifo];
}
}
......@@ -1273,16 +1410,17 @@ void ComputeDetFant(struct Net *net, double psi, double alpha, double sindelta,
{
int i, j;
REAL4 response[3][3];
double Fplus_temp, Fcross_temp;
for(i=0; i<3; i++)
{
for(j=0; j<3; j++) response[i][j] = (REAL4) net->response[id][i][j];
}
// for(i=0; i<3; i++) response[i] = (REAL4[]) net->response[id][i];
XLALComputeDetAMResponse(Fplus, Fcross, (const REAL4(*)[3])response, alpha, asin(sindelta), psi, gmst);
// printf("ComputeDetFant Fplus %f Fcross %f\n", Fplus, Fcross);
// printf("0 - ComputeDetFant Fplus %f Fcross %f\n", Fplus, Fcross);
F_ant(psi, alpha, sindelta, &Fplus_temp, &Fcross_temp, net->labels[id], gmst);
// printf("1 - ComputeDetFant Fplus %f Fcross %f\n", Fplus_temp, Fcross_temp);
}
......
......@@ -1020,6 +1020,13 @@ void print_chain_files(struct Data *data, struct Chain *chain, struct Model **mo
//print bayescbc parameters
if(data->cbcFlag)
{
for(int ifo = 0; ifo<data->NI; ifo++)
{
// Antenna patterns for detectors
bayescbc->net->Fplus[ifo] = model_x->projection->Fplus[ifo];
bayescbc->net->Fcross[ifo] = model_x->projection->Fcross[ifo];
}
print_cbc_model(chain->cbcChainFile[ic], bayescbc, chain->index[ic], data->Tobs);
}
......@@ -1146,9 +1153,9 @@ void print_cbc_model(FILE *fptr, struct bayesCBC *bayescbc, int ic, double Tobs)
// [0] log(Mc) [1] log(Mt) [2] chi1 [3] chi2 [4] phi0 [5] tp [6] log(DL) [7] alpha [8] sindelta [9] psi [10] ciota
// counter, log likelihood, chirp mass (DF), total mass (DF), effective spin, geocenter orbital phase, geocenter arrival time, distance, RA , sine of DEC,
// polarization angle, cos inclination, redshift, chirp mass (SF), total mass (SF), m1 (SF), m2 (SF), m1 (DF), m2 (DF), arrival time at each detector, SNR at each detector
fprintf(bayescbc->chainA,"%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e ", bayescbc->mxc[2], bayescbc->logLx[ic], Mchirp, Mtot, chieff, bayescbc->pallx[ic][4], \
fprintf(bayescbc->chainA,"%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e", bayescbc->mxc[2], bayescbc->logLx[ic], Mchirp, Mtot, chieff, bayescbc->pallx[ic][4], \
tc, DL, bayescbc->pallx[ic][bayescbc->NX], bayescbc->pallx[ic][bayescbc->NX+1], \
bayescbc->pallx[ic][bayescbc->NX+2], bayescbc->pallx[ic][bayescbc->NX+3], z, Mchirp/(1.0+z), Mtot/(1.0+z), m1/(1.0+z), m2/(1.0+z), m1, m2);
bayescbc->pallx[ic][bayescbc->NX+2], bayescbc->pallx[ic][bayescbc->NX+3], z, Mchirp/(1.0+z), Mtot/(1.0+z), m1/(1.0+z), m2/(1.0+z), m1, m2, bayescbc->net->Fplus[0], bayescbc->net->Fcross[0]);
fprintf(bayescbc->chainA,"\n");
......
......@@ -849,6 +849,7 @@ double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, doubl
//double sum = 0.0;
double logL = 0.0;
double **h, **d, **r, **t, **c;
double *dtimes;
if(extrinsic_checkrange(params)) return -1.0e60;
......@@ -886,7 +887,10 @@ double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, doubl
if(data->cbcFlag)
{
// for(i=0; i<NI; i++) for(n=0; n<N; n++) c[i][n] = 0.0;
projectCBCWaveform(amphase, N, NI, fmin, data->Tobs, params, c, projection->dtimes, projection->Fplus, projection->Fcross);
dtimes = double_vector(NI);
for(ifo=0; ifo<NI; ifo++) dtimes[ifo] = -1*projection->dtimes[ifo];
projectCBCWaveform(amphase, N, NI, fmin, data->Tobs, params, c, dtimes, projection->Fplus, projection->Fcross);
free_double_vector(dtimes);
}
//Form up residual
......
......@@ -1255,8 +1255,9 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
fprintf(stdout, " ======= Start CBC_burnin() ======\n");
fprintf(stdout, " Start #1 intrinsic burnin \n");
bayescbc->chainS = fopen("chains/cbc_searchchain.dat","w");
bayescbc->gmst = XLALGreenwichMeanSiderealTime(&data->epoch);
burnin_bayescbc(data->N, data->Tobs, data->trigtime, bayescbc->D, bayescbc->SN, chain->seed, filename, 0, bayescbc);
fclose(bayescbc->chainS);
// Quick search to find initial signal extrinisc parameters
// Currently is not executed with sky fixed flag (assuming sky paramters are given)
......@@ -1266,7 +1267,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
fprintf(stdout, " Start #2 (optional) extrinisic burnin \n");
extrinsic_burnin_cbc(bayescbc, data, model, chain);
}
fclose(bayescbc->chainS);
// Pass updated cbc parameters to model
for(ic=0; ic<chain->NC; ic++)
{
......@@ -1277,7 +1278,6 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
if (bayescbc->intrinsic_only == 0)
{
printf("updating model with burnin extrinisc params \n");
// Pass extrinsic parameters from BW model to BayesCBC
for(i = bayescbc->NX; i < bayescbc->NP; i++)
{
......@@ -2603,6 +2603,7 @@ void EvolveBayesCBCParameters(struct Data *data, struct Model **model, struct ba
int imin, imax;
double x, y;
double ciota, epsilon, scale;
double *dtimes;
char filename[MAXSTRINGSIZE];
int M = 500; // chain->cycle; // Nr. of local MCMC iterations
......@@ -2647,13 +2648,17 @@ void EvolveBayesCBCParameters(struct Data *data, struct Model **model, struct ba
model_x->cbcParams[i] = bayescbc->pallx[icbc][i];
}
dtimes = double_vector(NI);
for(ifo=0; ifo<NI; ifo++) dtimes[ifo] = -1 * projection->dtimes[ifo];
// Get CBC waveforms projected on the detectors
projectCBCWaveform(model_x->cbcamphase, N, NI, fmin, data->Tobs, model_x->extParams, model_x->cbctemplate, projection->dtimes, projection->Fplus, projection->Fcross);
projectCBCWaveform(model_x->cbcamphase, N, NI, fmin, data->Tobs, model_x->extParams, model_x->cbctemplate, dtimes, projection->Fplus, projection->Fcross);
free_double_vector(dtimes);
// Save 200 waveforms temporarily for checking residuals + template
//if (ic==0 && bayescbc->debug == 1 && chain->mc%(chain->count/200)==0) print_projected_cbc_waveform(bayescbc->SN, data->Tobs, data->trigtime, model_x->cbctemplate, bayescbc->D, N, bayescbc->mxc[0], bayescbc);
// Recompute likelihoods of current chain
model_x->logLnorm = 0.0;
model_x->logL = 0.0;
......@@ -2904,7 +2909,14 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
waveformProject(data, model_x->projection, model_x->extParams, model_x->response, model_x->h, data->fmin, data->fmax);
Shf_Geocenter_full(data, model_x->projection, model_x->Snf, model_x->SnGeo, model_x->extParams);
if (data->cbcFlag) projectCBCWaveform(model_x->cbcamphase, data->N, NI, data->fmin, data->Tobs, model_x->extParams, model_x->cbctemplate, model_x->projection->dtimes, model_x->projection->Fplus, model_x->projection->Fcross);
if (data->cbcFlag)
{
double *dtimes;
dtimes = double_vector(NI);
for(ifo=0; ifo<NI; ifo++) dtimes[ifo] = -1 * model_x->projection->dtimes[ifo];
projectCBCWaveform(model_x->cbcamphase, data->N, NI, data->fmin, data->Tobs, model_x->extParams, model_x->cbctemplate, dtimes, model_x->projection->Fplus, model_x->projection->Fcross);
free_double_vector(dtimes);
}
free(glitch);
......@@ -3073,6 +3085,12 @@ void extrinsic_burnin_cbc(struct bayesCBC *bayescbc, struct Data *data, struct M
}
printf("before skymcmc\n");
for(i=0; i<NC; i++)
{
double DL = exp(pallx[i][6])/(1.0e6*PC_SI);
printf("DL[%i][6] = %f, %f\n", i, DL, pallx[i][j]);
}
// intialize the local heat array
// run cold to force ML
for(i=0; i<NCC; i++) heat[i] = 0.25;
......@@ -3144,7 +3162,6 @@ void extrinsic_burnin_cbc(struct bayesCBC *bayescbc, struct Data *data, struct M
free_double_tensor(DHs,NC,net->Nifo);