diff --git a/src/BayesWave.c b/src/BayesWave.c
index ff818bfd53297713ed90c55432902ee9a6848e7b..44fe59fdbce986bf29ff6d1bfc5794c7c172f2ac 100644
--- a/src/BayesWave.c
+++ b/src/BayesWave.c
@@ -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);
diff --git a/src/BayesWave.h b/src/BayesWave.h
index 349d8b0640fd0feebc6878cf460f465e5e012de0..fb585e99d868e038afe03bdeef5ebf46f7c9fc45 100644
--- a/src/BayesWave.h
+++ b/src/BayesWave.h
@@ -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];
diff --git a/src/BayesWaveEvidence.c b/src/BayesWaveEvidence.c
index 5086d98d21db98af55b0047c9832941697d90d57..bd0a10d048ef5ea41bfa6005fc2fccc9c3fae2d0 100644
--- a/src/BayesWaveEvidence.c
+++ b/src/BayesWaveEvidence.c
@@ -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];
}
diff --git a/src/BayesWaveIO.c b/src/BayesWaveIO.c
index a5e8a1751dfc956ad7621e830ba61a124e4d7e9b..1f95c4448c14beee22d12b29cf7e55bb8fbeff9a 100644
--- a/src/BayesWaveIO.c
+++ b/src/BayesWaveIO.c
@@ -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]);
}
}
diff --git a/src/BayesWaveIO.h b/src/BayesWaveIO.h
index 83805df4516a4123fd23b0ff03af0e33f0a3b945..9dad20fe2dbcab8085787eba7959ee1c0a508d7e 100644
--- a/src/BayesWaveIO.h
+++ b/src/BayesWaveIO.h
@@ -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);
diff --git a/src/BayesWaveLikelihood.c b/src/BayesWaveLikelihood.c
index d972ccffe9104b7303da4932255dbb522eb5bd37..773a47ee82fd9f972f89174859841e783a8124b8 100644
--- a/src/BayesWaveLikelihood.c
+++ b/src/BayesWaveLikelihood.c
@@ -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);
diff --git a/src/BayesWaveLikelihood.h b/src/BayesWaveLikelihood.h
index 7389d3cf72c48c20b0018850e1c92c2092bf500c..3b9c3137919f23ed3274ef406f5c1c5716486624 100644
--- a/src/BayesWaveLikelihood.h
+++ b/src/BayesWaveLikelihood.h
@@ -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);
diff --git a/src/BayesWaveMCMC.c b/src/BayesWaveMCMC.c
index 6227dea7bc113ade72442c704ffa7c814166f7c2..e342d81596d761589171cba808db18cff2803469 100644
--- a/src/BayesWaveMCMC.c
+++ b/src/BayesWaveMCMC.c
@@ -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];
sprintf(filename,"waveforms/%s_frequency_residual_%s_%d.dat",modelname,data->ifos[ifo],frame);
- print_frequency_domain_waveforms(filename, data->r[ifo], N, model[chain->index[0]]->Snf[ifo], data->Tobs, model[chain->index[0]]->Sn[ifo], data->imin, data->imax);
+ print_frequency_domain_waveforms(filename, data->r[ifo], N, model[chain->index[0]]->Snf[ifo], data->Tobs, data->imin, data->imax);
sprintf(filename,"waveforms/%s_psd_%s_%d.dat", modelname, data->ifos[ifo], frame);
print_bayesline_spectrum(filename, bayesline[chain->index[0]][ifo]->freq, bayesline[chain->index[0]][ifo]->spow, bayesline[chain->index[0]][ifo]->Sbase, bayesline[chain->index[0]][ifo]->Sline, bayesline[chain->index[0]][ifo]->data->ncut);
@@ -1226,7 +1221,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
if(data->glitchFlag)
{
sprintf(filename,"waveforms/%s_glitch_%s_%d.dat", modelname, data->ifos[ifo], frame);
- print_time_domain_waveforms(filename, model[chain->index[0]]->glitch[ifo]->templates, 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, model[chain->index[0]]->glitch[ifo]->templates, N, model[chain->index[0]]->Snf[ifo], data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
sprintf(filename,"waveforms/%s_colored_glitch_%s_%d.dat", modelname, data->ifos[ifo], frame);
print_colored_time_domain_waveforms(filename, model[chain->index[0]]->glitch[ifo]->templates, N,data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
@@ -1235,7 +1230,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
if(data->signalFlag)
{
sprintf(filename,"waveforms/%s_wave_%s_%d.dat", modelname, data->ifos[ifo], frame);
- print_time_domain_waveforms(filename, model[chain->index[0]]->response[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, model[chain->index[0]]->response[ifo], N, model[chain->index[0]]->Snf[ifo], data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
sprintf(filename,"waveforms/%s_colored_wave_%s_%d.dat", modelname, data->ifos[ifo], frame);
print_colored_time_domain_waveforms(filename, model[chain->index[0]]->response[ifo], N, data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
@@ -1485,7 +1480,7 @@ void EvolveBayesLineParameters(struct Data *data, struct Model **model, struct B
if(data->glitchFlag) data->r[ifo][i] -= model_x->glitch[ifo]->templates[i];
}
- 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]);
if(!data->constantLogLFlag)
{
model_x->logL += model_x->detLogL[ifo];
@@ -1635,8 +1630,6 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
/* PRIOR */
int dmax = 1;
- double Snmin = prior->Snmin;
- double Snmax = prior->Snmax;
double **range = prior->range;
/* MODEL */
@@ -1645,9 +1638,6 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
model_x->size=model_x->signal[0]->size;
for(ifo=0; ifo<NI; ifo++) model_x->size+=model_x->glitch[ifo]->size;
-
- double *Snx = model_x->Sn;
-
struct Wavelet **wave_x;
int *larrayx;
double **paramsx;
@@ -1664,8 +1654,6 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
model_y->size=model_x->size;
- double *Sny = model_y->Sn;
-
struct Wavelet **wave_y;
int *larrayy;
double **paramsy;
@@ -1687,8 +1675,6 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
//acc gets set to 1 if model_y gets accepted (no need to copy at next iteration)
if(acc==0 && mc>1) copy_int_model(model_x, model_y, N, data->NI, det);
- if(data->psdFitFlag) for(ifo=0; ifo<NI; ifo++) Sny[ifo] = Snx[ifo];
- else for(ifo=0; ifo<NI; ifo++) Sny[ifo] = 1.0;
//reset acceptence counter
acc=0;
@@ -1968,12 +1954,12 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
larrayx = wave_x[n]->index;
larrayy = wave_y[n]->index;
- data->signal_amplitude_proposal(paramsy[ii], model_x->SnGeo, 1.0, seed, data->Tobs, prior->range, prior->sSNRpeak);
+ data->signal_amplitude_proposal(paramsy[ii], model_x->SnGeo, seed, data->Tobs, prior->range, prior->sSNRpeak);
logqy += ( data->signal_amplitude_prior(paramsy[ii], model_x->SnGeo, data->Tobs, prior->sSNRpeak) );
}
else
{
- data->glitch_amplitude_proposal(paramsy[ii], model_x->Snf[det], Sny[det], seed, data->Tobs, prior->range, prior->gSNRpeak);
+ data->glitch_amplitude_proposal(paramsy[ii], model_x->Snf[det], seed, data->Tobs, prior->range, prior->gSNRpeak);
logqy += (data->glitch_amplitude_prior(paramsy[ii], model_x->Snf[det], data->Tobs, prior->gSNRpeak) );
}
}
@@ -2000,11 +1986,6 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
{
typ = 1;
- for(ifo=0; ifo<NI; ifo++)
- {
- if(data->psdFitFlag) Sny[ifo] = Snx[ifo]+gaussian_draw(seed)/sqrt(0.5*(double)N);
- if(Sny[ifo] > Snmax || Sny[ifo] < Snmin) test = 1;
- }
//update stochastic background parameters
if(data->stochasticFlag) stochastic_background_proposal(model_x->background, model_y->background, seed, &test);
@@ -2061,14 +2042,14 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
{
if(det==-1)
{
- data->signal_amplitude_proposal(paramsy[ii], model_x->SnGeo, 1.0, seed, data->Tobs, prior->range, prior->sSNRpeak);
+ data->signal_amplitude_proposal(paramsy[ii], model_x->SnGeo, seed, data->Tobs, prior->range, prior->sSNRpeak);
logqy += (data->signal_amplitude_prior(paramsy[ii], model_x->SnGeo, data->Tobs, prior->sSNRpeak));
logqx += (data->signal_amplitude_prior(paramsx[ii], model_x->SnGeo, data->Tobs, prior->sSNRpeak));
}
else
{
- data->glitch_amplitude_proposal(paramsy[ii], model_x->Snf[det], Sny[det], seed, data->Tobs, prior->range, prior->gSNRpeak);
+ data->glitch_amplitude_proposal(paramsy[ii], model_x->Snf[det], seed, data->Tobs, prior->range, prior->gSNRpeak);
logqy += (data->glitch_amplitude_prior(paramsy[ii], model_x->Snf[det], data->Tobs, prior->gSNRpeak));
logqx += (data->glitch_amplitude_prior(paramsx[ii], model_x->Snf[det], data->Tobs, prior->gSNRpeak));
}
@@ -2130,12 +2111,12 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if(det==-1)
{
- intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->SnGeo, 1.0, &logqx, &logqy, seed, data->Tobs, &test, wave_y[n]->dimension, tfqFlag);
+ intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->SnGeo, &logqx, &logqy, seed, data->Tobs, &test, wave_y[n]->dimension, tfqFlag);
}
else
{
- intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->Snf[det], 1.0, &logqx, &logqy, seed, data->Tobs, &test, wave_y[n]->dimension, tfqFlag);
+ intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->Snf[det], &logqx, &logqy, seed, data->Tobs, &test, wave_y[n]->dimension, tfqFlag);
}
}
if(det==-1 && model_y->Npol > 1) constrain_hplus_hcross(wave_y, ii);
@@ -2157,12 +2138,12 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
{
if(det==-1)
{
- data->signal_amplitude_proposal(paramsy[ii], model_x->SnGeo, 1.0, seed, data->Tobs, prior->range, prior->sSNRpeak);
+ data->signal_amplitude_proposal(paramsy[ii], model_x->SnGeo, seed, data->Tobs, prior->range, prior->sSNRpeak);
logqy += ( data->signal_amplitude_prior(paramsy[ii], model_x->SnGeo, data->Tobs, prior->sSNRpeak) );
}
else
{
- data->glitch_amplitude_proposal(paramsy[ii], model_x->Snf[det], Sny[n], seed, data->Tobs, prior->range, prior->gSNRpeak);
+ data->glitch_amplitude_proposal(paramsy[ii], model_x->Snf[det], seed, data->Tobs, prior->range, prior->gSNRpeak);
logqy += ( data->glitch_amplitude_prior(paramsy[ii], model_x->Snf[n], data->Tobs, prior->gSNRpeak) );
}
}
@@ -2362,7 +2343,6 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
double logLx, logLy, logH;
double alpha;
double **glitch;
- double *Sn;
double *SnGeox;
double *SnGeoy;
double *paramsx, *paramsy;
@@ -2379,7 +2359,6 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
// extrinsic parameters
struct Model *model_x = model[chain->index[ic]];
struct Wavelet **smodel = model_x->signal;
- Sn = model_x->Sn;
paramsx = model_x->extParams;
SnGeox = model_x->SnGeo;
@@ -2426,7 +2405,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
double fmin = data->fmin;
double fmax = data->fmax;
- logLx = data->extrinsic_likelihood(model_x->projection, paramsx, model_x->invSnf, Sn, smodel, glitch, data, fmin, fmax);
+ logLx = data->extrinsic_likelihood(model_x->projection, paramsx, model_x->invSnf, smodel, glitch, data, fmin, fmax);
//Compute Fisher Matrix for current sky location
struct FisherMatrix *fisher = model_x->fisher;
@@ -2505,7 +2484,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if(test==0)
{
- logLy = data->extrinsic_likelihood(model_x->projection, paramsy, model_x->invSnf, Sn, smodel, glitch, data, fmin, fmax);
+ logLy = data->extrinsic_likelihood(model_x->projection, paramsy, model_x->invSnf, smodel, glitch, data, fmin, fmax);
//amplitude prior
if(data->amplitudePriorFlag)
diff --git a/src/BayesWaveMath.c b/src/BayesWaveMath.c
index e4c2b926279917b99f6159efb6259648a1224673..5fb4399fd3d6f691467f545c814f572e0089e9dd 100644
--- a/src/BayesWaveMath.c
+++ b/src/BayesWaveMath.c
@@ -86,18 +86,18 @@ double fourier_nwip(int imin, int imax, double *a, double *b, double *invSn)
return(2.0*arg);
}
-double network_nwip(int imin, int imax, double **a, double **b, double **invSn, double *eta, int NI)
+double network_nwip(int imin, int imax, double **a, double **b, double **invSn, int NI)
{
int i;
double sum=0.0;
- for(i=0; i<NI; i++) sum += fourier_nwip(imin,imax,a[i],b[i],invSn[i])/eta[i];
+ for(i=0; i<NI; i++) sum += fourier_nwip(imin,imax,a[i],b[i],invSn[i]);
return(sum);
}
-double network_snr(int imin, int imax, double **h, double **invpsd, double *eta, int NI )
+double network_snr(int imin, int imax, double **h, double **invpsd, int NI )
{
- return sqrt( network_nwip(imin,imax,h,h,invpsd,eta,NI) );
+ return sqrt( network_nwip(imin,imax,h,h,invpsd,NI) );
}
double detector_snr(int imin, int imax, double *g, double *invpsd, double eta)
diff --git a/src/BayesWaveMath.h b/src/BayesWaveMath.h
index e9a31fbb7454f0567065cfb1ba87fd098fb6de58..eb870b55d2694b16953fd7f0ad730cc13e4fbb7e 100644
--- a/src/BayesWaveMath.h
+++ b/src/BayesWaveMath.h
@@ -27,8 +27,8 @@ double gaussian_norm(double min, double max, double sigma);
void recursive_phase_evolution(double dre, double dim, double *cosPhase, double *sinPhase);
-double network_nwip(int imin, int imax, double **a, double **b, double **invSn, double *eta, int NI);
-double network_snr(int imin, int imax, double **h, double **invpsd, double *eta, int NI );
+double network_nwip(int imin, int imax, double **a, double **b, double **invSn, int NI);
+double network_snr(int imin, int imax, double **h, double **invpsd, int NI );
double detector_snr(int imin, int imax, double *g, double *invpsd, double eta);
double fourier_nwip(int imin, int imax, double *a, double *b, double *invSn);
diff --git a/src/BayesWaveModel.c b/src/BayesWaveModel.c
index a3a2b008d2a2f9f34acfd739bbc6d895722d42e5..92deb6dc57e432433f8d5eeebc98072550379c55 100644
--- a/src/BayesWaveModel.c
+++ b/src/BayesWaveModel.c
@@ -137,7 +137,7 @@ void OverlapReductionFunction(struct Data *data)
}
}
-void ComputeNoiseCorrelationMatrix(struct Data *data, double **Snf, double *Sn, struct Background *background)
+void ComputeNoiseCorrelationMatrix(struct Data *data, double **Snf, struct Background *background)//void ComputeNoiseCorrelationMatrix(struct Data *data, double **Snf, double *Sn, struct Background *background)
{
int i;
@@ -162,8 +162,8 @@ void ComputeNoiseCorrelationMatrix(struct Data *data, double **Snf, double *Sn,
Calculate noise correlation matrix
Combining instrument noise and stochastic backgrouns signal
*/
- C11 = background->spectrum[i] + Snf[0][i]*Sn[0];
- C22 = background->spectrum[i] + Snf[1][i]*Sn[1];
+ C11 = background->spectrum[i] + Snf[0][i];
+ C22 = background->spectrum[i] + Snf[1][i];
C12 = background->spectrum[i] * data->ORF[i];
C21 = C12;
@@ -175,9 +175,6 @@ void ComputeNoiseCorrelationMatrix(struct Data *data, double **Snf, double *Sn,
background->Cij[1][0][i] = -C21*invC;
background->Cij[1][1][i] = C11*invC;
-// printf("%lg %lg\n%lg %lg\n",background->Cij[0][0][i],background->Cij[0][1][i],background->Cij[1][0][i],background->Cij[1][1][i]);
-// printf("================\n");
-// printf("%lg %lg\n",Snf[0][i]*Sn[0],data->Snf[1][i]*Sn[1]);
}
}
@@ -708,7 +705,7 @@ void initialize_model(struct Model *model, struct Data *data, struct Prior *prio
//gaussian noise model
- model->Sn = double_vector(NI-1);
+ //model->Sn = double_vector(NI-1);
model->SnS = double_matrix(NI-1,halfN-1);
model->Snf = double_matrix(NI-1,halfN-1);
model->invSnf= double_matrix(NI-1,halfN-1);
@@ -1010,7 +1007,7 @@ void copy_int_model(struct Model *origin, struct Model *copy, int N, int NI, int
{
copy->detLogL[ifo] = origin->detLogL[ifo];
- copy->Sn[ifo] = origin->Sn[ifo]; //noise PSD scale
+ //copy->Sn[ifo] = origin->Sn[ifo]; //noise PSD scale
}
}
@@ -1054,7 +1051,7 @@ void free_model(struct Model *model, struct Data *data, struct Prior *prior)
int NI = data->NI;
int Npol = model->Npol;
- free_double_vector(model->Sn);
+ //free_double_vector(model->Sn);
free_double_vector(model->SnGeo);
free_double_matrix(model->SnS,NI-1);
free_double_matrix(model->Snf,NI-1);
diff --git a/src/BayesWaveModel.h b/src/BayesWaveModel.h
index 38cceafaa47aa6161cdbc77ae56cee620e323b89..53ce57d771efdb000b16b32c5af9252e5ebfe138 100644
--- a/src/BayesWaveModel.h
+++ b/src/BayesWaveModel.h
@@ -33,7 +33,7 @@ void Shf_Geocenter(struct Data *data, struct Model *model, double *SnGeo, double
void Shf_Geocenter_full(struct Data *data, struct Network *projection, double **Snf, double *SnGeo, double *params);
void OverlapReductionFunction(struct Data *data);
-void ComputeNoiseCorrelationMatrix(struct Data *data, double **Snf, double *Sn, struct Background *background);
+void ComputeNoiseCorrelationMatrix(struct Data *data, double **Snf, struct Background *background);//void ComputeNoiseCorrelationMatrix(struct Data *data, double **Snf, double *Sn, struct Background *background);
double symmetric_snr_ratio(struct Data *data, struct Network *projection, double *params);
diff --git a/src/BayesWaveProposal.c b/src/BayesWaveProposal.c
index 0934c5c9d2e6f71a27da87c81beec8566015efba..94087207ece9bf5c840eb0b05956d3f239870c95 100644
--- a/src/BayesWaveProposal.c
+++ b/src/BayesWaveProposal.c
@@ -52,7 +52,7 @@
/* */
/* ********************************************************************************** */
-void draw_glitch_amplitude(double *params, double *Snf, double Sa, gsl_rng *seed, double Tobs, double **range, double SNRpeak)
+void draw_glitch_amplitude(double *params, double *Snf, gsl_rng *seed, double Tobs, double **range, double SNRpeak)
{
int i, k;
double SNR, den=-1.0, alpha=1.0;
@@ -61,8 +61,6 @@ void draw_glitch_amplitude(double *params, double *Snf, double Sa, gsl_rng *seed
double dfac, dfac3;
- Sa = 1.0;
-
double SNR2 = 2.0*SNRpeak;
double SNRsq = 2.0*SNRpeak*SNRpeak;
@@ -72,8 +70,8 @@ void draw_glitch_amplitude(double *params, double *Snf, double Sa, gsl_rng *seed
invmax = 1./max;
i = (int)(params[1]*Tobs);
- double SNRmin = range[3][0]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
- double SNRmax = range[3][1]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
+ double SNRmin = range[3][0]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snf[i]*2.0/Tobs));
+ double SNRmax = range[3][1]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snf[i]*2.0/Tobs));
k = 0;
@@ -112,11 +110,11 @@ void draw_glitch_amplitude(double *params, double *Snf, double Sa, gsl_rng *seed
}
//SNR defined with Sn(f) but Snf array holdes <n_i^2>
- params[3] = SNR/sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
+ params[3] = SNR/sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snf[i]*2.0/Tobs));
}
-void draw_signal_amplitude(double *params, double *Snf, double Sa, gsl_rng *seed, double Tobs, double **range, double SNRpeak)
+void draw_signal_amplitude(double *params, double *Snf, gsl_rng *seed, double Tobs, double **range, double SNRpeak)
{
int i, k;
double SNR, den=-1.0, alpha=1.0;
@@ -125,8 +123,6 @@ void draw_signal_amplitude(double *params, double *Snf, double Sa, gsl_rng *seed
double dfac, dfac5;
- Sa = 1.0;
-
double SNR4 = 4.0*SNRpeak;
double SNRsq = 4.0*SNRpeak*SNRpeak;
@@ -137,8 +133,8 @@ void draw_signal_amplitude(double *params, double *Snf, double Sa, gsl_rng *seed
i = (int)(params[1]*Tobs);
- double SNRmin = range[3][0]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
- double SNRmax = range[3][1]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
+ double SNRmin = range[3][0]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snf[i]*2.0/Tobs));
+ double SNRmax = range[3][1]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snf[i]*2.0/Tobs));
k = 0;
@@ -179,7 +175,7 @@ void draw_signal_amplitude(double *params, double *Snf, double Sa, gsl_rng *seed
}
//SNR defined with Sn(f) but Snf array holdes <n_i^2>
- params[3] = SNR/sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
+ params[3] = SNR/sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snf[i]*2.0/Tobs));
// FILE *temp=fopen("prior.dat","a");
// fprintf(temp,"%lg\n",params[3]);
@@ -187,22 +183,20 @@ void draw_signal_amplitude(double *params, double *Snf, double Sa, gsl_rng *seed
}
-void draw_uniform_amplitude(double *params, double *Snf, double Sa, gsl_rng *seed, double Tobs, double **range)
+void draw_uniform_amplitude(double *params, double *Snf, gsl_rng *seed, double Tobs, double **range)
{
int i;
double SNR;
- Sa = 1.0;
-
i = (int)(params[1]*Tobs);
- double SNRmin = range[3][0];//*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
- double SNRmax = range[3][1];//*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
+ double SNRmin = range[3][0];
+ double SNRmax = range[3][1];
SNR = SNRmin + (SNRmax-SNRmin)*uniform_draw(seed);
//SNR defined with Sn(f) but Snf array holdes <n_i^2>
- params[3] = SNR/sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
+ params[3] = SNR/sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snf[i]*2.0/Tobs));
}
@@ -264,7 +258,7 @@ void intrinsic_halfcycle_proposal(double *x, double *y, gsl_rng *seed)
/* ********************************************************************************** */
-void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *params, double **invSnf, double *Sn, struct Model *model, double **glitch, struct Data *data)
+void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *params, double **invSnf, struct Model *model, double **glitch, struct Data *data)//void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *params, double **invSnf, double *Sn, struct Model *model, double **glitch, struct Data *data)
{
int i, j, k;
int NI,N;
@@ -276,7 +270,7 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
double *epsilon;
double **s;
double stab;
-
+
struct Network *projectionP = malloc(sizeof(struct Network));
struct Network *projectionM = malloc(sizeof(struct Network));
@@ -355,9 +349,8 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
paramsP[i] += epsilon[i];
paramsM[i] -= epsilon[i];
-
- logLP = EvaluateExtrinsicMarkovianLogLikelihood(projectionP, paramsP, invSnf, Sn, model->signal, glitch, data, data->fmin, data->fmax);
- logLM = EvaluateExtrinsicMarkovianLogLikelihood(projectionM, paramsM, invSnf, Sn, model->signal, glitch, data, data->fmin, data->fmax);
+ logLP = EvaluateExtrinsicMarkovianLogLikelihood(projectionP, paramsP, invSnf, model->signal, glitch, data, data->fmin, data->fmax);
+ logLM = EvaluateExtrinsicMarkovianLogLikelihood(projectionM, paramsM, invSnf, model->signal, glitch, data, data->fmin, data->fmax);
epsilon[i] = 0.1/sqrt(-(logLM+logLP)/(epsilon[i]*epsilon[i]));
if(epsilon[i]!=epsilon[i])epsilon[i]=1.0e-5;
}
@@ -408,7 +401,7 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
}
}
- for(i=0; i<NE; i++) for(j=0; j<NE; j++) fisher->matrix[i][j] = network_nwip(data->imin,data->imax,hdiv[i], hdiv[j], invSnf, Sn, NI);
+ for(i=0; i<NE; i++) for(j=0; j<NE; j++) fisher->matrix[i][j] = network_nwip(data->imin,data->imax,hdiv[i], hdiv[j], invSnf, NI);
// stabilizers
for(i=0; i<NE; i++) fisher->matrix[i][i] += stab;
@@ -439,7 +432,7 @@ void extrinsic_fisher_update(struct Data *data, struct Model *model)
int i,j;
double *params = model->extParams;
- double *Sn = model->Sn;
+ //double *Sn = model->Sn;
double **glitch = malloc(data->NI*sizeof(double*));
for(i=0; i<data->NI; i++)
@@ -451,7 +444,7 @@ void extrinsic_fisher_update(struct Data *data, struct Model *model)
struct FisherMatrix *fisher = model->fisher;
//calculate elements of Fisher
- extrinsic_fisher_information_matrix(fisher, params, model->invSnf, Sn, model, glitch, data);
+ extrinsic_fisher_information_matrix(fisher, params, model->invSnf, model, glitch, data);
matrix_eigenstuff(fisher->matrix, fisher->evector, fisher->evalue, fisher->N);
@@ -492,7 +485,7 @@ void fisher_matrix_proposal(struct FisherMatrix *fisher, double *params, gsl_rng
for(i=0; i<N; i++) y[i] = params[i] + gaussian_draw(seed)*jump[i];
}
-void intrinsic_fisher_update(double *params, double *dparams, double *Snf, double Snx, double Tobs, int NW, int TFQFLAG)
+void intrinsic_fisher_update(double *params, double *dparams, double *Snf, double Tobs, int NW, int TFQFLAG)
{
int i;
double SNR;//,SNR2;
@@ -500,7 +493,7 @@ void intrinsic_fisher_update(double *params, double *dparams, double *Snf, doubl
double betasq;
//SNR defined with Sn(f) but Snf array holdes <n_i^2>
- SNR = params[3]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snx*Snf[i]*2.0/Tobs));
+ SNR = params[3]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snf[i]*2.0/Tobs));
if(SNR < 5.0) SNR = 5.0; // this caps the size of the proposed jumps
@@ -637,7 +630,7 @@ void intrinsic_fisher_matrix(double *params, struct FisherMatrix *fisher, int NW
}
-void intrinsic_fisher_proposal(UNUSED struct Model *model, double **range, double *paramsx, double *paramsy, double *Snf, double Snx, double *logpx, double *logpy, gsl_rng *seed, double Tobs, int *test, int NW, int TFQFLAG)
+void intrinsic_fisher_proposal(UNUSED struct Model *model, double **range, double *paramsx, double *paramsy, double *Snf, double *logpx, double *logpy, gsl_rng *seed, double Tobs, int *test, int NW, int TFQFLAG)
{
int k;
double x,y;
@@ -648,7 +641,7 @@ void intrinsic_fisher_proposal(UNUSED struct Model *model, double **range, doubl
double scale = 1./sqrt((double)NW);//Number of intrinsic parameters
// this is a crude diagonal Fisher matrix for individual sine-gaussians
- intrinsic_fisher_update(paramsx,dparamsx,Snf,Snx,Tobs, NW, TFQFLAG);
+ intrinsic_fisher_update(paramsx,dparamsx,Snf, Tobs, NW, TFQFLAG);
// [0] t0 [1] f0 [2] Q [3] Amp [4] phi ([5] beta chirplet)
@@ -678,7 +671,7 @@ void intrinsic_fisher_proposal(UNUSED struct Model *model, double **range, doubl
qyx -= log(y);
//need fisher at new location to compute reverse probability
- intrinsic_fisher_update(paramsy,dparamsy,Snf,Snx,Tobs, NW, TFQFLAG);
+ intrinsic_fisher_update(paramsy,dparamsy,Snf,Tobs, NW, TFQFLAG);
//compute probabiltiy of reverse proposal
y = 1.0;
diff --git a/src/BayesWaveProposal.h b/src/BayesWaveProposal.h
index a602697159ea0e8b8d374eee4ab54e11710001fb..85915d3c088d690066783b62bcb79f412da8f924 100644
--- a/src/BayesWaveProposal.h
+++ b/src/BayesWaveProposal.h
@@ -33,9 +33,9 @@
void extrinsic_uniform_proposal(gsl_rng *seed, double *y);
-void draw_glitch_amplitude(double *params, double *Snf, double Sa, gsl_rng *seed, double Tobs, double **range, double SNRpeak);
+void draw_glitch_amplitude(double *params, double *Snf, gsl_rng *seed, double Tobs, double **range, double SNRpeak);
-void draw_signal_amplitude(double *params, double *Snf, double Sa, gsl_rng *seed, double Tobs, double **range, double SNRpeak);
+void draw_signal_amplitude(double *params, double *Snf, gsl_rng *seed, double Tobs, double **range, double SNRpeak);
void intrinsic_halfcycle_proposal(double *x, double *y, gsl_rng *seed);
@@ -47,13 +47,13 @@ void intrinsic_halfcycle_proposal(double *x, double *y, gsl_rng *seed);
void fisher_matrix_proposal(struct FisherMatrix *fisher, double *params, gsl_rng *seed, double *y);
-void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *params, double **invSnf, double *Sn, struct Model *model, double **glitch, struct Data *data);
+void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *params, double **invSnf, struct Model *model, double **glitch, struct Data *data);//void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *params, double **invSnf, double *Sn, struct Model *model, double **glitch, struct Data *data);
void extrinsic_fisher_update(struct Data *data, struct Model *model);
-void intrinsic_fisher_update(double *params, double *dparams, double *Snf, double Snx, double Tobs, int NW, int TFQFLAG);
+void intrinsic_fisher_update(double *params, double *dparams, double *Snf, double Tobs, int NW, int TFQFLAG);
-void intrinsic_fisher_proposal(UNUSED struct Model *model, double **range, double *paramsx, double *paramsy, double *Snf, double Snx, double *logpx, double *logpy, gsl_rng *seed, double Tobs, int *test, int NW, int TFQFLAG);
+void intrinsic_fisher_proposal(UNUSED struct Model *model, double **range, double *paramsx, double *paramsy, double *Snf, double *logpx, double *logpy, gsl_rng *seed, double Tobs, int *test, int NW, int TFQFLAG);
/* ********************************************************************************** */
/* */