Cheetah - SKA - PSS - Prototype Time Domain Search Pipeline
PulsarInjection.cpp
1 #include "cheetah/generators/PulsarInjection.h"
2 #include "cheetah/generators/pulse_profile/ProfileManager.h"
3 #include "cheetah/utils/ConvolvePlan.h"
4 #include "panda/Log.h"
5 #include <boost/math/special_functions/sinc.hpp>
6 #include <exception>
7 
8 #define M_PI_DBL 3.14159265358979323846264338327950288
9 
10 namespace {
11 
12 template<typename T>
13 class FftwAlloc {
14  public:
15  FftwAlloc(std::size_t size)
16  {
17  _ptr = static_cast<T*>( fftwf_malloc(size * sizeof(T)) );
18  }
19  operator T*() { return _ptr; }
20 
21  ~FftwAlloc()
22  {
23  fftwf_free(_ptr);
24  }
25  private:
26  T* _ptr;
27 };
28 
29 // constants for random number generator
30 static constexpr std::size_t IM1(2147483563);
31 static constexpr std::size_t IM2(2147483399);
32 static constexpr std::size_t IMM1(IM1 - 1);
33 static constexpr float AM(1.0/IM1);
34 static constexpr int IA1(40014);
35 static constexpr int IA2(40692);
36 static constexpr int IQ1(53668);
37 static constexpr int IQ2(52774);
38 static constexpr int IR1(12211);
39 static constexpr int IR2(3791);
40 static constexpr int NTAB(32);
41 static constexpr std::size_t NDIV(1 + IMM1/NTAB);
42 static constexpr float EPS(1.2e-7);
43 #define RNMX (1.0 - EPS)
44 #define MJK_RAND_R1279_SZ 2048
45 #define MJK_RAND_R1279_SZ1 2047
46 #define TWO_PI 6.2831853071795864769252866
47 
48 typedef struct MjkRand {
49  public:
50  MjkRand(uint64_t init_seed)
51  : gauss_buffer(nullptr)
52  {
53  long nseed = -(long) init_seed;
54  next = -1;
55  nthreadmax = 8;
56  buffer_len = 1024 * 16 * nthreadmax;
57  buffer = (double*) calloc(buffer_len, sizeof(double));
58  rmax = INT32_MAX/2;
59  gauss_next = -1;
60  gauss_buffer = NULL;
61  srand(init_seed);
62  ir = (int*) calloc(nthreadmax, sizeof(int));
63  seed = (uint64_t*) calloc(MJK_RAND_R1279_SZ * nthreadmax, sizeof(uint64_t));
64  uint64_t j, k;
65  for(std::size_t i = 0; i < (unsigned) (MJK_RAND_R1279_SZ * nthreadmax); ++i) {
66  j = (uint64_t) (UINT32_MAX * (double) nrran2(&nseed));
67  j = j << 8 * sizeof(uint32_t);
68  k = (uint64_t) (UINT32_MAX * (double) nrran2(&nseed));
69  seed[i] = j + k;
70  }
71  rand_fill();
72  }
73 
74  ~MjkRand()
75  {
76  free(buffer);
77  free(ir);
78  free(seed);
79  if (gauss_buffer != nullptr) free(gauss_buffer);
80  }
81 
82  void rand_fill()
83  {
84  uint64_t ibuf, iblk, blksize = buffer_len/nthreadmax;
85  for (iblk = 0; iblk < nthreadmax; iblk++) {
86  for (ibuf = 0; ibuf < blksize; ibuf++) {
87  buffer[ibuf + iblk * blksize] = a2281(seed + iblk * MJK_RAND_R1279_SZ, &(ir)[iblk]);
88  }
89  }
90  next = buffer_len - 1;
91  }
92 
93  void fill_gauss() {
94  rand_fill();
95  double o1, o2, i1, i2;
96  // const uint32_t chunk = 64;
97  if (gauss_buffer == NULL) {
98  gauss_buffer = (float*) calloc(buffer_len, sizeof(float));
99  }
100  //#pragma omp parallel private(i,o1,o2,i1,i2)
101  //#pragma omp for schedule(dynamic,chunk)
102  for(std::size_t i = 0; i < buffer_len; i += 2) {
103  i1 = buffer[i];
104  i2 = buffer[i + 1];
105  rand_gauss2(i1, i2, &o1, &o2);
106  gauss_buffer[i] = o1;
107  gauss_buffer[i + 1] = o2;
108  }
109  gauss_next = buffer_len - 1;
110  next = -1;
111  }
112 
113  double rand_double()
114  {
115  if (next < 0) {
116  rand_fill();
117  }
118  return buffer[next--];
119  }
120 
121  float rand_gauss()
122  {
123  if (gauss_next < 0) {
124  fill_gauss();
125  }
126  return gauss_buffer[gauss_next--];
127  }
128 
129  protected:
130  float nrran2(long *idum)
131  {
132  int j;
133  long k;
134  static long idum2 = 123456789;
135  static long iy = 0;
136  static long iv[NTAB];
137  float temp;
138 
139  if (*idum <= 0) {
140  if (-(*idum) < 1) *idum = 1;
141  else *idum = -(*idum);
142  idum2 = (*idum);
143  for (j = NTAB + 7; j >= 0; j--) {
144  k = (*idum)/IQ1;
145  *idum = IA1 * (*idum - k * IQ1) - k * IR1;
146  if (*idum < 0) *idum += IM1;
147  if (j < NTAB) iv[j] = *idum;
148  }
149  iy = iv[0];
150  }
151  k = (*idum)/IQ1;
152  *idum = IA1 * (*idum - k * IQ1) - k * IR1;
153  if (*idum < 0) *idum += IM1;
154  k = idum2/IQ2;
155  idum2 = IA2 * (idum2 - k * IQ2) - k * IR2;
156  if (idum2 < 0) idum2 += IM2;
157  j = iy/NDIV;
158  iy = iv[j] - idum2;
159  iv[j] = *idum;
160  if (iy < 1) iy += IMM1;
161  temp = AM * iy;
162  if (temp > RNMX) return RNMX;
163  else return temp;
164  }
165 
166  double a2281(uint64_t seed[MJK_RAND_R1279_SZ], int *i)
167  {
168  const int i1 = (*i - 2281) & MJK_RAND_R1279_SZ1;
169  const int i2 = (*i - 1252) & MJK_RAND_R1279_SZ1;
170  seed[*i] = seed[i1] +seed[i2];
171  double ret = seed[*i]/(double) UINT64_MAX;
172  (*i) = (*i + 1)&MJK_RAND_R1279_SZ1;
173  return ret;
174  }
175 
176  void rand_gauss2(double rand1, double rand2, double *o1, double *o2)
177  {
178  if (rand1 < 1e-100) rand1 = 1e-100;
179  rand1 = -2 * log(rand1);
180  rand2 = rand2 * TWO_PI;
181  *o1 = sqrt(rand1) * cos(rand2);
182  *o2 = sqrt(rand1) * sin(rand2);
183  }
184 
185  private:
186  uint32_t nthreadmax;
187  uint32_t buffer_len;
188  double *buffer;
189  int32_t next;
190  uint32_t rmax;
191  int32_t gauss_next;
192  float* gauss_buffer;
193  uint64_t* seed;
194  int32_t* ir;
195 } mjk_rand_t;
196 
197 
198 
199 inline void convolve(ska::cheetah::utils::ConvolvePlan& plan)
200 {
201  plan.convolve();
202 }
203 
204 struct IsmModels {
205  IsmModels(std::size_t number_of_models, std::size_t model_size)
206  : _size(number_of_models)
207  {
208  _models = (float**) malloc(sizeof(float*) * _size);
209  if(! _models) throw std::bad_alloc();
210  for(std::size_t i = 0; i < _size; ++i) {
211  _models[i] = (float*) malloc(sizeof(float) * model_size);
212  if(! _models[i]) {
213  clean_models(i);
214  throw std::bad_alloc();
215  }
216  }
217  }
218 
219  ~IsmModels()
220  {
221  clean_models(_size);
222  delete _models;
223  }
224 
225  inline
226  float* operator[](std::size_t i) {
227  return _models[i];
228  }
229 
230  private:
231  inline
232  void clean_models(std::size_t model_num) {
233  for(std::size_t i = 0; i < model_num; ++i) {
234  free(_models[i]);
235  }
236  }
237 
238  private:
239  std::size_t _size;
240  float** _models;
241 };
242 
243 } //namespace
244 
245 namespace ska {
246 namespace cheetah {
247 namespace generators {
248 
249 
250 template<typename T>
251 PulsarInjection<T>::PulsarInjection(PulsarInjectionConfig const& config, ProfileManager const& profile_manager)
252  : BaseT()
253  , _config(config)
254  , _profile_manager(profile_manager)
255  , _phase_model(config.phase_model())
256  , _last_number_of_channels(0)
257  , _last_number_of_spectra(0)
258 {
259  // these calls will throw if not found
260  _profile = profile_manager.profile(config.profile());
261 }
262 
263 template<typename T>
264 PulsarInjection<T>::~PulsarInjection()
265 {
266 }
267 
268 template<typename T>
269 void PulsarInjection<T>::next(DataType& data)
270 {
271  float spec_index = _config.spectral_index();
272  const float in_snr = _config.signal_to_noise(); // S/N
273  float pulse_energy_sigma = 0.2;
274  float frac;
275  uint32_t seed = time(NULL);
276  uint_fast32_t nsubpulse = 5;
277  uint_fast32_t n;
278 
279  // fftwf_init_threads();
280  // fftwf_plan_with_nthreads(1);
281 
282  const uint_fast32_t nchan_const = data.number_of_channels();;
283  const uint64_t nsamp = data.number_of_spectra();;
284  long double tsamp = static_cast<boost::units::quantity<boost::units::si::time, double>>(data.sample_interval()).value(); // ensure we get value in seconds
285  if( data.channel_frequencies().size() <=1 ) return;
286  long double const fch1 = static_cast<boost::units::quantity<data::MegaHertz, double>>(data.channel_frequencies()[0]).value();
287  utils::ModifiedJulianClock::time_point mjd(static_cast<cheetah::utils::ModifiedJulianClock::time_point>(data.start_time()));
288  utils::ModifiedJulianClock::duration const tsamp_mjd(tsamp/86400.0L);
289 
290  // TODO scale to the RMS of the existing data
291  const double noiseAmp = 24; // 8-bit only to match gaussian noise generated by the fake program
292 
293  // normalise the profile (if it hasn't been done already)
294  if(nsamp != _last_number_of_spectra || nchan_const != _last_number_of_channels)
295  {
296  // restore profile from original template and rescale
297  PANDA_LOG_DEBUG << "rescaling profile '" << static_cast<std::string>(_config.profile()) << "'";
298  _profile = _profile_manager.profile(_config.profile());
299  double sum = _profile.sum();
300  sum /= (double) _profile.size();
301  if(_profile.size() == 0) {
302  panda::Error e("profile is empty:");
303  e << static_cast<std::string>(_config.profile());
304  throw e;
305  }
306  double scale = noiseAmp * in_snr/sqrt(nchan_const * nsamp)/sum;
307  PANDA_LOG_DEBUG << "mean = " << sum << ", scale = " << scale;
308  // normalise to pseudo-S/N in 1s
309  std::for_each(_profile.begin(), _profile.end(), [scale](PulsarProfile::ProfileDataPoint& dp) { dp *= scale; });
310  _last_number_of_spectra = nsamp;
311  _last_number_of_channels = nchan_const;
312  }
313  const uint_fast32_t nprof = _profile.size();
314 
315  FftwAlloc<float> ism_convolution(nprof);
316  FftwAlloc<float> subpulse_map(nprof); // mask for the number of pulses actually present (nsubpulse out of n randomly selected)
317  FftwAlloc<float> subpulse_profile(nprof);
318  FftwAlloc<float> unsmeared_profile(nprof);
319  FftwAlloc<float> smeared_profile(nprof);
320 
321  utils::ConvolvePlan ism_convolution_plan(nprof, unsmeared_profile, ism_convolution, smeared_profile);
322  utils::ConvolvePlan subpulse_convolution_plan(nprof, subpulse_profile, subpulse_map, unsmeared_profile);
323 
324  for(std::size_t i = 0; i < nprof; ++i){
325  subpulse_profile[i] = 1.0;
326  }
327 
328  double sum = 0.0;
329  // normalise the subpulse profile
330  for(std::size_t i = 0; i < nprof; ++i) {
331  sum += subpulse_profile[i];
332  }
333  sum /= (double) nprof;
334  // we also normalise by the number of subpulses so that the final subpulse profile will have an area of 1.
335  double scale = 1.0/sum/(float) nsubpulse;
336  for(std::size_t i = 0; i < nprof; ++i) {
337  subpulse_profile[i] = subpulse_profile[i] * scale;
338  }
339 
340  long double phase;
341  int_fast32_t pbin;
342  uint_fast32_t ch = 0;
343  std::vector<boost::units::quantity<data::MegaHertz, double>> freq(nchan_const + 1);
344  std::vector<float> sidx(nchan_const);
345  std::vector<uint_fast32_t> ism_idx(nchan_const);
346  IsmModels ism_models(nchan_const, nprof);
347  uint_fast32_t nism = 0;
348  const uint_fast32_t NCOPY = nprof * sizeof(float);
349  std::vector<double> poff(nchan_const);
350 
351  const double pulse_energy_norm = exp(pow(pulse_energy_sigma, 2)/2.0);
352  double p0 = _phase_model(mjd, data.channel_frequencies()[0]);
353 
354  std::vector<int_fast32_t> prevbin(nchan_const);
355  int_fast32_t dbin;
356 
357  MjkRand rnd(seed);
358  for (ch = 0; ch < nchan_const; ++ch) {
359  freq[ch] = data.channel_frequencies()[ch];
360  poff[ch] = _phase_model(mjd, freq[ch]) - p0;
361  prevbin[ch] = -1;
362  sidx[ch] = pow(freq[ch].value()/fch1, spec_index);
363  }
364  freq[nchan_const] = freq[nchan_const -1] + (freq[nchan_const - 1] - freq[nchan_const -2]); // assume a const bin width
365 
366  // set up the ISM model
367  {
368  int_fast32_t dm_bin=-1;
369  int_fast32_t scatter_bin=-1;
370  float* dm_conv = ism_convolution; // re-use an existing convolution function
371  float* scatt_conv = unsmeared_profile;
372  float* out_conv = smeared_profile;
373 
374  double dm = _config.dm().value();
375  double const freq0 = boost::units::quantity<data::MegaHertz, double>(freq[0]).value();
376  auto phase0 = _phase_model(mjd, freq[0]) + dm/(freq0*freq0);
377  for (ch = 0; ch < nchan_const; ++ch) {
378  //phase = T2Predictor_GetPhase(&pred,mjd,freq[ch]) - T2Predictor_GetPhase(&pred,mjd,freq[ch+1]);
379  double const freq_ch = boost::units::quantity<data::MegaHertz, double>(freq[ch+1]).value();
380  auto phase1 = _phase_model(mjd, freq[ch+1]) + dm/(freq_ch*freq_ch);
381  phase = phase0 - phase1;
382  phase0 = phase1;
383  pbin = std::abs(phase) * nprof;
384  if (pbin < 1) pbin = 1;
385  else if (pbin%2 == 0) pbin -= 1;
386 
387  if ((pbin != dm_bin) != scatter_bin) {
388  // make new ISM model
389  sum = 0;
390  for(std::size_t i = 0; i < nprof; ++i) {
391  // create a dispersed sinc function profile modelling the telescope response function
392  // to use in convolution.
393  // The function sinc(l) == sin(pi*l)/(pi*l) is the Fourier transform of the unit rectangle function
394  // and is the electric-field pattern of a uniformly illuminated unit aperture.
395  int_fast32_t x = (i - pbin/2 + nprof)%nprof;
396  double y = 2. * M_PI_DBL * (i - pbin/2.0 + 0.5)/pbin;
397  if (i < (std::size_t)pbin) {
398  dm_conv[x] = boost::math::sinc_pi(y);
399  sum += dm_conv[x];
400  }
401  else dm_conv[x] = 0;
402  }
403  for(std::size_t i = 0; i < nprof; ++i){
404  dm_conv[i] /= sum;
405  }
406 
407  scatt_conv[0] = 1.0;
408  for(std::size_t i = 1; i < nprof; ++i){
409  scatt_conv[i] = 0;
410  }
411 
412  convolve(ism_convolution_plan); // dm_conv * scatt_conv => out_conv
413 
414  memcpy(ism_models[nism], out_conv, NCOPY);
415  ++nism;
416  dm_bin = pbin;
417  scatter_bin = 0;
418  }
419  ism_idx[ch] = nism - 1;
420  }
421  }
422  PANDA_LOG_DEBUG << "Generated " << nism << " ISM models";
423 
424  std::vector<std::vector<float>> grads(nism, std::vector<float>(nprof));
425  std::vector<std::vector<float>> output_prof(nism, std::vector<float>(nprof));
426 
427  PANDA_LOG_DEBUG << "Starting simulation";
428  uint64_t count = 0;
429  int64_t ip0;
430  int64_t prev_ip0 = INT64_MAX;
431  std::vector<float>rands(nchan_const);
432  for(std::size_t sample_index = 0; sample_index < nsamp; ++sample_index)
433  {
434  auto spectra = data.spectrum(sample_index);
435  auto channel_it = spectra.begin();
436 
437  for (ch = 0; ch < nchan_const; ++ch) {
438  rands[ch] = rnd.rand_double();
439  }
440  p0 = _phase_model(mjd, freq[0]);
441  ip0 = (int64_t) floor(p0);
442  //pragma omp paralell private(i,n,ch,phase,frac,pbin,A,dbin) shared(ip0,prev_ip0)
443  {
444  if (ip0 != prev_ip0) {
445  for(std::size_t i = 0; i < nprof; i++) {
446  subpulse_map[i] = 0;
447  }
448  if(nsubpulse==0) {
449  // take the profile as given
450  std::copy(_profile.begin(), _profile.end(), &unsmeared_profile[0]);
451  }
452  else {
453  // generate subpulses
454  for (n = 0; n < nsubpulse; ++n) {
455  // subpulse_map - emulate that only parts of the pulse are actually seen(<=nsubpulse)
456  std::size_t i = floor(rnd.rand_double() * nprof);
457  subpulse_map[i] += exp(rnd.rand_gauss() * pulse_energy_sigma)/pulse_energy_norm;
458  }
459  // this convolution gives us a masked profile (i.e. a profile with random bits missing)
460  // with the result in unsmeared_profile
461  convolve(subpulse_convolution_plan);
462 
463  // use this mask to modulate the profile and define our unsmeared_profile
464  std::transform(_profile.begin(), _profile.end(), &unsmeared_profile[0], &unsmeared_profile[0], std::multiplies<float>());
465  }
466 
467  //pragma omp single
468  {
469  for(std::size_t i = 0; i < nism; ++i) {
470  memcpy(static_cast<float*>(ism_convolution), ism_models[i], NCOPY);
471  convolve(ism_convolution_plan);
472  memcpy(static_cast<void*>(&output_prof[i][0]), static_cast<void*>(smeared_profile), NCOPY);
473  for (n = 0; n < nprof - 1; ++n) {
474  grads[i][n] = smeared_profile[n + 1] - smeared_profile[n];
475  }
476  grads[i][nprof - 1] = smeared_profile[0] - smeared_profile[nprof - 1];
477  }
478  prev_ip0 = ip0;
479  }
480  }
481  //pragma omp for schedule(dynamic,128)
482  for (ch = 0; ch < nchan_const; ++ch) {
483  phase = p0 + poff[ch];
484  phase = phase - floor(phase);
485  pbin = floor(phase * nprof);
486  frac = phase * nprof - pbin;
487  float final_signal = output_prof[ism_idx[ch]][pbin] + frac * grads[ism_idx[ch]][pbin];
488  if (prevbin[ch] < 0) prevbin[ch] = pbin;
489  dbin = pbin - prevbin[ch];
490  while (dbin < 0) dbin += nprof;
491  while (dbin > 1) {
492  final_signal += output_prof[ism_idx[ch]][prevbin[ch]];
493  prevbin[ch] += 1;
494  dbin = pbin - prevbin[ch];
495  while (dbin < 0) dbin += nprof;
496  }
497 
498  if (final_signal > 0) {
499  *channel_it += floor(sidx[ch] * final_signal + rands[ch]);
500  }
501  ++channel_it;
502  prevbin[ch] = pbin;
503  }
504  }
505 
506  mjd += tsamp_mjd;
507  count++;
508  }
509  PANDA_LOG_DEBUG << "Simulation ended.";
510 }
511 
512 } // namespace generators
513 } // namespace cheetah
514 } // namespace ska
Some limits and constants for FLDO.
Definition: Brdz.h:35
Perform a convolution using fft methods.
Definition: ConvolvePlan.h:39