From bca1d22d0e150ed465814b10403ac58d2fd72689 Mon Sep 17 00:00:00 2001
From: Tyson Littenberg <tyson.littenberg@gmail.com>
Date: Thu, 31 Jan 2019 22:55:31 -0600
Subject: [PATCH] A lot: +fixed mapping from burn in to psd parameters. +Now
apply PSD priors during cleaning phase +Allow user to skip cleaning phase
when using BayesLine +Tweaks to animation script +Undo change that reduced
number of samples in cleaning phase
---
src/BayesLine.c | 309 +++++++++++++++++++++-----------------------
src/BayesLine.h | 3 +-
src/BayesWave.c | 24 +---
src/BayesWaveIO.c | 19 ++-
src/BayesWaveMCMC.c | 18 ---
5 files changed, 166 insertions(+), 207 deletions(-)
diff --git a/src/BayesLine.c b/src/BayesLine.c
index 3e487128..3be82bb1 100644
--- a/src/BayesLine.c
+++ b/src/BayesLine.c
@@ -91,21 +91,11 @@ static double logprior(double *lower, double *upper, double *Snf, int ilow, int
lgp = 0;
for(i=ilow; i<ihigh; i++)
{
- /*
- x = (mean[i] - Snf[i])*invsigma[i];
- lgp -= x*x;
- */
-
-// 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;
+ lgp=-1e60;
}
}
-
- //return (0.5*lgp);
return (lgp);
}
@@ -564,7 +554,7 @@ void destroy_splineParams(splineParams *spline)
void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat, int steps, int focus, int priorFlag, double *dan)
{
- int nsy, nsx;
+ int nsy;
double logLx, logLy=0.0, logH;
int ilowx, ihighx, ilowy, ihighy;
int i, j, k=0, ki=0, ii=0, jj=0, mc;
@@ -593,7 +583,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
/* Make local pointers to BayesLineParams structure members */
dataParams *data = bayesline->data;
- lorentzianParams *lines_x = bayesline->lines_full;
+ lorentzianParams *lines_x = bayesline->lines_x;
splineParams *spline = bayesline->spline;
splineParams *spline_x = bayesline->spline_x;
BayesLinePriors *priors = bayesline->priors;
@@ -609,7 +599,6 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
double fhigh = data->fhigh;
int nspline = spline->n;
- int *nsplinex = &spline_x->n;
double *sdatax = spline_x->data;
@@ -645,8 +634,6 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
lSAmin = log(priors->SAmin);
lSAmax = log(priors->SAmax);
- nsx = *nsplinex;
-
dff = 0.01; // half-width of frequency focus region (used if focus == 1)
// this is the fractional error estimate on the noise level
@@ -654,16 +641,20 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
s3 = 0.5;
s4 = 0.5;
- for(i=0; i<lines_x->n; i++)
- {
- lines_x->larray[i] = i;
- lines_y->larray[i] = i;
- }
+ for(i=0; i<lines_x->n; i++)
+ {
+ lines_x->larray[i] = i;
+ lines_y->larray[i] = i;
+ }
+
+// FILE *temp=fopen("start_psd.dat","w");
+// for(i=0; i<ncut; i++) fprintf(temp,"%i %lg %lg %lg\n",i,Snf[i],priors->lower[i],priors->upper[i]);
+// fclose(temp);
baseav = 0.0;
//Interpolate {spointsx,sdatax} ==> {freq,xint}
- CubicSplineGSL(nsx,spointsx,sdatax,ncut,freq,xint);
+ CubicSplineGSL(spline_x->n,spointsx,sdatax,ncut,freq,xint);
for(i=0; i<ncut; i++)
{
@@ -677,7 +668,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
full_spectrum_spline(Sline, Sbase, freq, data, lines_x);
for(i=0; i< ncut; i++) Sn[i] = Sbase[i]+Sline[i];
- for(i=0; i<ncut; i++) Snx[i] = Sn[i];
+ for(i=0; i<ncut; i++) Snx[i] = Sn[i];
if(!bayesline->constantLogLFlag)
logLx = loglike(power, Sn, ncut);
@@ -692,8 +683,15 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
{
logPsx = logprior_gaussian(priors->mean, priors->sigma, Snx, 0, ncut);
}
-
-
+
+ if(logPsx<0)
+ {
+ printf("how did PSD come in outside of the prior?\n");
+ FILE *temp=fopen("failed_psd.dat","w");
+ for(i=0; i<ncut; i++) fprintf(temp,"%i %lg %lg %lg\n",i,Snx[i],priors->lower[i],priors->upper[i]);
+
+ exit(1);
+ }
// set up proposal for frequency jumps
xsm =0.0;
@@ -732,7 +730,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
//copy over current state
lines_y->n = lines_x->n;
- nsy = nsx;
+ nsy = spline_x->n;
for(i=0; i< tmax; i++)
{
lines_y->larray[i] = lines_x->larray[i];
@@ -760,14 +758,14 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
//decide between adding or removing spline point
alpha = gsl_rng_uniform(r);
- if(alpha > 0.5)// || nsx<3) // try and add a new term
+ if(alpha > 0.5)// || spline_x->n<3) // try and add a new term
{
- nsy = nsx+1;
+ nsy = spline_x->n+1;
typ = 5;
}
else // try and remove term
{
- nsy = nsx-1;
+ nsy = spline_x->n-1;
typ = 6;
}
@@ -775,12 +773,12 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
if(nsy > 0 && nsy < nspline)
{
//Spline death move
- if(nsy < nsx)
+ if(nsy < spline_x->n)
{
- ki=1+(int)(gsl_rng_uniform(r)*(double)(nsx-2)); // pick a term to try and kill - cant be first or last term
+ ki=1+(int)(gsl_rng_uniform(r)*(double)(spline_x->n-2)); // pick a term to try and kill - cant be first or last term
k = 0;
- for(j=0; j<nsx; j++)
+ for(j=0; j<spline_x->n; j++)
{
if(j != ki)
{
@@ -793,7 +791,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
}//end death moove
//Spline birth move
- if(nsy > nsx)
+ if(nsy > spline_x->n)
{
// have to randomly pick a new point that isn't already in use
@@ -803,7 +801,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
ii = 0;
// printf(" try adding point %i\n",ki);
- for(j=0; j<nsx; j++)
+ for(j=0; j<spline_x->n; j++)
{
if(fabs(spointsx[j]-spoints[ki]) < 1.0e-3) ii = 1; // this means the point is already in use
}
@@ -811,7 +809,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
ii = 0;
//slot new point into live array
- for(j=0; j<nsx; j++)
+ for(j=0; j<spline_x->n; j++)
{
if(spointsx[j] < spoints[ki])
{
@@ -840,10 +838,10 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
{
typ = 4;
- nsy = nsx;
+ nsy = spline_x->n;
//pick a term to update
- ii=(int)(gsl_rng_uniform(r)*(double)(nsx));
+ ii=(int)(gsl_rng_uniform(r)*(double)(spline_x->n));
// use a variety of jump sizes by using a sum of gaussians of different width
e1 = 0.0005;
@@ -1152,48 +1150,45 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
if(typ == 1) cc1++;
if(typ == 4) cc2++;
- if(!bayesline->constantLogLFlag)
+ if(typ > 3) // need to do a full update (slower)
{
- if(typ > 3) // need to do a full update (slower)
+ if(nsy>2)
{
- if(nsy>2)
- {
- //Interpolate {spointsy,sdatay} ==> {freq,xint}
- CubicSplineGSL(nsy,spointsy,sdatay,ncut,freq,xint);
+ //Interpolate {spointsy,sdatay} ==> {freq,xint}
+ CubicSplineGSL(nsy,spointsy,sdatay,ncut,freq,xint);
- for(i=0; i<ncut; i++) Sbase[i] = exp(xint[i]);
+ for(i=0; i<ncut; i++) Sbase[i] = exp(xint[i]);
- full_spectrum_spline(Sline, Sbase, freq, data, lines_y);
- for(i=0; i< ncut; i++) Sn[i] = Sbase[i]+Sline[i];
- logLy = loglike(power, Sn, ncut);
- }
- else logLy = -1e60;
- }
-
- if(typ == 1 || typ == 0) // fixed dimension MCMC of line ii
- {
- full_spectrum_single(Sn, Snx, Sbasex, freq, data, lines_x, lines_y, ii, &ilowx, &ihighx, &ilowy, &ihighy);
- logLy = logLx + loglike_single(power, Sn, Snx, ilowx, ihighx, ilowy, ihighy);
- }
-
- if(typ == 2) // add new line with label ii
- {
- full_spectrum_add_or_subtract(Sn, Snx, Sbasex, freq, data, lines_y, ii, &ilowy, &ihighy, 1);
- logLy = logLx + loglike_pm(power, Sn, Snx, ilowy, ihighy);
+ full_spectrum_spline(Sline, Sbase, freq, data, lines_y);
+ for(i=0; i< ncut; i++) Sn[i] = Sbase[i]+Sline[i];
+ logLy = loglike(power, Sn, ncut);
+ if(bayesline->constantLogLFlag) logLy = 1.0;
}
+ else logLy = -1e60;
+ }
- if(typ == 3) // remove line with label ki
- {
- full_spectrum_add_or_subtract(Sn, Snx, Sbasex, freq, data, lines_x, ki, &ilowx, &ihighx, -1);
- logLy = logLx + loglike_pm(power, Sn, Snx, ilowx, ihighx);
- }
+ if(typ == 1 || typ == 0) // fixed dimension MCMC of line ii
+ {
+ full_spectrum_single(Sn, Snx, Sbasex, freq, data, lines_x, lines_y, ii, &ilowx, &ihighx, &ilowy, &ihighy);
+ logLy = logLx + loglike_single(power, Sn, Snx, ilowx, ihighx, ilowy, ihighy);
+ if(bayesline->constantLogLFlag) logLy = 1.0;
+ }
+ if(typ == 2) // add new line with label ii
+ {
+ full_spectrum_add_or_subtract(Sn, Snx, Sbasex, freq, data, lines_y, ii, &ilowy, &ihighy, 1);
+ logLy = logLx + loglike_pm(power, Sn, Snx, ilowy, ihighy);
+ if(bayesline->constantLogLFlag) logLy = 1.0;
}
- else
+
+ if(typ == 3) // remove line with label ki
{
- logLy = 1.0;
+ full_spectrum_add_or_subtract(Sn, Snx, Sbasex, freq, data, lines_x, ki, &ilowx, &ihighx, -1);
+ logLy = logLx + loglike_pm(power, Sn, Snx, ilowx, ihighx);
+ if(bayesline->constantLogLFlag) logLy = 1.0;
}
+
// prior on line number e(-ZETA * n). (this is a prior, not a proposal, so opposite sign)
// effectively an SNR cut on lines
@@ -1220,7 +1215,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
alpha = log(gsl_rng_uniform(r));
if(logH > alpha)
- {
+ {
if(typ == 0) ac0++;
if(typ == 1) ac1++;
if(typ == 4) ac2++;
@@ -1228,7 +1223,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
//if(priorFlag!=0)
logPsx = logPsy;
lines_x->n = lines_y->n;
- nsx = nsy;
+ spline_x->n = nsy;
for(i=0; i< ncut; i++)
{
Snx[i] = Sn[i];
@@ -1241,7 +1236,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
lines_x->f[i] = lines_y->f[i];
lines_x->Q[i] = lines_y->Q[i];
}
- for(i=0; i<nspline; i++)
+ for(i=0; i<spline_x->n; i++)
{
sdatax[i] = sdatay[i];
spointsx[i] = spointsy[i];
@@ -1250,59 +1245,59 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
}//end prior check
- //Every 1000 steps update focus region
- if(mc%1000 == 0)
- {
- pmax = -1.0;
- for(i=0; i< ncut; i++)
- {
- x = power[i]/Snx[i];
- if(x > pmax)
- {
- pmax = x;
- k = i;
- }
- }
-
- // define the focus region (only used if focus flag = 1)
- fcl = freq[k]-dff;
- fch = freq[k]+dff;
- Ac = power[k];
- if(fcl < freq[0]) fcl = freq[0];
- if(fch > freq[ncut-1]) fch = freq[ncut-1];
-
- //if(focus==1)fprintf(stdout,"Focusing on [%f %f] Hz...\n", fcl, fch);
-
-
- // set up proposal for frequency jumps
- xsm =0.0;
- baseav = 0.0;
- for(i=0; i< ncut; i++)
- {
- x = power[i]/Snx[i];
- if(x < 16.0) x = 1.0;
- if(x >= 16.0) x = 100.0;
- fprop[i] = x;
- xsm += x;
- baseav += log(Sbasex[i]);
- }
- baseav /= (double)(ncut);
-
-
- pmax = -1.0;
- for(i=0; i< ncut; i++)
- {
- fprop[i] /= xsm;
- if(fprop[i] > pmax) pmax = fprop[i];
- }
-
- }//end focus update
-
+ //Every 1000 steps update focus region
+ if(mc%1000 == 0)
+ {
+ pmax = -1.0;
+ for(i=0; i< ncut; i++)
+ {
+ x = power[i]/Snx[i];
+ if(x > pmax)
+ {
+ pmax = x;
+ k = i;
+ }
+ }
+
+ // define the focus region (only used if focus flag = 1)
+ fcl = freq[k]-dff;
+ fch = freq[k]+dff;
+ Ac = power[k];
+ if(fcl < freq[0]) fcl = freq[0];
+ if(fch > freq[ncut-1]) fch = freq[ncut-1];
+
+ //if(focus==1)fprintf(stdout,"Focusing on [%f %f] Hz...\n", fcl, fch);
+
+
+ // set up proposal for frequency jumps
+ xsm =0.0;
+ baseav = 0.0;
+ for(i=0; i< ncut; i++)
+ {
+ x = power[i]/Snx[i];
+ if(x < 16.0) x = 1.0;
+ if(x >= 16.0) x = 100.0;
+ fprop[i] = x;
+ xsm += x;
+ baseav += log(Sbasex[i]);
+ }
+ baseav /= (double)(ncut);
+
+
+ pmax = -1.0;
+ for(i=0; i< ncut; i++)
+ {
+ fprop[i] /= xsm;
+ if(fprop[i] > pmax) pmax = fprop[i];
+ }
+
+ }//end focus update
+
}//End MCMC loop
//Interpolate {spointsx,sdatax} ==> {freq,xint}
- CubicSplineGSL(nsx,spointsx,sdatax,ncut,freq,xint);
+ CubicSplineGSL(spline_x->n,spointsx,sdatax,ncut,freq,xint);
for(i=0; i< ncut; i++) Sbase[i] = exp(xint[i]);
full_spectrum_spline(Sline, Sbase, freq, data, lines_x);
@@ -1315,7 +1310,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
// return updated spectral model
for(i=0; i< ncut; i++) Snf[i] = Snx[i];
-
+
// re-map the array to 0..mx ordering
for(i=0; i< lines_x->n; i++)
{
@@ -1342,8 +1337,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
}
*dan = pmax;
- *nsplinex = nsx;
-
+
free(foc);
free(xint);
free(fprop);
@@ -1365,7 +1359,6 @@ void BayesLineFree(struct BayesLineParams *bptr)
// 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);
@@ -1400,7 +1393,6 @@ void BayesLineSetup(struct BayesLineParams *bptr, double *freqData, double fmin,
bptr->data = malloc(sizeof(dataParams));
bptr->lines_x = malloc(sizeof(lorentzianParams));
- bptr->lines_full = malloc(sizeof(lorentzianParams));
bptr->spline = malloc(sizeof(splineParams));
bptr->spline_x = malloc(sizeof(splineParams));
@@ -1442,7 +1434,7 @@ void BayesLineSetup(struct BayesLineParams *bptr, double *freqData, double fmin,
create_lorentzianParams(bptr->lines_x,bptr->data->tmax);
// storage for master line model
- create_lorentzianParams(bptr->lines_full,bptr->data->tmax);
+ create_lorentzianParams(bptr->lines_x,bptr->data->tmax);
// start with a single line. This number gets updated and passed back
bptr->lines_x->n = 1;
@@ -1509,7 +1501,7 @@ void BayesLineInitialize(struct BayesLineParams *bayesline)
bayesline->data->tmax = 1000;
if(bayesline->data->tmax < 20) bayesline->data->tmax = 20;
- if(bayesline->lines_full->n < 1 ) bayesline->lines_full->n = 1;
+ if(bayesline->lines_x->n < 1 ) bayesline->lines_x->n = 1;
// Re-set dataParams structure to use full dataset
@@ -1651,7 +1643,6 @@ void copy_bayesline_params(struct BayesLineParams *origin, struct BayesLineParam
copy_splineParams(origin->spline, copy->spline);
copy_splineParams(origin->spline_x,copy->spline_x);
copy_lorentzianParams(origin->lines_x, copy->lines_x);
- copy_lorentzianParams(origin->lines_full,copy->lines_full);
}
@@ -1659,8 +1650,8 @@ void print_line_model(FILE *fptr, struct BayesLineParams *bayesline)
{
int j;
- fprintf(fptr,"%i ", bayesline->lines_full->n);
- 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,"%i ", bayesline->lines_x->n);
+ for(j=0; j< bayesline->lines_x->n; j++) fprintf(fptr,"%e %e %e ", bayesline->lines_x->f[j], bayesline->lines_x->A[j], bayesline->lines_x->Q[j]);
fprintf(fptr,"\n");
}
@@ -1677,8 +1668,8 @@ void parse_line_model(FILE *fptr, struct BayesLineParams *bayesline)
{
int j;
- fscanf(fptr,"%i",&bayesline->lines_full->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]);
+ fscanf(fptr,"%i",&bayesline->lines_x->n);
+ for(j=0; j< bayesline->lines_x->n; j++) fscanf(fptr,"%lg %lg %lg",&bayesline->lines_x->f[j],&bayesline->lines_x->A[j],&bayesline->lines_x->Q[j]);
}
void parse_spline_model(FILE *fptr, struct BayesLineParams *bayesline)
{
@@ -2584,12 +2575,12 @@ void blstart(double *data, double *residual, int N, double dt, double fmin, int
if(x < 5.0 || (x < xold && xnext > x))
{
flag = 0;
- j++;
linew[j] = (double)(k)/Tobs;
linef[j] = (double)(ii)/Tobs;
lineh[j] = (max-1.0)*sspecD[ii]*fac;
Abar = 0.5*(specD[ii+1]/sspecD[ii+1]+specD[ii-1]/sspecD[ii-1]);
lineQ[j] = sqrt((max/Abar-1.0))*linef[j]*Tobs;
+ j++;
//printf("%d %e %e %e\n", j, linef[j], linew[j], lineQ[j]);
}
}
@@ -2598,7 +2589,7 @@ void blstart(double *data, double *residual, int N, double dt, double fmin, int
xold = x;
}
- Nlines = j;
+ Nlines = j+1;
// printf("There are %d lines\n", Nlines);
@@ -2737,7 +2728,7 @@ void BayesLineBurnin(struct BayesLineParams *bayesline, double *timeData, double
//shortcuts to members of BayesLine structure
dataParams *data = bayesline->data;
- lorentzianParams *lines = bayesline->lines_full;
+ lorentzianParams *lines = bayesline->lines_x;
splineParams *spline = bayesline->spline_x;
int i,j;
@@ -2768,42 +2759,42 @@ void BayesLineBurnin(struct BayesLineParams *bayesline, double *timeData, double
bayesline->Sline[i] = 0.0;
}
- //interpolate cubic spline points
- CubicSplineGSL(spline->n,spline->points,spline->data,N/2-imin,f+imin,logSn+imin);
+ //interpolate cubic spline points
+ CubicSplineGSL(spline->n,spline->points,spline->data,N/2-imin,f+imin,logSn+imin);
//initialize work space for spline model
- bayesline->data->flow = bayesline->data->fmin;
-
- for(i=0; i< bayesline->data->ncut; i++)
- {
- j = i + imin - bayesline->data->nmin;
+ bayesline->data->flow = bayesline->data->fmin;
+
+ for(i=0; i< bayesline->data->ncut; i++)
+ {
+ j = i + imin - bayesline->data->nmin;
bayesline->power[j] = freqData[2*j]*freqData[2*j]+freqData[2*j+1]*freqData[2*j+1];
bayesline->spow[i] = bayesline->power[j];
- bayesline->sfreq[i] = bayesline->freq[j];
- }
+ bayesline->sfreq[i] = bayesline->freq[j];
+ }
- //form-up line model
- for(i=0; i<lines->n; i++)
- {
- lines->larray[i] = i;
- Lorentzian(Sl, Tobs, lines, i, N/2);
- }
+ //form-up line model
+ for(i=0; i<lines->n; i++)
+ {
+ lines->larray[i] = i;
+ Lorentzian(Sl, Tobs, lines, lines->larray[i], N);
+ }
- //combine spline and line model
- for(i=0; i<N/2-imin; i++)
- {
- bayesline->Sbase[i] = exp(logSn[i+imin]);
- bayesline->Sline[i] = Sl[i+imin];
- bayesline->Snf[i] = bayesline->Sbase[i] + bayesline->Sline[i];
- }
+ //combine spline and line model
+ for(i=0; i<N/2-imin; i++)
+ {
+ bayesline->Sbase[i] = exp(logSn[i+imin]);
+ bayesline->Sline[i] = Sl[i+imin];
+ bayesline->Snf[i] = bayesline->Sbase[i] + bayesline->Sline[i];
+ }
//Use initial estimate of PSD to set priors
for(i=0; i<N/2-imin; i++)
{
bayesline->priors->sigma[i] = bayesline->Snf[i];
bayesline->priors->mean[i] = bayesline->Snf[i];
- bayesline->priors->lower[i] = bayesline->Snf[i]/10.;
- bayesline->priors->upper[i] = bayesline->Snf[i]*10.;
+ bayesline->priors->lower[i] = bayesline->Sbase[i]/100.;
+ bayesline->priors->upper[i] = bayesline->Snf[i]*100.;
}
free(f);
diff --git a/src/BayesLine.h b/src/BayesLine.h
index eec86fd5..32553a6c 100644
--- a/src/BayesLine.h
+++ b/src/BayesLine.h
@@ -79,8 +79,7 @@ struct BayesLineParams
dataParams *data;
splineParams *spline;
splineParams *spline_x;
- lorentzianParams *lines_x;
- lorentzianParams *lines_full;
+ lorentzianParams *lines_x;
BayesLinePriors *priors;
double *Snf;
diff --git a/src/BayesWave.c b/src/BayesWave.c
index d7fbbf39..728f4955 100644
--- a/src/BayesWave.c
+++ b/src/BayesWave.c
@@ -380,6 +380,11 @@ int main(int argc, char *argv[])
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];
+ }
+ for(i=0; i< bayesline[chain->index[0]][ifo]->data->ncut; i++)
+ {
+ bayesline[chain->index[ic]][ifo]->priors->lower[i] = bayesline[chain->index[0]][ifo]->priors->lower[i];
+ bayesline[chain->index[ic]][ifo]->priors->upper[i] = bayesline[chain->index[0]][ifo]->priors->upper[i];
}
}
}
@@ -527,8 +532,8 @@ int main(int argc, char *argv[])
int CP = data->clusterPriorFlag;
data->clusterPriorFlag = 0;
chain->count = M;
- //chain->NC = NC/5;
- //if(chain->NC < 1)chain->NC=1;
+ chain->NC = NC/5;
+ if(chain->NC < 1)chain->NC=1;
//Initialize proposal ratios for glitch cleaning phase
chain->modelRate = 0.0; //even chance for each channel
@@ -557,12 +562,9 @@ int main(int argc, char *argv[])
printf("characterizing PSD+glitch model...\n\n");
//TODO disable checkpointing and adaptation during cleaning phase
- //int saveCheckpointFlag = data->checkpointFlag;
int saveAdaptationFlag = data->adaptTemperatureFlag;
- //data->checkpointFlag = 0;
data->adaptTemperatureFlag = 0;
RJMCMC(data, model, bayesline, chain, prior, &cleanLogZ, &cleanVarZ);
- //data->checkpointFlag = saveCheckpointFlag;
data->adaptTemperatureFlag = saveAdaptationFlag;
/******************************************************************************/
@@ -673,18 +675,6 @@ int main(int argc, char *argv[])
for(ifo=0; ifo<NI; ifo++)
{
for(i=0; i<data->N/2; i++) psd[ifo][i] = model[chain->index[0]]->Snf[ifo][i];
-
- //Reset priors on PSD
- imin = (int)(data->Tobs*data->fmin);
- for(i=0; i<(int)(data->Tobs*(data->fmax-data->fmin))+1; i++)
- {
- //bayesline[ifo]->priors->invsigma[i] = 1./(psd[ifo][i+(int)(Tobs*fmin)]*1.0);
- bayesline[chain->index[0]][ifo]->priors->sigma[i] = psd[ifo][i+imin];
- bayesline[chain->index[0]][ifo]->priors->mean[i] = psd[ifo][i+imin];
- bayesline[chain->index[0]][ifo]->priors->lower[i] = psd[ifo][i+imin]/10.;
- bayesline[chain->index[0]][ifo]->priors->upper[i] = psd[ifo][i+imin]*2.;
- }
-
for(ic=1; ic<chain->NC; ic++)
{
diff --git a/src/BayesWaveIO.c b/src/BayesWaveIO.c
index bbbbf6be..9210255b 100644
--- a/src/BayesWaveIO.c
+++ b/src/BayesWaveIO.c
@@ -710,7 +710,6 @@ void write_gnuplot_psd_animation(struct Data *data, char modelname[], int frame)
fprintf(framefile,"set logscale\n");
fprintf(framefile,"set tics scale 2 nomirror\n");
fprintf(framefile,"set border 3\n");
- fprintf(framefile,"unset key\n");
fprintf(framefile,"\n");
fprintf(framefile,"psdmin = 1e-48\n");
fprintf(framefile,"psdmax = 1e-38\n");
@@ -734,15 +733,15 @@ void write_gnuplot_psd_animation(struct Data *data, char modelname[], int frame)
fprintf(framefile," set format \"10^{%%T}\"\n");
fprintf(framefile," set xlabel 'f'\n");
fprintf(framefile," set ylabel '|h|^2'\n");
- fprintf(framefile," plot [*:*] [psdmin:psdmax] '%s_psd_%s_'.N.'.dat' u 1:2 w l lc 'gray', '' u 1:($3+$4) w l lc 'black' lw 2\n", modelname,ifos[n]);
+ fprintf(framefile," plot [*:*] [psdmin:psdmax] '%s_psd_%s_'.N.'.dat' u 1:2 w l lc 'gray' notitle, '' u 1:($3+$4) w l lc 'black' lw 2 notitle\n", modelname,ifos[n]);
fprintf(framefile,"\n");
fprintf(framefile," unset logscale\n");
fprintf(framefile," set format '%%g'\n");
fprintf(framefile," set xlabel 't'\n");
fprintf(framefile," set ylabel 'h'\n");
- fprintf(framefile," plot [tmin:tmax] [hmin:hmax] '%s_data_%s_'.N.'.dat' w l lc 'gray'",modelname,ifos[n]);
- if(data->glitchFlag) fprintf(framefile,", '%s_glitch_%s_'.N.'.dat' w l lc rgb '#a6611a' lw 2",modelname,ifos[n]);
- if(data->signalFlag) fprintf(framefile,", '%s_wave_%s_'.N.'.dat' w l lc rgb '#7b3294' lw 2",modelname,ifos[n]);
+ fprintf(framefile," plot [tmin:tmax] [hmin:hmax] '%s_data_%s_'.N.'.dat' w l lc 'gray' notitle",modelname,ifos[n]);
+ if(data->glitchFlag) fprintf(framefile,", '%s_glitch_%s_'.N.'.dat' w l lc rgb '#a6611a' lw 2 title 'glitch'",modelname,ifos[n]);
+ if(data->signalFlag) fprintf(framefile,", '%s_wave_%s_'.N.'.dat' w l lc rgb '#7b3294' lw 2 title 'signal'",modelname,ifos[n]);
fprintf(framefile,"\n");
fprintf(framefile,"\n");
}
@@ -751,7 +750,7 @@ void write_gnuplot_psd_animation(struct Data *data, char modelname[], int frame)
fprintf(framefile," unset multiplot\n");
fprintf(framefile,"}\n");
fprintf(framefile,"\n");
- fprintf(framefile,"CMD='ffmpeg -r 4 -f image2 -s '.xsize*rescale.'x'.rescale*nifo*ysize.' -start_number 0 -i animation/bayeswave%%03d.png -vframes '.nmax.' -vcodec libx264 -crf 25 -pix_fmt yuv420p %s'\n",filename);
+ fprintf(framefile,"CMD='ffmpeg -r 8 -f image2 -s '.xsize*rescale.'x'.rescale*nifo*ysize.' -start_number 0 -i animation/bayeswave%%03d.png -vframes '.nmax.' -vcodec libx264 -crf 25 -pix_fmt yuv420p %s'\n",filename);
fprintf(framefile,"\n");
fprintf(framefile,"system(CMD)\n");
fprintf(framefile,"\n");
@@ -1648,11 +1647,9 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
{
fprintf(stdout,"\n");
fprintf(stdout,"************************ WARNING ************************\n");
- fprintf(stdout,"BayesLine requires the cleaning phase \n");
- fprintf(stdout,"Ignoring --noClean argument \n");
- fprintf(stdout,"If you really want to skip that step remove --bayesLine \n");
+ fprintf(stdout,"Assuming clean data outside of the analysis window \n");
+ fprintf(stdout,"For reliable results remove --noClean argument \n");
fprintf(stdout,"*********************************************************\n");
- data->cleanModelFlag = 1;
}
//if doing a short run disable spline thermodynamic integration
@@ -1662,7 +1659,7 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
fprintf(stdout,"************************ WARNING ************************\n");
fprintf(stdout,"Too few iterations for spline integration (%i)\n",chain->count);
fprintf(stdout,"Reverting to trapezoid integration\n");
- fprintf(stdout,"Requst --Niter > 1000000 for spline integration\n");
+ fprintf(stdout,"Request --Niter > 1000000 for spline integration\n");
fprintf(stdout,"*********************************************************\n");
data->splineEvidenceFlag = 0;
}
diff --git a/src/BayesWaveMCMC.c b/src/BayesWaveMCMC.c
index e30acd48..86bcacfd 100644
--- a/src/BayesWaveMCMC.c
+++ b/src/BayesWaveMCMC.c
@@ -892,7 +892,6 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
int N = data->N;
int NI = data->NI;
int M = chain->count;
- if(data->cleanFlag) M/=5;
int NC = chain->NC;
int i, ic, ifo, n;
@@ -1146,22 +1145,6 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
recompute_residual(data, model, chain);
if(burnFlag==0 && chain->mc%1000==0 && data->signalFlag) TFprop_signal(data, prior->range, tf, model[chain->index[0]]->projection);
- //During cleaning phase track max and min BL PSDs to set up priors for model selection run
- if(data->bayesLineFlag && data->cleanFlag )
- {
- ic = chain->index[0];
- for(ifo=0; ifo<NI; ifo++)
- {
- for(i=0; i< bayesline[ic][ifo]->data->ncut; i++)
- {
- if(bayesline[ic][ifo]->Snf[i] > bayesline[ic][ifo]->priors->upper[i])
- bayesline[ic][ifo]->priors->upper[i] = bayesline[ic][ifo]->Snf[i]*10.;
- else if(bayesline[ic][ifo]->Snf[i] < bayesline[ic][ifo]->priors->lower[i])
- bayesline[ic][ifo]->priors->lower[i] = bayesline[ic][ifo]->Snf[i]/10.;
- }
- }
- }
-
//Parallel tempering
if(NC>1 && chain->mc>1)
{
@@ -1456,7 +1439,6 @@ void EvolveBayesLineParameters(struct Data *data, struct Model **model, struct B
struct BayesLineParams **bl_x = bayesline[chain->index[ic]];
int priorFlag = 1;
- if(data->cleanFlag) priorFlag = 0;
//update PSD for each interferometer
for(ifo=0; ifo<NI; ifo++)
--
GitLab