Merge branch 'ADoff' into 'master'

Add computation of the AD for different freq ranges

See merge request !181

src/BayesWavePost.c

@@ -3339,7 +3339,7 @@ void anderson_darling_test(char *outdir, int ifo, double fmin, double Tobs, char
   FILE *pFile;
   FILE *cFile;
 
-  int N,n,nn;
+  int N,n,nn, i;
   double f,dre,dim,Sn,Sn25,Sn75,Sn05,Sn95;
 
   int Nsamp;
@@ -3508,38 +3508,85 @@ void anderson_darling_test(char *outdir, int ifo, double fmin, double Tobs, char
   double S,A;
 
   // Anderson-Darling statistic
-  S=0.0;
   double lx,ly;
+
+  //We want to compute AD for different frequency ranges, from srate 256 up to the actual srate of the run.
   
-  for(n=0; n<Nsamp; n++)
+  double factorAD = ((Nsamp+2.*fmin*Tobs)/256./Tobs);// srate/256
+  int iAD = (int)(log(factorAD)/log(2.));
+  double p, c;
+  int newNsamp;
+
+  //We also need the unsorted samples again 
+
+  sprintf(filename,"%s/fourier_domain_%s_median_residual_%s.dat",outdir,model,ifos[ifo]);
+  dFile = fopen(filename,"r");
+  sprintf(filename,"%s/%s_median_PSD_%s.dat",outdir,model,ifos[ifo]);
+  nFile = fopen(filename,"r");  
+  nn=0;
+  for(n=0; n<N; n++)
   {
-    double x = gsl_cdf_ugaussian_P(samples[n]);
-    if(x<0.999)
-    {
-      lx = log(x);
-      ly = log(1.-x);
-    }
-    else
+    fscanf(dFile,"%lg %lg %lg",&f,&dre,&dim);
+    fscanf(nFile,"%lg %lg %lg %lg %lg %lg",&f,&Sn,&Sn25,&Sn75,&Sn05,&Sn95);
+    if(f>=fmin)
     {
-      double u = samples[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)));
+      //samples[2*nn]   = dre/sqrt(Sn/(Tobs/2.));
+      //samples[2*nn+1] = dim/sqrt(Sn/(Tobs/2.));
+      samples[2*nn]   = dre/sqrt(Sn/2.);//the factor of 2 works, but we need to understand why better
+      samples[2*nn+1] = dim/sqrt(Sn/2.);
+      nn++;
     }
-    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
-  double p = 1.-AD(Nsamp,A);
-  
-  free(samples);
+  fclose(dFile);
+  fclose(nFile);
+
   sprintf(filename,"%s/%s_anderson_darling_%s.dat",outdir,model,ifos[ifo]);
   FILE *afile = fopen(filename,"w");
-  fprintf(afile,"%lg %lg\n",A, p);
+  c=1;
+
+  for (i=0; i<=iAD; i++){
+    newNsamp = (int)(Nsamp+2.*fmin*Tobs)/(c)-2.*fmin*Tobs;// Samples only up to srate = (Nsamp+2.*fmin*Tobs)/(2**i)/Tobs, namely 256, 512, ...
+    S=0.0;
+
+    double *newsamples=NULL;
+    newsamples=malloc(newNsamp*sizeof(double));
+
+    for(n=0; n<newNsamp; n++) newsamples[n]=samples[n];
+
+    gsl_sort(newsamples,1,newNsamp);
+
+    for(n=0; n<newNsamp; n++)
+    {
+      double x = gsl_cdf_ugaussian_P(newsamples[n]);
+      if(x<0.999)
+      {
+        lx = log(x);
+        ly = log(1.-x);
+      }
+      else
+      {
+        double u = newsamples[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*newNsamp)-2.*(double)(n+1)+1.0)*ly)/(double)(newNsamp);
+    }
+    A = (double)(-1*newNsamp) - S;
+    
+    // p-value of A^2
+    p = 1.-AD(newNsamp,A);
+    fprintf(afile,"%g %lg %lg\n",(Nsamp+2.*fmin*Tobs)/(c)/Tobs, A, p);
+    c *=2;
+  
+    free(newsamples);
+  }
+
+  free(samples); 
   fclose(afile);
 
+
   // Besides the median PSD, we also compute the AD and p value with the burnin PSD
 
   double *samplesB=NULL;
@@ -3575,6 +3622,7 @@ void anderson_darling_test(char *outdir, int ifo, double fmin, double Tobs, char
 
   gsl_sort(samplesB,1,Nsamp);
 
+
   for(n=0; n<Nsamp; n++)
     {
        double x = gsl_cdf_ugaussian_P(samplesB[n]);