Commit 2465d034 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

working time-frequency proposal


git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@589 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent 17644733
......@@ -84,10 +84,6 @@ struct Wavelet
double *templates;
double **intParams;
double **TFhistory;
double **TFdensity;
double TFdensityMax;
char *name;
};
......@@ -124,10 +120,11 @@ struct Model
struct Network
{
double **expPhase;
double *deltaT;
double *Fplus;
double *Fcross;
double **expPhase;
double *deltaT;
double *Fplus;
double *Fcross;
double *dtimes;
};
struct Chain
......@@ -149,6 +146,7 @@ struct Chain
int fcount; // # of intrinsic fisher attempts
int pcount; // # of intrinsic phase attempts
int ccount; // # of cluster proposal attempts
int dcount; // # of density proposal attempts
int ucount; // # of uniform proposal attempts
int macc; // # of MCMC successess
int sacc; // # of RJMCMC successes
......@@ -156,6 +154,7 @@ struct Chain
int facc; // # of intrinsic fisher successes
int pacc; // # of intrinsic phase successes
int cacc; // # of cluster proposal successes
int dacc; // # of cluster proposal successes
int uacc; // # of uniform proposal successes
int burnFlag;// flag telling proposals if we are in burn-in phase
......
......@@ -1022,6 +1022,7 @@ void print_chain_status(struct Data *data, struct Chain *chain, struct Model **m
fprintf(stdout, " RJ: ");
fprintf(stdout, "rjacc=%.1e ", (double)chain->sacc/(double)(chain->scount));
fprintf(stdout, "cacc=%.1e " , (double)chain->cacc/(double)(chain->ccount));
fprintf(stdout, "dacc=%.1e " , (double)chain->dacc/(double)(chain->dcount));
fprintf(stdout, "uacc=%.1e " , (double)chain->uacc/(double)(chain->ucount));
fprintf(stdout, "\n");
fprintf(stdout," INT: ");
......@@ -1498,8 +1499,17 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
//store command line in data structure for printing later on
data->commandLine = commandLine;
//bayesLine takes precedent over PSD fitting
//if(data->bayesLineFlag) data->psdFitFlag = 0;
//BayesLine requires the cleaning phase
if(data->bayesLineFlag == 1 && data->cleanModelFlag == 0)
{
fprintf(stdout,"\n");
fprintf(stdout,"************************ WARNING ************************\n");
fprintf(stdout,"BayesLine requires the cleaning phase \n");
fprintf(stdout,"Ignoring --noClean argument \n");
fprintf(stdout,"If you really want to skip that step remove --bayesLine \n");
fprintf(stdout,"*********************************************************\n");
data->cleanModelFlag = 1;
}
//if doing a short run disable spline thermodynamic integration
if(chain->count < 1000000 && data->splineEvidenceFlag)
......
This diff is collapsed.
......@@ -8,7 +8,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
void PTMCMC(struct Model **model, struct Chain *chain, long *seed, int NC);
void EvolveBayesLineParameters(struct Data *data, struct Model **model, struct BayesLineParams ***bayesline, struct Chain *chain, struct Prior *prior, int ic);
void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Model **model, struct Chain *chain, long *seed, int ic);
void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Model **model, struct Chain *chain, struct TimeFrequencyMap *tf, long *seed, int ic);
void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Model **model, struct Chain *chain, long *seed, int ic);
/* ********************************************************************************** */
......
......@@ -178,25 +178,27 @@ void free_fisher(struct FisherMatrix *fisher)
void initialize_chain(struct Chain *chain, int flag)
{
chain->mod0 = 0; //Counter keeping track of noise model
chain->mod1 = 0; //Counter keeping track of signal model
chain->mod2 = 0; //Counter keeping track of noise model
chain->mod3 = 0; //Counter keeping track of signal model
chain->mcount = 1; //# of MCMC proposals
chain->scount = 1; //# of RJMCMC proposals
chain->xcount = 1; //# of extrinsic proposals
chain->mod0 = 0; //Counter keeping track of noise model
chain->mod1 = 0; //Counter keeping track of signal model
chain->mod2 = 0; //Counter keeping track of noise model
chain->mod3 = 0; //Counter keeping track of signal model
chain->mcount = 1; //# of MCMC proposals
chain->scount = 1; //# of RJMCMC proposals
chain->xcount = 1; //# of extrinsic proposals
chain->fcount = 1; //# of intrinsic fisher proposals
chain->pcount = 1; //# of intrinsic phase proposals
chain->ccount = 1; //# of cluster proposals
chain->dcount = 1; //# of density proposals
chain->ucount = 1; //# of uniform proposals
chain->macc = 0; //# of MCMC successes
chain->sacc = 0; //# of RJMCMC successes
chain->xacc = 0; //# of extrinsic successes
chain->macc = 0; //# of MCMC successes
chain->sacc = 0; //# of RJMCMC successes
chain->xacc = 0; //# of extrinsic successes
chain->facc = 0; //# of intrinsic fisher successes
chain->pacc = 0; //# of intrinsic phase successes
chain->cacc = 0; //# of cluster proposal successes
chain->dacc = 0; //# of density proposal successes
chain->uacc = 0; //# of uniform proposal successes
chain->burnFlag = flag; //Flag for burn-in (1 = yes)
chain->burnFlag = flag; //Flag for burn-in (1 = yes)
}
void free_chain(struct Chain *chain, int NI)
......@@ -919,10 +921,11 @@ void reset_params(struct Wavelet *wave_x, struct Wavelet *wave_y, double **range
void initialize_network(struct Network *projection, int N, int NI)
{
projection->expPhase = dmatrix(0,NI-1,0,N-1);
projection->deltaT = dvector(0,NI-1);
projection->Fplus = dvector(0,NI-1);
projection->Fcross = dvector(0,NI-1);
projection->expPhase = dmatrix(0,NI-1,0,N-1);
projection->deltaT = dvector(0,NI-1);
projection->Fplus = dvector(0,NI-1);
projection->Fcross = dvector(0,NI-1);
projection->dtimes = dvector(0,NI-1);
}
void free_network(struct Network *projection, int N, int NI)
......@@ -1068,10 +1071,6 @@ void initialize_wavelet(struct Wavelet *wave, int N, int smax, int size, int fsi
wave->cosy = 0.0;
wave->logp = 0.0;
//density proposal
wave->TFhistory = dmatrix(0,fsize-1,0,tsize-1);
wave->TFdensity = dmatrix(0,fsize-1,0,tsize-1);
}
void copy_wavelet(struct Wavelet *origin, struct Wavelet *copy, int N)
......@@ -1098,8 +1097,6 @@ void free_wavelet(struct Wavelet *wave, int N, int smax, int fsize, int tsize)
free_dmatrix(wave->intParams,1,smax,0,4);
free_ivector(wave->index,1,smax);
free_dvector(wave->templates,0,N-1);
free_dmatrix(wave->TFhistory,0,fsize-1,0,tsize-1);
free_dmatrix(wave->TFdensity,0,fsize-1,0,tsize-1);
}
void initialize_background(struct Background *background, int NI, int N)
......@@ -1146,53 +1143,61 @@ void free_background(struct Background *background, int NI, int N)
free_dvector(background->spectrum,0,N-1);
}
void initialize_TFproposal(struct Data *data, struct Prior *prior, struct Wavelet *wave)
void initialize_TF_proposal(struct Data *data, struct Prior *prior, struct TimeFrequencyMap *tf)
{
//Initialize TFdensity to uniform over allowed TF volume
int i,j;
int count = 0;
double t,f;
double dt = data->Tobs/(double)data->tsize;
double df = data->fmax/(double)data->fsize;
int i,ifo;
int N =data->N;
int NI=data->NI;
tf->N = data->N;
/*
prior->range[0][0] = tmin;
prior->range[0][1] = tmax;
prior->range[1][0] = fmin;
prior->range[1][1] = fmax;
*/
//1/2s grid for Q
tf->nQ = (int)(prior->range[2][1] - prior->range[2][0])*2;
double tmin = prior->range[0][0];//dt*(int)((data->Tobs/16.)/dt);
double tmax = prior->range[0][1];//data->Tobs-tmin;
double fmin = prior->range[1][0];//data->fmin;
double fmax = prior->range[1][1];//data->fmax;
//1 Hz grid for frequency
tf->nf = (int)(prior->range[1][1] - prior->range[1][0]);
double tmin = prior->range[0][0];
double tmax = prior->range[0][1];
for(i=0; i<data->fsize; i++)
{
f = i*df;
for(j=0; j<data->tsize; j++)
{
t = j*dt;
if(f>=fmin && f<fmax && t>=tmin && t<tmax)
{
wave->TFdensity[i][j] = 1.0;
count++;
}
else wave->TFdensity[i][j] = 0.0;
}
}
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++;
}
wave->TFdensityMax = 1./(double)count/0.5;
for(i=0; i<data->fsize; i++)
{
for(j=0; j<data->tsize;j++)
{
wave->TFhistory[i][j]=wave->TFdensity[i][j];
wave->TFdensity[i][j]/=(double)count/0.5;
}
}
tf->pdf = malloc((NI+1)*sizeof(double ***));
tf->snr = malloc((NI+1)*sizeof(double ***));
tf->max = malloc((NI+1)*sizeof(double *));
for(ifo=0; ifo<NI+1; ifo++)
{
tf->pdf[ifo] = d3tensor(0,tf->nQ-1,0,tf->nf-1,0,tf->nt-1);
tf->snr[ifo] = d3tensor(0,tf->nQ-1,0,tf->nf-1,0,tf->nt-1);
tf->max[ifo] = malloc(tf->nQ*sizeof(double));
}
}
void free_TF_proposal(struct Data *data, struct TimeFrequencyMap *tf)
{
int ifo;
for(ifo=0; ifo<data->NI+1; ifo++)
{
free_d3tensor(tf->pdf[ifo],0,tf->nQ-1,0,tf->nf-1,0,tf->nt-1);
free_d3tensor(tf->snr[ifo],0,tf->nQ-1,0,tf->nf-1,0,tf->nt-1);
free(tf->max[ifo]);
}
free(tf->pdf);
free(tf->snr);
free(tf->max);
}
int checkfile(char *name)
......
......@@ -50,11 +50,12 @@ void initialize_background(struct Background *background, int NI, int N);
void copy_background (struct Background *origin, struct Background *copy, int NI, int N);
void free_background (struct Background *background, int NI, int N);
void initialize_TFproposal(struct Data *data, struct Prior *prior, struct Wavelet *wave);
void initialize_network(struct Network *projection, int N, int NI);
void free_network(struct Network *projection, int N, int NI);
void initialize_TF_proposal(struct Data *data, struct Prior *prior, struct TimeFrequencyMap *tf);
void free_TF_proposal(struct Data *data, struct TimeFrequencyMap *tf);
int checkfile(char *name);
int *ivector(long nl, long nh);
......
......@@ -750,15 +750,6 @@ double wavlet_proposal_density(int ifo, double f, double t, struct Data *data)
i = (int)(f/df);
j = (int)(t/dt);
/*
ifo coming from EvolveIntrinsicParameters
ranges from [-1:NI-1], w/ -1 corresponding
to geocenter.
TFdensity ranges from [0:NI] w/ 0 corresponding
to geocenter
*/
//return(log(data->TFdensity[ifo+1][i][j]));
return(-log((double)(data->N)));
}
......@@ -1024,6 +1015,7 @@ void wavelet_sample(double *params, double **range, double **prop, int tsize, do
return;
}
//void TFdraw(double *params, double **range, double dt, double df, double **tfden, double tfmax, long *seed)
......@@ -1040,31 +1032,408 @@ void draw_wavelet_params(int ifo, double *params, struct Data *data, double **ra
return;
}
void draw_time_frequency(int ifo, int ii, struct Wavelet *wx, struct Wavelet *wy, struct Data *data, double **range, long *seed, double fr)
void draw_time_frequency(int ifo, int ii, struct Wavelet *wx, struct Wavelet *wy, struct Data *data, double **range, long *seed, double fr, struct TimeFrequencyMap *tf, int *prop)
{
if(ran2(seed)<fr)
wavelet_sample(wy->intParams[ii],range,wx->TFdensity,data->tsize, data->Tobs, wx->TFdensityMax,seed);
{
TFprop_draw(wy->intParams[ii], range, tf, ifo, seed);
*prop=0;
}
else
{
wavelet_proximity_new(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;
double tmin = range[0][0];
double tmax = range[0][1];
double fmin = range[1][0];
double fmax = range[1][1];
double Qmin = range[2][0];
double Qmax = range[2][1];
double dt = (tmax-tmin)/(double)tf->nt;
double df = (fmax-fmin)/(double)tf->nf;
double dQ = (Qmax-Qmin)/(double)tf->nQ;
int k = (int)((params[2]-Qmin)/dQ);
do
{
f = fmin + (fmax-fmin)*ran2(seed);
t = tmin + (tmax-tmin)*ran2(seed);
i = (int)((f-fmin)/df);
j = (int)((t)/dt);
alpha = tf->max[ifo][k]*ran2(seed);
} while(alpha > tf->pdf[ifo][k][i][j]);
double TFprop(double f, double t, double **prop, int tsize, int fsize, double Tobs)
params[0] = t;
params[1] = f;
}
double TFprop_density(double *params, double **range, struct TimeFrequencyMap *tf, int ifo)
{
int i, j;
double dt, df, alpha;
double x;
dt = Tobs/(double)tsize;
df = 0.5/dt;
double tmin = range[0][0];
double tmax = range[0][1];
double fmin = range[1][0];
double fmax = range[1][1];
double Qmin = range[2][0];
double Qmax = range[2][1];
i = (int)(f/df);
j = (int)(t/dt);
double dt = (tmax-tmin)/(double)tf->nt;
double df = (fmax-fmin)/(double)tf->nf;
double dQ = (Qmax-Qmin)/(double)tf->nQ;
int i = (int)((params[0])/dt);
int j = (int)((params[1]-fmin)/df);
int k = (int)((params[2]-Qmin)/dQ);
x = tf->pdf[ifo][k][j][i];
return(x);
}
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 *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);
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
// offset in time index between start of data and start of analysis region
ioff = (int)((tmin)/dt);
AC=dvector(0,n-1); AF=dvector(0,n-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];
}
double max;
double tmaxx;
int l,k;
for(m = 0; m < nQ; m++) // loop over Q
{
Q = Qmin + ((double)(m)+0.5)*dQx;
params[2] = Q;
x = 0.0;
for(j = 0; j < nf; j++)
{
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];
}
alpha = -1.0e60;
drealft(AC-1, n, -1);
drealft(AF-1, n, -1);
if((i > -1) && (j > -1) && (i < fsize) && (j < tsize)) alpha = prop[i][j];
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)
{
max=y;
tmaxx=ii*dt;
}
}
}
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);
}
void TFprop_plot(double **range, struct TimeFrequencyMap *tf, double dt, int ifo)
{
int nt= tf->nt;
int nf= tf->nf;
int nQ= tf->nQ;
double ***tfden = tf->pdf[ifo];
int m, j, i;
double Qmin, Qmax;
double tmin, tmax;
double fmin, fmax;
double f, t, dfx, dQx, a, x, y, u, Q, mx;
char filename[100];
FILE *out;
FILE *test;
tmin = range[0][0];
tmax = range[0][1];
fmin = range[1][0];
fmax = range[1][1];
Qmin = range[2][0];
Qmax = range[2][1];
dfx = (fmax-fmin)/(double)(nf); // spacing of grid in frequency
dQx = (Qmax-Qmin)/(double)(nQ); // spacing of grid in frequency
x = 0.0;
y = 1.0e60;
for(m = 0; m < nQ; m++) // loop over Q
{
for(j = 0; j < nf; j++)
{
for(i = 0; i < nt; i++)
{
a = tfden[m][j][i];
if(a > x) x = a;
if(a < y) y = a;
}
}
}
mx = x;
a = ceil(-log10(x));
y = x*pow(10.0, a);
u = floor(y);
if(y-u < 0.5) x = (u+0.5)*pow(10.0, -a);
if(y-u >= 0.5) x = (u+1.0)*pow(10.0, -a);
for(m = 0; m < nQ; m++) // loop over Q
{
Q = Qmin + ((double)(m)+0.5)*dQx;
return(alpha);
sprintf(filename,"waveforms/TF_map.dat");
out = fopen(filename,"w");
for(j = 0; j < nf; j++)
{