1 #include "cheetah/tdao/cuda/Tdao.cuh" 2 #include "cheetah/data/Candidate.h" 3 #include "cheetah/data/Units.h" 4 #include "cheetah/cuda_utils/nvtx.h" 20 inline __host__ __device__
21 bool operator()(thrust::tuple<int,T> t);
29 : _threshold(threshold)
34 inline __host__ __device__
37 return thrust::get<1>(t) > _threshold;
42 using boost::math::complement;
43 using boost::math::quantile;
44 using boost::math::cdf;
47 std::size_t Tdao::_execute(
typename thrust::device_vector<T>::const_iterator in, std::size_t size,
size_t offset,
float threshold)
50 typedef thrust::tuple<unsigned,T> PeakType;
51 typedef thrust::device_vector<unsigned>::iterator IdxIterator;
52 typedef typename thrust::device_vector<T>::iterator PowerIterator;
53 typedef thrust::tuple< thrust::counting_iterator<unsigned>, PowerIterator> PeakTupleInputIterator;
54 typedef typename thrust::tuple< IdxIterator, PowerIterator > PeakTupleOutputIterator;
57 using thrust::counting_iterator;
58 using thrust::zip_iterator;
60 counting_iterator<unsigned> index_counter(offset);
61 auto input_zip = make_zip_iterator(make_tuple(index_counter,in));
62 auto output_zip = make_zip_iterator(make_tuple(_idxs.begin(),_powers.begin()));
64 #ifdef CUDA_BROKEN_THRUST_WITH_ALLOCATOR 65 int num_copied = thrust::copy_if(thrust::cuda::par, input_zip, input_zip+size, output_zip,
68 int num_copied = thrust::copy_if(thrust::cuda::par(_allocator), input_zip, input_zip+size, output_zip,
71 return (std::size_t) num_copied;
74 template <
typename PowerSeriesType>
78 PowerSeriesType
const& input,
79 data::DedispersionMeasureType<float> dm,
81 std::size_t nharmonics)
83 double sigma = input.power_to_equiv_sigma(power);
84 auto period = data::Ccl::CandidateType::MsecTimeType(input.bin_to_period(idx));
85 auto pdot = data::Ccl::CandidateType::SecPerSecType(accel_fact * period.value()/1000.0);
86 auto width = data::Ccl::CandidateType::MsecTimeType(period/(2.0*nharmonics));
87 PANDA_LOG_DEBUG <<
"Candidate -> Period=" <<period<<
", Pdot="<<pdot<<
", Width="<<width<<
", Sigma="<<sigma;
91 template <
typename PowerSeriesType>
92 void Tdao::_filter_unique(PowerSeriesType
const& input,
94 std::size_t num_copied,
95 data::DedispersionMeasureType<float> dm,
96 data::AccelerationType acc,
97 std::size_t nharmonics)
105 thrust::copy(_idxs.begin(),_idxs.begin()+num_copied,_h_idxs.begin());
106 thrust::copy(_powers.begin(),_powers.begin()+num_copied,_h_powers.begin());
110 double accel_fact = (acc / (boost::units::si::constants::codata::c * input.frequency_step())).value();
112 std::size_t distance;
114 float cpeak = _h_powers[idx];
115 int cpeakidx = _h_idxs[idx];
117 std::size_t init_size = output.
size();
118 while(idx < num_copied)
120 distance = _h_idxs[idx]-_h_idxs[idx-1];
123 output.
push_back(_make_candidate(cpeakidx,cpeak,input,dm,accel_fact,nharmonics));
124 cpeak = _h_powers[idx];
125 cpeakidx = _h_idxs[idx];
127 else if (_h_powers[idx] > cpeak)
129 cpeak = _h_powers[idx];
130 cpeakidx = _h_idxs[idx];
134 output.
push_back(_make_candidate(cpeakidx,cpeak,input,dm,accel_fact,nharmonics));
135 std::size_t final_size = output.
size();
136 PANDA_LOG_DEBUG <<
"Found " << final_size-init_size <<
" candidates above threshold of " 137 << _algo_config.significance_threshold() <<
" (DM="<< dm <<
", Acc="<<acc<<
", filtered)";
140 template <
typename T,
typename Alloc>
144 data::DedispersionMeasureType<float> dm,
145 data::AccelerationType acc,
146 std::size_t nharmonics)
148 PUSH_NVTX_RANGE(
"cuda_Tdao_process",0);
149 PANDA_LOG_DEBUG <<
"GPU ID: "<<gpu.device_id();
150 _prepare(input.
size());
152 std::size_t end_bin = (std::size_t) (_algo_config.maximum_frequency()/input.
frequency_step()).value();
153 std::size_t start_bin = (std::size_t) (_algo_config.minimum_frequency()/input.
frequency_step()).value();
154 std::size_t size = std::min(end_bin,input.
size());
155 PANDA_LOG_DEBUG <<
"Bin range to be searched -> ("<<start_bin<<
", "<<size<<
")";
156 PUSH_NVTX_RANGE(
"cuda_Tdao_process_execute_kernels",1);
157 std::size_t num_copied = _execute<T>(input.
begin(),size,start_bin,power_threshold);
159 PANDA_LOG_DEBUG <<
"Found " << num_copied <<
" candidates above power threshold of " 160 << power_threshold <<
" (DM="<< dm <<
", Acc="<<acc<<
", unfiltered)";
161 PUSH_NVTX_RANGE(
"cuda_Tdao_process_filter_unique",2);
162 _filter_unique(input,output,num_copied,dm,acc,nharmonics);
FourierFrequencyType const & frequency_step() const
Retrieve the frequency step of the series.
void process(ResourceType &gpu, data::PowerSeries< Architecture, T, Alloc >const &input, data::Ccl &output, data::DedispersionMeasureType< float > dm, data::AccelerationType acc, std::size_t nharmonics)
Find peaks in a power series above a statistical threshold.
std::size_t size() const
Retrieve the size of the underlying vector.
Some limits and constants for FLDO.
ConstIterator begin() const
Iterators to device memory.
Class for power series (detected FrequencySeries).
void push_back(ValueType const &value)
Appends element to end of vector.
void reserve(std::size_t size)
Reserve space for this many elements.
float equiv_sigma_to_power(float sigma) const
Compute the power that corresponds to a Gaussian equivalent sigma.
std::size_t size() const
the size of the series
A simple record to hold 'candidate' proprerties.