BayesWaveProposal.c 46.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
//TODO: Amplitude prior and proposals need to call same functions!

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_statistics.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_blas.h>

#include "BayesWave.h"
16
#include "BayesLine.h"
17
18
19
20
21
22
23
#include "BayesWaveMath.h"
#include "BayesWavePrior.h"
#include "BayesWaveModel.h"
#include "BayesWaveWavelet.h"
#include "BayesWaveProposal.h"
#include "BayesWaveLikelihood.h"

24
25
26
27
28
29
#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

30
31
32
33
34
35
/* ********************************************************************************** */
/*                                                                                    */
/*                                    Prior proposals                                 */
/*                                                                                    */
/* ********************************************************************************** */

Tyson Littenberg's avatar
Tyson Littenberg committed
36
void draw_glitch_amplitude(double *params, double *Snf, double Sa, gsl_rng *seed, double Tobs, double **range, double SNRpeak)
37
{
Tyson Littenberg's avatar
Tyson Littenberg committed
38
39
40
41
  int i, k;
  double SNR, den=-1.0, alpha=1.0;
  double max;
  double invmax;
42

Tyson Littenberg's avatar
Tyson Littenberg committed
43
  double dfac, dfac3;
44

Tyson Littenberg's avatar
Tyson Littenberg committed
45
  Sa = 1.0;
46

Tyson Littenberg's avatar
Tyson Littenberg committed
47
48
  double SNR2 = 2.0*SNRpeak;
  double SNRsq = 2.0*SNRpeak*SNRpeak;
49

Tyson Littenberg's avatar
Tyson Littenberg committed
50
51
52
53
  dfac = 1.+SNRpeak/(SNR2);
  dfac3 = dfac*dfac*dfac;
  max = (SNRpeak)/(2.*SNRsq*dfac3);
  invmax = 1./max;
54

Tyson Littenberg's avatar
Tyson Littenberg committed
55
  i = (int)(params[1]*Tobs);
56

57
58
  double SNRmin = range[3][0]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
  double SNRmax = range[3][1]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
59
60


Tyson Littenberg's avatar
Tyson Littenberg committed
61
  k = 0;
Tyson Littenberg's avatar
Tyson Littenberg committed
62
  SNR = SNRmin + (SNRmax-SNRmin)*uniform_draw(seed);
63

Tyson Littenberg's avatar
Tyson Littenberg committed
64
65
  dfac = 1.+SNR/(SNR2);
  dfac3 = dfac*dfac*dfac;
66

Tyson Littenberg's avatar
Tyson Littenberg committed
67
  den = (SNR)/(SNRsq*dfac3);
68

Tyson Littenberg's avatar
Tyson Littenberg committed
69
  den *= invmax;
70

Tyson Littenberg's avatar
Tyson Littenberg committed
71
  alpha = uniform_draw(seed);
Tyson Littenberg's avatar
Tyson Littenberg committed
72
73
  while(alpha > den)
  {
74

Tyson Littenberg's avatar
Tyson Littenberg committed
75
    SNR = SNRmin + (SNRmax-SNRmin)*uniform_draw(seed);
76

Tyson Littenberg's avatar
Tyson Littenberg committed
77
78
    dfac = 1.+SNR/(SNR2);
    dfac3 = dfac*dfac*dfac;
79

Tyson Littenberg's avatar
Tyson Littenberg committed
80
    den = (SNR)/(2.*SNRsq*dfac3);
81

Tyson Littenberg's avatar
Tyson Littenberg committed
82
    den *= invmax;
83

Tyson Littenberg's avatar
Tyson Littenberg committed
84
    alpha = uniform_draw(seed);
85

Tyson Littenberg's avatar
Tyson Littenberg committed
86
    k++;
87
88
89
90
91
92
93

    //you had your chance
    if(k>10000)
    {
      SNR=0.0;
      break;
    }
Tyson Littenberg's avatar
Tyson Littenberg committed
94
  }
95

Tyson Littenberg's avatar
Tyson Littenberg committed
96
97
  //SNR defined with Sn(f) but Snf array holdes <n_i^2>
  params[3] = SNR/sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
98
99
100

}

Tyson Littenberg's avatar
Tyson Littenberg committed
101
void draw_signal_amplitude(double *params, double *Snf, double Sa, gsl_rng *seed, double Tobs, double **range, double SNRpeak)
102
{
Tyson Littenberg's avatar
Tyson Littenberg committed
103
104
105
106
  int i, k;
  double SNR, den=-1.0, alpha=1.0;
  double max;
  double invmax;
107

Tyson Littenberg's avatar
Tyson Littenberg committed
108
  double dfac, dfac5;
109

Tyson Littenberg's avatar
Tyson Littenberg committed
110
  Sa = 1.0;
111

Tyson Littenberg's avatar
Tyson Littenberg committed
112
113
  double SNR4 = 4.0*SNRpeak;
  double SNRsq = 4.0*SNRpeak*SNRpeak;
114

Tyson Littenberg's avatar
Tyson Littenberg committed
115
116
117
118
  dfac = 1.+SNRpeak/(SNR4);
  dfac5 = dfac*dfac*dfac*dfac*dfac;
  max = (3.*SNRpeak)/(SNRsq*dfac5);
  invmax = 1./max;
119

Tyson Littenberg's avatar
Tyson Littenberg committed
120
  i = (int)(params[1]*Tobs);
121

122
123
  double SNRmin = range[3][0]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
  double SNRmax = range[3][1]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
124
125


Tyson Littenberg's avatar
Tyson Littenberg committed
126
  k = 0;
Tyson Littenberg's avatar
Tyson Littenberg committed
127
  SNR = SNRmin + (SNRmax-SNRmin)*uniform_draw(seed);
128

Tyson Littenberg's avatar
Tyson Littenberg committed
129
130
  dfac = 1.+SNR/(SNR4);
  dfac5 = dfac*dfac*dfac*dfac*dfac;
131

Tyson Littenberg's avatar
Tyson Littenberg committed
132
  den = (3.*SNR)/(SNRsq*dfac5);
133

Tyson Littenberg's avatar
Tyson Littenberg committed
134
  den *= invmax;
135

Tyson Littenberg's avatar
Tyson Littenberg committed
136
  alpha = uniform_draw(seed);
Tyson Littenberg's avatar
Tyson Littenberg committed
137
138
  while(alpha > den)
  {
139

Tyson Littenberg's avatar
Tyson Littenberg committed
140
    SNR = SNRmin + (SNRmax-SNRmin)*uniform_draw(seed);
141

Tyson Littenberg's avatar
Tyson Littenberg committed
142
143
    dfac = 1.+SNR/(SNR4);
    dfac5 = dfac*dfac*dfac*dfac*dfac;
144

Tyson Littenberg's avatar
Tyson Littenberg committed
145
146
147
148
    den = (3.*SNR)/(SNRsq*dfac5);

    den *= invmax;

Tyson Littenberg's avatar
Tyson Littenberg committed
149
    alpha = uniform_draw(seed);
Tyson Littenberg's avatar
Tyson Littenberg committed
150
151
152

    k++;

153
154
155
156
157
158
159
160
    //you had your chance
    if(k>10000)
    {
      SNR=0.0;
      break;
    }


Tyson Littenberg's avatar
Tyson Littenberg committed
161
162
163
164
  }

  //SNR defined with Sn(f) but Snf array holdes <n_i^2>
  params[3] = SNR/sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
165
166
167
168
  
//  FILE *temp=fopen("prior.dat","a");
//  fprintf(temp,"%lg\n",params[3]);
//  fclose(temp);
169
170
171

}

Tyson Littenberg's avatar
Tyson Littenberg committed
172
void draw_uniform_amplitude(double *params, double *Snf, double Sa, gsl_rng *seed, double Tobs, double **range)
173
174
175
176
177
178
179
180
181
182
{
  int i;
  double SNR;

  Sa = 1.0;

  i = (int)(params[1]*Tobs);

  double SNRmin = range[3][0];//*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
  double SNRmax = range[3][1];//*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
Tyson Littenberg's avatar
Tyson Littenberg committed
183
  SNR = SNRmin + (SNRmax-SNRmin)*uniform_draw(seed);
184
185
186
187
188
189
190


  //SNR defined with Sn(f) but Snf array holdes <n_i^2>
  params[3] = SNR/sqrt((params[2]/(2.0*RT2PI*params[1]))/(Sa*Snf[i]*2.0/Tobs));
  
}

Tyson Littenberg's avatar
Tyson Littenberg committed
191
void extrinsic_uniform_proposal(gsl_rng *seed, double *y)
192
{
Tyson Littenberg's avatar
Tyson Littenberg committed
193
194
195
196
197
198
  y[0] = uniform_draw(seed)*LAL_TWOPI;
  y[1] = -1.0 + 2.0*uniform_draw(seed);
  y[2] = uniform_draw(seed)*LAL_PI_2;
  y[3] = -0.99 + 1.98*uniform_draw(seed);
  y[4] = uniform_draw(seed)*LAL_TWOPI;
  //y[4] = uniform_draw(seed)*LAL_PI;
199
200
201
  y[5] = 1.0;
}

Tyson Littenberg's avatar
Tyson Littenberg committed
202
void uniform_orientation_proposal(double *y, gsl_rng *seed)
203
204
205
206
{

  double newPsi,newPhi,newEcc;

Tyson Littenberg's avatar
Tyson Littenberg committed
207
208
209
  newPhi = uniform_draw(seed)*LAL_TWOPI;
  newPsi = uniform_draw(seed)*LAL_PI_2;
  newEcc = -1. + 2.*uniform_draw(seed);
210
211
212
213
214
215
216

  y[2] = newPsi;
  y[3] = newEcc;
  y[4] = newPhi;

}

Tyson Littenberg's avatar
Tyson Littenberg committed
217
void intrinsic_halfcycle_proposal(double *x, double *y, gsl_rng *seed)
218
{
219
  double n;
220
  double t0,f0,ph,dt,dp;
221
  
222
223
224
  t0 = x[0];
  f0 = x[1];
  ph = x[4];
225
226
  
  //choose number of half period shifts
Tyson Littenberg's avatar
Tyson Littenberg committed
227
  n = -2.0 + floor(5.0*uniform_draw(seed));
228
229
  
  //shift time by half period
Tyson Littenberg's avatar
Tyson Littenberg committed
230
  dt = (n/2.)/f0*(1. + 0.1*gaussian_draw(seed));
231
  dp = 0.0;
Tyson Littenberg's avatar
Tyson Littenberg committed
232
  if((int)n%2!=0) dp = LAL_PI*(1. + 0.1*gaussian_draw(seed));
233
234
235
236
  
  y[0] = t0 + dt;
  y[4] = ph + dp;
  
237
238
  while(y[4] > LAL_TWOPI) y[4] -=LAL_TWOPI;
  while(y[4] < LAL_TWOPI) y[4] +=LAL_TWOPI;
239
  
240
}
241
242
243
244
245
246
247
248

/* ********************************************************************************** */
/*                                                                                    */
/*                                Fisher matrix proposals                             */
/*                                                                                    */
/* ********************************************************************************** */


249
void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *params, double **invSnf, double *Sn, struct Model *model, double **glitch, struct Data *data)
250
{
Tyson Littenberg's avatar
Tyson Littenberg committed
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
  int i, j, k;
  int NI,N;

  double *paramsP, *paramsM;
  double logLP, logLM;
  double **hP, **hM, **hPP, **hMM;
  double ***hdiv;
  double *epsilon;
  double **s;
  double stab;

  struct Network *projectionP = malloc(sizeof(struct Network));
  struct Network *projectionM = malloc(sizeof(struct Network));

  NI   = data->NI;
  N    = data->N;

  initialize_network(projectionP,N, NI);
  initialize_network(projectionM,N, NI);

  stab = 10.0;

  // Plus and minus parameters:
Tyson Littenberg's avatar
Tyson Littenberg committed
274
275
  paramsP = double_vector(NE-1);
  paramsM = double_vector(NE-1);
Tyson Littenberg's avatar
Tyson Littenberg committed
276
277

  // Plus and minus templates for each detector:
Tyson Littenberg's avatar
Tyson Littenberg committed
278
279
280
281
282
  s     = double_matrix(NI-1,N-1);
  hP		= double_matrix(NI-1,N-1);
  hM		= double_matrix(NI-1,N-1);
  hPP		= double_matrix(NI-1,N-1);
  hMM		= double_matrix(NI-1,N-1);
Tyson Littenberg's avatar
Tyson Littenberg committed
283

Tyson Littenberg's avatar
Tyson Littenberg committed
284
  hdiv = double_tensor(NE-1,NI-1,N-1);
Tyson Littenberg's avatar
Tyson Littenberg committed
285

Tyson Littenberg's avatar
Tyson Littenberg committed
286
  epsilon = double_vector(NE-1);
Tyson Littenberg's avatar
Tyson Littenberg committed
287
288

  // Compute SNR of signal model
289
  combinePolarizations(data,model->signal,model->h,params,model->Npol);
290
  computeProjectionCoeffs(data, projectionP, params, data->fmin, data->fmax);
291
  waveformProject(data, projectionP, params, hP, model->h, data->fmin, data->fmax);
Tyson Littenberg's avatar
Tyson Littenberg committed
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314

  //store data
  for(i=0;i<NI;i++)
  {
    for(j=0;j<N;j++)
    {
      s[i][j]=data->s[i][j];
      data->s[i][j]=hP[i][j];
    }
  }

  double range;
  for(i = 0; i < NE; i++)
  {
    switch(i)
    {
      case 0:
        range = LAL_TWOPI;
        break;
      case 1:
        range = 2.0;
        break;
      case 2:
Tyson Littenberg's avatar
Tyson Littenberg committed
315
        range = LAL_PI_2;
Tyson Littenberg's avatar
Tyson Littenberg committed
316
317
318
319
320
        break;
      case 3:
        range = 0.99*2;
        break;
      case 4:
Tyson Littenberg's avatar
Tyson Littenberg committed
321
        range = LAL_TWOPI;
Tyson Littenberg's avatar
Tyson Littenberg committed
322
323
        break;
      case 5:
Tyson Littenberg's avatar
Tyson Littenberg committed
324
            //TODO:  HACK!  testing extrinsic amplitude update
325
        range = 2.0;
Tyson Littenberg's avatar
Tyson Littenberg committed
326
327
328
329
330
331
332
333
334
335
336
337
338
339
        break;
    }

    epsilon[i] = 1.0e-5*range;

    for(j = 0; j < NE; j++)
    {
      paramsP[j] = params[j];
      paramsM[j] = params[j];
    }

    paramsP[i] += epsilon[i];
    paramsM[i] -= epsilon[i];

340

341
342
    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);
Tyson Littenberg's avatar
Tyson Littenberg committed
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
    epsilon[i] = 0.1/sqrt(-(logLM+logLP)/(epsilon[i]*epsilon[i]));
    if(epsilon[i]!=epsilon[i])epsilon[i]=1.0e-5;
  }

  //restore data
  for(i=0;i<NI;i++)  for(j=0;j<N;j++) data->s[i][j]=s[i][j];

  /* assumes all the parameters are log or angles */
  for(i=0; i<NE; i++)
  {
    for(j=0; j<NE; j++)
    {
      paramsP[j] = params[j];
      paramsM[j] = params[j];
    }


    paramsP[i] += epsilon[i];
    paramsM[i] -= epsilon[i];

363
    combinePolarizations(data,model->signal,model->h,paramsP,model->Npol);
364
    computeProjectionCoeffs(data, projectionP, paramsP, data->fmin, data->fmax);
365
    waveformProject(data, projectionP, paramsP, hP, model->h, data->fmin, data->fmax);
Tyson Littenberg's avatar
Tyson Littenberg committed
366

367
368
369
    combinePolarizations(data,model->signal,model->h,paramsM,model->Npol);
    computeProjectionCoeffs(data, projectionM, paramsM, data->fmin, data->fmax);
    waveformProject(data, projectionM, paramsM, hM, model->h, data->fmin, data->fmax);
Tyson Littenberg's avatar
Tyson Littenberg committed
370
371
372
373

    paramsP[i] += epsilon[i];
    paramsM[i] -= epsilon[i];

374
375
376
377
378
379
380
    combinePolarizations(data,model->signal,model->h,paramsP,model->Npol);
    computeProjectionCoeffs(data, projectionP, paramsP, data->fmin, data->fmax);
    waveformProject(data, projectionP, paramsP, hPP, model->h, data->fmin, data->fmax);
    
    combinePolarizations(data,model->signal,model->h,paramsM,model->Npol);
    computeProjectionCoeffs(data, projectionM, paramsM, data->fmin, data->fmax);
    waveformProject(data, projectionM, paramsM, hMM, model->h, data->fmin, data->fmax);
Tyson Littenberg's avatar
Tyson Littenberg committed
381
382
383
384
385
386
387
388
389
390
391


    //Central differencing
    for(k=0; k<NI; k++)
    {
      for(j=0; j<N; j++)
      {
        hdiv[i][k][j] = (hMM[k][j] + 8.0*hP[k][j] - 8.0*hM[k][j] - hPP[k][j])/(12.0*epsilon[i]);
      }
    }
  }
392

Tyson Littenberg's avatar
Tyson Littenberg committed
393
  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);
394

Tyson Littenberg's avatar
Tyson Littenberg committed
395
396
  // stabilizers
  for(i=0; i<NE; i++)  fisher->matrix[i][i] += stab;
397
398


Tyson Littenberg's avatar
Tyson Littenberg committed
399
400
  free_double_vector(paramsP);
  free_double_vector(paramsM);
401

Tyson Littenberg's avatar
Tyson Littenberg committed
402
  free_double_tensor(hdiv,NE-1,NI-1);
403

Tyson Littenberg's avatar
Tyson Littenberg committed
404
405
406
407
408
  free_double_matrix(s,NI-1);
  free_double_matrix(hP,NI-1);
  free_double_matrix(hM,NI-1);
  free_double_matrix(hPP,NI-1);
  free_double_matrix(hMM,NI-1);
409

Tyson Littenberg's avatar
Tyson Littenberg committed
410
  free_double_vector(epsilon);
411

Tyson Littenberg's avatar
Tyson Littenberg committed
412
413
  free_network(projectionM,NI);
  free_network(projectionP,NI);
Tyson Littenberg's avatar
Tyson Littenberg committed
414
415
416
417
  free(projectionP);
  free(projectionM);

}
418

Tyson Littenberg's avatar
Tyson Littenberg committed
419
420
void extrinsic_fisher_update(struct Data *data, struct Model *model)
{
Tyson Littenberg's avatar
Tyson Littenberg committed
421
  int i,j;
422

Tyson Littenberg's avatar
Tyson Littenberg committed
423
424
  double *params = model->extParams;
  double *Sn     = model->Sn;
425

Tyson Littenberg's avatar
Tyson Littenberg committed
426
427
428
429
430
431
  double **glitch = malloc(data->NI*sizeof(double*));
  for(i=0; i<data->NI; i++)
  {
      glitch[i]=malloc(data->N*sizeof(double));
      for(j=0; j<data->N; j++)glitch[i][j] = model->glitch[i]->templates[j];
  }
432

Tyson Littenberg's avatar
Tyson Littenberg committed
433
  struct FisherMatrix *fisher = model->fisher;
434

Tyson Littenberg's avatar
Tyson Littenberg committed
435
  //calculate elements of Fisher
436
  extrinsic_fisher_information_matrix(fisher, params, model->invSnf, Sn, model, glitch, data);
437

438
  matrix_eigenstuff(fisher->matrix, fisher->evector, fisher->evalue, fisher->N);
439

Tyson Littenberg's avatar
Tyson Littenberg committed
440
  for(i=0; i<data->NI; i++) free(glitch[i]);
Tyson Littenberg's avatar
Tyson Littenberg committed
441
  free(glitch);
442

Tyson Littenberg's avatar
Tyson Littenberg committed
443
}
444

Tyson Littenberg's avatar
Tyson Littenberg committed
445
void fisher_matrix_proposal(struct FisherMatrix *fisher, double *params, gsl_rng *seed, double *y)
Tyson Littenberg's avatar
Tyson Littenberg committed
446
447
{
  int i, j;
448

449
450
451
452
  int N = fisher->N;

  double Amps[N];
  double jump[N];
453

Tyson Littenberg's avatar
Tyson Littenberg committed
454
455
    for(i=0; i<N; i++) Amps[i]=jump[i]=0.0;

Tyson Littenberg's avatar
Tyson Littenberg committed
456
  // not used since we jump along eigendirections
457
  double sqD = 1./sqrt((double)(N));
458

459
  //choose a single eigenvector to jump over
Tyson Littenberg's avatar
Tyson Littenberg committed
460
  i = (int)(uniform_draw(seed)*(double)N);
461

Tyson Littenberg's avatar
Tyson Littenberg committed
462
  //draw the eigen-jump amplitudes from N[0,1] scaled by evalue & dimension
463
464
//  for(i=0; i<N; i++)
//  {
465
    Amps[i] = 1.0/sqrt(fisher->evalue[i]);
Tyson Littenberg's avatar
Tyson Littenberg committed
466
    jump[i]	= 0.0;
467
//  }
468

Tyson Littenberg's avatar
Tyson Littenberg committed
469
  //decompose eigenjumps into parameter directions
470
471
  //for(i=0; i<N; i++)
  for (j=0; j<N; j++) jump[j] += Amps[i]*fisher->evector[j][i]*sqD;
472

Tyson Littenberg's avatar
Tyson Littenberg committed
473
  //jump from current position
Tyson Littenberg's avatar
Tyson Littenberg committed
474
  for(i=0; i<N; i++) y[i] = params[i] + gaussian_draw(seed)*jump[i];
Tyson Littenberg's avatar
Tyson Littenberg committed
475
}
476

477
void intrinsic_fisher_update(double *params, double *dparams, double *Snf, double Snx, double Tobs, int NW, int TFQFLAG)
Tyson Littenberg's avatar
Tyson Littenberg committed
478
479
480
481
{
  int i;
  double SNR;//,SNR2;
  i = (int)(params[1]*Tobs);
Meg Millhouse's avatar
Meg Millhouse committed
482
  double betasq;
483
484

  //SNR defined with Sn(f) but Snf array holdes <n_i^2>
Tyson Littenberg's avatar
Tyson Littenberg committed
485
  SNR = params[3]*sqrt((params[2]/(2.0*RT2PI*params[1]))/(Snx*Snf[i]*2.0/Tobs));
486

Tyson Littenberg's avatar
Tyson Littenberg committed
487
  if(SNR < 5.0) SNR = 5.0;  // this caps the size of the proposed jumps
488

Tyson Littenberg's avatar
Tyson Littenberg committed
489
490
  double sqrt3=1.73205080756888;
  double invSNR = 1./SNR;
Meg Millhouse's avatar
Meg Millhouse committed
491
492
493
  
  //double rtQ = sqrt(params[2]);
  betasq = 0.0;
494
495
  if (NW == 6) betasq = params[5]*params[5];

Meg Millhouse's avatar
Meg Millhouse committed
496
497
  double Qsq = params[2]*params[2];
  
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
  for(int k=0; k<NW; k++) dparams[k] = 0.0;

  // [0] t, [1] f, [2] Q, [3] A, [4] phi, [5] beta

  //always compute deltas for amplitude and phase
  dparams[3] = params[3]*invSNR; //amplitude
  dparams[4] = 1.0*invSNR;       //phase

  /*
   only update deltas for intrinsic parameters when TFQ flag is FALSE
   TFQ flag is TRUE when holding TFQ(beta) fixed, e.g. when updating hx
   */
  if(!TFQFLAG)
  {
    dparams[0] = params[2]/LAL_TWOPI/params[1]/sqrt(Qsq+1.0+LAL_PI_2*betasq)*invSNR; //time
    dparams[1] = 2.0*params[1]/sqrt(Qsq+3.0+3.0*LAL_PI_2*betasq)*invSNR;             //frequency
    dparams[2] = 2.0*params[2]/sqrt3/sqrt(1.0+LAL_PI_2*betasq)*invSNR;               //Q
    if (NW == 6) dparams[5] = 4.0/sqrt3/LAL_PI*invSNR;                               //beta
Meg Millhouse's avatar
Meg Millhouse committed
516
  }
Tyson Littenberg's avatar
Tyson Littenberg committed
517
518
}

Meg Millhouse's avatar
Meg Millhouse committed
519
void intrinsic_fisher_matrix(double *params, struct FisherMatrix *fisher, int NW)
520
{
521
522
523
524
  double **analytic_fisher = fisher->matrix;

  double t, f, Q, A, phi, beta = 0.0;
  //double f, Q, A, beta=0.0;
525
  int ii, jj, Nparams;
526
527

  //TODO: intrinsic_fisher_matrix has hard-coded D-wavelet!
528
  Nparams = NW;
529

Meg Millhouse's avatar
Meg Millhouse committed
530
  t = params[0];
531
532
533
  f = params[1];
  Q = params[2];
  A = params[3];
Meg Millhouse's avatar
Meg Millhouse committed
534
  phi = params[4];
535
536
537
538
539
540
541
542
543
544
545
546
  if(NW>5) beta = params[5];

  //precomput repeated constants
  double PI2 = LAL_PI*LAL_PI;
  double f2  = f*f;
  double Q2  = Q*Q;
  double b2  = beta*beta;
  double invf= 1./f;
  double invQ= 1./Q;
  double invA= 1./A;

  // Zero matrix out
547
548
  for (ii = 0; ii < Nparams; ii++) {
    for (jj = 0; jj < Nparams; jj++) {
549
      analytic_fisher[ii][jj] = 0.0;
550
551
    }
  }
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619

  // t t
  analytic_fisher[0][0] = (4.0*f2*PI2*(1.+Q2+PI2*b2))*invQ*invQ;

  // t f
  analytic_fisher[0][1] = (-2.0*PI2*beta);
  analytic_fisher[1][0] = analytic_fisher[0][1];

  // t Q
  analytic_fisher[0][2] = PI2*beta*f*invQ;
  analytic_fisher[2][0] = analytic_fisher[0][2];

  // t phi
  analytic_fisher[0][4] = (-2.0*f*LAL_PI);
  analytic_fisher[4][0] = analytic_fisher[0][4];

  // t beta
  analytic_fisher[0][5] = (-PI2*f)/2.0;
  analytic_fisher[5][0] = analytic_fisher[0][5];

  // f f
  analytic_fisher[1][1] = 0.25*((3.0+Q2+3.0*PI2*b2))*invf*invf;//((3.0+Q*Q+3.0*LAL_PI*LAL_PI*beta*beta)/(4.0*f*f));

  // f Q
  analytic_fisher[1][2] = 0.25*(-3.0*(1.0+PI2*b2))*invf*invQ;//(-3.0*(1.0+LAL_PI*LAL_PI*beta*beta)/(4.0*f*Q));
  analytic_fisher[2][1] = analytic_fisher[1][2];

  // f lnA
  analytic_fisher[1][3] = -0.5*invf*invA;//-(1.0/(2.0*f));
  analytic_fisher[3][1] = analytic_fisher[1][3];

  // f phi
  analytic_fisher[1][4] = 0.5*LAL_PI*beta*invf;//LAL_PI*beta/(2.0*f);
  analytic_fisher[4][1] = analytic_fisher[1][4];

  // f beta
  //analytic_fisher[1][5] = 3.0*beta*PI2/(8.0*f);
  //analytic_fisher[5][1] = analytic_fisher[1][5];

  // Q Q
  analytic_fisher[2][2] = 0.25*(3.0*(1.0+PI2*b2))*invQ*invQ;//(3.0*(1.0+LAL_PI*LAL_PI*beta*beta)/(4.0*Q*Q));

  // Q lnA
  analytic_fisher[2][3] = 0.5*invQ*invA;//(1.0/(2.0*Q));
  analytic_fisher[3][2] = analytic_fisher[2][3];

  // Q phi
  analytic_fisher[2][4] = -0.5*LAL_PI*beta*invQ;//-(LAL_PI*beta)/(2.0*Q);
  analytic_fisher[4][2] = analytic_fisher[2][4];

  // Q beta
  analytic_fisher[2][5] = -(3.0*beta*LAL_PI*LAL_PI)/(8.0*Q);
  analytic_fisher[5][2] = analytic_fisher[2][5];

  // lnA lnA
  analytic_fisher[3][3] = invA*invA;

  // phi phi
  analytic_fisher[4][4] = 1.0;

  // phi beta
  analytic_fisher[4][5] = 0.25*LAL_PI;//LAL_PI/4.0;
  analytic_fisher[5][4] = analytic_fisher[4][5];

  // beta beta
  analytic_fisher[5][5] = 3.0*LAL_PI*LAL_PI/16.0;

  // Get eigenvectors/values for Fisher matrix
Meg Millhouse's avatar
Meg Millhouse committed
620
621
  matrix_eigenstuff(fisher->matrix, fisher->evector, fisher->evalue, fisher->N);
  
622
623
}

Tyson Littenberg's avatar
Tyson Littenberg committed
624
void intrinsic_fisher_proposal(UNUSED struct Model *model, double **range, double *paramsx, double *paramsy, double *Snf, double Snx, double *logpx, double *logpy, gsl_rng *seed, double Tobs, int *test, int NW, int TFQFLAG)
Tyson Littenberg's avatar
Tyson Littenberg committed
625
626
627
628
629
{
  int k;
  double x,y;
  double qyx = *logpy;
  double qxy = *logpx;
Tyson Littenberg's avatar
Tyson Littenberg committed
630
631
  double *dparamsx = double_vector(NW-1);
  double *dparamsy = double_vector(NW-1);
Meg Millhouse's avatar
Meg Millhouse committed
632
  double scale = 1./sqrt((double)NW);//Number of intrinsic parameters
633
634
635
636
637

  // this is a crude diagonal Fisher matrix for individual sine-gaussians
  intrinsic_fisher_update(paramsx,dparamsx,Snf,Snx,Tobs, NW, TFQFLAG);

  // [0] t0 [1] f0 [2] Q [3] Amp [4] phi ([5] beta chirplet)
Meg Millhouse's avatar
Meg Millhouse committed
638
  
639
640
  //perturb current location by N[0,scale*sigma]
  for(k=0; k<NW; k++)
Tyson Littenberg's avatar
Tyson Littenberg committed
641
  {
Tyson Littenberg's avatar
Tyson Littenberg committed
642
    paramsy[k] = paramsx[k] + scale*dparamsx[k]*gaussian_draw(seed);
Tyson Littenberg's avatar
Tyson Littenberg committed
643
  }
644
645

  if(checkrange(paramsy, range, NW))
Tyson Littenberg's avatar
Tyson Littenberg committed
646
  {
Tyson Littenberg's avatar
Tyson Littenberg committed
647
648
    free_double_vector(dparamsx);
    free_double_vector(dparamsy);
649

Tyson Littenberg's avatar
Tyson Littenberg committed
650
651
652
    *test=1;
    return;
  }
653

Tyson Littenberg's avatar
Tyson Littenberg committed
654
  //compute probabiltiy of such a proposal (could just product the draws from gaussian_draw() above
Tyson Littenberg's avatar
Tyson Littenberg committed
655
  y = 1.0;
656
  for(k=0; k<NW; k++)
Tyson Littenberg's avatar
Tyson Littenberg committed
657
658
659
660
661
662
  {
    x = (paramsx[k]-paramsy[k])/(scale*dparamsx[k]);
    y *= scale*dparamsx[k];
    qyx +=  -0.5*x*x;
  }
  qyx -= log(y);
663
664
665
666
667

  //need fisher at new location to compute reverse probability
  intrinsic_fisher_update(paramsy,dparamsy,Snf,Snx,Tobs, NW, TFQFLAG);

  //compute probabiltiy of reverse proposal
Tyson Littenberg's avatar
Tyson Littenberg committed
668
  y = 1.0;
669
  for(k=0; k<NW; k++)
Tyson Littenberg's avatar
Tyson Littenberg committed
670
671
672
673
674
675
  {
    x = (paramsx[k]-paramsy[k])/(scale*dparamsy[k]);
    y *= scale*dparamsy[k];
    qxy +=  -0.5*x*x;
  }
  qxy -= log(y);
676
677
678
679
680
681
682
683
684

  // map phase onto [0,2pi]
  while(paramsy[4] > LAL_TWOPI) paramsy[4] -=LAL_TWOPI;
  while(paramsy[4] < 0.0)       paramsy[4] +=LAL_TWOPI;

  //update incoming pointers and return
  *logpy += qyx;
  *logpx += qxy;

Tyson Littenberg's avatar
Tyson Littenberg committed
685
686
  free_double_vector(dparamsx);
  free_double_vector(dparamsy);
687
}
688
689
690



Tyson Littenberg's avatar
Tyson Littenberg committed
691
692
693
694
695
/* ********************************************************************************** */
/*                                                                                    */
/*                       Custom intrinsic proposal distributions                      */
/*                                                                                    */
/* ********************************************************************************** */
696

Tyson Littenberg's avatar
Tyson Littenberg committed
697
double wavelet_proximity_density(double f, double t, double **params, int *larray, int cnt, struct Prior *prior)
698
{
Tyson Littenberg's avatar
Tyson Littenberg committed
699
700
701
702
703
704
  int jj, kk;
  double tau;
  double sigwt, signt, sigwf, signf;
  double sigwt2, signt2, sigwf2, signf2;
  double df, dt, df2, dt2;
  double pden;
705

Tyson Littenberg's avatar
Tyson Littenberg committed
706
707
  double **range = prior->range;
  double TFV = prior->TFV;
708

Tyson Littenberg's avatar
Tyson Littenberg committed
709
  // distribution we are drawing from is hollowed out 2-d Gaussian
710
711


Tyson Littenberg's avatar
Tyson Littenberg committed
712
713
  // we give equal weight to each wavelet (equally likely to draw from any of them)
  // so have to sum up densities for each wavelet and divide by the number of wavelets
714

Tyson Littenberg's avatar
Tyson Littenberg committed
715
  double *norm = NULL;
716
717


Tyson Littenberg's avatar
Tyson Littenberg committed
718
719
720
721
722
723
724
  // if no other active wavelets, the draw would have been from the full TF volume
  if(cnt == 0)
  {
    pden = 1.0/TFV;
  }
  else
  {
725

Tyson Littenberg's avatar
Tyson Littenberg committed
726
727
728
    //get normalization for prior
    norm = malloc( (prior->smax+1)*sizeof(double)  );
    proximity_normalization(params, range, norm, larray, cnt);
729
730


Tyson Littenberg's avatar
Tyson Littenberg committed
731
    pden = 0.0;
732

Tyson Littenberg's avatar
Tyson Littenberg committed
733
734
735
    // loop over the the existing wavelets
    for(jj=1; jj<= cnt; jj++)
    {
736

Tyson Littenberg's avatar
Tyson Littenberg committed
737
      kk = larray[jj];
738

Tyson Littenberg's avatar
Tyson Littenberg committed
739
740
      df = fabs(f-params[kk][1]);
      dt = fabs(t-params[kk][0]);
741

Tyson Littenberg's avatar
Tyson Littenberg committed
742
743
744
745
746
747
748
749
750
751
752
      tau = params[kk][2]/(LAL_TWOPI*params[kk][1]);

      sigwt = SPRD*tau;
      sigwf = SPRD/(LAL_PI*tau);

      signt = HLW*tau;
      signf = HLW/(LAL_PI*tau);


      if(dt < 4.0*sigwt && df < 4.0*sigwf)   // don't bother adding in contribution wavelets that are too far away
      {
753

Tyson Littenberg's avatar
Tyson Littenberg committed
754
755
        sigwt2 = sigwt*sigwt;
        signt2 = signt*signt;
756

Tyson Littenberg's avatar
Tyson Littenberg committed
757
758
        sigwf2 = sigwf*sigwf;
        signf2 = signf*signf;
759

Tyson Littenberg's avatar
Tyson Littenberg committed
760
761
        df2 = df*df;
        dt2 = dt*dt;
762

Tyson Littenberg's avatar
Tyson Littenberg committed
763
        pden += norm[kk]*(exp(-df2/(2.0*sigwf2)-dt2/(2.0*sigwt2))-exp(-df2/(2.0*signf2)-dt2/(2.0*signt2)));
764

Tyson Littenberg's avatar
Tyson Littenberg committed
765
766
767
768
769
770
771
772
773
774
775
      }

    }

    pden /= (double)(cnt);

    pden *= (1.0-UNIFORM);

    pden += UNIFORM/TFV;   // the uniform fraction

    free(norm);
776

Tyson Littenberg's avatar
Tyson Littenberg committed
777
778
779
  }

  return pden;
780
781
}

Tyson Littenberg's avatar
Tyson Littenberg committed
782
void wavelet_proximity(double *paramsy, double **params, int *larray, double **range, int cnt, gsl_rng *seed)
783
{
Tyson Littenberg's avatar
Tyson Littenberg committed
784
785
786
787
788
789
790
  int i, jj, kk;
  double alpha;
  double tau;
  double sigwt, signt, sigwf, signf;
  double sigwt2, signt2, sigwf2, signf2;
  double df, dt, df2, dt2;
  double rat=-1;
791

Tyson Littenberg's avatar
Tyson Littenberg committed
792
  alpha = uniform_draw(seed);
793

Tyson Littenberg's avatar
Tyson Littenberg committed
794
795
796
797
  if(alpha < UNIFORM || cnt == 0)  // 100% uniform if no other pixels lit
  {
    for(i=0; i< 2; i++)
    {
Tyson Littenberg's avatar
Tyson Littenberg committed
798
      paramsy[i] = range[i][0] + (range[i][1]-range[i][0])*uniform_draw(seed);
Tyson Littenberg's avatar
Tyson Littenberg committed
799
800
801
802
    }
  }
  else
  {
803

Tyson Littenberg's avatar
Tyson Littenberg committed
804
    // select one of the existing pixels and propose near it
Tyson Littenberg's avatar
Tyson Littenberg committed
805
    jj = (int)(uniform_draw(seed)*(double)(cnt))+1;
806

Tyson Littenberg's avatar
Tyson Littenberg committed
807
    kk = larray[jj];
808

Tyson Littenberg's avatar
Tyson Littenberg committed
809
    tau = params[kk][2]/(LAL_TWOPI*params[kk][1]);
810

Tyson Littenberg's avatar
Tyson Littenberg committed
811
    // scaling have been chosen to give distances of about 1 to 5 ish
812

Tyson Littenberg's avatar
Tyson Littenberg committed
813
814
815
816
    sigwt = SPRD*tau;
    signt = HLW*tau;
    sigwt2 = sigwt*sigwt;
    signt2 = signt*signt;
817

Tyson Littenberg's avatar
Tyson Littenberg committed
818
819
820
821
    sigwf = SPRD/(LAL_PI*tau);
    signf = HLW/(LAL_PI*tau);
    sigwf2 = sigwf*sigwf;
    signf2 = signf*signf;
822

Tyson Littenberg's avatar
Tyson Littenberg committed
823
    // distribution we are drawing from is hollowed out 2-D Gaussian
824

Tyson Littenberg's avatar
Tyson Littenberg committed
825
826
827
    // rejection sample to find t2, f2
    while(alpha < rat)
    {
Tyson Littenberg's avatar
Tyson Littenberg committed
828
      df = sigwf*gaussian_draw(seed);
Tyson Littenberg's avatar
Tyson Littenberg committed
829
      df2 = df*df;
830

Tyson Littenberg's avatar
Tyson Littenberg committed
831
      dt = sigwt*gaussian_draw(seed);
Tyson Littenberg's avatar
Tyson Littenberg committed
832
      dt2 = dt*dt;
833

Tyson Littenberg's avatar
Tyson Littenberg committed
834
      rat = (exp(-df2/(2.0*sigwf2)-dt2/(2.0*sigwt2))-exp(-df2/(2.0*signf2)-dt2/(2.0*signt2)));
835

Tyson Littenberg's avatar
Tyson Littenberg committed
836
      alpha = uniform_draw(seed);
Tyson Littenberg's avatar
Tyson Littenberg committed
837
838
839

      paramsy[1] = params[kk][1]+df;
      paramsy[0] = params[kk][0]+dt;
840

Tyson Littenberg's avatar
Tyson Littenberg committed
841
842
      if (paramsy[0]<range[0][0] || paramsy[0]>range[0][1]) rat = 0.0;
      if (paramsy[1]<range[1][0] || paramsy[1]>range[1][1]) rat = 0.0;
843

Tyson Littenberg's avatar
Tyson Littenberg committed
844
    }
845
846


Tyson Littenberg's avatar
Tyson Littenberg committed
847
848
  }
  return;
849
850
}

Tyson Littenberg's avatar
Tyson Littenberg committed
851
void wavelet_proximity_new(double *paramsy, double **params, int *larray, double **range, int cnt, gsl_rng *seed)
852
{
Tyson Littenberg's avatar
Tyson Littenberg committed
853
854
855
856
857
858
859
  int i, jj, kk;
  double alpha=1.0;
  double tau;
  double sigwt, signt, sigwf, signf;
  double sigwt2, signt2, sigwf2, signf2;
  double df, dt, df2, dt2;
  double rat=-1.0;
860
861


Tyson Littenberg's avatar
Tyson Littenberg committed
862
  alpha = uniform_draw(seed);
863

Tyson Littenberg's avatar
Tyson Littenberg committed
864
865
866
867
  if(alpha < UNIFORM || cnt == 0)  // 100% uniform if no other pixels lit
  {
    for(i=0; i< 2; i++)
    {
Tyson Littenberg's avatar
Tyson Littenberg committed
868
      paramsy[i] = range[i][0] + (range[i][1]-range[i][0])*uniform_draw(seed);
Tyson Littenberg's avatar
Tyson Littenberg committed
869
870
871
872
    }
  }
  else
  {
873

Tyson Littenberg's avatar
Tyson Littenberg committed
874
    // select one of the existing pixels and propose near it
Tyson Littenberg's avatar
Tyson Littenberg committed
875
    jj = (int)(uniform_draw(seed)*(double)(cnt))+1;
876

Tyson Littenberg's avatar
Tyson Littenberg committed
877
    kk = larray[jj];
878

Tyson Littenberg's avatar
Tyson Littenberg committed
879
    tau = params[kk][2]/(LAL_TWOPI*params[kk][1]);
880

Tyson Littenberg's avatar
Tyson Littenberg committed
881
    // scaling have been chosen to give distances of about 1 to 5 ish
882

Tyson Littenberg's avatar
Tyson Littenberg committed
883
884
885
886
    sigwt = SPRD*tau;
    signt = HLW*tau;
    sigwt2 = sigwt*sigwt;
    signt2 = signt*signt;
887

Tyson Littenberg's avatar
Tyson Littenberg committed
888
889
890
891
    sigwf = SPRD/(LAL_PI*tau);
    signf = HLW/(LAL_PI*tau);
    sigwf2 = sigwf*sigwf;
    signf2 = signf*signf;
892

Tyson Littenberg's avatar
Tyson Littenberg committed
893
894
    // distribution we are drawing from is hollowed out Gaussian
    // 1/(sqrt(2 Pi)(sigmaW-sigmaN))*(exp(-x^2/2 sigmaW^2) - exp(-x^2/2 sigmaN^2))
895

Tyson Littenberg's avatar
Tyson Littenberg committed
896
    // rejection sample to find t2, f2
897
898
899
    //TODO: Why were we rejection sampling??
//    while(alpha > rat)
//    {
Tyson Littenberg's avatar
Tyson Littenberg committed
900
      df = sigwf*gaussian_draw(seed);
Tyson Littenberg's avatar
Tyson Littenberg committed
901
      df2 = df*df;
902

Tyson Littenberg's avatar
Tyson Littenberg committed
903
      dt = sigwt*gaussian_draw(seed);
Tyson Littenberg's avatar
Tyson Littenberg committed
904
      dt2 = dt*dt;
905

Tyson Littenberg's avatar
Tyson Littenberg committed
906
      rat = 1.0-exp(df2*(1.0/(2.0*sigwf2) - 1.0/(2.0*signf2)))*exp(dt2*(1.0/(2.0*sigwt2) - 1.0/(2.0*signt2)));
907

Tyson Littenberg's avatar
Tyson Littenberg committed
908
      alpha = uniform_draw(seed);
909

Tyson Littenberg's avatar
Tyson Littenberg committed
910
911
      paramsy[1] = params[kk][1]+df;
      paramsy[0] = params[kk][0]+dt;
912

Tyson Littenberg's avatar
Tyson Littenberg committed
913
914
915
916
      if
        (paramsy[0]<range[0][0] || paramsy[0]>range[0][1]) rat = -1.0;
      else if
        (paramsy[1]<range[1][0] || paramsy[1]>range[1][1]) rat = -1.0;
917

918
//    }
Tyson Littenberg's avatar
Tyson Littenberg committed
919
920
  }
}
921

Tyson Littenberg's avatar
Tyson Littenberg committed
922
void wavelet_sample(double *params, double **range, double **prop, int tsize, double Tobs, double pmax, gsl_rng *seed)
Tyson Littenberg's avatar
Tyson Littenberg committed
923
924
925
926
{
  int i, j;
  double f, t, alpha;
  double df, dt;
927

Tyson Littenberg's avatar
Tyson Littenberg committed
928
929
  // Note: it is possible to speed the rejection sampling up by separating out the high density pixels
  // and the uniform piece.
930

Tyson Littenberg's avatar
Tyson Littenberg committed
931
932
  dt = Tobs/(double)tsize;
  df = 0.5/dt;
933

Tyson Littenberg's avatar
Tyson Littenberg committed
934
935
  f = range[1][0]+(range[1][1]-range[1][0])*uniform_draw(seed);
  t = range[0][0]+(range[0][1]-range[0][0])*uniform_draw(seed);
936

Tyson Littenberg's avatar
Tyson Littenberg committed
937
938
  i = (int)(f/df);
  j = (int)(t/dt);
939

Tyson Littenberg's avatar
Tyson Littenberg committed
940
  alpha = 2.*pmax*uniform_draw(seed);
941

Tyson Littenberg's avatar
Tyson Littenberg committed
942
943
  while(alpha > prop[i][j])
  {
Tyson Littenberg's avatar
Tyson Littenberg committed
944
945
    f = range[1][0]+(range[1][1]-range[1][0])*uniform_draw(seed);
    t = range[0][0]+(range[0][1]-range[0][0])*uniform_draw(seed);
946

Tyson Littenberg's avatar
Tyson Littenberg committed
947
948
    i = (int)(f/df);
    j = (int)(t/dt);
949

Tyson Littenberg's avatar
Tyson Littenberg committed
950
    alpha = 2.*pmax*uniform_draw(seed);
Tyson Littenberg's avatar
Tyson Littenberg committed
951
  }
952

Tyson Littenberg's avatar
Tyson Littenberg committed
953
954
  params[0] = t;
  params[1] = f;
955

Tyson Littenberg's avatar
Tyson Littenberg committed
956
957
  return;
}
958

Tyson Littenberg's avatar
Tyson Littenberg committed
959
void draw_wavelet_params(double *params, double **range, gsl_rng *seed, int NW)
Tyson Littenberg's avatar
Tyson Littenberg committed
960
961
962
963
{
  int i;

  // Draw the rest from uniform distribution
Meg Millhouse's avatar
Meg Millhouse committed
964
  for(i=0; i<NW; i++)
Tyson Littenberg's avatar
Tyson Littenberg committed
965
  {
Tyson Littenberg's avatar
Tyson Littenberg committed
966
    params[i] = range[i][0] + (range[i][1]-range[i][0])*uniform_draw(seed);
Tyson Littenberg's avatar
Tyson Littenberg committed
967
  }
968

Tyson Littenberg's avatar
Tyson Littenberg committed
969
  return;
970
971
}

Tyson Littenberg's avatar
Tyson Littenberg committed
972
void draw_time_frequency(int ifo, int ii, struct Wavelet *wx, struct Wavelet *wy, double **range, gsl_rng *seed, double fr, struct TimeFrequencyMap *tf, int *prop)
973
{
Tyson Littenberg's avatar
Tyson Littenberg committed
974
  if(uniform_draw(seed)<fr)
975
976
977
978
  {
    TFprop_draw(wy->intParams[ii], range, tf, ifo, seed);
    *prop=0;
  }
Tyson Littenberg's avatar
Tyson Littenberg committed
979
  else