Commit 64637a81 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

updates to fisher and proximity proposals


git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@575 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent c5b9ac16
......@@ -146,9 +146,17 @@ struct Chain
int mcount; // # of MCMC attempts
int scount; // # of RJMCMC attempts
int xcount; // # of extrinsic attempts
int fcount; // # of intrinsic fisher attempts
int pcount; // # of intrinsic phase attempts
int ccount; // # of cluster proposal attempts
int ucount; // # of uniform proposal attempts
int macc; // # of MCMC successess
int sacc; // # of RJMCMC successes
int xacc; // # of extrinsic successes
int facc; // # of intrinsic fisher successes
int pacc; // # of intrinsic phase successes
int cacc; // # of cluster proposal successes
int uacc; // # of uniform proposal successes
int burnFlag;// flag telling proposals if we are in burn-in phase
......
......@@ -1014,10 +1014,20 @@ void print_chain_status(struct Data *data, struct Chain *chain, struct Model **m
//print intrinsic parameter stats
fprintf(stdout, " INT: ");
fprintf(stdout, " DIM: ");
fprintf(stdout, "DGW=%d ", model[chain->index[0]]->signal->size);
for(ifo=0; ifo<NI; ifo++) fprintf(stdout, "D%i=%d ", ifo, model[chain->index[0]]->glitch[ifo]->size);
fprintf(stdout, "m0=%i m1=%i m2=%i m3=%i rjacc=%.1e macc=%.1e",chain->mod0, chain->mod1, chain->mod2, chain->mod3, (double)chain->sacc/(double)(chain->scount), (double)chain->macc/(double)(chain->mcount));
fprintf(stdout, "m0=%i m1=%i m2=%i m3=%i ", chain->mod0, chain->mod1, chain->mod2, chain->mod3);
fprintf(stdout, "\n");
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, "uacc=%.1e " , (double)chain->uacc/(double)(chain->ucount));
fprintf(stdout, "\n");
fprintf(stdout," INT: ");
fprintf(stdout, "macc=%.1e " , (double)chain->macc/(double)(chain->mcount));
fprintf(stdout, "facc=%.1e " , (double)chain->facc/(double)(chain->fcount));
fprintf(stdout, "pacc=%.1e " , (double)chain->pacc/(double)(chain->pcount));
fprintf(stdout, "\n");
......@@ -1392,7 +1402,7 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
ppt = LALInferenceGetProcParamVal(commandLine, "--clusterWeight");
if(ppt) data->TFtoProx = (double)atof(ppt->value);
else data->TFtoProx = 1.0;
else data->TFtoProx = 0.0;
ppt = LALInferenceGetProcParamVal(commandLine, "--runName");
if(ppt) sprintf(data->runName,"%s_",ppt->value);
......
......@@ -1638,14 +1638,17 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
double logH;
double Tobs = data->Tobs;
int fisherFlag;
int phaseFlag;
int clusterFlag;
int uniformFlag;
//ratio of TF proposals to proximity proposals
double TFtoProx = data->TFtoProx;
double FisherRate = 0.8;
if(data->constantLogLFlag) FisherRate = 0.2;
double scale = 1./sqrt(6.);//??
int burn = chain->burnFlag;
//Unpack structures and use convenient (and familiar) names
......@@ -1700,6 +1703,12 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
for(mc=1; mc<=M; mc++)
{
//Flags for which proposal is being used
uniformFlag=0;
clusterFlag=0;
fisherFlag=0;
phaseFlag=0;
//acc gets set to 1 if model_y gets accepted (no need to copy at next iteration)
if(acc==0) copy_int_model(model_x, model_y, N, data->NI, det);
if(data->psdFitFlag) for(ifo=0; ifo<NI; ifo++) Sny[ifo] = Snx[ifo];
......@@ -1849,8 +1858,14 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
draw_wavelet_params(det, paramsy[ii], data, range, seed);
// TF density for proposal
if(data->clusterProposalFlag)
if(data->clusterProposalFlag && (ran2(seed) < 0.9/chain->temperature[ic]) )
{
if(chain->index[ic]==0)
{
clusterFlag=1;
chain->ccount++;
}
// Draw new time and frequency from TFV density proposal
draw_time_frequency(det, ii, wave_x, wave_y, data, range, seed, TFtoProx);
......@@ -1860,6 +1875,12 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
// Uniform in TF proposal
else
{
if(chain->index[ic]==0)
{
uniformFlag=1;
chain->ucount++;
}
logqy += -prior->logTFV;
logqx += 0.0;
}
......@@ -1965,15 +1986,27 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
*/
if(ran2(seed)<0.1)
{
if(chain->index[ic]==0)
{
phaseFlag=1;
chain->pcount++;
}
intrinsic_halfcycle_proposal(paramsx[ii],paramsy[ii],seed);
}
else if(ran2(seed)<FisherRate)
{
if(chain->index[ic]==0)
{
fisherFlag=1;
chain->fcount++;
}
//TODO: Try new fisher proposal
if(det==-1) intrinsic_fisher_proposal(range, paramsx[ii], paramsy[ii], model_x->SnGeo, 1.0, &logqx, &logqy, scale, seed, Tobs, &test);
else intrinsic_fisher_proposal(range, paramsx[ii], paramsy[ii], model_x->Snf[det], 1.0, &logqx, &logqy, scale, seed, Tobs, &test);
// if(det==-1) intrinsic_fisher_proposal_2(model->x,range, paramsx[ii], paramsy[ii], &logqx, &logqy, seed, &test);
// else intrinsic_fisher_proposal_2(model->x,range, paramsx[ii], paramsy[ii], &logqx, &logqy, seed, &test);
if(det==-1) intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->SnGeo, 1.0, &logqx, &logqy, seed, Tobs, &test);
else intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->Snf[det], 1.0, &logqx, &logqy, seed, Tobs, &test);
//if(det==-1) intrinsic_fisher_proposal_2(model_x,range, paramsx[ii], paramsy[ii], model_x->SnGeo, 1.0, &logqx, &logqy, seed, Tobs, &test);
//else intrinsic_fisher_proposal_2(model_x,range, paramsx[ii], paramsy[ii], model_x->Snf[det], 1.0, &logqx, &logqy, seed, Tobs, &test);
}
else
......@@ -2127,6 +2160,11 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
{
if(wave_x->size>0)chain->macc++;
}
if(fisherFlag) chain->facc++;
if(phaseFlag) chain->pacc++;
if(clusterFlag)chain->cacc++;
if(uniformFlag)chain->uacc++;
}
wave_x->logp = wave_y->logp;
}
......@@ -2284,7 +2322,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
Don't bother for hot chains, which should largely be
sampling from the prior
*/
if(ran2(seed)< 0.5 && ic<10 && NI>1)
if(ran2(seed)< 0.5 && chain->index[ic]<10 && NI>1)
{
if(ic==0)cnt++;
if(ic==0)sky=1;
......
......@@ -185,9 +185,17 @@ void initialize_chain(struct Chain *chain, int flag)
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->ucount = 1; //# of uniform proposals
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->uacc = 0; //# of uniform proposal successes
chain->burnFlag = flag; //Flag for burn-in (1 = yes)
}
......
......@@ -21,13 +21,18 @@
#include "BayesWaveProposal.h"
#include "BayesWaveLikelihood.h"
#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif
/* ********************************************************************************** */
/* */
/* Prior proposals */
/* */
/* ********************************************************************************** */
void draw_glitch_amplitude(double *params, double *Snf, double Sa, long *seed, double Tobs, double **range, double SNRpeak)
{
int i, k;
......@@ -79,20 +84,6 @@ void draw_glitch_amplitude(double *params, double *Snf, double Sa, long *seed, d
alpha = ran2(seed);
k++;
// if(k>1e4)
// {
// printf("stuck in draw glitch amplitude\n");
// printf("SNR=%g\n",SNR);
// printf("SNRmin=%g\n",SNRmin);
// printf("SNRmax=%g\n",SNRmax);
// printf("Q=%g\n",params[2]);
// printf("f=%g\n",params[1]);
// printf("eta=%g\n",Sa);
// printf("Snf=%g\n",Snf[i]);
// exit(1);
// }
}
//SNR defined with Sn(f) but Snf array holdes <n_i^2>
......@@ -450,15 +441,19 @@ void fisher_matrix_proposal(struct FisherMatrix *fisher, double *params, long *s
// not used since we jump along eigendirections
double sqD = 1./sqrt((double)(N));
//choose a single eigenvector to jump over
i = (int)(ran2(seed)*(double)N);
//draw the eigen-jump amplitudes from N[0,1] scaled by evalue & dimension
for(i=0; i<N; i++)
{
// for(i=0; i<N; i++)
// {
Amps[i] = 1.0/sqrt(fisher->evalue[i]);
jump[i] = 0.0;
}
// }
//decompose eigenjumps into parameter directions
for(i=0; i<N; i++) for (j=0; j<N; j++) jump[j] += Amps[i]*fisher->evector[j][i]*sqD;
//for(i=0; i<N; i++)
for (j=0; j<N; j++) jump[j] += Amps[i]*fisher->evector[j][i]*sqD;
//jump from current position
for(i=0; i<N; i++) y[i] = params[i] + gasdev2(seed)*jump[i];
......@@ -509,6 +504,7 @@ void intrinsic_fisher_matrix(double *params, struct FisherMatrix *fisher)
double b2 = beta*beta;
double invf= 1./f;
double invQ= 1./Q;
double invA= 1./A;
// Zero matrix out
for (ii = 0; ii < Nparams; ii++) {
......@@ -544,7 +540,7 @@ void intrinsic_fisher_matrix(double *params, struct FisherMatrix *fisher)
analytic_fisher[2][1] = analytic_fisher[1][2];
// f lnA
analytic_fisher[1][3] = -0.5*invf;//-(1.0/(2.0*f));
analytic_fisher[1][3] = -0.5*invf*invA;//-(1.0/(2.0*f));
analytic_fisher[3][1] = analytic_fisher[1][3];
// f phi
......@@ -559,7 +555,7 @@ void intrinsic_fisher_matrix(double *params, struct FisherMatrix *fisher)
analytic_fisher[2][2] = 0.25*(3.0*(1.0+PI2*b2))*invQ*invQ;//(3.0*(1.0+LAL_PI*LAL_PI*beta*beta)/(4.0*Q*Q));
// Q lnA
analytic_fisher[2][3] = 0.5*invQ;//(1.0/(2.0*Q));
analytic_fisher[2][3] = 0.5*invQ*invA;//(1.0/(2.0*Q));
analytic_fisher[3][2] = analytic_fisher[2][3];
// Q phi
......@@ -571,7 +567,7 @@ void intrinsic_fisher_matrix(double *params, struct FisherMatrix *fisher)
//analytic_fisher[5][2] = analytic_fisher[2][5];
// lnA lnA
analytic_fisher[3][3] = 1.0;
analytic_fisher[3][3] = invA*invA;
// phi phi
analytic_fisher[4][4] = 1.0;
......@@ -585,7 +581,7 @@ void intrinsic_fisher_matrix(double *params, struct FisherMatrix *fisher)
}
void intrinsic_fisher_proposal(double **range, double *paramsx, double *paramsy, double *Snf, double Snx, double *logpx, double *logpy, double scale, long *seed, double Tobs, int *test)
void intrinsic_fisher_proposal(UNUSED struct Model *model, double **range, double *paramsx, double *paramsy, double *Snf, double Snx, double *logpx, double *logpy, long *seed, double Tobs, int *test)
{
int k;
double x,y;
......@@ -593,6 +589,7 @@ void intrinsic_fisher_proposal(double **range, double *paramsx, double *paramsy,
double qxy = *logpx;
double *dparamsx = dvector(0,5);
double *dparamsy = dvector(0,5);
double scale = 1./sqrt(5.);//Number of intrinsic parameters
// this is a crude diagonal Fisher matrix for individual sine-gaussians
intrinsic_fisher_update(paramsx,dparamsx,Snf,Snx,Tobs);
......@@ -649,9 +646,9 @@ void intrinsic_fisher_proposal(double **range, double *paramsx, double *paramsy,
free_dvector(dparamsy,0,5);
}
void intrinsic_fisher_proposal_2(struct Model *model, double **range, double *paramsx, double *paramsy, double *logpx, double *logpy, long *seed, int *test)
void intrinsic_fisher_proposal_2(struct Model *model, double **range, double *paramsx, double *paramsy, double *Snf, double Snx, double *logpx, double *logpy, long *seed, double Tobs, int *test)
{
int k;
int j,k;
double qyx = *logpy;
double qxy = *logpx;
double *dparams = dvector(0,model->intfisher->N);
......@@ -660,9 +657,21 @@ void intrinsic_fisher_proposal_2(struct Model *model, double **range, double *pa
struct FisherMatrix *fisher = model->intfisher;
//estimate SNR of current wavelet
double SNR;//,SNR2;
k = (int)(paramsx[1]*Tobs);
//SNR defined with Sn(f) but Snf array holdes <n_i^2>
SNR = paramsx[3]*sqrt((paramsx[2]/(2.0*RT2PI*paramsx[1]))/(Snx*Snf[k]*2.0/Tobs));
if(SNR < 5.0) SNR = 5.0; // this caps the size of the proposed jumps
//calculate fisher for current location
intrinsic_fisher_matrix(paramsx, fisher);
//rescale fisher matrix by 1/SNR^2
for(j=0; j<fisher->N; j++)for(k=0; k<fisher->N; k++) fisher->matrix[j][k]*=SNR*SNR;
//get eigenvectors of current fisher
matrix_eigenstuff(fisher->matrix, fisher->evector, fisher->evalue, fisher->N);
detJ = matrix_jacobian(fisher->matrix, fisher->N);
......@@ -687,11 +696,21 @@ void intrinsic_fisher_proposal_2(struct Model *model, double **range, double *pa
//deltax * Gamma * deltax
matrix_multiply(fisher->matrix,dparams,xdotF,fisher->N);
qyx += log(-0.5*dot_product(dparams,xdotF,fisher->N)) - 0.5*log(detJ);
qyx += -0.5*dot_product(dparams,xdotF,fisher->N) - 0.5*log(detJ);
//need fisher at new location to compute reverse probability
intrinsic_fisher_matrix(paramsy, fisher);
//SNR defined with Sn(f) but Snf array holdes <n_i^2>
k = (int)(paramsy[1]*Tobs);
SNR = paramsy[3]*sqrt((paramsy[2]/(2.0*RT2PI*paramsy[1]))/(Snx*Snf[k]*2.0/Tobs));
if(SNR < 5.0) SNR = 5.0; // this caps the size of the proposed jumps
//rescale fisher matrix by 1/SNR^2
for(j=0; j<fisher->N; j++)for(k=0; k<fisher->N; k++) fisher->matrix[j][k]*=SNR*SNR;
//get eigenvectors of current fisher
matrix_eigenstuff(fisher->matrix, fisher->evector, fisher->evalue, fisher->N);
detJ = matrix_jacobian(fisher->matrix, fisher->N);
......@@ -699,12 +718,13 @@ void intrinsic_fisher_proposal_2(struct Model *model, double **range, double *pa
//compute probabiltiy of reverse proposal
//deltax * Gamma * deltax
matrix_multiply(fisher->matrix,dparams,xdotF,fisher->N);
qxy += log(-0.5*dot_product(dparams,xdotF,fisher->N)) - 0.5*log(detJ);
qxy += -0.5*dot_product(dparams,xdotF,fisher->N) - 0.5*log(detJ);
// map phase onto [0,2pi]
while(paramsy[4] > LAL_TWOPI) paramsy[4] -=LAL_TWOPI;
while(paramsy[4] < 0.0) paramsy[4] +=LAL_TWOPI;
//update incoming pointers and return
*logpy += qyx;
*logpx += qxy;
......
#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif
......@@ -29,9 +34,9 @@ void extrinsic_fisher_update(struct Data *data, struct Model *model);
void intrinsic_fisher_update(double *params, double *dparams, double *Snf, double Snx, double Tobs);
void intrinsic_fisher_proposal(double **range, double *paramsx, double *paramsy, double *Snf, double Snx, double *logpx, double *logpy, double scale, long *seed, double Tobs, int *test);
void intrinsic_fisher_proposal(UNUSED struct Model *model, double **range, double *paramsx, double *paramsy, double *Snf, double Snx, double *logpx, double *logpy, long *seed, double Tobs, int *test);
void intrinsic_fisher_proposal_2(struct Model *model, double **range, double *paramsx, double *paramsy, double *logpx, double *logpy, long *seed, int *test);
void intrinsic_fisher_proposal_2(struct Model *model, double **range, double *paramsx, double *paramsy, double *Snf, double Snx, double *logpx, double *logpy, long *seed, double Tobs, int *test);
/* ********************************************************************************** */
/* */
......
......@@ -18,8 +18,13 @@
LIBS = lalsupport lalframe Frame lalinference lalinspiral lalpulsar lalmetaio lalsimulation metaio lal gsl fftw3 fftw3f gslcblas m
CCFLAGS = -g -Wall -O3 -std=gnu99
#CCFLAGS = -g -Wall
#CCFLAGS = -g -Wall -O3 -std=gnu99
#DEBUGGING OPTIONS
CCFLAGS = -g -Wall -O3 -march=native
#OPTIMIZATION OPTIONS
#CCFLAGS = -Wall -O3 -march=native
CC = gcc
......
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