Group the candidates using the fof algorithm.
103 std::vector<std::vector<size_t>> groups;
105 typedef std::pair<PointType, size_t> ValueType;
106 typedef bgi::rtree<ValueType, bgi::linear<16>> Tree;
108 std::vector< std::pair<PointType, size_t> > points;
109 points.reserve(cands.size());
111 typedef typename data::SpCcl<NumRepType>::Dm Dm;
112 std::pair<Dm, Dm> dm_range = cands.dm_range();
113 Dm dm_step=(dm_range.second - dm_range.first)/static_cast<typename Dm::value_type>(std::numeric_limits<std::size_t>::max());
114 if (dm_step == 0 * pss::astrotypes::units::parsecs_per_cube_cm)
115 dm_step = (1.0 * pss::astrotypes::units::parsecs_per_cube_cm );
117 if(cands.tf_blocks().size() == 0 || cands.size() == 0 )
return groups;
119 auto tsamp =
static_cast<boost::units::quantity<data::MilliSeconds,double>
>(((*(cands.tf_blocks().begin()))->sample_interval())).value();
121 for(
size_t i = 0 ; i < cands.size() ; ++i)
124 point.template set<0>(
static_cast<double>((cands[i].tstart()/tsamp)/(_config.
time_tolerance()/tsamp)));
127 point.template set<1>(
static_cast<double>((cands[i].dm()).value())/((_config.
dm_tolerance())).value());
129 point.template set<2>(
static_cast<double>(std::log2(cands[i].width().value()/tsamp)/std::log2(_config.
pulse_width_tolerance().value()/tsamp)));
132 points.emplace_back(point, i);
135 Tree tree(points.begin(), points.end());
137 while( !tree.empty())
139 std::vector<ValueType> to_add;
141 to_add.push_back( *tree.qbegin( bgi::satisfies([](ValueType
const &){ return true; })) );
142 tree.remove( to_add.begin(), to_add.end() );
143 std::vector<ValueType> added;
145 for( std::size_t to_add_i = 0L; to_add_i < to_add.size(); ++to_add_i )
147 auto it = to_add.begin() + to_add_i;
150 PointType lower = it->first;
151 add_scalar_to_point(lower, -_linking_length);
152 PointType upper = it->first;
153 add_scalar_to_point(upper, _linking_length);
155 boost::geometry::model::box<PointType> box( lower, upper );
157 auto within_ball = [
this, &it](ValueType
const &v)
160 d2_calc(it->first, v.first, d2);
161 return d2 < _linking_length_2;
165 tree.query( bgi::within(box) && bgi::satisfies(within_ball), std::back_inserter(added) );
167 to_add.reserve(to_add.size() + added.size());
169 std::copy(added.begin(), added.end(), std::back_inserter(to_add));
172 tree.remove( added.begin(), added.end() );
181 std::vector< size_t > group;
182 for(
auto p : to_add )
184 group.push_back(p.second);
186 groups.push_back(group);
MsecTimeType pulse_width_tolerance() const
: Threshold for pulse widths between two candidates to be clustered in a single group ...
Dm dm_tolerance() const
: Threshold for DM of two candidates to be clustered in a single group
MsecTimeType time_tolerance() const
: Threshold for time between two candidates to be clustered in a single group