source: orange/source/orange/hclust.cpp @ 11751:57be32d2ac3d

Revision 11751:57be32d2ac3d, 29.5 KB checked in by Ales Erjavec <ales.erjavec@…>, 5 months ago (diff)

Fix unordered_map include when using libc++ (clang).

Also omit the tr1 namespace if the compiler supports the c++11 standard.

(fixes #1341)

Line 
1
2#include <ciso646>
3#if defined(_LIBCPP_VERSION) || __cplusplus >= 201103L
4    // C++11 support or using libc++
5    #define HAVE_STD_UNORDERED_MAP
6#elif _MSC_VER
7    // At least on supported versions of MSVC (>= 2008)  we have
8    // std::unordered_map
9    #define HAVE_STD_UNORDERED_MAP
10#endif
11
12#ifdef HAVE_STD_UNORDERED_MAP
13    #include <unordered_map>
14    #define _HCLUST_UNORDERED_MAP std::unordered_map
15#else
16    // include and use from tr1 namespace
17    #include <tr1/unordered_map>
18    #define _HCLUST_UNORDERED_MAP std::tr1::unordered_map
19#endif
20
21
22#include "progress.hpp"
23
24#include "hclust.hpp"
25
26DEFINE_TOrangeVector_classDescription(PHierarchicalCluster, "THierarchicalClusterList", true, ORANGE_API)
27
28#include "hclust.ppp"
29
30
31class TClusterW {
32public:
33    TClusterW *next; // for cluster, this is next cluster, for elements next element
34    TClusterW *left, *right; // subclusters, if left==NULL, this is an element
35    int size;
36    int elementIndex;
37    float height;
38
39    float *distances; // distances to clusters before this one (lower left matrix)
40    float minDistance; // minimal distance
41    int rawIndexMinDistance; // index of minimal distance
42    int nDistances;
43
44    TClusterW(const int &elIndex, float *adistances, const int &anDistances)
45    : next(NULL),
46      left(NULL),
47      right(NULL),
48      size(1),
49      elementIndex(elIndex),
50      height(0.0),
51      distances(adistances),
52      minDistance(numeric_limits<float>::max()),
53      rawIndexMinDistance(-1),
54      nDistances(anDistances)
55    {
56      if (distances)
57        computeMinimalDistance();
58    }
59
60    void elevate(TClusterW **aright, const float &aheight)
61    {
62      left = mlnew TClusterW(*this);
63      right = *aright;
64
65      left->distances = NULL;
66      size = left->size + right->size;
67      elementIndex = -1;
68      height = aheight;
69
70      if (next == right)
71        next = right->next;
72      else
73        *aright = right->next;
74    }
75
76
77    void computeMinimalDistance()
78    {
79      float *dp = distances, *minp = dp++;
80      for(int i = nDistances; --i; dp++)
81        if ((*dp >= 0) && (*dp < *minp))
82          minp = dp;
83      minDistance = *minp;
84      rawIndexMinDistance = minp - distances;
85    }
86
87    TClusterW **clusterAt(int rawIndex, TClusterW **cluster)
88    {
89      for(float *dp = distances; rawIndex; dp++)
90        if (*dp>=0) {
91          cluster = &(*cluster)->next;
92          rawIndex--;
93        }
94      return cluster;
95    }
96
97    void clearDistances()
98    {
99      raiseErrorWho("TClusterW", "rewrite clean to adjust indexMinDistance");
100
101      float *dst = distances;
102      int i = nDistances;
103      for(; nDistances && (*dst >= 0); nDistances--, dst++);
104      if (!nDistances)
105        return;
106     
107      float *src = dst;
108      for(src++, nDistances--; nDistances && (*src >= 0); src++)
109        if (*src >= 0)
110          *dst++ = *src;
111
112      nDistances = dst - distances;
113    }
114};
115
116
117THierarchicalCluster::THierarchicalCluster()
118: branches(),
119  height(0.0),
120  mapping(),
121  first(0),
122  last(0)
123{}
124
125
126THierarchicalCluster::THierarchicalCluster(PIntList els, const int &elementIndex)
127: branches(),
128  height(0.0),
129  mapping(els),
130  first(elementIndex),
131  last(elementIndex+1)
132{}
133
134
135THierarchicalCluster::THierarchicalCluster(PIntList els, PHierarchicalCluster left, PHierarchicalCluster right, const float &h, const int &f, const int &l)
136: branches(new THierarchicalClusterList(2)),
137  height(h),
138  mapping(els),
139  first(f),
140  last(l)
141{ 
142  branches->at(0) = left;
143  branches->at(1) = right;
144}
145
146
147void THierarchicalCluster::swap()
148{
149  if (!branches || (branches->size()<2))
150    return;
151  if (branches->size() > 2)
152    raiseError("cannot swap multiple branches (use method 'permutation' instead)");
153
154  const TIntList::iterator beg0 = mapping->begin() + branches->at(0)->first;
155  const TIntList::iterator beg1 = mapping->begin() + branches->at(1)->first;
156  const TIntList::iterator end1 = mapping->begin() + branches->at(1)->last;
157
158  if ((branches->at(0)->first > branches->at(1)->first) || (branches->at(1)->first > branches->at(1)->last))
159    raiseError("internal inconsistency in clustering structure: invalid ordering of left's and right's 'first' and 'last'");
160
161  TIntList::iterator li0, li1;
162
163  int *temp = new int [beg1 - beg0], *t;
164  for(li0 = beg0, t = temp; li0 != beg1; *t++ = *li0++);
165  for(li0 = beg0, li1 = beg1; li1 != end1; *li0++ = *li1++);
166  for(t = temp; li0 != end1; *li0++ = *t++);
167  delete temp;
168
169  branches->at(0)->recursiveMove(end1 - beg1);
170  branches->at(1)->recursiveMove(beg0 - beg1);
171
172  PHierarchicalCluster tbr = branches->at(0);
173  branches->at(0) = branches->at(1);
174  branches->at(1) = tbr;
175}
176
177
178void THierarchicalCluster::permute(const TIntList &neworder)
179{
180  if ((!branches && neworder.size()) || (branches->size() != neworder.size()))
181    raiseError("the number of clusters does not match the lenght of the permutation vector");
182
183  int *temp = new int [last - first], *t = temp;
184  TIntList::const_iterator pi = neworder.begin();
185  THierarchicalClusterList::iterator bi(branches->begin()), be(branches->end());
186  THierarchicalClusterList newBranches;
187
188  for(; bi != be; bi++, pi++) {
189    PHierarchicalCluster branch = branches->at(*pi);
190    newBranches.push_back(branch);
191    TIntList::const_iterator bei(mapping->begin() + branch->first), bee(mapping->begin() + branch->last);
192    const int offset = (t - temp) - (branch->first - first);
193    for(; bei != bee; *t++ = *bei++);
194    if (offset)
195      branch->recursiveMove(offset);
196  }
197
198  TIntList::iterator bei(mapping->begin() + first), bee(mapping->begin() + last);
199  for(t = temp; bei!=bee; *bei++ = *t++);
200
201  bi = branches->begin();
202  THierarchicalClusterList::const_iterator nbi(newBranches.begin());
203  for(; bi != be; *bi++ = *nbi++);
204}
205
206
207void THierarchicalCluster::recursiveMove(const int &offset)
208{
209  first += offset;
210  last += offset;
211  if (branches)
212    PITERATE(THierarchicalClusterList, bi, branches)
213      (*bi)->recursiveMove(offset);
214}
215
216
217THierarchicalClustering::THierarchicalClustering()
218: linkage(Single),
219  overwriteMatrix(false)
220{}
221
222
223TClusterW **THierarchicalClustering::init(const int &dim, float *distanceMatrix)
224{
225  for(float *ddi = distanceMatrix, *dde = ddi + ((dim+1)*(dim+2))/2; ddi!=dde; ddi++)
226    if (*ddi < 0) {
227      int x, y;
228      TSymMatrix::index2coordinates(ddi-distanceMatrix, x, y);
229      raiseError("distance matrix contains negative element at (%i, %i)", x, y);
230    }
231
232  TClusterW **clusters = mlnew TClusterW *[dim];
233  TClusterW **clusteri = clusters;
234
235  *clusters = mlnew TClusterW(0, NULL, 0);
236  distanceMatrix++;
237 
238  for(int elementIndex = 1, e = dim; elementIndex < e; distanceMatrix += ++elementIndex) {
239    TClusterW *newcluster = mlnew TClusterW(elementIndex, distanceMatrix, elementIndex);
240    (*clusteri++)->next = newcluster;
241    *clusteri = newcluster; 
242  }
243
244  return clusters;
245}
246
247
248
249TClusterW *THierarchicalClustering::merge_SingleLinkage(TClusterW **clusters, float *milestones)
250{
251  float *milestone = milestones;
252
253  int step = 0;
254  while((*clusters)->next) {
255    if (milestone && (step++ ==*milestone))
256      progressCallback->call(*((++milestone)++));
257
258    TClusterW *cluster;
259    TClusterW **pcluster2;
260
261    float minDistance = numeric_limits<float>::max();
262    for(TClusterW **tcluster = &((*clusters)->next); *tcluster; tcluster = &(*tcluster)->next)
263      if ((*tcluster)->minDistance < minDistance) {
264        minDistance = (*tcluster)->minDistance;
265        pcluster2 = tcluster;
266      }
267
268    TClusterW *const cluster2 = *pcluster2;
269
270    const int rawIndex1 = cluster2->rawIndexMinDistance;
271    const int rawIndex2 = cluster2->nDistances;
272
273    TClusterW *const cluster1 = clusters[rawIndex1];
274
275    float *disti1 = cluster1->distances;
276    float *disti2 = cluster2->distances;
277
278    if (rawIndex1) { // not root - has no distances...
279      const float *minIndex1 = cluster1->distances + cluster1->rawIndexMinDistance;
280      for(int ndi = cluster1->nDistances; ndi--; disti1++, disti2++)
281        if (*disti2 < *disti1) // if one is -1, they both are
282          if ((*disti1 = *disti2) < *minIndex1)
283            minIndex1 = disti1;
284      cluster1->minDistance = *minIndex1;
285      cluster1->rawIndexMinDistance = minIndex1 - cluster1->distances;
286    }
287
288    while(*disti2 < 0)
289      disti2++;        // should have at least one more >=0  - the one corresponding to distance to cluster1
290
291    for(cluster = cluster1->next; cluster != cluster2; cluster = cluster->next) {
292      while(*++disti2 < 0); // should have more - the one corresponding to cluster
293      if (*disti2 < cluster->distances[rawIndex1])
294        if ((cluster->distances[rawIndex1] = *disti2) < cluster->minDistance) {
295          cluster->minDistance = *disti2;
296          cluster->rawIndexMinDistance = rawIndex1;
297        }
298    }
299
300    for(cluster = cluster->next; cluster; cluster = cluster->next) {
301      if (cluster->distances[rawIndex2] < cluster->distances[rawIndex1])
302        cluster->distances[rawIndex1] = cluster->distances[rawIndex2];
303
304      // don't nest this in the above if -- both distances can be equal, yet index HAS TO move
305      if (rawIndex2 == cluster->rawIndexMinDistance)
306        cluster->rawIndexMinDistance = rawIndex1; // the smallest element got moved
307
308      cluster->distances[rawIndex2] = -1;
309    }
310
311    cluster1->elevate(pcluster2, minDistance);
312  }
313
314  return *clusters;
315}
316
317// Also computes Ward's linkage
318TClusterW *THierarchicalClustering::merge_AverageLinkage(TClusterW **clusters, float *milestones)
319{
320  float *milestone = milestones;
321  bool ward = linkage == Ward;
322
323  int step = 0;
324  while((*clusters)->next) {
325    if (milestone && (step++ ==*milestone))
326      progressCallback->call(*((++milestone)++));
327
328    TClusterW *cluster;
329    TClusterW **pcluster2;
330
331    float minDistance = numeric_limits<float>::max();
332    for(TClusterW **tcluster = &((*clusters)->next); *tcluster; tcluster = &(*tcluster)->next)
333      if ((*tcluster)->minDistance < minDistance) {
334        minDistance = (*tcluster)->minDistance;
335        pcluster2 = tcluster;
336      }
337
338    TClusterW *const cluster2 = *pcluster2;
339
340    const int rawIndex1 = cluster2->rawIndexMinDistance;
341    const int rawIndex2 = cluster2->nDistances;
342
343    TClusterW *const cluster1 = clusters[rawIndex1];
344
345    float *disti1 = cluster1->distances;
346    float *disti2 = cluster2->distances;
347
348    const int size1 = cluster1->size;
349    const int size2 = cluster2->size;
350    const int sumsize = size1 + size2;
351    cluster = (*clusters)->next;
352
353    if (rawIndex1) { // not root - has no distances...
354      const float sizeK = cluster->size;
355      *disti1 = ward ? (*disti1 * (size1+sizeK) + *disti2 * (size2+sizeK) - minDistance * sizeK) / (sumsize+sizeK)
356                     : (*disti1 * size1 + *disti2 * size2) / sumsize;
357      const float *minIndex1 = disti1;
358      int ndi = cluster1->nDistances-1;
359      for(disti1++, disti2++, cluster = cluster->next; ndi--; disti1++, disti2++)
360        if (*disti1 >= 0) {
361          const float sizeK = cluster->size;
362          cluster = cluster->next;
363          *disti1 = ward ? (*disti1 * (size1+sizeK) + *disti2 * (size2+sizeK) - minDistance * sizeK) / (sumsize+sizeK)
364                         : (*disti1 * size1         + *disti2 * size2) / sumsize;
365          if (*disti1 < *minIndex1)
366            minIndex1 = disti1;
367        }
368      cluster1->minDistance = *minIndex1;
369      cluster1->rawIndexMinDistance = minIndex1 - cluster1->distances;
370    }
371
372    while(*disti2 < 0)
373      disti2++;        // should have at least one more >=0  - the one corresponding to distance to cluster1
374
375    for(cluster = cluster1->next; cluster != cluster2; cluster = cluster->next) {
376      while(*++disti2 < 0); // should have more - the one corresponding to cluster
377      float &distc = cluster->distances[rawIndex1];
378      const float sizeK = cluster->size;
379      distc = ward ? (distc * (size1+sizeK) + *disti2 * (size2+sizeK) - minDistance * sizeK) / (sumsize+sizeK)
380                   : (distc * size1         + *disti2 * size2) / sumsize;
381      if (distc < cluster->minDistance) {
382        cluster->minDistance = distc;
383        cluster->rawIndexMinDistance = rawIndex1;
384      }
385      else if ((distc > cluster->minDistance) && (cluster->rawIndexMinDistance == rawIndex1))
386        cluster->computeMinimalDistance();
387    }
388
389    for(cluster = cluster->next; cluster; cluster = cluster->next) {
390      float &distc = cluster->distances[rawIndex1];
391      const float sizeK = cluster->size;
392      distc = ward ? (distc * (size1+sizeK) + cluster->distances[rawIndex2] * (size2+sizeK) - minDistance * sizeK) / (sumsize+sizeK)
393                   : (distc * size1         + cluster->distances[rawIndex2] * size2) / sumsize;
394      cluster->distances[rawIndex2] = -1;
395      if (distc < cluster->minDistance) {
396        cluster->minDistance = distc;
397        cluster->rawIndexMinDistance = rawIndex1;
398      }
399      else if (   (distc > cluster->minDistance) && (cluster->rawIndexMinDistance == rawIndex1)
400               || (cluster->rawIndexMinDistance == rawIndex2))
401        cluster->computeMinimalDistance();
402    }
403
404    cluster1->elevate(pcluster2, minDistance);
405  }
406
407  return *clusters;
408}
409
410
411
412TClusterW *THierarchicalClustering::merge_CompleteLinkage(TClusterW **clusters, float *milestones)
413{
414  float *milestone = milestones;
415
416  int step = 0;
417  while((*clusters)->next) {
418    if (milestone && (step++ ==*milestone))
419      progressCallback->call(*((++milestone)++));
420
421    TClusterW *cluster;
422    TClusterW **pcluster2;
423
424    float minDistance = numeric_limits<float>::max();
425    for(TClusterW **tcluster = &((*clusters)->next); *tcluster; tcluster = &(*tcluster)->next)
426      if ((*tcluster)->minDistance < minDistance) {
427        minDistance = (*tcluster)->minDistance;
428        pcluster2 = tcluster;
429      }
430
431    TClusterW *const cluster2 = *pcluster2;
432
433    const int rawIndex1 = cluster2->rawIndexMinDistance;
434    const int rawIndex2 = cluster2->nDistances;
435
436    TClusterW *const cluster1 = clusters[rawIndex1];
437
438    float *disti1 = cluster1->distances;
439    float *disti2 = cluster2->distances;
440
441    if (rawIndex1) { // not root - has no distances...
442      if (*disti2 > *disti1)
443        *disti1 = *disti2;
444      const float *minIndex1 = disti1;
445      int ndi = cluster1->nDistances-1;
446      for(disti1++, disti2++; ndi--; disti1++, disti2++) {
447        if (*disti1 >= 0) {
448          if (*disti2 > *disti1) // if one is -1, they both are
449            *disti1 = *disti2;
450          if (*disti1 < *minIndex1)
451            minIndex1 = disti1;
452        }
453      }
454      cluster1->minDistance = *minIndex1;
455      cluster1->rawIndexMinDistance = minIndex1 - cluster1->distances;
456    }
457
458    while(*disti2 < 0)
459      disti2++;        // should have at least one more >=0  - the one corresponding to distance to cluster1
460
461    for(cluster = cluster1->next; cluster != cluster2; cluster = cluster->next) {
462      while(*++disti2 < 0); // should have more - the one corresponding to cluster
463      float &distc = cluster->distances[rawIndex1];
464      if (*disti2 > distc) {
465        distc = *disti2;
466        if (cluster->rawIndexMinDistance == rawIndex1)
467          if (distc <= cluster->minDistance)
468            cluster->minDistance = distc;
469          else
470            cluster->computeMinimalDistance();
471      }
472    }
473
474    for(cluster = cluster->next; cluster; cluster = cluster->next) {
475      float &distc = cluster->distances[rawIndex1];
476      if (cluster->distances[rawIndex2] > distc)
477        distc = cluster->distances[rawIndex2];
478
479      cluster->distances[rawIndex2] = -1;
480      if ((cluster->rawIndexMinDistance == rawIndex1) || (cluster->rawIndexMinDistance == rawIndex2))
481        cluster->computeMinimalDistance();
482    }
483
484    cluster1->elevate(pcluster2, minDistance);
485  }
486
487  return *clusters;
488}
489
490
491
492TClusterW *THierarchicalClustering::merge(TClusterW **clusters, float *milestones)
493{
494  switch(linkage) {
495    case Complete: return merge_CompleteLinkage(clusters, milestones);
496    case Single: return merge_SingleLinkage(clusters, milestones);
497    case Average: 
498    case Ward:
499    default: return merge_AverageLinkage(clusters, milestones);
500  }
501}
502
503PHierarchicalCluster THierarchicalClustering::restructure(TClusterW *root)
504{
505  PIntList elementIndices = new TIntList(root->size);
506  TIntList::iterator currentElement(elementIndices->begin());
507  int currentIndex = 0;
508
509  return restructure(root, elementIndices, currentElement, currentIndex);
510}
511
512
513PHierarchicalCluster THierarchicalClustering::restructure(TClusterW *node, PIntList elementIndices, TIntList::iterator &currentElement, int &currentIndex)
514{
515  PHierarchicalCluster cluster;
516
517  if (!node->left) {
518    *currentElement++ = node->elementIndex;
519    cluster = mlnew THierarchicalCluster(elementIndices, currentIndex++);
520  }
521  else {
522    PHierarchicalCluster left = restructure(node->left, elementIndices, currentElement, currentIndex);
523    PHierarchicalCluster right = restructure(node->right, elementIndices, currentElement, currentIndex);
524    cluster = mlnew THierarchicalCluster(elementIndices, left, right, node->height, left->first, right->last);
525  }
526
527  // no need to care about 'distances' - they've been removed during clustering (in 'elevate') :)
528  mldelete node;
529  return cluster;
530}
531
532
533PHierarchicalCluster THierarchicalClustering::operator()(PSymMatrix distanceMatrix)
534{
535  float *distanceMatrixElements = NULL;
536  TClusterW **clusters, *root;
537  float *callbackMilestones = NULL;
538 
539  try {
540    const int dim = distanceMatrix->dim;
541    const int size = ((dim+1)*(dim+2))/2;
542    float *distanceMatrixElements = overwriteMatrix ? distanceMatrix->elements : (float *)memcpy(new float[size], distanceMatrix->elements, size*sizeof(float));
543
544    clusters = init(dim, distanceMatrixElements);
545    callbackMilestones = (progressCallback && (distanceMatrix->dim >= 1000)) ? progressCallback->milestones(distanceMatrix->dim) : NULL;
546    root = merge(clusters, callbackMilestones);
547  }
548  catch (...) {
549    mldelete clusters;
550    mldelete callbackMilestones;
551    mldelete distanceMatrixElements;
552    throw;
553  }
554
555  mldelete clusters;
556  mldelete callbackMilestones;
557  mldelete distanceMatrixElements;
558
559  return restructure(root);
560}
561
562
563/*
564 *  Optimal leaf ordering.
565 */
566
567struct m_element {
568    THierarchicalCluster * cluster;
569    int left;
570    int right;
571
572    m_element(THierarchicalCluster * cluster, int left, int right);
573    m_element(const m_element & other);
574
575    inline bool operator< (const m_element & other) const;
576    inline bool operator== (const m_element & other) const;
577}; // cluster joined at left and right index
578
579struct ordering_element {
580    THierarchicalCluster * left;
581    unsigned int u; // the left most (outer) index of left luster
582    unsigned m; // the rightmost (inner) index of left cluster
583    THierarchicalCluster * right;
584    unsigned int w; // the right most (outer) index of the right cluster
585    unsigned int k; // the left most (inner) index of the right cluster
586
587    ordering_element();
588    ordering_element(THierarchicalCluster * left, unsigned int u, unsigned m,
589            THierarchicalCluster * right, unsigned int w, unsigned int k);
590    ordering_element(const ordering_element & other);
591};
592
593m_element::m_element(THierarchicalCluster * _cluster, int _left, int _right):
594    cluster(_cluster), left(_left), right(_right)
595{}
596
597m_element::m_element(const m_element & other):
598    cluster(other.cluster), left(other.left), right(other.right)
599{}
600
601bool m_element::operator< (const m_element & other) const
602{
603    if (cluster < other.cluster)
604        return true;
605    else
606        if (cluster == other.cluster)
607            if (left < other.left)
608                return true;
609            else
610                if (left == other.left)
611                    return right < other.right;
612                else
613                    return false;
614        else
615            return false;
616}
617
618bool m_element::operator== (const m_element & other) const
619{
620    return cluster == other.cluster && left == other.left && right == other.right;
621}
622
623struct m_element_hash
624{
625    inline size_t operator()(const m_element & m) const
626    {
627        size_t seed = 0;
628        hash_combine(seed, (size_t) m.cluster);
629        hash_combine(seed, (size_t) m.left);
630        hash_combine(seed, (size_t) m.right);
631        return seed;
632    }
633
634    // more or less taken from boost::hash_combine
635    inline void hash_combine(size_t &seed, size_t val) const
636    {
637        seed ^= val + 0x9e3779b9 + (seed << 6) + (seed >> 2);
638    }
639};
640
641ordering_element::ordering_element():
642    left(NULL), u(-1), m(-1),
643    right(NULL), w(-1), k(-1)
644{}
645
646ordering_element::ordering_element(THierarchicalCluster * _left,
647        unsigned int _u,
648        unsigned _m,
649        THierarchicalCluster * _right,
650        unsigned int _w,
651        unsigned int _k
652    ): left(_left), u(_u), m(_m),
653       right(_right), w(_w), k(_k)
654{}
655
656ordering_element::ordering_element(const ordering_element & other):
657       left(other.left), u(other.u), m(other.m),
658       right(other.right), w(other.w), k(other.k)
659{}
660
661
662
663typedef _HCLUST_UNORDERED_MAP<m_element, double, m_element_hash> join_scores;
664typedef _HCLUST_UNORDERED_MAP<m_element, ordering_element, m_element_hash> cluster_ordering;
665
666
667// Return the minimum distance between elements in matrix
668//
669float min_distance(
670        TIntList::iterator indices1_begin,
671        TIntList::iterator indices1_end,
672        TIntList::iterator indices2_begin,
673        TIntList::iterator indices2_end,
674        TSymMatrix & matrix)
675{
676//  TIntList::iterator iter1 = indices1.begin(), iter2 = indices2.begin();
677    float minimum = std::numeric_limits<float>::infinity();
678    TIntList::iterator indices2;
679    for (; indices1_begin != indices1_end; indices1_begin++)
680        for (indices2 = indices2_begin; indices2 != indices2_end; indices2++){
681            minimum = std::min(matrix.getitem(*indices1_begin, *indices2), minimum);
682        }
683    return minimum;
684}
685
686struct CompareByScores
687{
688    join_scores & scores;
689    const THierarchicalCluster & cluster;
690    const int & fixed;
691
692    CompareByScores(join_scores & _scores, const THierarchicalCluster & _cluster, const int & _fixed):
693        scores(_scores), cluster(_cluster), fixed(_fixed)
694    {}
695    bool operator() (int lhs, int rhs)
696    {
697        m_element left((THierarchicalCluster*)&cluster, fixed, lhs);
698        m_element right((THierarchicalCluster*)&cluster, fixed, rhs);
699        return scores[left] < scores[right];
700    }
701};
702
703
704//#include <iostream>
705//#include <cassert>
706
707// This needs to be called with all left, right pairs to
708// update all scores for cluster.
709void partial_opt_ordering(
710        THierarchicalCluster & cluster,
711        THierarchicalCluster & left,
712        THierarchicalCluster & right,
713        THierarchicalCluster & left_left,
714        THierarchicalCluster & left_right,
715        THierarchicalCluster & right_left,
716        THierarchicalCluster & right_right,
717        TSymMatrix &matrix,
718        join_scores & M,
719        cluster_ordering & ordering)
720{
721    int u = 0, w = 0;
722    TIntList & mapping = cluster.mapping.getReference();
723    for (TIntList::iterator u_iter = mapping.begin() + left_left.first;
724            u_iter != mapping.begin() + left_left.last;
725            u_iter++)
726        for (TIntList::iterator w_iter = mapping.begin() + right_right.first;
727                w_iter != mapping.begin() + right_right.last;
728                w_iter++)
729        {
730            u = *u_iter;
731            w = *w_iter;
732            float curr_min = std::numeric_limits<float>::infinity();
733            int curr_k = 0, curr_m = 0;
734            float C = min_distance(mapping.begin() + left_right.first,
735                    mapping.begin() + left_right.last,
736                    mapping.begin() + right_left.first,
737                    mapping.begin() + right_left.last,
738                    matrix);
739
740            vector<int> m_ordered(mapping.begin() + left_right.first,
741                    mapping.begin() + left_right.last);
742            vector<int> k_ordered(mapping.begin() + right_left.first,
743                    mapping.begin() + right_left.last);
744
745            // TODO: precompute the scores for m and k in an array and use a simpler
746            // comparison function
747            std::sort(m_ordered.begin(), m_ordered.end(), CompareByScores(M, left, u));
748            std::sort(k_ordered.begin(), k_ordered.end(), CompareByScores(M, right, w));
749
750
751            int k0 = k_ordered.front();
752            m_element m_right_k0(&right, w, k0);
753            int m = 0, k = 0;
754            for (vector<int>::iterator iter_m=m_ordered.begin(); iter_m != m_ordered.end(); iter_m++)
755            {
756                m = *iter_m;
757
758                m_element m_left(&left, u, m);
759
760                if (M[m_left] + M[m_right_k0] + C >= curr_min){
761                    break;
762                }
763                for (vector<int>::iterator iter_k = k_ordered.begin(); iter_k != k_ordered.end(); iter_k++)
764                {
765                    k = *iter_k;
766                    m_element m_right(&right, w, k);
767                    if (M[m_left] + M[m_right] + C >= curr_min)
768                    {
769                        break;
770                    }
771                    float test_val = M[m_left] + M[m_right] + matrix.getitem(m, k);
772                    if (curr_min > test_val)
773                    {
774                        curr_min = test_val;
775                        curr_k = k;
776                        curr_m = m;
777                    }
778                }
779
780            }
781
782            M[m_element(&cluster, u, w)] = curr_min;
783            M[m_element(&cluster, w, u)] = curr_min;
784
785//          assert(M[m_element(&cluster, u, w)] == M[m_element(&cluster, u, w)]);
786//          assert(M[m_element(&cluster, u, w)] == curr_min);
787
788//          assert(ordering.find(m_element(&cluster, u, w)) == ordering.end());
789//          assert(ordering.find(m_element(&cluster, w, u)) == ordering.end());
790
791            ordering[m_element(&cluster, u, w)] = ordering_element(&left, u, curr_m, &right, w, curr_k);
792            ordering[m_element(&cluster, w, u)] = ordering_element(&right, w, curr_k, &left, u, curr_m);
793        }
794}
795
796void order_clusters(
797        THierarchicalCluster & cluster,
798        TSymMatrix &matrix,
799        join_scores & M,
800        cluster_ordering & ordering,
801        TProgressCallback * callback)
802{
803    if (cluster.size() == 1)
804    {
805        M[m_element(&cluster, cluster.mapping->at(cluster.first), cluster.mapping->at(cluster.first))] = 0.0;
806        return;
807    }
808    else if (cluster.branches->size() == 2)
809    {
810        PHierarchicalCluster left = cluster.branches->at(0);
811        PHierarchicalCluster right = cluster.branches->at(1);
812
813        order_clusters(left.getReference(), matrix, M, ordering, callback);
814        order_clusters(right.getReference(), matrix, M, ordering, callback);
815
816        PHierarchicalCluster  left_left = (!left->branches) ? left : left->branches->at(0);
817        PHierarchicalCluster  left_right = (!left->branches) ? left : left->branches->at(1);
818        PHierarchicalCluster  right_left = (!right->branches) ? right : right->branches->at(0);
819        PHierarchicalCluster  right_right = (!right->branches) ? right : right->branches->at(1);
820
821        // 1.)
822        partial_opt_ordering(cluster,
823                left.getReference(), right.getReference(),
824                left_left.getReference(), left_right.getReference(),
825                right_left.getReference(), right_right.getReference(),
826                matrix, M, ordering);
827
828        if (right->branches)
829            // 2.) Switch right branches.
830            // (if there are no right branches the ordering has already been evaluated in 1.)
831            partial_opt_ordering(cluster,
832                    left.getReference(), right.getReference(),
833                    left_left.getReference(), left_right.getReference(),
834                    right_right.getReference(), right_left.getReference(),
835                    matrix, M, ordering);
836
837        if (left->branches)
838            // 3.) Switch left branches.
839            // (if there are no left branches the ordering has already been evaluated in 1. and 2.)
840            partial_opt_ordering(cluster,
841                    left.getReference(), right.getReference(),
842                    left_right.getReference(), left_left.getReference(),
843                    right_left.getReference(), right_right.getReference(),
844                    matrix, M, ordering);
845
846        if (left->branches && right->branches)
847            // 4.) Switch both branches.
848            partial_opt_ordering(cluster,
849                    left.getReference(), right.getReference(),
850                    left_right.getReference(), left_left.getReference(),
851                    right_right.getReference(), right_left.getReference(),
852                    matrix, M, ordering);
853    }
854    if (callback)
855        // TODO: count the number of already processed nodes.
856        callback->operator()(0.0, PHierarchicalCluster(&cluster));
857}
858
859/* Check if TIntList contains element.
860 */
861bool contains(TIntList::iterator iter_begin, TIntList::iterator iter_end, int element)
862{
863    return std::find(iter_begin, iter_end, element) != iter_end;
864}
865
866void optimal_swap(THierarchicalCluster * cluster, int u, int w, cluster_ordering & ordering)
867{
868    if (cluster->branches)
869    {
870        assert(ordering.find(m_element(cluster, u, w)) != ordering.end());
871        ordering_element ord = ordering[m_element(cluster, u, w)];
872
873        PHierarchicalCluster left_right = (ord.left->branches)? ord.left->branches->at(1) : PHierarchicalCluster(NULL);
874        PHierarchicalCluster right_left = (ord.right->branches)? ord.right->branches->at(0) : PHierarchicalCluster(NULL);
875
876        TIntList & mapping = cluster->mapping.getReference();
877        if (left_right && !contains(mapping.begin() + left_right->first,
878                mapping.begin() + left_right->last, ord.m))
879        {
880            assert(!contains(mapping.begin() + left_right->first,
881                             mapping.begin() + left_right->last, ord.m));
882            ord.left->swap();
883            left_right = ord.left->branches->at(1);
884            assert(contains(mapping.begin() + left_right->first,
885                            mapping.begin() + left_right->last, ord.m));
886        }
887        optimal_swap(ord.left, ord.u, ord.m, ordering);
888
889        assert(mapping.at(ord.left->first) == ord.u);
890        assert(mapping.at(ord.left->last - 1) == ord.m);
891
892        if (right_left && !contains(mapping.begin() + right_left->first,
893                mapping.begin() + right_left->last, ord.k))
894        {
895            assert(!contains(mapping.begin() + right_left->first,
896                             mapping.begin() + right_left->last, ord.k));
897            ord.right->swap();
898            right_left = ord.right->branches->at(0);
899            assert(contains(mapping.begin() + right_left->first,
900                            mapping.begin() + right_left->last, ord.k));
901        }
902        optimal_swap(ord.right, ord.k, ord.w, ordering);
903
904        assert(mapping.at(ord.right->first) == ord.k);
905        assert(mapping.at(ord.right->last - 1) == ord.w);
906
907        assert(mapping.at(cluster->first) == ord.u);
908        assert(mapping.at(cluster->last - 1) == ord.w);
909    }
910}
911
912PHierarchicalCluster THierarchicalClusterOrdering::operator() (
913        PHierarchicalCluster root,
914        PSymMatrix matrix)
915{
916    join_scores M; // scores
917    cluster_ordering ordering;
918    order_clusters(root.getReference(), matrix.getReference(), M, ordering,
919            progress_callback.getUnwrappedPtr());
920
921    int u = 0, w = 0;
922    int min_u = 0, min_w = 0;
923    float min_score = std::numeric_limits<float>::infinity();
924    TIntList & mapping = root->mapping.getReference();
925
926    for (TIntList::iterator u_iter = mapping.begin() + root->branches->at(0)->first;
927            u_iter != mapping.begin() + root->branches->at(0)->last;
928            u_iter++)
929        for (TIntList::iterator w_iter = mapping.begin() + root->branches->at(1)->first;
930                    w_iter != mapping.begin() + root->branches->at(1)->last;
931                    w_iter++)
932        {
933            u = *u_iter; w = *w_iter;
934            m_element el(root.getUnwrappedPtr(), u, w);
935            if (M[el] < min_score)
936            {
937                min_score = M[el];
938                min_u = u;
939                min_w = w;
940            }
941        }
942//  std::cout << "Min score "<< min_score << endl;
943
944    optimal_swap(root.getUnwrappedPtr(), min_u, min_w, ordering);
945    return root;
946}
Note: See TracBrowser for help on using the repository browser.