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

Two O2-related changes:

+less memory for TF proposal
+trying to generalize/tune BayesLine settings to supply PSDs for BNS triggers



git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@662 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent effd9f23
......@@ -1158,7 +1158,7 @@ void create_dataParams(dataParams *data, double *f, int n)
// size of segments in Hz
// If frequency snippets are too large need longer initial runs to get convergence
data->fstep = 256./data->Tobs;//30;//FSTEP;//9.0;//30.0;
data->fstep = data->fmax/16;//256./data->Tobs;//30;//FSTEP;//9.0;//30.0;
// This sets the maximum number of Lorentzian lines per segment.
// For segements shorter than 16 Hz this always seems to be enough
......@@ -1392,7 +1392,7 @@ void BayesLineNonMarkovianFit(struct BayesLineParams *bayesline, int *nj)
lines_full->A[kk] = lines_x->A[j];
kk++;
}
printf(" finished segment f = [%g,%g] Hz\n",data->flow,data->flow+data->fstep);
data->flow += data->fstep;
}while(data->fhigh < data->fmax);
......
......@@ -21,7 +21,7 @@
#define zeta 1.0
#define kappa_BL 0.8
#define FSTEP 10.0 // the stencil separation in Hz for the spline model. Going below 2 Hz is dangerous - will fight with line model
#define NSTEP 1000000
#define NSTEP 200000
// as we are getting down to the line width of the broadest lines
......
......@@ -1094,7 +1094,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
{
printf("Setting up time-frequency proposal for IFO%i...\n",ifo);
TFprop_setup(data, model[chain->index[0]], prior->range, tf, ifo);
if(data->gnuplotFlag) TFprop_plot(prior->range, tf, data->Tobs/(double)data->N, ifo);
if(data->gnuplotFlag) TFprop_plot(prior->range, tf, (prior->range[0][1]-prior->range[0][0])/(double)(tf->nt), ifo);
}
if(data->signalFlag)
{
......
......@@ -881,7 +881,7 @@ void initialize_bayesline(struct BayesLineParams **bayesline, struct Data *data,
bayesline[ifo]->priors->sigma = malloc(N*sizeof(double));
//Set BayesLine priors based on the data channel being used
set_bayesline_priors(data->channels[ifo], bayesline[ifo]);
set_bayesline_priors(data->channels[ifo], bayesline[ifo], data->Tobs);
// set default flags
bayesline[ifo]->constantLogLFlag = data->constantLogLFlag;
......@@ -1201,29 +1201,30 @@ void free_background(struct Background *background, int NI, int N)
void initialize_TF_proposal(struct Data *data, struct Prior *prior, struct TimeFrequencyMap *tf)
{
int i,ifo;
int N =data->N;
int ifo;
//int i,ifo;
//int N =data->N;
int NI=data->NI;
tf->N = data->N;
//4s grid for Q
tf->nQ = (int)(prior->range[2][1] - prior->range[2][0]);
tf->nQ = (int)(prior->range[2][1] - prior->range[2][0])/4;
//4 Hz grid for frequency
tf->nf = (int)(prior->range[1][1] - prior->range[1][0])/4;
double tmin = prior->range[0][0];
double tmax = prior->range[0][1];
tf->nt = 0;
double t;
for(i = 0; i < N; i++)
{
t = ((double)(i)+0.5)*(data->Tobs/data->N);
if(t > tmin && t < tmax) tf->nt++;
}
tf->nf = 128;//(int)(prior->range[1][1] - prior->range[1][0])/4;
// double tmin = prior->range[0][0];
// double tmax = prior->range[0][1];
//
// tf->nt = 0;
// double t;
// for(i = 0; i < N; i++)
// {
// t = ((double)(i)+0.5)*(8*data->Tobs/data->N);
// if(t > tmin && t < tmax) tf->nt++;
// }
tf->nt = 1024;
tf->pdf = malloc((NI+1)*sizeof(double ***));
tf->snr = malloc((NI+1)*sizeof(double ***));
......
......@@ -10,7 +10,7 @@
#include "BayesWavePrior.h"
#include "BayesWaveWavelet.h"
void set_bayesline_priors(char *channel, struct BayesLineParams *bayesline)
void set_bayesline_priors(char *channel, struct BayesLineParams *bayesline, double Tobs)
{
//S6 strain
if(!strcmp(channel,"LDAS-STRAIN"))
......@@ -59,11 +59,11 @@ void set_bayesline_priors(char *channel, struct BayesLineParams *bayesline)
//Use O1 strain as default
else
{
bayesline->priors->SAmin = 5.0e-47;
bayesline->priors->SAmax = 5.0e-45;
bayesline->priors->SAmin = 5.0e-47*Tobs/4.0;
bayesline->priors->SAmax = 5.0e-43*Tobs/4.0;
bayesline->priors->LQmin = 1.0e1;
bayesline->priors->LQmax = 1.0e7;
bayesline->priors->LAmin = 1.0e-40;
bayesline->priors->LAmin = 1.0e-45;
bayesline->priors->LAmax = 1.0e-18;
}
......
void set_bayesline_priors(char *channel, struct BayesLineParams *bayesline);
void set_bayesline_priors(char *channel, struct BayesLineParams *bayesline, double Tobs);
int checkrange(double *p, double **r);
......
......@@ -1124,7 +1124,6 @@ double TFprop_density(double *params, double **range, struct TimeFrequencyMap *t
void TFprop_setup(struct Data *data, struct Model *model, double **range, struct TimeFrequencyMap *tf, int ifo)
{
int n = tf->N;
int nt = tf->nt;
int nf = tf->nf;
int nQ = tf->nQ;
......@@ -1135,7 +1134,7 @@ void TFprop_setup(struct Data *data, struct Model *model, double **range, struct
double Tobs = data->Tobs;
int i, j, k, l, m, ii, ioff;
int i, j, k, l, m, ii, ioff,n;
double Qmin, Qmax;
double tmin, tmax;
double fmin, fmax;
......@@ -1144,10 +1143,10 @@ void TFprop_setup(struct Data *data, struct Model *model, double **range, struct
double tsnr;
double norm;
double *d = dvector(0,n-1);
double *h = dvector(0,n-1);
double *AC = dvector(0,n-1);
double *AF = dvector(0,n-1);
double *d = dvector(0,tf->N-1);
double *h = dvector(0,tf->N-1);
double *AC = dvector(0,tf->N-1);
double *AF = dvector(0,tf->N-1);
double *params = dvector(0,4);
......@@ -1163,7 +1162,7 @@ void TFprop_setup(struct Data *data, struct Model *model, double **range, struct
params[3] = 1.0;
params[4] = 0.0;
dt = data->Tobs/(double)data->N;
dt = (tmax-tmin)/(double)(nt);//8*data->Tobs/(double)data->N;
dfx = (fmax-fmin)/(double)(nf); // spacing of grid in frequency
dQx = (Qmax-Qmin)/(double)(nQ); // spacing of grid in frequency
......@@ -1171,10 +1170,13 @@ void TFprop_setup(struct Data *data, struct Model *model, double **range, struct
dV = (double)(nt*nf*nQ)/( (tmax-tmin) * (fmax-fmin) * (Qmax-Qmin) );
// offset in time index between start of data and start of analysis region
ioff = (int)((tmin)/dt);
ioff = (int)((tmin)/(data->Tobs/(double)data->N));
//how many time bins between tmin and tmax?
int nbin = (int)((tmax-tmin)/(data->Tobs/(double)data->N))/nt;
// Copy data into local array (Fourier transform will destroy contents)
for(i=0; i<n/2; i++)
// Copy data into local array (Fourier transform will destroy contents)
for(i=0; i<tf->N/2; i++)
{
d[2*i] = data->s[ifo][2*i];
d[2*i+1]=-data->s[ifo][2*i+1];
......@@ -1192,14 +1194,14 @@ void TFprop_setup(struct Data *data, struct Model *model, double **range, struct
// Step trhough different frequencies
params[1] = fmin + dfx*((double)(j)+0.5);
SineGaussianFourier(h, params, n, 0, Tobs);
SineGaussianFourier(h, params, tf->N, 0, Tobs);
fsnr = fourier_nwip(data->imin, data->imax, h, h, model->invSnf[ifo]);
//TODO: Call phase_blind_time_shift in TFprop_setup
AC[0]=AC[1]=AF[0]=AF[1]=0.0;
for (i=1; i < n/2; i++)
for (i=1; i < tf->N/2; i++)
{
l=2*i;
k=l+1;
......@@ -1211,17 +1213,22 @@ void TFprop_setup(struct Data *data, struct Model *model, double **range, struct
}
drealft(AC-1, n, -1);
drealft(AF-1, n, -1);
drealft(AC-1, tf->N, -1);
drealft(AF-1, tf->N, -1);
// Get filter weight at each time step
for(i = 0; i < nt; i++)
{
ii = i + ioff;
tsnr = (AC[ii]*AC[ii]+AF[ii]*AF[ii])/(fsnr);
// Store SNR for signal-model proposal
tsnr = 0.0;
//integrate snr over nbin
for(n=0; n<nbin; n++)
{
ii = nbin*i + ioff + n;
tsnr += (AC[ii]*AC[ii]+AF[ii]*AF[ii])/(fsnr);
}
// Store SNR for signal-model proposal
tfsnr[m][j][i] = tsnr;
// Proposal density scales as BF~exp(SNR^2/2)
......@@ -1253,10 +1260,10 @@ void TFprop_setup(struct Data *data, struct Model *model, double **range, struct
free_dvector(AC,0,n-1);
free_dvector(AF,0,n-1);
free_dvector(h, 0,n-1);
free_dvector(d, 0,n-1);
free_dvector(AC,0,tf->N-1);
free_dvector(AF,0,tf->N-1);
free_dvector(h, 0,tf->N-1);
free_dvector(d, 0,tf->N-1);
free_dvector(params,0,4);
......@@ -1426,7 +1433,7 @@ void TFprop_signal(struct Data *data, double **range, struct TimeFrequencyMap *t
for(i = 0; i < nt; i++)
{
ii = i - ioff;
if(ii < 0 || ii >= nt)
if(ii < 0 || ii > nt)
{
tsnr = tfsnr[0][m][j][i]; // leave out second detector if mapped out of time window (small edge effect)
}
......
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