Commit 9c91308a authored by Marcella Wijngaarden's avatar Marcella Wijngaarden

Update CBC model in same order as others + add burnin templates

In the EvolveParameters loop the CBC parameters of the same model
as updated in the intrinsic/extrinsic updates for that iteration
is updated. (i.e, model chain->index[ic] instead of model ic).
However, in doing this PTMCMC doesn't seem to like thing anymore.
Only running this with a single chain. And I am looking into
getting chain swaps working correctly again.

This commit also adds the CBC templates from the burnin so
they can be used in he first iteration of evolve extrinsic/
intrinsic parametes.
parent 87ac426b
Pipeline #159193 failed with stages
in 2 minutes and 2 seconds
......@@ -1204,7 +1204,7 @@ void CBC_start(struct Net *net, int *mxc, FILE *chainS, double **paramx, double
}
// Initial search to find signal
MCMC_intrinsic(net, 1, mxc, 10000, chainS, paramx, skyx, pallx, heat, history, global, freq, D, SN, N, Tobs, r, rundata);
MCMC_intrinsic(net, 1, mxc, 9500, chainS, paramx, skyx, pallx, heat, history, global, freq, D, SN, N, Tobs, r, rundata);
for(j = 0; j < NC; j++)
{
......@@ -3300,6 +3300,9 @@ void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int ista
RealVector *freq = rundata->freq;
double *max = rundata->max;
double *min = rundata->min;
int *who = rundata->who;
int icbc = who[ic];
av = int_vector(4);
cv = int_vector(4);
......@@ -3314,7 +3317,7 @@ void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int ista
lambda_type_version=rundata->lambda_type_version;
paramx = rundata->pallx[ic];
paramx = rundata->pallx[icbc];
// Determine which parameters to evolve in MCMC
evolveparams = int_vector(NP);
......@@ -3345,9 +3348,9 @@ void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int ista
}
// Initialize the cbc likelihood
rundata->logLx[ic] = log_likelihood_full(net, rundata->D, paramx, freq, rundata->SN, rhox, N, Tobs, rundata);
if(rundata->constantLogLFlag == 0) rundata->logLx[ic] = log_likelihood_full(net, rundata->D, paramx, freq, rundata->SN, rhox, N, Tobs, rundata);
// printf("logLx at start %i: %f \n", ic, rundata->logLx[ic]);
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]);
Fisher_All(net, fish, paramx, freq, rundata->SN, N, Tobs, rundata);
......@@ -3382,7 +3385,7 @@ void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int ista
}
// Actual parameter updates
update_chain(ic, net, rundata->logLx, rhox, paramx, paramy, min, max, rundata->heat[ic], rundata->history[ic], global, freq, rundata->D, rundata->SN, ejump, evec, N, Tobs, cv, av, r, evolveparams, rundata);
update_chain(icbc, net, rundata->logLx, rhox, paramx, paramy, min, max, rundata->heat[ic], rundata->history[ic], global, freq, rundata->D, rundata->SN, ejump, evec, N, Tobs, cv, av, r, evolveparams, rundata);
// add to the history file
if(mc%20 == 0)
......@@ -3394,9 +3397,9 @@ void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int ista
}
// update MAP parameters
if(rundata->logLx[ic] > mapL)
if(rundata->logLx[icbc] > mapL)
{
mapL = rundata->logLx[ic];
mapL = rundata->logLx[icbc];
for(i = 0; i < NP; i++) pmax[i] = paramx[i];
}
}
......@@ -3417,11 +3420,11 @@ void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int ista
lambda2 = paramx[8];
}
// Also print tidal parameters
printf("%d %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n", mc/100, rundata->logLx[ic], Mchirp, Mtot, paramx[4], paramx[5], paramx[6], paramx[7], paramx[8], paramx[NX], paramx[NX+1], paramx[NX+2], \
printf("%d %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n", mc/100, rundata->logLx[icbc], Mchirp, Mtot, paramx[4], paramx[5], paramx[6], paramx[7], paramx[8], paramx[NX], paramx[NX+1], paramx[NX+2], \
(double)(av[1])/(double)(cv[1]), (double)(av[2])/(double)(cv[2]), (double)(av[3])/(double)(cv[3]), (double)(av[4])/(double)(cv[4]),
lambda1, lambda2, lambdaT, dLambdaT);
} else {
printf("%d %f %f %f %f %f %f %f %f %f %f \n", mc/100, rundata->logLx[ic], Mchirp, Mtot, paramx[4], paramx[5], paramx[6], paramx[NX], paramx[NX+1], paramx[NX+2], paramx[NX+3]);
printf("%d %f %f %f %f %f %f %f %f %f %f \n", mc/100, rundata->logLx[icbc], Mchirp, Mtot, paramx[4], paramx[5], paramx[6], paramx[NX], paramx[NX+1], paramx[NX+2], paramx[NX+3]);
for (i = 0; i < NP; i++) printf("%f \t", paramx[i]);
printf("\n");
}
......@@ -3430,7 +3433,7 @@ void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int ista
if(ic==0 && rundata->verbose==1)
{
// Save waveforms to file
printwaveall(rundata->net, N, rundata->freq, rundata->pallx[ic], rundata->SN, Tobs, ttrig, rundata->mxc[0], rundata);
printwaveall(rundata->net, N, rundata->freq, rundata->pallx[icbc], rundata->SN, Tobs, ttrig, rundata->mxc[0], rundata);
}
// Save MAP likelihood. TODO: Do we want to have this?
......@@ -3438,9 +3441,9 @@ void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int ista
// Store geotemplate
// fulltemplates(rundata->net, template, rundata->freq, rundata->pallx[ic], N, rundata);
geotemplate(geocbc, rundata->freq, rundata->pallx[ic], N, rundata);
geotemplate(geocbc, rundata->freq, rundata->pallx[icbc], N, rundata);
ciota = rundata->pallx[ic][NX+3];
ciota = rundata->pallx[icbc][NX+3];
double epsilon =2*ciota/(1+ciota*ciota);
int re, im;
......@@ -4469,6 +4472,36 @@ void templates(struct Net *net, double **hwave, RealVector *freq, double *params
}
void cbcGeoTemplate(struct bayesCBC *rundata, int ic, int N, double **template)
{
int re, im, i;
int NX = rundata->NX;
double *geocbc;
geocbc = double_vector(N);
geotemplate(geocbc, rundata->freq, rundata->pallx[ic], N, rundata);
double ciota = rundata->pallx[ic][NX+3];
double epsilon =2*ciota/(1+ciota*ciota);
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];
}
free_double_vector(geocbc);
}
void geotemplate(double *gwave, RealVector *freq, double *params, int N, struct bayesCBC *rundata)
{
......
......@@ -116,6 +116,7 @@ void MCMC_all(struct Net *net, int *mxc, int M, FILE *chain, double **paramx, in
void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int istart, struct bayesCBC *rundata, int ic, double **template);
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);
double z_DL(double DL);
double DL_fit(double z);
......
......@@ -195,8 +195,8 @@ void recompute_residual(struct Data *data, struct Model **model, struct Chain *c
for(i=0; i<data->N; i++)
{
data->r[ifo][i] = data->s[ifo][i];
if(data->glitchFlag) data->r[ifo][i] -= model_x->glitch[ifo]->templates[i];
if(data->signalFlag) data->r[ifo][i] -= model_x->response[ifo][i];
if(data->glitchFlag && !data->cbcOnlyFlag) data->r[ifo][i] -= model_x->glitch[ifo]->templates[i];
if(data->signalFlag && !data->cbcOnlyFlag) data->r[ifo][i] -= model_x->response[ifo][i];
if(data->bayesCBCFlag) data->r[ifo][i] -= model_x->cbctemplate[ifo][i];
}
}
......@@ -880,11 +880,12 @@ double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, doubl
combinePolarizations(data, geo, t, params, data->Npol);
computeProjectionCoeffs(data, projection, params, fmin, fmax);
waveformProject(data, projection, params, h, t, fmin, fmax);
//compute instrument response from cbc waveform at geocenter with extrinsic params
if(data->bayesCBCFlag)
{
for(i=0; i<NI; i++) for(n=0; n<N; n++) c[i][n] = 0.0;
// 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);
}
......@@ -894,13 +895,14 @@ double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, doubl
for(n=0; n<N; n++)
{
r[ifo][n] = d[ifo][n];
if(data->signalFlag) r[ifo][n] -= h[ifo][n];
if(data->glitchFlag) r[ifo][n] -= g[ifo][n];
if(data->signalFlag && !data->cbcOnlyFlag) r[ifo][n] -= h[ifo][n];
if(data->glitchFlag && !data->cbcOnlyFlag) r[ifo][n] -= g[ifo][n];
if(data->bayesCBCFlag) r[ifo][n] -= c[ifo][n];
}
logL += loglike(imin, imax, r[ifo], invSnf[ifo]);
}
// if(data->bayesCBCFlag) print_projected_cbc_waveform(c[ifo][n]);
free_double_matrix(h,NI-1);
free_double_matrix(r,NI-1);
......@@ -956,8 +958,8 @@ double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *
for(n=0; n<N; n++)
{
r[ifo][n] = d[ifo][n];
if(data->signalFlag) r[ifo][n] -= h[ifo][n];
if(data->glitchFlag) r[ifo][n] -= g[ifo][n];
if(data->signalFlag && !data->cbcOnlyFlag) r[ifo][n] -= h[ifo][n];
if(data->glitchFlag && !data->cbcOnlyFlag) r[ifo][n] -= g[ifo][n];
if(data->bayesCBCFlag) r[ifo][n] -= c[ifo][n];
}
}
......@@ -1003,8 +1005,8 @@ double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *
for(n=0; n<N; n++)
{
r[ifo][n] = d[ifo][n];
if(data->signalFlag) r[ifo][n] -= h[ifo][n];
if(data->glitchFlag) r[ifo][n] -= g[ifo][n];
if(data->signalFlag && !data->cbcOnlyFlag) r[ifo][n] -= h[ifo][n];
if(data->glitchFlag && !data->cbcOnlyFlag) r[ifo][n] -= g[ifo][n];
if(data->bayesCBCFlag) r[ifo][n] -= c[ifo][n];
}
}
......
......@@ -266,12 +266,13 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
{
// TEMPORARY ADDITION MARCELLA TO KEEP ALL EXTRINISC FIXED FOR GW150914
if(data->fixExtrinsicFlag == 1)
// TEMPORARY ADDITION MARCELLA TO SET INITAL EXTRINISC FIXED FOR GW150914
// if(data->fixExtrinsicFlag == 1)
if(true)
{
model[ic]->extParams[2] = 1.5000;
model[ic]->extParams[3] = -0.975609756097561;
printf("IMPORTANT: FIXED EXTRINISC PARAMS TO CHOSEN VALUES\n");
printf("IMPORTANT: FIXED INITIAL EXTRINISC PARAMS TO CHOSEN VALUES\n");
}
else
{
......@@ -344,8 +345,8 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
for(i=0; i<data->N; i++)
{
data->r[ifo][i] = data->s[ifo][i];
if(data->signalFlag) data->r[ifo][i] -= model[ic]->response[ifo][i];
if(data->glitchFlag) data->r[ifo][i] -= glitch[ifo]->templates[i];
if(data->signalFlag && !data->cbcOnlyFlag) data->r[ifo][i] -= model[ic]->response[ifo][i];
if(data->glitchFlag && !data->cbcOnlyFlag) data->r[ifo][i] -= glitch[ifo]->templates[i];
if(data->bayesCBCFlag) data->r[ifo][i] -= model[ic]->cbctemplate[ifo][i];
}
if(!data->constantLogLFlag)
......@@ -1247,7 +1248,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
// Set BayesCBC rundata (this replaces the header file settings of QuickCBC)
struct bayesCBC *bayescbc = malloc(sizeof(struct bayesCBC));
if (data->runPhase && data->signalFlag && data->bayesCBCFlag)
if (data->runPhase && (data->signalFlag||data->glitchFlag||data->fullModelFlag) && data->bayesCBCFlag)
{
double x, y;
printf("Setting up BayesCBC in RJMCMC\n");
......@@ -1310,7 +1311,6 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
{
// bayescbc->rvec[i] = gsl_rng_alloc(T);
bayescbc->rvec[chain->index[i]] = chain->rvec[chain->index[i]];
model[chain->index[i]]->cbcParams = double_vector(bayescbc->NX);
}
// Setup antenna patterns and time delays for current extrinsic params
......@@ -1327,6 +1327,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
// Pass extrinsic params from BW to CBC
for(ic=0; ic<chain->NC; ic++)
{
model[chain->index[i]]->cbcParams = double_vector(bayescbc->NX);
for(i = bayescbc->NX; i < bayescbc->NP; i++)
{
// Convert Bayeswaves' ellipticity to BayesCBC's cos(inclination)
......@@ -1356,9 +1357,11 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
for(i = 0; i < bayescbc->NX; i++)
{
model[chain->index[ic]]->cbcParams[i] = bayescbc->pallx[chain->index[ic]][i];
// waveformProject(data, projection, model[chain->index[ic]]->extParams, model[chain->index[ic]]->cbctemplate, model[chain->index[ic]]->cbcgeotemplate, fmin, fmax);
}
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);
}
// Only evolve extrinsic params from now on
bayescbc->intrinsic_only = 1;
......@@ -1406,7 +1409,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
if(!data->fixIntrinsicFlag) EvolveIntrinsicParameters(data, prior, model, chain, tf, chain->seed, ic);
//This function executes [cycle] extrinsic parameter updates, common to all geocenter wavelets.
if(model[chain->index[ic]]->signal[0]->size>0 && !data->fixExtrinsicFlag) EvolveExtrinsicParameters(data, prior, model, chain, chain->seed, ic);
if((model[chain->index[ic]]->signal[0]->size>0 || data->bayesCBCFlag) && !data->fixExtrinsicFlag && chain->mc > 1) EvolveExtrinsicParameters(data, prior, model, chain, chain->seed, ic);
//update PSD with BayesLine
if(data->bayesLineFlag) EvolveBayesLineParameters(data, model, bayesline, chain, prior, ic);
......@@ -1503,16 +1506,15 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
logPost = model[chain->index[0]]->logL;
if(data->glitchFlag) for(ifo=0; ifo<NI; ifo++) logPost += model[chain->index[0]]->glitch[ifo]->logp;
if(data->signalFlag) logPost += model[chain->index[0]]->signal[0]->logp;
// TODO: Update LogP with bayesCBC MAP
// if(data->bayesCBCFlag) logPost += model[chain->index[0]]->signal[0]->logp;
// TODO: Update LogP with bayesCBC
if(logPost > logPostMap)
{
logPostMap = logPost;
copy_psd_model(model[chain->index[0]], model_MAP, N, NI);
copy_ext_model(model[chain->index[0]], model_MAP, N, NI);
copy_int_model(model[chain->index[0]], model_MAP, N, NI,-1);
copy_cbc_model(model[chain->index[0]], model_MAP, N, NI, data->NX);
for(ifo=0; ifo<NI; ifo++) copy_int_model(model[chain->index[0]], model_MAP, N, NI,ifo);
copy_cbc_model(model[chain->index[0]], model_MAP, N, NI, data->NX);
// LaplaceApproximation(data, model_MAP, chain, prior, &logZ);
}
}
......@@ -1700,7 +1702,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);
}
}
}
......@@ -1728,8 +1730,8 @@ void EvolveBayesLineParameters(struct Data *data, struct Model **model, struct B
for(i=0; i<N; i++)
{
data->r[ifo][i] = data->s[ifo][i];
if(data->signalFlag) data->r[ifo][i] -= model_x->response[ifo][i];
if(data->glitchFlag) data->r[ifo][i] -= model_x->glitch[ifo]->templates[i];
if(data->signalFlag && !data->cbcOnlyFlag) data->r[ifo][i] -= model_x->response[ifo][i];
if(data->glitchFlag && !data->cbcOnlyFlag) data->r[ifo][i] -= model_x->glitch[ifo]->templates[i];
if(data->bayesCBCFlag) data->r[ifo][i] -= model_x->cbctemplate[ifo][i];
}
......@@ -1749,8 +1751,8 @@ void EvolveBayesLineParameters(struct Data *data, struct Model **model, struct B
for(i=0; i< N; i++)
{
data->r[ifo][i] = data->s[ifo][i];
if(data->signalFlag) data->r[ifo][i] -= model_x->response[ifo][i];
if(data->glitchFlag) data->r[ifo][i] -= model_x->glitch[ifo]->templates[i];
if(data->signalFlag && !data->cbcOnlyFlag) data->r[ifo][i] -= model_x->response[ifo][i];
if(data->glitchFlag && !data->cbcOnlyFlag) data->r[ifo][i] -= model_x->glitch[ifo]->templates[i];
if(data->bayesCBCFlag) data->r[ifo][i] -= model_x->cbctemplate[ifo][i];
}
......@@ -2654,11 +2656,11 @@ void EvolveBayesCBCParameters(struct Data *data, struct Model **model, struct ba
int NI = data->NI;
double fmin = data->fmin;
double fmax = data->fmax;
int icbc = ic;
int icbc = chain->index[ic];
int istart = 0;
// Pointers to structures
struct Model *model_x = model[icbc];
struct Model *model_x = model[chain->index[ic]];
struct Network *projection = model_x->projection;
imin = fmin*data->Tobs;
......@@ -2672,15 +2674,16 @@ void EvolveBayesCBCParameters(struct Data *data, struct Model **model, struct ba
for(i=0; i<N; i++)
{
data->r[ifo][i] = data->s[ifo][i];
if(data->signalFlag) data->r[ifo][i] -= model_x->response[ifo][i];
if(data->glitchFlag) data->r[ifo][i] -= model_x->glitch[ifo]->templates[i];
if(data->signalFlag && !data->cbcOnlyFlag) data->r[ifo][i] -= model_x->response[ifo][i];
if(data->glitchFlag && !data->cbcOnlyFlag) data->r[ifo][i] -= model_x->glitch[ifo]->templates[i];
bayescbc->D[ifo][i] = data->r[ifo][i];
}
}
// Use latest estimates of PSD if using BayesLine
if (data->bayesLineFlag)
// if (data->bayesLineFlag)
if (true)
{
for(ifo=0; ifo<NI; ifo++)
{
......@@ -2734,12 +2737,14 @@ void EvolveBayesCBCParameters(struct Data *data, struct Model **model, struct ba
}
// Set heat
bayescbc->heat[icbc] = chain->temperature[icbc];
bayescbc->heat[ic] = chain->temperature[ic];
bayescbc->mxc[0] = chain->mc;
bayescbc->who = chain->index;
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, icbc, model_x->cbcgeotemplate);
BayesCBC_MCMC(M, N, data->Tobs, data->trigtime, chain->seed, istart, bayescbc, ic, model_x->cbcgeotemplate);
// Update CBC intrinsic parameters in chain model
for(i = 0; i < bayescbc->NX; i++)
......@@ -2747,6 +2752,7 @@ void EvolveBayesCBCParameters(struct Data *data, struct Model **model, struct ba
model_x->cbcParams[i] = bayescbc->pallx[icbc][i];
}
// Get CBC waveforms projected on the detectors
waveformProject(data, projection, model_x->extParams, model_x->cbctemplate, model_x->cbcgeotemplate, fmin, fmax);
// Recompute likelihoods of current chain
......@@ -2761,8 +2767,8 @@ void EvolveBayesCBCParameters(struct Data *data, struct Model **model, struct ba
for(i=0; i< N; i++)
{
data->r[ifo][i] = data->s[ifo][i];
if(data->signalFlag) data->r[ifo][i] -= model_x->response[ifo][i];
if(data->glitchFlag) data->r[ifo][i] -= model_x->glitch[ifo]->templates[i];
if(data->signalFlag && !data->cbcOnlyFlag) data->r[ifo][i] -= model_x->response[ifo][i];
if(data->glitchFlag && !data->cbcOnlyFlag) data->r[ifo][i] -= model_x->glitch[ifo]->templates[i];
data->r[ifo][i] -= model_x->cbctemplate[ifo][i];
}
......@@ -2862,13 +2868,13 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
double draw;
logpx=0.0;
if(data->amplitudePriorFlag)
if(data->amplitudePriorFlag )
{
//model_x->projection has be updated by data->extrinsic_likelihood, so it stores F+ and Fx for params_y.
Shf_Geocenter_full(data, model_x->projection, model_x->Snf, model_x->SnGeo, model_x->extParams);
for(i=1; i<ienddim; i++)
{
logpx += (data->signal_amplitude_prior(smodel[0]->intParams[smodel[0]->index[i]],SnGeox, data->Tobs, prior->sSNRpeak));
if (!data->cbcOnlyFlag) logpx += (data->signal_amplitude_prior(smodel[0]->intParams[smodel[0]->index[i]],SnGeox, data->Tobs, prior->sSNRpeak));
}
}
//logpx += dim*log(paramsx[5]);
......@@ -2941,7 +2947,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
for(i=1; i<ienddim; i++)
{
intParams[i][3] = smodel[0]->intParams[smodel[0]->index[i]][3]*paramsy[5];//amplitude
test += checkrange(intParams[i],prior->range, smodel[0]->dimension);
if (!data->cbcOnlyFlag) test += checkrange(intParams[i],prior->range, smodel[0]->dimension);
}
//check extrinsic parameters
......@@ -2956,7 +2962,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
{
//model_x->projection has be updated by data->extrinsic_likelihood, so it stores F+ and Fx for params_y.
Shf_Geocenter_full(data, model_x->projection, model_x->Snf, SnGeoy, paramsy);
for(i=1; i<ienddim; i++) logpy += (data->signal_amplitude_prior(intParams[i],SnGeoy, data->Tobs, prior->sSNRpeak));
if (!data->cbcOnlyFlag) for(i=1; i<ienddim; i++) logpy += (data->signal_amplitude_prior(intParams[i],SnGeoy, data->Tobs, prior->sSNRpeak));
}
logH = (logLy - logLx)*chain->beta + logpy - logpx + logJ;
......
......@@ -1247,7 +1247,7 @@ void copy_cbc_model(struct Model *origin, struct Model *copy, int N, int NI, int
for(ip=0; ip<origin->Npol; ip++) copy->cbcgeotemplate[ip][i] = origin->cbcgeotemplate[ip][i];
}
// Copy CBC model parameters
// Copy CBC (intrinsic) model parameters
for(i = 0; i < NX; i++) copy->cbcParams[i] = origin->cbcParams[i];
}
......
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