Cheetah - SKA - PSS - Prototype Time Domain Search Pipeline
Fof.cpp
1 /*
2  * The MIT License (MIT)
3  *
4  * Copyright (c) 2016 The SKA organisation
5  *
6  * Permission is hereby granted, free of charge, to any person obtaining a copy
7  * of this software and associated documentation files (the "Software"), to deal
8  * in the Software without restriction, including without limitation the rights
9  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
10  * copies of the Software, and to permit persons to whom the Software is
11  * furnished to do so, subject to the following conditions:
12  *
13  * The above copyright notice and this permission notice shall be included in all
14  * copies or substantial portions of the Software.
15  *
16  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
19  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
21  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
22  * SOFTWARE.
23  */
24 #include <boost/geometry/index/rtree.hpp>
25 #include <boost/mpl/range_c.hpp>
26 #include <boost/mpl/for_each.hpp>
27 #include <chrono>
28 
29 namespace ska {
30 namespace cheetah {
31 namespace sps_clustering {
32 
33 namespace bmpl = boost::mpl;
34 namespace bgi = boost::geometry::index;
35 
36 namespace
37 {
38  // Calculate the square of the euclidian distance between two points
39 
40  template<size_t I, typename PointType>
41  struct DCalc
42  {
43  inline static void calc_distance(size_t dim, const PointType& point1, const PointType& point2, double& d2)
44  {
45  d2 += pow( boost::geometry::get<I>(point1) - boost::geometry::get<I>(point2), 2);
46  DCalc<I-1, PointType>::calc_distance(dim-1, point1, point2, d2);
47  }
48  };
49 
50  template<typename PointType>
51  struct DCalc<0, PointType>
52  {
53  inline static void calc_distance(size_t, const PointType& point1, const PointType& point2, double& d2)
54  {
55  d2 += pow( boost::geometry::get<0>(point1) - boost::geometry::get<0>(point2), 2);
56  }
57  };
58 
59 
60  template<typename PointType>
61  void d2_calc(PointType const& point_1, PointType const& point_2, double& result)
62  {
63  static constexpr std::size_t Rank = boost::geometry::dimension<PointType>::value;
64  DCalc<Rank-1, PointType>::calc_distance(Rank-1, point_1, point_2, result);
65  }
66 
67  // Add a scaler to all the coordinates of a point
68 
69  template<size_t I, typename PointType>
70  struct AddScalar
71  {
72  inline static void add_scalar(size_t dim, PointType& point, double const& c)
73  {
74  double new_coord = boost::geometry::get<I>(point) + c;
75  boost::geometry::set<I>(point, new_coord);
76  AddScalar<I-1, PointType>::add_scalar(dim-1, point,c);
77  }
78  };
79 
80  template<typename PointType>
81  struct AddScalar<0, PointType>
82  {
83  inline static void add_scalar(size_t, PointType& point, double const& c)
84  {
85  double new_coord = boost::geometry::get<0>(point) + c;
86  boost::geometry::set<0>(point, new_coord);
87  }
88  };
89 
90 
91  template<typename PointType>
92  void add_scalar_to_point(PointType& point, double const& value)
93  {
94  static constexpr std::size_t Rank = boost::geometry::dimension<PointType>::value;
95  AddScalar<Rank-1, PointType>::add_scalar(Rank-1, point, value);
96  }
97 } // namespace
98 
99 template<typename NumRepType>
100 std::vector< std::vector<size_t>> Fof::operator()(data::SpCcl<NumRepType> const& cands)
101 {
102  //auto start = std::chrono::high_resolution_clock::now();
103  std::vector<std::vector<size_t>> groups;
104 
105  typedef std::pair<PointType, size_t> ValueType;
106  typedef bgi::rtree<ValueType, bgi::linear<16>> Tree;
107 
108  std::vector< std::pair<PointType, size_t> > points;
109  points.reserve(cands.size());
110 
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 );
116 
117  if(cands.tf_blocks().size() == 0 || cands.size() == 0 ) return groups;
118 
119  auto tsamp = static_cast<boost::units::quantity<data::MilliSeconds,double>>(((*(cands.tf_blocks().begin()))->sample_interval())).value();
120  // set up a point in the clustering search space for each candidate
121  for(size_t i = 0 ; i < cands.size() ; ++i)
122  {
123  PointType point;
124  point.template set<0>(static_cast<double>((cands[i].tstart()/tsamp)/(_config.time_tolerance()/tsamp)));
125 
126  //convert to appropriate scales for proper clustering
127  point.template set<1>(static_cast<double>((cands[i].dm()).value())/((_config.dm_tolerance())).value());
128 
129  point.template set<2>(static_cast<double>(std::log2(cands[i].width().value()/tsamp)/std::log2(_config.pulse_width_tolerance().value()/tsamp)));
130 
131  //set_point(point, point.begin() + i*D);
132  points.emplace_back(point, i);
133  }
134 
135  Tree tree(points.begin(), points.end());
136 
137  while( !tree.empty())
138  {
139  std::vector<ValueType> to_add;
140  // Grab a point from the tree.
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;
144 
145  for( std::size_t to_add_i = 0L; to_add_i < to_add.size(); ++to_add_i )
146  {
147  auto it = to_add.begin() + to_add_i;
148 
149  // Build box to query
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);
154 
155  boost::geometry::model::box<PointType> box( lower, upper );
156 
157  auto within_ball = [this, &it](ValueType const &v)
158  {
159  double d2 = 0.0;
160  d2_calc(it->first, v.first, d2);
161  return d2 < _linking_length_2;
162  };
163 
164  // Find all points within a linking length of the current point.
165  tree.query( bgi::within(box) && bgi::satisfies(within_ball), std::back_inserter(added) );
166 
167  to_add.reserve(to_add.size() + added.size());
168  // Add the found points to the list so we can find their "friends" as well
169  std::copy(added.begin(), added.end(), std::back_inserter(to_add));
170 
171  // Remove any points we find from the tree as they have been assigned.
172  tree.remove( added.begin(), added.end() );
173 
174  added.clear();
175  // Early exit when we have assigned all particles to a group
176  if (tree.empty())
177  {
178  break;
179  }
180  }
181  std::vector< size_t > group;
182  for( auto p : to_add )
183  {
184  group.push_back(p.second);
185  }
186  groups.push_back(group);
187  }
188  //auto stop = std::chrono::high_resolution_clock::now();
189  //auto duration = std::chrono::duration_cast<std::chrono::microseconds>(stop - start);
190  //PANDA_LOG << "Time taken by clustering: " << duration.count() << " us";
191 
192  return groups;
193 }
194 
195 } // namespace sps_clustering
196 } // namespace cheetah
197 } // namespace ska
std::pair< Dm, Dm > const & dm_range() const
the max and minimum Dispersion Measure values of the candidates
Definition: SpCcl.cpp:70
BlocksType const & tf_blocks() const
Access the TimeFrequency blocks associated with the single pulse candidtes.
Definition: SpCcl.cpp:64
std::size_t size() const
Retrieve the size of the underlying vector.
Definition: VectorLike.cpp:61
Some limits and constants for FLDO.
Definition: Brdz.h:35
SpCandidate list.
Definition: SpCcl.h:52
std::vector< std::vector< size_t > > operator()(data::SpCcl< NumRepType > const &cands)
Group the candidates using the fof algorithm.
Definition: Fof.cpp:100