Skip to content
GitLab
Commit 64c15c26 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

Merge branch 'no-psd-fit' into 'master' 

No psd fit

See merge request lscsoft/bayeswave!79
parents a61f3b04 6df706f3
Hide whitespace changes
Inline Side-by-side
src/BayesWave.c
View file @ 64c15c26
......@@ -605,7 +605,6 @@ int main(int argc, char *argv[])
data->signalFlag = 0;
data->runPhase = 0;
if(data->bayesLineFlag) data->psdFitFlag = 0;
initialize_chain(chain,1);
......@@ -753,8 +752,6 @@ int main(int argc, char *argv[])
//Shut off BayesLine in favor of PSD fitting
data->cleanModelFlag = 0;
//data->bayesLineFlag = 0;
//data->psdFitFlag = 1;
}
......@@ -943,7 +940,6 @@ int main(int argc, char *argv[])
data->runPhase = 1;
reset_priors(data, prior);
if(data->bayesLineFlag) data->psdFitFlag = 0;
//Initialize proposal ratios for GW search phase
chain->rjmcmcRate = 0.0;
......@@ -995,8 +991,6 @@ int main(int argc, char *argv[])
data->glitchFlag = 1;
data->signalFlag = 1;
if(data->bayesLineFlag) data->psdFitFlag = 0;
initialize_chain(chain,1);
reset_likelihood(data);
......
src/BayesWave.h
View file @ 64c15c26
......@@ -115,7 +115,7 @@ struct Model
int size;
int Npol;
double *Sn;
//double *Sn;
double **Snf;
double **invSnf;
double **SnS;
......@@ -278,7 +278,6 @@ struct Data
int Dmin;
int runPhase;
int psdFitFlag;
int noiseSimFlag;
int bayesLineFlag;
int stochasticFlag;
......@@ -358,14 +357,14 @@ struct Data
double varZfull;
void (*signal_amplitude_proposal)(double *, double *, double, gsl_rng *, double, double **, double);
void (*glitch_amplitude_proposal)(double *, double *, double, gsl_rng *, double, double **, double);
void (*signal_amplitude_proposal)(double *, double *, gsl_rng *, double, double **, double);
void (*glitch_amplitude_proposal)(double *, double *, gsl_rng *, double, double **, double);
double (*signal_amplitude_prior) (double *, double *, double, double);
double (*glitch_amplitude_prior) (double *, double *, double, double);
double (*intrinsic_likelihood)(int, int, int, struct Model*, struct Model*, struct Prior*, struct Chain*, struct Data*);
double (*extrinsic_likelihood)(struct Network*, double *, double **, double *, struct Wavelet **, double **, struct Data*, double, double);
double (*extrinsic_likelihood)(struct Network*, double *, double **, struct Wavelet **, double **, struct Data*, double, double);
char runName[100];
char outputDirectory[1000];
......
src/BayesWaveEvidence.c
View file @ 64c15c26
......@@ -2373,7 +2373,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
for(i=1; i<=model->signal[0]->size; i++)
{
ii = model->signal[0]->index[i];
intrinsic_fisher_update(model->signal[0]->intParams[ii], sigmas, model->SnGeo, 1.0, data->Tobs, data->NW,1 );
intrinsic_fisher_update(model->signal[0]->intParams[ii], sigmas, model->SnGeo, data->Tobs, data->NW,1 );
for(j=0; j<NW; j++) detC += log(sigmas[j]);
}
}
......@@ -2386,7 +2386,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
{
ii = model->glitch[ifo]->index[i];
//for(j=0; j<NW; j++)printf(" %i: %g\n",j,model->glitch[ifo]->intParams[ii][j]);
intrinsic_fisher_update(model->glitch[ifo]->intParams[ii], sigmas, model->Snf[ifo], 1.0, data->Tobs, data->NW, 1);
intrinsic_fisher_update(model->glitch[ifo]->intParams[ii], sigmas, model->Snf[ifo], data->Tobs, data->NW, 1);
//for(j=0; j<NW; j++)printf(" %i: %g\n",j,sigmas[j]);
for(j=0; j<NW; j++) detC += log(sigmas[j]);
}
......@@ -2428,7 +2428,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
for(i=0; i<100; i++)
{
draw_wavelet_params(params, prior->range, chain->seed, data->NW);
data->signal_amplitude_proposal(params, model->SnGeo, 1.0, chain->seed, data->Tobs, prior->range,prior->sSNRpeak);
data->signal_amplitude_proposal(params, model->SnGeo, chain->seed, data->Tobs, prior->range,prior->sSNRpeak);
if(params[3]<minAmp)minAmp=params[3];
if(params[3]>maxAmp)maxAmp=params[3];
}
......@@ -2455,7 +2455,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
for(i=0; i<100; i++)
{
draw_wavelet_params(params, prior->range, chain->seed, data->NW);
data->glitch_amplitude_proposal(params, model->Snf[ifo], 1.0, chain->seed, data->Tobs, prior->range, prior->gSNRpeak);
data->glitch_amplitude_proposal(params, model->Snf[ifo], chain->seed, data->Tobs, prior->range, prior->gSNRpeak);
if(params[3]<minAmp)minAmp=params[3];
if(params[3]>maxAmp)maxAmp=params[3];
}
......
src/BayesWaveIO.c
View file @ 64c15c26
......@@ -331,7 +331,7 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
{
data->glitchFlag = 1;
initialize_model(model, data, prior, psd, chain->seed);
for(ifo=0; ifo<NI; ifo++)model->Sn[ifo]=1.0;
//for(ifo=0; ifo<NI; ifo++)model->Sn[ifo]=1.0;
FILE **glitchParams = malloc(NI*sizeof(FILE *));
......@@ -367,7 +367,7 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
for(ifo=0; ifo<NI; ifo++)
{
GNR[ifo] = detector_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model->glitch[ifo]->templates, model->invSnf[ifo], model->Sn[ifo]);
GNR[ifo] = detector_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model->glitch[ifo]->templates, model->invSnf[ifo], 1.0);
fprintf(snrFile,"%g ",GNR[ifo]);
}
fprintf(snrFile,"\n");
......@@ -399,7 +399,7 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
{
data->signalFlag = 1;
initialize_model(model, data, prior, psd, chain->seed);
for(ifo=0; ifo<NI; ifo++)model->Sn[ifo]=1.0;
//for(ifo=0; ifo<NI; ifo++)model->Sn[ifo]=1.0;
FILE **signalParams = malloc(data->Npol*sizeof(FILE *));
......@@ -429,8 +429,8 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
double SNRnet;
double SNRifo[NI];
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]);
SNRnet = network_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model->response, model->invSnf, 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], 1.0);
if(SNRnet>100.) NCflag++;
if(SNRnet>200.) NCflag++;
sprintf(filename,"%ssnr/snr.dat",data->runName);
......@@ -835,11 +835,6 @@ void print_run_flags(FILE *fptr, struct Data *data, struct Prior *prior)
else fprintf(fptr, "DISABLED");
fprintf(fptr, "\n");
fprintf(fptr, " PSD fitting is ........... ");
if(data->psdFitFlag) fprintf(fptr, "ENABLED");
else fprintf(fptr, "DISABLED");
fprintf(fptr, "\n");
fprintf(fptr, " BayesLine PSD model is ... ");
if(data->bayesLineFlag) fprintf(fptr, "ENABLED");
else fprintf(fptr, "DISABLED");
......@@ -1013,7 +1008,7 @@ void print_model(FILE *fptr, struct Data *data, struct Chain *chain, struct Mode
fprintf(fptr,"%d %f ", chain->mc, model->logL-data->logLc+model->logLnorm);
fprintf(fptr,"%d ", model->signal[0]->size);
for(ifo=0; ifo<NI; ifo++) fprintf(fptr,"%d ", model->glitch[ifo]->size);
for(ifo=0; ifo<NI; ifo++) fprintf(fptr, "%f ",model->Sn[ifo]);
//for(ifo=0; ifo<NI; ifo++) fprintf(fptr, "%f ",model->Sn[ifo]);
if(data->stochasticFlag) fprintf(fptr,"%f %f ",model->background->logamp,model->background->index);
fprintf(fptr,"%f %f ",model->logL,model->logLnorm);
fprintf(fptr,"\n");
......@@ -1061,12 +1056,12 @@ void print_chain_status(struct Data *data, struct Chain *chain, struct Model **m
SNR = 0.0;
if(data->signalFlag) SNR = network_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model[chain->index[0]]->response, model[chain->index[0]]->invSnf, model[chain->index[0]]->Sn, data->NI);
if(data->signalFlag) SNR = network_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model[chain->index[0]]->response, model[chain->index[0]]->invSnf, data->NI);
fprintf(stdout, " logL=%f hSNR=%.0f ", model[chain->index[0]]->logL - data->logLc, SNR);
for(ifo=0; ifo<NI; ifo++)
{
SNR = 0.0;
if(data->glitchFlag) SNR = detector_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model[chain->index[0]]->glitch[ifo]->templates, model[chain->index[0]]->invSnf[ifo], model[chain->index[0]]->Sn[ifo]);
if(data->glitchFlag) SNR = detector_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model[chain->index[0]]->glitch[ifo]->templates, model[chain->index[0]]->invSnf[ifo], 1.0);
fprintf(stdout, " g%iSNR=%.0f ", ifo, SNR);
}
fprintf(stdout, "\n");
......@@ -1104,7 +1099,7 @@ void print_chain_status(struct Data *data, struct Chain *chain, struct Model **m
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_waveforms(char filename[], double *h, int N, double *Snf, double Tobs, int imin, int imax, double tmin, double tmax)
{
int i;
double x,t;
......@@ -1123,7 +1118,7 @@ void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf,
{
if(i>imin && i<imax)
{
x = sqrt(Snf[i]*eta);
x = sqrt(Snf[i]);
hf[i][0] = h[2*i]/x;
hf[i][1] = h[2*i+1]/x;
}
......@@ -1211,7 +1206,7 @@ void print_colored_time_domain_waveforms(char filename[], double *h, int N, doub
}
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_time_domain_hdot(char filename[], double *h, int N, double *Snf, double Tobs, int imin, int imax, double tmin, double tmax)
{
int i;
double x,t;
......@@ -1239,7 +1234,7 @@ void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, doub
{
if(i>imin && i<imax)
{
x = sqrt(Snf[i]*eta);
x = sqrt(Snf[i]);
// Need to switch real and imaginary parts [ie -i*(x+iy)=-ix+y]
......@@ -1277,7 +1272,7 @@ void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, doub
fclose(waveout);
}
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_waveforms(char filename[], double *h, int N, double *Snf, double Tobs, int imin, int imax)
{
int i;
double x,f;
......@@ -1288,7 +1283,7 @@ void print_frequency_domain_waveforms(char filename[], double *h, int N, double
for(i=imin; i<imax; i++)
{
x = sqrt(Snf[i]*eta);
x = sqrt(Snf[i]);
hf[2*i] /= x;
hf[2*i+1] /= x;
}
......@@ -1568,10 +1563,6 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
if(ppt) sprintf(data->outputDirectory,"%s/",ppt->value);
else sprintf(data->outputDirectory,"./");
ppt = LALInferenceGetProcParamVal(commandLine, "--noPSDfit");
if(ppt) data->psdFitFlag = 0;
else data->psdFitFlag = 1;
ppt = LALInferenceGetProcParamVal(commandLine, "--bayesLine");
if(ppt) data->bayesLineFlag = 1;
else data->bayesLineFlag = 0;
......@@ -1685,26 +1676,6 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
fprintf(stdout,"*********************************************************\n");
data->splineEvidenceFlag = 0;
}
if(data->psdFitFlag && !data->noiseSimFlag)
{
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;
}
/*if(data->waveletPriorFlag)
{
fprintf(stdout,"************************ WARNING ************************\n");
fprintf(stdout,"waveletPrior flag overrides --Dmax \n");
fprintf(stdout,"Setting Dmax = 100 \n");
fprintf(stdout,"Remove --waveletPrior if you really want to cap N \n");
fprintf(stdout,"*********************************************************\n");
data->Dmax = 100;
}*/
}
void parse_bayesline_parameters(struct Data *data, struct Model *model, FILE **splinechain, FILE **linechain, double **psd)
......@@ -1908,7 +1879,7 @@ void dump_time_domain_waveform(struct Data *data, struct Model *model, struct Pr
for(ifo=0; ifo<data->NI; ifo++)
{
sprintf(filename,"%s.%i",root,ifo);
print_time_domain_waveforms(filename, model->response[ifo], data->N, model->Snf[ifo], 1, data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
print_time_domain_waveforms(filename, model->response[ifo], data->N, model->Snf[ifo], data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
}
}
......
src/BayesWaveIO.h
View file @ 64c15c26
......@@ -61,9 +61,9 @@ void print_model(FILE *fptr, struct Data *data, struct Chain *chain, struct Mode
void print_signal_model(FILE *fptr, struct Model *model, int n);
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);
void print_frequency_domain_waveforms(char filename[], double *h, int N, double *Snf, double Tobs, double eta, int imin, int imax);
void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf, double Tobs, int imin, int imax, double tmin, double tmax);
void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, 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, int imin, int imax);
void print_frequency_domain_data(char filename[], double *h, int N, double Tobs, int imin, int imax);
void print_colored_time_domain_waveforms(char filename[], double *h, int N, double Tobs, int imin, int imax, double tmin, double tmax);
......
src/BayesWaveLikelihood.c
View file @ 64c15c26
......@@ -109,12 +109,12 @@ void phase_blind_time_shift(double *corr, double *corrf, double *data1, double *
}
double loglike(int imin, int imax, double *r, double Sn, double *invSnf)
double loglike(int imin, int imax, double *r, double *invSnf)
{
/*
logL = -(r|r)/2 - Nlog(Sn)
logL = -(r|r)/2
*/
return -0.5*(fourier_nwip(imin, imax, r, r, invSnf)/Sn) - (double)(imax-imin)*log(Sn);
return -0.5*fourier_nwip(imin, imax, r, r, invSnf);
}
double loglike_stochastic(int NI, int imin, int imax, double **r, double ***invCij, double *detC)
......@@ -206,14 +206,14 @@ void recompute_residual(struct Data *data, struct Model **model, struct Chain *c
if(data->stochasticFlag)
{
//TODO: No support for glitch model in stochastic likelihood
ComputeNoiseCorrelationMatrix(data, model_x->Snf, model_x->Sn, model_x->background);
ComputeNoiseCorrelationMatrix(data, model_x->Snf, model_x->background);
model_x->logL = loglike_stochastic(data->NI, data->imin, data->imax, data->r, model_x->background->Cij, model_x->background->detCij);
}
else
{
for(ifo=0; ifo<NI; ifo++)
{
model_x->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], model_x->Sn[ifo], model_x->invSnf[ifo]);
model_x->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], model_x->invSnf[ifo]);
model_x->logL += model_x->detLogL[ifo];
for(i=0; i<data->N/2; i++)
{
......@@ -306,7 +306,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
}
//proposed model
double *Sny = my->Sn;
//double *Sny = my->Sn;
double **ry = my->response;
double **hy = my->h;
......@@ -588,7 +588,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
if(data->stochasticFlag)
{
//TODO: No support for glitch model in stochastic likelihood
ComputeNoiseCorrelationMatrix(data, mx->Snf, my->Sn, my->background);
ComputeNoiseCorrelationMatrix(data, mx->Snf, my->background);
logLy = loglike_stochastic(NI, imin, imax, r, my->background->Cij, my->background->detCij);
}
......@@ -607,7 +607,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
*/
for(ifo=0; ifo<NI; ifo++)
{
my->detLogL[ifo] = (mx->detLogL[ifo] + (imax-imin)*log(mx->Sn[ifo]))*(mx->Sn[ifo]/Sny[ifo]) -(imax-imin)*log(Sny[ifo]);
my->detLogL[ifo] = mx->detLogL[ifo];
logLy += my->detLogL[ifo];
}
......@@ -645,7 +645,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
for(ifo=0; ifo<NI; ifo++)
{
//remove contribution between imin and imax of old model
logLseg = loglike(jmin, jmax, r[ifo], Sny[ifo], invSnf[ifo]);
logLseg = loglike(jmin, jmax, r[ifo], invSnf[ifo]);
my->detLogL[ifo] -= logLseg;
logLy -= logLseg;
......@@ -667,7 +667,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
r[ifo][k] -= ry[ifo][k];
}
}
logLseg = loglike(jmin, jmax, r[ifo], Sny[ifo], invSnf[ifo]);
logLseg = loglike(jmin, jmax, r[ifo], invSnf[ifo]);
my->detLogL[ifo] += logLseg;
logLy += logLseg;
}
......@@ -675,7 +675,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
else
{
//remove contribution between imin and imax of old model
logLseg = loglike(jmin, jmax, r[det], Sny[det], invSnf[det]);
logLseg = loglike(jmin, jmax, r[det], invSnf[det]);
my->detLogL[det] -= logLseg;
logLy -= logLseg;
......@@ -687,7 +687,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
r[det][j] = s[det][j] - ry[det][j] - gy[det][j];
r[det][k] = s[det][k] - ry[det][k] - gy[det][k];
}
logLseg = loglike(jmin, jmax, r[det], Sny[det], invSnf[det]);
logLseg = loglike(jmin, jmax, r[det], invSnf[det]);
my->detLogL[det] += logLseg;
logLy += logLseg;
}
......@@ -716,7 +716,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
my->detLogL[ifo] = mx->detLogL[ifo];
//remove contribution between imin and imax of old model
logLseg = loglike(jmin, jmax, r[ifo], mx->Sn[ifo], invSnf[ifo]);
logLseg = loglike(jmin, jmax, r[ifo], invSnf[ifo]);
my->detLogL[ifo] -= logLseg;
logLy -= logLseg;
......@@ -738,7 +738,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
r[ifo][k] -= ry[ifo][k];
}
}
logLseg = loglike(jmin, jmax, r[ifo], Sny[ifo], invSnf[ifo]);
logLseg = loglike(jmin, jmax, r[ifo], invSnf[ifo]);
my->detLogL[ifo] += logLseg;
logLy += logLseg;
}
......@@ -748,7 +748,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
my->detLogL[det] = mx->detLogL[det];
//remove contribution between imin and imax of old model
logLseg = loglike(jmin, jmax, r[det], mx->Sn[det], invSnf[det]);
logLseg = loglike(jmin, jmax, r[det], invSnf[det]);
my->detLogL[det] -= logLseg;
logLy -= logLseg;
......@@ -760,7 +760,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
r[det][j] = s[det][j] - ry[det][j] - gy[det][j];
r[det][k] = s[det][k] - ry[det][k] - gy[det][k];
}
logLseg = loglike(jmin, jmax, r[det], Sny[det], invSnf[det]);
logLseg = loglike(jmin, jmax, r[det], invSnf[det]);
my->detLogL[det] += logLseg;
logLy += logLseg;
}
......@@ -782,7 +782,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
r[ifo][i] -= ry[ifo][i];
}
}
my->detLogL[ifo] = loglike(imin, imax, r[ifo], Sny[ifo], invSnf[ifo]);
my->detLogL[ifo] = loglike(imin, imax, r[ifo], invSnf[ifo]);
logLy += my->detLogL[ifo];
}
}
......@@ -794,14 +794,14 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
}
double EvaluateExtrinsicConstantLogLikelihood(UNUSED struct Network *projection, double *params, UNUSED double **invSnf, UNUSED double *Sn, UNUSED struct Wavelet **geo, UNUSED double **g, UNUSED struct Data *data, UNUSED double fmin, UNUSED double fmax)
double EvaluateExtrinsicConstantLogLikelihood(UNUSED struct Network *projection, double *params, UNUSED double **invSnf, UNUSED struct Wavelet **geo, UNUSED double **g, UNUSED struct Data *data, UNUSED double fmin, UNUSED double fmax)
{
if(extrinsic_checkrange(params)) return -1.0e60;
else return 0.0;
}
double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, double *params, double **invSnf, double *Sn, struct Wavelet **geo, double **g, struct Data *data, double fmin, double fmax)
double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, double *params, double **invSnf, struct Wavelet **geo, double **g, struct Data *data, double fmin, double fmax)
{
int i, n;
int NI,NP,N;
......@@ -852,7 +852,7 @@ double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, doubl
if(data->signalFlag) r[ifo][n] -= h[ifo][n];
if(data->glitchFlag) r[ifo][n] -= g[ifo][n];
}
logL += loglike(imin, imax, r[ifo], Sn[ifo], invSnf[ifo]);
logL += loglike(imin, imax, r[ifo], invSnf[ifo]);
}
......@@ -863,7 +863,7 @@ double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, doubl
return logL;
}
double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *params, double **invSnf, double *Sn, struct Wavelet **geo, double **g, struct Data *data, double fmin, double fmax)
double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *params, double **invSnf, UNUSED double *Sn, struct Wavelet **geo, double **g, struct Data *data, double fmin, double fmax)
{
int i, n;
int NI,N;
......@@ -945,7 +945,7 @@ double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *
}
}
logL = 0.0;
for(ifo=0; ifo<NI; ifo++) logL += loglike(imin, imax, r[ifo], Sn[ifo], invSnf[ifo]);
for(ifo=0; ifo<NI; ifo++) logL += loglike(imin, imax, r[ifo], invSnf[ifo]);
free_double_matrix(h,NI-1);
free_double_matrix(r,NI-1);
......
src/BayesWaveLikelihood.h
View file @ 64c15c26
......@@ -29,7 +29,7 @@
/* */
/* ********************************************************************************** */
double loglike(int imin, int imax, double *r, double Sn, double *Snf);
double loglike(int imin, int imax, double *r, double *Snf);
double loglike_stochastic(int NI, int imin, int imax, double **r, double ***invCij, double *detC);
void recompute_residual(struct Data *data, struct Model **model, struct Chain *chain);
......@@ -37,9 +37,9 @@ void recompute_residual(struct Data *data, struct Model **model, struct Chain *c
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, struct Wavelet **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 struct Wavelet **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, struct Wavelet **geo, double **g, struct Data *data, double fmin, double fmax);
double EvaluateExtrinsicSearchLogLikelihood (struct Network *projection, double *params, double **invSnf, UNUSED double *Sn, struct Wavelet **geo, double **g, struct Data *data, double fmin, double fmax);
double EvaluateExtrinsicConstantLogLikelihood (UNUSED struct Network *projection, double *params, UNUSED double **invSnf, UNUSED struct Wavelet **geo, UNUSED double **g, UNUSED struct Data *data, UNUSED double fmin, UNUSED double fmax);
double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, double *params, double **invSnf, struct Wavelet **geo, double **g, struct Data *data, double fmin, double fmax);
......
src/BayesWaveMCMC.c
View file @ 64c15c26
......@@ -169,7 +169,7 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
for(j=1; j<=glitch[ifo]->size; j++)
{
draw_wavelet_params(glitch[ifo]->intParams[j], prior->range, chain->seed, data->NW);
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);
if(data->amplitudePriorFlag) data->glitch_amplitude_proposal(glitch[ifo]->intParams[j], model[ic]->Snf[ifo], chain->seed, data->Tobs, prior->range, prior->gSNRpeak);
model[ic]->wavelet(glitch[ifo]->templates, glitch[ifo]->intParams[j], data->N, 1, Tobs);
}
......@@ -205,7 +205,7 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
{
draw_wavelet_params(signal[n]->intParams[j], prior->range, chain->seed, data->NW);
if(data->amplitudePriorFlag)
data->signal_amplitude_proposal(signal[n]->intParams[j], model[ic]->SnGeo, 1.0, chain->seed, data->Tobs, prior->range, prior->sSNRpeak);
data->signal_amplitude_proposal(signal[n]->intParams[j], model[ic]->SnGeo, chain->seed, data->Tobs, prior->range, prior->sSNRpeak);
model[ic]->wavelet(signal[n]->templates, signal[n]->intParams[j], data->N, 1, Tobs);
}
......@@ -215,15 +215,12 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
waveformProject(data, model[ic]->projection, model[ic]->extParams, model[ic]->response, model[ic]->h, data->fmin, data->fmax);
}
// initialize the noise parameters
for(ifo=0; ifo<NI; ifo++) model[ic]->Sn[ifo] = 1.0;
// form up residual and compute initial likelihood
if(data->stochasticFlag)
{
//TODO: No support for glitch model in stochastic likelihood
ComputeNoiseCorrelationMatrix(data, model[ic]->Snf, model[ic]->Sn, model[ic]->background);
ComputeNoiseCorrelationMatrix(data, model[ic]->Snf, model[ic]->background);
model[ic]->logL = loglike_stochastic(data->NI, data->imin, data->imax, data->r, model[ic]->background->Cij, model[ic]->background->detCij);
}
......@@ -244,7 +241,7 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
}
if(!data->constantLogLFlag)
{
model[ic]->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], model[ic]->Sn[ifo], model[ic]->invSnf[ifo]);
model[ic]->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], model[ic]->invSnf[ifo]);
model[ic]->logL += model[ic]->detLogL[ifo];
for(i=0; i<data->N/2; i++)
{
......@@ -469,7 +466,6 @@ static void resume_sampler(struct Data *data, struct Chain *chain, struct Model
fscanf(fptr,"%i %lg",&i,&logL); //iteration and logL
fscanf(fptr,"%i",&NS); //signal model dimension
for(i=0; i<NI; i++)fscanf(fptr,"%i",&NG[i]); //glitch model dimension
for(i=0; i<NI; i++)fscanf(fptr,"%lg",&model[ic]->Sn[i]); //PSD scale parameters
fclose(fptr);
signal = model[ic]->signal;
......@@ -628,7 +624,7 @@ static void resume_sampler(struct Data *data, struct Chain *chain, struct Model
}
if(!data->constantLogLFlag)
{
model[ic]->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], model[ic]->Sn[ifo], model[ic]->invSnf[ifo]);
model[ic]->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], model[ic]->invSnf[ifo]);
model[ic]->logL += model[ic]->detLogL[ifo];
for(i=0; i<data->N/2; i++)
{
......@@ -818,7 +814,6 @@ static void save_sampler(struct Data *data, struct Chain *chain, struct Model **
fprintf(fptr,"%i %.12g ",i,model[ic]->logL); //iteration and logL
fprintf(fptr,"%i ",signal[0]->size); //signal model dimension
for(i=0; i<NI; i++)fprintf(fptr,"%i ",glitch[i]->size); //glitch model dimension
for(i=0; i<NI; i++)fprintf(fptr,"%lg ",model[ic]->Sn[i]); //PSD scale parameters
fprintf(fptr,"\n");
fclose(fptr);
......@@ -1210,14 +1205,14 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
for(ifo=0; ifo<NI; ifo++)
{
sprintf(filename,"waveforms/%s_data_%s_%d.dat", modelname, data->ifos[ifo], frame);
print_time_domain_waveforms(filename, data->s[ifo], N, model[chain->index[0]]->Snf[ifo], model[chain->index[0]]->Sn[ifo],data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
print_time_domain_waveforms(filename, data->s[ifo], N, model[chain->index[0]]->Snf[ifo], data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
if(data->bayesLineFlag)
{
for(i=0; i<N; i++) data->r[ifo][i] = data->s[ifo][i] - model[chain->index[0]]->glitch[ifo]->templates[i];