Commit 542dfb34 authored by Marcella Wijngaarden's avatar Marcella Wijngaarden

Fix projected templates from BayesCBC

There was a mismatch in the templates used internally in
bayesCBC and the one used in BayesWave to compute the residuals
(and likelihood). This commit updates the way the projection
is calculated (and this fixes the previous issue with parallel
tempering because now the likelihood is computed correctly..)
parent 9c91308a
Pipeline #159227 failed with stages
in 41 seconds
......@@ -424,6 +424,125 @@ void wavemax(struct Net *net, int N, double **tp, RealVector *freq, double **par
}
void print_projected_cbc_waveform(double **SN, double Tobs, double ttrig, double **projected, double **data, int N, int iter, struct bayesCBC *rundata)
{
char command[1024];
FILE *out;
double fs;
int i, id;
double fac, x, y, f, fr, dt;
double fmin, fmax;
double **twave, **dwave;
double **ww, **dw;
struct Net *net = rundata->net;
int Nifo = net->Nifo;
RealVector *freq = rundata->freq;
fmin = rundata->fmin;
fmax = rundata->fmax;
dt = Tobs/(double)(N);
twave = double_matrix(net->Nifo,N);
ww = double_matrix(net->Nifo,N); // Template/model whitened
dwave = double_matrix(net->Nifo,N);
dw = double_matrix(net->Nifo,N); // Data whitened
fs = fbegin(rundata->pallx[rundata->who[0]]);
// Convert indices for FFT
for (id = 0; id < Nifo; ++id)
{
for (i = 1; i < N/2; ++i)
{
twave[id][i] = projected[id][2*i];
twave[id][N-i] = projected[id][2*i+1];
dwave[id][i] = data[id][2*i];
dwave[id][N-i] = data[id][2*i+1];
}
twave[id][0] = 0.0;
twave[id][N/2] = 0.0;
dwave[id][0] = 0.0;
dwave[id][N/2] = 0.0;
}
// MW: For plotting was fr/3.0, but this also has better overlap at low frequencies
fr = fs/fmin;
for (i = 1; i < N/2; ++i)
{
f = freq->data[i];
x = 0.5*(1.0+tanh((f-(fs+fr))/fr));
for (id = 0; id < net->Nifo; ++id)
{
y = 1.0/sqrt(SN[id][i]*2.0);
twave[id][i] *= x;
twave[id][N-i] *= x;
dwave[id][i] *= x;
dwave[id][N-i] *= x;
ww[id][i] = y*twave[id][i];
ww[id][N-i] = y*twave[id][N-i];
dw[id][i] = y*dwave[id][i];
dw[id][N-i] = y*dwave[id][N-i];
}
}
x = sqrt(2.0)/dt;
fac = sqrt((double)(N/2));
for (id = 0; id < Nifo; ++id)
{
ww[id][0] = 0.0;
ww[id][N/2] = 0.0;
dw[id][0] = 0.0;
dw[id][N/2] = 0.0;
gsl_fft_halfcomplex_radix2_inverse(ww[id], 1, N);
gsl_fft_halfcomplex_radix2_inverse(dw[id], 1, N);
for (i = 0; i < N; ++i)
{
ww[id][i] *= fac;
dw[id][i] *= fac;
}
}
fac = ceil(sqrt(1.0/4.0));
for (id = 0; id < net->Nifo; ++id)
{
sprintf(command, "wavesCBC/projected_wavewhite_%d_%d_%d_%d.dat", iter, (int)(Tobs), (int)ttrig, net->labels[id]);
out = fopen(command,"w");
for (i = 0; i < N; ++i)
{
fprintf(out,"%e %e\n", (double)(i)*dt-Tobs+2.0, ww[id][i]/fac);
}
fclose(out);
sprintf(command, "wavesCBC/projected_wavedatawhite_%d_%d_%d_%d.dat", iter, (int)(Tobs), (int)ttrig, net->labels[id]);
out = fopen(command,"w");
for (i = 0; i < N; ++i)
{
fprintf(out,"%e %.16e\n", (double)(i)*dt-Tobs+2.0, dw[id][i]/fac);
}
fclose(out);
}
free_double_matrix(dwave,net->Nifo);
free_double_matrix(twave,net->Nifo);
free_double_matrix(ww,net->Nifo);
free_double_matrix(dw,net->Nifo);
}
void printwaveall(struct Net *net, int N, RealVector *freq, double *paramx, double **SN, double Tobs, double ttrig, int iter, struct bayesCBC *rundata)
......@@ -637,6 +756,7 @@ void printwaveall(struct Net *net, int N, RealVector *freq, double *paramx, doub
free_double_matrix(thold,net->Nifo);
free_double_matrix(twave,net->Nifo);
free_double_matrix(dwave,net->Nifo);
free_double_matrix(twave_tmp,net->Nifo);
free_double_matrix(Dtime,net->Nifo);
free_double_matrix(tf,net->Nifo);
......@@ -3238,9 +3358,9 @@ void set_bayescbc_priors(struct Net *net, struct bayesCBC *rundata)
void MCMC_all_wrapper(int M, int N, double Tobs, double ttrig, gsl_rng *r, int istart, struct bayesCBC *rundata)
{
int i, k;
double **cbctemplate;
double *amphase;
cbctemplate = double_matrix(rundata->net->Nifo-1,N-1);
amphase = double_vector(N-1);
//open MP modifications
omp_set_num_threads(rundata->NC);
......@@ -3249,17 +3369,17 @@ void MCMC_all_wrapper(int M, int N, double Tobs, double ttrig, gsl_rng *r, int i
for(k=0; k < rundata->NC; k++)
{
printf("Updating chain %i/%i\n", k, rundata->NC);
BayesCBC_MCMC(M, N, Tobs, ttrig, rundata->rvec[k], istart, rundata, k, cbctemplate);
BayesCBC_MCMC(M, N, Tobs, ttrig, rundata->rvec[k], istart, rundata, k, amphase);
}
free(cbctemplate);
free(amphase);
return;
}
/**
* Run MCMC_all for a single chain with the full parameter array (NP)
*/
void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int istart, struct bayesCBC *rundata, int ic, double **template)
void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int istart, struct bayesCBC *rundata, int ic, double *amphase)
{
int i, j, q, k, mc, flag;
int m, *evolveparams;
......@@ -3270,7 +3390,6 @@ void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int ista
double *rhox;
double mapL;
double *pmax;
double *geocbc;
// Tidal parameters
double lambda1, lambda2, lambdaT, dLambdaT;
......@@ -3313,7 +3432,6 @@ void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int ista
evec = double_matrix(NP,NP);
pmax = double_vector(NP);
paramy = double_vector(NP);
geocbc = double_vector(N);
lambda_type_version=rundata->lambda_type_version;
......@@ -3350,7 +3468,7 @@ void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int ista
// Initialize the cbc likelihood
rundata->logLx[icbc] = log_likelihood_full(net, rundata->D, paramx, freq, rundata->SN, rhox, N, Tobs, rundata);
if(rundata->constantLogLFlag == 0) rundata->logLx[icbc] = log_likelihood_full(net, rundata->D, paramx, freq, rundata->SN, rhox, N, Tobs, rundata);
printf("logLx at start %i: %f \n", icbc, rundata->logLx[icbc]);
// printf("logLx at start %i: %f \n", icbc, rundata->logLx[icbc]);
Fisher_All(net, fish, paramx, freq, rundata->SN, N, Tobs, rundata);
......@@ -3439,28 +3557,9 @@ void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int ista
// Save MAP likelihood. TODO: Do we want to have this?
// for(i = 0; i < NP; i++) rundata->pallx[ic][i] = pmax[i];
// Store geotemplate
// Store Amplitude and Phase for geotemplate
// fulltemplates(rundata->net, template, rundata->freq, rundata->pallx[ic], N, rundata);
geotemplate(geocbc, rundata->freq, rundata->pallx[icbc], N, rundata);
ciota = rundata->pallx[icbc][NX+3];
double epsilon =2*ciota/(1+ciota*ciota);
int re, im;
for (i=0; i<N/2; ++i) {
re = 2*i;
im = re+1;
//h+
template[0][re] = geocbc[re];
template[0][im] = geocbc[im];
//hx = -i * epsilon * h+
template[1][re] = -epsilon*template[0][im];
template[1][im] = epsilon*template[0][re];
}
// printf("logLx at end update %i: %f \n", ic, rundata->logLx[ic]);
geowave(amphase, rundata->freq, rundata->pallx[icbc], N, rundata);
free_int_vector(av);
free_int_vector(cv);
......@@ -3472,7 +3571,6 @@ void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int ista
free_double_matrix(evec,NP);
free_double_vector(paramy);
free_double_vector(geocbc);
}
......@@ -4502,6 +4600,151 @@ void cbcGeoTemplate(struct bayesCBC *rundata, int ic, int N, double **template)
free_double_vector(geocbc);
}
void projectCBCWaveform(double *ampphase, int N, int NI, double fmin, double Tobs, double *extParams, double **response, double *dtimes, double *Fplus, double *Fcross)
{
int i,id;
double Fs, lambda;
double A, x, td, f;
double epsilon = extParams[3];
double ciota =(1.0-sqrt(1.0-epsilon*epsilon))/epsilon;
double Ap = (1.0+ciota*ciota)/2.0;
double Ac = ciota;
double df = 1./Tobs;
for (id=0; id<NI; id++)
{
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]);
td = dtimes[id];
if(lambda < 0.0) lambda += TPI;
response[id][0] = 0.0;
response[id][1] = 0.0;
for (i=1; i< N/2; i++)
{
f = i*df;
if (f < fmin) {
response[id][2*i] = 0.0;
response[id][2*i+1] = 0.0;
}
else
{
A = ampphase[2*i]*Fs;
x = ampphase[2*i+1]+lambda-TPI*i*df*td;
response[id][2*i] = A*cos(x);
response[id][2*i+1] = A*sin(x);
}
}
}
printf("done the loop!\n");
}
void geowave(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[1] = 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[2*i] = A;
gwave[2*i+1] = x;
}
else
{
gwave[2*i] = 0.0;
gwave[2*i+1] = 0.0;
}
}
DestroyAmpPhaseFDWaveform(ap);
free(extraParams);
}
void geotemplate(double *gwave, RealVector *freq, double *params, int N, struct bayesCBC *rundata)
{
......
......@@ -114,9 +114,11 @@ void CBC_initialize(struct Net *net, int *mxc, FILE *chainS, double **paramx, do
void MCMC_all(struct Net *net, int *mxc, int M, FILE *chain, double **paramx, int *who, double *heat, double ***history, double ***global, RealVector *freq, double **D, double **SN, int N, double Tobs, double ttrig, gsl_rng *r, int istart, struct bayesCBC *rundata);
void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int istart, struct bayesCBC *rundata, int ic, double **template);
void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int istart, struct bayesCBC *rundata, int ic, double *amphase);
void MCMC_all_wrapper(int M, int N, double Tobs, double ttrig, gsl_rng *r, int istart, struct bayesCBC *rundata);
void cbcGeoTemplate(struct bayesCBC *rundata, int ic, int N, double **template);
void projectCBCWaveform(double *AmpPhase, int N, int NI, double fmin, double Tobs, double *extParams, double **response, double *dtimes, double *Fplus, double *Fcross);
void print_projected_cbc_waveform(double **SN, double Tobs, double ttrig, double **projected, double **data, int N, int iter, struct bayesCBC *rundata);
double z_DL(double DL);
double DL_fit(double z);
......
......@@ -112,6 +112,7 @@ void extrinsictemplates(struct Net *net, double **hwave, RealVector *freq, doubl
void templates(struct Net *net, double **hwave, RealVector *freq, double *params, int N, struct bayesCBC *rundata);
void geotemplate(double *gwave, RealVector *freq, double *params, int N, struct bayesCBC *rundata);
void geowave(double *gwave, RealVector *freq, double *params, int N, struct bayesCBC *rundata);
void fulltemplates(struct Net *net, double **hwave, RealVector *freq, double *params, int N, struct bayesCBC *rundata);
void fullphaseamp(struct Net *net, double **amp, double **phase, RealVector *freq, double *params, int N, struct bayesCBC *rundata);
void Fisher_All(struct Net *net, double **fish, double *params, RealVector *freq, double **SN, int N, double Tobs, struct bayesCBC *rundata);
......
......@@ -129,7 +129,7 @@ struct Model
double **h; // reference waveform at one interferometer (NpolS times N)
double **response; // response in each detector (NI times N)
double **cbctemplate; // CBC model for each detector (NI times N) TODO: rename to cbcresponse
double **cbcgeotemplate; // CBC model for each detector (NP times N)
double *cbcamphase; // CBC geocenter amplitude & phase
int chirpletFlag; // whether to use chirplets
int NW; // number of intrinsic wavelet parameters
......@@ -388,7 +388,7 @@ struct Data
double (*glitch_amplitude_prior) (double *, double *, double, double);
double (*intrinsic_likelihood)(int, int, int, struct Model*, struct Model*, struct Prior*, struct Chain*, struct Data*);
double (*extrinsic_likelihood)(struct Network*, double *, double **, struct Wavelet **, double **, double **, struct Data*, double, double);
double (*extrinsic_likelihood)(struct Network*, double *, double **, struct Wavelet **, double **, double *, struct Data*, double, double);
char runName[100];
char outputDirectory[1000];
......
......@@ -832,14 +832,14 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
}
double EvaluateExtrinsicConstantLogLikelihood(UNUSED struct Network *projection, double *params, UNUSED double **invSnf, UNUSED struct Wavelet **geo, UNUSED double **g, UNUSED double **geocbc, UNUSED struct Data *data, UNUSED double fmin, UNUSED double fmax)
double EvaluateExtrinsicConstantLogLikelihood(UNUSED struct Network *projection, double *params, UNUSED double **invSnf, UNUSED struct Wavelet **geo, UNUSED double **g, UNUSED double *amphase, UNUSED struct Data *data, UNUSED double fmin, UNUSED double fmax)
{
if(extrinsic_checkrange(params)) return -1.0e60;
else return 0.0;
}
double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, double *params, double **invSnf, struct Wavelet **geo, double **g, double **geocbc, struct Data *data, double fmin, double fmax)
double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, double *params, double **invSnf, struct Wavelet **geo, double **g, double *amphase, struct Data *data, double fmin, double fmax)
{
int i, n;
int NI,NP,N;
......@@ -886,7 +886,7 @@ double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, doubl
if(data->bayesCBCFlag)
{
// for(i=0; i<NI; i++) for(n=0; n<N; n++) c[i][n] = 0.0;
waveformProject(data, projection, params, c, geocbc, fmin, fmax);
projectCBCWaveform(amphase, N, NI, fmin, data->Tobs, extParams, c, projection->dtimes, projection->Fplus, projection->Fcross);
}
//Form up residual
......@@ -912,7 +912,7 @@ double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, doubl
return logL;
}
double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *params, double **invSnf, UNUSED double *Sn, struct Wavelet **geo, double **g, double **geocbc, struct Data *data, double fmin, double fmax)
double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *params, double **invSnf, UNUSED double *Sn, struct Wavelet **geo, double **g, double *amphase, struct Data *data, double fmin, double fmax)
{
int i, n;
int NI,N;
......@@ -949,7 +949,7 @@ double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *
//compute instrument response from cbc waveform at geocenter with extrinsic params
if(data->bayesCBCFlag)
{
waveformProject(data, projection, params, c, geocbc, data->fmin, data->fmax);
projectCBCWaveform(amphase, N, NI, fmin, data->Tobs, extParams, c, projection->dtimes, projection->Fplus, projection->Fcross);
}
//Form up residual
......@@ -993,10 +993,10 @@ double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *
computeProjectionCoeffs(data, projection, params, data->fmin, data->fmax);
waveformProject(data, projection, params, h, t, data->fmin, data->fmax);
//compute instrument response from cbc waveform at geocenter with extrinsic params
compute instrument response from cbc waveform at geocenter with extrinsic params
if(data->bayesCBCFlag)
{
waveformProject(data, projection, params, c, geocbc, data->fmin, data->fmax);
projectCBCWaveform(amphase, N, NI, fmin, data->Tobs, extParams, c, projection->dtimes, projection->Fplus, projection->Fcross);
}
//Form up residual
......
......@@ -37,8 +37,8 @@ void recompute_residual(struct Data *data, struct Model **model, struct Chain *c
double EvaluateConstantLogLikelihood (int typ, int ii, int det, UNUSED struct Model *mx, struct Model *my, struct Prior *prior, UNUSED struct Chain *chain, UNUSED struct Data *data);
double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx, struct Model *my, struct Prior *prior, struct Chain *chain, struct Data *data);
double EvaluateExtrinsicSearchLogLikelihood (struct Network *projection, double *params, double **invSnf, UNUSED double *Sn, struct Wavelet **geo, double **g, double **geocbc, struct Data *data, double fmin, double fmax);
double EvaluateExtrinsicConstantLogLikelihood (UNUSED struct Network *projection, double *params, UNUSED double **invSnf, UNUSED struct Wavelet **geo, UNUSED double **g, UNUSED double **geocbc, UNUSED struct Data *data, UNUSED double fmin, UNUSED double fmax);
double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, double *params, double **invSnf, struct Wavelet **geo, double **g, double **geocbc, struct Data *data, double fmin, double fmax);
double EvaluateExtrinsicSearchLogLikelihood (struct Network *projection, double *params, double **invSnf, UNUSED double *Sn, struct Wavelet **geo, double **g, double *amphase, struct Data *data, double fmin, double fmax);
double EvaluateExtrinsicConstantLogLikelihood (UNUSED struct Network *projection, double *params, UNUSED double **invSnf, UNUSED struct Wavelet **geo, UNUSED double **g, UNUSED double *amphase, UNUSED struct Data *data, UNUSED double fmin, UNUSED double fmax);
double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, double *params, double **invSnf, struct Wavelet **geo, double **g, double *amphase, struct Data *data, double fmin, double fmax);
void phase_blind_time_shift(double *corr, double *corrf, double *data1, double *data2, double *invpsd, int n);
......@@ -1358,8 +1358,8 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
{
model[chain->index[ic]]->cbcParams[i] = bayescbc->pallx[chain->index[ic]][i];
}
cbcGeoTemplate(bayescbc, chain->index[ic], data->N, model[chain->index[ic]]->cbcgeotemplate);
waveformProject(data, model[chain->index[ic]]->projection, model[chain->index[ic]]->extParams, model[chain->index[ic]]->cbctemplate, model[chain->index[ic]]->cbcgeotemplate, data->fmin, data->fmax);
// cbcGeoTemplate(bayescbc, chain->index[ic], data->N, model[chain->index[ic]]->cbcamphase);
// waveformProject(data, model[chain->index[ic]]->projection, model[chain->index[ic]]->extParams, model[chain->index[ic]]->cbctemplate, model[chain->index[ic]]->cbcamphase, data->fmin, data->fmax);
}
// Only evolve extrinsic params from now on
......@@ -1368,6 +1368,8 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
// Set MAP model to initial state
copy_cbc_model(model[chain->index[0]], model_MAP, N, NI, data->NX);
fprintf(stdout, " ======= End of CBC_burnin() ======\n");
for(ic=0; ic<chain->NC; ic++) printf("chain %i -> index %i\n", ic, chain->index[ic]);
}
/******************************************************************************/
......@@ -1421,7 +1423,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
//After so many iterations recompute the residuals and likelihood (prevent accumulation of roundoff error)
recompute_residual(data, model, chain);
// recompute_residual(data, model, chain);
if(burnFlag==0 && data->signalFlag && !data->fixIntrinsicFlag)
TFprop_signal(data, prior->range, tf, model[chain->index[0]]->projection);
......@@ -1702,7 +1704,7 @@ void PTMCMC(struct Model **model, struct Chain *chain, gsl_rng *seed, int NC)
chain->index[a] = oldb;
chain->index[b] = olda;
chain->A[a]=1;
printf("\nDEBUG: In PTMCMC() Swapping %i and %i (real index a: %i).\n\n", olda, oldb, a);
// printf("\nDEBUG: In PTMCMC() Swapping %i and %i (real index a: %i).\n\n", olda, oldb, a);
}
}
}
......@@ -2744,7 +2746,7 @@ void EvolveBayesCBCParameters(struct Data *data, struct Model **model, struct ba
if (bayescbc->debug==1) printf("before BayesCBC_MCMC chain %i, %i, model_x->logL = %f, logLnorm = %f, logLx = %f \n", chain->index[ic], ic, model_x->logL, model_x->logLnorm, bayescbc->logLx[icbc]);
// Do updates
BayesCBC_MCMC(M, N, data->Tobs, data->trigtime, chain->seed, istart, bayescbc, ic, model_x->cbcgeotemplate);
BayesCBC_MCMC(M, N, data->Tobs, data->trigtime, chain->seed, istart, bayescbc, ic, model_x->cbcamphase);
// Update CBC intrinsic parameters in chain model
for(i = 0; i < bayescbc->NX; i++)
......@@ -2753,7 +2755,9 @@ void EvolveBayesCBCParameters(struct Data *data, struct Model **model, struct ba
}
// Get CBC waveforms projected on the detectors
waveformProject(data, projection, model_x->extParams, model_x->cbctemplate, model_x->cbcgeotemplate, fmin, fmax);
projectCBCWaveform(model_x->cbcamphase, N, NI, fmin, data->Tobs, model_x->extParams, model_x->cbctemplate, projection->dtimes, projection->Fplus, projection->Fcross);
if (ic==0 && bayescbc->debug == 1) 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;
......@@ -2860,7 +2864,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
double fmin = data->fmin;
double fmax = data->fmax;
logLx = data->extrinsic_likelihood(model_x->projection, paramsx, model_x->invSnf, smodel, glitch, model_x->cbcgeotemplate, data, fmin, fmax);
logLx = data->extrinsic_likelihood(model_x->projection, paramsx, model_x->invSnf, smodel, glitch, model_x->cbcamphase, data, fmin, fmax);
//Compute Fisher Matrix for current sky location
struct FisherMatrix *fisher = model_x->fisher;
......@@ -2955,7 +2959,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if(test==0)
{
logLy = data->extrinsic_likelihood(model_x->projection, paramsy, model_x->invSnf, smodel, glitch, model_x->cbcgeotemplate, data, fmin, fmax);
logLy = data->extrinsic_likelihood(model_x->projection, paramsy, model_x->invSnf, smodel, glitch, model_x->cbcamphase, data, fmin, fmax);
//amplitude prior
if(data->amplitudePriorFlag)
......
......@@ -722,7 +722,7 @@ void initialize_model(struct Model *model, struct Data *data, struct Prior *prio
//instrument cbc template model
model->cbctemplate = double_matrix(NI-1,N-1);
// model->cbcgeotemplate = double_vector(N-1);
model->cbcgeotemplate = double_matrix(Npol, N-1);
model->cbcamphase = double_vector(N-1);
//store sum of 1/psd in 1/SnGeo
for(i=0; i<halfN; i++)
......@@ -1243,8 +1243,8 @@ void copy_cbc_model(struct Model *origin, struct Model *copy, int N, int NI, int
// Projected
for(ifo=0; ifo<NI; ifo++) copy->cbctemplate[ifo][i] = origin->cbctemplate[ifo][i];
// Geocenter
for(ip=0; ip<origin->Npol; ip++) copy->cbcgeotemplate[ip][i] = origin->cbcgeotemplate[ip][i];
// Geocenter Amp/Phase
copy->cbcamphase[i] = origin->cbcamphase[i];
}
// Copy CBC (intrinsic) model parameters
......@@ -1269,7 +1269,7 @@ void free_model(struct Model *model, struct Data *data, struct Prior *prior)
free_double_matrix(model->response,NI-1);
free_double_matrix(model->h,Npol);
free_double_matrix(model->cbctemplate,NI-1);
free_double_matrix(model->cbcgeotemplate, Npol);
free_double_vector(model->cbcamphase);
for(i=0; i<Npol; i++)
{
......
......@@ -349,8 +349,8 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
paramsP[i] += epsilon[i];
paramsM[i] -= epsilon[i];
logLP = EvaluateExtrinsicMarkovianLogLikelihood(projectionP, paramsP, invSnf, model->signal, glitch, model->cbcgeotemplate, data, data->fmin, data->fmax);
logLM = EvaluateExtrinsicMarkovianLogLikelihood(projectionM, paramsM, invSnf, model->signal, glitch, model->cbcgeotemplate, data, data->fmin, data->fmax);
logLP = EvaluateExtrinsicMarkovianLogLikelihood(projectionP, paramsP, invSnf, model->signal, glitch, model->cbcamphase, data, data->fmin, data->fmax);
logLM = EvaluateExtrinsicMarkovianLogLikelihood(projectionM, paramsM, invSnf, model->signal, glitch, model->cbcamphase, data, data->fmin, data->fmax);
epsilon[i] = 0.1/sqrt(-(logLM+logLP)/(epsilon[i]*epsilon[i]));
if(epsilon[i]!=epsilon[i])epsilon[i]=1.0e-5;
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment