Commit a3a05add authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

automatically rescale NC for loud injections

tweak to time/phase proposal for review



git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@611 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent 18ec9a97
......@@ -268,7 +268,7 @@ int main(int argc, char *argv[])
if(LALInferenceGetProcParamVal(runState->commandLine,"--BW-inject"))
{
fprintf(stdout, " ==== BayesWaveInjection(): started. ====\n");
BayesWaveInjection(runState->commandLine, data, chain, prior, psd);
BayesWaveInjection(runState->commandLine, data, chain, prior, psd, &NC);
fprintf(stdout, " ==== BayesWaveInjection(): finished. ====\n\n");
}
......@@ -555,7 +555,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)>50.0) chain->NC+=5;
if(sqrt(rSNR)>100.0) chain->NC+=5;
/******************************************************************************/
/* */
......@@ -653,6 +653,15 @@ int main(int argc, char *argv[])
fprintf(snrFile,"%.12g %.12g\n",SNRmin,SNRinj[NI]);
fclose(snrFile);
if(data->cleanModelFlag==0 && SNRinj[NI]>100.)
{
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);
}
......@@ -926,7 +935,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," --varyExtrinsicAmplitude update wavelet amplitudes in extrinsic update (BROKEN)\n");
fprintf(stdout," --fixExtrinsicAmplitude 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");
......
......@@ -245,10 +245,12 @@ void InjectFromMDC(ProcessParamsTable *commandLine, LALInferenceIFOData *IFOdata
}
void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, struct Chain *chain, struct Prior *prior, double **psd)
void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, struct Chain *chain, struct Prior *prior, double **psd, int *NC)
{
int i,n,ifo;
int NCflag=0;
int N = data->N;
int NI = data->NI;
......@@ -324,7 +326,6 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
for(ifo=0; ifo<NI; ifo++)
{
printf("ifo=%i\n",ifo);
for(n=0; n<Nsample+1; n++) fgets(burnrows,100000,glitchParams[ifo]);
......@@ -350,7 +351,14 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
fprintf(snrFile,"\n");
fclose(snrFile);
for(ifo=0; ifo<NI; ifo++)fprintf(stdout, " g%iSNR=%.0f\n", ifo, GNR[ifo]);
for(ifo=0; ifo<NI; ifo++)
{
fprintf(stdout, "g%iSNR=%.0f\n", ifo, GNR[ifo]);
if(GNR[ifo]>100.) NCflag++;
}
//Add model to data
......@@ -394,7 +402,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++;
sprintf(filename,"%ssnr/snr.dat",data->runName);
snrFile=fopen(filename,"w");
......@@ -406,9 +414,9 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
fprintf(snrFile,"%g\n",SNRnet);
fclose(snrFile);
fprintf(stdout, " SNR=%.0f\n",SNRnet);
fprintf(stdout, "SNR=%.0f\n",SNRnet);
for(ifo=0; ifo<NI; ifo++)
fprintf(stdout, " IFO%iSNR=%.0f\n", ifo, SNRifo[ifo]);
fprintf(stdout, "IFO%iSNR=%.0f\n", ifo, SNRifo[ifo]);
//Add model to data
......@@ -418,6 +426,7 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
free_model(model, data, prior);
}
fprintf(stdout,"\n");
for(ifo=0; ifo<NI; ifo++)
{
......@@ -427,6 +436,14 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
free(hinj);
free(ginj);
if(data->cleanModelFlag==0 && NCflag>0)
{
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;
}
data->signalFlag = signalFlag;
data->glitchFlag = glitchFlag;
}
......
......@@ -14,7 +14,7 @@
REAL8TimeSeries *readTseries(CHAR *cachefile, CHAR *channel, LIGOTimeGPS start, REAL8 length);
void InjectFromMDC(ProcessParamsTable *commandLine, LALInferenceIFOData *IFOdata, double *SNR);
void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, struct Chain *chain, struct Prior *prior, double **psd);
void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, struct Chain *chain, struct Prior *prior, double **psd, int *NC);
void getChannels(ProcessParamsTable *commandLine, char **channel);
int outputFrameFile(ProcessParamsTable *commandLine, struct Data *data, struct Model *model);
......
......@@ -2346,7 +2346,6 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
}
}
if(!data->extrinsicAmplitudeFlag) paramsy[5]=1.0;
//else paramsy[5] = 1.0 + 0.1*gasdev2(seed);
//Rescale intrinsic amplitudes and phases with extrinsic parameters
for(i=1; i<ienddim; i++)
......@@ -2366,7 +2365,6 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if(data->geocenterPSDFlag)Shf_Geocenter(data, model_x, SnGeoy, paramsy);
for(i=1; i<ienddim; i++) logpy += (data->signal_amplitude_prior(intParams[i],SnGeoy, 1.0,data->Tobs, prior->sSNRpeak));
}
//logpy += dim*log(paramsy[5])/2.;
logH = (logLy - logLx)*chain->beta + logpy - logpx + logJ;
......
......@@ -384,7 +384,7 @@ int main(int argc, char *argv[])
if(LALInferenceGetProcParamVal(runState->commandLine,"--BW-inject"))
{
fprintf(stdout, " ==== BayesWaveInjection(): started. ====\n");
BayesWaveInjection(runState->commandLine, data, chain, prior, psd);
BayesWaveInjection(runState->commandLine, data, chain, prior, psd, &chain->NC);
fprintf(stdout, " ==== BayesWaveInjection(): finished. ====\n\n");
injectionFlag = 1;
}
......
......@@ -208,10 +208,10 @@ void intrinsic_halfcycle_proposal(double *x, double *y, long *seed)
ph = x[4];
//shift time by half period
dt = 0.5/f0;
dp = LAL_PI;
dt = 0.5/f0*gasdev2(seed);
dp = LAL_PI*gasdev2(seed);
if(ran2(seed)<0.5) dt*=-1;
//if(ran2(seed)<0.5) dt*=-1;
y[0] = t0 + dt;
y[4] = ph + dp;
......@@ -305,7 +305,7 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
break;
case 5:
//TODO: HACK! testing extrinsic amplitude update
range = 2.0;
range = 0.0;
break;
}
......@@ -320,9 +320,6 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
paramsP[i] += epsilon[i];
paramsM[i] -= epsilon[i];
computeProjectionCoeffs(data, projectionP, paramsP, data->fmin, data->fmax);
computeProjectionCoeffs(data, projectionM, paramsM, data->fmin, data->fmax);
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