Commit 8796780c authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

BayesLine and Paralle Tempering, together at last


git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@565 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent dc4c66fb
This diff is collapsed.
......@@ -23,16 +23,6 @@
#define FSTEP 10.0 // the stencil separation in Hz for the spline model. Going below 2 Hz is dangerous - will fight with line model
// as we are getting down to the line width of the broadest lines
typedef struct
{
int gnuplot;
int verbose;
int zerologL;
char ifile[64];
char ofile[64];
}bayesline_opt;
typedef struct
{
......@@ -78,33 +68,21 @@ typedef struct
typedef struct
{
int n;
double *f;
double *pow;
double *Sn;
double *Sbase;
double *Sline;
}psdParams;
typedef struct
{
double SAmin;// 1.0e-48
double SAmax;// 1.0e-30
double LQmin;// 1.0e2 // allowing Q's below ~ 1e2 starts a fight with the spline model as lines get too fat
double LQmax;// 1.0e8
double LAmin;// 1.0e-44
double LAmax;// 1.0e-30
double *invsigma; //variances for each frequency bin
double SAmin;
double SAmax;
double LQmin;
double LQmax;
double LAmin;
double LAmax;
//double *invsigma; //variances for each frequency bin
double *sigma; //variances for each frequency bin
double *mean; //means for each frequency bin
}BayesLinePriors;
struct BayesLineParams
{
bayesline_opt opts;
dataParams *data;
splineParams *spline;
splineParams *spline_x;
......@@ -121,37 +99,23 @@ struct BayesLineParams
double *sfreq;
double *Sbase;
double *Sline;
double Snorm;
double RJrate;
int *lines_histogram;
int constantLogLFlag;
double TwoDeltaT;
gsl_rng *r;
FILE *pipe;
FILE *splineChainFile;
FILE *lineChainFile;
};
void BayesLineRJMCMC(struct BayesLineParams *bayesline, double *freqData, double *psd, double *invpsd, double *splinePSD, int N, int cycle);
void BayesLineSetup(struct BayesLineParams *bptr, double *freqData, double fmin, double fmax, double deltaT, double Tobs);
void BayesLineRJMCMC(struct BayesLineParams *bayesline, double *freqData, double *psd, double *invpsd, double *splinePSD, int N, int cycle, double beta);
void BayesLineSearch(struct BayesLineParams *bptr, double *freqData, double fmin, double fmax, double deltaT, double Tobs);
void BayesLineNonMarkovianFit (struct BayesLineParams *bayesline, int *nj);
void BayesLineLorentzSplineMCMC (struct BayesLineParams *bayesline, int steps, int focus, double *dan);
void BayesLineLorentzSplineMCMC (struct BayesLineParams *bayesline, double heat, int steps, int focus, double *dan);
void BayesLineMarkovianSplineOnly (struct BayesLineParams *bayesline, int nspline, int jj);
void BayesLineMarkovianFocusedAnalysis (struct BayesLineParams *bayesline);
//void KILL(char* Message);
int bayesline_defaults(bayesline_opt *opts);
//double loglike (double *respow, double *Snf, int ncut);
double loglike_fit_spline(double *respow, double *Snf, int ncut);
double loglike_pm (double *respow, double *Sn, double *Snx, int ilow, int ihigh);
......@@ -166,23 +130,24 @@ void full_spectrum_spline(double *Sline, double *Sbase, double *sfreq, dataParam
void spectrum_spline(double *Sn, double *Sbase, double *sfreq, dataParams *data, lorentzianParams *lines, splineParams *spline);
void SpecFitSpline (BayesLinePriors *priors, bayesline_opt opts, int steps, double *freq, double *power, splineParams *spline, double *Snf, int ncut, gsl_rng *r);
void LorentzSplineFit (BayesLinePriors *priors, bayesline_opt opts, int steps, dataParams *data, lorentzianParams *lines, splineParams *spline, double *sfreq, double *spow, gsl_rng *r);
void SpecFitSpline (BayesLinePriors *priors, int zeroLogL, int steps, double *freq, double *power, splineParams *spline, double *Snf, int ncut, gsl_rng *r);
void LorentzSplineFit (BayesLinePriors *priors, int zeroLogL, int steps, dataParams *data, lorentzianParams *lines, splineParams *spline, double *sfreq, double *spow, gsl_rng *r);
void CubicSplineGSL(int N, double *x, double *y, int Nint, double *xint, double *yint);
void create_dataParams(dataParams *data, double *f, int n);
void create_lorentzianParams(lorentzianParams *lines, int size);
void copy_lorentzianParams(lorentzianParams *origin, lorentzianParams *copy);
void destroy_lorentzianParams(lorentzianParams *lines);
void create_splineParams(splineParams *spline, int size);
void copy_splineParams(splineParams *origin, splineParams *copy);
void destroy_splineParams(splineParams *spline);
void create_psdParams(psdParams *psd, int size);
void destroy_psdParams(psdParams *psd);
void gnuplot_full_spectrum(int n, double *f, double *p, double *Sn, double *Sbase, double *Sline, FILE *pipe);
void copy_bayesline_params(struct BayesLineParams *origin, struct BayesLineParams *copy);
void print_line_model(FILE *fptr, struct BayesLineParams *bayesline);
void print_spline_model(FILE *fptr, struct BayesLineParams *bayesline);
......@@ -313,53 +313,79 @@ int main(int argc, char *argv[])
/* */
/******************************************************************************/
/*
Setup BayesLine structure
*/
struct BayesLineParams **bayesline = malloc(NI*sizeof(struct BayesLineParams *));
data->bayesline = bayesline;
struct BayesLineParams ***bayesline = malloc(chain->NC*sizeof(struct BayesLineParams **));
/*
BayesLine spectral estimation search-phase
*/
if(data->bayesLineFlag)
{
fprintf(stdout,"\n ============ BayesLine ==============\n");
for(ifo=0; ifo<NI; ifo++)
/*
Setup BayesLine structure
*/
for(ic=0; ic<chain->NC; ic++)
{
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));
bayesline[ic] = malloc(data->NI*sizeof(struct BayesLineParams *));
for(ifo=0; ifo<NI; ifo++)
{
bayesline[ic][ifo] = malloc(sizeof(struct BayesLineParams));
bayesline[ic][ifo]->priors = malloc(sizeof(BayesLinePriors));
//bayesline[ic][ifo]->priors->invsigma = malloc((int)(Tobs*(fmax-fmin))*sizeof(double));
bayesline[ic][ifo]->priors->sigma = malloc((int)(Tobs*(fmax-fmin))*sizeof(double));
bayesline[ic][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");
sprintf(filename,"chains/%ssearch_lorentzchain_ifo%i.dat",data->runName,ifo);
data->bayesline[ifo]->lineChainFile = fopen(filename,"w");
//Set BayesLine priors based on the data channel being used
set_bayesline_priors(data->channels[ifo], bayesline[ic][ifo]);
//Set BayesLine priors based on the data channel being used
set_bayesline_priors(data->channels[ifo], data->bayesline[ifo]);
// set default flags
bayesline[ic][ifo]->constantLogLFlag = data->constantLogLFlag;
//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];
//Use initial estimate of PSD to set priors
for(i=0; i<(int)(Tobs*(fmax-fmin)); i++)
{
//bayesline[ic][ifo]->priors->invsigma[i] = 1./(psd[ifo][i+(int)(Tobs*fmin)]*1.0);
bayesline[ic][ifo]->priors->sigma[i] = psd[ifo][i+(int)(Tobs*fmin)]*0.1;
bayesline[ic][ifo]->priors->mean[i] = psd[ifo][i+(int)(Tobs*fmin)];//*bayesline[ic][ifo]->priors->invsigma[i];
}
//Allocates arrays and sets constants for for BayesLine
BayesLineSetup(bayesline[ic][ifo], freqData[ifo], fmin, fmax, deltaT, Tobs);
}
}
/*
BayesLine spectral estimation search-phase
*/
for(ifo=0; ifo<NI; ifo++)
{
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], 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);
fprintf(stdout,"\n");
for(ic=0; ic<chain->NC; ic++)
{
BayesLineSearch(bayesline[ic][ifo], freqData[ifo], fmin, fmax, deltaT, Tobs);
}
}
for(ifo=0; ifo<NI; ifo++)
{
for(ic=0; ic<chain->NC; ic++)
{
if(ic==0)
{
fprintf(stdout,"BayesLine characterization phase for IFO %i\n", ifo);
BayesLineRJMCMC(bayesline[ic][ifo], freqData[ifo], model[ic]->Snf[ifo], model[ic]->invSnf[ifo], model[ic]->SnS[ifo], N, 100000, 1.0);
}
else
{
copy_bayesline_params(bayesline[0][ifo], bayesline[ic][ifo]);
for(i=0; i<data->N/2; i++)
{
model[ic]->Snf[ifo][i] = model[0]->Snf[ifo][i];
model[ic]->SnS[ifo][i] = model[0]->SnS[ifo][i];
model[ic]->invSnf[ifo][i] = model[0]->invSnf[ifo][i];
}
}
}
}
fprintf(stdout,"\n");
}
......@@ -484,7 +510,7 @@ int main(int argc, char *argv[])
double cleanLogZ, cleanVarZ;
printf("characterizing PSD+glitch model...\n\n");
RJMCMC(data, model, chain, prior, &cleanLogZ, &cleanVarZ);
RJMCMC(data, model, bayesline, chain, prior, &cleanLogZ, &cleanVarZ);
/******************************************************************************/
/* */
......@@ -566,8 +592,8 @@ int main(int argc, char *argv[])
if(data->bayesLineFlag) export_cleaned_data(data, model[chain->index[0]]);
//Shut off BayesLine in favor of PSD fitting
data->bayesLineFlag = 0;
data->psdFitFlag = 1;
//data->bayesLineFlag = 0;
//data->psdFitFlag = 1;
}
......@@ -653,6 +679,7 @@ 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;
......@@ -667,7 +694,7 @@ int main(int argc, char *argv[])
reset_priors(data, prior);
reset_model(data, chain, prior, model);
RJMCMC(data, model, chain, prior, &data->logZnoise, &data->varZnoise);
RJMCMC(data, model, bayesline, chain, prior, &data->logZnoise, &data->varZnoise);
//shut off noise model for checkpointing
data->noiseModelFlag = 0;
......@@ -700,7 +727,7 @@ int main(int argc, char *argv[])
reset_priors(data, prior);
reset_model(data, chain, prior, model);
RJMCMC(data, model, chain, prior, &data->logZglitch, &data->varZglitch);
RJMCMC(data, model, bayesline, chain, prior, &data->logZglitch, &data->varZglitch);
/*if(data->glitchOnlyFlag==1 && strcmp(data->channels[0],"AdLIGO"))
outputFrameFile(runState->commandLine, data, model[chain->index[0]]);*/
......@@ -736,7 +763,7 @@ int main(int argc, char *argv[])
reset_priors(data, prior);
reset_model(data, chain, prior, model);
RJMCMC(data, model, chain, prior, &data->logZsignal, &data->varZsignal);
RJMCMC(data, model, bayesline, chain, prior, &data->logZsignal, &data->varZsignal);
//shut off signal model for checkpointing
data->signalModelFlag = 0;
......@@ -771,7 +798,7 @@ int main(int argc, char *argv[])
reset_priors(data, prior);
reset_model(data, chain, prior, model);
RJMCMC(data, model, chain, prior, &data->logZfull, &data->varZfull);
RJMCMC(data, model, bayesline, chain, prior, &data->logZfull, &data->varZfull);
fprintf(evidence,"full %lg %lg\n",data->logZfull,data->varZfull);
fflush(evidence);
......
......@@ -132,52 +132,54 @@ struct Network
struct Chain
{
int mc;
int NC; // Number of parallel chains
int burn; // Number of burn in samples
int mod0; // Counter keeping track of noise model
int mod1; // Counter keeping track of glitch model
int mod2; // Counter keeping track of signal model
int mod3; // Counter keeping track of glitch+signal model
int cycle; // # of model updates per BurstRJMCMC iteration
int count; // # of BurstRJMCMC iterations
int mcount; // # of MCMC attempts
int scount; // # of RJMCMC attempts
int xcount; // # of extrinsic attempts
int macc; // # of MCMC successess
int sacc; // # of RJMCMC successes
int xacc; // # of extrinsic successes
int burnFlag;// flag telling proposals if we are in burn-in phase
int *index; // keep track of which chain is at which temperature
int *ptacc;
int *ptprop;
double *A;
long seed;
double modelRate; //fraction of updates to signal model
double rjmcmcRate; //fraction of trans- vs. inter-dimensional moves
double intPropRate; //fraction of prior vs. fisher intrinsic moves
double *dT;
double beta;
double tempMin;
double tempStep;
double *temperature;
double *avgLogLikelihood;
double *varLogLikelihood;
int NC; // Number of parallel chains
int burn; // Number of burn in samples
int mod0; // Counter keeping track of noise model
int mod1; // Counter keeping track of glitch model
int mod2; // Counter keeping track of signal model
int mod3; // Counter keeping track of glitch+signal model
int cycle; // # of model updates per BurstRJMCMC iteration
int count; // # of BurstRJMCMC iterations
int mcount; // # of MCMC attempts
int scount; // # of RJMCMC attempts
int xcount; // # of extrinsic attempts
int macc; // # of MCMC successess
int sacc; // # of RJMCMC successes
int xacc; // # of extrinsic successes
int burnFlag;// flag telling proposals if we are in burn-in phase
int *index; // keep track of which chain is at which temperature
int *ptacc;
int *ptprop;
double *A;
long seed;
double modelRate; //fraction of updates to signal model
double rjmcmcRate; //fraction of trans- vs. inter-dimensional moves
double intPropRate; //fraction of prior vs. fisher intrinsic moves
double *dT;
double beta;
double tempMin;
double tempStep;
double *temperature;
double *avgLogLikelihood;
double *varLogLikelihood;
int nPoints;// = 3*chain->count/4/chain->cycle;
double **logLchain;//=dmatrix(0,chain->NC-1,0,nPoints-1);
FILE *runFile;
FILE **intChainFile;
FILE **intWaveChainFile;
FILE ***intGlitchChainFile;
FILE *runFile;
FILE **intChainFile;
FILE **intWaveChainFile;
FILE ***intGlitchChainFile;
FILE ***lineChainFile;
FILE ***splineChainFile;
};
......@@ -313,8 +315,6 @@ struct Data
char outputDirectory[1000];
char **channels;
struct BayesLineParams **bayesline;
ProcessParamsTable *commandLine;
LALDetector **detector;
LIGOTimeGPS epoch;
......
#include "BayesLine.h"
#include "BayesWave.h"
#include "BayesWaveIO.h"
......@@ -862,12 +861,13 @@ void print_run_flags(FILE *fptr, struct Data *data, struct Prior *prior)
fprintf(fptr, "\n");
}
void print_chain_files(struct Data *data, struct Chain *chain, struct Model **model, int ic, double logZ, double varZ)
void print_chain_files(struct Data *data, struct Chain *chain, struct Model **model, struct BayesLineParams ***bayesline, int ic, double logZ, double varZ)
{
int ifo;
int NI = data->NI;
struct Model *model_x = model[chain->index[ic]];
struct BayesLineParams **psd_x = bayesline[chain->index[ic]];
//TODO "model" chain is very redundent
//print "model" chain (logL, PSD parameters, etc.)
......@@ -879,6 +879,16 @@ void print_chain_files(struct Data *data, struct Chain *chain, struct Model **mo
//print glitch model parameters
if(data->glitchFlag) for(ifo=0; ifo<NI; ifo++) print_glitch_model(chain->intGlitchChainFile[ifo][ic], model_x->glitch[ifo]);
//print PSD parameters
if(data->bayesLineFlag)
{
for(ifo=0; ifo<NI; ifo++)
{
print_line_model(chain->lineChainFile[ic][ifo],psd_x[ifo]);
print_spline_model(chain->splineChainFile[ic][ifo],psd_x[ifo]);
}
}
if(data->verboseFlag)fflush(chain->intChainFile[ic]);
}
......
......@@ -27,7 +27,7 @@ void print_bayesline_spectrum(char filename[], double *f, double *power, double
void print_version(FILE *fptr);
void print_run_stats(FILE *fptr, struct Data *data, struct Chain *chain);
void print_run_flags(FILE *fptr, struct Data *data, struct Prior *prior);
void print_chain_files(struct Data *data, struct Chain *chain, struct Model **model, int ic, double logZ, double varZ);
void print_chain_files(struct Data *data, struct Chain *chain, struct Model **model, struct BayesLineParams ***bayesline, int ic, double logZ, double varZ);
void print_chain_status(struct Data *data, struct Chain *chain, struct Model **model, int searchFlag, int NC);
void print_model(FILE *fptr, struct Data *data, struct Chain *chain, struct Model *model);
......
......@@ -3,6 +3,7 @@
#include <math.h>
#include "BayesWave.h"
#include "BayesLine.h"
#include "BayesWaveIO.h"
#include "BayesWaveMath.h"
#include "BayesWavePrior.h"
......
......@@ -65,7 +65,7 @@ static void check_interrupt(volatile sig_atomic_t __bayeswave_exitFlag)
}
}
static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior *prior, struct Model **model, char modelname[])
static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior *prior, struct Model **model, struct BayesLineParams ***bayesline, char modelname[])
{
int i,j;
int ic,ifo;
......@@ -253,15 +253,15 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
}
//Parameters for reproducing PSD model
if(data->bayesLineFlag && ic==0)
if(data->bayesLineFlag)
{
for(ifo=0; ifo<NI; ifo++)
{
sprintf(filename,"chains/%s_spline_ifo%i.dat.0",modelname,ifo);
data->bayesline[ifo]->splineChainFile = fopen(filename,"w");
sprintf(filename,"chains/%s_spline_ifo%i.dat.%i",modelname,ifo,ic);
chain->splineChainFile[ic][ifo] = fopen(filename,"w");
sprintf(filename,"chains/%s_lorentz_ifo%i.dat.0",modelname,ifo);
data->bayesline[ifo]->lineChainFile = fopen(filename,"w");
sprintf(filename,"chains/%s_lorentz_ifo%i.dat.%i",modelname,ifo,ic);
chain->lineChainFile[ic][ifo] = fopen(filename,"w");
}
}
......@@ -269,7 +269,7 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
}
}
static void resume_sampler(struct Data *data, struct Chain *chain, struct Prior *prior, struct Model **model)
static void resume_sampler(struct Data *data, struct Chain *chain, struct Prior *prior, struct Model **model, struct BayesLineParams ***bayesline)
{
//TODO: set internal mcmc flags right based on chain->mc
int i;
......@@ -510,15 +510,15 @@ static void resume_sampler(struct Data *data, struct Chain *chain, struct Prior
}
//Parameters for reproducing PSD model
if(data->bayesLineFlag && ic==0)
if(data->bayesLineFlag)
{
for(ifo=0; ifo<NI; ifo++)
{
sprintf(filename,"chains/%s_spline_ifo%i.dat.0",modelname,ifo);
data->bayesline[ifo]->splineChainFile = fopen(filename,"a");
sprintf(filename,"chains/%s_spline_ifo%i.dat.%i",modelname,ifo,ic);
chain->splineChainFile[ic][ifo] = fopen(filename,"a");
sprintf(filename,"chains/%s_lorentz_ifo%i.dat.0",modelname,ifo);
data->bayesline[ifo]->lineChainFile = fopen(filename,"a");
sprintf(filename,"chains/%s_lorentz_ifo%i.dat.%i",modelname,ifo,ic);
chain->lineChainFile[ic][ifo] = fopen(filename,"a");
}
}
......@@ -720,7 +720,7 @@ static void check_save(volatile sig_atomic_t __bayeswave_saveStateFlag, struct D
/* */
/* ********************************************************************************** */
void RJMCMC(struct Data *data, struct Model **model, struct Chain *chain, struct Prior *prior, double *logEvidence, double *varLogEvidence)
void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***bayesline, struct Chain *chain, struct Prior *prior, double *logEvidence, double *varLogEvidence)
{
print_run_stats(stdout, data, chain);
......@@ -865,7 +865,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct Chain *chain, struct
if(data->checkpointFlag && data->resumeFlag && !data->cleanModelFlag)
{
printf("try resuming...\n");
resume_sampler(data, chain, prior, model);
resume_sampler(data, chain, prior, model, bayesline);
data->resumeFlag = 0;
/******************************************************************************/
......@@ -919,7 +919,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct Chain *chain, struct
burnFlag=1;
}
}
else restart_sampler(data, chain, prior, model, modelname);
else restart_sampler(data, chain, prior, model, bayesline, modelname);
//Set MAP model to initial state
......@@ -1070,11 +1070,11 @@ void RJMCMC(struct Data *data, struct Model **model, struct Chain *chain, struct
}
//print state of cold chain to file
print_chain_files(data, chain, model, 0, logZ, varZ);
print_chain_files(data, chain, model, bayesline, 0, logZ, varZ);
if(data->verboseFlag)
{
//only print hot chains if asked to
for(ic=1; ic<NC; ic++)print_chain_files(data, chain, model, ic, logZ, varZ);
for(ic=1; ic<NC; ic++)print_chain_files(data, chain, model, bayesline, ic, logZ, varZ);
}
......@@ -1130,11 +1130,10 @@ void RJMCMC(struct Data *data, struct Model **model, struct Chain *chain, struct
//This function executes [cycle] extrinsic parameter updates, common to all geocenter wavelets.
if(model[chain->index[ic]]->signal->size>0 && !data->fixExtrinsicFlag) EvolveExtrinsicParameters(data, prior, model, chain, &chain->seed, ic);
//update PSD with BayesLine
if(data->bayesLineFlag) EvolveBayesLineParameters(data, model, bayesline, chain, prior, ic);
}
//update PSD with BayesLine
if(data->bayesLineFlag) EvolveBayesLineParameters(data, model, chain, prior);
//After so many iterations recompute the residuals and likelihood (prevent accumulation of roundoff error)
recompute_residual(data, model, chain);
......@@ -1210,14 +1209,16 @@ void RJMCMC(struct Data *data, struct Model **model, struct Chain *chain, struct
if(data->bayesLineFlag)
{
sprintf(filename,"waveforms/%s_psd_%d.dat.%i", modelname, frame,ifo);
print_bayesline_spectrum(filename, data->bayesline[ifo]->freq, data->bayesline[ifo]->spow, data->bayesline[ifo]->Sbase, data->bayesline[ifo]->Sline, data->bayesline[ifo]->data->ncut);
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);
sprintf(filename,"waveforms/%s_psd_%d.dat.%i.hot", modelname, frame,ifo);
print_bayesline_spectrum(filename, bayesline[chain->index[chain->NC-1]][ifo]->freq, bayesline[chain->index[chain->NC-1]][ifo]->spow, bayesline[chain->index[chain->NC-1]][ifo]->Sbase, bayesline[chain->index[chain->NC-1]][ifo]->Sline, bayesline[chain->index[chain->NC-1]][ifo]->data->ncut);
}
}
write_gnuplot_script_frame(script, modelname, data->runName, data->Tobs, data->runPhase, data->signalFlag, data->glitchFlag, frame, data->NI);
if(data->bayesLineFlag)
{
write_bayesline_gnuplot_script_frame(psdscript, modelname, data->Tobs, frame);