source: orange/source/orange/hclust.cpp @ 10888:4296caed2ecd

Revision 10888:4296caed2ecd, 30.1 KB checked in by Ales Erjavec <ales.erjavec@…>, 2 years ago (diff)

Moved unordered_map include to the top (see #1188)

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