Commit 7f0c5338 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

removed psd fit parameters from functions

parent 134c6f03
......@@ -358,14 +358,14 @@ struct Data
double varZfull;
void (*signal_amplitude_proposal)(double *, double *, double, gsl_rng *, double, double **, double);
void (*glitch_amplitude_proposal)(double *, double *, double, gsl_rng *, double, double **, double);
void (*signal_amplitude_proposal)(double *, double *, gsl_rng *, double, double **, double);
void (*glitch_amplitude_proposal)(double *, double *, gsl_rng *, 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 *, struct Wavelet **, double **, struct Data*, double, double);
double (*extrinsic_likelihood)(struct Network*, double *, double **, struct Wavelet **, double **, struct Data*, double, double);
char runName[100];
char outputDirectory[1000];
......
......@@ -2373,7 +2373,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
for(i=1; i<=model->signal[0]->size; i++)
{
ii = model->signal[0]->index[i];
intrinsic_fisher_update(model->signal[0]->intParams[ii], sigmas, model->SnGeo, 1.0, data->Tobs, data->NW,1 );
intrinsic_fisher_update(model->signal[0]->intParams[ii], sigmas, model->SnGeo, data->Tobs, data->NW,1 );
for(j=0; j<NW; j++) detC += log(sigmas[j]);
}
}
......@@ -2386,7 +2386,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
{
ii = model->glitch[ifo]->index[i];
//for(j=0; j<NW; j++)printf(" %i: %g\n",j,model->glitch[ifo]->intParams[ii][j]);
intrinsic_fisher_update(model->glitch[ifo]->intParams[ii], sigmas, model->Snf[ifo], 1.0, data->Tobs, data->NW, 1);
intrinsic_fisher_update(model->glitch[ifo]->intParams[ii], sigmas, model->Snf[ifo], data->Tobs, data->NW, 1);
//for(j=0; j<NW; j++)printf(" %i: %g\n",j,sigmas[j]);
for(j=0; j<NW; j++) detC += log(sigmas[j]);
}
......@@ -2428,7 +2428,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
for(i=0; i<100; i++)
{
draw_wavelet_params(params, prior->range, chain->seed, data->NW);
data->signal_amplitude_proposal(params, model->SnGeo, 1.0, chain->seed, data->Tobs, prior->range,prior->sSNRpeak);
data->signal_amplitude_proposal(params, model->SnGeo, chain->seed, data->Tobs, prior->range,prior->sSNRpeak);
if(params[3]<minAmp)minAmp=params[3];
if(params[3]>maxAmp)maxAmp=params[3];
}
......@@ -2455,7 +2455,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
for(i=0; i<100; i++)
{
draw_wavelet_params(params, prior->range, chain->seed, data->NW);
data->glitch_amplitude_proposal(params, model->Snf[ifo], 1.0, chain->seed, data->Tobs, prior->range, prior->gSNRpeak);
data->glitch_amplitude_proposal(params, model->Snf[ifo], chain->seed, data->Tobs, prior->range, prior->gSNRpeak);
if(params[3]<minAmp)minAmp=params[3];
if(params[3]>maxAmp)maxAmp=params[3];
}
......
......@@ -1104,7 +1104,7 @@ 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)
void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf, double Tobs, int imin, int imax, double tmin, double tmax)
{
int i;
double x,t;
......@@ -1123,7 +1123,7 @@ void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf,
{
if(i>imin && i<imax)
{
x = sqrt(Snf[i]*eta);
x = sqrt(Snf[i]);
hf[i][0] = h[2*i]/x;
hf[i][1] = h[2*i+1]/x;
}
......@@ -1211,7 +1211,7 @@ void print_colored_time_domain_waveforms(char filename[], double *h, int N, doub
}
void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, double eta, double Tobs, int imin, int imax, double tmin, double tmax)
void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, double Tobs, int imin, int imax, double tmin, double tmax)
{
int i;
double x,t;
......@@ -1239,7 +1239,7 @@ void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, doub
{
if(i>imin && i<imax)
{
x = sqrt(Snf[i]*eta);
x = sqrt(Snf[i]);
// Need to switch real and imaginary parts [ie -i*(x+iy)=-ix+y]
......@@ -1277,7 +1277,7 @@ void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, doub
fclose(waveout);
}
void print_frequency_domain_waveforms(char filename[], double *h, int N, double *Snf, double Tobs, double eta, int imin, int imax)
void print_frequency_domain_waveforms(char filename[], double *h, int N, double *Snf, double Tobs, int imin, int imax)
{
int i;
double x,f;
......@@ -1288,7 +1288,7 @@ void print_frequency_domain_waveforms(char filename[], double *h, int N, double
for(i=imin; i<imax; i++)
{
x = sqrt(Snf[i]*eta);
x = sqrt(Snf[i]);
hf[2*i] /= x;
hf[2*i+1] /= x;
}
......@@ -1908,7 +1908,7 @@ void dump_time_domain_waveform(struct Data *data, struct Model *model, struct Pr
for(ifo=0; ifo<data->NI; ifo++)
{
sprintf(filename,"%s.%i",root,ifo);
print_time_domain_waveforms(filename, model->response[ifo], data->N, model->Snf[ifo], 1, data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
print_time_domain_waveforms(filename, model->response[ifo], data->N, model->Snf[ifo], data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
}
}
......
......@@ -61,9 +61,9 @@ void print_model(FILE *fptr, struct Data *data, struct Chain *chain, struct Mode
void print_signal_model(FILE *fptr, struct Model *model, int n);
void print_glitch_model(FILE *fptr, struct Wavelet *glitch);
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);
void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, double eta, double Tobs, int imin, int imax, double tmin, double tmax);
void print_frequency_domain_waveforms(char filename[], double *h, int N, double *Snf, double Tobs, double eta, int imin, int imax);
void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf, double Tobs, int imin, int imax, double tmin, double tmax);
void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, double Tobs, int imin, int imax, double tmin, double tmax);
void print_frequency_domain_waveforms(char filename[], double *h, int N, double *Snf, double Tobs, int imin, int imax);
void print_frequency_domain_data(char filename[], double *h, int N, double Tobs, int imin, int imax);
void print_colored_time_domain_waveforms(char filename[], double *h, int N, double Tobs, int imin, int imax, double tmin, double tmax);
......
......@@ -109,12 +109,12 @@ void phase_blind_time_shift(double *corr, double *corrf, double *data1, double *
}
double loglike(int imin, int imax, double *r, double Sn, double *invSnf)
double loglike(int imin, int imax, double *r, double *invSnf)
{
/*
logL = -(r|r)/2 - Nlog(Sn)
logL = -(r|r)/2
*/
return -0.5*(fourier_nwip(imin, imax, r, r, invSnf)/Sn) - (double)(imax-imin)*log(Sn);
return -0.5*fourier_nwip(imin, imax, r, r, invSnf);
}
double loglike_stochastic(int NI, int imin, int imax, double **r, double ***invCij, double *detC)
......@@ -213,7 +213,7 @@ void recompute_residual(struct Data *data, struct Model **model, struct Chain *c
{
for(ifo=0; ifo<NI; ifo++)
{
model_x->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], 1.0, model_x->invSnf[ifo]);
model_x->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], model_x->invSnf[ifo]);
model_x->logL += model_x->detLogL[ifo];
for(i=0; i<data->N/2; i++)
{
......@@ -607,7 +607,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
*/
for(ifo=0; ifo<NI; ifo++)
{
my->detLogL[ifo] = (mx->detLogL[ifo] + (imax-imin)*log(1.0))*(1.0/1.0) -(imax-imin)*log(1.0);
my->detLogL[ifo] = mx->detLogL[ifo];
logLy += my->detLogL[ifo];
}
......@@ -645,7 +645,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
for(ifo=0; ifo<NI; ifo++)
{
//remove contribution between imin and imax of old model
logLseg = loglike(jmin, jmax, r[ifo], 1.0, invSnf[ifo]);
logLseg = loglike(jmin, jmax, r[ifo], invSnf[ifo]);
my->detLogL[ifo] -= logLseg;
logLy -= logLseg;
......@@ -667,7 +667,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
r[ifo][k] -= ry[ifo][k];
}
}
logLseg = loglike(jmin, jmax, r[ifo], 1.0, invSnf[ifo]);
logLseg = loglike(jmin, jmax, r[ifo], invSnf[ifo]);
my->detLogL[ifo] += logLseg;
logLy += logLseg;
}
......@@ -675,7 +675,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
else
{
//remove contribution between imin and imax of old model
logLseg = loglike(jmin, jmax, r[det], 1.0, invSnf[det]);
logLseg = loglike(jmin, jmax, r[det], invSnf[det]);
my->detLogL[det] -= logLseg;
logLy -= logLseg;
......@@ -687,7 +687,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
r[det][j] = s[det][j] - ry[det][j] - gy[det][j];
r[det][k] = s[det][k] - ry[det][k] - gy[det][k];
}
logLseg = loglike(jmin, jmax, r[det], 1.0, invSnf[det]);
logLseg = loglike(jmin, jmax, r[det], invSnf[det]);
my->detLogL[det] += logLseg;
logLy += logLseg;
}
......@@ -716,7 +716,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
my->detLogL[ifo] = mx->detLogL[ifo];
//remove contribution between imin and imax of old model
logLseg = loglike(jmin, jmax, r[ifo], 1.0, invSnf[ifo]);
logLseg = loglike(jmin, jmax, r[ifo], invSnf[ifo]);
my->detLogL[ifo] -= logLseg;
logLy -= logLseg;
......@@ -738,7 +738,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
r[ifo][k] -= ry[ifo][k];
}
}
logLseg = loglike(jmin, jmax, r[ifo], 1.0, invSnf[ifo]);
logLseg = loglike(jmin, jmax, r[ifo], invSnf[ifo]);
my->detLogL[ifo] += logLseg;
logLy += logLseg;
}
......@@ -748,7 +748,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
my->detLogL[det] = mx->detLogL[det];
//remove contribution between imin and imax of old model
logLseg = loglike(jmin, jmax, r[det], 1.0, invSnf[det]);
logLseg = loglike(jmin, jmax, r[det], invSnf[det]);
my->detLogL[det] -= logLseg;
logLy -= logLseg;
......@@ -760,7 +760,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
r[det][j] = s[det][j] - ry[det][j] - gy[det][j];
r[det][k] = s[det][k] - ry[det][k] - gy[det][k];
}
logLseg = loglike(jmin, jmax, r[det], 1.0, invSnf[det]);
logLseg = loglike(jmin, jmax, r[det], invSnf[det]);
my->detLogL[det] += logLseg;
logLy += logLseg;
}
......@@ -782,7 +782,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
r[ifo][i] -= ry[ifo][i];
}
}
my->detLogL[ifo] = loglike(imin, imax, r[ifo], 1.0, invSnf[ifo]);
my->detLogL[ifo] = loglike(imin, imax, r[ifo], invSnf[ifo]);
logLy += my->detLogL[ifo];
}
}
......@@ -794,14 +794,14 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
}
double EvaluateExtrinsicConstantLogLikelihood(UNUSED struct Network *projection, double *params, UNUSED double **invSnf, UNUSED double *Sn, UNUSED struct Wavelet **geo, UNUSED double **g, UNUSED struct Data *data, UNUSED double fmin, UNUSED double fmax)
double EvaluateExtrinsicConstantLogLikelihood(UNUSED struct Network *projection, double *params, UNUSED double **invSnf, UNUSED struct Wavelet **geo, UNUSED double **g, UNUSED struct Data *data, UNUSED double fmin, UNUSED double fmax)
{
if(extrinsic_checkrange(params)) return -1.0e60;
else return 0.0;
}
double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, double *params, double **invSnf, UNUSED double *Sn, struct Wavelet **geo, double **g, struct Data *data, double fmin, double fmax)
double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, double *params, double **invSnf, struct Wavelet **geo, double **g, struct Data *data, double fmin, double fmax)
{
int i, n;
int NI,NP,N;
......@@ -852,7 +852,7 @@ double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, doubl
if(data->signalFlag) r[ifo][n] -= h[ifo][n];
if(data->glitchFlag) r[ifo][n] -= g[ifo][n];
}
logL += loglike(imin, imax, r[ifo], 1.0, invSnf[ifo]);
logL += loglike(imin, imax, r[ifo], invSnf[ifo]);
}
......@@ -945,7 +945,7 @@ double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *
}
}
logL = 0.0;
for(ifo=0; ifo<NI; ifo++) logL += loglike(imin, imax, r[ifo], 1.0, invSnf[ifo]);
for(ifo=0; ifo<NI; ifo++) logL += loglike(imin, imax, r[ifo], invSnf[ifo]);
free_double_matrix(h,NI-1);
free_double_matrix(r,NI-1);
......
......@@ -29,7 +29,7 @@
/* */
/* ********************************************************************************** */
double loglike(int imin, int imax, double *r, double Sn, double *Snf);
double loglike(int imin, int imax, double *r, double *Snf);
double loglike_stochastic(int NI, int imin, int imax, double **r, double ***invCij, double *detC);
void recompute_residual(struct Data *data, struct Model **model, struct Chain *chain);
......@@ -38,8 +38,8 @@ double EvaluateConstantLogLikelihood (int typ, int ii, int det, UNUSED struct Mo
double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx, struct Model *my, struct Prior *prior, struct Chain *chain, struct Data *data);
double EvaluateExtrinsicSearchLogLikelihood (struct Network *projection, double *params, double **invSnf, UNUSED double *Sn, struct Wavelet **geo, double **g, struct Data *data, double fmin, double fmax);
double EvaluateExtrinsicConstantLogLikelihood (UNUSED struct Network *projection, double *params, UNUSED double **invSnf, UNUSED double *Sn, UNUSED struct Wavelet **geo, UNUSED double **g, UNUSED struct Data *data, UNUSED double fmin, UNUSED double fmax);
double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, double *params, double **invSnf, double *Sn, struct Wavelet **geo, double **g, struct Data *data, double fmin, double fmax);
double EvaluateExtrinsicConstantLogLikelihood (UNUSED struct Network *projection, double *params, UNUSED double **invSnf, UNUSED struct Wavelet **geo, UNUSED double **g, UNUSED struct Data *data, UNUSED double fmin, UNUSED double fmax);
double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, double *params, double **invSnf, struct Wavelet **geo, double **g, struct Data *data, double fmin, double fmax);
......
......@@ -169,7 +169,7 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
for(j=1; j<=glitch[ifo]->size; j++)
{
draw_wavelet_params(glitch[ifo]->intParams[j], prior->range, chain->seed, data->NW);
if(data->amplitudePriorFlag) data->glitch_amplitude_proposal(glitch[ifo]->intParams[j], model[ic]->Snf[ifo], 1.0, chain->seed, data->Tobs, prior->range, prior->gSNRpeak);
if(data->amplitudePriorFlag) data->glitch_amplitude_proposal(glitch[ifo]->intParams[j], model[ic]->Snf[ifo], chain->seed, data->Tobs, prior->range, prior->gSNRpeak);
model[ic]->wavelet(glitch[ifo]->templates, glitch[ifo]->intParams[j], data->N, 1, Tobs);
}
......@@ -205,7 +205,7 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
{
draw_wavelet_params(signal[n]->intParams[j], prior->range, chain->seed, data->NW);
if(data->amplitudePriorFlag)
data->signal_amplitude_proposal(signal[n]->intParams[j], model[ic]->SnGeo, 1.0, chain->seed, data->Tobs, prior->range, prior->sSNRpeak);
data->signal_amplitude_proposal(signal[n]->intParams[j], model[ic]->SnGeo, chain->seed, data->Tobs, prior->range, prior->sSNRpeak);
model[ic]->wavelet(signal[n]->templates, signal[n]->intParams[j], data->N, 1, Tobs);
}
......@@ -215,9 +215,6 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
waveformProject(data, model[ic]->projection, model[ic]->extParams, model[ic]->response, model[ic]->h, data->fmin, data->fmax);
}
// initialize the noise parameters
//for(ifo=0; ifo<NI; ifo++) model[ic]->Sn[ifo] = 1.0;
// form up residual and compute initial likelihood
if(data->stochasticFlag)
......@@ -244,7 +241,7 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
}
if(!data->constantLogLFlag)
{
model[ic]->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], 1.0, model[ic]->invSnf[ifo]);
model[ic]->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], model[ic]->invSnf[ifo]);
model[ic]->logL += model[ic]->detLogL[ifo];
for(i=0; i<data->N/2; i++)
{
......@@ -469,7 +466,6 @@ static void resume_sampler(struct Data *data, struct Chain *chain, struct Model
fscanf(fptr,"%i %lg",&i,&logL); //iteration and logL
fscanf(fptr,"%i",&NS); //signal model dimension
for(i=0; i<NI; i++)fscanf(fptr,"%i",&NG[i]); //glitch model dimension
//for(i=0; i<NI; i++)fscanf(fptr,"%lg",&model[ic]->Sn[i]); //PSD scale parameters
fclose(fptr);
signal = model[ic]->signal;
......@@ -628,7 +624,7 @@ static void resume_sampler(struct Data *data, struct Chain *chain, struct Model
}
if(!data->constantLogLFlag)
{
model[ic]->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], 1.0, model[ic]->invSnf[ifo]);
model[ic]->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], model[ic]->invSnf[ifo]);
model[ic]->logL += model[ic]->detLogL[ifo];
for(i=0; i<data->N/2; i++)
{
......@@ -818,7 +814,6 @@ static void save_sampler(struct Data *data, struct Chain *chain, struct Model **
fprintf(fptr,"%i %.12g ",i,model[ic]->logL); //iteration and logL
fprintf(fptr,"%i ",signal[0]->size); //signal model dimension
for(i=0; i<NI; i++)fprintf(fptr,"%i ",glitch[i]->size); //glitch model dimension
//for(i=0; i<NI; i++)fprintf(fptr,"%lg ",model[ic]->Sn[i]); //PSD scale parameters
fprintf(fptr,"\n");
fclose(fptr);
......@@ -1210,14 +1205,14 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
for(ifo=0; ifo<NI; ifo++)
{
sprintf(filename,"waveforms/%s_data_%s_%d.dat", modelname, data->ifos[ifo], frame);
print_time_domain_waveforms(filename, data->s[ifo], N, model[chain->index[0]]->Snf[ifo], 1.0, data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
print_time_domain_waveforms(filename, data->s[ifo], N, model[chain->index[0]]->Snf[ifo], data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
if(data->bayesLineFlag)
{
for(i=0; i<N; i++) data->r[ifo][i] = data->s[ifo][i] - model[chain->index[0]]->glitch[ifo]->templates[i];
sprintf(filename,"waveforms/%s_frequency_residual_%s_%d.dat",modelname,data->ifos[ifo],frame);
print_frequency_domain_waveforms(filename, data->r[ifo], N, model[chain->index[0]]->Snf[ifo], data->Tobs, 1.0, data->imin, data->imax);
print_frequency_domain_waveforms(filename, data->r[ifo], N, model[chain->index[0]]->Snf[ifo], data->Tobs, data->imin, data->imax);
sprintf(filename,"waveforms/%s_psd_%s_%d.dat", modelname, data->ifos[ifo], frame);
print_bayesline_spectrum(filename, bayesline[chain->index[0]][ifo]->freq, bayesline[chain->index[0]][ifo]->spow, bayesline[chain->index[0]][ifo]->Sbase, bayesline[chain->index[0]][ifo]->Sline, bayesline[chain->index[0]][ifo]->data->ncut);
......@@ -1226,7 +1221,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
if(data->glitchFlag)
{
sprintf(filename,"waveforms/%s_glitch_%s_%d.dat", modelname, data->ifos[ifo], frame);
print_time_domain_waveforms(filename, model[chain->index[0]]->glitch[ifo]->templates, N, model[chain->index[0]]->Snf[ifo], 1.0, data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
print_time_domain_waveforms(filename, model[chain->index[0]]->glitch[ifo]->templates, N, model[chain->index[0]]->Snf[ifo], data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
sprintf(filename,"waveforms/%s_colored_glitch_%s_%d.dat", modelname, data->ifos[ifo], frame);
print_colored_time_domain_waveforms(filename, model[chain->index[0]]->glitch[ifo]->templates, N,data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
......@@ -1235,7 +1230,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
if(data->signalFlag)
{
sprintf(filename,"waveforms/%s_wave_%s_%d.dat", modelname, data->ifos[ifo], frame);
print_time_domain_waveforms(filename, model[chain->index[0]]->response[ifo], N, model[chain->index[0]]->Snf[ifo], 1.0,data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
print_time_domain_waveforms(filename, model[chain->index[0]]->response[ifo], N, model[chain->index[0]]->Snf[ifo], data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
sprintf(filename,"waveforms/%s_colored_wave_%s_%d.dat", modelname, data->ifos[ifo], frame);
print_colored_time_domain_waveforms(filename, model[chain->index[0]]->response[ifo], N, data->Tobs, data->imin, data->imax, prior->range[0][0], prior->range[0][1]);
......@@ -1485,7 +1480,7 @@ void EvolveBayesLineParameters(struct Data *data, struct Model **model, struct B
if(data->glitchFlag) data->r[ifo][i] -= model_x->glitch[ifo]->templates[i];
}
model_x->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], 1.0, model_x->invSnf[ifo]);
model_x->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], model_x->invSnf[ifo]);
if(!data->constantLogLFlag)
{
model_x->logL += model_x->detLogL[ifo];
......@@ -1635,8 +1630,6 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
/* PRIOR */
int dmax = 1;
// double Snmin = prior->Snmin;
// double Snmax = prior->Snmax;
double **range = prior->range;
/* MODEL */
......@@ -1645,9 +1638,6 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
model_x->size=model_x->signal[0]->size;
for(ifo=0; ifo<NI; ifo++) model_x->size+=model_x->glitch[ifo]->size;
//double *Snx = model_x->Sn;
struct Wavelet **wave_x;
int *larrayx;
double **paramsx;
......@@ -1664,8 +1654,6 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
model_y->size=model_x->size;
//double *Sny = model_y->Sn;
struct Wavelet **wave_y;
int *larrayy;
double **paramsy;
......@@ -1687,8 +1675,6 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
//acc gets set to 1 if model_y gets accepted (no need to copy at next iteration)
if(acc==0 && mc>1) copy_int_model(model_x, model_y, N, data->NI, det);
//if(data->psdFitFlag) for(ifo=0; ifo<NI; ifo++) Sny[ifo] = Snx[ifo];
//else for(ifo=0; ifo<NI; ifo++) Sny[ifo] = 1.0;
//reset acceptence counter
acc=0;
......@@ -1968,12 +1954,12 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
larrayx = wave_x[n]->index;
larrayy = wave_y[n]->index;
data->signal_amplitude_proposal(paramsy[ii], model_x->SnGeo, 1.0, seed, data->Tobs, prior->range, prior->sSNRpeak);
data->signal_amplitude_proposal(paramsy[ii], model_x->SnGeo, seed, data->Tobs, prior->range, 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], 1.0, seed, data->Tobs, prior->range, prior->gSNRpeak);
data->glitch_amplitude_proposal(paramsy[ii], model_x->Snf[det], seed, data->Tobs, prior->range, prior->gSNRpeak);
logqy += (data->glitch_amplitude_prior(paramsy[ii], model_x->Snf[det], data->Tobs, prior->gSNRpeak) );
}
}
......@@ -2000,11 +1986,6 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
{
typ = 1;
for(ifo=0; ifo<NI; ifo++)
{
//if(data->psdFitFlag) Sny[ifo] = Snx[ifo]+gaussian_draw(seed)/sqrt(0.5*(double)N);
//if(Sny[ifo] > Snmax || Sny[ifo] < Snmin) test = 1;
}
//update stochastic background parameters
if(data->stochasticFlag) stochastic_background_proposal(model_x->background, model_y->background, seed, &test);
......@@ -2061,14 +2042,14 @@ 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);
data->signal_amplitude_proposal(paramsy[ii], model_x->SnGeo, seed, data->Tobs, prior->range, 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], 1.0, seed, data->Tobs, prior->range, prior->gSNRpeak);
data->glitch_amplitude_proposal(paramsy[ii], model_x->Snf[det], seed, data->Tobs, prior->range, 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));
}
......@@ -2130,12 +2111,12 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if(det==-1)
{
intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->SnGeo, 1.0, &logqx, &logqy, seed, data->Tobs, &test, wave_y[n]->dimension, tfqFlag);
intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->SnGeo, &logqx, &logqy, seed, data->Tobs, &test, wave_y[n]->dimension, tfqFlag);
}
else
{
intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->Snf[det], 1.0, &logqx, &logqy, seed, data->Tobs, &test, wave_y[n]->dimension, tfqFlag);
intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->Snf[det], &logqx, &logqy, seed, data->Tobs, &test, wave_y[n]->dimension, tfqFlag);
}
}
if(det==-1 && model_y->Npol > 1) constrain_hplus_hcross(wave_y, ii);
......@@ -2157,12 +2138,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);
data->signal_amplitude_proposal(paramsy[ii], model_x->SnGeo, seed, data->Tobs, prior->range, 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], 1.0, seed, data->Tobs, prior->range, prior->gSNRpeak);
data->glitch_amplitude_proposal(paramsy[ii], model_x->Snf[det], seed, data->Tobs, prior->range, prior->gSNRpeak);
logqy += ( data->glitch_amplitude_prior(paramsy[ii], model_x->Snf[n], data->Tobs, prior->gSNRpeak) );
}
}
......@@ -2362,7 +2343,6 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
double logLx, logLy, logH;
double alpha;
double **glitch;
double *Sn=NULL;
double *SnGeox;
double *SnGeoy;
double *paramsx, *paramsy;
......@@ -2379,7 +2359,6 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
// extrinsic parameters
struct Model *model_x = model[chain->index[ic]];
struct Wavelet **smodel = model_x->signal;
//Sn = model_x->Sn;
paramsx = model_x->extParams;
SnGeox = model_x->SnGeo;
......@@ -2426,7 +2405,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
double fmin = data->fmin;
double fmax = data->fmax;
logLx = data->extrinsic_likelihood(model_x->projection, paramsx, model_x->invSnf, Sn, smodel, glitch, data, fmin, fmax);
logLx = data->extrinsic_likelihood(model_x->projection, paramsx, model_x->invSnf, smodel, glitch, data, fmin, fmax);
//Compute Fisher Matrix for current sky location
struct FisherMatrix *fisher = model_x->fisher;
......@@ -2505,7 +2484,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if(test==0)
{
logLy = data->extrinsic_likelihood(model_x->projection, paramsy, model_x->invSnf, Sn, smodel, glitch, data, fmin, fmax);
logLy = data->extrinsic_likelihood(model_x->projection, paramsy, model_x->invSnf, smodel, glitch, data, fmin, fmax);
//amplitude prior
if(data->amplitudePriorFlag)
......
......@@ -86,16 +86,16 @@ double fourier_nwip(int imin, int imax, double *a, double *b, double *invSn)
return(2.0*arg);
}
double network_nwip(int imin, int imax, double **a, double **b, double **invSn, int NI)//double network_nwip(int imin, int imax, double **a, double **b, double **invSn, double *eta, int NI)
double network_nwip(int imin, int imax, double **a, double **b, double **invSn, int NI)
{
int i;
double sum=0.0;
for(i=0; i<NI; i++) sum += fourier_nwip(imin,imax,a[i],b[i],invSn[i]);///eta[i];
for(i=0; i<NI; i++) sum += fourier_nwip(imin,imax,a[i],b[i],invSn[i]);
return(sum);
}
double network_snr(int imin, int imax, double **h, double **invpsd, int NI )//double network_snr(int imin, int imax, double **h, double **invpsd, double *eta, int NI )
double network_snr(int imin, int imax, double **h, double **invpsd, int NI )
{
return sqrt( network_nwip(imin,imax,h,h,invpsd,NI) );
}
......
......@@ -27,8 +27,8 @@ double gaussian_norm(double min, double max, double sigma);