Commit 134c6f03 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

commented/set to 1 all psdFit params

parent a61f3b04
......@@ -115,7 +115,7 @@ struct Model
int size;
int Npol;
double *Sn;
//double *Sn;
double **Snf;
double **invSnf;
double **SnS;
......
......@@ -331,7 +331,7 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
{
data->glitchFlag = 1;
initialize_model(model, data, prior, psd, chain->seed);
for(ifo=0; ifo<NI; ifo++)model->Sn[ifo]=1.0;
//for(ifo=0; ifo<NI; ifo++)model->Sn[ifo]=1.0;
FILE **glitchParams = malloc(NI*sizeof(FILE *));
......@@ -367,7 +367,7 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
for(ifo=0; ifo<NI; ifo++)
{
GNR[ifo] = detector_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model->glitch[ifo]->templates, model->invSnf[ifo], model->Sn[ifo]);
GNR[ifo] = detector_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model->glitch[ifo]->templates, model->invSnf[ifo], 1.0);
fprintf(snrFile,"%g ",GNR[ifo]);
}
fprintf(snrFile,"\n");
......@@ -399,7 +399,7 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
{
data->signalFlag = 1;
initialize_model(model, data, prior, psd, chain->seed);
for(ifo=0; ifo<NI; ifo++)model->Sn[ifo]=1.0;
//for(ifo=0; ifo<NI; ifo++)model->Sn[ifo]=1.0;
FILE **signalParams = malloc(data->Npol*sizeof(FILE *));
......@@ -429,8 +429,8 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
double SNRnet;
double SNRifo[NI];
SNRnet = network_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model->response, model->invSnf, model->Sn, data->NI);
for(ifo=0; ifo<NI; ifo++) SNRifo[ifo] = detector_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model->response[ifo], model->invSnf[ifo], model->Sn[ifo]);
SNRnet = network_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model->response, model->invSnf, data->NI);
for(ifo=0; ifo<NI; ifo++) SNRifo[ifo] = detector_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model->response[ifo], model->invSnf[ifo], 1.0);
if(SNRnet>100.) NCflag++;
if(SNRnet>200.) NCflag++;
sprintf(filename,"%ssnr/snr.dat",data->runName);
......@@ -1013,7 +1013,7 @@ void print_model(FILE *fptr, struct Data *data, struct Chain *chain, struct Mode
fprintf(fptr,"%d %f ", chain->mc, model->logL-data->logLc+model->logLnorm);
fprintf(fptr,"%d ", model->signal[0]->size);
for(ifo=0; ifo<NI; ifo++) fprintf(fptr,"%d ", model->glitch[ifo]->size);
for(ifo=0; ifo<NI; ifo++) fprintf(fptr, "%f ",model->Sn[ifo]);
//for(ifo=0; ifo<NI; ifo++) fprintf(fptr, "%f ",model->Sn[ifo]);
if(data->stochasticFlag) fprintf(fptr,"%f %f ",model->background->logamp,model->background->index);
fprintf(fptr,"%f %f ",model->logL,model->logLnorm);
fprintf(fptr,"\n");
......@@ -1061,12 +1061,12 @@ void print_chain_status(struct Data *data, struct Chain *chain, struct Model **m
SNR = 0.0;
if(data->signalFlag) SNR = network_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model[chain->index[0]]->response, model[chain->index[0]]->invSnf, model[chain->index[0]]->Sn, data->NI);
if(data->signalFlag) SNR = network_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model[chain->index[0]]->response, model[chain->index[0]]->invSnf, data->NI);
fprintf(stdout, " logL=%f hSNR=%.0f ", model[chain->index[0]]->logL - data->logLc, SNR);
for(ifo=0; ifo<NI; ifo++)
{
SNR = 0.0;
if(data->glitchFlag) SNR = detector_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model[chain->index[0]]->glitch[ifo]->templates, model[chain->index[0]]->invSnf[ifo], model[chain->index[0]]->Sn[ifo]);
if(data->glitchFlag) SNR = detector_snr((int)floor(data->fmin*data->Tobs), (int)floor(data->fmax*data->Tobs), model[chain->index[0]]->glitch[ifo]->templates, model[chain->index[0]]->invSnf[ifo], 1.0);
fprintf(stdout, " g%iSNR=%.0f ", ifo, SNR);
}
fprintf(stdout, "\n");
......
......@@ -206,14 +206,14 @@ void recompute_residual(struct Data *data, struct Model **model, struct Chain *c
if(data->stochasticFlag)
{
//TODO: No support for glitch model in stochastic likelihood
ComputeNoiseCorrelationMatrix(data, model_x->Snf, model_x->Sn, model_x->background);
ComputeNoiseCorrelationMatrix(data, model_x->Snf, model_x->background);
model_x->logL = loglike_stochastic(data->NI, data->imin, data->imax, data->r, model_x->background->Cij, model_x->background->detCij);
}
else
{
for(ifo=0; ifo<NI; ifo++)
{
model_x->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], model_x->Sn[ifo], model_x->invSnf[ifo]);
model_x->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], 1.0, model_x->invSnf[ifo]);
model_x->logL += model_x->detLogL[ifo];
for(i=0; i<data->N/2; i++)
{
......@@ -306,7 +306,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
}
//proposed model
double *Sny = my->Sn;
//double *Sny = my->Sn;
double **ry = my->response;
double **hy = my->h;
......@@ -588,7 +588,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
if(data->stochasticFlag)
{
//TODO: No support for glitch model in stochastic likelihood
ComputeNoiseCorrelationMatrix(data, mx->Snf, my->Sn, my->background);
ComputeNoiseCorrelationMatrix(data, mx->Snf, my->background);
logLy = loglike_stochastic(NI, imin, imax, r, my->background->Cij, my->background->detCij);
}
......@@ -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(mx->Sn[ifo]))*(mx->Sn[ifo]/Sny[ifo]) -(imax-imin)*log(Sny[ifo]);
my->detLogL[ifo] = (mx->detLogL[ifo] + (imax-imin)*log(1.0))*(1.0/1.0) -(imax-imin)*log(1.0);
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], Sny[ifo], invSnf[ifo]);
logLseg = loglike(jmin, jmax, r[ifo], 1.0, 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], Sny[ifo], invSnf[ifo]);
logLseg = loglike(jmin, jmax, r[ifo], 1.0, 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], Sny[det], invSnf[det]);
logLseg = loglike(jmin, jmax, r[det], 1.0, 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], Sny[det], invSnf[det]);
logLseg = loglike(jmin, jmax, r[det], 1.0, 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], mx->Sn[ifo], invSnf[ifo]);
logLseg = loglike(jmin, jmax, r[ifo], 1.0, 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], Sny[ifo], invSnf[ifo]);
logLseg = loglike(jmin, jmax, r[ifo], 1.0, 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], mx->Sn[det], invSnf[det]);
logLseg = loglike(jmin, jmax, r[det], 1.0, 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], Sny[det], invSnf[det]);
logLseg = loglike(jmin, jmax, r[det], 1.0, 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], Sny[ifo], invSnf[ifo]);
my->detLogL[ifo] = loglike(imin, imax, r[ifo], 1.0, invSnf[ifo]);
logLy += my->detLogL[ifo];
}
}
......@@ -801,7 +801,7 @@ double EvaluateExtrinsicConstantLogLikelihood(UNUSED struct Network *projection,
else return 0.0;
}
double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, double *params, double **invSnf, double *Sn, struct Wavelet **geo, double **g, struct Data *data, double fmin, double fmax)
double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, double *params, double **invSnf, UNUSED double *Sn, 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], Sn[ifo], invSnf[ifo]);
logL += loglike(imin, imax, r[ifo], 1.0, invSnf[ifo]);
}
......@@ -863,7 +863,7 @@ double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, doubl
return logL;
}
double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *params, double **invSnf, double *Sn, struct Wavelet **geo, double **g, struct Data *data, double fmin, double fmax)
double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *params, double **invSnf, UNUSED double *Sn, struct Wavelet **geo, double **g, struct Data *data, double fmin, double fmax)
{
int i, n;
int NI,N;
......@@ -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], Sn[ifo], invSnf[ifo]);
for(ifo=0; ifo<NI; ifo++) logL += loglike(imin, imax, r[ifo], 1.0, invSnf[ifo]);
free_double_matrix(h,NI-1);
free_double_matrix(r,NI-1);
......
......@@ -37,7 +37,7 @@ void recompute_residual(struct Data *data, struct Model **model, struct Chain *c
double EvaluateConstantLogLikelihood (int typ, int ii, int det, UNUSED struct Model *mx, struct Model *my, struct Prior *prior, UNUSED struct Chain *chain, UNUSED struct Data *data);
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, double *Sn, struct Wavelet **geo, double **g, struct Data *data, double fmin, double fmax);
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);
......
......@@ -216,14 +216,14 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
}
// initialize the noise parameters
for(ifo=0; ifo<NI; ifo++) model[ic]->Sn[ifo] = 1.0;
//for(ifo=0; ifo<NI; ifo++) model[ic]->Sn[ifo] = 1.0;
// form up residual and compute initial likelihood
if(data->stochasticFlag)
{
//TODO: No support for glitch model in stochastic likelihood
ComputeNoiseCorrelationMatrix(data, model[ic]->Snf, model[ic]->Sn, model[ic]->background);
ComputeNoiseCorrelationMatrix(data, model[ic]->Snf, model[ic]->background);
model[ic]->logL = loglike_stochastic(data->NI, data->imin, data->imax, data->r, model[ic]->background->Cij, model[ic]->background->detCij);
}
......@@ -244,7 +244,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], model[ic]->Sn[ifo], model[ic]->invSnf[ifo]);
model[ic]->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], 1.0, model[ic]->invSnf[ifo]);
model[ic]->logL += model[ic]->detLogL[ifo];
for(i=0; i<data->N/2; i++)
{
......@@ -469,7 +469,7 @@ 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
//for(i=0; i<NI; i++)fscanf(fptr,"%lg",&model[ic]->Sn[i]); //PSD scale parameters
fclose(fptr);
signal = model[ic]->signal;
......@@ -628,7 +628,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], model[ic]->Sn[ifo], model[ic]->invSnf[ifo]);
model[ic]->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], 1.0, model[ic]->invSnf[ifo]);
model[ic]->logL += model[ic]->detLogL[ifo];
for(i=0; i<data->N/2; i++)
{
......@@ -818,7 +818,7 @@ 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
//for(i=0; i<NI; i++)fprintf(fptr,"%lg ",model[ic]->Sn[i]); //PSD scale parameters
fprintf(fptr,"\n");
fclose(fptr);
......@@ -1210,14 +1210,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], model[chain->index[0]]->Sn[ifo],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], 1.0, 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, model[chain->index[0]]->Sn[ifo], data->imin, data->imax);
print_frequency_domain_waveforms(filename, data->r[ifo], N, model[chain->index[0]]->Snf[ifo], data->Tobs, 1.0, 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 +1226,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],model[chain->index[0]]->Sn[ifo],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], 1.0, 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 +1235,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],model[chain->index[0]]->Sn[ifo],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], 1.0,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 +1485,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], model_x->Sn[ifo], model_x->invSnf[ifo]);
model_x->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], 1.0, model_x->invSnf[ifo]);
if(!data->constantLogLFlag)
{
model_x->logL += model_x->detLogL[ifo];
......@@ -1635,8 +1635,8 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
/* PRIOR */
int dmax = 1;
double Snmin = prior->Snmin;
double Snmax = prior->Snmax;
// double Snmin = prior->Snmin;
// double Snmax = prior->Snmax;
double **range = prior->range;
/* MODEL */
......@@ -1646,7 +1646,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
for(ifo=0; ifo<NI; ifo++) model_x->size+=model_x->glitch[ifo]->size;
double *Snx = model_x->Sn;
//double *Snx = model_x->Sn;
struct Wavelet **wave_x;
int *larrayx;
......@@ -1664,7 +1664,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
model_y->size=model_x->size;
double *Sny = model_y->Sn;
//double *Sny = model_y->Sn;
struct Wavelet **wave_y;
int *larrayy;
......@@ -1687,8 +1687,8 @@ 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;
//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;
......@@ -1973,7 +1973,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
}
else
{
data->glitch_amplitude_proposal(paramsy[ii], model_x->Snf[det], Sny[det], seed, data->Tobs, prior->range, prior->gSNRpeak);
data->glitch_amplitude_proposal(paramsy[ii], model_x->Snf[det], 1.0, seed, data->Tobs, prior->range, prior->gSNRpeak);
logqy += (data->glitch_amplitude_prior(paramsy[ii], model_x->Snf[det], data->Tobs, prior->gSNRpeak) );
}
}
......@@ -2002,8 +2002,8 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
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;
//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
......@@ -2068,7 +2068,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
}
else
{
data->glitch_amplitude_proposal(paramsy[ii], model_x->Snf[det], Sny[det], seed, data->Tobs, prior->range, prior->gSNRpeak);
data->glitch_amplitude_proposal(paramsy[ii], model_x->Snf[det], 1.0, 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));
}
......@@ -2162,7 +2162,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
}
else
{
data->glitch_amplitude_proposal(paramsy[ii], model_x->Snf[det], Sny[n], seed, data->Tobs, prior->range, prior->gSNRpeak);
data->glitch_amplitude_proposal(paramsy[ii], model_x->Snf[det], 1.0, seed, data->Tobs, prior->range, prior->gSNRpeak);
logqy += ( data->glitch_amplitude_prior(paramsy[ii], model_x->Snf[n], data->Tobs, prior->gSNRpeak) );
}
}
......@@ -2362,7 +2362,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
double logLx, logLy, logH;
double alpha;
double **glitch;
double *Sn;
double *Sn=NULL;
double *SnGeox;
double *SnGeoy;
double *paramsx, *paramsy;
......@@ -2379,7 +2379,7 @@ 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;
//Sn = model_x->Sn;
paramsx = model_x->extParams;
SnGeox = model_x->SnGeo;
......
......@@ -86,18 +86,18 @@ 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, double *eta, int NI)
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)
{
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]);///eta[i];
return(sum);
}
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 )//double network_snr(int imin, int imax, double **h, double **invpsd, double *eta, int NI )
{
return sqrt( network_nwip(imin,imax,h,h,invpsd,eta,NI) );
return sqrt( network_nwip(imin,imax,h,h,invpsd,NI) );
}
double detector_snr(int imin, int imax, double *g, double *invpsd, double eta)
......
......@@ -27,8 +27,8 @@ double gaussian_norm(double min, double max, double sigma);
void recursive_phase_evolution(double dre, double dim, double *cosPhase, double *sinPhase);
double network_nwip(int imin, int imax, double **a, double **b, double **invSn, double *eta, int NI);
double network_snr(int imin, int imax, double **h, double **invpsd, double *eta, int NI );
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_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 detector_snr(int imin, int imax, double *g, double *invpsd, double eta);
double fourier_nwip(int imin, int imax, double *a, double *b, double *invSn);
......
......@@ -137,7 +137,7 @@ void OverlapReductionFunction(struct Data *data)
}
}
void ComputeNoiseCorrelationMatrix(struct Data *data, double **Snf, double *Sn, struct Background *background)
void ComputeNoiseCorrelationMatrix(struct Data *data, double **Snf, struct Background *background)//void ComputeNoiseCorrelationMatrix(struct Data *data, double **Snf, double *Sn, struct Background *background)
{
int i;
......@@ -162,8 +162,8 @@ void ComputeNoiseCorrelationMatrix(struct Data *data, double **Snf, double *Sn,
Calculate noise correlation matrix
Combining instrument noise and stochastic backgrouns signal
*/
C11 = background->spectrum[i] + Snf[0][i]*Sn[0];
C22 = background->spectrum[i] + Snf[1][i]*Sn[1];
C11 = background->spectrum[i] + Snf[0][i];
C22 = background->spectrum[i] + Snf[1][i];
C12 = background->spectrum[i] * data->ORF[i];
C21 = C12;
......@@ -175,9 +175,6 @@ void ComputeNoiseCorrelationMatrix(struct Data *data, double **Snf, double *Sn,
background->Cij[1][0][i] = -C21*invC;
background->Cij[1][1][i] = C11*invC;
// printf("%lg %lg\n%lg %lg\n",background->Cij[0][0][i],background->Cij[0][1][i],background->Cij[1][0][i],background->Cij[1][1][i]);
// printf("================\n");
// printf("%lg %lg\n",Snf[0][i]*Sn[0],data->Snf[1][i]*Sn[1]);
}
}
......@@ -708,7 +705,7 @@ void initialize_model(struct Model *model, struct Data *data, struct Prior *prio
//gaussian noise model
model->Sn = double_vector(NI-1);
//model->Sn = double_vector(NI-1);
model->SnS = double_matrix(NI-1,halfN-1);
model->Snf = double_matrix(NI-1,halfN-1);
model->invSnf= double_matrix(NI-1,halfN-1);
......@@ -1010,7 +1007,7 @@ void copy_int_model(struct Model *origin, struct Model *copy, int N, int NI, int
{
copy->detLogL[ifo] = origin->detLogL[ifo];
copy->Sn[ifo] = origin->Sn[ifo]; //noise PSD scale
//copy->Sn[ifo] = origin->Sn[ifo]; //noise PSD scale
}
}
......@@ -1054,7 +1051,7 @@ void free_model(struct Model *model, struct Data *data, struct Prior *prior)
int NI = data->NI;
int Npol = model->Npol;
free_double_vector(model->Sn);
//free_double_vector(model->Sn);
free_double_vector(model->SnGeo);
free_double_matrix(model->SnS,NI-1);
free_double_matrix(model->Snf,NI-1);
......
......@@ -33,7 +33,7 @@ void Shf_Geocenter(struct Data *data, struct Model *model, double *SnGeo, double
void Shf_Geocenter_full(struct Data *data, struct Network *projection, double **Snf, double *SnGeo, double *params);
void OverlapReductionFunction(struct Data *data);
void ComputeNoiseCorrelationMatrix(struct Data *data, double **Snf, double *Sn, struct Background *background);
void ComputeNoiseCorrelationMatrix(struct Data *data, double **Snf, struct Background *background);//void ComputeNoiseCorrelationMatrix(struct Data *data, double **Snf, double *Sn, struct Background *background);
double symmetric_snr_ratio(struct Data *data, struct Network *projection, double *params);
......
......@@ -264,7 +264,7 @@ void intrinsic_halfcycle_proposal(double *x, double *y, gsl_rng *seed)
/* ********************************************************************************** */
void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *params, double **invSnf, double *Sn, struct Model *model, double **glitch, struct Data *data)
void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *params, double **invSnf, struct Model *model, double **glitch, struct Data *data)//void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *params, double **invSnf, double *Sn, struct Model *model, double **glitch, struct Data *data)
{
int i, j, k;
int NI,N;
......@@ -276,6 +276,8 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
double *epsilon;
double **s;
double stab;
double *DUMMY=NULL;
struct Network *projectionP = malloc(sizeof(struct Network));
struct Network *projectionM = malloc(sizeof(struct Network));
......@@ -355,9 +357,8 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
paramsP[i] += epsilon[i];
paramsM[i] -= epsilon[i];
logLP = EvaluateExtrinsicMarkovianLogLikelihood(projectionP, paramsP, invSnf, Sn, model->signal, glitch, data, data->fmin, data->fmax);
logLM = EvaluateExtrinsicMarkovianLogLikelihood(projectionM, paramsM, invSnf, Sn, model->signal, glitch, data, data->fmin, data->fmax);
logLP = EvaluateExtrinsicMarkovianLogLikelihood(projectionP, paramsP, invSnf, DUMMY, model->signal, glitch, data, data->fmin, data->fmax);
logLM = EvaluateExtrinsicMarkovianLogLikelihood(projectionM, paramsM, invSnf, DUMMY, model->signal, glitch, data, data->fmin, data->fmax);
epsilon[i] = 0.1/sqrt(-(logLM+logLP)/(epsilon[i]*epsilon[i]));
if(epsilon[i]!=epsilon[i])epsilon[i]=1.0e-5;
}
......@@ -408,7 +409,7 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
}
}