1 #include "cheetah/generators/PulsarInjection.h" 2 #include "cheetah/generators/pulse_profile/ProfileManager.h" 3 #include "cheetah/utils/ConvolvePlan.h" 5 #include <boost/math/special_functions/sinc.hpp> 8 #define M_PI_DBL 3.14159265358979323846264338327950288 15 FftwAlloc(std::size_t size)
17 _ptr =
static_cast<T*
>( fftwf_malloc(size *
sizeof(T)) );
19 operator T*() {
return _ptr; }
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 48 typedef struct MjkRand {
50 MjkRand(uint64_t init_seed)
51 : gauss_buffer(nullptr)
53 long nseed = -(long) init_seed;
56 buffer_len = 1024 * 16 * nthreadmax;
57 buffer = (
double*) calloc(buffer_len,
sizeof(
double));
62 ir = (
int*) calloc(nthreadmax,
sizeof(
int));
63 seed = (uint64_t*) calloc(MJK_RAND_R1279_SZ * nthreadmax,
sizeof(uint64_t));
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));
79 if (gauss_buffer !=
nullptr) free(gauss_buffer);
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]);
90 next = buffer_len - 1;
95 double o1, o2, i1, i2;
97 if (gauss_buffer == NULL) {
98 gauss_buffer = (
float*) calloc(buffer_len,
sizeof(
float));
102 for(std::size_t i = 0; i < buffer_len; i += 2) {
105 rand_gauss2(i1, i2, &o1, &o2);
106 gauss_buffer[i] = o1;
107 gauss_buffer[i + 1] = o2;
109 gauss_next = buffer_len - 1;
118 return buffer[next--];
123 if (gauss_next < 0) {
126 return gauss_buffer[gauss_next--];
130 float nrran2(
long *idum)
134 static long idum2 = 123456789;
136 static long iv[NTAB];
140 if (-(*idum) < 1) *idum = 1;
141 else *idum = -(*idum);
143 for (j = NTAB + 7; j >= 0; j--) {
145 *idum = IA1 * (*idum - k * IQ1) - k * IR1;
146 if (*idum < 0) *idum += IM1;
147 if (j < NTAB) iv[j] = *idum;
152 *idum = IA1 * (*idum - k * IQ1) - k * IR1;
153 if (*idum < 0) *idum += IM1;
155 idum2 = IA2 * (idum2 - k * IQ2) - k * IR2;
156 if (idum2 < 0) idum2 += IM2;
160 if (iy < 1) iy += IMM1;
162 if (temp > RNMX)
return RNMX;
166 double a2281(uint64_t seed[MJK_RAND_R1279_SZ],
int *i)
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;
176 void rand_gauss2(
double rand1,
double rand2,
double *o1,
double *o2)
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);
205 IsmModels(std::size_t number_of_models, std::size_t model_size)
206 : _size(number_of_models)
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);
214 throw std::bad_alloc();
226 float* operator[](std::size_t i) {
232 void clean_models(std::size_t model_num) {
233 for(std::size_t i = 0; i < model_num; ++i) {
247 namespace generators {
251 PulsarInjection<T>::PulsarInjection(PulsarInjectionConfig
const& config, ProfileManager
const& profile_manager)
254 , _profile_manager(profile_manager)
255 , _phase_model(config.phase_model())
256 , _last_number_of_channels(0)
257 , _last_number_of_spectra(0)
260 _profile = profile_manager.profile(config.profile());
264 PulsarInjection<T>::~PulsarInjection()
269 void PulsarInjection<T>::next(DataType& data)
271 float spec_index = _config.spectral_index();
272 const float in_snr = _config.signal_to_noise();
273 float pulse_energy_sigma = 0.2;
275 uint32_t seed = time(NULL);
276 uint_fast32_t nsubpulse = 5;
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();
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);
291 const double noiseAmp = 24;
294 if(nsamp != _last_number_of_spectra || nchan_const != _last_number_of_channels)
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());
306 double scale = noiseAmp * in_snr/sqrt(nchan_const * nsamp)/sum;
307 PANDA_LOG_DEBUG <<
"mean = " << sum <<
", scale = " << scale;
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;
313 const uint_fast32_t nprof = _profile.size();
315 FftwAlloc<float> ism_convolution(nprof);
316 FftwAlloc<float> subpulse_map(nprof);
317 FftwAlloc<float> subpulse_profile(nprof);
318 FftwAlloc<float> unsmeared_profile(nprof);
319 FftwAlloc<float> smeared_profile(nprof);
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);
324 for(std::size_t i = 0; i < nprof; ++i){
325 subpulse_profile[i] = 1.0;
330 for(std::size_t i = 0; i < nprof; ++i) {
331 sum += subpulse_profile[i];
333 sum /= (double) nprof;
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;
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);
351 const double pulse_energy_norm = exp(pow(pulse_energy_sigma, 2)/2.0);
352 double p0 = _phase_model(mjd, data.channel_frequencies()[0]);
354 std::vector<int_fast32_t> prevbin(nchan_const);
358 for (ch = 0; ch < nchan_const; ++ch) {
359 freq[ch] = data.channel_frequencies()[ch];
360 poff[ch] = _phase_model(mjd, freq[ch]) - p0;
362 sidx[ch] = pow(freq[ch].value()/fch1, spec_index);
364 freq[nchan_const] = freq[nchan_const -1] + (freq[nchan_const - 1] - freq[nchan_const -2]);
368 int_fast32_t dm_bin=-1;
369 int_fast32_t scatter_bin=-1;
370 float* dm_conv = ism_convolution;
371 float* scatt_conv = unsmeared_profile;
372 float* out_conv = smeared_profile;
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) {
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;
383 pbin = std::abs(phase) * nprof;
384 if (pbin < 1) pbin = 1;
385 else if (pbin%2 == 0) pbin -= 1;
387 if ((pbin != dm_bin) != scatter_bin) {
390 for(std::size_t i = 0; i < nprof; ++i) {
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);
403 for(std::size_t i = 0; i < nprof; ++i){
408 for(std::size_t i = 1; i < nprof; ++i){
412 convolve(ism_convolution_plan);
414 memcpy(ism_models[nism], out_conv, NCOPY);
419 ism_idx[ch] = nism - 1;
422 PANDA_LOG_DEBUG <<
"Generated " << nism <<
" ISM models";
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));
427 PANDA_LOG_DEBUG <<
"Starting simulation";
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)
434 auto spectra = data.spectrum(sample_index);
435 auto channel_it = spectra.begin();
437 for (ch = 0; ch < nchan_const; ++ch) {
438 rands[ch] = rnd.rand_double();
440 p0 = _phase_model(mjd, freq[0]);
441 ip0 = (int64_t) floor(p0);
444 if (ip0 != prev_ip0) {
445 for(std::size_t i = 0; i < nprof; i++) {
450 std::copy(_profile.begin(), _profile.end(), &unsmeared_profile[0]);
454 for (n = 0; n < nsubpulse; ++n) {
456 std::size_t i = floor(rnd.rand_double() * nprof);
457 subpulse_map[i] += exp(rnd.rand_gauss() * pulse_energy_sigma)/pulse_energy_norm;
461 convolve(subpulse_convolution_plan);
464 std::transform(_profile.begin(), _profile.end(), &unsmeared_profile[0], &unsmeared_profile[0], std::multiplies<float>());
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];
476 grads[i][nprof - 1] = smeared_profile[0] - smeared_profile[nprof - 1];
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;
492 final_signal += output_prof[ism_idx[ch]][prevbin[ch]];
494 dbin = pbin - prevbin[ch];
495 while (dbin < 0) dbin += nprof;
498 if (final_signal > 0) {
499 *channel_it += floor(sidx[ch] * final_signal + rands[ch]);
509 PANDA_LOG_DEBUG <<
"Simulation ended.";
Some limits and constants for FLDO.
Perform a convolution using fft methods.