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

split Subroutines.c int smaller source files

changed executable names to match Unix style guide



git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@539 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent a953c25e
"""
This script generates a webpage to display the results of a BayesWaveBurst run.
This script generates a webpage to display the results of a BayesWave run.
This script was tested on ldas-pcdev1.ligo.caltech.edu only.
"""
......@@ -99,43 +99,6 @@ htmlDir = 'html/'
# Read in run data and info
#
######################################################################################################################
# ---------------------------------------------------------------------
# Determine name of job and check for special flags and for mdc status
# ---------------------------------------------------------------------
def get_jobname():
# -- Determine name of job
bayeswaverunfile = glob.glob('*bayeswave.run')
if not len(bayeswaverunfile):
sys.exit("\n bayeswave.run not found! \n")
bayeswaverunfile = bayeswaverunfile[0]
bayeswave = open(bayeswaverunfile, 'r')
jobName = bayeswaverunfile.split('_')
del jobName[-1]
if not len(jobName) == 0:
jobName = '_'.join(jobName)+'_'
else:
jobName = ''
# -- Determine if the job was glitchOnly, noiseOnly, or signalOnly
restrictModel = ''
mdc = False
injFlag = False
for line in bayeswave:
if (line.find('glitchOnly')>-1):
print '\nThis run was executed with the --glitchOnly flag\n'
restrictModel = 'glitch'
elif (line.find('noiseOnly')>-1):
print '\nThis run was executed with the --noiseOnly flag\n'
restrictModel = 'noise'
elif (line.find('signalOnly')>-1):
print '\nThis run was executed with the --signalOnly flag\n'
restrictModel = 'signal'
if (line.find('--inj')>-1) or (line.find('MDC')>-1):
mdc = True
if line.find('--inj')>-1 and line.find('--lalinspiralinjection')<0:
injFlag = True
# -- Report to the user
print "The job name is: {0}".format(jobName)
return(jobName, restrictModel, mdc, injFlag)
# ----------------------------------
# Extract information about the run
......@@ -157,8 +120,11 @@ def readbwb():
sys.exit("\n *bayeswave.run not found! \n")
bayeswaverunfile = bayeswaverunfile[0]
bayeswave = open(bayeswaverunfile, 'r')
# -- Parse command line arguments
cmdline = bayeswave.readline().split(' ')
lines = bayeswave.readlines()
cmdline = lines[7].split(' ')
for index, arg in enumerate(cmdline):
# Determine the IFOs involved
if arg=='--ifo':
......@@ -319,7 +285,7 @@ def snrfunctime(median_waveform):
def get_axes(jobName, postDir, ifoList, worc, model, time):
axisList = []
for ifo in ifoList[0:-1]:
print ifo
# -- Read Signal model
try:
filename = str(jobName)+postDir+'{1}_recovered_{2}_waveform.dat.{0}'.format(ifo, model, worc)
......@@ -1017,7 +983,8 @@ def make_mode_table(page, subpage, imdc, momentsPlus, momentsCross, ifoList, tab
subpage.write(' </tr>\n')
# -- Loop over whole signal, + and x polarizations
for polarization in ['', '_plus', '_cross']:
for ifo in ifoList: # Single detector
shortlist = ifoList[:2]
for ifo in shortlist: # Single detector
plotsrc = './'+plotsDir+html_dict[page][2]+polarization+'_{1}_{0}.png'.format(ifo,'lin')
if os.path.exists(plotsrc): # Indirect way of making sure this overlap was calculated
mode, low, high, inmode = mode_values('overlap'+polarization, ifo, mdc, momentsPlus, momentsCross, ifoList, tablesDir, worc)
......@@ -1031,7 +998,7 @@ def make_mode_table(page, subpage, imdc, momentsPlus, momentsCross, ifoList, tab
for value in [mode, low, high]:
subpage.write(' <td>{0}</td>\n'.format(value))
subpage.write(' </tr>\n')
if len(ifoList)>1: # Network overlap
if len(shortlist)>1: # Network overlap
mode, low, high, inmode = mode_values('network_overlap'+polarization, '0', mdc, momentsPlus, momentsCross, ifoList, tablesDir, worc)
subpage.write(' <tr>\n')
if polarization == '':
......@@ -1385,20 +1352,20 @@ def make_index(htmlDir, plotsDir, tablesDir, model, gps, ifoList, ifoNames, snrL
def make_webpage(htmlDir, model, mdc, gps, ifoList, ifoNames, modelList, snrList, sig_gl, sig_noise, momentsPlus, momentsCross):
# TODO: move in separate function
# -- Find out the path to the BayesWaveBurst executable
p = subprocess.Popen(["which","BayesWaveBurst"],stdout=subprocess.PIPE,stderr=subprocess.STDOUT)
# -- Find out the path to the BayesWave executable
p = subprocess.Popen(["which","bayeswave"],stdout=subprocess.PIPE,stderr=subprocess.STDOUT)
for line in iter(p.stdout.readline, b''):
BayesWaveBurst = line.rstrip()
if not os.path.isfile(BayesWaveBurst):
sys.exit("\nMake sure the BayesWaveBurst command is in your path before running this script!\n")
bayeswave = line.rstrip()
if not os.path.isfile(bayeswave):
sys.exit("\nMake sure the bayeswave command is in your path before running this script!\n")
# -- Use the path to copy over the css and js files
BayesWaveBurst = BayesWaveBurst.split('src')
del BayesWaveBurst[-1]
BayesWaveBurst.append('postprocess/')
BayesWaveBurst = ''.join(BayesWaveBurst)
os.system('cp '+BayesWaveBurst+'BWBweb.css '+htmlDir+'.')
os.system('cp '+BayesWaveBurst+'secure_ajax.js '+htmlDir+'.')
os.system('cp '+BayesWaveBurst+'navigate.js '+htmlDir+'.')
bayeswave = bayeswave.split('src')
del bayeswave[-1]
bayeswave.append('postprocess/')
bayeswave = ''.join(bayeswave)
os.system('cp '+bayeswave+'BWBweb.css '+htmlDir+'.')
os.system('cp '+bayeswave+'secure_ajax.js '+htmlDir+'.')
os.system('cp '+bayeswave+'navigate.js '+htmlDir+'.')
# -- Write the index
......@@ -1430,7 +1397,7 @@ def make_webpage(htmlDir, model, mdc, gps, ifoList, ifoNames, modelList, snrList
#
######################################################################################################################
# -- Parse command line arguments
parser = argparse.ArgumentParser(description='Produce html page for a BayesWaveBurst run.')
parser = argparse.ArgumentParser(description='Produce html page for a bayeswave run.')
# -- Get basic info on the run
jobName, restrictModel, mdc, injFlag, bayeswaverunfile, ifoList, ifoNames, gps, snrList, info = readbwb()
......
......@@ -26,7 +26,7 @@ import lalinference.cmap
# Works with single output directory from
# BayesWaveBurst
# -------------------------------------------
def make_skyview(directory='.', mdc=None, NSIDE=64, ra=None, dec=None, results=None, npost=5000):
def make_skyview(directory='.', mdc=None, NSIDE=128, ra=None, dec=None, results=None, npost=5000):
# -- Close any open figures
plt.close('all')
......@@ -226,7 +226,7 @@ if __name__ == "__main__":
# -- Set default argument values
directory='.'
mdc=None
NSIDE=64
NSIDE=128
ra=None
dec=None
geo=False
......
......@@ -300,9 +300,9 @@ void LorentzSplineFit(BayesLinePriors *priors, bayesline_opt opts, int steps, da
}
else
{
z = 2.0*kappa*gsl_ran_gaussian_pdf(y, lAwidth);
z = 2.0*kappa_BL*gsl_ran_gaussian_pdf(y, lAwidth);
}
logpx += log((1.0-kappa)/(lAmax-lAmin) + z);
logpx += log((1.0-kappa_BL)/(lAmax-lAmin) + z);
logpy += -log((lAmax-lAmin));
}
......@@ -342,7 +342,7 @@ void LorentzSplineFit(BayesLinePriors *priors, bayesline_opt opts, int steps, da
logpx = -log((double)(ncut)); // corresponds to uniform density - just as likely to kill from alines_y->nwhere
alpha = gsl_rng_uniform(r);
if(alpha < kappa)
if(alpha < kappa_BL)
{
y = fabs(gsl_ran_gaussian(r, lAwidth));
lines_y->A[i] = exp(baseav+y);
......@@ -358,9 +358,9 @@ void LorentzSplineFit(BayesLinePriors *priors, bayesline_opt opts, int steps, da
}
else
{
z = 2.0*kappa*gsl_ran_gaussian_pdf(y, lAwidth);
z = 2.0*kappa_BL*gsl_ran_gaussian_pdf(y, lAwidth);
}
logpy += log((1.0-kappa)/(lAmax-lAmin) + z);
logpy += log((1.0-kappa_BL)/(lAmax-lAmin) + z);
logpx += -log((lAmax-lAmin));
if(lines_y->A[i] > Amaxx) check = 1;
......@@ -407,7 +407,7 @@ void LorentzSplineFit(BayesLinePriors *priors, bayesline_opt opts, int steps, da
alpha = gsl_rng_uniform(r);
if(alpha < kappa)
if(alpha < kappa_BL)
{
y = fabs(gsl_ran_gaussian(r, lAwidth));
lines_y->A[ii] = exp(baseav+y);
......@@ -423,9 +423,9 @@ void LorentzSplineFit(BayesLinePriors *priors, bayesline_opt opts, int steps, da
}
else
{
z = 2.0*kappa*gsl_ran_gaussian_pdf(y, lAwidth);
z = 2.0*kappa_BL*gsl_ran_gaussian_pdf(y, lAwidth);
}
logpy += log((1.0-kappa)/(lAmax-lAmin) + z);
logpy += log((1.0-kappa_BL)/(lAmax-lAmin) + z);
y = log(lines_x->A[ii])-baseav;
......@@ -435,9 +435,9 @@ void LorentzSplineFit(BayesLinePriors *priors, bayesline_opt opts, int steps, da
}
else
{
z = 2.0*kappa*gsl_ran_gaussian_pdf(y, lAwidth);
z = 2.0*kappa_BL*gsl_ran_gaussian_pdf(y, lAwidth);
}
logpx += log((1.0-kappa)/(lAmax-lAmin) + z);
logpx += log((1.0-kappa_BL)/(lAmax-lAmin) + z);
}
else
{
......@@ -1715,9 +1715,9 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, int steps, in
}
else
{
z = 2.0*kappa*gsl_ran_gaussian_pdf(y, lAwidth);
z = 2.0*kappa_BL*gsl_ran_gaussian_pdf(y, lAwidth);
}
logpx += log((1.0-kappa)/(lAmax-lAmin) + z);
logpx += log((1.0-kappa_BL)/(lAmax-lAmin) + z);
logpy += -log((lAmax-lAmin));
}
......@@ -1757,7 +1757,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, int steps, in
logpx = -log((double)(ncut)); // corresponds to uniform density - just as likely to kill from anywhere
alpha = gsl_rng_uniform(r);
if(alpha < kappa)
if(alpha < kappa_BL)
{
y = fabs(gsl_ran_gaussian(r, lAwidth));
lines_y->A[i] = exp(baseav+y);
......@@ -1773,9 +1773,9 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, int steps, in
}
else
{
z = 2.0*kappa*gsl_ran_gaussian_pdf(y, lAwidth);
z = 2.0*kappa_BL*gsl_ran_gaussian_pdf(y, lAwidth);
}
logpy += log((1.0-kappa)/(lAmax-lAmin) + z);
logpy += log((1.0-kappa_BL)/(lAmax-lAmin) + z);
logpx += -log((lAmax-lAmin));
......@@ -1856,7 +1856,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, int steps, in
lines_y->Q[ii] = exp(lQmin+(lQmax-lQmin)*gsl_rng_uniform(r));
alpha = gsl_rng_uniform(r);
if(alpha < kappa)
if(alpha < kappa_BL)
{
y = fabs(gsl_ran_gaussian(r, lAwidth));
lines_y->A[ii] = exp(baseav+y);
......@@ -1872,9 +1872,9 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, int steps, in
}
else
{
z = 2.0*kappa*gsl_ran_gaussian_pdf(y, lAwidth);
z = 2.0*kappa_BL*gsl_ran_gaussian_pdf(y, lAwidth);
}
logpy += log((1.0-kappa)/(lAmax-lAmin) + z);
logpy += log((1.0-kappa_BL)/(lAmax-lAmin) + z);
logpx += -log((lAmax-lAmin));
......@@ -2479,3 +2479,5 @@ void BayesLineRJMCMC(struct BayesLineParams *bayesline, double *freqData, double
invpsd[i] = 1./psd[i];
}
}
......@@ -19,7 +19,7 @@
//#define LAmax 1.0e-30
#define lAwidth 2.3 // 2.3 is one decade
#define zeta 1.0
#define kappa 0.8
#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
// as we are getting down to the line width of the broadest lines
......
This diff is collapsed.
/* ********************************************************************************** */
/* */
/* LAL Dependencies */
/* */
/* ********************************************************************************** */
#include <lal/LALInspiral.h>
#include <lal/LALCache.h>
#include <lal/LALFrStream.h>
#include <lal/TimeFreqFFT.h>
#include <lal/LALDetectors.h>
#include <lal/ResampleTimeSeries.h>
#include <lal/TimeSeries.h>
#include <lal/FrequencySeries.h>
#include <lal/Units.h>
#include <lal/Date.h>
#include <lal/StringInput.h>
#include <lal/VectorOps.h>
#include <lal/Random.h>
#include <lal/LALNoiseModels.h>
#include <lal/XLALError.h>
#include <lal/GenerateInspiral.h>
#include <lal/LIGOLwXMLRead.h>
#include <lal/LIGOLwXMLInspiralRead.h>
#include <lal/LALConstants.h>
#include <lal/SeqFactories.h>
#include <lal/DetectorSite.h>
#include <lal/GenerateInspiral.h>
#include <lal/GeneratePPNInspiral.h>
#include <lal/SimulateCoherentGW.h>
#include <lal/Inject.h>
#include <lal/LIGOMetadataTables.h>
#include <lal/LIGOMetadataUtils.h>
#include <lal/LIGOMetadataInspiralUtils.h>
#include <lal/LIGOMetadataRingdownUtils.h>
#include <lal/LALInspiralBank.h>
#include <lal/FindChirp.h>
#include <lal/LALInspiralBank.h>
#include <lal/GenerateInspiral.h>
#include <lal/NRWaveInject.h>
#include <lal/GenerateInspRing.h>
#include <lal/LALError.h>
#include <lal/LALInference.h>
#include <lal/LALInferenceReadData.h>
/* ********************************************************************************** */
/* */
/* Data structures */
/* */
/* ********************************************************************************** */
struct Wavelet
{
int size;
int smax;
int *index;
//clustering prior
int Ncluster;
double cosy;
double logp;
double *templates;
double **intParams;
double **TFhistory;
double **TFdensity;
double TFdensityMax;
char *name;
};
struct Model
{
int size;
double *Sn;
double logL;
double *detLogL;
double *extParams;
double **response;
struct Wavelet *signal;
struct Wavelet **glitch;
struct Network *projection;
struct FisherMatrix *fisher;
};
struct Network
{
double **expPhase;
double *deltaT;
double *Fplus;
double *Fcross;
};
struct Chain
{
int NC; // Number of parallel chains
// int NCmin; // Minimum number of parallel chains
// int NCmax; // Maximum number of parallel chains
int burn; // Number of burn in samples
int mod0; // Counter keeping track of noise model
int mod1; // Counter keeping track of glitch model
int mod2; // Counter keeping track of signal model
int mod3; // Counter keeping track of glitch+signal model
int cycle; // # of model updates per BurstRJMCMC iteration
int count; // # of BurstRJMCMC iterations
int mcount; // # of MCMC attempts
int scount; // # of RJMCMC attempts
int xcount; // # of extrinsic attempts
int macc; // # of MCMC successess
int sacc; // # of RJMCMC successes
int xacc; // # of extrinsic successes
int burnFlag;// flag telling proposals if we are in burn-in phase
int *index; // keep track of which chain is at which temperature
int *ptacc;
int *ptprop;
double *A;
long seed;
double modelRate; //fraction of updates to signal model
double rjmcmcRate; //fraction of trans- vs. inter-dimensional moves
double intPropRate; //fraction of prior vs. fisher intrinsic moves
double *dT;
double beta;
double tempMin;
double tempStep;
double *temperature;
double *avgLogLikelihood;
double *varLogLikelihood;
FILE *runFile;
FILE **intChainFile;
FILE **extChainFile;
FILE **intWaveChainFile;
FILE ***intGlitchChainFile;
};
struct Prior
{
int smin;
int smax;
int *gmin;
int *gmax;
int omax;
double Snmin;
double Snmax;
double TFV;
double logTFV;
double *bias;
double **range;
//background prior vectors
char bkgName[1024];
double *bkgFreq;
double *bkgCnts;
int bkgBins;
//cluster prior parameters
double alpha;
double beta;
double gamma;
//amplitude prior parameters
double sSNRpeak;
//glitch prior parameters
double gSNRpeak;
//dimension prior parameters
double Dtau;
char path[100];
};
struct Data
{
int N;
int NI;
int tsize;
int fsize;
int imin;
int imax;
int Dmax;
int Dmin;
int runPhase;
int psdFitFlag;
int noiseSimFlag;
int bayesLineFlag;
int fullOnlyFlag;
int noiseOnlyFlag;
int glitchOnlyFlag;
int signalOnlyFlag;
int cleanOnlyFlag;
int skipCleanFlag;
int skipSignalFlag;
int glitchFlag;
int signalFlag;
int noiseFlag;
int fixIntrinsicFlag;
int fixExtrinsicFlag;
int clusterPriorFlag;
int backgroundPriorFlag;
int amplitudePriorFlag;
int amplitudeProposalFlag;
int signalAmplitudePriorFlag;
int extrinsicAmplitudeFlag;
int orientationProposalFlag;
int clusterProposalFlag;
int adaptTemperatureFlag;
int splineEvidenceFlag;
int constantLogLFlag;
int gnuplotFlag;
int verboseFlag;
int rjFlag;
double dt;
double df;
double fmin;
double fmax;
double Tobs;
double Twin;
double Qmin;
double Qmax;
double gmst;
double logLc;
double TwoDeltaToverN;
double TFtoProx;
double **r;
double **s;
double **rh;
double *h;
double **Snf;
double **invSnf;
double **SnS;
double *SnGeo;
void (*signal_amplitude_proposal)(double *, double *, double, long *, double, double **, double);
void (*glitch_amplitude_proposal)(double *, double *, double, long *, double, double **, double);
double (*signal_amplitude_prior) (double *, double *, double, double, double);
double (*glitch_amplitude_prior) (double *, double *, double, double, double);
double (*intrinsic_likelihood)(int, int, int, struct Model*, struct Model*, struct Prior*, struct Chain*, struct Data*);
double (*extrinsic_likelihood)(struct Network*, double *, double *, double *, double **, struct Data*, double, double);
char runName[100];
char outputDirectory[1000];
char **channels;
struct BayesLineParams **bayesline;
ProcessParamsTable *commandLine;
LALDetector **detector;
LIGOTimeGPS epoch;
};
struct FisherMatrix
{
double *sigma;
double *evalue;
double **evector;
double **matrix;
double *aamps;
double **count;
};
\ No newline at end of file
This diff is collapsed.
void average_log_likelihood_via_direct_downsampling(int M, double **logLchain, struct Chain *chain, double *ACL, int nPoints);
void average_log_likelihood_via_recursion_downsampling(int M, double **logLchain, struct Chain *chain);
void average_log_likelihood_via_ACF(int M, double **logLchain, struct Chain *chain);
void SplineThermodynamicIntegration(double *temperature, double *avgLogLikelihood, double *varLogLikelihood, int NC, double *logEvidence, double *varLogEvidence);
void CovarianceThermodynamicIntegration(double *temperature, double **loglikelihood, int M, int NC, char modelname[], double *logEvidence, double *varLogEvidence);
void ThermodynamicIntegration(double *temperature, double **loglikelihood, int M, int NC, char modelname[], double *logEvidence, double *varLogEvidence);
void TrapezoidIntegration(struct Chain *chain, double **logLchain, int M, int counts, char modelname[], double *logEvidence, double *varLogEvidence);
void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *chain, struct Prior *prior, double *logEvidence);
This diff is collapsed.
/* ********************************************************************************** */
/* */
/* Data handling and injection routines */
/* */
/* ********************************************************************************** */
REAL8TimeSeries *readTseries(CHAR *cachefile, CHAR *channel, LIGOTimeGPS start, REAL8 length);
void InjectFromMDC(ProcessParamsTable *commandLine, LALInferenceIFOData *IFOdata, double *SNR);
void BayesWaveBurstInject(ProcessParamsTable *commandLine, struct Data *data, struct Chain *chain, struct Prior *prior);
void getChannels(ProcessParamsTable *commandLine, char **channel);
int outputFrameFile(ProcessParamsTable *commandLine, struct Data *data, struct Model *model);
/* ********************************************************************************** */
/* */
/* Output Files */