Commit a953c25e authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

Merging O1_offline tag with trunk, cleaning up directories



git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@536 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent 59407e23
This diff is collapsed.
......@@ -165,15 +165,6 @@ def make_skyview(directory='.', mdc=None, NSIDE=64, ra=None, dec=None, results=N
# Shift the declination to match hp.projscatter conventions
decrad = np.pi*0.5 - decrad
if mdc is not None:
# Convert the right ascension to a reference angle from -pi to pi
# for mdc injections.
while injra > np.pi:
injra = injra - 2*np.pi
while injra < -np.pi:
injra = injra + 2*np.pi
# -- Plot the skymap and injection location
skymap = kde.as_healpix(NSIDE, nest=True)
......
......@@ -1049,11 +1049,18 @@ void BurstRJMCMC(struct Data *data, struct Model **model, struct Chain *chain, s
chain->dT[0] = pow(chain->temperature[chain->NC-1],(1./(double)(NC-1)));
for(ic=0; ic<chain->NC-1; ic++)
{
//chain->temperature[ic]=chain->temperature[ic-1]*chain->dT[0];
chain->temperature[ic]=exp((double)(ic*ic)/(double)(((chain->NC-1)*(chain->NC-1))/(log(1.0e6)))); //empirically determined initial guess for temp spacing
if(chain->NC==15)
{
//chain->temperature[ic]=chain->temperature[ic-1]*chain->dT[0];
chain->temperature[ic]=exp((double)(ic*ic)/(double)(((chain->NC-1)*(chain->NC-1))/(log(1.0e6)))); //empirically determined initial guess for temp spacing
chain->temperature[ic] = 1./(1.0 - 0.07*(double)ic);
if(ic>14)chain->temperature[ic]=0.5*(1.0e6-chain->temperature[ic-1]);
chain->temperature[ic] = 1./(1.0 - 0.07*(double)ic);
}
else
{
chain->temperature[ic] = pow(chain->dT[0],ic);
}
printf("ic=%i temp=%g\n",ic,chain->temperature[ic]);
}
}
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -416,6 +416,8 @@ void nrerror(char error_text[]);
/* ********************************************************************************** */
void Shf_Geocenter(struct Data *data, struct Network *projection, double *SnGeo, double *params);
double symmetric_snr_ratio(struct Data *data, struct Network *projection, double *params);
/* ********************************************************************************** */
/* */
......@@ -444,6 +446,7 @@ void getChannels(ProcessParamsTable *commandLine, char **channel);
int outputFrameFile(ProcessParamsTable *commandLine, struct Data *data, struct Model *model);
void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf, double eta, double Tobs, double fix, int imin, int imax, double tmin, double tmax);
void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, double eta, double Tobs, double fix, int imin, int imax, double tmin, double tmax, double fmax);
void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *prior, ProcessParamsTable *commandLine);
void parse_glitch_parameters(struct Data *data, struct Model *model, FILE **paramsfile, double **grec);
......
......@@ -4353,7 +4353,7 @@ double TFprop(double f, double t, double **prop, int tsize, int fsize, double To
return(alpha);
}
/*
double glitch_amplitude_prior(double *params, double *Snf, double Sa, double Tobs, double SNRpeak)
{
int i;
......@@ -4374,6 +4374,34 @@ double glitch_amplitude_prior(double *params, double *Snf, double Sa, double Tob
return(alpha);
}
*/
double glitch_amplitude_prior(double *params, double *Snf, double Sa, double Tobs, double SNRpeak)
{
int i;
double SNR, alpha;
double dfac, dfac3;
Sa = 1.0;
// x/a^2 exp(-x/a) prior on SNR. Peaks at x = a. Good choice is a=5
i = (int)(params[1]*Tobs);
//SNR defined with Sn(f) but Snf array holdes <n_i^2>
SNR = analytic_snr(params, Snf, Tobs);
dfac = 1.+SNR/(2.*SNRpeak);
dfac3 = dfac*dfac*dfac;
alpha = log((SNR)/(2.*SNRpeak*SNRpeak*dfac3)) + log(SNR/params[3]);
return(alpha);
}
double signal_amplitude_prior(double *params, double *Snf, double Sa, double Tobs, double SNRpeak)
{
......@@ -4395,13 +4423,13 @@ double signal_amplitude_prior(double *params, double *Snf, double Sa, double Tob
dfac5 = dfac*dfac*dfac*dfac*dfac;
alpha = log((3.*SNR)/(4.*SNRpeak*SNRpeak*dfac5)) + log(SNR/params[3]);
alpha = log((3.*SNR)/(4.*SNRpeak*SNRpeak*dfac5)) + log(SNR/params[3]);
return(alpha);
}
/*
void draw_glitch_amplitude(double *params, double *Snf, double Sa, long *seed, double Tobs, double **range, double SNRpeak)
{
int i, k;
......@@ -4449,6 +4477,65 @@ void draw_glitch_amplitude(double *params, double *Snf, double Sa, long *seed, d
//SNR defined with Sn(f) but Snf array holdes <n_i^2>
params[3] = SNR/sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
}
*/
void draw_glitch_amplitude(double *params, double *Snf, double Sa, long *seed, double Tobs, double **range, double SNRpeak)
{
int i, k;
double SNR, den=-1.0, alpha=1.0;
double max;
double invmax;
double dfac, dfac3;
Sa = 1.0;
double SNR2 = 2.0*SNRpeak;
double SNRsq = 2.0*SNRpeak*SNRpeak;
dfac = 1.+SNRpeak/(SNR2);
dfac3 = dfac*dfac*dfac;
max = (SNRpeak)/(2.*SNRsq*dfac3);
invmax = 1./max;
i = (int)(params[1]*Tobs);
double SNRmin = range[3][0]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
double SNRmax = range[3][1]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
k = 0;
SNR = SNRmin + (SNRmax-SNRmin)*ran2(seed);
dfac = 1.+SNR/(SNR2);
dfac3 = dfac*dfac*dfac;
den = (SNR)/(SNRsq*dfac3);
den *= invmax;
alpha = ran2(seed);
while(alpha > den)
{
SNR = SNRmin + (SNRmax-SNRmin)*ran2(seed);
dfac = 1.+SNR/(SNR2);
dfac3 = dfac*dfac*dfac;
den = (SNR)/(2.*SNRsq*dfac3);
den *= invmax;
alpha = ran2(seed);
k++;
}
//SNR defined with Sn(f) but Snf array holdes <n_i^2>
params[3] = SNR/sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
}
void draw_signal_amplitude(double *params, double *Snf, double Sa, long *seed, double Tobs, double **range, double SNRpeak)
......@@ -6088,6 +6175,26 @@ void Shf_Geocenter(struct Data *data, struct Network *projection, double *SnGeo,
}
double symmetric_snr_ratio(struct Data *data, struct Network *projection, double *params)
{
int ifo;
int NI = data->NI;
double AntennaPattern;
double ecc = params[3];
double SNR[NI];
for(ifo=0; ifo<NI; ifo++)
{
AntennaPattern = projection->Fplus[ifo]*projection->Fplus[ifo] + ecc*ecc*projection->Fcross[ifo]*projection->Fcross[ifo];
SNR[ifo] = AntennaPattern;
}
return SNR[0]*SNR[1]/(SNR[0]+SNR[1])/(SNR[0]+SNR[1]);
}
/* ********************************************************************************** */
/* */
/* Fourier Routines */
......@@ -6769,6 +6876,64 @@ void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf,
fclose(waveout);
}
void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, double eta, double Tobs, double fix, int imin, int imax, double tmin, double tmax, double fmax)
{
/*
TODO: invFFT comes out time-reversed.
LAL uses exp(-i2pift) for their Fourier transforms
Why doesn't drealft(ty-1,N,1) work for inverse FFT?
In the meantime, hacked invFFT by changing sign of imaginary part
*/
int i;
double x,t;
double xx, yy;
double deltaF = 1./Tobs;
FILE *waveout = fopen(filename,"w");
// hdot = IFT(-2pi*i*f*h(f))
double *ht = malloc(N*sizeof(double));
for(i=0; i<N; i++)ht[i]= -TPI*h[i]; // Take care of the -2pi part
ht[0] = 0.0;
ht[1] = 0.0;
for(i=1; i< N/2; i++)
{
if(i>imin && i<imax)
{
x = sqrt(Snf[i]*eta);
// Need to switch real and imaginary parts [ie -i*(x+iy)=-ix+y]
xx = ht[2*i]; // Real part
yy = ht[2*i+1]; // Imaginary part
ht[2*i] = yy*(double)i*deltaF/x; // f = i*deltaf
ht[2*i+1] = xx*(double)i*deltaF/x; // Change sign of new imaginary part so it doesn't come out time reversed
}
else
{
ht[2*i]=ht[2*i+1]=0.0;
}
}
drealft(ht-1,N,-1);
double t0 = Tobs-2.0;
for(i=0; i<N; i++)
{
t = (double)(i)/(double)(N)*Tobs;
if(t>=tmin && t<tmax)fprintf(waveout,"%e %e\n", t-t0, ht[i]/fix);
}
free(ht);
fclose(waveout);
}
void print_frequency_domain_waveforms(char filename[], double *h, int N, double *Snf, double Tobs, double eta, int imin, int imax)
{
/*
......@@ -7270,7 +7435,7 @@ void initialize_priors(struct Data *data, struct Prior *prior, int omax)
if(!strcmp(data->channels[0],"LDAS-STRAIN"))
{
Amin = 5.e-25;
Amax = 1.e-19;
Amax = 1.e-18;
}
//VSR2 strain
......@@ -7291,7 +7456,7 @@ void initialize_priors(struct Data *data, struct Prior *prior, int omax)
else
{
Amin = 1.e-26;
Amax = 1.e-19;
Amax = 1.e-18;
}
//double dt = data->Tobs/(double)data->tsize;
......@@ -7427,7 +7592,7 @@ void reset_priors(struct Data *data, struct Prior *prior)
if(!strcmp(data->channels[0],"LDAS-STRAIN"))
{
Amin = 5.e-25;
Amax = 1.e-19;
Amax = 1.e-18;
}
//VSR2 strain
......@@ -7448,7 +7613,7 @@ void reset_priors(struct Data *data, struct Prior *prior)
else
{
Amin = 1.e-26;
Amax = 1.e-19;
Amax = 1.e-18;
}
double Qmin = data->Qmin;
......
This diff is collapsed.
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