Commit 6b37b212 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

Almost ready for O2 review:

-Clean up warnings
-Optimizations in Extrinsic Update
-New flags to toggle new features
-New defaults for RJ proposals
-Added support for smaller time windows




git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@601 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent 226da92f
......@@ -2425,7 +2425,7 @@ void BayesLineSetup(struct BayesLineParams *bptr, double *freqData, double fmin,
}
void BayesLineSearch(struct BayesLineParams *bptr, double *freqData, double fmin, double fmax, double deltaT, double Tobs)
void BayesLineSearch(struct BayesLineParams *bayesline)
{
int i;
......@@ -2441,20 +2441,20 @@ void BayesLineSearch(struct BayesLineParams *bptr, double *freqData, double fmin
/* */
/******************************************************************************/
BayesLineNonMarkovianFit(bptr, &jj);
BayesLineNonMarkovianFit(bayesline, &jj);
// maximum number of terms in the Lorentzian model
bptr->data->tmax = 4*bptr->lines_full->n;
bayesline->data->tmax = 4*bayesline->lines_full->n;
if(bptr->data->tmax < 20) bptr->data->tmax = 20;
if(bptr->lines_full->n < 1 ) bptr->lines_full->n = 1;
if(bayesline->data->tmax < 20) bayesline->data->tmax = 20;
if(bayesline->lines_full->n < 1 ) bayesline->lines_full->n = 1;
// Re-set dataParams structure to use full dataset
bptr->data->flow = bptr->data->fmin;
imax = (int)(floor(bptr->data->fhigh*bptr->data->Tobs));
imin = (int)(floor(bptr->data->flow*bptr->data->Tobs));
bptr->data->ncut = imax-imin;
bayesline->data->flow = bayesline->data->fmin;
imax = (int)(floor(bayesline->data->fhigh*bayesline->data->Tobs));
imin = (int)(floor(bayesline->data->flow*bayesline->data->Tobs));
bayesline->data->ncut = imax-imin;
/******************************************************************************/
/* */
......@@ -2462,23 +2462,23 @@ void BayesLineSearch(struct BayesLineParams *bptr, double *freqData, double fmin
/* */
/******************************************************************************/
BayesLineMarkovianSplineOnly(bptr, bptr->spline->n, jj);
BayesLineMarkovianSplineOnly(bayesline, bayesline->spline->n, jj);
for(j=0; j<bptr->spline->n; j++)
for(j=0; j<bayesline->spline->n; j++)
{
bptr->spline_x->points[j] = bptr->spline->points[j];
bptr->spline_x->data[j] = bptr->spline->data[j];
bayesline->spline_x->points[j] = bayesline->spline->points[j];
bayesline->spline_x->data[j] = bayesline->spline->data[j];
}
bptr->data->flow = bptr->data->fmin;
bayesline->data->flow = bayesline->data->fmin;
imin = (int)(floor(bptr->data->flow*bptr->data->Tobs));
imin = (int)(floor(bayesline->data->flow*bayesline->data->Tobs));
for(i=0; i< bptr->data->ncut; i++)
for(i=0; i< bayesline->data->ncut; i++)
{
j = i + imin - bptr->data->nmin;
bptr->spow[i] = bptr->power[j];
bptr->sfreq[i] = bptr->freq[j];
j = i + imin - bayesline->data->nmin;
bayesline->spow[i] = bayesline->power[j];
bayesline->sfreq[i] = bayesline->freq[j];
}
/******************************************************************************/
......@@ -2487,11 +2487,11 @@ void BayesLineSearch(struct BayesLineParams *bptr, double *freqData, double fmin
/* */
/******************************************************************************/
BayesLineMarkovianFocusedAnalysis(bptr);
BayesLineMarkovianFocusedAnalysis(bayesline);
double dan;
BayesLineLorentzSplineMCMC(bptr, 1.0, 200000, 0, 0, &dan);
BayesLineLorentzSplineMCMC(bayesline, 1.0, 200000, 0, 0, &dan);
}
......
......@@ -112,7 +112,7 @@ void BayesLineFree(struct BayesLineParams *bptr);
void BayesLineSetup(struct BayesLineParams *bptr, double *freqData, double fmin, double fmax, double deltaT, double Tobs);
void BayesLineRJMCMC(struct BayesLineParams *bayesline, double *freqData, double *psd, double *invpsd, double *splinePSD, int N, int cycle, double beta, int priorFlag);
void BayesLineSearch(struct BayesLineParams *bptr, double *freqData, double fmin, double fmax, double deltaT, double Tobs);
void BayesLineSearch(struct BayesLineParams *bayesline);
void BayesLineNonMarkovianFit (struct BayesLineParams *bayesline, int *nj);
void BayesLineLorentzSplineMCMC (struct BayesLineParams *bayesline, double heat, int steps, int focus, int priorFlag, double *dan);
......
......@@ -142,13 +142,8 @@ int main(int argc, char *argv[])
//Output Fourier domain data for BayesLine
dataPtr = runState->data;
//TODO: We have to pick frequency resolution for TF proposal. What's optimal?
//TODO: A continuous density would be better
// and do-able -- store (fraction of) burnin wavelets, ~KDE (has Q information)
double fmin = dataPtr->fLow;
double fmax = ceil(dataPtr->fHigh);
double deltaT = dataPtr->timeData->deltaT;
//set frequency resolution to be 32 Hz because why not?
double df,dt;
......@@ -208,25 +203,6 @@ int main(int argc, char *argv[])
print_run_stats(chain->runFile, data, chain);
fclose(chain->runFile);
//likelhood is a pointer to a function which returns logL
if(data->constantLogLFlag)
{
// logL=const to test detailed balance
data->intrinsic_likelihood = EvaluateConstantLogLikelihood;
data->extrinsic_likelihood = EvaluateExtrinsicConstantLogLikelihood;
}
else
{
// Start with search-phase log likelihood
data->intrinsic_likelihood = EvaluateSearchLogLikelihood;
// Jump right into characterization likelihood if noise-only run
if(!data->glitchFlag && !data->signalFlag) data->intrinsic_likelihood = EvaluateMarkovianLogLikelihood;
// No ned (or good implementatino) of extrinsic search likelihood
data->extrinsic_likelihood = EvaluateExtrinsicMarkovianLogLikelihood;
}
data->detector = malloc(NI*sizeof(LALDetector*));
data->epoch = runState->data->epoch;
......@@ -347,9 +323,8 @@ int main(int argc, char *argv[])
bayesline[0][ifo]->power[i] = (asd[2*j]*asd[2*j]+asd[2*j+1]*asd[2*j+1]);
}
//TODO: Is passing d(f) to BayesLineSearch etc. redundant?
fprintf(stdout,"BayesLine search phase for IFO %i\n", ifo);
BayesLineSearch(bayesline[0][ifo], asd, fmin, fmax, deltaT, Tobs);
BayesLineSearch(bayesline[0][ifo]);
fprintf(stdout,"BayesLine characterization phase for IFO %i\n", ifo);
BayesLineRJMCMC(bayesline[0][ifo], asd, model[0]->Snf[ifo], model[0]->invSnf[ifo], model[0]->SnS[ifo], N, 100000, 1.0, 2);
......@@ -899,16 +874,17 @@ void print_help_message(void)
fprintf(stdout," ----------------------------------------------------------------------------------\n");
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," --Ncycle number of model updates per RJMCMC iteration (100)\n");
fprintf(stdout," --Nburnin number of burn-in iterations (10000)\n");
fprintf(stdout," --Nburnin number of burn-in iterations (50000)\n");
fprintf(stdout," --maxLogL use maximized likelihood during burnin\n");
fprintf(stdout," --chainseed random number seed for Markov chain (1234)\n");
fprintf(stdout," --runName run name for output files\n");
fprintf(stdout," --outputDir absolute path of where the output should be written\n");
fprintf(stdout," --0noise no noise realization\n");
fprintf(stdout," --zeroLogL logL = constant test\n");
fprintf(stdout," --restart restart run using psd and residual from file\n");
fprintf(stdout," --prior sample from prior using logL = constant test\n");
fprintf(stdout," --gnuplot output files for gnuplot animations\n");
fprintf(stdout," --verbose output hot chains\n");
fprintf(stdout," --window duration of time window for model characterization runs\n");
......@@ -944,13 +920,13 @@ void print_help_message(void)
fprintf(stdout," --clusterBeta cluster prior exp(-beta K) (4)\n");
fprintf(stdout," --clusterGamma occam penalty exp(gamma*J) (0)\n");
fprintf(stdout," --backgroundPrior name of 2-column bkg frequency distribution file\n");
fprintf(stdout," --orientationProposal enable MCMC proposal for psi/ecc\n");
fprintf(stdout," --noOrientationProposal disable MCMC proposal for psi/ecc\n");
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," --clusterProposal use clustering & TF density proposal (BROKEN)\n");
fprintf(stdout," --clusterWeight fractional weight for TF to proximity proposal (1.0)\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");
fprintf(stdout," --signalPriorPeak SNR where signal amplitude prior peaks (5)\n");
fprintf(stdout," --dimensionDecayRate e-folding time for number of wavelets (1000)\n");
......
......@@ -264,6 +264,7 @@ struct Data
int fixIntrinsicFlag;
int fixExtrinsicFlag;
int geocenterPSDFlag;
int clusterPriorFlag;
int backgroundPriorFlag;
int amplitudePriorFlag;
......@@ -272,6 +273,7 @@ struct Data
int extrinsicAmplitudeFlag;
int orientationProposalFlag;
int clusterProposalFlag;
int logClusterProposalFlag;
int adaptTemperatureFlag;
int splineEvidenceFlag;
int constantLogLFlag;
......@@ -279,6 +281,7 @@ struct Data
int verboseFlag;
int checkpointFlag;
int resumeFlag;
int searchFlag;
int rjFlag;
double dt;
......
......@@ -2408,7 +2408,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
maxAmp=-1.e60;
for(i=0; i<100; i++)
{
draw_wavelet_params(-1, params, data, prior->range, &chain->seed);
draw_wavelet_params(params, prior->range, &chain->seed);
data->signal_amplitude_proposal(params, model->SnGeo, 1.0, &chain->seed, data->Tobs, prior->range,prior->sSNRpeak);
if(params[3]<minAmp)minAmp=params[3];
if(params[3]>maxAmp)maxAmp=params[3];
......@@ -2435,7 +2435,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
maxAmp=-1.e60;
for(i=0; i<100; i++)
{
draw_wavelet_params(ifo, params, data, prior->range, &chain->seed);
draw_wavelet_params(params, prior->range, &chain->seed);
data->glitch_amplitude_proposal(params, model->Snf[ifo], 1.0, &chain->seed, data->Tobs, prior->range, prior->gSNRpeak);
if(params[3]<minAmp)minAmp=params[3];
if(params[3]>maxAmp)maxAmp=params[3];
......
......@@ -5,6 +5,11 @@
#include "BayesWaveModel.h"
#include "BayesWaveWavelet.h"
#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif
/* ********************************************************************************** */
/* */
......@@ -172,8 +177,9 @@ void InjectFromMDC(ProcessParamsTable *commandLine, LALInferenceIFOData *IFOdata
COMPLEX16FrequencySeries* injF=(COMPLEX16FrequencySeries *)XLALCreateCOMPLEX16FrequencySeries("injF",&IFOdata->timeData->epoch,0.0,IFOdata->freqData->deltaF,&lalDimensionlessUnit, IFOdata->freqData->data->length);
if(!injF) {
XLALPrintError("Unable to allocate memory for injection buffer\n");
XLAL_ERROR_VOID(XLAL_EFUNC);
fprintf(stdout,"Unable to allocate memory for injection buffer\n");
fflush(stdout);
exit(1);
}
REAL4 WinNorm = sqrt(IFOdata->window->sumofsquares/IFOdata->window->data->length);
......@@ -587,7 +593,7 @@ int outputFrameFile(ProcessParamsTable *commandLine, struct Data *data, struct M
/* */
/* ********************************************************************************** */
void write_gnuplot_script_header(FILE *script, char modelname[], double Tobs)
void write_gnuplot_script_header(FILE *script, char modelname[])
{
char filename[100];
sprintf(filename,"waveforms/%s_waveforms_animate.gpi",modelname);
......@@ -604,7 +610,7 @@ void write_gnuplot_script_header(FILE *script, char modelname[], double Tobs)
fprintf(script," \n");
}
void write_gnuplot_script_frame (FILE *script, char modelname[], char runname[], double Tobs, int phase, int signal, int glitch, int cycle, int NI)
void write_gnuplot_script_frame (FILE *script, char modelname[], int signal, int glitch, int cycle, int NI)
{
int i;
......@@ -632,7 +638,7 @@ void write_gnuplot_script_frame (FILE *script, char modelname[], char runname[],
fflush(script);
}
void write_bayesline_gnuplot_script_frame(FILE *script, char modelname[], double Tobs, int cycle)
void write_bayesline_gnuplot_script_frame(FILE *script, char modelname[], int cycle)
{
fprintf(script,"set multiplot \n");
fprintf(script,"set size 1,0.5 \n");
......@@ -654,7 +660,7 @@ void write_bayesline_gnuplot_script_frame(FILE *script, char modelname[], double
fflush(script);
}
void write_bayeswave_gnuplot_script_frame(FILE *script, char modelname[], char runname[], double Tobs, double fmin, double fmax, double SAmin, double SAmax, int phase, int signal, int glitch, int cycle, int NI)
void write_bayeswave_gnuplot_script_frame(FILE *script, char modelname[], double Tobs, double fmin, double fmax, int phase, int signal, int glitch, int cycle, int NI)
{
int i;
......@@ -809,6 +815,11 @@ void print_run_flags(FILE *fptr, struct Data *data, struct Prior *prior)
else fprintf(fptr, "DISABLED");
fprintf(fptr, "\n");
fprintf(fptr, " maximized logL is ........ ");
if(data->searchFlag) fprintf(fptr, "ENABLED");
else fprintf(fptr, "DISABLED");
fprintf(fptr, "\n");
fprintf(fptr, " PSD fitting is ........... ");
if(data->psdFitFlag) fprintf(fptr, "ENABLED");
else fprintf(fptr, "DISABLED");
......@@ -867,7 +878,11 @@ void print_run_flags(FILE *fptr, struct Data *data, struct Prior *prior)
fprintf(fptr, "\n");
fprintf(fptr, " density proposal is ...... ");
if(data->clusterProposalFlag) fprintf(fptr, "ENABLED");
if(data->clusterProposalFlag)
{
if(data->logClusterProposalFlag) fprintf(fptr, "SNR^2/2");
else fprintf(fptr, "exp(SNR^2/2)-1");
}
else fprintf(fptr, "DISABLED");
fprintf(fptr, "\n");
......@@ -904,7 +919,7 @@ void print_run_flags(FILE *fptr, struct Data *data, struct Prior *prior)
fprintf(fptr, "\n");
}
void print_chain_files(struct Data *data, struct Chain *chain, struct Model **model, struct BayesLineParams ***bayesline, int ic, double logZ, double varZ)
void print_chain_files(struct Data *data, struct Chain *chain, struct Model **model, struct BayesLineParams ***bayesline, int ic)
{
int ifo;
int NI = data->NI;
......@@ -980,7 +995,7 @@ void print_glitch_model(FILE *fptr, struct Wavelet *glitch)
fprintf(fptr,"\n");
}
void print_chain_status(struct Data *data, struct Chain *chain, struct Model **model, int searchFlag, int NC)
void print_chain_status(struct Data *data, struct Chain *chain, struct Model **model, int searchFlag)
{
int ifo;
int NI = data->NI;
......@@ -1082,7 +1097,7 @@ void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf,
fclose(waveout);
}
void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, double eta, double Tobs, int imin, int imax, double tmin, double tmax, double fmax)
void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, double eta, double Tobs, int imin, int imax, double tmin, double tmax)
{
/*
TODO: invFFT comes out time-reversed.
......@@ -1228,6 +1243,10 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
if(ppt) chain->burn = atoi(ppt->value);
else chain->burn = 50000;
ppt = LALInferenceGetProcParamVal(commandLine, "--maxLogL");
if(ppt) data->searchFlag = 1;
else data->searchFlag = 0;
ppt = LALInferenceGetProcParamVal(commandLine, "--chainseed");
if(ppt) chain->seed = atoi(ppt->value);
else chain->seed = 1234;
......@@ -1355,15 +1374,19 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
if(ppt) data->fixIntrinsicFlag = 1;
else data->fixIntrinsicFlag = 0;
ppt = LALInferenceGetProcParamVal(commandLine, "--orientationProposal");
if(ppt) data->orientationProposalFlag = 1;
else data->orientationProposalFlag = 0;
ppt = LALInferenceGetProcParamVal(commandLine, "--noOrientationProposal");
if(ppt) data->orientationProposalFlag = 0;
else data->orientationProposalFlag = 1;
//TODO: Disabling the extrinsic amplitude updates (severly biasing Q, A)
ppt = LALInferenceGetProcParamVal(commandLine, "--varyExtrinsicAmplitude");
if(ppt) data->extrinsicAmplitudeFlag = 1;
else data->extrinsicAmplitudeFlag = 0;
ppt = LALInferenceGetProcParamVal(commandLine, "--fixGeocenterPSD");
if(ppt) data->geocenterPSDFlag = 1;
else data->geocenterPSDFlag = 0;
ppt = LALInferenceGetProcParamVal(commandLine, "--uniformAmplitudePrior");
if(ppt) data->amplitudePriorFlag = 0;
else data->amplitudePriorFlag = 1;
......@@ -1388,14 +1411,17 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
if(ppt) prior->Dtau = (double)atof(ppt->value);
else prior->Dtau=1000.0;
//TODO: Cluster proposal is off by default (problems in normalization?)
ppt = LALInferenceGetProcParamVal(commandLine, "--clusterProposal");
if(ppt) data->clusterProposalFlag = 1;
else data->clusterProposalFlag = 0;
ppt = LALInferenceGetProcParamVal(commandLine, "--noClusterProposal");
if(ppt) data->clusterProposalFlag = 0;
else data->clusterProposalFlag = 1;
ppt = LALInferenceGetProcParamVal(commandLine, "--noLogClusterProposal");
if(ppt) data->logClusterProposalFlag = 0;
else data->logClusterProposalFlag = 1;
ppt = LALInferenceGetProcParamVal(commandLine, "--clusterWeight");
if(ppt) data->TFtoProx = 1.0-(double)atof(ppt->value);
else data->TFtoProx = 0.0;
else data->TFtoProx = 0.5;
ppt = LALInferenceGetProcParamVal(commandLine, "--runName");
if(ppt) sprintf(data->runName,"%s_",ppt->value);
......@@ -1417,7 +1443,7 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
if(ppt) data->noiseSimFlag = 0;
else data->noiseSimFlag = 1;
ppt = LALInferenceGetProcParamVal(commandLine, "--zeroLogL");
ppt = LALInferenceGetProcParamVal(commandLine, "--prior");
if(ppt) data->constantLogLFlag = 1;
else data->constantLogLFlag = 0;
......@@ -1513,6 +1539,15 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
data->splineEvidenceFlag = 0;
}
if(data->psdFitFlag == 1 && data->noiseSimFlag == 0)
{
fprintf(stdout,"************************ WARNING ************************\n");
fprintf(stdout,"0noise runs do not get along with PSD fitting \n");
fprintf(stdout,"Disabling PSD fitting \n");
fprintf(stdout,"Silence warning with --noPSDfit \n");
fprintf(stdout,"*********************************************************\n");
data->psdFitFlag = 0;
}
}
void parse_bayesline_parameters(struct Data *data, struct Model *model, FILE **splinechain, FILE **linechain, double **psd)
......
#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif
/* ********************************************************************************** */
/* */
/* Data handling and injection routines */
......@@ -17,25 +24,25 @@ int outputFrameFile(ProcessParamsTable *commandLine, struct Data *data, struct M
/* */
/* ********************************************************************************** */
void write_gnuplot_script_header(FILE *script, char modelname[], double Tobs);
void write_gnuplot_script_frame (FILE *script, char modelname[], char runname[], double Tobs, int phase, int signal, int glitch, int cycle, int NI);
void write_bayesline_gnuplot_script_frame(FILE *script, char modelname[], double Tobs, int cycle);
void write_bayeswave_gnuplot_script_frame(FILE *script, char modelname[], char runname[], double Tobs, double fmin, double fmax, double SAmin, double SAmax, int phase, int signal, int glitch, int cycle, int NI);
void write_gnuplot_script_header(FILE *script, char modelname[]);
void write_gnuplot_script_frame (FILE *script, char modelname[], int signal, int glitch, int cycle, int NI);
void write_bayesline_gnuplot_script_frame(FILE *script, char modelname[], int cycle);
void write_bayeswave_gnuplot_script_frame(FILE *script, char modelname[], double Tobs, double fmin, double fmax, int phase, int signal, int glitch, int cycle, int NI);
void write_evidence_gnuplot_script(FILE *script, char modelname[]);
void print_bayesline_spectrum(char filename[], double *f, double *power, double *Sbase, double *Sline, int N);
void print_version(FILE *fptr);
void print_run_stats(FILE *fptr, struct Data *data, struct Chain *chain);
void print_run_flags(FILE *fptr, struct Data *data, struct Prior *prior);
void print_chain_files(struct Data *data, struct Chain *chain, struct Model **model, struct BayesLineParams ***bayesline, int ic, double logZ, double varZ);
void print_chain_status(struct Data *data, struct Chain *chain, struct Model **model, int searchFlag, int NC);
void print_chain_files(struct Data *data, struct Chain *chain, struct Model **model, struct BayesLineParams ***bayesline, int ic);
void print_chain_status(struct Data *data, struct Chain *chain, struct Model **model, int searchFlag);
void print_model(FILE *fptr, struct Data *data, struct Chain *chain, struct Model *model);
void print_signal_model(FILE *fptr, struct Model *model);
void print_glitch_model(FILE *fptr, struct Wavelet *glitch);
void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf, double eta, double Tobs, int imin, int imax, double tmin, double tmax);
void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, double eta, double Tobs, int imin, int imax, double tmin, double tmax, double fmax);
void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, double eta, double Tobs, int imin, int imax, double tmin, double tmax);
void print_frequency_domain_waveforms(char filename[], double *h, int N, double *Snf, double Tobs, double eta, int imin, int imax);
void print_frequency_domain_data(char filename[], double *h, int N, double Tobs, int imin, int imax);
......
......@@ -11,6 +11,12 @@
#include "BayesWaveWavelet.h"
#include "BayesWaveLikelihood.h"
#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif
/* ********************************************************************************** */
/* */
/* Likelihood Functions */
......@@ -641,8 +647,7 @@ double EvaluateSearchLogLikelihood(int typ, int ii, int det, struct Model *mx, s
return logLy;
}
//double EvaluateConstantLogLikelihood(int typ, int ii, int det, struct Model *mx, struct Model *my, double **range, struct Data *data)
double EvaluateConstantLogLikelihood(int typ, int ii, int det, struct Model *mx, struct Model *my, struct Prior *prior, struct Chain *chain, struct Data *data)
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)
{
// Rejection sample on prior range
struct Wavelet *wave = NULL;
......@@ -1120,7 +1125,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
}
double EvaluateExtrinsicConstantLogLikelihood(struct Network *projection, double *params, double **invSnf, double *Sn, double *geo, double **g, struct Data *data, double fmin, double fmax)
double EvaluateExtrinsicConstantLogLikelihood(UNUSED struct Network *projection, double *params, UNUSED double **invSnf, UNUSED double *Sn, UNUSED double *geo, UNUSED double **g, UNUSED struct Data *data, UNUSED double fmin, UNUSED double fmax)
{
if(extrinsic_checkrange(params)) return -1.0e60;
......
#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif
/* ********************************************************************************** */
/* */
/* Likelihood Functions */
......@@ -13,13 +19,13 @@ void recompute_residual(struct Data *data, struct Model **model, struct Chain *c
double Sum_Extreme(double **a, double **b, double **Sn, int n, double *delt, double *pshift, double Tobs, int NI, int imin, int imax, double t0);
double EvaluateSearchLogLikelihood (int typ, int ii, int det, struct Model *mx, struct Model *my, struct Prior *prior, struct Chain *chain, struct Data *data);
double EvaluateConstantLogLikelihood (int typ, int ii, int det, struct Model *mx, struct Model *my, struct Prior *prior, struct Chain *chain, struct Data *data);
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, double *Sn, double *geo, double **g, struct Data *data, double fmin, double fmax);
double EvaluateExtrinsicConstantLogLikelihood (struct Network *projection, double *params, double **invSnf, double *Sn, double *geo, double **g, struct Data *data, double fmin, double fmax);
double EvaluateExtrinsicConstantLogLikelihood (UNUSED struct Network *projection, double *params, UNUSED double **invSnf, UNUSED double *Sn, UNUSED double *geo, UNUSED double **g, UNUSED struct Data *data, UNUSED double fmin, UNUSED double fmax);
double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, double *params, double **invSnf, double *Sn, double *geo, double **g, struct Data *data, double fmin, double fmax);
void phase_blind_time_shift(double *corr, double *corrf, double *data1, double *data2, double *psd, int n);
\ No newline at end of file
void phase_blind_time_shift(double *corr, double *corrf, double *data1, double *data2, double *psd, int n);
......@@ -52,7 +52,7 @@ static void catch_alarm(UNUSED int sig, UNUSED siginfo_t *siginfo,UNUSED void *c
}
static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior *prior, struct Model **model, struct BayesLineParams ***bayesline, char modelname[])
static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior *prior, struct Model **model, char modelname[])
{
int i,j;
int ic,ifo;
......@@ -143,7 +143,7 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
model[ic]->size++;
for(j=1; j<=glitch[ifo]->size; j++)
{
draw_wavelet_params(ifo, glitch[ifo]->intParams[j], data, prior->range, &chain->seed);
draw_wavelet_params(glitch[ifo]->intParams[j], prior->range, &chain->seed);
if(data->amplitudePriorFlag) data->glitch_amplitude_proposal(glitch[ifo]->intParams[j], model[ic]->Snf[ifo], 1.0, &chain->seed, data->Tobs, prior->range, prior->gSNRpeak);
model[ic]->wavelet(glitch[ifo]->templates, glitch[ifo]->intParams[j], data->N, 1, Tobs);
......@@ -166,11 +166,11 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
current geocenter waveform and the detector response function.
*/
computeProjectionCoeffs(data, model[ic]->projection, model[ic]->extParams, data->fmin, data->fmax);
Shf_Geocenter_inv(data, model[ic]->projection, model[ic]->invSnf, model[ic]->SnGeo, model[ic]->extParams);
Shf_Geocenter_full(data, model[ic]->projection, model[ic]->Snf, model[ic]->SnGeo, model[ic]->extParams);
for(j=1; j<=signal->size; j++)
{
draw_wavelet_params(ifo, signal->intParams[j], data, prior->range, &chain->seed);
draw_wavelet_params(signal->intParams[j], prior->range, &chain->seed);
if(data->amplitudePriorFlag) data->signal_amplitude_proposal(signal->intParams[j], model[ic]->SnGeo, 1.0, &chain->seed, data->Tobs, prior->range, prior->sSNRpeak);
model[ic]->wavelet(signal->templates, signal->intParams[j], data->N, 1, Tobs);
}
......@@ -263,7 +263,7 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
}
}
static void resume_sampler(struct Data *data, struct Chain *chain, struct Prior *prior, struct Model **model, struct BayesLineParams ***bayesline)
static void resume_sampler(struct Data *data, struct Chain *chain, struct Model **model, struct BayesLineParams ***bayesline)
{
//TODO: set internal mcmc flags right based on chain->mc
int i;
......@@ -642,7 +642,7 @@ static void resume_sampler(struct Data *data, struct Chain *chain, struct Prior
free(h);
}
static void save_sampler(struct Data *data, struct Chain *chain, struct Prior *prior, struct Model **model, struct BayesLineParams ***bayesline, char modelname[])
static void save_sampler(struct Data *data, struct Chain *chain, struct Model **model, struct BayesLineParams ***bayesline, char modelname[])
{
int i;
int ic,ifo;
......@@ -974,7 +974,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
if(data->checkpointFlag && data->resumeFlag)
{
printf("try resuming...\n");
resume_sampler(data, chain, prior, model, bayesline);
resume_sampler(data, chain, model, bayesline);
data->resumeFlag = 0;
/******************************************************************************/
......@@ -1028,7 +1028,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b