Skip to content

GitLab

bayeswave
bayeswave
Commit 1c409a39 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files
Options

Merge branch 'psd-prior' into 'master' 

Psd prior

See merge request lscsoft/bayeswave!90
parents b8f7a62b 61770292
Hide whitespace changes
Inline Side-by-side
Showing with 74 additions and 107 deletions
src/BayesLine.c
View file @ 1c409a39
......@@ -525,7 +525,7 @@ void copy_lorentzianParams(lorentzianParams *origin, lorentzianParams *copy)
copy->size = origin->size;
int n;
for(n=0; n<origin->size; n++)
for(n=0; n<origin->n; n++)
{
copy->larray[n] = origin->larray[n];
......@@ -708,6 +708,8 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
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]);
print_line_model(stdout, bayesline);
exit(1);
}
......@@ -746,6 +748,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
for(mc=0; mc < steps; mc++)
{
typ=-1;
check = 0;
//copy over current state
lines_y->n = lines_x->n;
......@@ -883,8 +886,6 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
}
check = 0;
if(nsy < 5 || nsy > nspline-1) check = 1;
for(i=0; i<nsy; i++)
......@@ -912,10 +913,8 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
typ = 3;
}
check = 0;
if(lines_y->n < 0 || lines_y->n > tmax) check = 1;
if(check == 0)
{
if(lines_y->n < lines_x->n)
......@@ -1145,8 +1144,6 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
logpx = logpy = 0.0;
}
check =0;
if(lines_y->A[ii] > priors->LAmax) check = 1;
if(lines_y->A[ii] < priors->LAmin) check = 1;
if(lines_y->f[ii] < flow) check = 1;
......@@ -1401,7 +1398,7 @@ void BayesLineFree(struct BayesLineParams *bptr)
}
void BayesLineSetup(struct BayesLineParams *bptr, double *freqData, double fmin, double fmax, double deltaT, double Tobs, int nstep)
void BayesLineSetup(struct BayesLineParams *bptr, double *freqData, double fmin, double fmax, double deltaT, double Tobs)
{
int i;
int n;
......@@ -1418,8 +1415,6 @@ void BayesLineSetup(struct BayesLineParams *bptr, double *freqData, double fmin,
bptr->TwoDeltaT = deltaT*2.0;
bptr->nstep = nstep;
/* set up GSL random number generator */
const gsl_rng_type *T = gsl_rng_default;
bptr->r = gsl_rng_alloc (T);
......@@ -1449,9 +1444,6 @@ void BayesLineSetup(struct BayesLineParams *bptr, double *freqData, double fmin,
// storage for data meta parameters
create_dataParams(bptr->data,bptr->freq,n);
// storage for line model for a segment. These get updated and passed back from the fitting routine
create_lorentzianParams(bptr->lines_x,bptr->data->tmax);
// storage for master line model
create_lorentzianParams(bptr->lines_x,bptr->data->tmax);
......@@ -1516,13 +1508,8 @@ void BayesLineInitialize(struct BayesLineParams *bayesline)
/* */
/******************************************************************************/
// maximum number of terms in the Lorentzian model
bayesline->data->tmax = 1000;
if(bayesline->data->tmax < 20) bayesline->data->tmax = 20;
if(bayesline->lines_x->n < 1 ) bayesline->lines_x->n = 1;
// Re-set dataParams structure to use full dataset
bayesline->data->flow = bayesline->data->fmin;
imax = (int)(floor(bayesline->data->fhigh*bayesline->data->Tobs));
......@@ -2395,7 +2382,7 @@ static void clean(double *D, double *Draw, double *sqf, double *freqs, double *S
}
void blstart(double *data, double *residual, int N, double dt, double fmin, int *Nsp, int *Nl, double *dspline, double *pspline, double *linef, double *lineh, double *lineQ)
void blstart(double *data, double *residual, int N, double dt, double fmin, int *Nsp, int *Nl, double *dspline, double *pspline, double *linef, double *lineh, double *lineQ, char *ifo)
{
int i, j, k, Nf, ii;
int Nlines;
......@@ -2416,11 +2403,9 @@ void blstart(double *data, double *residual, int N, double dt, double fmin, int
double t_rise, s1, s2, fac;
double *linew;
int modelprint;
char command[1024];
char filename[1024];
FILE *out;
FILE *outFile;
Q = Qs; // Q of transform
Tobs = (double)(N)*dt; // duration
......@@ -2441,7 +2426,6 @@ void blstart(double *data, double *residual, int N, double dt, double fmin, int
// Set the range of the spectrogram. fmin, fmax must be a power of 2
fmax = fny;
modelprint = 1;
D = (double*)malloc(sizeof(double)* (N));
Draw = (double*)malloc(sizeof(double)* (N));
......@@ -2527,17 +2511,15 @@ void blstart(double *data, double *residual, int N, double dt, double fmin, int
residual[2*i+1] = D[N-i] * sqrt(fac);
}
if(modelprint == 1)
{
sprintf(command, "PSD_%d_%d.dat", (int)(Q), (int)(Tobs));
out = fopen(command,"w");
for (i = 1; i < N/2; ++i)
{
fprintf(out,"%.15e %.15e %.15e %.15e\n", (double)(i)/Tobs, Sn[i]*fac, specD[i]*fac, sspecD[i]*fac);
}
fclose(out);
}
/*
sprintf(filename, "PSD_%d_%d.dat", (int)(Q), (int)(Tobs));
outFile = fopen(filename,"w");
for (i = 1; i < N/2; ++i)
{
fprintf(outFile,"%.15e %.15e %.15e %.15e\n", (double)(i)/Tobs, Sn[i]*fac, specD[i]*fac, sspecD[i]*fac);
}
fclose(outFile);
*/
// set up the spline model
int Nspline;
......@@ -2618,7 +2600,7 @@ void blstart(double *data, double *residual, int N, double dt, double fmin, int
xold = x;
}
Nlines = j+1;
Nlines = j;
// printf("There are %d lines\n", Nlines);
......@@ -2634,52 +2616,45 @@ void blstart(double *data, double *residual, int N, double dt, double fmin, int
pspline // reference grid of spline control points
dspline // values at reference points
*/
if(modelprint == 1)
for (i = 1; i < N/2; ++i) Sn[i] = 0.0;
for (i = imin; i < N/2; ++i)
{
for (i = 1; i < N/2; ++i) Sn[i] = 0.0;
for (i = imin; i < N/2; ++i)
f = (double)(i)/Tobs;
for (j = 0; j < Nlines; ++j)
{
f = (double)(i)/Tobs;
for (j = 0; j < Nlines; ++j)
x = 0.25*fabs(f - linef[j]);
if(x < linew[j])
{
x = 0.25*fabs(f - linef[j]);
if(x < linew[j])
{
Sn[i] += lineh[j]*pow(linef[j],4.0)/(pow(f*linef[j],2.0)+lineQ[j]*lineQ[j]*pow((pow(linef[j],2.0)-pow(f,2.0)),2.0));
}
Sn[i] += lineh[j]*pow(linef[j],4.0)/(pow(f*linef[j],2.0)+lineQ[j]*lineQ[j]*pow((pow(linef[j],2.0)-pow(f,2.0)),2.0));
}
}
int Nint;
double *xint, *yint;
Nint = (int)((fmax-fmin)*Tobs)+1;
xint = (double*)malloc(sizeof(double)*(Nint));
yint = (double*)malloc(sizeof(double)*(Nint));
for (i = 0; i < Nint; ++i)
{
xint[i] = fmin+(double)(i)/Tobs;
}
CubicSplineGSL(Nspline, pspline, dspline, Nint, xint, yint);
k = (int)(fmin*Tobs)-1;
sprintf(command, "PSD_check_%d_%d.dat", (int)(Q), (int)(Tobs));
out = fopen(command,"w");
for (i = 0; i < Nint; ++i)
{
fprintf(out,"%.15e %.15e %.15e\n", (double)(i+k)/Tobs, Sn[k+i]+exp(yint[i]), specD[k+i]*fac);
}
fclose(out);
free(xint);
free(yint);
}
int Nint;
double *xint, *yint;
Nint = (int)((fmax-fmin)*Tobs)+1;
xint = (double*)malloc(sizeof(double)*(Nint));
yint = (double*)malloc(sizeof(double)*(Nint));
for (i = 0; i < Nint; ++i)
{
xint[i] = fmin+(double)(i)/Tobs;
}
CubicSplineGSL(Nspline, pspline, dspline, Nint, xint, yint);
k = (int)(fmin*Tobs)-1;
sprintf(filename, "%s_burnin_psd.dat", ifo);
outFile = fopen(filename,"w");
for (i = 0; i < Nint; ++i)
{
fprintf(outFile,"%.15e %.15e\n", (double)(i+k)/Tobs, (Sn[k+i]+exp(yint[i]))/(Tobs/2.));
}
fclose(outFile);
free(xint);
free(yint);
free(D);
free(Draw);
free(linew);
......@@ -2741,7 +2716,7 @@ static void Lorentzian(double *Slines, double Tobs, lorentzianParams *lines, int
}
void BayesLineBurnin(struct BayesLineParams *bayesline, double *timeData, double *freqData)
void BayesLineBurnin(struct BayesLineParams *bayesline, double *timeData, double *freqData, char *ifo)
{
//Initialize BayesLine data structures
BayesLineInitialize(bayesline);
......@@ -2769,7 +2744,7 @@ void BayesLineBurnin(struct BayesLineParams *bayesline, double *timeData, double
int N = 2*(int)floor(fmax*Tobs);
blstart(timeData, freqData, N, dt, data->fmin, &spline->n, &lines->n, spline->data, spline->points, lines->f, lines->A, lines->Q);
blstart(timeData, freqData, N, dt, data->fmin, &spline->n, &lines->n, spline->data, spline->points, lines->f, lines->A, lines->Q, ifo);
//Work space for assembling PSD model
double *f = malloc(sizeof(double)*(N+1));
......@@ -2806,6 +2781,7 @@ void BayesLineBurnin(struct BayesLineParams *bayesline, double *timeData, double
for(i=0; i<lines->n; i++)
{
lines->larray[i] = i;
//printf("Line %i: %g %g\n",i,lines->f[i], lines->A[i]);
Lorentzian(Sl, Tobs, lines, lines->larray[i], N);
}
......
src/BayesLine.h
View file @ 1c409a39
......@@ -113,14 +113,12 @@ struct BayesLineParams
int constantLogLFlag;
int nstep;
double TwoDeltaT;
gsl_rng *r;
};
void BayesLineFree(struct BayesLineParams *bptr);
void BayesLineSetup(struct BayesLineParams *bptr, double *freqData, double fmin, double fmax, double deltaT, double Tobs, int nstep);
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, int priorFlag);
void BayesLineInitialize(struct BayesLineParams *bayesline);
......@@ -160,5 +158,5 @@ void parse_line_model(FILE *fptr, struct BayesLineParams *bayesline);
void parse_spline_model(FILE *fptr, struct BayesLineParams *bayesline);
void BayesLineBurnin(struct BayesLineParams *bayesline, double *timeData, double *freqData);
void BayesLineBurnin(struct BayesLineParams *bayesline, double *timeData, double *freqData, char *ifo);
src/BayesWave.c
View file @ 1c409a39
......@@ -204,16 +204,16 @@ int main(int argc, char *argv[])
//read in frequency domain data from LALInference data structure
for(i=0; i<N/2; i++)
{
if(simDataFlag)
// if(simDataFlag)
// {
// dataPtr->freqData->data->data[imin] = dataPtr->freqData->data->data[imin+1];
// dataPtr->oneSidedNoisePowerSpectrum->data->data[imin] = dataPtr->oneSidedNoisePowerSpectrum->data->data[imin+1];
// }
if(i<imin+1)
{
dataPtr->freqData->data->data[imin] = dataPtr->freqData->data->data[imin+1];
dataPtr->oneSidedNoisePowerSpectrum->data->data[imin] = dataPtr->oneSidedNoisePowerSpectrum->data->data[imin+1];
}
if(i<imin)
{
freqData[ifo][2*i] = creal(dataPtr->freqData->data->data[imin]);
freqData[ifo][2*i+1] = cimag(dataPtr->freqData->data->data[imin]);
psd[ifo][i] = dataPtr->oneSidedNoisePowerSpectrum->data->data[imin]*Tobs/2.0;
freqData[ifo][2*i] = creal(dataPtr->freqData->data->data[imin+1]);
freqData[ifo][2*i+1] = cimag(dataPtr->freqData->data->data[imin+1]);
psd[ifo][i] = dataPtr->oneSidedNoisePowerSpectrum->data->data[imin+1]*Tobs/2.0;
}
else
{
......@@ -352,7 +352,7 @@ int main(int argc, char *argv[])
for(ic=0; ic<chain->NC; ic++)
{
bayesline[ic] = malloc(data->NI*sizeof(struct BayesLineParams *));
initialize_bayesline(bayesline[ic], data, psd, chain->blcount);
initialize_bayesline(bayesline[ic], data, psd);
}
/*
......@@ -396,7 +396,7 @@ int main(int argc, char *argv[])
for(i=0; i<N; i++) timeData[ifo][i]/=data->Tobs;
fprintf(stdout,"BayesLine search phase for IFO %s\n", data->ifos[ifo]);
BayesLineBurnin(bayesline[0][ifo], timeData[ifo], asd);
BayesLineBurnin(bayesline[0][ifo], timeData[ifo], asd, data->ifos[ifo]);
//Copy BayesLine model into BayesWave model
for(i=0; i<data->N/2; i++)
......@@ -1093,7 +1093,6 @@ void print_help_message(void)
fprintf(stdout," --Nchain number of parallel chains (20)\n");
fprintf(stdout," --Ncycle number of model updates per RJMCMC iteration (100)\n");
fprintf(stdout," --Nburnin number of burn-in iterations (50000)\n");
fprintf(stdout," --Nbayesline number of burn-in iterations for BayesLine (200000)\n");
fprintf(stdout," --maxLogL use maximized likelihood during burnin\n");
fprintf(stdout," --chainseed random number seed for Markov chain (1234)\n");
fprintf(stdout," --runName run name for output files\n");
......
src/BayesWave.h
View file @ 1c409a39
......@@ -167,7 +167,6 @@ struct Chain
int cycle; // # of model updates per BurstRJMCMC iteration
int count; // # of BurstRJMCMC iterations
int blcount;// # of BayesLine burnin iterations
int zcount; // # of logL samples for thermodynamic integration
int mcount; // # of MCMC attempts
int scount; // # of RJMCMC attempts
......
src/BayesWaveIO.c
View file @ 1c409a39
......@@ -800,7 +800,6 @@ void print_run_stats(FILE *fptr, struct Data *data, struct Chain *chain)
fprintf(fptr, " number of samples %i\n", data->N);
fprintf(fptr, " number of iterations %i\n", chain->count);
fprintf(fptr, " number of burn-in iterations %i\n", chain->burn);
fprintf(fptr, " number of BL burn-in iterations %i\n", chain->blcount);
fprintf(fptr, " max N wavelets %i\n", data->Dmax);
fprintf(fptr, " max Q for wavelets %.0f\n", data->Qmax);
fprintf(fptr, " number of chains %i\n", chain->NC);
......@@ -1361,10 +1360,6 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
ppt = LALInferenceGetProcParamVal(commandLine, "--Nburnin");
if(ppt) chain->burn = atoi(ppt->value);
else chain->burn = 50000;
ppt = LALInferenceGetProcParamVal(commandLine, "--Nbayesline");
if(ppt) chain->blcount = atoi(ppt->value);
else chain->blcount = 200000;
ppt = LALInferenceGetProcParamVal(commandLine, "--maxLogL");
if(ppt) data->searchFlag = 1;
......@@ -1896,11 +1891,11 @@ void export_cleaned_data(struct Data *data, struct Model *model)
for(ifo=0; ifo<NI; ifo++)
{
sprintf(filename,"%s%s_res.dat",data->runName,data->ifos[ifo]);
sprintf(filename,"%s%s_fairdraw_res.dat",data->runName,data->ifos[ifo]);
resFile = fopen(filename,"w");
sprintf(filename,"%s%s_psd.dat",data->runName,data->ifos[ifo]);
sprintf(filename,"%s%s_fairdraw_psd.dat",data->runName,data->ifos[ifo]);
psdFile = fopen(filename,"w");
sprintf(filename,"%s%s_asd.dat",data->runName,data->ifos[ifo]);
sprintf(filename,"%s%s_fairdraw_asd.dat",data->runName,data->ifos[ifo]);
asdFile = fopen(filename,"w");
imin = data->imin;//(int)(floor(data->bayesline[0]->data->flow * data->bayesline[0]->data->Tobs));
......
src/BayesWaveModel.c
View file @ 1c409a39
......@@ -342,7 +342,7 @@ void resize_model(struct Data *data, struct Chain *chain, struct Prior *prior, s
struct BayesLineParams **bptr;
bptr=malloc(data->NI*sizeof(struct BayesLineParams *));
initialize_bayesline(bptr,data,psd,chain->blcount);
initialize_bayesline(bptr,data,psd);
for(ifo=0; ifo<data->NI; ifo++) copy_bayesline_params(bayesline[0][ifo], bptr[ifo]);
......@@ -371,7 +371,7 @@ void resize_model(struct Data *data, struct Chain *chain, struct Prior *prior, s
initialize_model(model[ic], data, prior, psd, chain->seed);
bayesline[ic] = malloc(data->NI*sizeof(struct BayesLineParams *));
initialize_bayesline(bayesline[ic], data, psd, chain->blcount);
initialize_bayesline(bayesline[ic], data, psd);
for(ifo=0; ifo<data->NI; ifo++) copy_bayesline_params(bptr[ifo], bayesline[ic][ifo]);
}
......@@ -889,7 +889,7 @@ void free_bayesline(struct BayesLineParams **bayesline, struct Data *data)
}
void initialize_bayesline(struct BayesLineParams **bayesline, struct Data *data, double **psd, int nstep)
void initialize_bayesline(struct BayesLineParams **bayesline, struct Data *data, double **psd)
{
int i,ifo;
......@@ -924,7 +924,7 @@ void initialize_bayesline(struct BayesLineParams **bayesline, struct Data *data,
}
//Allocates arrays and sets constants for for BayesLine
BayesLineSetup(bayesline[ifo], data->s[ifo], data->fmin, data->fmax, data->dt, data->Tobs, nstep);
BayesLineSetup(bayesline[ifo], data->s[ifo], data->fmin, data->fmax, data->dt, data->Tobs);
}
// FILE *fptr=fopen("psdprior.dat","w");
......
src/BayesWaveModel.h
View file @ 1c409a39
......@@ -64,7 +64,7 @@ void copy_ext_model (struct Model *origin, struct Model *copy, int N, int NI)
void copy_psd_model (struct Model *origin, struct Model *copy, int N, int NI);
void free_bayesline(struct BayesLineParams **bayesline, struct Data *data);
void initialize_bayesline(struct BayesLineParams **bayesline, struct Data *data, double **psd, int nstep);
void initialize_bayesline(struct BayesLineParams **bayesline, struct Data *data, double **psd);
void reset_likelihood(struct Data *data);
void reset_params(struct Wavelet *wave_x, struct Wavelet *wave_y, double **range);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment