diff --git a/src/BayesWave.h b/src/BayesWave.h index 349d8b0640fd0feebc6878cf460f465e5e012de0..b11c3c93142ddbb9cb6dd5bbf6865aad78468ccb 100644 --- a/src/BayesWave.h +++ b/src/BayesWave.h @@ -115,7 +115,7 @@ struct Model int size; int Npol; - double *Sn; + //double *Sn; double **Snf; double **invSnf; double **SnS; diff --git a/src/BayesWaveIO.c b/src/BayesWaveIO.c index a5e8a1751dfc956ad7621e830ba61a124e4d7e9b..c3eedd23b2cd5ebed1aa4b0fae96d29a15f2a828 100644 --- a/src/BayesWaveIO.c +++ b/src/BayesWaveIO.c @@ -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"); diff --git a/src/BayesWaveLikelihood.c b/src/BayesWaveLikelihood.c index d972ccffe9104b7303da4932255dbb522eb5bd37..4903ffe57d6a27cb77adbe97dec08072d7f2ec62 100644 --- a/src/BayesWaveLikelihood.c +++ b/src/BayesWaveLikelihood.c @@ -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); diff --git a/src/BayesWaveLikelihood.h b/src/BayesWaveLikelihood.h index 7389d3cf72c48c20b0018850e1c92c2092bf500c..692a353151916a6ea01fceb7b902421e9953d1e5 100644 --- a/src/BayesWaveLikelihood.h +++ b/src/BayesWaveLikelihood.h @@ -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); diff --git a/src/BayesWaveMCMC.c b/src/BayesWaveMCMC.c index 6227dea7bc113ade72442c704ffa7c814166f7c2..60e4103bf268199612af6f3c784ecf04ee967765 100644 --- a/src/BayesWaveMCMC.c +++ b/src/BayesWaveMCMC.c @@ -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; diff --git a/src/BayesWaveMath.c b/src/BayesWaveMath.c index e4c2b926279917b99f6159efb6259648a1224673..e8f8f478d0853a842e4517dc92bb3c11542ac2ed 100644 --- a/src/BayesWaveMath.c +++ b/src/BayesWaveMath.c @@ -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) diff --git a/src/BayesWaveMath.h b/src/BayesWaveMath.h index e9a31fbb7454f0567065cfb1ba87fd098fb6de58..6bf28988234f1122106daa134b3a3d41da2ed900 100644 --- a/src/BayesWaveMath.h +++ b/src/BayesWaveMath.h @@ -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); diff --git a/src/BayesWaveModel.c b/src/BayesWaveModel.c index a3a2b008d2a2f9f34acfd739bbc6d895722d42e5..92deb6dc57e432433f8d5eeebc98072550379c55 100644 --- a/src/BayesWaveModel.c +++ b/src/BayesWaveModel.c @@ -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); diff --git a/src/BayesWaveModel.h b/src/BayesWaveModel.h index 38cceafaa47aa6161cdbc77ae56cee620e323b89..53ce57d771efdb000b16b32c5af9252e5ebfe138 100644 --- a/src/BayesWaveModel.h +++ b/src/BayesWaveModel.h @@ -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); diff --git a/src/BayesWaveProposal.c b/src/BayesWaveProposal.c index 0934c5c9d2e6f71a27da87c81beec8566015efba..1903041e95b10248c0d591c4ac410d8cb3a4edad 100644 --- a/src/BayesWaveProposal.c +++ b/src/BayesWaveProposal.c @@ -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 } } - for(i=0; i<NE; i++) for(j=0; j<NE; j++) fisher->matrix[i][j] = network_nwip(data->imin,data->imax,hdiv[i], hdiv[j], invSnf, Sn, NI); + for(i=0; i<NE; i++) for(j=0; j<NE; j++) fisher->matrix[i][j] = network_nwip(data->imin,data->imax,hdiv[i], hdiv[j], invSnf, NI); // stabilizers for(i=0; i<NE; i++) fisher->matrix[i][i] += stab; @@ -439,7 +440,7 @@ void extrinsic_fisher_update(struct Data *data, struct Model *model) int i,j; double *params = model->extParams; - double *Sn = model->Sn; + //double *Sn = model->Sn; double **glitch = malloc(data->NI*sizeof(double*)); for(i=0; i<data->NI; i++) @@ -451,7 +452,7 @@ void extrinsic_fisher_update(struct Data *data, struct Model *model) struct FisherMatrix *fisher = model->fisher; //calculate elements of Fisher - extrinsic_fisher_information_matrix(fisher, params, model->invSnf, Sn, model, glitch, data); + extrinsic_fisher_information_matrix(fisher, params, model->invSnf, model, glitch, data); matrix_eigenstuff(fisher->matrix, fisher->evector, fisher->evalue, fisher->N); diff --git a/src/BayesWaveProposal.h b/src/BayesWaveProposal.h index a602697159ea0e8b8d374eee4ab54e11710001fb..52a5ed46c49502c824e8391ca94513ea3d63b26b 100644 --- a/src/BayesWaveProposal.h +++ b/src/BayesWaveProposal.h @@ -47,7 +47,7 @@ void intrinsic_halfcycle_proposal(double *x, double *y, gsl_rng *seed); void fisher_matrix_proposal(struct FisherMatrix *fisher, double *params, gsl_rng *seed, double *y); -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); void extrinsic_fisher_update(struct Data *data, struct Model *model);