PROBI  1.0
 All Classes Functions
Weiszfeld.hpp
1 #ifndef WEISZFELD_H
2 #define WEISZFELD_H
3 
4 #include <functional>
5 #include <algorithm>
6 #include <limits>
7 
8 #include "Point.hpp"
9 #include "Metric.hpp"
10 
16 class Weiszfeld
17 {
18 public:
19  class IterationFailed {};
20 
21  Weiszfeld(std::function<Metric<Point>*() > createMetric);
22 
26  template<typename ForwardIterator>
27  Point approximateOneMedian(ForwardIterator begin, ForwardIterator end, int max_iteration = 15);
28 
29 private:
30  Metric<Point>* metric;
31 };
32 
33 template<typename ForwardIterator>
34 Point Weiszfeld::approximateOneMedian(ForwardIterator begin, ForwardIterator end, int max_iteration)
35 {
36  size_t n = 0;
37  double eps = std::numeric_limits<double>::epsilon();
38  int dimension = begin->getDimension();
39 
40  // Compute starting point
41  Point cog(dimension);
42  double weight = 0;
43  for (ForwardIterator it = begin; it != end; ++it)
44  {
45  ++n;
46  weight += it->getWeight();
47  cog += weight * (*it);
48  }
49  cog = (1.0 / weight) * cog;
50 
51  // Check for sets with only one element
52  if(n == 1)
53  return *begin;
54 
55  bool running = true;
56  bool failed = false;
57  Point y(cog);
58  double lastDist = std::numeric_limits<double>::infinity();
59  int iteration = 0;
60  while (running && iteration < max_iteration)
61  {
62  ++iteration;
63  // Compute new iteration point
64  Point numerator(dimension);
65  double denominator = 0;
66  for (ForwardIterator it = begin; it != end; ++it)
67  {
68  double dist = metric->distance(*it, y);
69  if (dist <= eps)
70  {
71  failed = true;
72  break;
73  }
74  double invdist = (1.0 / dist);
75  numerator += it->getWeight() * invdist * (*it);
76  denominator += it->getWeight() * invdist;
77  }
78  Point newPoint((1.0 / denominator) * numerator);
79 
80  // Check for break condition
81  double currentDist = metric->distance(y, newPoint);
82  if (currentDist / lastDist > 0.10 || currentDist <= eps)
83  running = false;
84 
85  // Update iteration data
86  lastDist = currentDist;
87  y = newPoint;
88  }
89 
90  if(failed)
91  throw IterationFailed();
92  else
93  return y;
94 }
95 
96 #endif