ceph/src/boost/libs/geometry/test/strategies/haversine.cpp

   1 // Boost.Geometry (aka GGL, Generic Geometry Library)
   2 // Unit Test
   3
   4 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
   5 // Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
   6 // Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
   7
   8 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
   9 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
  10
  11 // Use, modification and distribution is subject to the Boost Software License,
  12 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  13 // http://www.boost.org/LICENSE_1_0.txt)
  14
  15
  16 #include <geometry_test_common.hpp>
  17
  18 #include <boost/concept/requires.hpp>
  19 #include <boost/concept_check.hpp>
  20 #include <boost/core/ignore_unused.hpp>
  21
  22 #include <boost/geometry/algorithms/assign.hpp>
  23 #include <boost/geometry/strategies/spherical/distance_haversine.hpp>
  24 #include <boost/geometry/strategies/concepts/distance_concept.hpp>
  25
  26
  27 #include <boost/geometry/geometries/point.hpp>
  28
  29 double const average_earth_radius = 6372795.0;
  30
  31
  32 template <typename Point, typename LatitudePolicy>
  33 struct test_distance
  34 {
  35     typedef bg::strategy::distance::haversine<double> haversine_type;
  36     typedef typename bg::strategy::distance::services::return_type<haversine_type, Point, Point>::type return_type;
  37
  38     BOOST_CONCEPT_ASSERT
  39         (
  40             (bg::concepts::PointDistanceStrategy<haversine_type, Point, Point>)
  41         );
  42
  43
  44     static void test(double lon1, double lat1, double lon2, double lat2,
  45                        double radius, double expected, double tolerance)
  46     {
  47         haversine_type strategy(radius);
  48
  49         Point p1, p2;
  50         bg::assign_values(p1, lon1, LatitudePolicy::apply(lat1));
  51         bg::assign_values(p2, lon2, LatitudePolicy::apply(lat2));
  52         return_type d = strategy.apply(p1, p2);
  53
  54         BOOST_CHECK_CLOSE(d, expected, tolerance);
  55     }
  56 };
  57
  58 template <typename Point, typename LatitudePolicy>
  59 void test_all()
  60 {
  61     // earth to unit-sphere -> divide by earth circumference, then it is from 0-1,
  62     // then multiply with 2 PI, so effectively just divide by earth radius
  63     double e2u = 1.0 / average_earth_radius;
  64
  65     // ~ Amsterdam/Paris, 467 kilometers
  66     double const a_p = 467.2704 * 1000.0;
  67     test_distance<Point, LatitudePolicy>::test(4, 52, 2, 48, average_earth_radius, a_p, 1.0);
  68     test_distance<Point, LatitudePolicy>::test(2, 48, 4, 52, average_earth_radius, a_p, 1.0);
  69     test_distance<Point, LatitudePolicy>::test(4, 52, 2, 48, 1.0, a_p * e2u, 0.001);
  70
  71     // ~ Amsterdam/Barcelona
  72     double const a_b = 1232.9065 * 1000.0;
  73     test_distance<Point, LatitudePolicy>::test(4, 52, 2, 41, average_earth_radius, a_b, 1.0);
  74     test_distance<Point, LatitudePolicy>::test(2, 41, 4, 52, average_earth_radius, a_b, 1.0);
  75     test_distance<Point, LatitudePolicy>::test(4, 52, 2, 41, 1.0, a_b * e2u, 0.001);
  76 }
  77
  78
  79 template <typename P1, typename P2, typename CalculationType, typename LatitudePolicy>
  80 void test_services()
  81 {
  82     namespace bgsd = bg::strategy::distance;
  83     namespace services = bg::strategy::distance::services;
  84
  85     {
  86
  87         // Compile-check if there is a strategy for this type
  88         typedef typename services::default_strategy
  89             <
  90                 bg::point_tag, bg::point_tag, P1, P2
  91             >::type haversine_strategy_type;
  92
  93         boost::ignore_unused<haversine_strategy_type>();
  94     }
  95
  96     P1 p1;
  97     bg::assign_values(p1, 4, 52);
  98
  99     P2 p2;
 100     bg::assign_values(p2, 2, 48);
 101
 102     // ~ Amsterdam/Paris, 467 kilometers
 103     double const km = 1000.0;
 104     double const a_p = 467.2704 * km;
 105     double const expected = a_p;
 106
 107     double const expected_lower = 460.0 * km;
 108     double const expected_higher = 470.0 * km;
 109
 110     // 1: normal, calculate distance:
 111
 112     typedef bgsd::haversine<double, CalculationType> strategy_type;
 113     typedef typename bgsd::services::return_type<strategy_type, P1, P2>::type return_type;
 114
 115     strategy_type strategy(average_earth_radius);
 116     return_type result = strategy.apply(p1, p2);
 117     BOOST_CHECK_CLOSE(result, return_type(expected), 0.001);
 118
 119     // 2: the strategy should return the same result if we reverse parameters
 120     result = strategy.apply(p2, p1);
 121     BOOST_CHECK_CLOSE(result, return_type(expected), 0.001);
 122
 123
 124     // 3: "comparable" to construct a "comparable strategy" for P1/P2
 125     //    a "comparable strategy" is a strategy which does not calculate the exact distance, but
 126     //    which returns results which can be mutually compared (e.g. avoid sqrt)
 127
 128     // 3a: "comparable_type"
 129     typedef typename services::comparable_type<strategy_type>::type comparable_type;
 130
 131     // 3b: "get_comparable"
 132     comparable_type comparable = bgsd::services::get_comparable<strategy_type>::apply(strategy);
 133
 134     // Check vice versa:
 135     // First the result of the comparable strategy
 136     return_type c_result = comparable.apply(p1, p2);
 137     // Second the comparable result of the expected distance
 138     return_type c_expected = services::result_from_distance<comparable_type, P1, P2>::apply(comparable, expected);
 139     // And that one should be equa.
 140     BOOST_CHECK_CLOSE(c_result, return_type(c_expected), 0.001);
 141
 142     // 4: the comparable_type should have a distance_strategy_constructor as well,
 143     //    knowing how to compare something with a fixed distance
 144     return_type c_dist_lower = services::result_from_distance<comparable_type, P1, P2>::apply(comparable, expected_lower);
 145     return_type c_dist_higher = services::result_from_distance<comparable_type, P1, P2>::apply(comparable, expected_higher);
 146
 147     // If this is the case:
 148     BOOST_CHECK(c_dist_lower < c_result && c_result < c_dist_higher);
 149
 150     // Calculate the Haversine by hand here:
 151     return_type c_check = return_type(2.0) * asin(sqrt(c_result)) * average_earth_radius;
 152     BOOST_CHECK_CLOSE(c_check, expected, 0.001);
 153
 154     // This should also be the case
 155     return_type dist_lower = services::result_from_distance<strategy_type, P1, P2>::apply(strategy, expected_lower);
 156     return_type dist_higher = services::result_from_distance<strategy_type, P1, P2>::apply(strategy, expected_higher);
 157     BOOST_CHECK(dist_lower < result && result < dist_higher);
 158 }
 159
 160 /****
 161 template <typename P, typename Strategy>
 162 void time_compare_s(int const n)
 163 {
 164     boost::timer t;
 165     P p1, p2;
 166     bg::assign_values(p1, 1, 1);
 167     bg::assign_values(p2, 2, 2);
 168     Strategy strategy;
 169     typename Strategy::return_type s = 0;
 170     for (int i = 0; i < n; i++)
 171     {
 172         for (int j = 0; j < n; j++)
 173         {
 174             s += strategy.apply(p1, p2);
 175         }
 176     }
 177     std::cout << "s: " << s << " t: " << t.elapsed() << std::endl;
 178 }
 179
 180 template <typename P>
 181 void time_compare(int const n)
 182 {
 183     time_compare_s<P, bg::strategy::distance::haversine<double> >(n);
 184     time_compare_s<P, bg::strategy::distance::comparable::haversine<double> >(n);
 185 }
 186
 187 #include <time.h>
 188 double time_sqrt(int n)
 189 {
 190     clock_t start = clock();
 191
 192     double v = 2.0;
 193     double s = 0.0;
 194     for (int i = 0; i < n; i++)
 195     {
 196         for (int j = 0; j < n; j++)
 197         {
 198             s += sqrt(v);
 199             v += 1.0e-10;
 200         }
 201     }
 202     clock_t end = clock();
 203     double diff = double(end - start) / CLOCKS_PER_SEC;
 204
 205     std::cout << "Check: " << s << " Time: " << diff << std::endl;
 206     return diff;
 207 }
 208
 209 double time_normal(int n)
 210 {
 211     clock_t start = clock();
 212
 213     double v = 2.0;
 214     double s = 0.0;
 215     for (int i = 0; i < n; i++)
 216     {
 217         for (int j = 0; j < n; j++)
 218         {
 219             s += v;
 220             v += 1.0e-10;
 221         }
 222     }
 223     clock_t end = clock();
 224     double diff = double(end - start) / CLOCKS_PER_SEC;
 225
 226     std::cout << "Check: " << s << " Time: " << diff << std::endl;
 227     return diff;
 228 }
 229 ***/
 230
 231 int test_main(int, char* [])
 232 {
 233     test_all<bg::model::point<int, 2, bg::cs::spherical_equatorial<bg::degree> >, geographic_policy>();
 234     test_all<bg::model::point<float, 2, bg::cs::spherical_equatorial<bg::degree> >, geographic_policy>();
 235     test_all<bg::model::point<double, 2, bg::cs::spherical_equatorial<bg::degree> >, geographic_policy>();
 236
 237     // NYI: haversine for mathematical spherical coordinate systems
 238     // test_all<bg::model::point<double, 2, bg::cs::spherical<bg::degree> >, mathematical_policy>();
 239
 240     //double t1 = time_sqrt(20000);
 241     //double t2 = time_normal(20000);
 242     //std::cout << "Factor: " << (t1 / t2) << std::endl;
 243     //time_compare<bg::model::point<double, 2, bg::cs::spherical<bg::radian> > >(10000);
 244
 245     test_services
 246         <
 247             bg::model::point<double, 2, bg::cs::spherical_equatorial<bg::degree> >,
 248             bg::model::point<double, 2, bg::cs::spherical_equatorial<bg::degree> >,
 249             double,
 250             geographic_policy
 251         >();
 252
 253     return 0;
 254 }