Commit 80b5cad3 authored by Marcella Wijngaarden's avatar Marcella Wijngaarden

Cleanup mcmc_intrinsic, add data waves output, initialize who arrat

parent 10775237
......@@ -429,9 +429,9 @@ void wavemax(struct Net *net, int N, double **tp, RealVector *freq, double **par
void printwaveall(struct Net *net, int N, RealVector *freq, double *paramx, double **SN, double Tobs, double ttrig, int iter, struct bayesCBC *rundata)
{
double **twave, **tf, **thold, **twave_tmp;
double **twave, **tf, **thold, **twave_tmp, **dwave;
double **Dtime;
double **ww, **wf;
double **ww, **wf, **dw;
double fmin, fmax;
double fr, f, x, y, dt, fac, fs;
int i, j, id;
......@@ -449,6 +449,9 @@ void printwaveall(struct Net *net, int N, RealVector *freq, double *paramx, doub
tf = double_matrix(net->Nifo,N);
ww = double_matrix(net->Nifo,N);
wf = double_matrix(net->Nifo,N);
dwave = double_matrix(net->Nifo,N);
dw = double_matrix(net->Nifo,N); // Data whitened
fulltemplates(net, twave_tmp, freq, paramx, N, rundata);
......@@ -457,6 +460,24 @@ void printwaveall(struct Net *net, int N, RealVector *freq, double *paramx, doub
// printf("%d %f %f\n", iter, exp(paramx[0])/MSUN_SI, fs);
for (id = 0; id < net->Nifo; ++id)
{
sprintf(command, "wavesCBC/wavefcheck_%d_%d_%d_%d.dat", iter, (int)(Tobs), (int)ttrig, net->labels[id]);
out = fopen(command,"w");
for (i = 1; i < N/2; ++i)
{
f = (double)(i)/Tobs;
if(f > fs)
{
fprintf(out,"%e %.16e %.16e %.16e %.16e\n", f, twave_tmp[id][2*i]/(2.0*SN[id][i]), twave_tmp[id][2*i+1]/(2.0*SN[id][i]), rundata->D[id][2*i]/(2.0*SN[id][i]), rundata->D[id][2*i+1])/(2.0*SN[id][i]);
}
else
{
fprintf(out,"%e %.16e %.16e %.16e %.16e\n", f, 0.0, 0.0, 0.0, 0.0);
}
}
fclose(out);
}
// Convert indices for FFT
for (id = 0; id < net->Nifo; ++id)
......@@ -467,11 +488,15 @@ void printwaveall(struct Net *net, int N, RealVector *freq, double *paramx, doub
twave[id][N-i] = twave_tmp[id][2*i+1];
Dtime[id][i] = rundata->D[id][2*i];
Dtime[id][N-i] = rundata->D[id][2*i+1];
dwave[id][i] = rundata->D[id][2*i];
dwave[id][N-i] = rundata->D[id][2*i+1];
}
twave[id][0] = 0.0;
twave[id][N/2] = 0.0;
Dtime[id][0] = 0.0;
Dtime[id][N/2] = 0.0;
dwave[id][0] = 0.0;
dwave[id][N/2] = 0.0;
}
for (id = 0; id < net->Nifo; ++id)
......@@ -495,6 +520,8 @@ void printwaveall(struct Net *net, int N, RealVector *freq, double *paramx, doub
twave[id][i] *= x;
twave[id][N-i] *= x;
dwave[id][i] *= x;
dwave[id][N-i] *= x;
tf[id][i] = twave[id][N-i];
tf[id][N-i] = -twave[id][i];
......@@ -502,6 +529,9 @@ void printwaveall(struct Net *net, int N, RealVector *freq, double *paramx, doub
ww[id][N-i] = y*twave[id][N-i];
wf[id][i] = ww[id][N-i];
wf[id][N-i] = -ww[id][i];
dw[id][i] = y*dwave[id][i];
dw[id][N-i] = y*dwave[id][N-i];
}
}
......@@ -521,12 +551,15 @@ void printwaveall(struct Net *net, int N, RealVector *freq, double *paramx, doub
ww[id][N/2] = 0.0;
wf[id][0] = 0.0;
wf[id][N/2] = 0.0;
dw[id][0] = 0.0;
dw[id][N/2] = 0.0;
gsl_fft_halfcomplex_radix2_inverse(twave[id], 1, N);
gsl_fft_halfcomplex_radix2_inverse(tf[id], 1, N);
gsl_fft_halfcomplex_radix2_inverse(ww[id], 1, N);
gsl_fft_halfcomplex_radix2_inverse(wf[id], 1, N);
gsl_fft_halfcomplex_radix2_inverse(Dtime[id], 1, N);
gsl_fft_halfcomplex_radix2_inverse(dw[id], 1, N);
for (i = 0; i < N; ++i)
{
......@@ -535,6 +568,7 @@ void printwaveall(struct Net *net, int N, RealVector *freq, double *paramx, doub
tf[id][i] *= x;
ww[id][i] *= fac;
wf[id][i] *= fac;
dw[id][i] *= fac;
}
}
......@@ -579,6 +613,7 @@ void printwaveall(struct Net *net, int N, RealVector *freq, double *paramx, doub
// Why do we need this for segments > 4s to scale back the amplitude???
//fac = ceil(sqrt(Tobs/4.0));
fac = ceil(sqrt(1.0/4.0));
// double fac2 = sqrt((double)(N/2));
// printf(" amplitude factor: %f\n\n",fac);
for (id = 0; id < net->Nifo; ++id)
{
......@@ -589,6 +624,14 @@ void printwaveall(struct Net *net, int N, RealVector *freq, double *paramx, doub
fprintf(out,"%e %e %e\n", (double)(i)*dt-Tobs+2.0, ww[id][i]/fac, sqrt(ww[id][i]*ww[id][i]+wf[id][i]*wf[id][i]));
}
fclose(out);
sprintf(command, "wavesCBC/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);
}
......@@ -887,7 +930,7 @@ void burnin_bayescbc(int N, double Tobs, double ttrig, double **D, double **SN,
// Read in injection parameters from injection XML file.
// Note: assuming always starting from XML when injecting
if (start_at_injection > 0) {
if (true) {
printf("Reading XML injection parameters from file: %s\n", xml_file_path);
sprintf(filename, xml_file_path);
......@@ -951,54 +994,7 @@ void burnin_bayescbc(int N, double Tobs, double ttrig, double **D, double **SN,
sky[4] = 1.0;
sky[5] = 0.0;
sky[6] = 0.0;
// Else set initial parameters from ConstCBC.h (TODO: Make them cmd arguments)
} else if (start_at_injection == 100) {
// Set initial chain params to given values
params[0] = log(init_Mc*MSUN_SI); // log(Mc)
q = init_q; // q = m1/m2
mc = exp(params[0]);
m2 = mc*pow(1.0+q, 0.2)/pow(q,0.6);
m1 = q*m2;
mt = m1+m2;
printf("\n%f %f %f %f %f %f %f %f\n", m1/MSUN_SI, m2/MSUN_SI, q, mc/MSUN_SI, init_phi0, init_distance, rundata->init_tidal1, rundata->init_tidal2);
params[1] = log(mt); // log(Mt)
params[2] = init_chi1; // chi1
params[3] = init_chi2; // chi2
params[4] = init_phi0; // init_phi0;
params[5] = Tobs-2.0000;
params[6] = log(init_distance* PC_SI);
// Add tidal parameters
lambda_type_version = rundata->lambda_type_version;
if (NX > 7) {
params[7] = rundata->init_tidal1;
params[8] = rundata->init_tidal2;
}
// Parameters for computing scaling factor distance
printf(" Fixing sky parameters");
sky[0] = longitude; // alpha RA
sky[1] = sin(latitude); // sin(delta) sin(declination)
sky[2] = polarization; // polarization phase psi
sky[3] = cos(inclination); // ciota cos inclination angle orbital plane
// Initialize [4] scale [5] dphi [6] dt
// Scale changes the overall distance reference at geocenter
// dt shifts the geocenter time
// dphi shifts the geocenter phase
sky[4] = 1.0;
sky[5] = 0.0;
sky[6] = 0.0;
printf(" Fixing inital MCMC parameters");
printf("\n%f %f %f %f %f %f %f %f\n", init_Mc, init_q, init_chi1, init_chi2, init_phi0, init_distance, rundata->init_tidal1, rundata->init_tidal2);
}
}
printf("\n");
......@@ -1524,12 +1520,11 @@ void CBC_start(struct Net *net, int *mxc, FILE *chainS, double **paramx, double
for(k = 0; k < 4; k++)
{
skyx[j][k] = pallx[j][NX+k];
printf("sky[%i][%i] = %f\n", j, k, skyx[j][k]);
}
}
// Initial search to find signal
MCMC_intrinsic(net, 1, mxc, 10000, chainS, paramx, skyx, pallx, who, heat, history, global, freq, D, SN, N, Tobs, r, rundata);
MCMC_intrinsic(net, 1, mxc, 10000, chainS, paramx, skyx, pallx, heat, history, global, freq, D, SN, N, Tobs, r, rundata);
for(j = 0; j < NC; j++)
{
......@@ -1694,7 +1689,7 @@ void CBC_update(struct Net *net, int *mxc, FILE *chainI, FILE *chainE, double **
for(k = 0; k < 5; k++)
{
MCMC_intrinsic(net, 0, mxc, 200, chainI, paramx, skyx, pallx, who, heat, history, global, freq, D, SN, N, Tobs, r, rundata);
MCMC_intrinsic(net, 0, mxc, 200, chainI, paramx, skyx, pallx, heat, history, global, freq, D, SN, N, Tobs, r, rundata);
// Set up whitend signal arrays
......@@ -2490,43 +2485,31 @@ double globeden(double ***global, double *max, double *min, double Tobs, double
void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, double **paramx, double **skyx, double **pallx, int *who, double *heat, double ***history, double ***global, RealVector *freq, double **D, double **SN, int N, double Tobs, gsl_rng *r, struct bayesCBC *rundata)
void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, double **paramx, double **skyx, double **pallx, double *heat, double ***history, double ***global, RealVector *freq, double **D, double **SN, int N, double Tobs, gsl_rng *r, struct bayesCBC *rundata)
{
Lambda_type lambda_type_version;
int i, j, q, k, mc, flag;
int *m, *evolveparams;
double *logLx, logLy, logL, x;
double **paramy;
double *pref;
double *max, *min;
double fmin, fmax;
double *href;
double *pml;
double alpha, beta, H;
double Mchirp, Mtot, M1, M2, ciota;
double eta, dm, m1, m2, chieff, DL, f;
double a, b, c, d;
double sqh;
double Lmax, mch;
double logpx, logpy;
double pMcMx, pMcMy;
double dA, dTD, dP, z, y;
double leta;
int bin;
double eta, dm, m1, m2, chieff, DL, f, z;
int **av, **cv;
int typ, hold;
int scount, sacc, mcount;
double *Fscale;
double pDx, pDy, qx, qy;
double qx;
double maxL;
double *pmax;
double lambda1, lambda2, lambdaT, dLambdaT;
FILE *like;
double ***fish, ***evec;
double **ejump;
FILE *like;
// Local pointers to run settings
int NX = rundata->NX;
int NP = rundata->NP;
int NS = rundata->NS;
......@@ -2535,25 +2518,21 @@ void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, dou
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
int *who = rundata->who;
fmin = rundata->fmin;
fmax = rundata->fmax;
av = int_matrix(4,NC);
cv = int_matrix(4,NC);
m = int_vector(NC);
ejump = double_matrix(NC,NX);
fish = double_tensor(NC,NX,NX);
evec = double_tensor(NC,NX,NX);
pmax = double_vector(NX+3*net->Nifo);
paramy = double_matrix(NC,NX+3*net->Nifo);
logLx = double_vector(NC);
Fscale = double_vector(NC);
fmin = rundata->fmin;
fmax = rundata->fmax;
// Determine which parameters to evolve in MCMC
evolveparams = int_vector(NX);
for (i=0; i<NX; i++)
......@@ -2561,7 +2540,6 @@ void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, dou
evolveparams[i] = 1;
}
// reference antenna pattern scaling used to convert reference detector frame amplitude to a physical distance
if(lmax == 0)
{
......@@ -2611,7 +2589,7 @@ void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, dou
if (NX > 7) {
if (lambda_type_version == lambdaTilde) {
min[7] = rundata->lambdaTmin; // lambdaT
min[8] = rundata->dlambdaTmin; //-0.1; // dlambdaT
min[8] = rundata->dlambdaTmin; // dlambdaT
} else {
min[7] = rundata->lambda1min; // lambda1
min[8] = rundata->lambda2min; // lambda2
......@@ -2666,11 +2644,8 @@ void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, dou
for(k=0; k < NC; k++) m[k] = 0;
q = who[0];
printf("%f %f %f %f %f %f\n", paramx[q][4], paramx[q][5], paramx[q][6], paramx[q][7], paramx[q][8], paramx[q][9]);
q = who[0];
maxL = -1.0e20;
for(mc = 1; mc <= M; mc++)
{
......@@ -2721,29 +2696,6 @@ void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, dou
}
alpha = gsl_rng_uniform(r);
// if((NC > 1) && (alpha < 0.2)) // decide if we are doing a MCMC update of all the chains or a PT swap
// {
// // chain swap
// scount++;
// alpha = (double)(NC-1)*gsl_rng_uniform(r);
// j = (int)(alpha) + 1;
// beta = exp((logLx[who[j]]-logLx[who[j+1]])/heat[j+1] - (logLx[who[j]]-logLx[who[j+1]])/heat[j]);
// alpha = gsl_rng_uniform(r);
// if(beta > alpha)
// {
// hold = who[j];
// who[j] = who[j+1];
// who[j+1] = hold;
// sacc++;
// }
// }
// else // MCMC update
// {
mcount++;
for(j = 0; j < NC; j++)
......@@ -2771,9 +2723,7 @@ void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, dou
m[k]++;
}
}
// }
// update maxL parameters
if(logLx[who[0]] > maxL)
{
......@@ -2797,12 +2747,10 @@ void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, dou
// pall [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
pmap(net, pallx[q], paramx[q], skyx[q], NX, rundata->gmst);
DL = Fscale[q]*exp(paramx[q][6])/(1.0e6*PC_SI); // includes conversion to Mpc
z = z_DL(DL);
// printf("%e %e\n", DL, exp(pallx[q][6])/(1.0e6*PC_SI));
// counter, log likelihood, chirp mass, total mass, effective spin, geocenter GW phase, geocenter arrival time, distance, RA , sine of DEC,
// polarization angle, cos inclination
......@@ -2850,41 +2798,15 @@ void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, dou
printf("%d %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n\n", k, logLx[q], Mchirp, Mtot, f, paramx[q][5], paramx[q][6], paramx[q][NX], paramx[q][NX+1], paramx[q][NX+2], \
(double)(sacc)/(double)(scount), (double)(av[1][1])/(double)(cv[1][1]), (double)(av[2][1])/(double)(cv[2][1]), (double)(av[3][1])/(double)(cv[3][1]), (double)(av[4][1])/(double)(cv[4][1]));
//printf("%d %f %f %f %f %f\n", mc/100, paramx[q][5], paramx[q][6], paramx[q][7], paramx[q][8], paramx[q][9]);
for(k=0; k < NC; k++)
{
q = who[k];
Mchirp = exp(paramx[q][0])/MSUN_SI;
printf("%i %f %f %f %f %d \n", k, logLx[q], Mchirp, Mtot, DL, q);
}
printf("\n");
/* printf("%d ", mc/100);
for(k=1; k <= NC; k++)
{
q = who[k];
x = log_likelihood(net, D, paramx[q], freq, SN, N, Tobs);
printf("%f ", x);
}
printf("\n"); */
/* printf("%d ", mc/100);
for(k=1; k <= NC; k++)
{
q = who[k];
printf("%f ", paramx[q][5]);
}
printf("\n"); */
}
//printf("%g\n",maxL);
printf("\n");
}
}
......@@ -3076,14 +2998,11 @@ void updatei(int k, struct Net *net, int lmax, double *logLx, double **paramx,
if(flag == 1) logLy = 0.0;
// DLx = Fscale[q]*exp(paramx[q][6]);
// DLy = Fscale[q]*exp(paramy[q][6]);
// since the Fscale and sky parameters are held fixed in the extrinsic MCMC can
// use the expresion below for the DL^2 prior
// variable in MCMC is x=logD, so p(x) = dD/dx p(D) = D p(D) = D^3
pDx = 3.0*paramx[q][6]; // unet->Nifoform in volume prior
......@@ -3096,11 +3015,8 @@ void updatei(int k, struct Net *net, int lmax, double *logLx, double **paramx,
qy = qx = 0.0;
}
H = (logLy-logLx[q])/heat[k] + pMcMy + pDy - qy - pDx - pMcMx + qx;
alpha = log(gsl_rng_uniform(r));
if(H > alpha && flag == 0)
......@@ -3705,7 +3621,7 @@ void MCMC_all_wrapper(int M, int N, double Tobs, double ttrig, gsl_rng *r, int i
void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int istart, struct bayesCBC *rundata, int ic, double **template)
{
int i, j, q, k, mc, flag;
int m, *evolveparams, *who;
int m, *evolveparams;
double logLy, logL, x;
double *paramy, *paramx;
double Mchirp, Mtot, M1, M2, ciota, tc;
......@@ -3744,12 +3660,6 @@ void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int ista
av = int_vector(4);
cv = int_vector(4);
// This is for compatibility with mcmc_all loop updates, but
// parallel tempering/chain swaps is not actually done here
who = int_vector(ic);
who[ic] = ic; // The current chain
rhox = double_vector(net->Nifo);
ejump = double_vector(NP);
......@@ -3764,7 +3674,6 @@ void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int ista
paramx = rundata->pallx[ic];
// Determine which parameters to evolve in MCMC
// Currently hardcoded here for testing
evolveparams = int_vector(NP);
for (i=0; i<NP; i++)
{
......@@ -3819,9 +3728,9 @@ void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int ista
// if(mc > 0 && mc%500==0) printwaveall(net, N, freq, paramx, rundata->SN, Tobs, ttrig, mc/500, rundata);
// if((rundata->mxc[0]+mc)%5000==0 )
if(mc%5000==0 )
if(mc%200==0 )
{
if (ic == 0) printf("\nupdating fisher matric at i: %i\n", rundata->mxc[0]+mc);
if (ic == 0) printf("\nupdating fisher matric at mc: %i\n", mc);
// update the Fisher matrices
Fisher_All(net, fish, paramx, freq, rundata->SN, N, Tobs, rundata);
FisherEvec(fish, ejump, evec, NP);
......@@ -3853,7 +3762,7 @@ void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int ista
}
if(ic==0)
if(rundata->who[ic]==0)
{
Mchirp = exp(paramx[0])/MSUN_SI;
Mtot = exp(paramx[1])/MSUN_SI;
......@@ -3878,7 +3787,7 @@ void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int ista
}
// Print waveform
printwaveall(rundata->net, N, rundata->freq, paramx, rundata->SN, Tobs, ttrig, rundata->mxc[0], rundata);
printwaveall(rundata->net, N, rundata->freq, rundata->pallx[ic], rundata->SN, Tobs, ttrig, rundata->mxc[0], rundata);
}
......
......@@ -189,7 +189,7 @@ void CBC_start(struct Net *net, int *mxc, FILE *chainI, double **paramx, double
void CBC_update(struct Net *net, int *mxc, FILE *chainI, FILE *chainE, double **paramx, double **skyx, double **pallx, int *who, double *heat, double ***history, double ***global, RealVector *freq, double **D, double **SN, int N, double Tobs, gsl_rng *r, struct bayesCBC *rundata);
void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, double **paramx, double **skyx, double **pallx, int *who, double *heat, double ***history, double ***global, RealVector *freq, double **D, double **SN, int N, double Tobs, gsl_rng *r, struct bayesCBC *rundata);
void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, double **paramx, double **skyx, double **pallx, double *heat, double ***history, double ***global, RealVector *freq, double **D, double **SN, int N, double Tobs, gsl_rng *r, struct bayesCBC *rundata);
void skymcmc(struct Net *net, int MCX, int *mxc, FILE *chain, double **paramx, double **skyx, double **pallx, int *who, double *heat, double dtx, int nt, double *DD, double **WW, double ***DHc, double ***DHs, double ***HH, double Tobs, gsl_rng * r, struct bayesCBC *rundata);
......
......@@ -1260,7 +1260,8 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
if(LALInferenceGetProcParamVal(data->commandLine, "--inj"))
{
sprintf(filename, LALInferenceGetProcParamVal(data->commandLine,"--inj")->value);
start_at_injection = 1;
// If injecting tidal parameters, read values from command line and setup version + sampler type
if(LALInferenceGetProcParamVal(data->commandLine, "--inj-lambda1")||LALInferenceGetProcParamVal(data->commandLine, "--inj-lambdaT"))
{
......@@ -2733,6 +2734,7 @@ void EvolveBayesCBCParameters(struct Data *data, struct Model **model, struct ba
// Set heat
bayescbc->heat[icbc] = chain->temperature[icbc];
bayescbc->mxc[0] = chain->mc;
bayescbc->who = chain->index;
// Do updates
BayesCBC_MCMC(M, N, data->Tobs, data->trigtime, chain->seed, istart, bayescbc, icbc, model_x->cbcgeotemplate);
......
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