Commit 93b84304 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

create O2_online tag

flag high SNR assymetric data as glitch
clean up build warnings



git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@642 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent 0da623f0
#!/usr/bin/python
"""
This script generates a webpage to display the results of a BayesWave run.
......@@ -2358,4 +2359,4 @@ if not os.path.exists(htmlDir):
os.makedirs(htmlDir)
# -- Make webpage
make_webpage(htmlDir, model, mdc, gps, ifoList, ifoNames, modelList, snrList, sig_gl, sig_noise, momentsPlus, momentsCross, postprocesspath)
\ No newline at end of file
make_webpage(htmlDir, model, mdc, gps, ifoList, ifoNames, modelList, snrList, sig_gl, sig_noise, momentsPlus, momentsCross, postprocesspath)
#!/usr/bin/python
import matplotlib
matplotlib.use('Agg')
import numpy as np
......
......@@ -158,13 +158,14 @@ static double logprior_gaussian(double *mean, double *sigma, double *Snf, int il
//leaving off normalizations since they cancel in Hastings ratio
lgp = 0;
for(i=ilow; i<ihigh; i++)
{
x = (mean[i] - Snf[i])/sigma[i];
lgp -= x*x;
}
return (0.5*lgp);
return (0.0*lgp);
}
static double loglike(double *respow, double *Snf, int ncut)
......@@ -193,7 +194,7 @@ void LorentzSplineFit(BayesLinePriors *priors, int zeroLogL, int steps, dataPara
double SAmaxx, SAminn, lSAmax, lSAmin;
double lQmin, lQmax, QQ;
double Aminn, Amaxx, lAmin, lAmax;
int ac0, ac1, cnt;
int ac0, ac1;
int cc0, cc1, cc2;
double *Sn, *Sbase;
double e1, e2, e3, e4;
......@@ -297,7 +298,6 @@ void LorentzSplineFit(BayesLinePriors *priors, int zeroLogL, int steps, dataPara
if(!zeroLogL) logLx = loglike(spow, Sn, ncut);
else logLx = 1.0;
cnt = 0;
passes = 0;
spass = steps/4;
......@@ -1158,7 +1158,7 @@ void create_dataParams(dataParams *data, double *f, int n)
// size of segments in Hz
// If frequency snippets are too large need longer initial runs to get convergence
data->fstep = 30;//FSTEP;//9.0;//30.0;
data->fstep = 128./data->Tobs;//30;//FSTEP;//9.0;//30.0;
// This sets the maximum number of Lorentzian lines per segment.
// For segements shorter than 16 Hz this always seems to be enough
......@@ -1174,7 +1174,7 @@ void create_dataParams(dataParams *data, double *f, int n)
data->nmin = (int)(f[0]*data->Tobs);
// the stencil separation in Hz for the spline model. Going below 2 Hz is dangerous - will fight with line model
data->fgrid = 15.0;//FSTEP;///4;//4.0;//FSTEP;
data->fgrid = data->fstep/2.;//15.0;//FSTEP;///4;//4.0;//FSTEP;
}
void create_lorentzianParams(lorentzianParams *lines, int size)
......@@ -1253,7 +1253,7 @@ void BayesLineNonMarkovianFit(struct BayesLineParams *bayesline, int *nj)
int j=0, jj=0;
int kk=0;
double mdn,sm;
double mdn;
double frs,deltafmax,spread,fsq;
double x,z;
......@@ -1317,8 +1317,6 @@ void BayesLineNonMarkovianFit(struct BayesLineParams *bayesline, int *nj)
sfreq = malloc((size_t)(sizeof(double)*(data->ncut)));
wndw = malloc((size_t)(sizeof(double)*(data->ncut)));
sm = data->ncut/2;
for(i=0; i<data->ncut; i++)
{
j = i+imin-data->nmin;
......@@ -1426,12 +1424,12 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
double *xint;
double e1, e2, e3, e4;
double x2, x3, x4;
double s1, s2, s3, s4;
double s2, s3, s4;
int typ=0;
double xsm, pmax, fcl, fch, dff;
double baseav;
double logpx=0.0, logpy=0.0, x, y, z, beta;
double logPsy,logPsx;
double logPsy=0.0,logPsx=0.0;
double Ac;
double *fprop;
double *sdatay;
......@@ -1497,8 +1495,6 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
dff = 0.01; // half-width of frequency focus region (used if focus == 1)
// this is the fractional error estimate on the noise level
s1 = 1.0/sqrt((double)(ncut));
s2 = 0.01;
s3 = 0.5;
s4 = 0.5;
......@@ -2252,7 +2248,7 @@ void BayesLineMarkovianSplineOnly(struct BayesLineParams *bayesline, int nspline
}
// produce an initial spline fit to the smooth part of the spectrum
SpecFitSpline(priors, bayesline->constantLogLFlag, 500000, fa, Sna, spline, Snf, jj, r);
SpecFitSpline(priors, bayesline->constantLogLFlag, 1000000, fa, Sna, spline, Snf, jj, r);
}
......@@ -2268,15 +2264,15 @@ void BayesLineMarkovianFocusedAnalysis(struct BayesLineParams *bayesline)
int j;
double dan;
BayesLineLorentzSplineMCMC(bayesline, 1.0, 500000, 0, 0, &dan);
BayesLineLorentzSplineMCMC(bayesline, 1.0, 1000000, 0, 0, &dan);
//alternate between targeting outliers and MCMCing full solution until outliers are gone or max iterations are reached
j = 0;
do
{
BayesLineLorentzSplineMCMC(bayesline, 1.0, 100000, 1, 0, &dan);
BayesLineLorentzSplineMCMC(bayesline, 1.0, 200000, 1, 0, &dan);
BayesLineLorentzSplineMCMC(bayesline, 1.0, 100000, 0, 0, &dan);
BayesLineLorentzSplineMCMC(bayesline, 1.0, 200000, 0, 0, &dan);
j++;
......
......@@ -404,6 +404,7 @@ int main(int argc, char *argv[])
fscanf(fptr,"%i",&data->signalModelFlag);
fscanf(fptr,"%i",&data->fullModelFlag);
fscanf(fptr,"%i",&data->cleanOnlyFlag);
fscanf(fptr,"%i",&data->loudGlitchFlag);
fclose(fptr);
/****************************************************/
......@@ -513,6 +514,9 @@ int main(int argc, char *argv[])
//SNR of signal remaining in data
double rSNR=0.0;
double ifoSNR[NI];
double snrRatio;
for(ifo=0; ifo<NI; ifo++) ifoSNR[ifo]=0.0;
struct Model *m = NULL;
struct Wavelet *g = NULL;
......@@ -534,10 +538,33 @@ int main(int argc, char *argv[])
//regress all wavelets from residual array for plotting
m->wavelet(r[ifo], g->intParams[i], N, -1, data->Tobs);
}
rSNR += fourier_nwip(data->imin,data->imax,g->templates,g->templates,m->invSnf[ifo]);
ifoSNR[ifo] += fourier_nwip(data->imin,data->imax,g->templates,g->templates,m->invSnf[ifo]);
rSNR += ifoSNR[ifo];
fprintf(stdout, " Residual SNR in IFO%i = %g\n",ifo,sqrt(ifoSNR[ifo]));
}
fprintf(stdout, " Residual SNR = %g\n",sqrt(rSNR));
/*
BayesWave can not handle very loud, very assymetric SNRs.
loudGlitchFlag is set to 1 if the SNR > 100 and the SNRratio is > 5:1
This hack is only implemented for 2-detector networks.
*/
if(NI==2)
{
snrRatio = ifoSNR[0]/ifoSNR[1];
if(snrRatio < 1.0) snrRatio = 1./snrRatio;
if(sqrt(rSNR)>100. && snrRatio > 5. )
{
fprintf(stdout,"\n");
fprintf(stdout,"************************ WARNING ************************\n");
fprintf(stdout,"Cleaning phase found loud, assymetric residual power \n");
fprintf(stdout,"Model selection cannot handle this trigger \n");
fprintf(stdout,"Assume it is a glitch but proceed with analysis \n");
fprintf(stdout,"*********************************************************\n");
data->loudGlitchFlag = 1;
}
}
// Print residual to file
for(ifo=0; ifo<data->NI; ifo++)
{
......@@ -687,7 +714,8 @@ int main(int argc, char *argv[])
data->Dmax = Dmax;
FILE *evidence;
sprintf(filename,"%sevidence.dat",data->runName);
if(data->loudGlitchFlag)sprintf(filename,"%sevidence_ignore.dat",data->runName);
else sprintf(filename,"%sevidence.dat",data->runName);
evidence = fopen(filename,"w");
/******************************************************************************/
......@@ -857,13 +885,30 @@ int main(int argc, char *argv[])
fclose(evidence);
//output file format for stacked histogram
sprintf(filename,"%sevidence_stacked.dat",data->runName);
if(data->loudGlitchFlag) sprintf(filename,"%sevidence_stacked_ignore.dat",data->runName);
else sprintf(filename,"%sevidence_stacked.dat",data->runName);
evidence = fopen(filename,"w");
fprintf(evidence,"%.12g %.12g %.12g %lg %lg %lg\n",data->logZnoise,data->logZglitch,data->logZsignal,data->varZnoise,data->varZglitch,data->varZsignal);
fclose(evidence);
//output evidence files if loudGlitchFlag is triggered
if(data->loudGlitchFlag)
{
sprintf(filename,"%sevidence_stacked.dat",data->runName);
evidence = fopen(filename,"w");
fprintf(evidence,"%.12g %.12g %.12g %lg %lg %lg\n",0.0,10000.0,0.0,0.0,1.0,0.0);
fclose(evidence);
sprintf(filename,"%sevidence.dat",data->runName);
evidence = fopen(filename,"w");
fprintf(evidence,"signal %.12g %.12g\n",0.0,0.0);
fprintf(evidence,"glitch %.12g %.12g\n",10000.0,1.0);
fprintf(evidence,"noise %.12g %.12g\n",0.0,0.0);
fclose(evidence);
}
//write gnuplot script to make single-run evidence histograms
sprintf(filename,"%sevidence.gpi",data->runName);
evidence = fopen(filename,"w");
......@@ -908,7 +953,7 @@ void print_help_message(void)
fprintf(stdout," --- Run parameters -------------------------------------------------------------\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --segment-start GPS start time of segment (trigtime + 2 - seglen) \n");
fprintf(stdout," --Niter number of iterations (2000000)\n");
fprintf(stdout," --Niter number of iterations (4000000)\n");
fprintf(stdout," --Nchain number of parallel chains (20)\n");
fprintf(stdout," --Ncycle number of model updates per RJMCMC iteration (100)\n");
fprintf(stdout," --Nburnin number of burn-in iterations (50000)\n");
......
......@@ -281,6 +281,7 @@ struct Data
int adaptTemperatureFlag;
int splineEvidenceFlag;
int constantLogLFlag;
int loudGlitchFlag;
int gnuplotFlag;
int verboseFlag;
int checkpointFlag;
......@@ -324,8 +325,8 @@ struct Data
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 (*signal_amplitude_prior) (double *, double *, double, double);
double (*glitch_amplitude_prior) (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 *, double **, struct Data*, double, double);
......
......@@ -35,16 +35,13 @@ void average_log_likelihood_via_direct_downsampling(struct Chain *chain, double
double cumACF;
double ACF;
double ACLmin;
double norm;
double mean;
double *temp = malloc(sizeof(double)*nPoints);
double *lagged = malloc(sizeof(double)*nPoints);
double **logLchain = chain->logLchain;
ACLmin = (double)nPoints/100.;
double **logLchain = chain->logLchain;
// Compute autocorrelation length of each logL chain
for(ic=0; ic<NC; ic++)
......@@ -187,7 +184,6 @@ void average_log_likelihood_via_ACL_downsampling(int M, struct Chain *chain)
}
//convert sumLogL to average logL
double beta[NC];
for(ic=0; ic<NC; ic++)
{
count=0;
......@@ -203,7 +199,6 @@ void average_log_likelihood_via_ACL_downsampling(int M, struct Chain *chain)
}
}
beta[ic] = 1./chain->temperature[ic];
chain->avgLogLikelihood[ic] /= (double)count;
chain->varLogLikelihood[ic] /= (double)count;
......@@ -240,8 +235,6 @@ void average_log_likelihood_via_recursion_downsampling(int M, struct Chain *chai
double *temp = malloc(sizeof(double)*nPoints);
double *lagged = malloc(sizeof(double)*nPoints);
double variance;
double beta[NC];
double norm;
for(ic=0; ic<NC; ic++)
{
......@@ -264,7 +257,6 @@ void average_log_likelihood_via_recursion_downsampling(int M, struct Chain *chai
if(lag==1)
{
mean = gsl_stats_mean(temp, 1, nPoints-lag);
variance = gsl_stats_variance(temp, 1, nPoints-lag);
}
for (i=0; i<nPoints-lag; i++)
{
......@@ -296,7 +288,6 @@ void average_log_likelihood_via_recursion_downsampling(int M, struct Chain *chai
/******************************************************************************/
//convert sumLogL to average logL
//double beta[NC];
for(ic=0; ic<NC; ic++)
{
count=0;
......@@ -313,7 +304,6 @@ void average_log_likelihood_via_recursion_downsampling(int M, struct Chain *chai
}
printf("number of samples in average = %i\n",count);
beta[ic] = 1./chain->temperature[ic];
chain->avgLogLikelihood[ic] /= (double)count;
chain->varLogLikelihood[ic] /= (double)count;
......@@ -340,11 +330,8 @@ void average_log_likelihood_via_ACF(struct Chain *chain)
int i,j,ic;
int NC = chain->NC;
int lag;
int imax;
int count=0;
double ACL[NC];
double ACF;
double s;
double mean;
int nPoints = chain->nPoints;
......@@ -353,9 +340,6 @@ void average_log_likelihood_via_ACF(struct Chain *chain)
double *temp = malloc(sizeof(double)*nPoints);
double norm;
//convert sumLogL to average logL
double beta[NC];
double correction;// = 0.0;
for(ic=0; ic<NC; ic++)
{
......@@ -377,11 +361,8 @@ void average_log_likelihood_via_ACF(struct Chain *chain)
for (i=0; i<nPoints; i++) temp[i] -= mean;
norm=0.0;
for(i=0; i<nPoints; i++)norm+=temp[i]*temp[i];
imax=nPoints/2;
lag=1;
ACL[ic]=1.0;
ACF=1.0;
s=1.0;
for(lag=1; lag<nPoints; lag++)
{
ACF=0.0;
......@@ -391,7 +372,6 @@ void average_log_likelihood_via_ACF(struct Chain *chain)
}
printf("number of samples in average = %i\n",count);
beta[ic] = 1./chain->temperature[ic];
chain->avgLogLikelihood[ic] /= (double)count;
chain->varLogLikelihood[ic] /= (double)count;
......@@ -841,9 +821,9 @@ void SplineThermodynamicIntegration(double *temperature, double *avgLogLikelihoo
void bootcorrposdef(double **data, int NC, int N, double **cmat, double *mu, long *seed)
{
double av, cc, vc;
double av;
double x, y, shrink;
int i, j, k, cmin, icm, cmax, icx;
int i, j, k, cmin, icm, cmax;
int ii, jj, kk, nn;
int *cyc;
int cycles;
......@@ -899,7 +879,6 @@ void bootcorrposdef(double **data, int NC, int N, double **cmat, double *mu, lon
// if b > 1, then what I'm calling cycles are actually chunks, made up of b cycles
icm = 0;
icx = NC;
do
{
kk = 0;
......@@ -936,7 +915,6 @@ void bootcorrposdef(double **data, int NC, int N, double **cmat, double *mu, lon
if(cycles > cmax)
{
cmax = cycles;
icx = ii;
}
kk += cycles;
......@@ -944,9 +922,6 @@ void bootcorrposdef(double **data, int NC, int N, double **cmat, double *mu, lon
}
cc = ((double)(kk)/(double)(NC));
vc = sqrt((double)(nn)/(double)(NC)-(double)(kk*kk)/(double)(NC*NC));
// looking for between 100 and 1000 chunks
// This assumes that we have a decent amount of data to begin with
......@@ -1276,15 +1251,6 @@ void CovarianceThermodynamicIntegration(double *temperature, double **loglikelih
// Invert the matrix
gsl_linalg_LU_invert (matrix, perm, inverse);
double val=0;
for(i=0; i<ND; i++)
{
for(j=0; j<ND; j++)
{
val=gsl_matrix_get(inverse,i,j);
}
}
// print out data with error bars
sprintf(filename,"%s_bootstrap_evidence.dat",modelname);
......@@ -1759,7 +1725,6 @@ void ThermodynamicIntegration(double *temperature, double **loglikelihood, int M
double smooth;
double a1, a2, base;
double avr = 0.0;
int nA;
N = 10000000; // number of MCMC steps
......@@ -1780,8 +1745,6 @@ void ThermodynamicIntegration(double *temperature, double **loglikelihood, int M
// maximum number of spline points - we initially overcover
NSP = 2*ND-1;
nA = ND;
datax = dvector(0,ND-1);
chains = dmatrix(0,ND-1,0,NL-1);
datay = dvector(0,ND-1);
......
......@@ -507,7 +507,7 @@ int outputFrameFile(ProcessParamsTable *commandLine, struct Data *data, struct M
REAL8TimeSeries *ht; /* input strain data */
REAL8TimeSeries *rt; /* output strain data */
char cache_name_1[256]; /* input frame cache name */
//char cache_name_1[256]; /* input frame cache name */
LALCache *cache; /* input frame cache */
LALFrStream *stream=NULL; /* input stream */
......@@ -527,7 +527,7 @@ int outputFrameFile(ProcessParamsTable *commandLine, struct Data *data, struct M
gps_end.gpsSeconds = GPStrig.gpsSeconds+2;
gps_end.gpsNanoSeconds = 0;
cache_name_1[0]='\0';
//cache_name_1[0]='\0';
filename[0]='\0';
ht=NULL;
rt=NULL;
......@@ -1257,6 +1257,8 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
data->fullModelFlag = 0;
data->cleanOnlyFlag = 0;
//assume the data is analyzable
data->loudGlitchFlag = 0;
ppt = LALInferenceGetProcParamVal(commandLine, "--window");
if(ppt) data->Twin = atof(ppt->value);
......@@ -1264,11 +1266,11 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
ppt = LALInferenceGetProcParamVal(commandLine, "--Nchain");
if(ppt) chain->NC = atoi(ppt->value);
else chain->NC = 15;
else chain->NC = 20;
ppt = LALInferenceGetProcParamVal(commandLine, "--Niter");
if(ppt) chain->count = atoi(ppt->value);
else chain->count = 2000000;
else chain->count = 4000000;
ppt = LALInferenceGetProcParamVal(commandLine, "--Ncycle");
if(ppt) chain->cycle = atoi(ppt->value);
......@@ -1599,10 +1601,10 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
{
fprintf(stdout,"************************ WARNING ************************\n");
fprintf(stdout,"waveletPrior flag overrides --Dmax \n");
fprintf(stdout,"Setting Dmax = 200 \n");
fprintf(stdout,"Setting Dmax = 100 \n");
fprintf(stdout,"Remove --waveletPrior if you really want to cap N \n");
fprintf(stdout,"*********************************************************\n");
data->Dmax = 200;
data->Dmax = 100;
}
}
......
......@@ -321,6 +321,7 @@ static void resume_sampler(struct Data *data, struct Chain *chain, struct Model
fscanf(fptr,"%i",&data->signalModelFlag);
fscanf(fptr,"%i",&data->fullModelFlag);
fscanf(fptr,"%i",&data->cleanOnlyFlag);
fscanf(fptr,"%i",&data->loudGlitchFlag);
fclose(fptr);
......@@ -657,8 +658,6 @@ static void save_sampler(struct Data *data, struct Chain *chain, struct Model **
FILE *fptr = NULL;
char command[128];
/****************************************************/
/* */
/* STATE VECTOR */
......@@ -1512,7 +1511,7 @@ void EvolveBayesLineParameters(struct Data *data, struct Model **model, struct B
wave->logp = 0.0;
for(i=1; i<=wave->size; i++)
{
wave->logp += (data->glitch_amplitude_prior(wave->intParams[wave->index[i]],model_x->Snf[ifo], model_x->Sn[ifo], data->Tobs, prior->gSNRpeak));
wave->logp += (data->glitch_amplitude_prior(wave->intParams[wave->index[i]],model_x->Snf[ifo], data->Tobs, prior->gSNRpeak));
}
}
......@@ -1520,7 +1519,7 @@ void EvolveBayesLineParameters(struct Data *data, struct Model **model, struct B
wave->logp = 0.0;
for(i=1; i<=wave->size; i++)
{
wave->logp += (data->glitch_amplitude_prior(wave->intParams[wave->index[i]],model_x->SnGeo, 1.0, data->Tobs, prior->sSNRpeak));
wave->logp += (data->glitch_amplitude_prior(wave->intParams[wave->index[i]],model_x->SnGeo, data->Tobs, prior->sSNRpeak));
}
}
......@@ -1567,7 +1566,6 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
int NI = data->NI;
/* PRIOR */
int dmin = 0;
int dmax = 1;
double Snmin = prior->Snmin;
......@@ -1640,7 +1638,6 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
{
det = -1;
dmin = prior->smin;
dmax = prior->smax;
wave_x = model_x->signal;
......@@ -1652,7 +1649,6 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
{
det = (int)floor(ran2(seed)*(float)NI);
dmin = 0;//prior->gmin[det];
dmax = prior->gmax[det];
wave_x = model_x->glitch[det];
......@@ -1746,8 +1742,8 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if(data->amplitudePriorFlag)
{
if(det==-1) logqx += ( data->signal_amplitude_prior(paramsx[ii], model_x->SnGeo, 1.0, data->Tobs, prior->sSNRpeak) );
else logqx += ( data->glitch_amplitude_prior(paramsx[ii], model_x->Snf[det], Snx[det], data->Tobs, prior->gSNRpeak) );
if(det==-1) logqx += ( data->signal_amplitude_prior(paramsx[ii], model_x->SnGeo, data->Tobs, prior->sSNRpeak) );
else logqx += ( data->glitch_amplitude_prior(paramsx[ii], model_x->Snf[det], data->Tobs, prior->gSNRpeak) );
}
}//end death move
......@@ -1834,12 +1830,12 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if(det==-1)
{
data->signal_amplitude_proposal(paramsy[ii], model_x->SnGeo, 1.0, seed, data->Tobs, prior->range, prior->sSNRpeak);
logqy += ( data->signal_amplitude_prior(paramsy[ii], model_x->SnGeo, 1.0, data->Tobs, prior->sSNRpeak) );
logqy += ( data->signal_amplitude_prior(paramsy[ii], model_x->SnGeo, data->Tobs, prior->sSNRpeak) );
}
else
{
data->glitch_amplitude_proposal(paramsy[ii], model_x->Snf[det], Sny[det], seed, data->Tobs, prior->range, prior->gSNRpeak);
logqy += (data->glitch_amplitude_prior(paramsy[ii], model_x->Snf[det], Sny[det], data->Tobs, prior->gSNRpeak) );
logqy += (data->glitch_amplitude_prior(paramsy[ii], model_x->Snf[det], data->Tobs, prior->gSNRpeak) );
}
}
else logqy += -log(prior->range[3][1]-prior->range[3][0]);
......@@ -1917,15 +1913,15 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if(det==-1)
{
data->signal_amplitude_proposal(paramsy[ii], model_x->SnGeo, 1.0, seed, data->Tobs, prior->range, prior->sSNRpeak);
logqy += (data->signal_amplitude_prior(paramsy[ii], model_x->SnGeo, 1.0, data->Tobs, prior->sSNRpeak));
logqx += (data->signal_amplitude_prior(paramsx[ii], model_x->SnGeo, 1.0, data->Tobs, prior->sSNRpeak));
logqy += (data->signal_amplitude_prior(paramsy[ii], model_x->SnGeo, data->Tobs, prior->sSNRpeak));
logqx += (data->signal_amplitude_prior(paramsx[ii], model_x->SnGeo, data->Tobs, prior->sSNRpeak));
}
else
{
data->glitch_amplitude_proposal(paramsy[ii], model_x->Snf[det], Sny[det], seed, data->Tobs, prior->range, prior->gSNRpeak);
logqy += (data->glitch_amplitude_prior(paramsy[ii], model_x->Snf[det], Sny[det], data->Tobs, prior->gSNRpeak));
logqx += (data->glitch_amplitude_prior(paramsx[ii], model_x->Snf[det], Snx[det], data->Tobs, prior->gSNRpeak));
logqy += (data->glitch_amplitude_prior(paramsy[ii], model_x->Snf[det], data->Tobs, prior->gSNRpeak));
logqx += (data->glitch_amplitude_prior(paramsx[ii], model_x->Snf[det], data->Tobs, prior->gSNRpeak));
}
}
......@@ -1973,12 +1969,12 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if(det==-1)
{
data->signal_amplitude_proposal(paramsy[ii], model_x->SnGeo, 1.0, seed, data->Tobs, prior->range, prior->sSNRpeak);
logqy += ( data->signal_amplitude_prior(paramsy[ii], model_x->SnGeo, 1.0, data->Tobs, prior->sSNRpeak) );
logqy += ( data->signal_amplitude_prior(paramsy[ii], model_x->SnGeo, data->Tobs, prior->sSNRpeak) );
}
else
{
data->glitch_amplitude_proposal(paramsy[ii], model_x->Snf[det], Sny[det], seed, data->Tobs, prior->range, prior->gSNRpeak);
logqy += (data->glitch_amplitude_prior(paramsy[ii], model_x->Snf[det], Sny[det], data->Tobs, prior->gSNRpeak) );
logqy += (data->glitch_amplitude_prior(paramsy[ii], model_x->Snf[det], data->Tobs, prior->gSNRpeak) );
}
}
}
......@@ -2030,16 +2026,16 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
for(i=1; i<=wave_y->size; i++)
{
if(det==-1)
wave_y->logp += (data->signal_amplitude_prior(paramsy[larrayy[i]], model_x->SnGeo, 1.0, data->Tobs, prior->sSNRpeak));
wave_y->logp += (data->signal_amplitude_prior(paramsy[larrayy[i]], model_x->SnGeo, data->Tobs, prior->sSNRpeak));