Commit 97c15098 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

review-related tweaks to proposals, temp ladder


git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@612 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent a3a05add
......@@ -223,7 +223,7 @@ int main(int argc, char *argv[])
//give chain buffer in case we need more
int NC=chain->NC;
chain->NC+=5;
chain->NC+=10;
allocate_chain(chain,NI);
for(ic=0; ic<chain->NC; ic++) chain->index[ic]=ic;
......@@ -556,6 +556,7 @@ int main(int argc, char *argv[])
//if(sqrt(rSNR)>0.0) resize_model(data, chain, prior, model, bayesline, psd, chain->NC+5);
if(sqrt(rSNR)>100.0) chain->NC+=5;
if(sqrt(rSNR)>200.0) chain->NC+=5;
/******************************************************************************/
/* */
......@@ -660,8 +661,15 @@ int main(int argc, char *argv[])
printf(" Increasing number of chains to %i\n",chain->NC+5);
chain->NC+=5;
}
if(data->cleanModelFlag==0 && SNRinj[NI]>200.)
{
printf("MDC/XML injection with SNR>100\n");
printf(" Cleaning phase disabled\n");
printf(" Increasing number of chains to %i\n",chain->NC+5);
chain->NC+=5;
}
free(SNRinj);
}
......@@ -679,99 +687,97 @@ int main(int argc, char *argv[])
sprintf(filename,"%sevidence.dat",data->runName);
evidence = fopen(filename,"w");
/******************************************************************************/
/* */
/* Gaussian Noise Characterization phase */
/* */
/******************************************************************************/
if(data->noiseModelFlag)
{
/******************************************************************************/
/* */
/* Signal Characterization phase */
/* */
/******************************************************************************/
data->runPhase = 1;
reset_priors(data, prior);
if(data->bayesLineFlag) data->psdFitFlag = 0;
if(data->signalModelFlag)
{
//Initialize proposal ratios for GW search phase
chain->rjmcmcRate = 0.0;
data->runPhase = 1;
reset_priors(data, prior);
//set flags for no wavelets in model
data->cleanFlag = 0;
data->glitchFlag = 0;
data->signalFlag = 0;
//Initialize proposal ratios for GW search phase
chain->modelRate = 1; //always do signal model updates
chain->rjmcmcRate = 0; //model fixed to one wavelet for extrinsic search
initialize_chain(chain,1);
data->cleanFlag = 0;
data->glitchFlag = 0;
data->signalFlag = 1;
reset_likelihood(data);
reset_priors(data, prior);
reset_model(data, chain, prior, model);
initialize_chain(chain,1);
printf("characterizing Gaussian noise model...\n\n");
reset_likelihood(data);
reset_priors(data, prior);
reset_model(data, chain, prior, model);
RJMCMC(data, model, bayesline, chain, prior, &data->logZnoise, &data->varZnoise);
RJMCMC(data, model, bayesline, chain, prior, &data->logZsignal, &data->varZsignal);
//shut off noise model for checkpointing
data->noiseModelFlag = 0;
}
fprintf(evidence,"noise %lg %lg\n",data->logZnoise,data->varZnoise);
fflush(evidence);
//shut off signal model for checkpointing
data->signalModelFlag = 0;
}
fprintf(evidence,"signal %lg %lg\n",data->logZsignal,data->varZsignal);
fflush(evidence);
/******************************************************************************/
/* */
/* Glitch Characterization phase */
/* */
/******************************************************************************/
/******************************************************************************/
/* */
/* Glitch Characterization phase */
/* */
/******************************************************************************/
if(data->glitchModelFlag)
{
if(data->glitchModelFlag)
{
data->runPhase = 1;
reset_priors(data, prior);
data->runPhase = 1;
reset_priors(data, prior);
//Initialize proposal ratios for GW search phase
chain->modelRate = 1; //even chance for each channel
chain->rjmcmcRate = 0; //model fixed to one wavelet
//Initialize proposal ratios for GW search phase
chain->modelRate = 1; //even chance for each channel
chain->rjmcmcRate = 0; //model fixed to one wavelet
data->cleanFlag = 0;
data->glitchFlag = 1;
data->signalFlag = 0;
data->cleanFlag = 0;
data->glitchFlag = 1;
data->signalFlag = 0;
initialize_chain(chain,1);
initialize_chain(chain,1);
reset_likelihood(data);
reset_priors(data, prior);
reset_model(data, chain, prior, model);
reset_likelihood(data);
reset_priors(data, prior);
reset_model(data, chain, prior, model);
RJMCMC(data, model, bayesline, chain, prior, &data->logZglitch, &data->varZglitch);
RJMCMC(data, model, bayesline, chain, prior, &data->logZglitch, &data->varZglitch);
/*if(data->glitchOnlyFlag==1 && strcmp(data->channels[0],"AdLIGO"))
outputFrameFile(runState->commandLine, data, model[chain->index[0]]);*/
/*if(data->glitchOnlyFlag==1 && strcmp(data->channels[0],"AdLIGO"))
outputFrameFile(runState->commandLine, data, model[chain->index[0]]);*/
//shut off glitch model for checkpointing
data->glitchModelFlag = 0;
}
fprintf(evidence,"glitch %lg %lg\n",data->logZglitch,data->varZglitch);
fflush(evidence);
//shut off glitch model for checkpointing
data->glitchModelFlag = 0;
}
fprintf(evidence,"glitch %lg %lg\n",data->logZglitch,data->varZglitch);
fflush(evidence);
/******************************************************************************/
/* */
/* Signal Characterization phase */
/* Gaussian Noise Characterization phase */
/* */
/******************************************************************************/
if(data->signalModelFlag)
if(data->noiseModelFlag)
{
data->runPhase = 1;
reset_priors(data, prior);
if(data->bayesLineFlag) data->psdFitFlag = 0;
//Initialize proposal ratios for GW search phase
chain->modelRate = 1; //always do signal model updates
chain->rjmcmcRate = 0; //model fixed to one wavelet for extrinsic search
chain->rjmcmcRate = 0.0;
//set flags for no wavelets in model
data->cleanFlag = 0;
data->glitchFlag = 0;
data->signalFlag = 1;
data->signalFlag = 0;
initialize_chain(chain,1);
......@@ -779,12 +785,14 @@ int main(int argc, char *argv[])
reset_priors(data, prior);
reset_model(data, chain, prior, model);
RJMCMC(data, model, bayesline, chain, prior, &data->logZsignal, &data->varZsignal);
printf("characterizing Gaussian noise model...\n\n");
//shut off signal model for checkpointing
data->signalModelFlag = 0;
RJMCMC(data, model, bayesline, chain, prior, &data->logZnoise, &data->varZnoise);
//shut off noise model for checkpointing
data->noiseModelFlag = 0;
}
fprintf(evidence,"signal %lg %lg\n",data->logZsignal,data->varZsignal);
fprintf(evidence,"noise %lg %lg\n",data->logZnoise,data->varZnoise);
fflush(evidence);
/******************************************************************************/
......@@ -884,8 +892,8 @@ void print_help_message(void)
fprintf(stdout," --- Run parameters -------------------------------------------------------------\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --segment-start GPS start time of segment (trigtime + 2 - seglen) \n");
fprintf(stdout," --Niter number of iterations (4000000)\n");
fprintf(stdout," --Nchain number of parallel chains (15)\n");
fprintf(stdout," --Niter number of iterations (2000000)\n");
fprintf(stdout," --Nchain number of parallel chains (20)\n");
fprintf(stdout," --Ncycle number of model updates per RJMCMC iteration (100)\n");
fprintf(stdout," --Nburnin number of burn-in iterations (50000)\n");
fprintf(stdout," --maxLogL use maximized likelihood during burnin\n");
......@@ -935,7 +943,7 @@ void print_help_message(void)
fprintf(stdout," --uniformAmplitudePrior don't use SNR-dependent amplitude prior\n");
fprintf(stdout," --noSignalAmplitudePrior use same SNR-dependent prior for signal & glitch model\n");
fprintf(stdout," --noAmplitudeProposal don't draw from SNR-dependent amplitude prior\n");
fprintf(stdout," --fixExtrinsicAmplitude update wavelet amplitudes in extrinsic update (FIXED?)\n");
fprintf(stdout," --varyExtrinsicAmplitude update wavelet amplitudes in extrinsic update (FIXED?)\n");
fprintf(stdout," --noClusterProposal disable clustering & TF density proposal \n");
fprintf(stdout," --clusterWeight fractional weight for TF to proximity proposal (0.5)\n");
fprintf(stdout," --ampPriorPeak SNR where amplitude prior peaks (5)\n");
......
......@@ -356,6 +356,7 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
{
fprintf(stdout, "g%iSNR=%.0f\n", ifo, GNR[ifo]);
if(GNR[ifo]>100.) NCflag++;
if(GNR[ifo]>200.) NCflag++;
}
......@@ -403,6 +404,7 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
SNRnet = network_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model->response, model->invSnf, model->Sn, data->NI);
for(ifo=0; ifo<NI; ifo++) SNRifo[ifo] = detector_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model->response[ifo], model->invSnf[ifo], model->Sn[ifo]);
if(SNRnet>100.) NCflag++;
if(SNRnet>200.) NCflag++;
sprintf(filename,"%ssnr/snr.dat",data->runName);
snrFile=fopen(filename,"w");
......@@ -441,7 +443,9 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
fprintf(stdout,"BayesWave internal injection with SNR>100\n");
fprintf(stdout,"Cleaning phase already disabled\n");
fprintf(stdout,"Manually increasing number of chains to %i\n\n",*NC+5);
*NC+=5;
if(NCflag>0)*NC+=5;
if(NCflag>1)*NC+=5;
}
data->signalFlag = signalFlag;
......@@ -477,7 +481,7 @@ void getChannels(ProcessParamsTable *commandLine, char **channel)
{
channel[nifo] = (char*)malloc(strlen(ppt->value) + 1 * sizeof(char));
strcpy(channel[nifo], ppt->value);
printf("here\n");
//This removes the ifo name from the channel string
memmove(channel[nifo], channel[nifo]+3, strlen(channel[nifo]));
}
......@@ -1260,11 +1264,11 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
ppt = LALInferenceGetProcParamVal(commandLine, "--Nchain");
if(ppt) chain->NC = atoi(ppt->value);
else chain->NC = 15;
else chain->NC = 20;
ppt = LALInferenceGetProcParamVal(commandLine, "--Niter");
if(ppt) chain->count = atoi(ppt->value);
else chain->count = 4000000;
else chain->count = 2000000;
ppt = LALInferenceGetProcParamVal(commandLine, "--Ncycle");
if(ppt) chain->cycle = atoi(ppt->value);
......@@ -1316,7 +1320,7 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
ppt = LALInferenceGetProcParamVal(commandLine, "--tempSpacing");
if(ppt) chain->tempStep = (double)atof(ppt->value);
else chain->tempStep = 1.5;
else chain->tempStep = 1.2;
ppt = LALInferenceGetProcParamVal(commandLine, "--noAdaptTemperature");
if(ppt) data->adaptTemperatureFlag = 0;
......
......@@ -91,21 +91,24 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
if(data->adaptTemperatureFlag)// && data->runPhase!=0)
{
chain->temperature[chain->NC-1] = 1.0e6;
chain->dT[0] = pow(chain->temperature[chain->NC-1],(1./(double)(NC-1)));
//chain->dT[0] = pow(chain->temperature[chain->NC-1],(1./(double)(NC-1)));
//chain->dT[0] = pow(chain->temperature[chain->NC-1],(1./(double)(NC-1))*(1./(double)(NC-1))*(1./(double)(NC-1)));
chain->dT[0] = pow(chain->temperature[chain->NC-1],(1./(double)(NC-1))*(1./(double)(NC-1)));
for(ic=0; ic<chain->NC-1; ic++)
{
if(chain->NC==15)
{
//if(chain->NC==15)
//{
//chain->temperature[ic]=chain->temperature[ic-1]*chain->dT[0];
//TODO: DO better initializing chain temperature ladder (hard-wired for NCmax=15!)
chain->temperature[ic]=exp((double)(ic*ic)/(double)(((chain->NC-1)*(chain->NC-1))/(log(1.0e6)))); //empirically determined initial guess for temp spacing
chain->temperature[ic] = 1./(1.0 - 0.07*(double)ic);
}
else
{
chain->temperature[ic] = pow(chain->dT[0],ic);
}
//chain->temperature[ic]=exp((double)(ic*ic)/(double)(((chain->NC-1)*(chain->NC-1))/(log(1.0e6)))); //empirically determined initial guess for temp spacing
//chain->temperature[ic] = 1./(1.0 - 0.07*(double)ic);
//}
//else
//{
//chain->temperature[ic] = pow(chain->dT[0],ic*ic*ic);
chain->temperature[ic] = pow(chain->dT[0],ic*ic);
//}
}
}
......@@ -1270,7 +1273,18 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
}
//adapt PTMCMC spacing and reset various counters
if(chain->mc<M/2 && data->adaptTemperatureFlag) adapt_temperature_ladder(chain, NC);
if(chain->mc<M/2 && data->adaptTemperatureFlag)
{
adapt_temperature_ladder(chain, NC);
if(chain->mc%100000==0)
{
for(ic=0; ic<NC; ic++)
{
chain->ptprop[ic] = 1;
chain->ptacc[ic] = 1;
}
}
}
}//end NC>1 condition for parallel tempering
......@@ -1398,8 +1412,16 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
printf("Trapezoid Rule logZ = %g +/- %g\n",logZ,sqrt(varZ));
if(data->splineEvidenceFlag)
//check for convergence of PTMCMC ladder
int convergence_check=0;
for(ic=0; ic<NC-1; ic++)
{
if((double)chain->ptacc[ic]/(double)chain->ptprop[ic]<1.e-1) convergence_check++;
}
if(data->splineEvidenceFlag && convergence_check==0)
{
/******************************************************************************/
/* */
/* Monte Carlo integration of <logL> v beta */
......@@ -1410,6 +1432,11 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
printf("Thermodynamic Integration logZ = %g +/- %g\n",logZ,sqrt(varZ));
}
else
{
printf("Warning: Trapezoid Rule integration used for logZ\n");
printf(" Spline integration skipped due to flag or poor mixing\n");
}
}
/******************************************************************************/
......
......@@ -1212,7 +1212,7 @@ void initialize_TF_proposal(struct Data *data, struct Prior *prior, struct TimeF
tf->N = data->N;
//4s grid for Q
tf->nQ = (int)(prior->range[2][1] - prior->range[2][0])/4;
tf->nQ = (int)(prior->range[2][1] - prior->range[2][0]);
//4 Hz grid for frequency
tf->nf = (int)(prior->range[1][1] - prior->range[1][0])/4;
......
......@@ -201,21 +201,33 @@ void uniform_orientation_proposal(double *y, long *seed)
void intrinsic_halfcycle_proposal(double *x, double *y, long *seed)
{
double n;
double t0,f0,ph,dt,dp;
t0 = x[0];
f0 = x[1];
ph = x[4];
//shift time by half period
dt = 0.5/f0*gasdev2(seed);
dp = LAL_PI*gasdev2(seed);
//if(ran2(seed)<0.5) dt*=-1;
y[0] = t0 + dt;
y[4] = ph + dp;
// if(ran2(seed)<0.9)
// {
//choose number of half period shifts
n = -2.0 + floor(5.0*ran2(seed));
//shift time by half period
dt = (n/2.)/f0*(1. + 0.1*gasdev2(seed));
dp = 0.0;
if((int)n%2!=0) dp = LAL_PI*(1. + 0.1*gasdev2(seed));
// }
//
// else{
// dt = 0.5/f0*gasdev2(seed);
// dp = LAL_PI*gasdev2(seed);
//
// }
y[0] = t0 + dt;
y[4] = ph + dp;
while(y[4] > LAL_TWOPI) y[4] -=LAL_TWOPI;
while(y[4] < LAL_TWOPI) y[4] +=LAL_TWOPI;
......@@ -305,7 +317,7 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
break;
case 5:
//TODO: HACK! testing extrinsic amplitude update
range = 0.0;
range = 2.0;
break;
}
......@@ -320,6 +332,7 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
paramsP[i] += epsilon[i];
paramsM[i] -= epsilon[i];
logLP = EvaluateExtrinsicMarkovianLogLikelihood(projectionP, paramsP, invSnf, Sn, geo, glitch, data, data->fmin, data->fmax);
logLM = EvaluateExtrinsicMarkovianLogLikelihood(projectionM, paramsM, invSnf, Sn, geo, glitch, data, data->fmin, data->fmax);
epsilon[i] = 0.1/sqrt(-(logLM+logLP)/(epsilon[i]*epsilon[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