Merge branch 'Qbump' into 'master'

Qbump

See merge request lscsoft/bayeswave!182

src/BayesWaveMCMC.c

@@ -1950,7 +1950,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
             }
           }
 
-          if(data->amplitudePriorFlag)
+          if(data->amplitudeProposalFlag)
           {
             //loop over Npol for signal model, single IFO (nmax=nmin+1) for glitch model
             for(n=nmin; n<nmax; n++)
@@ -2279,28 +2279,31 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
           */
           else
           {
+            alpha=uniform_draw(seed);
             for(n=nmin; n<nmax; n++)
             {
               paramsx = wave_x[n]->intParams;
               paramsy = wave_y[n]->intParams;
 
               draw_wavelet_params(paramsy[ii], range, seed, wave_y[n]->dimension);
-
-              if(uniform_draw(seed)<0.5 && data->amplitudeProposalFlag)
+              if(det==-1 && model_y->Npol > 1) constrain_hplus_hcross(wave_y, ii);
+              
+              if(alpha<0.5 && data->amplitudeProposalFlag)
               {
                 if(det==-1)
                 {
                   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) );
+                  logqy += ( data->signal_amplitude_prior(paramsy[ii],  model_y->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], seed, data->Tobs, prior->range, prior->gSNRpeak);
-                  logqy += ( data->glitch_amplitude_prior(paramsy[ii],  model_x->Snf[n], data->Tobs, prior->gSNRpeak) );
+                  logqy += ( data->glitch_amplitude_prior(paramsy[ii],  model_y->Snf[n], data->Tobs, prior->gSNRpeak) );
+                  logqx += ( data->glitch_amplitude_prior(paramsx[ii],  model_x->Snf[n], data->Tobs, prior->gSNRpeak) );
                 }
               }
             }
-            if(det==-1 && model_y->Npol > 1) constrain_hplus_hcross(wave_y, ii);
           }
         }//end perturb step
         //check priors

src/BayesWaveProposal.c

@@ -69,8 +69,8 @@ void draw_glitch_amplitude(double *params, double *Snf, gsl_rng *seed, double To
   invmax = 1./max;
 
   i = (int)(params[1]*Tobs);
-  double SNRmin = range[3][0]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snf[i]*2.0/Tobs));
-  double SNRmax = range[3][1]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snf[i]*2.0/Tobs));
+  double SNRmin = range[3][0]*SineGaussianSNR(params, Snf, Tobs);
+  double SNRmax = range[3][1]*SineGaussianSNR(params, Snf, Tobs);
 
 
   k = 0;
@@ -109,7 +109,9 @@ void draw_glitch_amplitude(double *params, double *Snf, gsl_rng *seed, double To
   }
 
   //SNR defined with Sn(f) but Snf array holdes <n_i^2>
-  params[3] = SNR/sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snf[i]*2.0/Tobs));
+  //params[3] = SNR/sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snf[i]*2.0/Tobs));
+
+  params[3] *= SNR/SineGaussianSNR(params, Snf, Tobs);
 
 }
 
@@ -132,10 +134,10 @@ void draw_signal_amplitude(double *params, double *Snf, gsl_rng *seed, double To
 
   i = (int)(params[1]*Tobs);
 
-  double SNRmin = range[3][0]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snf[i]*2.0/Tobs));
-  double SNRmax = range[3][1]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snf[i]*2.0/Tobs));
-
+  double SNRmin = range[3][0]*SineGaussianSNR(params, Snf, Tobs)/params[3];
+  double SNRmax = range[3][1]*SineGaussianSNR(params, Snf, Tobs)/params[3];
 
+  
   k = 0;
   SNR = SNRmin + (SNRmax-SNRmin)*uniform_draw(seed);
 
@@ -174,7 +176,9 @@ void draw_signal_amplitude(double *params, double *Snf, gsl_rng *seed, double To
   }
 
   //SNR defined with Sn(f) but Snf array holdes <n_i^2>
-  params[3] = SNR/sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snf[i]*2.0/Tobs));
+  //params[3] = SNR/sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snf[i]*2.0/Tobs));
+
+  params[3] *= SNR/SineGaussianSNR(params, Snf, Tobs);
   
 //  FILE *temp=fopen("prior.dat","a");
 //  fprintf(temp,"%lg\n",params[3]);
@@ -195,7 +199,7 @@ void draw_uniform_amplitude(double *params, double *Snf, gsl_rng *seed, double T
 
 
   //SNR defined with Sn(f) but Snf array holdes <n_i^2>
-  params[3] = SNR/sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snf[i]*2.0/Tobs));
+  params[3] *= SNR/SineGaussianSNR(params, Snf, Tobs);
   
 }
 
@@ -492,7 +496,7 @@ void intrinsic_fisher_update(double *params, double *dparams, double *Snf, doubl
   double betasq;
 
   //SNR defined with Sn(f) but Snf array holdes <n_i^2>
-  SNR = params[3]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snf[i]*2.0/Tobs));
+  SNR = SineGaussianSNR(params, Snf, Tobs);
 
   if(SNR < 5.0) SNR = 5.0;  // this caps the size of the proposed jumps