Commit 231962d2 authored by Tyson Littenberg's avatar Tyson Littenberg Committed by James Clark
Browse files

License branch

parent b2340762
......@@ -662,7 +662,8 @@ def plot_tf_tracks(jobName, postDir, ifo, plotsDir, worc, mdc, model, axwinner,
elif model == 'signal':
colour = scolor
t = time - 2.0
#t = time - 2.0
t = time;
# -- Kludgy fix to keep megaplot from crashing when data is bad
if np.amax(high_50) > 1e5: # Something is wrong
......@@ -833,7 +834,6 @@ def bayesogram(filename, plotsDir, worc, time, twodrange, median_waveform, axwin
# --------------------------------------
def Q_scan(model, Q, t, f, ifo, axwinner, f_axwinner, climwinner=[1,1]):
data = np.genfromtxt(postDir+'/{0}_spectrogram_{1}_{2}.dat'.format(model, Q, ifoNames[int(ifo)]))
data = data[...,::-1] # gets printed out time reversed, need to fix it
#plt.clim(np.min(data)/np.sum(data),np.max(data)/np.sum(data))
......@@ -2178,7 +2178,8 @@ for mod in modelList:
filename = str(jobName)+postDir+'{1}_median_tf_{0}.dat'.format(ifoNames[int(ifo)],mod)
freqsamp, median_powerspec, high_50, low_50, high_90, low_90 = get_waveform(filename)
plot_tf_tracks(jobName, postDir, ifo, plotsDir, worc, mdc, mod, axwinner, f_axwinner, freqsamp, low_50, high_50, low_90, high_90, median_powerspec);
#plot_tf_tracks(jobName, postDir, ifo, plotsDir, worc, mdc, mod, axwinner, f_axwinner, freqsamp, low_50, high_50, low_90, high_90, median_powerspec);
plot_tf_tracks(jobName, postDir, ifo, plotsDir, worc, mdc, mod, axwinner, f_axwinner, time, low_50, high_50, low_90, high_90, median_powerspec);
#if injFlag:
......
......@@ -99,13 +99,17 @@ int main(int argc, char *argv[])
}
NI=ifo;
double **freqData = dmatrix(0,NI-1,0,N-1);
double **psd = dmatrix(0,NI-1,0,N/2-1);
double **freqData = double_matrix(NI-1,N-1);
double **psd = double_matrix(NI-1,N/2-1);
struct Data *data = malloc(sizeof(struct Data));
struct Chain *chain = malloc(sizeof(struct Chain));
struct Prior *prior = malloc(sizeof(struct Prior));
const gsl_rng_type *T = gsl_rng_default;
chain->seed = gsl_rng_alloc(T);
gsl_rng_env_setup();
parse_command_line(data, chain, prior, runState->commandLine);
if (data->chirpletFlag) {
......@@ -267,7 +271,7 @@ int main(int argc, char *argv[])
for(ic=0; ic<chain->NC; ic++)
{
model[ic] = malloc(sizeof(struct Model));
initialize_model(model[ic], data, prior, psd, &chain->seed);
initialize_model(model[ic], data, prior, psd, chain->seed);
}
......@@ -431,7 +435,8 @@ int main(int argc, char *argv[])
fscanf(fptr,"%i",&chain->mc);
fscanf(fptr,"%i",&chain->NC);
fscanf(fptr,"%i",&chain->zcount);
fscanf(fptr,"%li",&chain->seed);
gsl_rng_fread(fptr,chain->seed);
fclose(fptr);
/****************************************************/
......@@ -844,9 +849,6 @@ int main(int argc, char *argv[])
RJMCMC(data, model, bayesline, chain, prior, &data->logZglitch, &data->varZglitch);
/*if(data->glitchOnlyFlag==1 && strcmp(data->channels[0],"AdLIGO"))
outputFrameFile(runState->commandLine, data, model[chain->index[0]]);*/
//shut off glitch model for checkpointing
data->glitchModelFlag = 0;
}
......
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
/* ********************************************************************************** */
/* */
......@@ -171,7 +173,7 @@ struct Chain
double *A;
long seed;
gsl_rng *seed;
double modelRate; //fraction of updates to signal model
double rjmcmcRate; //fraction of trans- vs. inter-dimensional moves
......@@ -335,8 +337,8 @@ struct Data
double varZfull;
void (*signal_amplitude_proposal)(double *, double *, double, long *, double, double **, double);
void (*glitch_amplitude_proposal)(double *, double *, double, long *, double, double **, double);
void (*signal_amplitude_proposal)(double *, double *, double, gsl_rng *, double, double **, double);
void (*glitch_amplitude_proposal)(double *, double *, double, gsl_rng *, double, double **, double);
double (*signal_amplitude_prior) (double *, double *, double, double);
double (*glitch_amplitude_prior) (double *, double *, double, double);
......
......@@ -9,6 +9,7 @@
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <fftw3.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
......@@ -25,9 +26,6 @@
/************* PROTOTYPE DECLARATIONS FOR INTERNAL FUNCTIONS **************/
double swap, tempr;
#define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}
#define REQARG 12
struct Data
......@@ -65,8 +63,7 @@ static REAL8TimeSeries *readTseries(CHAR *cachefile, CHAR *channel, LIGOTimeGPS
static void output_frame(REAL8TimeSeries *timeData, REAL8TimeSeries *timeRes, REAL8TimeSeries *timeGlitch, CHAR *frameType, CHAR *ifo);
void dfour1(double data[], unsigned long nn, int isign);
void drealft(double data[], unsigned long n, int isign);
void fftw_wrapper(double *data, int N, int flag);
void tukey_scale(double *s1, double *s2, double alpha, int N);
void tukey(double *data, double alpha, int N);
......@@ -246,7 +243,7 @@ int main(int argc, char *argv[])
tukey(glitch_model, alpha, data->bw_N);
// FFT
drealft(glitch_model-1, data->bw_N, 1);
fftw_wrapper(glitch_model, data->bw_N, 1);
// FFT scaling
norm = 2.0/(double)data->bw_N;
......@@ -263,7 +260,7 @@ int main(int argc, char *argv[])
}
// Inverse FFT
drealft(glitch_model_resampled-1, Nup, -1);
fftw_wrapper(glitch_model_resampled, Nup, -1);
/***********************************************************
......@@ -312,73 +309,6 @@ int main(int argc, char *argv[])
output_frame(timeData, timeRes, timeGlitch, outframeType, data->ifo);
XLALDestroyREAL8TimeSeries(timeData);
// printf("whitening a 32 segment surrounding the glitch and producing a spectrogram\n");
// FILE *outfile;
// double Tprint = 32.0;
//// Used to check cleaning
// N32 = (int)(Tprint/dT);
// double *G32 = malloc(N32*sizeof(double));
// double *S32 = malloc(N32*sizeof(double));
// double *L = malloc(N*sizeof(double));
// const gsl_rng_type * P;
// gsl_rng * r;
// gsl_rng_env_setup();
// P = gsl_rng_default;
// r = gsl_rng_alloc (P);
//
// double f;
// j -= N32/2;
//
// if(j < 0)
// {
// j = 0;
// }
//
// if(j > N-N32)
// {
// j = N-N32;
// }
//
// for(i=0; i< N32; i++)
// {
// G32[i] = Lclean[i+j];
// }
//
// alpha = (2.0*t_rise/32.0);
//
// tukey(G32, alpha, N32);
// drealft(G32-1, N32, 1);
//
// double *SL32, *SLn32;
//
// SL32 = dvector(0,N32/2-1);
// SLn32 = dvector(0,N32/2-1);
//
// // Power spectra
// for(i=0; i< N32/2; i++) SL32[i] = 2.0*(G32[2*i]*G32[2*i]+G32[2*i+1]*G32[2*i+1]);
//
// // Form spectral model for whitening data (lines plus a smooth component)
// spectrum(SL32, SLn32, 1.0/Tprint, r, N32);
//
// printf("spectrum computed\n");
//
// outfile=fopen("spec32.dat","w");
// for(i=1;i< N32/2;i++)
// {
// f = (double)(i)/32.0;
// fprintf(outfile,"%e %e\n", f, SLn32[i]);
// }
// fclose(outfile);
//
// printf("computing Q transform\n");
// QTransfrom(20.0, G32, SLn32, 8.0, 1024.0, Tprint, N32, 4);
//
// printf("plotting Q transform\n");
// system ("gnuplot specto.gnu");
// system ("open Qspec.png");
//
fprintf(stdout,"\n");
return 0;
}
......@@ -435,98 +365,77 @@ void tukey(double *data, double alpha, int N)
/* */
/* ********************************************************************************** */
void dfour1(double data[], unsigned long nn, int isign)
void fftw_wrapper(double *data, int N, int flag)
{
unsigned long n,mmax,m,j,istep,i;
double wtemp,wr,wpr,wpi,wi,theta;
double tempr,tempi, swap;
n=nn << 1;
j=1;
for (i=1;i<n;i+=2) {
if (j > i) {
SWAP(data[j],data[i]);
SWAP(data[j+1],data[i+1]);
}
m=n >> 1;
while (m >= 2 && j > m) {
j -= m;
m >>= 1;
}
j += m;
}
mmax=2;
while (n > mmax) {
istep=mmax << 1;
theta=isign*(6.28318530717959/mmax);
wtemp=sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi=sin(theta);
wr=1.0;
wi=0.0;
for (m=1;m<mmax;m+=2) {
for (i=m;i<=n;i+=istep) {
j=i+mmax;
tempr=wr*data[j]-wi*data[j+1];
tempi=wr*data[j+1]+wi*data[j];
data[j]=data[i]-tempr;
data[j+1]=data[i+1]-tempi;
data[i] += tempr;
data[i+1] += tempi;
int n;
double *timeData = (double *)malloc(N*sizeof(double));
fftw_complex *freqData = (fftw_complex *) fftw_malloc(sizeof(fftw_complex)*N);
switch (flag)
{
//Reverse transform
case -1:
//fill freqData with contents of data
for(n=0; n<N/2; n++)
{
freqData[n][0] = data[2*n];
freqData[n][1] = data[2*n+1];
}
wr=(wtemp=wr)*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;
}
mmax=istep;
}
}
//get DC and Nyquist where FFTW wants them
freqData[N/2][0] = freqData[0][1];
freqData[N/2][1] = freqData[0][1] = 0.0;
//setup FFTW plan
fftw_plan reverse = fftw_plan_dft_c2r_1d(N, freqData, timeData, FFTW_MEASURE);
//The FFT
fftw_execute(reverse);
fftw_destroy_plan(reverse);
//Copy output back into data
for(n=0; n<N; n++) data[n] = timeData[n];
break;
//Forward transform
case 1:
//fill timeData with contents of data
for(n=0; n<N; n++) timeData[n] = data[n];
//setup FFTW plan
fftw_plan forward = fftw_plan_dft_r2c_1d(N, timeData, freqData, FFTW_MEASURE);
//The FFT
fftw_execute(forward);
fftw_destroy_plan(forward);
//Copy output back into data
for(n=0; n<N/2; n++)
{
data[2*n] = freqData[n][0];
data[2*n+1] = freqData[n][1];
}
//get DC and Nyquist where FFTW wants them
data[1] = freqData[N/2][1];
void drealft(double data[], unsigned long n, int isign)
{
void dfour1(double data[], unsigned long nn, int isign);
unsigned long i,i1,i2,i3,i4,np3;
double c1=0.5,c2,h1r,h1i,h2r,h2i;
double wr,wi,wpr,wpi,wtemp,theta;
theta=3.141592653589793/(double) (n>>1);
if (isign == 1) {
c2 = -0.5;
dfour1(data,n>>1,1);
} else {
c2=0.5;
theta = -theta;
}
wtemp=sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi=sin(theta);
wr=1.0+wpr;
wi=wpi;
np3=n+3;
for (i=2;i<=(n>>2);i++) {
i4=1+(i3=np3-(i2=1+(i1=i+i-1)));
h1r=c1*(data[i1]+data[i3]);
h1i=c1*(data[i2]-data[i4]);
h2r = -c2*(data[i2]+data[i4]);
h2i=c2*(data[i1]-data[i3]);
data[i1]=h1r+wr*h2r-wi*h2i;
data[i2]=h1i+wr*h2i+wi*h2r;
data[i3]=h1r-wr*h2r+wi*h2i;
data[i4] = -h1i+wr*h2i+wi*h2r;
wr=(wtemp=wr)*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;
}
if (isign == 1) {
data[1] = (h1r=data[1])+data[2];
data[2] = h1r-data[2];
} else {
data[1]=c1*((h1r=data[1])+data[2]);
data[2]=c1*(h1r-data[2]);
dfour1(data,n>>1,-1);
break;
default:
fprintf(stdout,"Error: unsupported type in fftw_wrapper()\n");
exit(1);
break;
}
free(timeData);
fftw_free(freqData);
}
// for n odd, median given by k = (n+1)/2. for n even, median is average of k=n/2, k=n/2+1 elements
/* ********************************************************************************** */
/* */
/* Frame I/O */
......
This diff is collapsed.
This diff is collapsed.
#include <fftw3.h>
#include "BayesLine.h"
#include "BayesWave.h"
#include "BayesWaveIO.h"
......@@ -309,7 +311,7 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
if(strcmp(injmodel, "glitch") == 0)
{
data->glitchFlag = 1;
initialize_model(model, data, prior, psd, &chain->seed);
initialize_model(model, data, prior, psd, chain->seed);
for(ifo=0; ifo<NI; ifo++)model->Sn[ifo]=1.0;
FILE **glitchParams = malloc(NI*sizeof(FILE *));
......@@ -377,7 +379,7 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
if(strcmp(injmodel, "signal") == 0)
{
data->signalFlag = 1;
initialize_model(model, data, prior, psd, &chain->seed);
initialize_model(model, data, prior, psd, chain->seed);
for(ifo=0; ifo<NI; ifo++)model->Sn[ifo]=1.0;
FILE **signalParams = malloc(data->Npol*sizeof(FILE *));
......@@ -505,126 +507,6 @@ void getChannels(ProcessParamsTable *commandLine, char **channel)
}
}
int outputFrameFile(ProcessParamsTable *commandLine, struct Data *data, struct Model *model)
{
int N = data->N, halfN=N/2;
LALFrameH *frame; /* output frame */
char filename[256]; /* output file name */
REAL8TimeSeries *ht; /* input strain data */
REAL8TimeSeries *rt; /* output strain data */
//char cache_name_1[256]; /* input frame cache name */
LALCache *cache; /* input frame cache */
LALFrStream *stream=NULL; /* input stream */
/*** Get GMST for the trigger time ***/
LALStatus status;
char *chartmp=NULL;
LIGOTimeGPS GPStrig, gps_start, gps_end;
memset(&status,0,sizeof(status));
LALStringToGPS(&status,&GPStrig,LALInferenceGetProcParamVal(commandLine,"--trigtime")->value,&chartmp);
int Tobs = atoi(LALInferenceGetProcParamVal(commandLine,"--seglen")->value);
/* define time segment */
gps_start.gpsSeconds = GPStrig.gpsSeconds+2-Tobs;
gps_start.gpsNanoSeconds = 0;
gps_end.gpsSeconds = GPStrig.gpsSeconds+2;
gps_end.gpsNanoSeconds = 0;
//cache_name_1[0]='\0';
filename[0]='\0';
ht=NULL;
rt=NULL;
int i,n;
char **ifos=NULL;
ProcessParamsTable *this=commandLine;
char tmp[128];
/* Construct a list of IFOs */
int NI=0;
for(this=commandLine;this;this=this->next)
{
if(!strcmp(this->param,"--ifo")||!strcmp(this->param,"--IFO"))
{
NI++;
ifos=XLALRealloc(ifos,NI*sizeof(char *));
ifos[NI-1]=XLALStringDuplicate(this->value);
}
}
char **caches=XLALCalloc(NI,sizeof(char *));
char **channels=XLALCalloc(NI,sizeof(char *));
double *r = malloc(N*sizeof(double));
double fix = (sqrt((double)Tobs))*sqrt(0.5*(double)N)/8.0;
/* For each IFO, fetch the other options if available */
for(i=0;i<NI;i++)
{
for(n=0; n<halfN; n++)
{
r[2*n] = (data->s[i][2*n] - model->glitch[i]->templates[2*n] );
r[2*n+1] = -(data->s[i][2*n+1]- model->glitch[i]->templates[2*n+1]);
}
drealft(r-1,N,-1);
for(n=0; n<N; n++) r[n] /= fix;
/* Cache */
sprintf(tmp,"--%s-cache",ifos[i]);
this=LALInferenceGetProcParamVal(commandLine,tmp);
caches[i]=XLALStringDuplicate(this->value);
/* Channel */
sprintf(tmp,"--%s-channel",ifos[i]);
this=LALInferenceGetProcParamVal(commandLine,tmp);
channels[i]=XLALStringDuplicate(this?this->value:"Unknown channel");
cache = XLALCacheImport(caches[i]);
stream = XLALFrStreamCacheOpen(cache);
sprintf(filename,"%s_CLEAN",channels[i]);
/* read raw data */
ht = XLALFrStreamReadREAL8TimeSeries(stream, channels[i], &gps_start, XLALGPSDiff(&gps_end, &gps_start), 0);
/* units */
ht->sampleUnits = lalStrainUnit;
rt = XLALCreateREAL8TimeSeries(filename, &gps_start, 0, (double)N/data->Tobs, &ht->sampleUnits, data->N);
/* fill residual data vector */
for(n=0;n<data->N; n++) rt->data->data[n]=r[n];
/* save time series into a frame */
frame = XLALFrameNew(&gps_start, XLALGPSDiff(&gps_end, &gps_start), "BayesWaveResidual", 0, 1, 0);
XLALFrameAddREAL8TimeSeriesProcData(frame, rt);
sprintf(filename,"%c-%s_BAYESWAVE_RESIDUAL-%d-%d.gwf", ifos[i][0],ifos[i],gps_start.gpsSeconds,(int)Tobs);
printf("writing %s\n",filename);
XLALFrameWrite(frame, filename);
/* cleaning */
XLALDestroyREAL8TimeSeries(ht);
XLALDestroyREAL8TimeSeries(rt);
XLALFrameFree(frame);
}
/* close streams */
XLALFrStreamClose(stream);
free(r);
return 0;
}
/* ********************************************************************************** */
/* */
/* Output Files */
......@@ -826,7 +708,6 @@ void print_run_stats(FILE *fptr, struct Data *data, struct Chain *chain)
fprintf(fptr, " number of iterations %i\n", chain->count);
fprintf(fptr, " number of burn-in iterations %i\n", chain->burn);
fprintf(fptr, " number of BL burn-in iterations %i\n", chain->blcount);
fprintf(fptr, " RNG seed for chain %li\n", chain->seed);
fprintf(fptr, " max N wavelets %i\n", data->Dmax);
fprintf(fptr, " max Q for wavelets %.0f\n", data->Qmax);
fprintf(fptr, " number of chains %i\n", chain->NC);
......@@ -1133,37 +1014,44 @@ void print_chain_status(struct Data *data, struct Chain *chain, struct Model **m
void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf, double eta, double Tobs, int imin, int imax, double tmin, double tmax)
{
/*
TODO: invFFT comes out time-reversed.
LAL uses exp(-i2pift) for their Fourier transforms
Why doesn't drealft(ty-1,N,1) work for inverse FFT?
In the meantime, hacked invFFT by changing sign of imaginary part
*/
int i;
double x,t;
FILE *waveout = fopen(filename,"w");
double *ht = malloc(N*sizeof(double));
for(i=0; i<N; i++)ht[i]=h[i];
ht[0] = 0.0;
ht[1] = 0.0;
fftw_complex *hf = (fftw_complex *) fftw_malloc(sizeof(fftw_complex)*N);
fftw_plan reverse = fftw_plan_dft_c2r_1d(N, hf, ht, FFTW_MEASURE);
for (i=0; i < N; i++)
{
ht[i] = 0.0;
}
hf[0][0] = 0.0;
hf[0][1] = 0.0;
for(i=1; i< N/2; i++)
{
if(i>imin && i<imax)
{
x = sqrt(Snf[i]*eta);
ht[2*i] /= x;
ht[2*i+1] /= -x;
hf[i][0] = h[2*i]/x;
hf[i][1] = h[2*i+1]/x;
}
else
{
ht[2*i]=ht[2*i+1]=0.0;
hf[i][0] = 0.0;
hf[i][1] = 0.0;
}
}
//get DC and Nyquist where FFTW wants them
hf[N/2][0] = hf[0][1];
hf[N/2][1] = hf[0][1] = 0.0;
drealft(ht-1,N,-1);
double norm = 0.5*sqrt((double)N);
fftw_execute(reverse);