diff --git a/src/BayesLine.c b/src/BayesLine.c
index f4783284b9d6c7d7c09925704e09f9fa2a7c8b15..c2133e3432af5bceee6b023c62104d1b7a94dd81 100644
--- a/src/BayesLine.c
+++ b/src/BayesLine.c
@@ -113,6 +113,7 @@ static double logprior(double *lower, double *upper, double *Snf, int ilow, int
if(Snf[i]>upper[i] || Snf[i]<lower[i])
{
lgp=-1e60;
+ //printf("The PSD is outside the prior. Offending point: %i %lg %lg %lg\n",i,Snf[i],lower[i],upper[i]);
}
}
return (lgp);
@@ -706,9 +707,20 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
if(logPsx<0)
{
printf("how did PSD come in outside of the prior?\n");
+ printf("%lg\n",logPsx);
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]);
+ for(i=0; i<ncut; i++)
+ {
+ if(Snx[i]>priors->upper[i] || Snx[i]<priors->lower[i])
+ {
+
+ printf("The PSD is outside the prior. Offending point: %i %lg %lg %lg\n",i,Snx[i],priors->lower[i],priors->upper[i]);
+ }
+ }
+
+
print_line_model(stdout, bayesline);
exit(1);
@@ -1657,7 +1669,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_x->n; j++) fprintf(fptr,"%e %e %e ", bayesline->lines_x->f[j], bayesline->lines_x->A[j], bayesline->lines_x->Q[j]);
+ for(j=0; j< bayesline->lines_x->n; j++) fprintf(fptr,"%lg %lg %lg ", bayesline->lines_x->f[j], bayesline->lines_x->A[j], bayesline->lines_x->Q[j]);
fprintf(fptr,"\n");
}
@@ -2646,9 +2658,14 @@ void blstart(double *data, double *residual, int N, double dt, double fmin, int
k = (int)(fmin*Tobs)-1;
sprintf(filename, "%s_burnin_psd.dat", ifo);
outFile = fopen(filename,"w");
- for (i = 0; i < Nint; ++i)
+ for (i = 0; i < Nint+k; ++i)
{
- fprintf(outFile,"%.15e %.15e\n", (double)(i+k)/Tobs, (Sn[k+i]+exp(yint[i]))/(Tobs/2.));
+ if (i<k){
+ fprintf(outFile,"%.15e %.15e\n", (double)(i)/Tobs, 1.);
+ }
+ else{
+ fprintf(outFile,"%.15e %.15e\n", (double)(i)/Tobs, (Sn[i]+exp(yint[i-k]))/(Tobs/2.));
+ }
}
fclose(outFile);
diff --git a/src/BayesWavePost.c b/src/BayesWavePost.c
index 076ce8cfb5f3d4d5eab14da2a175351661ef7c1a..e0a6377b4abd06cb0e6b8ad675cc2c67f3616a4d 100644
--- a/src/BayesWavePost.c
+++ b/src/BayesWavePost.c
@@ -3375,7 +3375,7 @@ void anderson_darling_test(char *outdir, int ifo, double fmin, double Tobs, char
sprintf(filename,"%s/fourier_domain_%s_median_residual_%s.dat",outdir,model,ifos[ifo]);
dFile = fopen(filename,"r");
- //get PSD file
+ //get median PSD file
sprintf(filename,"%s/%s_median_PSD_%s.dat",outdir,model,ifos[ifo]);
nFile = fopen(filename,"r");
@@ -3536,13 +3536,79 @@ void anderson_darling_test(char *outdir, int ifo, double fmin, double Tobs, char
// p-value of A^2
double p = 1.-AD(Nsamp,A);
-
+
free(samples);
sprintf(filename,"%s/%s_anderson_darling_%s.dat",outdir,model,ifos[ifo]);
FILE *afile = fopen(filename,"w");
fprintf(afile,"%lg %lg\n",A, p);
fclose(afile);
+ // Besides the median PSD, we also compute the AD and p value with the burnin PSD
+
+ double *samplesB=NULL;
+ samplesB=malloc(Nsamp*sizeof(double));
+ nn=0;
+ //get residual file
+ sprintf(filename,"%s/fourier_domain_%s_median_residual_%s.dat",outdir,model,ifos[ifo]);
+ dFile = fopen(filename,"r");
+
+ //get burnin PSD file
+ sprintf(filename,"%s_burnin_psd.dat",ifos[ifo]);
+ nFile = fopen(filename,"r");
+
+ //whiten residual with the burnin PSD
+ for(n=0; n<N; n++)
+ {
+ fscanf(dFile,"%lg %lg %lg",&f,&dre,&dim);
+ fscanf(nFile,"%lg %lg",&f,&Sn);
+ if(f>=fmin)
+ {
+ //samples[2*nn] = dre/sqrt(Sn/(Tobs/2.));
+ //samples[2*nn+1] = dim/sqrt(Sn/(Tobs/2.));
+ samplesB[2*nn] = dre/sqrt(Sn*(Tobs/2.)/2.);//I am totally winging this factor here
+ samplesB[2*nn+1] = dim/sqrt(Sn*(Tobs/2.)/2.);
+ nn++;
+ }
+ }
+ fclose(dFile);
+ fclose(nFile);
+
+ //Compute and print the Anderson-Darling statistic for the residual whitened by the burnin PSD
+ S=0.0;
+
+ gsl_sort(samplesB,1,Nsamp);
+
+ for(n=0; n<Nsamp; n++)
+ {
+ double x = gsl_cdf_ugaussian_P(samplesB[n]);
+ if(x<0.999)
+ {
+ lx = log(x);
+ ly = log(1.-x);
+ }
+ else
+ {
+ double u = samplesB[n];
+ double z = 1.0-1.0/(u*u)+3.0/(u*u*u*u)-15.0/(u*u*u*u*u*u)+105.0/(u*u*u*u*u*u*u*u);
+ double y = z*exp(-u*u/2.0)/(u*sqrt(M_PI));
+ lx = -y;
+ ly = -u*u/2.0 +log(z/(u*sqrt(M_PI)));
+ }
+ S += ((2.*(double)(n+1)-1.0)*lx + ((double)(2*Nsamp)-2.*(double)(n+1)+1.0)*ly)/(double)(Nsamp);
+ }
+ A = (double)(-1*Nsamp) - S;
+
+ // p-value of A^2
+ p = 1.-AD(Nsamp,A);
+
+ free(samplesB);
+ sprintf(filename,"%s/%s_anderson_darling_burnin_%s.dat",outdir,model,ifos[ifo]);
+ FILE *afileB = fopen(filename,"w");
+ fprintf(afileB,"%lg %lg\n",A, p);
+ fclose(afileB);
+
+
+
}
void get_time_domain_envelope(double *envelope, double *hdata_sample, int N, int imin, int imax)