Source code for finishing

"""
OpenEnsembles is a resource for performing and analyzing ensemble clustering

Copyright (C) 2017 Naegle Lab

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
										   
"""
import numpy as np
import pandas as pd
import pprint as pp
import copy
import matplotlib.pyplot as plt
import sklearn.cluster
import sys
import os
import networkx as nx
from collections import defaultdict

[docs]class mixture_model: """ Implementation of the article Mixture Models for Ensemble CLustering Topchy, Jain, and Punch, "A mixture model for clustering ensembles Proc. SIAM Int. Conf. Data Mining (2004) Written by Pedro da Silva Tavares and adapted by Kristen M. Naegle Parameters ---------- parg: list of lists Solutions of assignments of objects to clusters across an ensemble N: int Number of data points to cluster nEnsCluster: int Number of clusters to create in mixture model Default is 2 iterations: Number of expectation maximization iterations Default is 10 Attributes ---------- labels: list of ints Final Mixture Model partitions of objects See Also -------- openensembles.cluster.mixture_model References ---------- """ def __init__(self, parg, N, nEnsCluster=2, iterations=10): self.parg = parg #list of lists of solutions self.N = N# number of data points self.nEnsCluster = nEnsCluster #number of clusters to make from ensemble self.iterations = iterations self.y = self.gatherPartitions() self.y self.K = self.genKj() self.alpha, self.v, self.ExpZ = self.initParameters() self.labels = [] self.piFinishing = {}
[docs] def gatherPartitions(self): ''' Returns the y vector. parg: list of H-labeling solutions nElem: number of features/objects ''' H = len(self.parg) listParts = np.concatenate(self.parg).reshape(H,self.N) #print listParts[:,0] y = [] [y.append(listParts[:,i]) for i in range(self.N)] y = pd.DataFrame(y, columns= np.arange(H)) y.index.name = 'objs' y.columns.name = 'partition' return y
[docs] def genKj(self): ''' Generates the K(j) H-array that contains the tuples of unique clusters of each j-th partition, eg: K = [(X,Y), (A,B)] ''' #K = np.zeros(y.shape[1], dtype= int) K = [] aux = [] for i in range(self.y.shape[1]): if 'NaN' in np.unique(self.y.iloc[:,i].values): aux = copy.copy(self.y.iloc[:,i].values) aux = [x for x in aux if x != 'NaN'] K.append(aux) else: K.append(tuple(np.unique(self.y.iloc[:,i].values))) return K
[docs] def initParameters(self): ''' The function initializes the parameters of the mixture model. ''' def initAlpha(): return np.ones(self.nEnsCluster) / self.nEnsCluster def initV(): v = [] [v.append([]) for j in range(self.y.shape[1])] #[v[j].append(list(np.ones(len(self.K[j])) / len(self.K[j]))) for j in range(self.y.shape[1]) for m in range(self.nEnsCluster)] for j in range(self.y.shape[1]): for m in range(self.nEnsCluster): aux = abs(np.random.randn(len(self.K[j]))) v[j].append( aux / sum(aux) ) return v def initExpZ(): return np.zeros(self.y.shape[0] * self.nEnsCluster).reshape(self.y.shape[0],self.nEnsCluster) alpha = initAlpha() v = initV() ExpZ = initExpZ() return alpha, v, ExpZ
[docs] def expectation(self): ''' Compute the Expectation (ExpZ) according to parameters. Obs: y(N,H) Kj(H) alpha(M) v(H,M,K(j)) ExpZ(N,M) ''' def sigma(a,b): return 1 if a == b else 0 M = self.ExpZ.shape[1] nElem = self.y.shape[0] for m in range(M): for i in range(nElem): prod1 = 1 num = 0 for j in range(self.y.shape[1]): ix1 = 0 for k in self.K[j]: prod1 *= ( self.v[j][m][ix1] ** sigma(self.y.iloc[i][j],k) ) ix1 += 1 num += self.alpha[m] * prod1 den = 0 for n in range(M): prod2 = 1 for j2 in range(self.y.shape[1]): ix2 = 0 for k in self.K[j2]: prod2 *= ( self.v[j2][n][ix2] ** sigma(self.y.iloc[i][j2],k) ) ix2 += 1 den += self.alpha[n] * prod2 self.ExpZ[i][m] = float(num) / float(den) return self.ExpZ
[docs] def maximization(self): ''' Update the parameters taking into account the ExpZ computed in the Expectation (ExpZ) step. Obs: y(N,H) Kj(H) alpha(M) v(H,M,K(j)) ExpZ(N,M) ''' def vecSigma(vec, k): ''' Compare i-th elements of vector to k assigining to a vector 1 if i-th == k, 0 otherwise. ''' aux = [] for i in vec: if i == k: aux.append(1) else: aux.append(0) return np.asarray(aux) def updateAlpha(): for m in range(self.alpha.shape[0]): self.alpha[m] = float(sum(self.ExpZ[:,m])) / float(sum(sum(self.ExpZ))) return self.alpha def updateV(): for j in range(self.y.shape[1]): for m in range(self.alpha.shape[0]): ix = 0 for k in self.K[j]: num = sum(vecSigma(self.y.iloc[:,j],k) * np.array(self.ExpZ[:,m])) den = 0 for kx in self.K[j]: den += sum(vecSigma(self.y.iloc[:,j],kx) * np.asarray(self.ExpZ[:,m])) self.v[j][m][ix] = float(num) / float(den) ix += 1 return self.v self.alpha = updateAlpha() self.v = updateV() return self.alpha, self.v
def emProcess(self): def piConsensus(): ''' The function outputs the final ensemble solution based on ExpZ values. ''' maxExpZValues = {} piFinishing = {} labels = [] for i in range(self.ExpZ.shape[0]): maxExpZValues[i] = max(self.ExpZ[i,:]) piFinishing[i] = [] for j in range(self.ExpZ.shape[1]): if (maxExpZValues[i] == self.ExpZ[i,j]): piFinishing[i].append(j + 1) labels.append(j+1) # choose randomly in the case of same values of ExpZ[i,:] #[piFinishing[i].delete(random.choice(piFinishing[i])) for i in piFinishing.keys() if (len(piFinishing[i]) > 1)] return piFinishing, labels i = 0 while(i<self.iterations): self.ExpZ = self.expectation() self.alpha, self.v = self.maximization() i += 1 piFinishing, labels = piConsensus() self.piFinishing = piFinishing self.labels = np.asarray(labels) # return piFinishing, labels return self.labels
[docs]class co_occurrence_linkage: """ Returns a final solution to the ensemble that is the agglomerative clustering of the co_occurrence matrix according to the linkage passed by linkage (default is Average). Parameters ---------- co_occ_object: openensembles.coMat object The co-occurrence object to operate on threshold: float The threshold to cut the linkage at to create hard partitions linkage: string linkage type. See `scipy.cluster.hierarchy <https://docs.scipy.org/doc/scipy/reference/generated/scipy.cluster.hierarchy.linkage.html#scipy.cluster.hierarchy.linkage>`_ Returns ------- labels: list of ints The final solution See Also -------- openensembles.cluster.finish_co_occ_linkage() """ def __init__(self, co_occ_object, threshold, linkage='average'): self.coMat = co_occ_object self.N = len(self.coMat.co_matrix) self.labels = np.empty(self.N) self.K= 0 #number of clusters made by the cut, to be replaced self.linkage = linkage self.threshold = threshold
[docs] def finish(self): """ This is the function that is called ona co_occurrence_linkage object that agglomeratively clusters the co_occurrence_matrix (self.coMat.co_matrix) According to self.linkage (linkage parameter set in initialization of object) with clusters equal to self.K (also set in intialization) """ #first get linkage, then cut lnk = self.coMat.link(linkage=self.linkage) labels = self.coMat.cut(lnk, self.threshold) self.K = len(np.unique(labels)) self.labels = labels return self.labels
[docs]class majority_vote: """ Based on Ana Fred's 2001 paper, this algorithm assingns clusters to the same class if they co-cluster at least 50% of the time. It greedily joins clusters with the evidence that at least one pair of items from two different clusters co-cluster a majority of the time. Outliers will get their own cluster. The original algorithm chose 0.5 as the threshold for a majority vote, but here you can override the default parameter. Parameters ---------- co_occ_object: openensembles.coMat object The co-occurrence object to operate on threshold: float The threshold, fraction of times any pair of objects clusters together that is considred for majority. For example threshold=0.5 means 50% of the time or more. Attributes ---------- K: int Number of clusters created labels: list of ints Solution of cluster assignments References ---------- Fred, Ana. “Finding Consistent Clusters in Data Partitions.” In Multiple Classifier Systems, edited by Josef Kittler and Fabio Roli, LNCS 2096., 309–18. Springer, 2001. See Also -------- openensembles.cluster.finish_majority_vote() """ def __init__(self, co_occ_matrix, threshold=0.5): self.co_matrix = co_occ_matrix self.K= 0 #number of clusters made self.N = len(co_occ_matrix) self.labels = np.empty(self.N) self.threshold=threshold
[docs] def finish(self): """ Finish the ensemble by majority vote """ labels = np.zeros(self.N).astype(int) currCluster = 1 x = self.co_matrix.as_matrix() for i in range(0, self.N): for j in range(i+1, self.N): if x[i,j] > self.threshold: # the clusters should be set to the same value and if both belong to an existing cluster, all # members of the clusters should be joined if labels[i] and labels[j]: cluster_num = min(labels[i] , labels[j]) cluster_toChange = max(labels[i] , labels[j]) indices = [k for k, x in enumerate(labels) if x == cluster_toChange] labels[indices] = cluster_num elif not labels[i] and not labels[j]: #a new cluster labels[i] = currCluster labels[j] = currCluster currCluster += 1 else: #one of them is in a cluster and one is not, one will be assigned the same thing, but saves an if cluster_num = max(labels[i] , labels[j]) labels[i] = cluster_num labels[j] = cluster_num else: #else don't join them and give them an assignment the first time they are traversed if not labels[i]: labels[i] = currCluster currCluster += 1 if not labels[j]: labels[j] = currCluster currCluster += 1 #Add a cleanup state, if there are any cluster numbers that jump, move them down clusters = np.sort(np.unique(labels)) for ind in range(0,len(clusters)-1): #operating on ind+1 if clusters[ind+1] != clusters[ind]+1: cluster_num = clusters[ind] + 1 #print("updating cluster num %d to %d")%(clusters[ind+1], cluster_num) indices = [k for k, x in enumerate(labels) if x == clusters[ind+1]] labels[indices] = cluster_num clusters[ind+1] = cluster_num self.labels = labels
[docs]class graph_closure: """ Returns a final solution of the ensemble based on treating the co-occurrence matrix as a weighted graph whose solution is found from identifying network components within the graph Finds k-percolated cliques in G, e.g, Unless the cliques argument evaluates to True, this algorithm first enumerates all cliques in G. These are stored in memory, which in large graphs can consume large amounts of memory. Returns a generator object. To return a list of percolated k-cliques, Notes ----- References ---------- * Based on the method outlined in Palla et. al., Nature 435, 814-818 (2005) * Based on Code for `Percolation From ConradLee <https://gist.github.com/conradlee/1341985>`_ See Also -------- openensembles.cluster.finish_graph_closure() """ def __init__(self, co_occ_matrix, threshold, clique_size = 3): self.co_matrix = co_occ_matrix self.K= 0 #number of clusters made self.N = len(co_occ_matrix) self.labels = np.empty(self.N) self.threshold = threshold self.coMat_binary = np.array(self.co_matrix >= threshold).astype(int) self.clique_size = clique_size
[docs] def finish(self): """ Finishes the ensemble by taking a binary adjacency matrix, defined in initilization according to the threshold given and percolates the cliques """ # From ConradLee on GitHUB def get_percolated_cliques(G, k): perc_graph = nx.Graph() cliques = [frozenset(c) for c in nx.find_cliques(G) if len(c) >= k] perc_graph.add_nodes_from(cliques) # First index which nodes are in which cliques membership_dict = defaultdict(list) for clique in cliques: for node in clique: membership_dict[node].append(clique) # For each clique, see which adjacent cliques percolate for clique in cliques: for adj_clique in get_adjacent_cliques(clique, membership_dict): if len(clique.intersection(adj_clique)) >= (k - 1): perc_graph.add_edge(clique, adj_clique) # Connected components of clique graph with perc edges # are the percolated cliques for component in nx.connected_components(perc_graph): yield(frozenset.union(*component)) def get_adjacent_cliques(clique, membership_dict): adjacent_cliques = set() for n in clique: for adj_clique in membership_dict[n]: if clique != adj_clique: adjacent_cliques.add(adj_clique) return adjacent_cliques G = nx.from_numpy_matrix(self.coMat_binary) y = get_percolated_cliques(G, self.clique_size) z = list(y) clusterNum = 0 while z: l = list(z.pop()) self.labels[l] = int(clusterNum) clusterNum+=1 self.labels = self.labels.astype(int) return self.labels