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

Tested TF proposals for RJMCMC

Generalized Makefile


git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@595 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent 0c4da0d3
......@@ -612,7 +612,7 @@ int main(int argc, char *argv[])
/******************************************************************************/
//if(sqrt(rSNR)>0.0) resize_model(data, chain, prior, model, bayesline, psd, chain->NC+5);
if(sqrt(rSNR)>80.0) chain->NC+=5;
if(sqrt(rSNR)>50.0) chain->NC+=5;
/******************************************************************************/
/* */
......
......@@ -814,11 +814,6 @@ void print_run_flags(FILE *fptr, struct Data *data, struct Prior *prior)
else fprintf(fptr, "DISABLED");
fprintf(fptr, "\n");
fprintf(fptr, " noise realization is ..... ");
if(data->noiseSimFlag) fprintf(fptr, "INCLUDED");
else fprintf(fptr, "NOT INCLUDED");
fprintf(fptr, "\n");
fprintf(fptr, " PSD fitting is ........... ");
if(data->psdFitFlag) fprintf(fptr, "ENABLED");
else fprintf(fptr, "DISABLED");
......@@ -1404,7 +1399,7 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
else data->clusterProposalFlag = 0;
ppt = LALInferenceGetProcParamVal(commandLine, "--clusterWeight");
if(ppt) data->TFtoProx = (double)atof(ppt->value);
if(ppt) data->TFtoProx = 1.0-(double)atof(ppt->value);
else data->TFtoProx = 0.0;
ppt = LALInferenceGetProcParamVal(commandLine, "--runName");
......
......@@ -990,7 +990,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
if( (data->glitchFlag || data->signalFlag) && data->rjFlag)
{
if(data->constantLogLFlag)
chain->rjmcmcRate = 0.1;
chain->rjmcmcRate = 0.01;
else
chain->rjmcmcRate = 0.5;
}
......@@ -1006,7 +1006,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
if( (data->glitchFlag || data->signalFlag) && data->rjFlag)
{
if(data->constantLogLFlag)
chain->rjmcmcRate = 0.1;
chain->rjmcmcRate = 0.01;
else
chain->rjmcmcRate = 0.5;
}
......@@ -1139,8 +1139,11 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
TFprop_setup(data, model[0], prior->range, tf, ifo);
if(data->gnuplotFlag) TFprop_plot(prior->range, tf, data->Tobs/(double)data->N, ifo);
}
printf("Setting up time-frequency proposal for signal model...\n");
TFprop_signal(data, tf, model[0]->projection);
if(data->signalFlag)
{
printf("Setting up time-frequency proposal for signal model...\n");
TFprop_signal(data, prior->range, tf, model[0]->projection);
}
}
......@@ -1160,7 +1163,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
if( (data->glitchFlag || data->signalFlag) && data->rjFlag)
{
if(data->constantLogLFlag)
chain->rjmcmcRate = 0.1;
chain->rjmcmcRate = 0.01;
else
chain->rjmcmcRate = 0.5;
}
......@@ -1176,7 +1179,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
if( (data->glitchFlag || data->signalFlag) && data->rjFlag)
{
if(data->constantLogLFlag)
chain->rjmcmcRate = 0.1;
chain->rjmcmcRate = 0.01;
else
chain->rjmcmcRate = 0.5;
}
......@@ -1234,7 +1237,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
//After so many iterations recompute the residuals and likelihood (prevent accumulation of roundoff error)
recompute_residual(data, model, chain);
if(burnFlag==0 && chain->mc%1000==0 && data->signalFlag) TFprop_signal(data, tf, model[chain->index[0]]->projection);
if(burnFlag==0 && chain->mc%1000==0 && data->signalFlag) TFprop_signal(data, prior->range, tf, model[chain->index[0]]->projection);
//During cleaning phase track max and min BL PSDs to set up priors for model selection run
......@@ -1341,7 +1344,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
for(ic=0; ic<NC; ic++) chain->logLchain[ic][chain->zcount] = model[chain->index[ic]]->logL + model[chain->index[ic]]->logLnorm - data->logLc;
chain->zcount++;
if(chain->mc%(M/100)==0)TrapezoidIntegration(chain, chain->zcount, modelname, &logZ, &varZ);
if(chain->mc%(M/100)==0 && !data->constantLogLFlag)TrapezoidIntegration(chain, chain->zcount, modelname, &logZ, &varZ);
}
chain->mc+=chain->cycle;
......@@ -1623,9 +1626,14 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
//ratio of TF proposals to proximity proposals
double TFtoProx = data->TFtoProx;
double ClusterRate = 0.9;
double FisherRate = 0.8;
if(data->constantLogLFlag) FisherRate = 0.2;
//optimize ratio for sampling prior
if(data->constantLogLFlag || ic > 10)
{
FisherRate = 0.2;
ClusterRate = 0.2;
}
int burn = chain->burnFlag;
......@@ -2090,6 +2098,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
wave_y->logp += -(double)(wave_y->size)*prior->logTFV;
}
//amplitude prior
if(data->amplitudePriorFlag)
{
......@@ -2126,12 +2135,13 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
// and finally, the Hastings ratio
logH = (model_y->logL - model_x->logL)*chain->beta + wave_y->logp - wave_x->logp - logqy + logqx;
/*
if(model_x->size < model_y->size)
/*
if(model_x->size != model_y->size)
{
printf("%i: %i->%i logLy=%g logLx=%g py=%g px=%g qy=%g qx=%g, dx=%i,dy=%i,num=%g, den=%g, logH=%g\n",typ,model_x->size,model_y->size,model_y->logL,model_x->logL,wave_y->logp,wave_x->logp,logqy,logqx,wave_x->size,wave_y->size,wave_y->logp+ logqx,- wave_x->logp - logqy, logH);
}
*/
*/
//if(chain->index[ic]==0 && data->stochasticFlag)printf("logLy=%g logLx=%g Ay=%g Ax=%g detC=%g logH=%g\n",model_y->logL-data->logLc,model_x->logL-data->logLc,model_x->background->logamp,model_y->background->logamp, model_y->background->detCij[data->N/4],logH);
......@@ -2139,7 +2149,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
with the maximization on the sytsem effectively has fewer degrees of freedom per SG
compensate with a prior penalty
*/
if(burn)
if(burn && !data->constantLogLFlag)
{
if(det==-1) logH += 10.0*(double)(wave_x->size-wave_y->size);//*chain->beta;
else logH += 10.0*(double)(wave_x->size-wave_y->size);//*chain->beta;
......
......@@ -402,7 +402,7 @@ void initialize_priors(struct Data *data, struct Prior *prior, int omax)
prior->range[4][0] = 0.0;
prior->range[4][1] = LAL_TWOPI;
prior->TFV = (tmax-tmin)*(fmax-fmin);
prior->TFV = (tmax-tmin)*(fmax-fmin)*(Qmax-Qmin);
prior->logTFV = log(prior->TFV);
......@@ -589,7 +589,7 @@ void reset_priors(struct Data *data, struct Prior *prior)
prior->range[4][0] = 0.0;
prior->range[4][1] = LAL_TWOPI;
prior->TFV = (tmax-tmin)*(fmax-fmin);
prior->TFV = (tmax-tmin)*(fmax-fmin)*(Qmax-Qmin);
prior->logTFV = log(prior->TFV);
}
......@@ -1151,11 +1151,11 @@ void initialize_TF_proposal(struct Data *data, struct Prior *prior, struct TimeF
tf->N = data->N;
//1/2s grid for Q
tf->nQ = (int)(prior->range[2][1] - prior->range[2][0])*2;
//4s grid for Q
tf->nQ = (int)(prior->range[2][1] - prior->range[2][0])/4;
//1 Hz grid for frequency
tf->nf = (int)(prior->range[1][1] - prior->range[1][0]);
//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];
......
......@@ -202,7 +202,9 @@ void proximity_normalization(double **params, double **range, double *norm, int
double tau;
double tmin, tmax, sigma_t_thin, sigma_t_fat;
double fmin, fmax, sigma_f_thin, sigma_f_fat;
double QV;
QV = range[2][1]-range[2][0];
for(i=1; i<= cnt; i++)
{
......@@ -222,6 +224,7 @@ void proximity_normalization(double **params, double **range, double *norm, int
fmax = (range[1][1]-params[j][1]);
norm[j] = 1./(gaussian_norm(tmin,tmax,sigma_t_fat)*gaussian_norm(fmin,fmax,sigma_f_fat) - gaussian_norm(tmin,tmax,sigma_t_thin)*gaussian_norm(fmin,fmax,sigma_f_thin));
norm[j] /= QV;
}
}
......
......@@ -1042,14 +1042,15 @@ void draw_time_frequency(int ifo, int ii, struct Wavelet *wx, struct Wavelet *wy
else
{
wavelet_proximity_new(wy->intParams[ii],wx->intParams,wy->index,range,wx->size, seed);
//wavelet_proximity(wy->intParams[ii],wx->intParams,wy->index,range,wx->size, seed);
*prop=1;
}
}
void TFprop_draw(double *params, double **range, struct TimeFrequencyMap *tf, int ifo, long *seed)
{
int j, i;
double f, t, alpha;
int j, i, k;
double f, t, Q, alpha;
double tmin = range[0][0];
double tmax = range[0][1];
......@@ -1062,16 +1063,15 @@ void TFprop_draw(double *params, double **range, struct TimeFrequencyMap *tf, in
double df = (fmax-fmin)/(double)tf->nf;
double dQ = (Qmax-Qmin)/(double)tf->nQ;
int k = (int)((params[2]-Qmin)/dQ);
do
{
Q = Qmin + (Qmax-Qmin)*ran2(seed);
f = fmin + (fmax-fmin)*ran2(seed);
t = tmin + (tmax-tmin)*ran2(seed);
i = (int)((f-fmin)/df);
j = (int)((t)/dt);
k = (int)((Q-Qmin)/dQ);
alpha = tf->max[ifo][k]*ran2(seed);
......@@ -1079,6 +1079,7 @@ void TFprop_draw(double *params, double **range, struct TimeFrequencyMap *tf, in
params[0] = t;
params[1] = f;
params[2] = Q;
}
double TFprop_density(double *params, double **range, struct TimeFrequencyMap *tf, int ifo)
......@@ -1110,143 +1111,138 @@ 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;
double ***tfden = tf->pdf[ifo];
double ***tfsnr = tf->snr[ifo];
double *tfmax = tf->max[ifo];
double Tobs = data->Tobs;
int i, j, k, l, m, ii, ioff;
double Qmin, Qmax;
double tmin, tmax;
double fmin, fmax;
double dt, df, dfx, dQx, dV;
double fsnr;
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 *params = dvector(0,4);
tmin = range[0][0];
tmax = range[0][1];
fmin = range[1][0];
fmax = range[1][1];
Qmin = range[2][0];
Qmax = range[2][1];
//Fix time, amplitude, and phse of wavelet filter
params[0] = 0.0;
params[3] = 1.0;
params[4] = 0.0;
dt = data->Tobs/(double)data->N;
df = data->df;
dfx = (fmax-fmin)/(double)(nf); // spacing of grid in frequency
dQx = (Qmax-Qmin)/(double)(nQ); // spacing of grid in frequency
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);
// Copy data into local array (Fourier transform will destroy contents)
for(i=0; i<n/2; i++)
{
d[2*i] = data->s[ifo][2*i];
d[2*i+1]=-data->s[ifo][2*i+1];
}
norm = 0.0;
int n = tf->N;
int nt = tf->nt;
int nf = tf->nf;
int nQ = tf->nQ;
double ***tfden = tf->pdf[ifo];
double ***tfsnr = tf->snr[ifo];
double *tfmax = tf->max[ifo];
double *a = dvector(0,n-1);//data->s[ifo];
double Tobs = data->Tobs;
int i, ii, j, m, ioff;
double Qmin, Qmax;
double tmin, tmax;
double fmin, fmax;
double f, dt, df, dfx, dQx, x, y, Q;
double *AC, *AF;
double *b;
double *params;
double bsq;
tmin = range[0][0];
tmax = range[0][1];
fmin = range[1][0];
fmax = range[1][1];
Qmin = range[2][0];
Qmax = range[2][1];
// [0] t0 [1] f0 [2] Q [3] Amp [4] phi
params= dvector(0,4);
// Loop over different Q layers
for(m = 0; m < nQ; m++)
{
// Set Q of wavelet filter to median of layer
params[2] = Qmin + ((double)(m)+0.5)*dQx;
params[0] = 0.0;
params[3] = 1.0;
params[4] = 0.0;
for(j = 0; j < nf; j++)
{
// Step trhough different frequencies
params[1] = fmin + dfx*((double)(j)+0.5);
dt = data->Tobs/(double)data->N;
df = data->df;
SineGaussianFourier(h, params, n, 0, Tobs);
dfx = (fmax-fmin)/(double)(nf); // spacing of grid in frequency
dQx = (Qmax-Qmin)/(double)(nQ); // spacing of grid in frequency
fsnr = fourier_nwip(data->imin, data->imax, h, h, model->invSnf[ifo]);
// offset in time index between start of data and start of analysis region
ioff = (int)((tmin)/dt);
//TODO: Call phase_blind_time_shift in TFprop_setup
AC[0]=AC[1]=AF[0]=AF[1]=0.0;
AC=dvector(0,n-1); AF=dvector(0,n-1);
for (i=1; i < n/2; i++)
{
l=2*i;
k=l+1;
b = dvector(0,n-1);
for(i=0; i<n/2; i++)
{
a[2*i] = data->s[ifo][2*i];
a[2*i+1]=-data->s[ifo][2*i+1];
}
AC[l] = ( d[l]*h[l] + d[k]*h[k]) * model->invSnf[ifo][i];
AC[k] = ( d[k]*h[l] - d[l]*h[k]) * model->invSnf[ifo][i];
AF[l] = ( d[l]*h[k] - d[k]*h[l]) * model->invSnf[ifo][i];
AF[k] = ( d[k]*h[k] + d[l]*h[l]) * model->invSnf[ifo][i];
double max;
double tmaxx;
int l,k;
for(m = 0; m < nQ; m++) // loop over Q
{
}
Q = Qmin + ((double)(m)+0.5)*dQx;
drealft(AC-1, n, -1);
drealft(AF-1, n, -1);
params[2] = Q;
// 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);
x = 0.0;
// Store SNR for signal-model proposal
tfsnr[m][j][i] = tsnr;
for(j = 0; j < nf; j++)
// Proposal density scales as BF~exp(SNR^2/2)
tfden[m][j][i] = exp(tsnr/2.0)-1.0;
norm += tfden[m][j][i];
}
}
} // end Q loop
//Normalize and find max density for rejection sampling
tfmax[0] = 0;
for(m = 0; m < nQ; m++)
{
max=0.0;
f = fmin + dfx*((double)(j)+0.5);
params[1] = f;
SineGaussianFourier(b, params, n, 0, Tobs);
bsq = fourier_nwip(data->imin, data->imax, b, b, 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++)
{
l=2*i;
k=l+1;
AC[l] = ( a[l]*b[l] + a[k]*b[k]) / model->Snf[ifo][i];
AC[k] = ( a[k]*b[l] - a[l]*b[k]) / model->Snf[ifo][i];
AF[l] = ( a[l]*b[k] - a[k]*b[l]) / model->Snf[ifo][i];
AF[k] = ( a[k]*b[k] + a[l]*b[l]) / model->Snf[ifo][i];
}
drealft(AC-1, n, -1);
drealft(AF-1, n, -1);
for(i = 0; i < nt; i++)
{
ii = i + ioff;
y = (AC[ii]*AC[ii]+AF[ii]*AF[ii])/(bsq);
tfsnr[m][j][i] = y;
y = exp(y/2.0)-1.0; // BF scales as exp(SNR^2/2)
tfden[m][j][i] = y;
x += y;
if(y>max)
for(j = 0; j < nf; j++)
{
max=y;
tmaxx=ii*dt;
for(i = 0; i < nt; i++)
{
tfden[m][j][i] /= norm;
tfden[m][j][i] *= dV;
if(tfden[m][j][i] > tfmax[0]) tfmax[0] = tfden[m][j][i];
}
}
}
}
y = 0.0;
for(j = 0; j < nf; j++)
{
for(i = 0; i < nt; i++)
{
tfden[m][j][i] /= x;
if(tfden[m][j][i] > y) y = tfden[m][j][i];
}
}
tfmax[m] = y;
} // end Q loop
free_dvector(AC,0,n-1);
free_dvector(AF,0,n-1);
free_dvector(b, 0,n-1);
free_dvector(a, 0,n-1);
for(m = 1; m < nQ; m++) tfmax[m] = tfmax[0];
free_dvector(AC,0,n-1);
free_dvector(AF,0,n-1);
free_dvector(h, 0,n-1);
free_dvector(d, 0,n-1);
}
void TFprop_plot(double **range, struct TimeFrequencyMap *tf, double dt, int ifo)
......@@ -1361,79 +1357,89 @@ void TFprop_plot(double **range, struct TimeFrequencyMap *tf, double dt, int ifo
}
void TFprop_signal(struct Data *data, struct TimeFrequencyMap *tf, struct Network *net)
void TFprop_signal(struct Data *data, double **range, struct TimeFrequencyMap *tf, struct Network *net)
{
int n = tf->N;
int nt = tf->nt;
int nf = tf->nf;
int nQ = tf->nQ;
int n = tf->N;
int nt = tf->nt;
int nf = tf->nf;
int nQ = tf->nQ;
double ****tfden = tf->pdf;
double ****tfsnr = tf->snr;
double **tfmax = tf->max;
double ****tfden = tf->pdf;
double ****tfsnr = tf->snr;
double **tfmax = tf->max;
/*
double ****pdf;
double ****snr;
double **max;
*/
double Tobs = data->Tobs;
double Tobs = data->Tobs;
int i, ii, j, m, ioff;
double dt,dV;
double norm;
double tsnr;
int i, ii, j, m, ioff;
double dt, x, y;
dt = Tobs/(double)n;
dt = Tobs/(double)n;
double tmin, tmax;
double fmin, fmax;
double Qmin, Qmax;
tmin = range[0][0];
tmax = range[0][1];
fmin = range[1][0];
fmax = range[1][1];
Qmin = range[2][0];
Qmax = range[2][1];
// Here we combine the TF maps from two detectors allowing for the time delay tdelay = t1 -t2
// It is assumed that the signal model times are referenced to detector 1
// This code should be generalized to handle N detectors
dV = (double)(nt*nf*nQ)/( (tmax-tmin) * (fmax-fmin) * (Qmax-Qmin) );
int ifo;
// Here we combine the TF maps from two detectors allowing for the time delay tdelay = t1 -t2
// It is assumed that the signal model times are referenced to detector 1
// This code should be generalized to handle N detectors
for(m = 0; m < nQ; m++) // loop over Q
{
int ifo;
x = 0.0;
for(ifo=1; ifo<data->NI; ifo++)
norm=0.0;
for(m = 0; m < nQ; m++) // loop over Q
{
ioff = (int)((net->dtimes[ifo]-net->dtimes[0])/dt);
for(j = 0; j < nf; j++)
{
for(i = 0; i < nt; i++)
// x = 0.0;
for(ifo=1; ifo<data->NI; ifo++)
{
ii = i - ioff;
if(ii < 0 || ii >= nt)
{
y = tfsnr[0][m][j][i]; // leave out second detector if mapped out of time window (small edge effect)
}
else
{
y = tfsnr[0][m][j][i] + tfsnr[ifo][m][j][ii];
}