MDL network population clustering +++++++++ Source Code ------------ .. code-block:: python from collections import Counter import numpy as np from scipy.special import loggamma import random import time .. _generate-synthetic: .. code-block:: python def generate_synthetic(S,N,modes,alphas,betas,pis): """ generate synthetic networks from the heterogeneous population model uses fast binomial sampling for false positive edges """ def ind2ij(ind,N): i = N - 2 - np.floor(np.sqrt(-8*ind + 4*N*(N-1)-7)/2.0 - 0.5) j = ind + i + 1 - N*(N-1)/2 + (N-i)*((N-i)-1)/2 return int(i),int(j) K = len(modes) NC2 = int(N*(N-1)/2) nets,cluster_labels = [],[] for t in range(S): k = np.random.choice(range(K),p=pis) Mk = len(modes[k]) net = set() for e in modes[k]: if np.random.rand() < alphas[k]: net.add(e) num_fps = np.random.binomial(NC2-Mk,betas[k]) while num_fps > 0: ind = np.random.randint(NC2-1) i,j = ind2ij(ind,N) if not((i,j) in modes[k]): net.add((i,j)) num_fps -= 1 nets.append(net) cluster_labels.append(k) return nets,cluster_labels .. _remap-keys: .. code-block:: python def remap_keys(Dict): """ remap dict keys to first K integers """ sorted_keys = sorted(list(Dict.keys())) for i,u in enumerate(sorted_keys): Dict[i] = Dict.pop(u) return Dict .. _mdl-populations-init: .. code-block:: python class MDL_populations(): """ MDL population clustering class Inputs: edgesets: list of sets. the s-th set contains all the edges (i,j) in the s-th network in the sample (do not include the other direction (j,i) if network is undirected). the order of edgesets within D only matters for contiguous clustering, where we want the edgesets to be in order of the samples in time N: number of nodes in each network K0: initial number of clusters (for discontiguous clustering, usually K0 = 1 works well; for contiguous clustering it doesn't matter) n_fails: number of failed reassign/merge/split/merge-split moves before terminating algorithm bipartite: 'None' for unipartite network populations, array [# of nodes of type 1, # of nodes of type 2] otherwise directed: boolean indicating whether edgesets contain directed edges max_runs: maximum number of allowed moves, regardless of number of fails Outputs of 'run_sims' (unconstrained description length optimization) and 'dynamic_contiguous' (restriction to contiguous clusters): C: dictionary with items (cluster label):(set of indices corresponding to networks in cluster) A: dictionary with items (cluster label):(set of edges corresponding to mode of cluster) L: inverse compression ratio (description length after clustering)/(description length of naive transmission) """ def __init__(self, edgesets, N, K0 = 1, n_fails = 100, bipartite = None, directed = False, max_runs = np.inf): """ initialize class attributes """ self.edgesets = edgesets self.K0 = K0 self.n_fails = n_fails self.S = len(self.edgesets) self.N = N self.max_runs = max_runs if bipartite is not None: self.NC2 = bipartite[0]*bipartite[1] #bipartite networks only differentiated from unipartite ones through this term if directed: self.NC2 = self.N*(self.N-1) #directed networks only differentiated from undirected ones through this term else: self.NC2 = self.N*(self.N-1)/2 self.C,self.E,self.A = {},{},{} self.attmerges,self.attsplits,self.attmergesplits = set(),set(),set() .. _mdl-populations-initialize-clusters: .. code-block:: python def initialize_clusters(self): """ initialize K0 random clusters and find their modes as well as the total description length of this configuration """ for i,s in enumerate(np.random.permutation(range(self.S))): k = str(i % self.K0) if k in self.C: self.C[k].add(s) else: self.C[k] = set() self.C[k].add(s) for k in range(self.K0): k = str(k) self.E[k] = self.generate_Ek(self.C[k]) self.A[k] = self.update_mode(self.E[k],len(self.C[k])) self.L = sum([self.Lk(self.A[k],self.E[k],len(self.C[k])) for k in self.C]) .. _mdl-populations-random_key: .. code-block:: python def random_key(self): """generate random key for new cluster""" return str(np.random.randint(1000000000000)) .. _mdl-populations-logchoose: .. code-block:: python def logchoose(self,N,K): """logarithm of binomial coefficient""" return loggamma(N+1) - loggamma(N-K+1) - loggamma(K+1) .. _mdl-populations-logmult: .. code-block:: python def logmult(self,Ns): """logarithm of multinomial coefficient with denominator Ns[0]!Ns[1]!...""" return loggamma(sum(Ns)+1) - sum(loggamma(i+1) for i in Ns) .. _mdl-populations-generate-ek: .. code-block:: python def generate_Ek(self,cluster): """ tally edge counts for networks in cluster """ Ek = {} for s in cluster: for e in self.edgesets[s]: if e in Ek: Ek[e] += 1 else: Ek[e] = 1 return Ek .. _mdl-populations-update-mode: .. code-block:: python def update_mode(self,Ek,Sk): """ generate mode from cluster edge counts by greedily removing least common edges in cluster from mode of all in-cluster edges """ Ek_vals = set(Ek.values()) Acomplete = set(Ek.keys()) if (Sk == 1): return Acomplete # return network itself if cluster only has one network elif not(Ek): return set() # return empty mode if networks in cluster are empty elif (len(Ek_vals) == 1) and (next(iter(Ek_vals)) > 1): return Acomplete # return network itself if networks are all duplicates Etil = sorted(Ek.items(),key=lambda x:x[1]) r,tk,fk,Mk,Ak = 0,sum(Ek.values()),0,len(Ek),Acomplete.copy() Mmax = len(Ek) best_mode,deltaL,deltaL_best = 0,0,0 while (r < Mmax): e,Xij = Etil[r] Lafter = self.logchoose(self.NC2,Mk-1) + Sk*np.log(self.S/Sk) + self.logchoose(Sk*(Mk-1),tk-Xij) \ + self.logchoose(Sk*(self.NC2-(Mk-1)),fk+Xij) Lbefore = self.logchoose(self.NC2,Mk) + Sk*np.log(self.S/Sk) + self.logchoose(Sk*Mk,tk) + self.logchoose(Sk*(self.NC2-Mk),fk) Ak.discard(e) r += 1 tk -= Xij fk += Xij Mk -= 1 deltaL += (Lafter - Lbefore) if deltaL < deltaL_best: deltaL_best = deltaL best_mode = r Ak = Acomplete.copy() for r in range(best_mode): e,Xij = Etil[r] Ak.discard(e) return Ak .. _mdl-populations-lk: .. code-block:: python def Lk(self,Ak,Ek,Sk): """ cluster description length as function of mode, edge counts, and size of cluster """ if Sk == 0: return 0. Mk = len(Ak) tk,fk = 0,sum(Ek.values()) for e in Ak: if e in Ek: tk += Ek[e] fk -= Ek[e] return self.logchoose(self.NC2,Mk) + Sk*np.log(self.S/Sk) + self.logchoose(Sk*Mk,tk) + self.logchoose(Sk*(self.NC2-Mk),fk) .. _mdl-populations-move1: .. code-block:: python def move1(self,k=None): """ move type 1: reassign randomly chosen network to best cluster """ ks = list(self.C.keys()) if k is None: k = random.choice(ks) if len(self.C) == 1: return self.move3() #try splitting if only one cluster else: s = np.random.choice(list(self.C[k])) Ek_after = self.E[k].copy() for e in self.edgesets[s]: Ek_after[e] -= 1 if Ek_after[e] == 0: del Ek_after[e] deltaL1s = {} L_kbefore = self.Lk(self.A[k],self.E[k],len(self.C[k])) ks = list(self.C.keys()) for kp in set(ks) - set({k}): L_kpbefore = self.Lk(self.A[kp],self.E[kp],len(self.C[kp])) Ekp_after = self.E[kp].copy() for e in self.edgesets[s]: if e in Ekp_after: Ekp_after[e] += 1 else: Ekp_after[e] = 1 deltaL1s[kp] = self.Lk(self.A[kp],Ekp_after,len(self.C[kp])+1) \ + self.Lk(self.A[k],Ek_after,len(self.C[k])-1) - L_kbefore - L_kpbefore if min(deltaL1s.values()) < 0: min_kp = min(deltaL1s, key=deltaL1s.get) self.C[k].discard(s) self.C[kp].add(s) for e in self.edgesets[s]: self.E[k][e] -= 1 if self.E[k][e] == 0: self.E[k].pop(e, None) if e in self.E[kp]: self.E[kp][e] += 1 else: self.E[kp][e] = 1 knew,kpnew = self.random_key(),self.random_key() self.C[knew] = self.C.pop(k) self.C[kpnew] = self.C.pop(kp) self.E[knew] = self.E.pop(k) self.E[kpnew] = self.E.pop(kp) self.A[knew] = self.A.pop(k) self.A[kpnew] = self.A.pop(kp) self.A[knew] = self.update_mode(self.E[knew],len(self.C[knew])) self.A[kpnew] = self.update_mode(self.E[kpnew],len(self.C[kpnew])) if not(self.C[knew]): del self.C[knew] del self.A[knew] del self.E[knew] return True, deltaL1s[min_kp] else: return False, 0 .. _mdl-populations-move2: .. code-block:: python def move2(self): """ move type 2: merge two randomly chosen clusters """ if len(self.C) == 1: #try splitting if only one cluster return self.move3() ks = list(self.C.keys()) kp,kpp = np.random.choice(ks,size=2,replace=False) if ((kp,kpp) in self.attmerges) or ((kpp,kp) in self.attmerges): #check if merge already has been tried and failed return False,0 Ek = self.E[kp].copy() for e in self.E[kpp]: if e in Ek: Ek[e] += self.E[kpp][e] else: Ek[e] = self.E[kpp][e] Skp,Skpp = len(self.C[kp]),len(self.C[kpp]) Sk = Skp + Skpp Ak = self.update_mode(Ek,Sk) deltaL2 = self.Lk(Ak,Ek,Sk) - self.Lk(self.A[kp],self.E[kp],Skp) - self.Lk(self.A[kpp],self.E[kpp],Skpp) if deltaL2 < 0: k = self.random_key() self.C[k] = self.C[kp].union(self.C[kpp]) del self.C[kp] del self.C[kpp] self.E[k] = Ek.copy() del self.E[kp] del self.E[kpp] self.A[k] = Ak.copy() del self.A[kp] del self.A[kpp] return True, deltaL2 else: self.attmerges.add((kp,kpp)) #add to attempted merges if move fails return False, 0 .. _mdl-populations-move3: .. code-block:: python def move3(self): """ move type 3: split randomly chosen cluster in two and perform K-means type algorithm to get these clusters and modes """ ks = list(self.C.keys()) k = random.choice(ks) if len(self.C[k]) == 1: #if only one network in cluster, try move 1 with this cluster return self.move1(k) if k in self.attsplits: #if split already tried and failed, exit return False,0 Sk = len(self.C[k]) localC = {0:set(),1:set()} for i,s in enumerate(np.random.permutation(list(self.C[k]))): localC[i % 2].add(s) localS,localE,localA,localL = {0:None,1:None},{0:None,1:None},{0:None,1:None},{0:None,1:None} for kl in [0,1]: localS[kl] = len(localC[kl]) localE[kl] = self.generate_Ek(localC[kl]) localA[kl] = self.update_mode(localE[kl],localS[kl]) localL[kl] = self.Lk(localA[kl],localE[kl],localS[kl]) #local 2-means type algorithm for identifying clusters C[k] will split into movement = True num_iters,max_2means = 0,10 while (movement == True) and (num_iters < max_2means): movement = False to_move = [] localEafter = {} for s in self.C[k]: for kl in [0,1]: localEafter[kl] = localE[kl].copy() if s in localC[0]: old,new = 0,1 else: old,new = 1,0 for e in self.edgesets[s]: if e in localEafter[new]: localEafter[new][e] += 1 else: localEafter[new][e] = 1 localEafter[old][e] -= 1 if localEafter[old][e] == 0: del localEafter[old][e] deltaLmove = self.Lk(localA[new],localEafter[new],len(localC[new])+1) \ + self.Lk(localA[old],localEafter[old],len(localC[old])-1) \ - localL[0] - localL[1] if deltaLmove < 0: to_move.append((s,old,new)) movement = True for tup in to_move: s,old,new = tup localC[new].add(s) localC[old].discard(s) localS[new] += 1 localS[old] -= 1 if localS[old] == 0: return self.move2() for e in self.edgesets[s]: if e in localE[new]: localE[new][e] += 1 else: localE[new][e] = 1 localE[old][e] -= 1 if localE[old][e] == 0: del localE[old][e] for kl in [0,1]: localA[kl] = self.update_mode(localE[kl],localS[kl]) localL[kl] = self.Lk(localA[kl],localE[kl],localS[kl]) num_iters += 1 deltaL3 = self.Lk(localA[0],localE[0],localS[0]) + self.Lk(localA[1],localE[1],localS[1]) - self.Lk(self.A[k],self.E[k],Sk) if deltaL3 < 0: kp,kpp = self.random_key(),self.random_key() self.C[kp] = localC[0].copy() self.C[kpp] = localC[1].copy() del self.C[k] self.E[kp] = localE[0].copy() self.E[kpp] = localE[1].copy() del self.E[k] self.A[kp] = localA[0].copy() self.A[kpp] = localA[1].copy() del self.A[k] return True,deltaL3 else: self.attsplits.add(k) return False, 0 .. _mdl-populations-move4: .. code-block:: python def move4(self): """ move type 4: merge two randomly chosen clusters then split them (perform moves 2 and 3 in a row) """ if len(self.C) == 1: return self.move3() # try split move if only a single cluster exists ks = list(self.C.keys()) k1,k2 = np.random.choice(ks,size=2,replace=False) if ((k1,k2) in self.attmergesplits) or ((k2,k1) in self.attmergesplits): #check if merge-split combo already tried with these clusters return False,0 Ck = self.C[k1].union(self.C[k2]) Sk = len(Ck) localC = {0:set(),1:set()} for i,s in enumerate(np.random.permutation(list(Ck))): localC[i % 2].add(s) localS,localE,localA,localL = {0:None,1:None},{0:None,1:None},{0:None,1:None},{0:None,1:None} for kl in [0,1]: localS[kl] = len(localC[kl]) localE[kl] = self.generate_Ek(localC[kl]) localA[kl] = self.update_mode(localE[kl],localS[kl]) localL[kl] = self.Lk(localA[kl],localE[kl],localS[kl]) movement = True num_iters,max_2means = 0,10 while (movement == True) and (num_iters < max_2means): movement = False to_move = [] localEafter = {} for s in Ck: for kl in [0,1]: localEafter[kl] = localE[kl].copy() if s in localC[0]: old,new = 0,1 else: old,new = 1,0 for e in self.edgesets[s]: if e in localEafter[new]: localEafter[new][e] += 1 else: localEafter[new][e] = 1 localEafter[old][e] -= 1 if localEafter[old][e] == 0: del localEafter[old][e] deltaLmove = self.Lk(localA[new],localEafter[new],len(localC[new])+1) \ + self.Lk(localA[old],localEafter[old],len(localC[old])-1) \ - localL[0] - localL[1] if deltaLmove < 0: to_move.append((s,old,new)) movement = True for tup in to_move: s,old,new = tup localC[new].add(s) localC[old].discard(s) localS[new] += 1 localS[old] -= 1 if localS[old] == 0: #if all networks go into one cluster, try merge move instead return self.move2() for e in self.edgesets[s]: if e in localE[new]: localE[new][e] += 1 else: localE[new][e] = 1 localE[old][e] -= 1 if localE[old][e] == 0: del localE[old][e] for kl in [0,1]: localA[kl] = self.update_mode(localE[kl],localS[kl]) localL[kl] = self.Lk(localA[kl],localE[kl],localS[kl]) num_iters += 1 deltaL4 = self.Lk(localA[0],localE[0],localS[0]) + self.Lk(localA[1],localE[1],localS[1]) \ - self.Lk(self.A[k1],self.E[k1],len(self.C[k1])) - self.Lk(self.A[k2],self.E[k2],len(self.C[k2])) if deltaL4 < 0: kp,kpp = self.random_key(),self.random_key() self.C[kp] = localC[0].copy() self.C[kpp] = localC[1].copy() del self.C[k1] del self.C[k2] self.E[kp] = localE[0].copy() self.E[kpp] = localE[1].copy() del self.E[k1] del self.E[k2] self.A[kp] = localA[0].copy() self.A[kpp] = localA[1].copy() del self.A[k1] del self.A[k2] return True,deltaL4 else: self.attmergesplits.add((k1,k2)) return False, 0 .. _mdl-populations-run-sims: .. code-block:: python def run_sims(self): """ run discontiguous (unconstrained) merge split simulations to identify modes and clusters that minimize the description length """ nf,runs = 0,0 move_times,move_types = [],[] while (nf < self.n_fails) and (runs < self.max_runs): start = time.time() move = np.random.choice([1,2,3,4]) accepted,deltaL = eval('self.move'+str(move)+'()') if accepted: nf = 0 else: nf += 1 self.L += deltaL runs += 1 move_times.append(time.time()-start) move_types.append(move) M = sum([len(D) for D in self.edgesets]) self.L = sum([self.Lk(self.A[k],self.E[k],len(self.C[k])) for k in self.C]) self.L /= self.logchoose(self.S*self.NC2,M) #return (minimum description length)/(naive code length transmitting all networks separately) self.move_times = np.array(move_times) self.move_types = np.array(move_types) return remap_keys(self.C),remap_keys(self.A),self.L .. _mdl-populations-dynamic-contiguous: .. code-block:: python def dynamic_contiguous(self): """ miimize description length while constraining clusters to be contiguous in time (according to order of networks in 'edgesets') uses dynamic programming approach for exact optimization, and ignores cluster label entropy terms Sk*log(S/Sk) """ self.LMDL = {} self.clusters = {} self.modes = {} self.LMDL[-1] = 0 self.clusters[-1] = {} self.modes[-1] = {} self.LMDL[0] = self.logchoose(self.NC2,len(self.edgesets[0])) + np.log(self.S) key0 = self.random_key() self.clusters[0] = {key0:set([0])}.copy() self.modes[0] = {key0:self.edgesets[0].copy()}.copy() start = time.time() for j in range(1,self.S): jkey = self.random_key() Lj = self.LMDL[j-1] + self.logchoose(self.NC2,len(self.edgesets[j])) + np.log(self.S) Cj = self.clusters[j-1].copy() Cj[jkey] = set([j]) Aj = self.modes[j-1].copy() Aj[jkey] = self.edgesets[j].copy() localE = Counter(list(Aj[jkey])) for i in np.arange(j-1,-1,-1): Lprop = self.LMDL[i-1] Cprop = self.clusters[i-1].copy() Cprop[jkey] = set(range(i,j+1)) Aprop = self.modes[i-1].copy() for e in self.edgesets[i]: if e in localE: localE[e] += 1 else: localE[e] = 1 Aprop[jkey] = self.update_mode(localE,len(Cprop[jkey])) Lprop += self.Lk(Aprop[jkey],localE,len(Cprop[jkey])) - len(Cprop[jkey])*np.log(self.S/len(Cprop[jkey])) + np.log(self.S) if Lprop < Lj: Lj = Lprop Cj = Cprop.copy() Aj = Aprop.copy() self.LMDL[j] = Lj self.clusters[j] = Cj.copy() self.modes[j] = Aj.copy() self.C = self.clusters[self.S-1].copy() self.A = self.modes[self.S-1].copy() M = sum([len(D) for D in self.edgesets]) self.L = self.LMDL[self.S-1]/self.logchoose(self.S*self.NC2,M) self.runtime = time.time() - start return remap_keys(self.C),remap_keys(self.A),self.L .. _mdl-populations-evaluate-partition: .. code-block:: python def evaluate_partition(self,partition,contiguous=False): """ evaluate description length of partition. contiguous option removes cluster label entropy term from description length """ for s,k in enumerate(partition): if k in self.C: self.C[k].add(s) else: self.C[k] = set() self.C[k].add(s) K = len(self.C) for k in range(K): self.E[k] = self.generate_Ek(self.C[k]) self.A[k] = self.update_mode(self.E[k],len(self.C[k])) self.L = sum([self.Lk(self.A[k],self.E[k],len(self.C[k])) for k in self.C]) if contiguous: self.L -= sum([len(self.C[k])*np.log(self.S/len(self.C[k])) for k in self.C]) M = sum([len(D) for D in self.edgesets]) self.L /= self.logchoose(self.S*self.NC2,M) return self.L