Skip to content
GitLab
Commit 3e85e9d3 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

Merge branch 'odds-n-ends' into 'master' 

A lot:

See merge request !55
parents f0a48c6a bca1d22d
Hide whitespace changes
Inline Side-by-side
src/BayesLine.c
View file @ 3e85e9d3
......@@ -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);
......
src/BayesLine.h
View file @ 3e85e9d3
......@@ -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;
......
src/BayesWave.c
View file @ 3e85e9d3
......@@ -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;