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

BayesLine priors and checkpointing


git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@568 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent c71d2670
......@@ -18,6 +18,13 @@
#include "BayesLine.h"
static void system_pause()
{
printf("Press Any Key to Continue\n");
getchar();
}
double sample(double *fprop, double pmax, dataParams *data, gsl_rng *r)
{
int i;
......@@ -72,7 +79,8 @@ double loglike_fit_spline(double *respow, double *Snf, int ncut)
}
//static double logprior(double *invsigma, double *mean, double *Snf, int ilow, int ihigh)
static double logprior(double *sigma, double *mean, double *Snf, int ilow, int ihigh)
//static double logprior(double *sigma, double *mean, double *Snf, int ilow, int ihigh)
static double logprior(double *lower, double *upper, double *Snf, int ilow, int ihigh)
{
//double x;
double lgp;
......@@ -87,13 +95,77 @@ static double logprior(double *sigma, double *mean, double *Snf, int ilow, int i
lgp -= x*x;
*/
if( fabs(mean[i] - Snf[i]) > 3.*sigma[i]) lgp += -1e60;
// if( fabs(mean[i] - Snf[i]) > 3.*sigma[i]) lgp += -1e60;
if(Snf[i]>upper[i] || Snf[i]<lower[i])
{
//printf("PSD out of range: f=%g, Sn=%g, min=%g, max=%g\n",(double)i/4., Snf[i],lower[i],upper[i]);
lgp=-1e60;
}
}
//return (0.5*lgp);
return (lgp);
}
static double logprior_gaussian_model(double *mean, double *sigma, double *Snf, double *spline_f, int spline_n, double *lines_f, int lines_n, dataParams *data)
{
double x,f;
double lgp;
int i,n;
double flow = data->flow;
double Tobs = data->Tobs;
lgp = 0.0;
//Lorentzian model
for(n=0; n<lines_n; n++)
{
f = lines_f[n];
i = (int)((f-flow)*Tobs);
x = (mean[i] - Snf[i])/sigma[i];
lgp -= x*x;
// printf(" line %i: f=%g, PSD=%g, , P=%g, sigma=%g, logP=%g\n",n,f,Snf[i],mean[i],fabs(x),lgp);
}
//Spline model
for(n=0; n<spline_n; n++)
{
f = spline_f[n];
i = (int)((f-flow)*Tobs);
x = (mean[i] - Snf[i])/sigma[i];
lgp -= x*x;
// printf(" spline %i: f=%g, PSD=%g, , P=%g, sigma=%g, logP=%g\n",n,f,Snf[i],mean[i],fabs(x),lgp);
}
// system_pause();
return (0.5*lgp);
}
static double logprior_gaussian(double *mean, double *sigma, double *Snf, int ilow, int ihigh)
{
double x;
double 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])/sigma[i];
lgp -= x*x;
}
return (0.5*lgp);
}
static double loglike(double *respow, double *Snf, int ncut)
{
double lgl, x;
......@@ -1265,11 +1337,12 @@ void BayesLineNonMarkovianFit(struct BayesLineParams *bayesline, int *nj)
spline->data[j] = mdn;
}
LorentzSplineFit(priors, bayesline->constantLogLFlag, 40000, data, lines_x, spline, sfreq, spow, r);
LorentzSplineFit(priors, bayesline->constantLogLFlag, 400000, data, lines_x, spline, sfreq, spow, r);
//Interpoloate {spoints,sdata} ==> {sfreq,y}
CubicSplineGSL(spline->n,spline->points,spline->data,data->ncut,sfreq,y);
//FILE *temp=fopen("temp2.dat","w");
for(i=0; i< data->ncut; i++)
{
Sbase[i] = exp(y[i]);
......@@ -1287,7 +1360,13 @@ void BayesLineNonMarkovianFit(struct BayesLineParams *bayesline, int *nj)
if(x > deltafmax) z = exp(-(x-deltafmax)/deltafmax);
Sn[i] += z*lines_x->A[j]*frs*frs/(frs*fsq+lines_x->Q[j]*lines_x->Q[j]*(fsq-frs)*(fsq-frs));
}
//at frequencies below fmin output SAmax (killing lower frequencies in any future inner products)
//fprintf(temp,"%lg %lg %lg %lg\n",sfreq[i],spow[i],Sbase[i],Sn[i]);
}
//fclose(temp);
//system_pause();
free(y);
......@@ -1329,7 +1408,7 @@ void BayesLineNonMarkovianFit(struct BayesLineParams *bayesline, int *nj)
}
void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat, int steps, int focus, double *dan)
void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat, int steps, int focus, int priorFlag, double *dan)
{
int nsy, nsx;
double logLx, logLy=0.0, logH;
......@@ -1454,7 +1533,29 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
logLx = 1.0;
// logPsx = logprior(priors->invsigma, priors->mean, Snx, 0, ncut);
logPsx = logprior(priors->sigma, priors->mean, Snx, 0, ncut);
// if(priorFlag)logPsx = logprior(priors->sigma, priors->mean, Snx, 0, ncut);
if(priorFlag==1)
{
logPsx = logprior(priors->lower, priors->upper, Snx, 0, ncut);
// printf("starting position for PSD: %g\n",logPsx);
// if(logPsx < -1e50)
// {
// FILE *temp = fopen("temp.dat","w");
// for(i=0; i<ncut; i++)
// {
// fprintf(temp,"%lg %lg %lg %lg\n",(double)i/4.,priors->lower[i],Snx[i],priors->upper[i]);
// }
// fclose(temp);
// system_pause();
// }
// logPsx = logprior_gaussian_model(priors->mean, priors->sigma, Snx, spline_x->points, spline_x->n, lines_x->f, lines_x->n, data);
}
if(priorFlag==2)
{
logPsx = logprior_gaussian(priors->mean, priors->sigma, Snx, 0, ncut);
}
// set up proposal for frequency jumps
xsm =0.0;
......@@ -1622,7 +1723,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
check = 0;
if(nsy < 1 || nsy > nspline-1) check = 1;
if(nsy < 5 || nsy > nspline-1) check = 1;
for(i=0; i<nsy; i++)
{
......@@ -1956,9 +2057,21 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
logpx += zeta*(double)(lines_x->n);
//logPsy = logprior(priors->invsigma, priors->mean, Sn, 0, ncut);
logPsy = logprior(priors->sigma, priors->mean, Sn, 0, ncut);
//if(priorFlag)logPsy = logprior(priors->sigma, priors->mean, Sn, 0, ncut);
if(priorFlag==1)
{
logPsy = logprior(priors->lower, priors->upper, Sn, 0, ncut);
//logPsy = logprior_gaussian_model(priors->mean, priors->sigma, Sn, spointsy, nsy, lines_y->f, lines_y->n, data);
}
if(priorFlag==2)
{
logPsy = logprior_gaussian(priors->mean, priors->sigma, Sn, 0, ncut);
}
logH = (logLy - logLx)*heat - logpy + logpx + logPsy - logPsx;
logH = (logLy - logLx)*heat - logpy + logpx;
if(priorFlag!=0) logH += logPsy - logPsx;
alpha = log(gsl_rng_uniform(r));
......@@ -1968,6 +2081,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
if(typ == 1) ac1++;
if(typ == 4) ac2++;
logLx = logLy;
//if(priorFlag!=0)
logPsx = logPsy;
lines_x->n = lines_y->n;
nsx = nsy;
......@@ -2151,7 +2265,7 @@ void BayesLineMarkovianSplineOnly(struct BayesLineParams *bayesline, int nspline
}
// produce an initial spline fit to the smooth part of the spectrum
SpecFitSpline(priors, bayesline->constantLogLFlag, 10000, fa, Sna, spline, Snf, jj, r);
SpecFitSpline(priors, bayesline->constantLogLFlag, 100000, fa, Sna, spline, Snf, jj, r);
}
......@@ -2167,15 +2281,15 @@ void BayesLineMarkovianFocusedAnalysis(struct BayesLineParams *bayesline)
int j;
double dan;
BayesLineLorentzSplineMCMC(bayesline, 1.0, 20000, 0, &dan);
BayesLineLorentzSplineMCMC(bayesline, 1.0, 200000, 0, 0, &dan);
//alternate between targeting outliers and MCMCing full solution until outliers are gone or max iterations are reached
j = 0;
do
{
BayesLineLorentzSplineMCMC(bayesline, 1.0, 5000, 1, &dan);
BayesLineLorentzSplineMCMC(bayesline, 1.0, 50000, 1, 0, &dan);
BayesLineLorentzSplineMCMC(bayesline, 1.0, 5000, 0, &dan);
BayesLineLorentzSplineMCMC(bayesline, 1.0, 50000, 0, 0, &dan);
j++;
......@@ -2183,6 +2297,36 @@ void BayesLineMarkovianFocusedAnalysis(struct BayesLineParams *bayesline)
}
void BayesLineFree(struct BayesLineParams *bptr)
{
free(bptr->data);
// storage for line model for a segment. These get updated and passed back from the fitting routine
destroy_lorentzianParams(bptr->lines_x);
destroy_lorentzianParams(bptr->lines_full);
destroy_splineParams(bptr->spline);
destroy_splineParams(bptr->spline_x);
free(bptr->fa);
free(bptr->Sna);
free(bptr->freq);
free(bptr->power);
free(bptr->Snf);
free(bptr->Sbase);
free(bptr->Sline);
free(bptr->spow);
free(bptr->sfreq);
/* set up GSL random number generator */
gsl_rng_free(bptr->r);
free(bptr);
}
void BayesLineSetup(struct BayesLineParams *bptr, double *freqData, double fmin, double fmax, double deltaT, double Tobs)
{
int i;
......@@ -2247,17 +2391,46 @@ void BayesLineSetup(struct BayesLineParams *bptr, double *freqData, double fmin,
bptr->Sna = malloc((size_t)(sizeof(double)*(smodn+1)));
//start with ~4 Hz resolution for the spline
nspline = (int)((bptr->data->fmax-bptr->data->fmin)/bptr->data->fgrid)+2;
create_splineParams(bptr->spline,nspline);
create_splineParams(bptr->spline_x,nspline);
// Re-set dataParams structure to use full dataset
bptr->data->flow = bptr->data->fmin;
bptr->data->fhigh = bptr->data->fmax;
imax = (int)(floor(bptr->data->fhigh*bptr->data->Tobs));
imin = (int)(floor(bptr->data->flow*bptr->data->Tobs));
bptr->data->ncut = imax-imin;
/*
The spow and sfreq arrays hold the frequency and power of the full section of data to be whitened
The model parameters are the number of Lorentzians, nx, and their frequency ff, amplitude AA and
quality factor QQ, as well as the number of terms in the power law fit to the smooth part of the
spectrum, segs, and their amplitudes at 100 Hz, SA, and their spectral slopes SP. The maximum number
of terms in the Lorentzian model is tmax (around 200-300 should cover most spectra), and the
maximum number of terms in the power law fit is segmax (aound 12-20 should be enough). When called
from another MCMC code, the LorentzMCMC code needs to be passed the most recent values for the
noise model. Arrays to hold them have to be declared somewhere in advance. The updated values are
passed back. The first entry in the call to LoretnzMCMC are the number of iterations. It probably
makes sense to do 100 or so each time it is called.
*/
bptr->spow = malloc((size_t)(sizeof(double)*(bptr->data->ncut)));
bptr->sfreq = malloc((size_t)(sizeof(double)*(bptr->data->ncut)));
}
void BayesLineSearch(struct BayesLineParams *bptr, double *freqData, double fmin, double fmax, double deltaT, double Tobs)
{
int i;
int imin, imax;
int imin,imax;
int j;
//start with ~4 Hz resolution for the spline
int nspline = (int)(bptr->data->fstep/bptr->data->fgrid)+2;
int jj = 0;
......@@ -2275,12 +2448,12 @@ void BayesLineSearch(struct BayesLineParams *bptr, double *freqData, double fmin
if(bptr->data->tmax < 20) bptr->data->tmax = 20;
if(bptr->lines_full->n < 1 ) bptr->lines_full->n = 1;
//start with ~4 Hz resolution for the spline
nspline = (int)((bptr->data->fhigh-bptr->data->fmin)/bptr->data->fgrid)+2;
bptr->spline = malloc( (size_t)sizeof(splineParams) );
create_splineParams(bptr->spline,nspline);
// Re-set dataParams structure to use full dataset
bptr->data->flow = bptr->data->fmin;
imax = (int)(floor(bptr->data->fhigh*bptr->data->Tobs));
imin = (int)(floor(bptr->data->flow*bptr->data->Tobs));
bptr->data->ncut = imax-imin;
/******************************************************************************/
/* */
......@@ -2288,43 +2461,17 @@ void BayesLineSearch(struct BayesLineParams *bptr, double *freqData, double fmin
/* */
/******************************************************************************/
BayesLineMarkovianSplineOnly(bptr, nspline, jj);
BayesLineMarkovianSplineOnly(bptr, bptr->spline->n, jj);
//copy spineParams structure spline to spline_x
/*
spline_x will be varied (and updated) by BayesLineLorentzSplineMCMC.
spline contains the original grid used in the RJ proposals
*/
create_splineParams(bptr->spline_x,nspline);
for(j=0; j<nspline; j++)
for(j=0; j<bptr->spline->n; j++)
{
bptr->spline_x->points[j] = bptr->spline->points[j];
bptr->spline_x->data[j] = bptr->spline->data[j];
}
// Re-set dataParams structure to use full dataset
bptr->data->flow = bptr->data->fmin;
imax = (int)(floor(bptr->data->fhigh*bptr->data->Tobs));
imin = (int)(floor(bptr->data->flow*bptr->data->Tobs));
bptr->data->ncut = imax-imin;
/*
The spow and sfreq arrays hold the frequency and power of the full section of data to be whitened
The model parameters are the number of Lorentzians, nx, and their frequency ff, amplitude AA and
quality factor QQ, as well as the number of terms in the power law fit to the smooth part of the
spectrum, segs, and their amplitudes at 100 Hz, SA, and their spectral slopes SP. The maximum number
of terms in the Lorentzian model is tmax (around 200-300 should cover most spectra), and the
maximum number of terms in the power law fit is segmax (aound 12-20 should be enough). When called
from another MCMC code, the LorentzMCMC code needs to be passed the most recent values for the
noise model. Arrays to hold them have to be declared somewhere in advance. The updated values are
passed back. The first entry in the call to LoretnzMCMC are the number of iterations. It probably
makes sense to do 100 or so each time it is called.
*/
bptr->spow = malloc((size_t)(sizeof(double)*(bptr->data->ncut)));
bptr->sfreq = malloc((size_t)(sizeof(double)*(bptr->data->ncut)));
imin = (int)(floor(bptr->data->flow*bptr->data->Tobs));
for(i=0; i< bptr->data->ncut; i++)
{
......@@ -2341,12 +2488,13 @@ void BayesLineSearch(struct BayesLineParams *bptr, double *freqData, double fmin
BayesLineMarkovianFocusedAnalysis(bptr);
double dan;
BayesLineLorentzSplineMCMC(bptr, 1.0, 20000, 0, &dan);
BayesLineLorentzSplineMCMC(bptr, 1.0, 200000, 0, 0, &dan);
}
void BayesLineRJMCMC(struct BayesLineParams *bayesline, double *freqData, double *psd, double *invpsd, double *splinePSD, int N, int cycle, double beta)
void BayesLineRJMCMC(struct BayesLineParams *bayesline, double *freqData, double *psd, double *invpsd, double *splinePSD, int N, int cycle, double beta, int priorFlag)
{
int i,j;
int imin,imax;
......@@ -2354,8 +2502,8 @@ void BayesLineRJMCMC(struct BayesLineParams *bayesline, double *freqData, double
double dan;
//at frequencies below fmin output SAmax (killing lower frequencies in any future inner products)
imax = (int)(floor(bayesline->data->fhigh * bayesline->data->Tobs));
imin = (int)(floor(bayesline->data->flow * bayesline->data->Tobs));
imax = (int)(floor(bayesline->data->fmax * bayesline->data->Tobs));
imin = (int)(floor(bayesline->data->fmin * bayesline->data->Tobs));
for(i=0; i< bayesline->data->ncut; i++)
{
......@@ -2367,7 +2515,7 @@ void BayesLineRJMCMC(struct BayesLineParams *bayesline, double *freqData, double
bayesline->spow[i] = (reFreq*reFreq+imFreq*imFreq);
}
BayesLineLorentzSplineMCMC(bayesline, beta, cycle, 0, &dan);
BayesLineLorentzSplineMCMC(bayesline, beta, cycle, 0, priorFlag, &dan);
/******************************************************************************/
/* */
......@@ -2394,40 +2542,25 @@ void BayesLineRJMCMC(struct BayesLineParams *bayesline, double *freqData, double
void copy_bayesline_params(struct BayesLineParams *origin, struct BayesLineParams *copy)
{
/*
dataParams *data;
splineParams *spline;
splineParams *spline_x;
lorentzianParams *lines_x;
lorentzianParams *lines_full;
BayesLinePriors *priors;
double *Snf;
double *Sna;
double *fa;
double *freq;
double *power;
double *spow;
double *sfreq;
double *Sbase;
double *Sline;
int constantLogLFlag;
double TwoDeltaT;
gsl_rng *r;
FILE *splineChainFile;
FILE *lineChainFile;
*/
copy->data->tmax = origin->data->tmax;
int i;
for(i=0; i<origin->data->ncut; i++)
{
copy->spow[i] = origin->spow[i];
copy->sfreq[i] = origin->sfreq[i];
}
// copy->spline = malloc( (size_t)sizeof(splineParams) );
// copy->spline_x = malloc( (size_t)sizeof(splineParams) );
// create_splineParams(copy->spline,origin->spline->n);
// create_splineParams(copy->spline_x,origin->spline_x->n);
//
// copy->spow = malloc((size_t)(sizeof(double)*(origin->data->ncut)));
// copy->sfreq = malloc((size_t)(sizeof(double)*(origin->data->ncut)));
int imax = (int)(floor(origin->data->fmax * origin->data->Tobs));
int imin = (int)(floor(origin->data->fmin * origin->data->Tobs));
for(i=0; i<imax-imin; i++)
{
copy->Snf[i] = origin->Snf[i];
copy->spow[i] = origin->spow[i];
copy->sfreq[i] = origin->sfreq[i];
copy->Sbase[i] = origin->Sbase[i];
copy->Sline[i] = origin->Sline[i];
}
copy_splineParams(origin->spline, copy->spline);
copy_splineParams(origin->spline_x,copy->spline_x);
......@@ -2441,7 +2574,7 @@ void print_line_model(FILE *fptr, struct BayesLineParams *bayesline)
int j;
fprintf(fptr,"%i ", bayesline->lines_x->n);
for(j=0; j< bayesline->lines_full->n; j++) fprintf(fptr,"%e %e %e ", bayesline->lines_full->f[ bayesline->lines_full->larray[j]], bayesline->lines_full->A[ bayesline->lines_full->larray[j]], bayesline->lines_full->Q[ bayesline->lines_full->larray[j]]);
for(j=0; j< bayesline->lines_full->n; j++) fprintf(fptr,"%e %e %e ", bayesline->lines_full->f[j], bayesline->lines_full->A[j], bayesline->lines_full->Q[j]);
fprintf(fptr,"\n");
}
......@@ -2454,4 +2587,19 @@ void print_spline_model(FILE *fptr, struct BayesLineParams *bayesline)
fprintf(fptr,"\n");
}
void parse_line_model(FILE *fptr, struct BayesLineParams *bayesline)
{
int j;
fscanf(fptr,"%i",&bayesline->lines_x->n);
for(j=0; j< bayesline->lines_full->n; j++) fscanf(fptr,"%lg %lg %lg",&bayesline->lines_full->f[j],&bayesline->lines_full->A[j],&bayesline->lines_full->Q[j]);
}
void parse_spline_model(FILE *fptr, struct BayesLineParams *bayesline)
{
int j;
fscanf(fptr,"%i",&bayesline->spline_x->n);
for(j=0; j<bayesline->spline_x->n; j++)fscanf(fptr,"%lg %lg",&bayesline->spline_x->points[j],&bayesline->spline_x->data[j]);
}
......@@ -77,6 +77,8 @@ typedef struct
//double *invsigma; //variances for each frequency bin
double *sigma; //variances for each frequency bin
double *upper; //variances for each frequency bin
double *lower; //variances for each frequency bin
double *mean; //means for each frequency bin
}BayesLinePriors;
......@@ -106,13 +108,14 @@ struct BayesLineParams
gsl_rng *r;
};
void BayesLineFree(struct BayesLineParams *bptr);
void BayesLineSetup(struct BayesLineParams *bptr, double *freqData, double fmin, double fmax, double deltaT, double Tobs);
void BayesLineRJMCMC(struct BayesLineParams *bayesline, double *freqData, double *psd, double *invpsd, double *splinePSD, int N, int cycle, double beta);
void BayesLineRJMCMC(struct BayesLineParams *bayesline, double *freqData, double *psd, double *invpsd, double *splinePSD, int N, int cycle, double beta, int priorFlag);
void BayesLineSearch(struct BayesLineParams *bptr, double *freqData, double fmin, double fmax, double deltaT, double Tobs);
void BayesLineNonMarkovianFit (struct BayesLineParams *bayesline, int *nj);
void BayesLineLorentzSplineMCMC (struct BayesLineParams *bayesline, double heat, int steps, int focus, double *dan);
void BayesLineLorentzSplineMCMC (struct BayesLineParams *bayesline, double heat, int steps, int focus, int priorFlag, double *dan);
void BayesLineMarkovianSplineOnly (struct BayesLineParams *bayesline, int nspline, int jj);
void BayesLineMarkovianFocusedAnalysis (struct BayesLineParams *bayesline);
......@@ -148,6 +151,8 @@ void destroy_splineParams(splineParams *spline);
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);
void parse_line_model(FILE *fptr, struct BayesLineParams *bayesline);
void parse_spline_model(FILE *fptr, struct BayesLineParams *bayesline);
......@@ -33,7 +33,7 @@ int main(int argc, char *argv[])
{
/* Variable declaration */
int i, ic, ifo, imin;
int i, j, ic, ifo, imin,imax;
char filename[100];
......@@ -160,6 +160,7 @@ int main(int argc, char *argv[])
//Copy Fourier domain data & PSD from LAL data types into simple arrays
ifo=0;
imin = (int)(fmin*Tobs);
imax = (int)(fmax*Tobs);
while(dataPtr!=NULL)
{
for(i=0; i<N/2; i++)
......@@ -242,6 +243,10 @@ int main(int argc, char *argv[])
/*
Setup CHAIN structure
*/
//give chain buffer in case we need more
int NC=chain->NC;
chain->NC+=5;
allocate_chain(chain,NI);
for(ic=0; ic<chain->NC; ic++) chain->index[ic]=ic;
......@@ -325,69 +330,87 @@ int main(int argc, char *argv[])
for(ic=0; ic<chain->NC; ic++)
{
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));
//Set BayesLine priors based on the data channel being used
set_bayesline_priors(data->channels[ifo], bayesline[ic][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[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);
}
initialize_bayesline(bayesline[ic], data, psd);
}
/*
BayesLine spectral estimation search-phase
*/