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

+Fixed extrinsic proposal

+PSD stored in Model structure
+First pass at BayesLine priors



git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@562 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent 3778f4b7
/*********************************************************************************/
/* */
/* */
/* BayesLine v1.0. Written by Neil J. Cornish between 5/20/13 - 6/5/13 */
/* */
/* NJC thanks Tyson Littenberg and Will Farr for comments and suggestions */
/* */
/* */
/* */
/* BayesLine fits the LIGO/Virgo power spectra using a model made up */
/* of N Lorentzian lines (described by central frequency f, quality */
/* factor Q and amplitude A) and cubic spline with M control points. */
......@@ -20,7 +13,6 @@
/* transdimensional (and very free in its movement between dimensions) */
/* it will not "over-fit" the spectral model. */
/* */
/* */
/*********************************************************************************/
......@@ -80,7 +72,21 @@ double loglike_fit_spline(double *respow, double *Snf, int ncut)
return(lgl);
}
static double logprior(double *invsigma, double *mean, double *Snf, int ilow, int ihigh)
{
double x,lgp;
int i;
//leaving off normalizations since they cancel in Hastings ratio
lgp = 0;
for(i=ilow; i<ihigh; i++)
{
x = mean[i] - Snf[i]*invsigma[i];
lgp -= x*x;
}
return (0.5*lgp);
}
static double loglike(double *respow, double *Snf, int ncut)
{
......@@ -1367,6 +1373,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, int steps, in
double xsm, pmax, fcl, fch, dff;
double baseav;
double logpx=0.0, logpy=0.0, x, y, z, beta;
double logPsy,logPsx;
double Ac;
double *fprop;
double *sdatay;
......@@ -1472,6 +1479,8 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, int steps, in
if(!opts.zerologL) logLx = loglike(power, Sn, ncut);
else logLx = 1.0;
logPsx = logprior(priors->invsigma, priors->mean, Snx, 0, ncut);
// set up proposal for frequency jumps
xsm =0.0;
pmax = -1.0;
......@@ -1968,7 +1977,9 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, int steps, in
logpy += zeta*(double)(lines_y->n);
logpx += zeta*(double)(lines_x->n);
logH = (logLy - logLx) - logpy + logpx;
logPsy = logprior(priors->invsigma, priors->mean, Sn, 0, ncut);
logH = (logLy - logLx) - logpy + logpx + logPsy - logPsx;
alpha = log(gsl_rng_uniform(r));
if(logH > alpha)
......@@ -1977,6 +1988,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, int steps, in
if(typ == 1) ac1++;
if(typ == 4) ac2++;
logLx = logLy;
logPsx = logPsy;
lines_x->n = lines_y->n;
nsx = nsy;
for(i=0; i< ncut; i++)
......@@ -2453,13 +2465,13 @@ void BayesLineRJMCMC(struct BayesLineParams *bayesline, double *freqData, double
{
if(i>=imin && i<imax)
{
psd[i] = bayesline->Snf[i-imin];///bayesline->data->Tobs/(bayesline->TwoDeltaT/(double)(N));
splinePSD[i] = bayesline->Sbase[i-imin];///bayesline->data->Tobs/(bayesline->TwoDeltaT/(double)(N));
psd[i] = bayesline->Snf[i-imin];
splinePSD[i] = bayesline->Sbase[i-imin];
}
else
{
psd[i] = bayesline->priors->SAmax;///bayesline->data->Tobs/(bayesline->TwoDeltaT/(double)(N));
splinePSD[i] = bayesline->priors->SAmax;///bayesline->data->Tobs/(bayesline->TwoDeltaT/(double)(N));
psd[i] = bayesline->priors->SAmax;
splinePSD[i] = bayesline->priors->SAmax;
}
invpsd[i] = 1./psd[i];
}
......
......@@ -97,6 +97,9 @@ typedef struct
double LAmin;// 1.0e-44
double LAmax;// 1.0e-30
double *invsigma; //variances for each frequency bin
double *mean; //means for each frequency bin
}BayesLinePriors;
struct BayesLineParams
......
......@@ -196,7 +196,7 @@ int main(int argc, char *argv[])
Setup DATA structure
*/
initialize_data(data,freqData,N,tsize,psd,Tobs,NI,fmin,fmax);
initialize_data(data,freqData,N,tsize,Tobs,NI,fmin,fmax);
//print run summary file
sprintf(filename,"%sbayeswave.run",data->runName);
......@@ -285,7 +285,7 @@ int main(int argc, char *argv[])
for(ic=0; ic<chain->NC; ic++)
{
model[ic] = malloc(sizeof(struct Model));
initialize_model(model[ic], data, prior, &chain->seed);
initialize_model(model[ic], data, prior, psd, &chain->seed);
}
......@@ -330,6 +330,9 @@ int main(int argc, char *argv[])
{
bayesline[ifo] = malloc(sizeof(struct BayesLineParams));
bayesline[ifo]->priors = malloc(sizeof(BayesLinePriors));
bayesline[ifo]->priors->invsigma = malloc((int)(Tobs*(fmax-fmin))*sizeof(double));
bayesline[ifo]->priors->mean = malloc((int)(Tobs*(fmax-fmin))*sizeof(double));
sprintf(filename,"chains/%ssearch_splinechain_ifo%i.dat",data->runName,ifo);
data->bayesline[ifo]->splineChainFile = fopen(filename,"w");
......@@ -339,11 +342,18 @@ int main(int argc, char *argv[])
//Set BayesLine priors based on the data channel being used
set_bayesline_priors(data->channels[ifo], data->bayesline[ifo]);
//Use initial estimate of PSD to set priors
for(i=0; i<(int)(Tobs*(fmax-fmin)); i++)
{
bayesline[ifo]->priors->invsigma[i] = 1./(psd[ifo][i+(int)(Tobs*fmin)]*0.2);
bayesline[ifo]->priors->mean[i] = psd[ifo][i+(int)(Tobs*fmin)]*bayesline[ifo]->priors->invsigma[i];
}
fprintf(stdout,"BayesLine search phase for IFO %i\n", ifo);
BayesLineSearch(bayesline[ifo], freqData[ifo], fmin, fmax, deltaT, Tobs);
fprintf(stdout,"BayesLine characterization phase for IFO %i\n", ifo);
BayesLineRJMCMC(bayesline[ifo], freqData[ifo], data->Snf[ifo], data->invSnf[ifo], data->SnS[ifo], N, 1000000);
BayesLineRJMCMC(bayesline[ifo], freqData[ifo], model[chain->index[0]]->Snf[ifo], model[chain->index[0]]->invSnf[ifo], model[chain->index[0]]->SnS[ifo], N, 1000000);
fclose(data->bayesline[ifo]->splineChainFile);
fclose(data->bayesline[ifo]->lineChainFile);
......@@ -520,7 +530,7 @@ int main(int argc, char *argv[])
//regress all wavelest from residual array for plotting
m->wavelet(r[ifo], g->intParams[i], N, -1, data->Tobs);
}
rSNR += fourier_nwip(data->imin,data->imax,g->templates,g->templates,data->invSnf[ifo]);
rSNR += fourier_nwip(data->imin,data->imax,g->templates,g->templates,m->invSnf[ifo]);
}
fprintf(stdout, " Residual SNR = %g\n",sqrt(rSNR));
......@@ -531,7 +541,7 @@ int main(int argc, char *argv[])
/* */
/******************************************************************************/
if(sqrt(rSNR)>100.0) resize_chain(data, chain, prior, model, chain->NC+5);
if(sqrt(rSNR)>100.0) resize_chain(data, chain, prior, model, psd, chain->NC+5);
/******************************************************************************/
/* */
......@@ -553,7 +563,7 @@ int main(int argc, char *argv[])
free(r);
//Export BayesLine PSD and cleaned data so we can re-start runs
if(data->bayesLineFlag) export_cleaned_data(data);
if(data->bayesLineFlag) export_cleaned_data(data, model[chain->index[0]]);
//Shut off BayesLine in favor of PSD fitting
data->bayesLineFlag = 0;
......@@ -577,7 +587,7 @@ int main(int argc, char *argv[])
dataPtr = runState->data;
while(dataPtr!=NULL)
{
for(i=0; i<N/2; i++) dataPtr->oneSidedNoisePowerSpectrum->data->data[i] = 2.0/Tobs/data->invSnf[ifo][i];
for(i=0; i<N/2; i++) dataPtr->oneSidedNoisePowerSpectrum->data->data[i] = 2.0/Tobs/model[chain->index[0]]->invSnf[ifo][i];
dataPtr = dataPtr->next;
ifo++;
}
......@@ -597,7 +607,7 @@ int main(int argc, char *argv[])
while(dataPtr!=NULL)
{
// Fill runState PSD with results from BayesLine
for(i=0; i<N/2; i++) dataPtr->oneSidedNoisePowerSpectrum->data->data[i] = 2.0/Tobs/data->invSnf[ifo][i];
for(i=0; i<N/2; i++) dataPtr->oneSidedNoisePowerSpectrum->data->data[i] = 2.0/Tobs/model[chain->index[0]]->invSnf[ifo][i];
dataPtr = dataPtr->next;
ifo++;
}
......
......@@ -96,6 +96,11 @@ struct Model
int size;
double *Sn;
double **Snf;
double **invSnf;
double **SnS;
double *SnGeo;
double logL;
double *detLogL;
double *extParams;
......@@ -219,74 +224,70 @@ struct Prior
struct Data
{
int N;
int NI;
int tsize;
int fsize;
int imin;
int imax;
int Dmax;
int Dmin;
int runPhase;
int psdFitFlag;
int noiseSimFlag;
int bayesLineFlag;
int stochasticFlag;
int fullModelFlag;
int noiseModelFlag;
int glitchModelFlag;
int signalModelFlag;
int cleanModelFlag;
int cleanOnlyFlag;
int glitchFlag;
int signalFlag;
int noiseFlag;
int fixIntrinsicFlag;
int fixExtrinsicFlag;
int clusterPriorFlag;
int backgroundPriorFlag;
int amplitudePriorFlag;
int amplitudeProposalFlag;
int signalAmplitudePriorFlag;
int extrinsicAmplitudeFlag;
int orientationProposalFlag;
int clusterProposalFlag;
int adaptTemperatureFlag;
int splineEvidenceFlag;
int constantLogLFlag;
int gnuplotFlag;
int verboseFlag;
int N;
int NI;
int tsize;
int fsize;
int imin;
int imax;
int Dmax;
int Dmin;
int runPhase;
int psdFitFlag;
int noiseSimFlag;
int bayesLineFlag;
int stochasticFlag;
int fullModelFlag;
int noiseModelFlag;
int glitchModelFlag;
int signalModelFlag;
int cleanModelFlag;
int cleanOnlyFlag;
int glitchFlag;
int signalFlag;
int noiseFlag;
int fixIntrinsicFlag;
int fixExtrinsicFlag;
int clusterPriorFlag;
int backgroundPriorFlag;
int amplitudePriorFlag;
int amplitudeProposalFlag;
int signalAmplitudePriorFlag;
int extrinsicAmplitudeFlag;
int orientationProposalFlag;
int clusterProposalFlag;
int adaptTemperatureFlag;
int splineEvidenceFlag;
int constantLogLFlag;
int gnuplotFlag;
int verboseFlag;
int checkpointFlag;
int resumeFlag;
int rjFlag;
double dt;
double df;
double fmin;
double fmax;
double Tobs;
double Twin;
double Qmin;
double Qmax;
double gmst;
double logLc;
double TwoDeltaToverN;
double TFtoProx;
double **r;
double **s;
double **rh;
double *h;
double **Snf;
double **invSnf;
double **SnS;
double *SnGeo;
double *ORF;
int rjFlag;
double dt;
double df;
double fmin;
double fmax;
double Tobs;
double Twin;
double Qmin;
double Qmax;
double gmst;
double logLc;
double TwoDeltaToverN;
double TFtoProx;
double **r;
double **s;
double **rh;
double *h;
double *ORF;
double logZglitch;
double logZnoise;
......@@ -300,23 +301,23 @@ struct Data
void (*signal_amplitude_proposal)(double *, double *, double, long *, double, double **, double);
void (*glitch_amplitude_proposal)(double *, double *, double, long *, double, double **, double);
double (*signal_amplitude_prior) (double *, double *, double, double, double);
double (*glitch_amplitude_prior) (double *, 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 *, double **, struct Data*, double, double);
char runName[100];
char outputDirectory[1000];
char **channels;
struct BayesLineParams **bayesline;
ProcessParamsTable *commandLine;
LALDetector **detector;
LIGOTimeGPS epoch;
void (*glitch_amplitude_proposal)(double *, double *, double, long *, double, double **, double);
double (*signal_amplitude_prior) (double *, double *, double, double, double);
double (*glitch_amplitude_prior) (double *, 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 *, double *, double **, struct Data*, double, double);
char runName[100];
char outputDirectory[1000];
char **channels;
struct BayesLineParams **bayesline;
ProcessParamsTable *commandLine;
LALDetector **detector;
LIGOTimeGPS epoch;
};
......
......@@ -2354,7 +2354,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
for(i=1; i<=model->signal->size; i++)
{
ii = model->signal->index[i];
intrinsic_fisher_update(model->signal->intParams[ii], sigmas, data->SnGeo, 1.0, data->Tobs);
intrinsic_fisher_update(model->signal->intParams[ii], sigmas, model->SnGeo, 1.0, data->Tobs);
for(j=0; j<NW; j++) detC += log(sigmas[j]);
}
}
......@@ -2367,7 +2367,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, data->Snf[ifo], 1.0, data->Tobs);
intrinsic_fisher_update(model->glitch[ifo]->intParams[ii], sigmas, model->Snf[ifo], 1.0, data->Tobs);
//for(j=0; j<NW; j++)printf(" %i: %g\n",j,sigmas[j]);
for(j=0; j<NW; j++) detC += log(sigmas[j]);
}
......@@ -2409,7 +2409,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
for(i=0; i<100; i++)
{
draw_wavelet_params(-1, params, data, prior->range, &chain->seed);
data->signal_amplitude_proposal(params, data->SnGeo, 1.0, &chain->seed, data->Tobs, prior->range,prior->sSNRpeak);
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];
}
......@@ -2436,7 +2436,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
for(i=0; i<100; i++)
{
draw_wavelet_params(ifo, params, data, prior->range, &chain->seed);
data->glitch_amplitude_proposal(params, data->Snf[ifo], 1.0, &chain->seed, data->Tobs, prior->range, prior->gSNRpeak);
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];
}
......
......@@ -212,11 +212,14 @@ void BayesWaveBurstInject(ProcessParamsTable *commandLine, struct Data *data, st
//User-friendly arrays for PSD, signal, and injected waveform
double **hinj = malloc(NI*sizeof(double*));
double **ginj = malloc(NI*sizeof(double*));
double **psd = malloc(NI*sizeof(double *));
for(ifo=0; ifo<NI; ifo++)
{
hinj[ifo] = malloc(N*sizeof(double));
ginj[ifo] = malloc(N*sizeof(double));
psd[ifo] = malloc(N/2*sizeof(double));
for(i=0; i<data->N/2; i++) psd[ifo][i]=1.0;
}
......@@ -253,7 +256,7 @@ void BayesWaveBurstInject(ProcessParamsTable *commandLine, struct Data *data, st
if(strcmp(injmodel, "glitch") == 0)
{
data->glitchFlag = 1;
initialize_model(model, data, prior, &chain->seed);
initialize_model(model, data, prior, psd, &chain->seed);
for(ifo=0; ifo<NI; ifo++)model->Sn[ifo]=1.0;
FILE **glitchParams = malloc(NI*sizeof(FILE *));
......@@ -291,7 +294,7 @@ void BayesWaveBurstInject(ProcessParamsTable *commandLine, struct Data *data, st
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, data->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], model->Sn[ifo]);
fprintf(snrFile,"%g ",GNR[ifo]);
}
fprintf(snrFile,"\n");
......@@ -314,7 +317,7 @@ void BayesWaveBurstInject(ProcessParamsTable *commandLine, struct Data *data, st
if(strcmp(injmodel, "signal") == 0)
{
data->signalFlag = 1;
initialize_model(model, data, prior, &chain->seed);
initialize_model(model, data, prior, psd, &chain->seed);
for(ifo=0; ifo<NI; ifo++)model->Sn[ifo]=1.0;
//File containing signal-model waveform parameters
......@@ -339,8 +342,8 @@ void BayesWaveBurstInject(ProcessParamsTable *commandLine, struct Data *data, st
double SNRnet;
double SNRifo[NI];
SNRnet = network_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model->response, data->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], data->invSnf[ifo], model->Sn[ifo]);
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]);
sprintf(filename,"%ssnr/snr.dat",data->runName);
snrFile=fopen(filename,"w");
......@@ -370,9 +373,11 @@ void BayesWaveBurstInject(ProcessParamsTable *commandLine, struct Data *data, st
{
free(hinj[ifo]);
free(ginj[ifo]);
free(psd[ifo]);
}
free(hinj);
free(ginj);
free(psd);
data->signalFlag = signalFlag;
data->glitchFlag = glitchFlag;
......@@ -934,12 +939,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, data->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, model[chain->index[0]]->Sn, 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, data->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], model[chain->index[0]]->Sn[ifo]);
fprintf(stdout, " g%iSNR=%.0f ", ifo, SNR);
}
fprintf(stdout, "\n");
......@@ -1528,11 +1533,11 @@ 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, data->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], 1, data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
}
}
void export_cleaned_data(struct Data *data)
void export_cleaned_data(struct Data *data, struct Model *model)
{
int NI = data->NI;
double Tobs = data->Tobs;
......@@ -1557,17 +1562,17 @@ void export_cleaned_data(struct Data *data)
imax = data->imax;//(int)(floor(data->bayesline[0]->data->fhigh * data->bayesline[0]->data->Tobs));
//output PSD (padded for ease of use in LALInferenceReadData())
fprintf(psdFile, "%.12g %.12g\n",(double)(imin-1)/data->Tobs, data->Snf[ifo][imin]/(Tobs/2.0));
fprintf(asdFile, "%.12g %.12g\n",(double)(imin-1)/data->Tobs, sqrt(data->Snf[ifo][imin]/(Tobs/2.0)));
fprintf(psdFile, "%.12g %.12g\n",(double)(imin-1)/data->Tobs, model->Snf[ifo][imin]/(Tobs/2.0));
fprintf(asdFile, "%.12g %.12g\n",(double)(imin-1)/data->Tobs, sqrt(model->Snf[ifo][imin]/(Tobs/2.0)));
fprintf(resFile, "%.12g %.12g %.12g\n",(double)(imin-1)/data->Tobs,data->s[ifo][2*imin],data->s[ifo][2*imin+1]);
for(i=imin; i<imax; i++)
{
fprintf(psdFile, "%.12g %.12g\n",(double)i/data->Tobs, data->Snf[ifo][i]/(Tobs/2.0));
fprintf(asdFile, "%.12g %.12g\n",(double)i/data->Tobs, sqrt(data->Snf[ifo][i]/(Tobs/2.0)));
fprintf(psdFile, "%.12g %.12g\n",(double)i/data->Tobs, model->Snf[ifo][i]/(Tobs/2.0));
fprintf(asdFile, "%.12g %.12g\n",(double)i/data->Tobs, sqrt(model->Snf[ifo][i]/(Tobs/2.0)));
fprintf(resFile, "%.12g %.12g %.12g\n",(double)i/data->Tobs,data->s[ifo][2*i],data->s[ifo][2*i+1]);
}
fprintf(psdFile, "%.12g %.12g\n",(double)imax/data->Tobs, data->Snf[ifo][imax-1]/(Tobs/2.0));
fprintf(asdFile, "%.12g %.12g\n",(double)imax/data->Tobs, sqrt(data->Snf[ifo][imax-1]/(Tobs/2.0)));
fprintf(psdFile, "%.12g %.12g\n",(double)imax/data->Tobs, model->Snf[ifo][imax-1]/(Tobs/2.0));
fprintf(asdFile, "%.12g %.12g\n",(double)imax/data->Tobs, sqrt(model->Snf[ifo][imax-1]/(Tobs/2.0)));
fprintf(resFile, "%.12g %.12g %.12g\n",(double)imax/data->Tobs,data->s[ifo][2*(imax-1)],data->s[ifo][2*(imax-1)+1]);
fclose(psdFile);
......
......@@ -44,7 +44,7 @@ void parse_glitch_parameters(struct Data *data, struct Model *model, FILE **para
void parse_signal_parameters(struct Data *data, struct Model *model, FILE *paramsfile, double **hrec, double **hrecPlus, double **hrecCross);
void dump_time_domain_waveform(struct Data *data, struct Model *model, struct Prior *prior, char root[]);
void export_cleaned_data(struct Data *data);
void export_cleaned_data(struct Data *data, struct Model *model);
void system_pause();
......@@ -3,6 +3,7 @@
#include <math.h>
#include "BayesWave.h"
#include "BayesWaveIO.h"
#include "BayesWaveMath.h"
#include "BayesWavePrior.h"
#include "BayesWaveModel.h"
......@@ -88,15 +89,13 @@ double loglike_stochastic(int NI, int imin, int imax, double **r, double ***invC
return logL;
}
double loglike_maxNET(struct Data *data, struct Prior *prior, double **sNET, double **hNET, double *SnNET, double t0, double *dt, double *dphi, double *dA)
double loglike_maxNET(struct Data *data, struct Prior *prior, double **sNET, double **hNET, double **SnfNET, double **invSnfNET, double *SnNET, double t0, double *dt, double *dphi, double *dA)
{
int N = data->N;
int halfN = N/2;
int NI = data->NI;
int imin = data->imin;
int imax = data->imax;
double **invSnfNET = data->invSnf;
double **SnfNET = data->Snf;
double Tobs = data->Tobs;
int tmin = (int)floor(prior->range[0][0]/(Tobs/(double)N));
......@@ -311,14 +310,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->Sn, model_x->background);
ComputeNoiseCorrelationMatrix(data, model_x->Snf, model_x->Sn, 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++)