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

new orientation proposal

rearranged extrinsic parameter arrays
continued cleaning code, making memory and constants more consistent



git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@549 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent 5bc42499
......@@ -89,7 +89,7 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, ra=None, dec=None, results=
# num, L, ralist, sin_dec, psi, e, dA, dphi, dt = np.loadtxt(filename, unpack=True)
print "Extracting RA/DEC samples"
filename = './chains/' + jobname + 'signal_params.dat.0'
data = np.loadtxt(filename, unpack=True)
data = np.loadtxt(filename, unpack=True,usecols=(0,1,2))
ralist = data[1]
sin_dec = data[2]
print "Total samples are {0}".format(ralist.size)
......
This diff is collapsed.
......@@ -45,6 +45,26 @@
#include <lal/LALInference.h>
#include <lal/LALInferenceReadData.h>
/* ********************************************************************************** */
/* */
/* Globally defined parameters */
/* */
/* ********************************************************************************** */
#define RTPI 1.772453850905516 // sqrt(Pi)
#define RT2PI 2.5066282746310005024 // sqrt(2 Pi)
#define CLIGHT 2.99792458e8 // Speed of light (m/s)
#define NE 6 //Number of extrinsic parameters
#define NW 5 //Number of intrinsic parameters
//Tuning parameters for wavelet proximity proposal
//TODO: These should not be global!
#define UNIFORM 0.2 // fraction of proximity prior that is uniform
#define SPRD 4 // width of outer part of proximity proposal
#define HLW 1 // width of hollowed-out region
/* ********************************************************************************** */
/* */
/* Data structures */
......
//Globally defined parameters.
#define TSUN 4.92569043916e-6 // mass to seconds conversion
#define PI 3.1415926535897931159979634685442 // Pi
#define PI2 9.869604401089357992304940125905 // Pi^2
#define TPI 6.2831853071795862319959269370884 // 2 Pi
#define RTPI 1.772453850905516 // sqrt(Pi)
#define RT2PI 2.5066282746310005024 // sqrt(2 Pi)
#define EULER 0.577215664901532860606512 // Euler's constant
#define E 2.71828183 // e
#define MPC 3.08568025e22 // MegaParsec in meters
#define CLIGHT 2.99792458e8 // Speed of light in m/s
#define AF0 0.1063937550698317
#define NS 100 // maximum number of sine-gaussian
#define NE 7
#define NW 5
#define NP 57 // number of parameters = NE + NS*NG
#define LIGO3 0 // use LIGO 3 noise curve when set equal to 1
#define skyfix 0 // can fix sky location
#define anglefix 0 // fix sky and orientation
//Tuning parameters for wavelet proximity proposal
#define UNIFORM 0.2 // fraction of proximity prior that is uniform
#define SPRD 4 // width of outer part of proximity proposal
#define HLW 1 // width of hollowed-out region
#define APN 3.4774624728292402883 // normalization factor in signal amplitude prior
//Runge-Kutta parameters
#define SAFETY 0.9
#define PGROW -0.2
#define PSHRNK -0.25
#define ERRCON 1.89e-4
#define MAXSTP 1000000
//#define BURNIN 10000 // number of samples in the burnin phase
#define derate 0.05 // fraction of the time we add a point to the diff_ev array
......@@ -10,7 +10,6 @@
#include "BayesWaveWavelet.h"
#include "BayesWaveProposal.h"
#include "BayesWaveEvidence.h"
#include "BayesWaveConstants.h"
#include "BayesWaveLikelihood.h"
......
This diff is collapsed.
......@@ -8,7 +8,6 @@
#include "BayesWaveModel.h"
#include "BayesWaveWavelet.h"
#include "BayesWaveLikelihood.h"
#include "BayesWaveConstants.h"
/* ********************************************************************************** */
/* */
......@@ -124,25 +123,14 @@ double loglike_maxNET(struct Data *data, struct Prior *prior, double **sNET, dou
logL = -0.5*(SSNET - HSNET*HSNET/rhoNET);
//Amplitude
//params[3] *= norm;
*dA = norm;
//Time
//params[0] += delt;
*dt = delt;
//Phase
//params[4] -= pshift;
//*dphi = -pshift;
*dphi = pshift;
/*
if(params[4] < 0.0) params[4] += TPI;
if(params[4] > TPI) params[4] -= TPI;
if(params[0] < 0.0) params[0] += Tobs;
if(params[0] > Tobs) params[0] -= Tobs;
*/
for(ifo=0; ifo<NI; ifo++) logL -= (double)(imax-imin)*log(SnNET[ifo]);
if(logL!=logL)
......@@ -575,10 +563,10 @@ double EvaluateSearchLogLikelihood(int typ, int ii, int det, struct Model *mx, s
// fclose(prop);
//abort();
if(paramsy[4] < 0.0) paramsy[4] += TPI;
if(paramsy[4] > TPI) paramsy[4] -= TPI;
if(paramsy[4] < 0.0) paramsy[4] += LAL_TWOPI;
if(paramsy[4] > LAL_TWOPI) paramsy[4] -= LAL_TWOPI;
if(paramsy[0] < 0.0) paramsy[0] += Tobs;
if(paramsy[0] < 0.0) paramsy[0] += Tobs;
if(paramsy[0] > Tobs) paramsy[0] -= Tobs;
// check that we haven't mapped out of range
......@@ -861,7 +849,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
{
wBirth->intParams[jj][0] += my->projection->deltaT[ifo];
wBirth->intParams[jj][3] *= sqrt(my->projection->Fplus[ifo]*my->projection->Fplus[ifo]+my->projection->Fcross[ifo]*my->projection->Fcross[ifo]);
wBirth->intParams[jj][4] = ran2(&chain->seed)*TPI;
wBirth->intParams[jj][4] = ran2(&chain->seed)*LAL_TWOPI;
}
else wBirth->intParams[jj][0] -= my->projection->deltaT[ifo];
......@@ -1206,10 +1194,10 @@ double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *
//Phase
params[6] += dphi;
if(params[6] < 0.0) params[6] += TPI;
if(params[6] > TPI) params[6] -= TPI;
if(params[6] < 0.0) params[6] += LAL_TWOPI;
if(params[6] > LAL_TWOPI) params[6] -= LAL_TWOPI;
if(params[5] < 0.0) params[5] += Tobs;
if(params[5] < 0.0) params[5] += Tobs;
if(params[5] > Tobs) params[5] -= Tobs;
// check that we haven't mapped out of range
......
This diff is collapsed.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "BayesWave.h"
#define IM1 2147483563
#define IM2 2147483399
......@@ -28,9 +29,8 @@ double swap, tempr;
double gaussian_norm(double min, double max, double sigma)
{
double sq2 = 1.4142135623731; //sqrt(2)
double sqp2 = 1.2533141373155; //sqrt(pi/2)
return sqp2*sigma*(erfc(min/sq2/sigma) - erfc(max/sq2/sigma) );
return sqp2*sigma*(erfc(min/LAL_SQRT2/sigma) - erfc(max/LAL_SQRT2/sigma) );
}
void recursive_phase_evolution(double dre, double dim, double *cosPhase, double *sinPhase)
......
......@@ -7,7 +7,6 @@
#include "BayesWavePrior.h"
#include "BayesWaveWavelet.h"
#include "BayesWaveProposal.h"
#include "BayesWaveConstants.h"
#include "BayesWaveLikelihood.h"
#define NR_END 1
......@@ -258,7 +257,7 @@ void initialize_priors(struct Data *data, struct Prior *prior, int omax)
prior->range[3][0] = Amin;
prior->range[3][1] = Amax;
prior->range[4][0] = 0.0;
prior->range[4][1] = TPI;
prior->range[4][1] = LAL_TWOPI;
prior->TFV = (tmax-tmin)*(fmax-fmin);
prior->logTFV = log(prior->TFV);
......@@ -445,7 +444,7 @@ void reset_priors(struct Data *data, struct Prior *prior)
prior->range[3][0] = Amin;
prior->range[3][1] = Amax;
prior->range[4][0] = 0.0;
prior->range[4][1] = TPI;
prior->range[4][1] = LAL_TWOPI;
prior->TFV = (tmax-tmin)*(fmax-fmin);
prior->logTFV = log(prior->TFV);
......@@ -538,13 +537,11 @@ void initialize_model(struct Model *model, struct Data *data, struct Prior *prio
model->extParams[6] = 0.0; //geocenter phase shift
*/
//draw extrinsic parameters from prior
model->extParams[5] = 0.0; //geocenter time shift
extrinsic_uniform_proposal(seed,model->extParams);
//set "delta" parameters to not modify signal model
model->extParams[4] = 1.0; //amplitude scale
model->extParams[5] = 0.0; //geocenter time shift
model->extParams[6] = 0.0; //geocenter phase shift
model->extParams[4] = 0.0; //geocenter phase shift
model->extParams[5] = 1.0; //amplitude scale
//allocate memory for signal wavelet structure
......
......@@ -26,7 +26,6 @@
#include "BayesWaveWavelet.h"
#include "BayesWaveProposal.h"
#include "BayesWaveEvidence.h"
#include "BayesWaveConstants.h"
#include "BayesWaveLikelihood.h"
/************* PROTOTYPE DECLARATIONS FOR INTERNAL FUNCTIONS **************/
......@@ -2637,7 +2636,7 @@ void get_time_domain_hdot(double *hdot, double *h, int N, double *Snf, int imin,
int i;
double xx, yy;
for(i=0; i<N; i++)hdot[i]= -TPI*h[i]; // Take care of the -2pi part
for(i=0; i<N; i++)hdot[i]= -LAL_TWOPI*h[i]; // Take care of the -2pi part
hdot[0] = 0.0;
hdot[1] = 0.0;
......
......@@ -10,8 +10,6 @@
#include "BayesWavePrior.h"
#include "BayesWaveWavelet.h"
#include "BayesWaveConstants.h"
//data->channels[ifo]
void set_bayesline_priors(char *channel, struct BayesLineParams *bayesline)
{
......@@ -86,27 +84,28 @@ int checkrange(double *p, double **r)
int extrinsic_checkrange(double *p)
{
if(p[0] < 0.0) p[0] += TPI;
if(p[0] >= TPI) p[0] -= TPI;
if(p[0] < 0.0) p[0] += LAL_TWOPI;
if(p[0] >= LAL_TWOPI) p[0] -= LAL_TWOPI;
if(p[1] < -1.0 || p[1] >= 1.0) return 1;
while(p[2] < 0.0)
p[2] += 0.5*LAL_PI;
while(p[2] >= 0.5*LAL_PI)
p[2] -= 0.5*LAL_PI;
p[2] += LAL_PI_2; //PI/2
while(p[2] >= LAL_PI_2)
p[2] -= LAL_PI_2;
if(p[3] < -0.99 || p[3] >= 0.99) return 1;
if(p[4] < 0.0) p[4] += LAL_TWOPI;
if(p[4] >=LAL_TWOPI) p[4] -= LAL_TWOPI;
//Each wavelet's amplitude is checked against intrinsic prior during extrinsic updates
if(p[4] < 0) return 1;
if(p[5] < 0) return 1;
//Each wavelet's time is checked against intrinsic prior duirng extrinsic updates
//if(p[5] < -0.1 || p[5] >=0.1) return 1;//{printf("***********5 %g,%g\n",p[5],Tobs);return 1;}//return 1;
if(p[6] < 0.0) p[6] += TPI;
if(p[6] >=TPI) p[6] -= TPI;
return 0;
}
......@@ -149,10 +148,10 @@ static double proximity_prior_density(double t, double f, double **params, int *
df = fabs(f-params[kk][1]);
dt = fabs(t-params[kk][0]);
tau = params[kk][2]/(TPI*params[kk][1]);
tau = params[kk][2]/(LAL_TWOPI*params[kk][1]);
sigwt = SPRD*tau;
sigwf = SPRD/(PI*tau);
sigwf = SPRD/(LAL_PI*tau);
if(dt < 6.0*sigwt && df < 6.0*sigwf) // don't bother adding in contribution from wavelets that are too far away
......@@ -162,7 +161,7 @@ static double proximity_prior_density(double t, double f, double **params, int *
sigwt2 = sigwt*sigwt;
signt2 = signt*signt;
signf = HLW/(PI*tau);
signf = HLW/(LAL_PI*tau);
sigwf2 = sigwf*sigwf;
signf2 = signf*signf;
......
This diff is collapsed.
......@@ -59,6 +59,10 @@ double TFprop(double f, double t, double **prop, int tsize, int fsize, double To
/* */
/* ********************************************************************************** */
void sky_ring_proposal(double *x, long *seed, double *y, struct Data *data, double *snr, double *logq, int ic);
void sky_ring_proposal(double *x, double *y, struct Data *data, long *seed, double *logq, int ic);
void uniform_orientation_proposal(double *x, double *y, struct Data *data, long *seed);
void network_orientation_proposal(double *paramsx, double *paramsy, struct Data *data, double *logJ);
void uniform_orientation_proposal(double *x, long *seed, double *y, struct Data *data);
......@@ -5,7 +5,6 @@
#include "BayesWave.h"
#include "BayesWaveMath.h"
#include "BayesWaveWavelet.h"
#include "BayesWaveConstants.h"
/* ********************************************************************************** */
......@@ -28,7 +27,7 @@ void SineGaussianBandwidth(double *sigpar, double *fmin, double *fmax)
double f0 = sigpar[1];
double Q = sigpar[2];
double tau = Q/(TPI*f0);
double tau = Q/(LAL_TWOPI*f0);
*fmax = f0 + 3.0/(tau); // no point evaluating waveform past this time (many efolds down)
*fmin = f0 - 3.0/(tau); // no point evaluating waveform before this time (many efolds down)
......@@ -54,7 +53,7 @@ void SineGaussianFourier(double *hs, double *sigpar, int N, int flag, double Tob
SineGaussianBandwidth(sigpar,&fmin,&fmax);
tau = Q/(TPI*f0);
tau = Q/(LAL_TWOPI*f0);
//Make sure frequency is in range
if(fmin < 0.0)
......@@ -81,19 +80,19 @@ void SineGaussianFourier(double *hs, double *sigpar, int N, int flag, double Tob
/* Use recursion relationship for cos(Phase) & sin(Phase) */
//initial values of exp(iPhase)
double phase = TPI*fmin*t0;
double phase = LAL_TWOPI*fmin*t0;
double cosPhase_m = cos(phase-phi);
double sinPhase_m = sin(phase-phi);
double cosPhase_p = cos(phase+phi);
double sinPhase_p = sin(phase+phi);
//incremental values of exp(iPhase)
double dim = -sin(TPI*t0*invTobs);
double dre = sin(0.5*(TPI*t0*invTobs));
double dim = -sin(LAL_TWOPI*t0*invTobs);
double dre = sin(0.5*(LAL_TWOPI*t0*invTobs));
dre = -2.0*dre*dre;
double amplitude = 0.5*(Amp)*RTPI*tau;
double pi2tau2 = PI*PI*tau*tau;
double pi2tau2 = LAL_PI*LAL_PI*tau*tau;
double Q2 = Q*Q/f0;
for(i = istart; i < imin; i++)
......@@ -175,13 +174,12 @@ void computeProjectionCoeffs(struct Data *data, struct Network *projection, doub
double Fplus, Fcross;
// common extrinsic parameters
//[0] alpha, [1] sin(delta) [2] psi [3] ellipticity [4] scale [5] dt [6] phi0
//[0] alpha, [1] sin(delta) [2] psi [3] ellipticity [4] phi0 [5] scale
double alpha = extParams[0];
double sindelta = extParams[1];
double psi = extParams[2];
double dt = extParams[5];
double phi0 = extParams[6];
double phi0 = extParams[4];
double cosPhase;
double sinPhase;
......@@ -199,8 +197,8 @@ void computeProjectionCoeffs(struct Data *data, struct Network *projection, doub
XLALComputeDetAMResponse(&Fplus, &Fcross, (const REAL4(*)[3])data->detector[ifo]->response, alpha, asin(sindelta), psi, XLALGreenwichMeanSiderealTime(&data->epoch));
//initialize times
toff = dt + dtimes[ifo]; //offset time for detector ifo
toff -= dtimes[0]; //reference time is TOA at IFO0
toff = dtimes[ifo]; //offset time for detector ifo
toff -= dtimes[0]; //reference time is TOA at IFO0
//store the timedelay to each IFO for use elsewhere
projection->deltaT[ifo] = toff;
......@@ -251,7 +249,7 @@ void waveformProject(struct Data *data, struct Network *projection, double *extP
int i, id, even, odd;
double ecc = extParams[3];
double scale = extParams[4];
double scale = extParams[5];
int n = data->N;
int iend = n/2;
......@@ -309,7 +307,7 @@ void compute_reconstructed_plus_and_cross(struct Data *data, struct Network *pro
int i, id;
double ecc = extParams[3];
double scale = extParams[4];
double scale = extParams[5];
int n = data->N;
double Tobs = data->Tobs;
......
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