PROBI  1.0
 All Classes Functions
LloydProbMedian.hpp
1 #ifndef LLOYDPROBMEDIAN_HPP
2 #define LLOYDPROBMEDIAN_HPP
3 
4 #include <algorithm>
5 #include <utility>
6 #include <memory>
7 
8 #include "Metric.hpp"
9 #include "WeightedPoint.hpp"
10 #include "AdaptiveSampling.hpp"
11 #include "Weiszfeld.hpp"
12 #include "CenterOfGravity.hpp"
13 #include "KumarMedian.hpp"
14 #include "ProbabilisticPoint.hpp"
15 
22 {
23 private:
24  AdaptiveSampling sampling;
25  Weiszfeld weiszfeld;
26  CenterOfGravity centerOfGravity;
27  KumarMedian kumar;
28  Metric<Point>* metric;
29 
30  int weiszfeldMaxIt;
31  int kumarMedianIterations;
32 public:
33  LloydProbMedian(std::function<Metric<Point>*() > createMetric, std::function<Norm<Point>*() > createNorm);
34 
44  template<typename ForwardIterator, typename OutputIterator>
45  void computeCenterSet(ForwardIterator begin, ForwardIterator end, OutputIterator output, size_t k, size_t maxIterations, size_t n = 0);
46 
47  virtual ~LloydProbMedian();
48 };
49 
50 template<typename ForwardIterator, typename OutputIterator>
51 void LloydProbMedian::computeCenterSet(ForwardIterator begin, ForwardIterator end, OutputIterator output, size_t k, size_t maxIterations, size_t n)
52 {
53  //std::cout << "LloydProbMedian" << std::endl;
54  if (n == 0)
55  for (ForwardIterator it = begin; it != end; ++it)
56  ++n;
57  int dimension = (*begin)[0].getDimension();
58 
59  //std::cout << "Bastele Realisationsvektor (" << n << ")" << std::endl;
60  std::vector<WeightedPoint> allRealizations;
61  int z = 0;
62  for (auto it = begin; it != end; ++it)
63  {
64  ++z;
65  //std::cout << z << ", ";
66  ProbabilisticPoint const & pp = *it;
67  for (auto wpit = pp.cbegin(); wpit != pp.cend(); ++wpit)
68  {
69  WeightedPoint wp = *wpit;
70  wp.setWeight(wp.getWeight() * pp.getWeight());
71  allRealizations.push_back(wp);
72  }
73  }
74  int m = allRealizations.size();
75  std::unique_ptr < std::vector < Point >> initialCenters = sampling.computeCenterSet(allRealizations.begin(), allRealizations.end(), k, m);
76 
77  std::vector<Point> centers(*initialCenters);
78  std::vector<size_t> centerAssignmentIndices(n, 0);
79  bool assignmentChanged = true;
80  for (size_t i = 0; i < maxIterations && assignmentChanged; ++i)
81  {
82  //std::cout << "Iteration " << i << std::endl;
83  assignmentChanged = false;
84  std::vector < std::vector < WeightedPoint >> centerAssignments(k, std::vector<WeightedPoint>());
85  size_t p = 0;
86  for (ForwardIterator it = begin; it != end; ++it)
87  {
88  ProbabilisticPoint const & pp = *it;
89  double minDist = std::numeric_limits<double>::infinity();
90  size_t minCenter = 0;
91  for (size_t c = 0; c < centers.size(); ++c)
92  {
93  double tmpDist = 0;
94  for (auto wpit = pp.cbegin(); wpit != pp.cend(); ++wpit)
95  {
96  WeightedPoint const & wp = *wpit;
97  tmpDist += wp.getWeight() * metric->distance(centers[c], *wpit);
98  }
99  if (tmpDist < minDist)
100  {
101  minDist = tmpDist;
102  minCenter = c;
103  }
104  }
105  if (centerAssignmentIndices[p] != minCenter)
106  {
107  centerAssignmentIndices[p] = minCenter;
108  assignmentChanged = true;
109  }
110  for (auto wpit = pp.cbegin(); wpit != pp.cend(); ++wpit)
111  centerAssignments[minCenter].push_back(*wpit);
112  ++p;
113  }
114 
115  for (size_t c = 0; c < centers.size(); ++c)
116  {
117  if (centerAssignments[c].size() > 0)
118 #ifdef KMEANS
119  centers[c] = centerOfGravity.cog(centerAssignments[c].begin(), centerAssignments[c].end());
120 #else
121  try
122  {
123  centers[c] = weiszfeld.approximateOneMedian(centerAssignments[c].begin(), centerAssignments[c].end(), weiszfeldMaxIt);
124  }
125  catch (Weiszfeld::IterationFailed err)
126  {
127  centers[c] = kumar.approximateOneMedianRounds(centerAssignments[c].begin(), centerAssignments[c].end(), 0.9999999999, weiszfeldMaxIt);
128  }
129 #endif
130 else
131  centers[c] = Point(dimension);
132  }
133  }
134 
135  for (size_t i = 0; i < centers.size(); ++i)
136  {
137  *output = centers[i];
138  ++output;
139  }
140 }
141 
142 #endif /* LLOYDMEDIAN_HPP */
143