Commit 947c598c authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

small memory fix in bootstrap. could that be it?

Anderson-Darling in BWP


git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@632 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent 030346f8
......@@ -1127,7 +1127,7 @@ void bootcorrposdef(double **data, int NC, int N, double **cmat, double *mu, lon
free_dmatrix(rmat,0,NC,0,NC);
free_dmatrix(rvar,0,NC,0,NC);
free_ivector(cyc,0,N);
free_ivector(cyc,0,N/2);
free_ivector(cnts,0,NC);
free_dmatrix(fmean,0,NC,0,bsteps);
free_dvector(var,0,NC);
......
......@@ -631,11 +631,7 @@ static void resume_sampler(struct Data *data, struct Chain *chain, struct Model
/* CLEAR OUT CHECKPOINT DIRECTORY */
/* */
/****************************************************/
char command[128];
sprintf(command,"rm -f checkpoint/*");
fprintf(stdout,"...resuming model %s at iteration %i\n",modelname,chain->mc);
fprintf(stdout,"...clearing out checkpoint directory\n");
system(command);
fflush(stdout);
for(ifo=0; ifo<NI; ifo++)
......@@ -661,7 +657,10 @@ static void save_sampler(struct Data *data, struct Chain *chain, struct Model **
FILE *fptr = NULL;
printf("...checkpoint model %s at iteration %i\n",modelname,chain->mc);
char command[128];
sprintf(command,"rm -f checkpoint/*");
fprintf(stdout,"...clearing out checkpoint directory...\n");
system(command);
/****************************************************/
......@@ -831,6 +830,7 @@ static void save_sampler(struct Data *data, struct Chain *chain, struct Model **
fprintf(fptr,"%.12g %lg\n",data->logZfull,data->varZfull);
fclose(fptr);
fprintf(stdout,"...saved model %s at iteration %i\n",modelname,chain->mc);
}
/* ********************************************************************************** */
......
......@@ -10,6 +10,8 @@
#include <string.h>
#include <math.h>
#include <sys/stat.h>
#include <gsl/gsl_cdf.h>
#include <lal/LALInferenceTemplate.h>
#include <lal/LALSimIMR.h>
......@@ -102,6 +104,8 @@ void median_and_CIs_fourierdom(double **waveforminfo, int N, int Nwaveforms, dou
void Q_scan(struct Data *data, double *median_waveform, double Q, double deltaT, double deltaF, double df_resolution, double *invPSD, FILE *outfile);
void Transform(double *hoff, double Q, double dt, double df_resolution, double fmin, double *invPSD, int n, int m, FILE *outfile);
void anderson_darling_test(char *outdir, int ifo, double fmin, double Tobs, char *model);
/* ============================ MAIN PROGRAM ============================ */
......@@ -1420,6 +1424,11 @@ int main(int argc, char *argv[])
print_stats(data, GPStrig, injectionFlag, outdir, "signal", dimension);
for(ifo=0; ifo<NI; ifo++)
{
anderson_darling_test(outdir, ifo, data->fmin, data->Tobs, "signal");
}
} //End post processing of signal model
......@@ -1740,6 +1749,12 @@ int main(int argc, char *argv[])
print_stats(data, GPStrig, injectionFlag, outdir, "glitch", dimension);
for(ifo=0; ifo<NI; ifo++)
{
anderson_darling_test(outdir, ifo, data->fmin, data->Tobs, "glitch");
}
} //End post processing of glitch model
......@@ -1885,6 +1900,12 @@ int main(int argc, char *argv[])
print_stats(data, GPStrig, injectionFlag, outdir, "clean", dimension);
for(ifo=0; ifo<NI; ifo++)
{
anderson_darling_test(outdir, ifo, data->fmin, data->Tobs, "clean");
}
} //End post processing of cleaning phase
......@@ -1977,9 +1998,12 @@ int main(int argc, char *argv[])
free_dmatrix(noise_median, 0,NI-1,0,N/2-1);
free_d3tensor(noise_model, 0, NI-1, 0, N_sample_noise_model-1, 0, N/2-1);
for(ifo=0; ifo<NI; ifo++)
{
anderson_darling_test(outdir, ifo, data->fmin, data->Tobs, "noise");
}
}//End post processing of cleaning phase
}//End post processing of noise phase
......@@ -4111,3 +4135,86 @@ int acor(double *mean, double *sigma, double *tau, double X[], int L, int maxlag
free(C);
return 0;
}
void anderson_darling_test(char *outdir, int ifo, double fmin, double Tobs, char *model)
{
char filename[1024];
FILE *dFile;
FILE *nFile;
int N,n,nn;
double f,dre,dim,Sn,Sn25,Sn75,Sn05,Sn95;
int Nsamp;
double *samples=NULL;
// how many samples in the data?
N=0;
Nsamp=0;
//get residual file
sprintf(filename,"%s/fourier_domain_%s_residual_ifo%i.dat",outdir,model,ifo);
dFile = fopen(filename,"r");
while(!feof(dFile))
{
fscanf(dFile,"%lg %lg %lg",&f,&dre,&dim);
N++;
if(f>fmin) Nsamp++;
}
N--;
fclose(dFile);
// get memory to ingest file
Nsamp *= 2; //re and im parts
samples=malloc(Nsamp*sizeof(double));
nn=0;
// combine Fourier coefficints
//get residual file
sprintf(filename,"%s/fourier_domain_%s_residual_ifo%i.dat",outdir,model,ifo);
dFile = fopen(filename,"r");
//get PSD file
sprintf(filename,"%s/%s_median_PSD.dat.%i",outdir,model,ifo);
nFile = fopen(filename,"r");
//whiten residual
for(n=0; n<N; n++)
{
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)
{
samples[2*nn] = dre/sqrt(Sn/(Tobs/2.));
samples[2*nn+1] = dim/sqrt(Sn/(Tobs/2.));
nn++;
}
}
fclose(dFile);
fclose(nFile);
gsl_sort(samples,1,Nsamp);
double S,A;
// Anderson-Darling statistics
S=0.0;
for(n=0; n<Nsamp; n++)
{
S += ( log(gsl_cdf_ugaussian_P(samples[n])) + log(1 - gsl_cdf_ugaussian_P(samples[Nsamp-1-n])) )*(2.*(double)(n+1)-1.0)/(double)(Nsamp);
}
A = (double)(-1*Nsamp) - S;
free(samples);
sprintf(filename,"%s/%s_anderson_darling.dat.%i",outdir,model,ifo);
FILE *afile = fopen(filename,"w");
fprintf(afile,"%lg\n",A);
fclose(afile);
}
......@@ -84,6 +84,13 @@ void draw_glitch_amplitude(double *params, double *Snf, double Sa, long *seed, d
alpha = ran2(seed);
k++;
//you had your chance
if(k>10000)
{
SNR=0.0;
break;
}
}
//SNR defined with Sn(f) but Snf array holdes <n_i^2>
......@@ -143,6 +150,14 @@ void draw_signal_amplitude(double *params, double *Snf, double Sa, long *seed, d
k++;
//you had your chance
if(k>10000)
{
SNR=0.0;
break;
}
}
//SNR defined with Sn(f) but Snf array holdes <n_i^2>
......
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