Changeset 7846:ed4f33ee23ce in orange


Ignore:
Timestamp:
04/15/11 13:30:55 (3 years ago)
Author:
ales_erjavec <ales.erjavec@…>
Branch:
default
Convert:
90b83ddd1380d2ce245edbf5cdd7f21413c15078
Message:

Added utility functions (iterative traversal of clusters, dendrogram layout function ...).

File:
1 edited

Legend:

Unmodified
Added
Removed
  • orange/Orange/clustering/hierarchical.py

    r7809 r7846  
    1717""" 
    1818import orange 
     19import Orange 
    1920from Orange.core import HierarchicalClustering, \ 
    2021                        HierarchicalCluster, \ 
    2122                        HierarchicalClusterList 
     23                         
     24from Orange.misc import progressBarMilestones 
     25                         
     26import sys 
    2227 
    2328 
     
    8691        progressCallback --function used to report progress 
    8792    """ 
     93#    from Orange.misc import recursion_limit 
     94     
    8895    objects = getattr(tree.mapping, "objects", None) 
    8996    tree.mapping.setattr("objects", range(len(tree))) 
     
    9198    ordering = {} 
    9299    visitedClusters = set() 
    93     def _optOrdering(tree): 
     100     
     101    def _optOrderingRecursive(tree): 
    94102        if len(tree)==1: 
    95103            for leaf in tree: 
    96104                M[tree, leaf, leaf] = 0 
    97 ##                print "adding:", tree, leaf, leaf 
    98105        else: 
    99             _optOrdering(tree.left) 
    100             _optOrdering(tree.right) 
    101 ##            print "ordering", [i for i in tree] 
     106            _optOrderingRecursive(tree.left) 
     107            _optOrderingRecursive(tree.right) 
     108             
    102109            Vl = set(tree.left) 
    103110            Vr = set(tree.right) 
     
    110117            for u in Vl: 
    111118                for w in Vr: 
    112                     if True: #Improved search 
     119#                    if True: #Improved search 
     120                    C = min([matrix[m, k] for m in other(u, Vll, Vlr) for k in other(w, Vrl, Vrr)]) 
     121                    orderedMs = sorted(other(u, Vll, Vlr), key=lambda m: M[tree_left, u, m]) 
     122                    orderedKs = sorted(other(w, Vrl, Vrr), key=lambda k: M[tree_right, w, k]) 
     123                    k0 = orderedKs[0] 
     124                    curMin = 1e30000  
     125                    curM = curK = None 
     126                    for m in orderedMs: 
     127                        if M[tree_left, u, m] + M[tree_right, w, k0] + C >= curMin: 
     128                            break 
     129                        for k in  orderedKs: 
     130                            if M[tree_left, u, m] + M[tree_right, w, k] + C >= curMin: 
     131                                break 
     132                            testMin = M[tree_left, u, m] + M[tree_right, w, k] + matrix[m, k] 
     133                            if curMin > testMin: 
     134                                curMin = testMin 
     135                                curM = m 
     136                                curK = k 
     137                    M[tree, u, w] = M[tree, w, u] = curMin 
     138                    ordering[tree, u, w] = (tree_left, u, curM, tree_right, w, curK) 
     139                    ordering[tree, w, u] = (tree_right, w, curK, tree_left, u, curM) 
     140#                    else: 
     141#                    def MFunc((m, k)): 
     142#                        return M[tree_left, u, m] + M[tree_right, w, k] + matrix[m, k] 
     143#                    m, k = min([(m, k) for m in other(u, Vll, Vlr) for k in other(w, Vrl, Vrr)], key=MFunc) 
     144#                    M[tree, u, w] = M[tree, w, u] = MFunc((m, k)) 
     145#                    ordering[tree, u, w] = (tree_left, u, m, tree_right, w, k) 
     146#                    ordering[tree, w, u] = (tree_right, w, k, tree_left, u, m) 
     147 
     148            if progressCallback: 
     149                progressCallback(100.0 * len(visitedClusters) / len(tree.mapping)) 
     150                visitedClusters.add(tree) 
     151     
     152    if False: 
     153        with recursion_limit(sys.getrecursionlimit() + len(tree)): 
     154            _optOrderingRecursive(tree) 
     155         
     156    def _optOrderingIterative(tree): 
     157        if len(tree)==1: 
     158            for leaf in tree: 
     159                M[tree, leaf, leaf] = 0 
     160        else: 
     161#            _optOrdering(tree.left) 
     162#            _optOrdering(tree.right) 
     163             
     164            Vl = set(tree.left) 
     165            Vr = set(tree.right) 
     166            Vlr = set(tree.left.right or tree.left) 
     167            Vll = set(tree.left.left or tree.left) 
     168            Vrr = set(tree.right.right or tree.right) 
     169            Vrl = set(tree.right.left or tree.right) 
     170            other = lambda e, V1, V2: V2 if e in V1 else V1 
     171            tree_left, tree_right = tree.left, tree.right 
     172            for u in Vl: 
     173                for w in Vr: 
     174#                    if True: #Improved search 
    113175                        C = min([matrix[m, k] for m in other(u, Vll, Vlr) for k in other(w, Vrl, Vrr)]) 
    114176                        orderedMs = sorted(other(u, Vll, Vlr), key=lambda m: M[tree_left, u, m]) 
     
    131193                        ordering[tree, u, w] = (tree_left, u, curM, tree_right, w, curK) 
    132194                        ordering[tree, w, u] = (tree_right, w, curK, tree_left, u, curM) 
    133                     else: 
    134                         def MFunc((m, k)): 
    135                             return M[tree_left, u, m] + M[tree_right, w, k] + matrix[m, k] 
    136                         m, k = min([(m, k) for m in other(u, Vll, Vlr) for k in other(w, Vrl, Vrr)], key=MFunc) 
    137                         M[tree, u, w] = M[tree, w, u] = MFunc((m, k)) 
    138                         ordering[tree, u, w] = (tree_left, u, m, tree_right, w, k) 
    139                         ordering[tree, w, u] = (tree_right, w, k, tree_left, u, m) 
    140  
    141             if progressCallback: 
    142                 progressCallback(100.0 * len(visitedClusters) / len(tree.mapping)) 
    143                 visitedClusters.add(tree) 
    144          
    145     _optOrdering(tree) 
    146  
    147     def _order(tree, u, w): 
    148         if len(tree)==1: 
    149             return 
    150         left, u, m, right, w, k = ordering[tree, u, w] 
    151         if len(left)>1 and m not in left.right: 
    152             left.swap() 
    153         _order(left, u, m) 
    154 ##        if u!=left[0] or m!=left[-1]: 
    155 ##            print "error 4:", u, m, list(left) 
    156         if len(right)>1 and k not in right.left: 
    157             right.swap() 
    158         _order(right, k, w) 
    159 ##        if k!=right[0] or w!=right[-1]: 
    160 ##            print "error 5:", k, w, list(right) 
     195#                    else: 
     196#                        def MFunc((m, k)): 
     197#                            return M[tree_left, u, m] + M[tree_right, w, k] + matrix[m, k] 
     198#                        m, k = min([(m, k) for m in other(u, Vll, Vlr) for k in other(w, Vrl, Vrr)], key=MFunc) 
     199#                        M[tree, u, w] = M[tree, w, u] = MFunc((m, k)) 
     200#                        ordering[tree, u, w] = (tree_left, u, m, tree_right, w, k) 
     201#                        ordering[tree, w, u] = (tree_right, w, k, tree_left, u, m) 
     202 
     203#            if progressCallback: 
     204#                progressCallback(100.0 * len(visitedClusters) / len(tree.mapping)) 
     205#                visitedClusters.add(tree) 
     206    from Orange.misc import progressBarMilestones 
     207     
     208    subtrees = postorder(tree) 
     209    milestones = progressBarMilestones(len(subtrees), 1000) 
     210     
     211    for i, subtree in enumerate(subtrees): 
     212        _optOrderingIterative(subtree) 
     213        if progressCallback and i in milestones: 
     214            progressCallback(100.0 * i / len(subtrees)) 
     215 
     216#    def _orderRecursive(tree, u, w): 
     217#        """ Order the tree based on the computed optimal ordering.  
     218#        """ 
     219#        if len(tree)==1: 
     220#            return 
     221#        left, u, m, right, w, k = ordering[tree, u, w] 
     222#        if len(left)>1 and m not in left.right: 
     223#            left.swap() 
     224#        _orderRecursive(left, u, m) 
     225#         
     226#        if len(right)>1 and k not in right.left: 
     227#            right.swap() 
     228#        _orderRecursive(right, k, w) 
     229         
     230    def _orderIterative(tree, u, w): 
     231        """ Order the tree based on the computed optimal ordering.  
     232        """ 
     233        opt_uw = {tree: (u, w)} 
     234        for subtree in preorder(tree): 
     235            if subtree.branches: 
     236                u, w = opt_uw[subtree] 
     237                left, u, m, right, w, k = ordering[subtree, u, w] 
     238                opt_uw[left] = (u, m) 
     239                opt_uw[right] = (k, w) 
     240                 
     241                if left.branches and m not in left.right: 
     242                    left.swap() 
     243                 
     244                if right.branches and k not in right.left: 
     245                    right.swap() 
     246             
     247#        if len(tree) == 1: 
     248#            return 
     249#        left, u, m, right, w, k = ordering[tree, u, w] 
     250#        if len(left) > 1 and m not in left.right: 
     251#            left.swap() 
     252#        _orderIterative(left, u, m) 
     253#         
     254#        if len(right) > 1 and k not in right.left: 
     255#            right.swap() 
     256#        _orderIterative(right, k, w) 
    161257     
    162258    u, w = min([(u, w) for u in tree.left for w in tree.right], key=lambda (u, w): M[tree, u, w]) 
     
    164260##    print "M(v) =", M[tree, u, w] 
    165261     
    166     _order(tree, u, w) 
    167  
    168     def _check(tree, u, w): 
    169         if len(tree)==1: 
    170             return 
    171         left, u, m, right, w, k = ordering[tree, u, w] 
    172         if tree[0] == u and tree[-1] == w: 
    173             _check(left, u, m) 
    174             _check(right, k, w) 
    175         else: 
    176             print "Error:", u, w, tree[0], tree[-1] 
    177  
    178     _check(tree, u ,w) 
    179      
     262    if False: 
     263        with recursion_limit(sys.getrecursionlimit() + len(tree)): 
     264            _orderRecursive(tree, u, w) 
     265             
     266    _orderIterative(tree, u, w) 
     267             
     268 
     269#    def _check(tree, u, w): 
     270#        if len(tree)==1: 
     271#            return 
     272#        left, u, m, right, w, k = ordering[tree, u, w] 
     273#        if tree[0] == u and tree[-1] == w: 
     274#            _check(left, u, m) 
     275#            _check(right, k, w) 
     276#        else: 
     277#            print "Error:", u, w, tree[0], tree[-1] 
     278# 
     279#    with recursion_limit(sys.getrecursionlimit() + len(tree)): 
     280#        _check(tree, u ,w) 
    180281 
    181282    if objects: 
     
    563664    d.plot(filename) 
    564665     
     666     
     667""" 
     668Utility functions 
     669================= 
     670 
     671""" 
     672     
     673def postorder(cluster): 
     674    """ Return a post order list of clusters. 
     675     
     676    :param cluster: Cluster 
     677    :type cluster: :class:`HierarchicalCluster` 
     678     
     679    """ 
     680    order = [] 
     681    visited = set() 
     682    stack = [cluster] 
     683    while stack: 
     684        cluster = stack.pop(0) 
     685         
     686        if cluster.branches: 
     687            if cluster in visited: 
     688                order.append(cluster) 
     689            else: 
     690                stack = cluster.branches + [cluster] + stack 
     691                visited.add(cluster) 
     692        else: 
     693            order.append(cluster) 
     694            visited.add(cluster) 
     695    return order 
     696     
     697     
     698def preorder(cluster): 
     699    """ Return a pre order list of clusters. 
     700     
     701    :param cluster: Cluster 
     702    :type cluster: :class:`HierarchicalCluster` 
     703     
     704    """ 
     705    order = [] 
     706    stack = [cluster] 
     707    while stack: 
     708        cluster = stack.pop(0) 
     709        order.append(cluster) 
     710        if cluster.branches: 
     711            stack = cluster.branches + stack 
     712    return order 
     713     
     714     
     715def dendrogram_layout(root_cluster, expand_leaves=False): 
     716    """ Return a layout of the cluster dendrogram on a 2D plane. The return  
     717    value if a list of (subcluster, (start, center, end)) tuples where 
     718    subcluster is an instance of :class:`HierarchicalCluster` and start, 
     719    end are the two end points for the cluster branch. The list is sorted 
     720    in post-order. 
     721     
     722    :param root_cluster: Cluster to layout. 
     723    :type root_cluster: :class:`HierarchicalCluster` 
     724     
     725    :param expand_leaves: If True leaves will span proportional to the number 
     726        of items they map, else all leaves will be the same width.  
     727      
     728    """ 
     729    result = [] 
     730    cluster_geometry = {} 
     731    leaf_idx = 0 
     732    for cluster in postorder(root_cluster): 
     733        if not cluster.branches: 
     734            if expand_leaves: 
     735                start, end = float(cluster.first), float(cluster.last) 
     736            else: 
     737                start = end = leaf_idx 
     738                leaf_idx += 1 
     739            center = (start + end) / 2.0 
     740            cluster_geometry[cluster] = (start, center, end) 
     741            result.append((cluster, (start, center, end))) 
     742        else: 
     743            left = cluster.branches[0] 
     744            right = cluster.branches[1] 
     745            left_center = cluster_geometry[left][1] 
     746            right_center = cluster_geometry[right][1] 
     747            start, end = left_center, right_center 
     748            center = (start + end) / 2.0 
     749            cluster_geometry[cluster] = (start, center, end) 
     750            result.append((cluster, (start, center, end))) 
     751             
     752    return result 
     753     
     754     
     755def pruned(root_cluster, level=None, height=None, condition=None): 
     756    """ Return a new pruned clustering instance. 
     757     
     758    .. note:: This uses `copy.deepcopy` to create a copy of the root_cluster 
     759        instance. 
     760     
     761    :param cluster: Cluster to prune. 
     762    :type cluster: :class:`HierarchicalCluster` 
     763     
     764    :param level: If not None prune all clusters deeper then `level`. 
     765    :type level: int 
     766     
     767    :param height: If not None prune all clusters with height lower then 
     768        `height`. 
     769    :type height: float 
     770     
     771    :param condition: If not None condition must be a function taking a 
     772        single :class:`HierarchicalCluster` instance and returning a True  
     773        or False value indicating if the cluster should be pruned. 
     774    :type condition: function  
     775     
     776    """ 
     777    import copy 
     778     
     779    # XXX This is unsafe HierarchicalCluster should take care of copying 
     780    if hasattr(root_cluster.mapping, "objects"): 
     781        objects = root_cluster.mapping.objects 
     782        root_cluster.mapping.objects = None 
     783        has_objects = True 
     784    else: 
     785        has_objects = False 
     786         
     787    root_cluster = copy.deepcopy(root_cluster) 
     788     
     789    prune(root_cluster, level, height, condition) 
     790     
     791    if has_objects: 
     792        root_cluster.mapping.objects = objects 
     793         
     794    return root_cluster 
     795     
     796     
     797def prune(root_cluster, level=None, height=None, condition=None): 
     798    """ Prune clustering instance in place 
     799     
     800    :param cluster: Cluster to prune. 
     801    :type cluster: :class:`HierarchicalCluster` 
     802     
     803    :param level: If not None prune all clusters deeper then `level`. 
     804    :type level: int 
     805     
     806    :param height: If not None prune all clusters with height lower then 
     807        `height`. 
     808    :type height: float 
     809     
     810    :param condition: If not None condition must be a function taking a 
     811        single :class:`HierarchicalCluster` instance and returning a True  
     812        or False value indicating if the cluster should be pruned. 
     813    :type condition: function 
     814     
     815    """ 
     816    if not any(arg is not None for arg in [level, height, condition]): 
     817        raise ValueError("At least one pruning argument must be supplied") 
     818     
     819    level_check = height_check = condition_check = lambda cl: False 
     820    cluster_depth = cluster_depths(root_cluster) 
     821     
     822    if level is not None: 
     823        level_check = lambda cl: cluster_depth[cl] > level 
     824         
     825    if height is not None: 
     826        height_check = lambda cl: cl.height < height 
     827 
     828    if condition is not None: 
     829        condition_check = condition 
     830         
     831    pruned_set = set() 
     832     
     833    def check_all(cl): 
     834        return any([check(cl) for check in [level_check, height_check, 
     835                                            condition_check]]) 
     836         
     837    for cluster in preorder(root_cluster): 
     838        if cluster not in pruned_set: 
     839            if check_all(cluster): 
     840                cluster.branches = None 
     841                pruned_set.update(set(preorder(cluster))) 
     842            else: 
     843                pass 
     844     
     845     
     846def cluster_depths(root_cluster): 
     847    """ Return a dictionary mapping :class:`HierarchicalCluster` instances to 
     848    their depths in the `root_cluster` hierarchy. 
     849     
     850    """ 
     851    depths = {} 
     852    depths[root_cluster] = 0 
     853    for cluster in preorder(root_cluster): 
     854        cl_depth = depths[cluster] 
     855        if cluster.branches: 
     856            depths.update(dict.fromkeys(cluster.branches, cl_depth + 1)) 
     857    return depths 
     858 
     859 
     860def instance_distance_matrix(data, 
     861            distance_constructor=orange.ExamplesDistanceConstructor_Euclidean, 
     862            progress_callback=None): 
     863    """ A helper function that computes an :class:`Orange.core.SymMatrix` of all 
     864    pairwise distances between instances in `data`. 
     865     
     866    :param data: A data table 
     867    :type data: :class:`Orange.data.Table` 
     868     
     869    :param distance_constructor: An ExamplesDistance_Constructor instance. 
     870    :type distance_constructor: :class:`Orange.distances.ExampleDistConstructor` 
     871     
     872    """ 
     873    matrix = orange.SymMatrix(len(data)) 
     874    dist = distance_constructor(data) 
     875     
     876    iter_count = matrix.dim * (matrix.dim - 1) / 2 
     877    milestones = progressBarMilestones(iter_count, 100) 
     878     
     879    for count, ((i, ex1), (j, ex2)) in enumerate(_pairs(enumerate(data))): 
     880        matrix[i, j] = dist(ex1, ex2) 
     881        if progress_callback and count in milestones: 
     882            progress_callback(100.0 * count / iter_count) 
     883             
     884    return matrix  
     885 
     886 
     887def feature_distance_matrix(data, distance=None, progress_callback=None): 
     888    """ A helper function that computes an :class:`Orange.core.SymMatrix` of 
     889    all pairwise distances between features in `data`. 
     890     
     891    :param data: A data table 
     892    :type data: :class:`Orange.data.Table` 
     893    :param distance: a function taking two lists and returning the distance. 
     894    :type distance: function 
     895      
     896    """ 
     897    attributes = data.domain.attributes 
     898    matrix = orange.SymMatrix(len(attributes)) 
     899    iter_count = matrix.dim * (matrix.dim - 1) / 2 
     900    milestones = progressBarMilestones(iter_count, 100) 
     901     
     902    for count, ((i, a1), (j, a2)) in enumerate(_pairs(enumerate(attributes))): 
     903        matrix[i, j] = (1.0 - orange.PearsonCorrelation(a1, a2, data, 0).r) / 2.0 
     904        if progress_callback and count in milestones: 
     905            progress_callback(100.0 * count / iter_count) 
     906             
     907    return matrix 
     908 
     909 
     910def _pairs(seq, same = False): 
     911    """ Return all pairs from elements of `seq`. 
     912    """ 
     913    seq = list(seq) 
     914    same = 0 if same else 1 
     915    for i in range(len(seq)): 
     916        for j in range(i + same, len(seq)): 
     917            yield seq[i], seq[j] 
     918     
     919     
     920def joining_cluster(root_cluster, item1, item2): 
     921    """ Return the cluster where `item1` and `item2` are first joined 
     922     
     923    :param root_cluster: Clustering. 
     924    :type root_cluster: :class:`HierarchicalCluster` 
     925    :param item1: Index of the item or an element of `root_cluster.objects` 
     926    :param item2: Index of the item or an element of `root_cluster.objects` 
     927     
     928    """ 
     929    cluster = root_cluster 
     930    while True: 
     931        if cluster.branches: 
     932            for branch in cluster.branches: 
     933                if item1 in branch and item2 in branch: 
     934                    cluster = branch 
     935                    break 
     936            else: 
     937                return cluster 
     938        else: 
     939            return cluster 
     940         
     941 
     942def cophenetic_distances(root_cluster): 
     943    """ Return the cophenetic distance matrix between items in clustering. 
     944    Cophenetic distance is defined as the height of the cluster where the  
     945    two items are first joined. 
     946     
     947    :param root_cluster: Clustering. 
     948    :type root_cluster: :class:`HierarchicalCluster` 
     949      
     950    """ 
     951 
     952    mapping = root_cluster.mapping   
     953    matrix = Orange.core.SymMatrix(len(mapping)) 
     954    for cluster in postorder(root_cluster): 
     955        if cluster.branches: 
     956            for branch1, branch2 in _pairs(cluster.branches): 
     957                for idx1 in mapping[branch1.first: branch1.last]: 
     958                    for idx2 in mapping[branch2.first: branch2.last]: 
     959                        matrix[idx1, idx2] = cluster.height 
     960                 
     961        else: 
     962            for ind1, ind2 in _pairs(mapping[cluster.first: cluster.last]): 
     963                matrix[ind1, ind2] = cluster.height 
     964     
     965    return matrix 
     966     
     967def cophenetic_correlation(root_cluster, matrix): 
     968    """ Return the `cophenetic correlation coefficient 
     969    <http://en.wikipedia.org/wiki/Cophenetic_correlation>`_ of the given 
     970    clustering. 
     971     
     972    :param root_cluster: Clustering. 
     973    :type root_cluster: :class:`HierarchicalCluster` 
     974     
     975    :param matrix: The distance matrix from which the `root_cluster` was 
     976        derived. 
     977      
     978    """ 
     979    import numpy 
     980    cophenetic = cophenetic_distances(root_cluster) 
     981    cophenetic = [list(row) for row in cophenetic] 
     982    original = [list(row) for row in matrix] 
     983    cophenetic = numpy.ravel(cophenetic) 
     984    original = numpy.ravel(original) 
     985    return numpy.corrcoef(cophenetic, original)[0, 1] 
     986     
     987     
    565988if __name__=="__main__": 
    566989    data = orange.ExampleTable("doc//datasets//brown-selected.tab") 
Note: See TracChangeset for help on using the changeset viewer.