Fréchet View  1.6.0
A Tool for Exploring Fréchet Distance Algorithms
kalgorithm.cpp
Go to the documentation of this file.
1 
2 #include <k_frechet/kalgorithm.h>
3 
4 using namespace frechet;
5 using namespace data;
6 using namespace k;
7 using namespace app;
8 using namespace fs;
9 
10 kAlgorithm::kAlgorithm(FreeSpace::ptr fs)
11  : components(fs->components()),
12  domain(),
13  working(), greedy(), bf()
14 {
15  domain.H = Interval(0,fs->n-1);
16  domain.V = Interval(0,fs->m-1);
17  //setupMapsAndResults();
18 }
19 
20 
32 {
33  Q_ASSERT(domain.H.lower()==0.0 && domain.V.lower()==0.0);
34 
36  size_t count = *components.last()+1;
37  working.set = BitSet((int)count);
38  greedy.result = BitSet((int)count);
39  bf.result = BitSet((int)count);
40 
41 
42  count = components.count();
43  working.hmap.clear();
44  working.vmap.clear();
45  working.hmap.reserve(count);
46  working.vmap.reserve(count);
47 
48  for(i = components.begin(); i.valid(); ++i)
49  {
50  // i.offset() is the index in the Array = ID of component
51  component_id_t compID = *i;
52  const IntervalPair& ival = components.intervals() [compID];
53 
54  working.hmap.push_back(MappedInterval(ival.H,compID));
55  working.vmap.push_back(MappedInterval(ival.V,compID));
56  }
57 
58  std::sort(working.hmap.begin(),working.hmap.end());
59  std::sort(working.vmap.begin(),working.vmap.end());
60 
63 }
64 
79 {
80  QMutexLocker guard(&working.lock);
81 
83 
84  if (greedy_xy()==NO_SOLUTION) {
85  greedy.k_yx =
87 
88  emit greedyFinished();
89  return NO_SOLUTION;
90  }
91 
93  greedy_yx();
94  Q_ASSERT(greedy.k_xy>0 && greedy.k_yx>0);
95 
96  if (greedy.k_yx < greedy.k_xy)
98 
99  /* The results of the Greedy algorithm
100  * serve as lower and upper bounds for the Brute Force algorithm
101  * */
104 
105  if (bf.k_min == bf.k_max) {
106  // the Greedy result is already optimal. Fine.
107  bf.k_optimal = bf.k_max;
109  }
110  else {
111  bf.k_optimal = UNKNOWN;
112  }
113 
114  emit greedyFinished();
115  return bf.k_max;
116 }
117 
118 
120 {
121  working.set.clear();
122 
125  working.hmap, domain.H,
126  working.vmap, domain.V,
127  greedy.k_x);
128  Q_ASSERT(greedy.k_xy==NO_SOLUTION || greedy.k_xy==working.set.count());
129  return greedy.k_xy;
130 }
131 
133 {
134  working.set.clear();
135 
138  working.vmap, domain.V,
139  working.hmap, domain.H,
140  greedy.k_y);
141  Q_ASSERT(greedy.k_yx==NO_SOLUTION || greedy.k_yx==working.set.count());
142  return greedy.k_yx;
143 }
144 
146  const IntervalMap& map1, const Interval& range1,
147  const IntervalMap& map2, const Interval& range2,
148  int& k_1)
149 {
150  k_1 = greedy_1(map1,range1); // first axis
151  if (k_1 <= 0) return NO_SOLUTION;
152 
153  int k_2 = greedy_2(map2,range2); // second axis
154  // note: k_2 can be 0, if all intervals are already contained in k_1
155  if (k_2 >= 0)
156  return k_1 + k_2;
157  else
158  return NO_SOLUTION;
159 }
160 
161 
162 int kAlgorithm::greedy_1(const IntervalMap& map, const Interval& range)
163 {
164  // scan first axis
165  int k = 0;
166  int i = 0;
167  double x = range.lower();
168 
169  while((x < range.upper()) && (i < map.size()))
170  { // find the largest interval containing x
171  if (map[i].lower() > x)
172  return NO_SOLUTION;
173 
174  double x_next = x;
175  component_id_t comp_next;
176 
177  for( ; (i < map.size()) && (map[i].lower() <= x); i++)
178  {
179  double xupper = map[i].upper();
180  if (xupper > x_next) {
181  x_next = xupper;
182  comp_next = map[i].componentID;
183  }
184  }
185 
186  if (x_next > x) {
187  x = x_next;
188  working_add(comp_next);
189  ++k;
190  }
191  else {
192  return NO_SOLUTION; // we have made no progress
193  }
194  }
195 
196  if (x < range.upper())
197  return NO_SOLUTION;
198  else
199  return k;
200 }
201 
202 int kAlgorithm::greedy_2(const IntervalMap& map, const Interval& range)
203 {
204  // Scan second axis
205  // in addition, we must look for already selected intervals
206  int k = 0;
207  int i = 0;
208  double y = range.lower();
209 
210  while((y < range.upper()) && (i < map.size()))
211  {
212  if (map[i].lower() > y)
213  return NO_SOLUTION; // a gap
214 
215  double y0 = y;
216  double y_next = y;
217  size_t comp_next;
218 
219  for( ; (i < map.size()) && (map[i].lower() <= y); ++i)
220  {
221  double yupper = map[i].upper();
222  if (working.set.contains(map[i].componentID))
223  {
224  if (yupper > y)
225  y = yupper;
226  }
227  else if (yupper > y_next) {
228  y_next = yupper;
229  comp_next = map[i].componentID;
230  }
231  }
232 
233  if (y_next > y) {
234  y = y_next;
235  working_add(comp_next);
236  ++k;
237  }
238  else if (y == y0) {
239  return NO_SOLUTION; // we have made no progress
240  }
241  }
242 
243  if (y < range.upper())
244  return NO_SOLUTION;
245  else
246  return k;
247 }
248 
267 int kAlgorithm::runBruteForce(volatile bool* cancelFlag) throw (InterruptedException)
268 {
269  QMutexLocker guard(&working.lock);
270 
271  if (bf.k_max == NO_SOLUTION)
272  return NO_SOLUTION; // Greedy was unsucessful. No reason to continue.
273 
274  Q_ASSERT(bf.k_min > 0 && bf.k_max > 0);
275  // these must have been calculated by the Greedy algorithm.
276 
277  if (bf.k_min == bf.k_max) {
278  // Greedy is already optimal. That was easy.
279  bf.k_current = bf.k_iteration = bf.k_optimal = bf.k_max;
280  bf.result = greedy.result;
281  return bf.k_optimal;
282  }
283 
284  // TODO delegate to thread and return RUNNING immediately.
285  bf.stack1.reserve(bf.k_max);
286  bf.stack2.reserve(bf.k_max);
287 
288  bf.k_optimal = RUNNING;
289  for(bf.k_iteration = bf.k_min; bf.k_iteration <= bf.k_max; bf.k_iteration++)
290  {
291  // TODO parallel: we could run two searches in parallel (x-y and y-x)
292  // the fastest one wins
293  // for start, search only one way. (or: change direction every iteration?)
294  //working.clear();
295  //bf.k_current=0;
296 
297  emit bruteForceIteration(bf.k_iteration);
298  if (bf_search_1(cancelFlag)) {
299  // found a result
300  Q_ASSERT((bf.stack1.size()+bf.stack2.size()) == bf.k_iteration);
301  // copy stack2 to working set
302  for( ; !bf.stack2.empty(); bf.stack2.pop_back()) {
303  component_id_t comp = working.vmap[bf_top2().first].componentID;
304  Q_ASSERT(!working.set.contains(comp));
305  working.set += comp;
306  }
307  // move working set to bf.result
308  bf.result.swap(working.set);
309  bf.k_optimal = bf.k_iteration;
310  Q_ASSERT(bf.result.count() == bf.k_optimal);
311 
312  emit bruteForceFinished();
313  return bf.k_optimal;
314  }
315  }
316  // We MUST have a result !!
317  Q_ASSERT(false);
318  return NO_SOLUTION;
319 }
320 
322  working_add(ival.componentID);
323 }
324 
327 }
328 
330 {
331  Q_ASSERT(!working.set.contains(comp));
332  working.set += comp;
333 }
334 
336 {
337  Q_ASSERT(working.set.contains(comp));
338  working.set -= comp;
339 }
340 
341 std::pair<int, double> &kAlgorithm::bf_top1() {
342  Q_ASSERT(!bf.stack1.empty());
343  return bf.stack1[bf.stack1.size()-1];
344 }
345 
346 std::pair<int, double> &kAlgorithm::bf_top2() {
347  Q_ASSERT(!bf.stack1.empty());
348  return bf.stack2[bf.stack2.size()-1];
349 }
350 
351 void kAlgorithm::bf_push1(int i, double x) {
352  bf.stack1.push_back(std::make_pair(i,x));
353  Q_ASSERT(working.set.count() == (bf.stack1.size()-1));
354 }
355 
357  Q_ASSERT(bf.stack1.size() > 0);
358  bf.stack1.pop_back();
359  return !bf.stack1.empty();
360 }
361 
362 void kAlgorithm::bf_push2(int i, double x) {
363  bf.stack2.push_back(std::make_pair(i,x));
364  // don't maintain the working set
365 }
366 
368  bf.stack2.clear();
369  return false;
370 }
371 
372 bool kAlgorithm::bf_search_1(volatile bool* cancelFlag) throw (InterruptedException)
373 {
374  // search horizontal range first
375  const IntervalMap& map = working.hmap;
376  Interval range = domain.H;
377 
378  bf.stack1.clear();
379  bf.stack2.clear();
380  working.set.clear();
381  bf.stack1.push_back(std::make_pair(0,domain.H.lower()));
382 
383  for(;;)
384  {
385 outer_loop:
386  if (cancelFlag && *cancelFlag)
387  {
388  bf.k_current = bf.k_iteration = bf.k_optimal = UNKNOWN;
389  throw InterruptedException();
390  }
391 
392  Q_ASSERT(!bf.stack1.empty());
393  int& i = bf_top1().first;
394  double& x = bf_top1().second;
395 
396  for( ; (i < map.size()) && (map[i].lower() <= x); ++i)
397  {
398  if (map[i].upper() <= x) continue; // already covered, not interesting
399  // TODO sort choices?
400  // By interval length -> the first solution is the greedy solution
401  // By partial solution (stored from previous iteration) ?
402 
403  // recurse
404  if (map[i].upper() >= range.upper()) {
405  working_add(map[i]);
406 
407  if (bf_search_2()) return true; // switch to vertical
408 
409  working_remove(map[i]);
410  }
411  else if (bf.stack1.size() < bf.k_iteration) {
412  // push
413  working_add(map[i]);
414  bf_push1(i+1,map[i].upper());
415  goto outer_loop;
416  }
417  }
418 
419  // pop
420  if (!bf_pop1()) return false;
421  working_remove(map[bf_top1().first]);
422  // continue
423  bf_top1().first++;
424  }
425  // never come here
426  Q_ASSERT(false);
427  return false;
428 }
429 
437 bool kAlgorithm::bf_search_2() throw (InterruptedException)
438 {
439 // if (abortRequested)
440 // throw InterruptedException();
441 
442  const IntervalMap& map = working.vmap;
443  Interval range = domain.V;
444 
445  // Scan second axis
446  int i = 0;
447  double y = range.lower();
448 
449  while((y < range.upper()) && (i < map.size()))
450  {
451  if (map[i].lower() > y)
452  return bf_backtrack2(); // a gap
453 
454  double y0 = y;
455  double y_next = y;
456  int i_next = -1;
457 
458  for( ; (i < map.size()) && (map[i].lower() <= y); ++i)
459  {
460  double yupper = map[i].upper();
461  if (working.set.contains(map[i].componentID))
462  {
463  if (yupper > y)
464  y = yupper;
465  }
466  else if (yupper > y_next) {
467  y_next = yupper;
468  i_next = i;
469  }
470  }
471 
472  if (y_next > y) {
473  if ((bf.stack1.size()+bf.stack2.size()) >= bf.k_iteration)
474  return bf_backtrack2(); // reached k-limit
475 
476  Q_ASSERT(i_next >= 0);
477  Q_ASSERT(!working.set.contains(map[i_next].componentID));
478  bf_push2(i_next,y_next);
479  y = y_next;
480  }
481  else if (y == y0)
482  return bf_backtrack2(); // we have made no progress
483  }
484 
485  if (y < range.upper())
486  return bf_backtrack2();
487  else
488  return true;
489 }
490 
491 
492 #define SORT_TIDY
493 
494 bool MappedInterval::operator<(const MappedInterval &that) const {
495  if (lower() != that.lower())
496  return lower() < that.lower();
497 #ifdef SORT_TIDY
498  if (upper() != that.upper())
499  return upper() > that.upper();
500  return componentID < that.componentID;
501 #else
502  // sort arbitrarily (but deterministic, of course)
503  std::hash<const MappedInterval*> hsh;
504  return hsh(this) < hsh(&that);
505 #endif
506 }
507 
509  return (lower()==that.lower())
510  && (upper()==that.upper())
511  && (componentID==that.componentID);
512 }
void bf_push1(int i, double x)
push an value to first stack
Definition: kalgorithm.cpp:351
initial status; the algorithm has not ye been run
Definition: kalgorithm.h:132
Result upperBound() const
get an upper bound for the optimal result. Observe that k_optimal <= k_xy,k_yx
Definition: kalgorithm.h:233
void greedyFinished()
raised when the 'greedy' algorithm finishes Connected to UI components to update their states.
an iterator over a BitSet
Definition: bitset.h:93
int k_x
number of components that is needed to cover the x-axis. 0 if the x-axis can not be covered.
Definition: kalgorithm.h:209
const fs::Components & components
connected components of the free-space diagram
Definition: kalgorithm.h:147
IntervalMap vmap
map to vertical (y-)axis
Definition: kalgorithm.h:162
data::BitSet set
subset of currently selected interval
Definition: kalgorithm.h:164
int k_optimal
optimal value so far
Definition: kalgorithm.h:271
bool bf_search_1(volatile bool *cancelFlag=0)
run the brute-force algorithm on the x-axis. Results are placed in this->bf.
Definition: kalgorithm.cpp:372
void working_remove(const MappedInterval &ival)
remove a mapped interval from the current working-set
Definition: kalgorithm.cpp:325
struct BruteForce bf
result of the brute-force algorithm
Definition: kalgorithm.h:293
Interval H
horizontal interval
Definition: interval.h:459
global definitions for all algorithms.
data::IntervalPair domain
the x-axis and the y-axis of the free-space diagram (from Grid)
Definition: kalgorithm.h:151
std::vector< MappedInterval > IntervalMap
a vector of MappedInterval objects. Usually these are sorted according to the natural ordering impose...
Definition: kalgorithm.h:71
void bf_push2(int i, double y)
push an value to second stack
Definition: kalgorithm.cpp:362
indicates that the algorithm is currently running
Definition: kalgorithm.h:133
Result lowerBound() const
get a lower bound for the optimal result. Observe that k_x,k_y <= k_optimal
Definition: kalgorithm.h:227
int runBruteForce(volatile bool *cancelFlag=0)
run the brute force algorithm. The greedy algorithm must have been run before!
Definition: kalgorithm.cpp:267
a pair of horizonal / vertical intervals.
Definition: interval.h:456
IntervalMap hmap
map to horizontal (x-)axis
Definition: kalgorithm.h:160
int k_iteration
optimal value of last iteration
Definition: kalgorithm.h:269
int k_current
current number of selected components
Definition: kalgorithm.h:267
attaches a component ID to an interval
Definition: kalgorithm.h:32
Interval & clear()
make this an invalid interval
Definition: interval.h:197
void clear()
assign false to all bits
Definition: bitset.cpp:67
int greedy_1(const IntervalMap &map, const data::Interval &range)
run the greedy algorithm on the first axis
Definition: kalgorithm.cpp:162
int runGreedy()
run the greedy algorithm
Definition: kalgorithm.cpp:78
void setupMapsAndResults()
initialize working-set data structures
Definition: kalgorithm.cpp:31
double upper() const
Definition: interval.h:96
component_id_t count() const
Definition: components.h:86
data::BitSet result
subsset of selected components
Definition: kalgorithm.h:274
int greedy_yx()
run greedy algorithm (y-axis first)
Definition: kalgorithm.cpp:132
const IntervalArray & intervals() const
Definition: components.h:99
Interval V
vertical interval
Definition: interval.h:460
int k_yx
number of components, when scanning the y-axis first. 0 if the axes can no be covered.
Definition: kalgorithm.h:215
bool operator<(const MappedInterval &that) const
comparison operator, needed for sorting a list of MappedInterval objects. We first compare the lower ...
Definition: kalgorithm.cpp:494
struct Greedy greedy
results of the greedy algorithm
Definition: kalgorithm.h:246
bool contains(int offset) const
Definition: bitset.cpp:71
bool bf_pop1()
pop from first stack
Definition: kalgorithm.cpp:356
void working_add(const MappedInterval &ival)
add a mapped interval to the current working-set
Definition: kalgorithm.cpp:321
std::pair< int, double > & bf_top2()
Definition: kalgorithm.cpp:346
int k_max
upper bound derived from greedy results. k_max = min{greedy.k_xy,greedy.k_yx}
Definition: kalgorithm.h:265
bool bf_search_2()
run the brute-force algorithm on the y-axis. Results are placed in this->bf.
Definition: kalgorithm.cpp:437
double lower() const
Definition: interval.h:92
data::BitSet::iterator last() const
Definition: components.h:96
data::BitSet result
subset of selected components
Definition: kalgorithm.h:217
data::BitSet::iterator begin() const
Definition: components.h:89
A simple bit vector of fixed size.
Definition: bitset.h:16
bool bf_backtrack2()
clear second stack
Definition: kalgorithm.cpp:367
BitSet & swap(BitSet &that)
swap contents with other object
Definition: bitset.cpp:44
int greedy_2(const IntervalMap &map, const data::Interval &range)
run the greedy algorithm on the second axis
Definition: kalgorithm.cpp:202
bool operator==(const MappedInterval &that) const
equality operator, needed for sorting a list of MappedInterval objects.
Definition: kalgorithm.cpp:508
Stack stack1
stacks used for backtracking
Definition: kalgorithm.h:276
int k_min
lower bound derived from greedy results. k_min = max{greedy.k_x,greedy.k_y}
Definition: kalgorithm.h:263
int count() const
This methods is expensive, as it scans the whole array.
Definition: bitset.cpp:58
std::pair< int, double > & bf_top1()
Definition: kalgorithm.cpp:341
struct frechet::k::kAlgorithm::WorkingSet working
int k_y
number of components that is needed to cover the y-axis. 0 if the y-axis can not be covered.
Definition: kalgorithm.h:211
unsigned int component_id_t
used as identifier for free-space components
Definition: types.h:74
an interval of two double values.
Definition: interval.h:31
int k_xy
number of components, when scanning the x-axis first. 0 if the axes can no be covered.
Definition: kalgorithm.h:213
there is no solution at all
Definition: kalgorithm.h:134
int greedy_12(const IntervalMap &map1, const data::Interval &range1, const IntervalMap &map2, const data::Interval &range2, int &k_1)
run the greedy algorithm on tow axes
Definition: kalgorithm.cpp:145
boost::shared_ptr< FreeSpace > ptr
smart pointer to FreeSpace object
Definition: freespace.h:92
int greedy_xy()
run the greedy algorithm (x-axis first)
Definition: kalgorithm.cpp:119