major refactor including analysis scripts

PSD/analysis.py

+__author__ = 'paul'
+
+import numpy as np
+import pandas
+import os, fnmatch
+
+try:
+    from pyunicorn.core import resistive_network
+except:
+    print "ERROR: pyunicorn not availiable."
+    exit(1)
+
+class ResNet(resistive_network.ResNetwork):
+    def __init__(self, netfile):
+        import networkx as nx
+        from scipy.sparse import dok_matrix
+        from re import findall
+
+        G = nx.read_graphml(netfile)
+        mapping = dict(zip(sorted(G.nodes()),range(G.number_of_nodes())))
+        G = nx.relabel_nodes(G, mapping, copy=False)
+
+        weights = dok_matrix(nx.to_numpy_matrix(G))
+        for key in weights.keys():
+            weights[key] = self._get_link_length(G, key)
+
+        self.nodes = np.array(G.nodes(data=False))
+        self.run = int(findall(r"\d+", netfile)[0])
+
+
+        super(ResNet, self).__init__(resistances=weights.todense())
+
+    def _get_link_length(self, G, link):
+        """
+        return euclidean distance between x and y
+        """
+        x = np.array([G.node[link[0]]['lat'], G.node[link[0]]['lon']])
+        y = np.array([G.node[link[1]]['lat'], G.node[link[1]]['lon']])
+        return np.sqrt(sum((x-y)**2))
+
+class BasicMeasures(object):
+    def __init__(self):
+        self.data = pandas.DataFrame()
+
+    def load(self, rn):
+        assert isinstance(rn, ResNet)
+        self.net = rn
+
+        if not self.data.empty:
+            # add space for new network
+            for row in xrange(rn.N):
+                self.data.loc[row + rn.run * rn.N] = np.repeat(np.NaN, len(self.data.columns))
+
+        self.add(rn.nodes, "idx")
+        self.add(np.repeat(rn.run, rn.N), "run")
+
+        # initialise with saved resilience measures
+        self.add(self._get_mu_b(), "mu_b(2)")
+        self.add(self._get_mu_b(model="fourth"), "mu_b(4)")
+        self.add(self._get_mu_s(), "mu_s(2)")
+        self.add(self._get_mu_s(model="fourth"), "mu_s(4)")
+
+        # add a set of standard unweighted measures
+        self.add(rn.degree(), "deg")
+        self.add(rn.average_neighbors_degree(), "avn_deg")
+        self.add(rn.betweenness() / rn.N, "spb[N]")
+
+        # add a set of standard weighted measures
+        self.add(self._get_min_clust(), "minC")
+        self.add(self.net.admittive_degree(), "strength")
+        self.add(self.net.average_neighbors_admittive_degree(), "avn_strength")
+        # use different norm to compare with spb
+        self.add([1. * self.net.vertex_current_flow_betweenness(a) * ((self.net.N - 1) / 2)  for a in xrange(self.net.N)], "vcfb[N]")
+        self.add([self.net.effective_resistance_closeness_centrality(a) for a in xrange(self.net.N)], "ercc")
+
+
+        # add motif information
+        self.add(np.array([int(i) for i in np.array(rn.degree()) == 1]), "dead_end")
+        dead = np.zeros(rn.N)
+        for i, b in enumerate(rn.betweenness()):
+            if np.isclose(b, [rn.N - 2, 2 * rn.N - 5, 2 * rn.N - 6, 3 * rn.N - 10, 4 * rn.N - 17, 5 * rn.N - 26], atol=0.01).any():
+                dead[i] = 1
+            else:
+                dead[i] = 0
+        self.add(dead, "dead_root")
+
+    def add(self, measure, label):
+        measure = np.array(measure)
+        #assert measure.shape == (self.net.N, )
+        if not label in self.data:
+            self.data[label] = measure
+        else:
+            self.data[label][self.net.N * self.net.run: self.net.N * (self.net.run + 1)] = measure
+
+    def report(self):
+        print self.data.describe()
+
+    def report_to_csv(self, filename="aggregated_data.csv"):
+        self.data.to_csv(filename)
+
+    def _get_mu_b(self, model="swing", threshold=5):
+        conv = np.load("random_net" + str(self.net.run) + "_" + model + "_converged.npy")
+        conv /= 1. * np.load("random_net" + str(self.net.run) + "_" + model + "_total.npy")
+        freq_filter = np.load("random_net" + str(self.net.run) + "_freq_filter.npy")
+        index = (np.abs(freq_filter - threshold)).argmin()
+        return conv[index, :]
+
+    def _get_mu_s(self, model="swing", threshold=5):
+        surv = np.load("random_net" + str(self.net.run) + "_" + model + "_survived.npy")
+        surv /= 1. * np.load("random_net" + str(self.net.run) + "_" + model + "_total.npy")
+        freq_filter = np.load("random_net" + str(self.net.run) + "_freq_filter.npy")
+        index = (np.abs(freq_filter - threshold)).argmin()
+        return surv[index, :]
+
+    def _snbs_predictor(self, threshold=0.15):
+
+        '''predicted probability to have poor basin stability'''
+        # %Coefficients:
+        # %              Estimate Std. Error z value Pr(>|z|)
+        # %(Intercept)   3.325612   0.143404   23.19   <2e-16 ***
+        # %ED            0.088743   0.003832   23.16   <2e-16 ***
+        # %ANED         -0.348762   0.011514  -30.29   <2e-16 ***
+        # %minC        -10.389121   0.271047  -38.33   <2e-16 ***
+        # %VCFB         -0.107782   0.007994  -13.48   <2e-16 ***
+        # %ERCC         -1.515226   0.033785  -44.85   <2e-16 ***
+        # %dead          4.925139   0.084272   58.44   <2e-16 ***
+
+        g = 3.325612 + \
+            0.088743 * self.data["strength"] - \
+            0.348762 * self.data["avn_strength"] -\
+            10.389121 * self.data["minC"] -\
+            0.107782 * self.data["vcfb[N]"] -\
+            1.515226 * self.data["ercc"] + \
+            4.925139 * self.data["dead_root"]
+        prob = 1.0 / (1.0 + np.exp(- g))
+        t = threshold
+        poor_bs = np.zeros(self.net.N)
+        for i in xrange(self.net.N):
+            if prob[i]>t:
+                poor_bs[i]=1
+        return prob
+
+    def _get_min_clust(self):
+        ''' W is the admittance matrix'''
+        W = self.net.get_admittance()
+        N = len(W)
+        C = np.zeros(N)
+        for i in xrange(N):
+            norm = 0
+            for j in xrange(N):
+                for k in xrange(j):
+                    if j!=k:
+                        norm += W[i,k]*W[i,j]
+                        if W[i,j]*W[i,k]*W[j,k]>0:
+                            C[i] += min( W[i,k],W[j,k] ) * min( W[i,j],W[j,k] )
+                            #print i,norm[i],C[i]
+            C[i] /=  max(norm,1)
+            #print C[i]
+        return C
+
+
+def main():
+
+    net = BasicMeasures()
+
+
+    for count in xrange(50):
+        print count
+        rn = ResNet(netfile="random_net" + str(count) + ".gml")
+        net.load(rn)
+
+    # add secondary data
+    net.add(net._snbs_predictor(), "prob_poor")
+
+    net.report()
+    net.report_to_csv()
+
+
+
+if __name__ == "__main__":
+    main()
+

PSD/node_classification.py

+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jul 28 15:40:53 2016
+
+@author: jannitzbon
+
+This is an implementation of a node classification algorithm which aims at
+systematically classifying nodes in  tree-like structures of networks.
+"""
+
+from __future__ import print_function
+
+import numpy as np
+import networkx as nx
+
+
+# functions
+
+# auxiliary function which returns the full branch of a node x given its children C
+def branch(x, C):
+    b = [x]
+    for c in C[x]:
+        for i in branch(c, C):
+            b.append(i)
+    return b
+
+# procedure which does the full node classification of a networkx graph G
+def full_node_classification(G):
+    # return values
+    Gs = [] # 0 list of all graphs for each iteration step
+
+    N = []  # 1 list of non-root nodes in trees
+    R = []  # 2 list of root nodes
+    X = []  # 3 list of nodes not in tree-shaped parts
+    
+    B = {}  # 4 dict of branches for all nodes
+    C = {}  # 5 dict of children for all nodes
+    p = {}  # 6 dict of parents for all tree-nodes and roots
+    
+    delta = {}  # 7 dict of depth levels for all tree-nodes and roots
+    eta = {}    # 8 dict of heights for all nodes
+    beta = {}   # 9 dict of branch sizes for all nodes
+    
+    # special case for tree graphs with odd diameter
+    maxIter = np.infty
+    if nx.is_tree(G):
+        diam = nx.diameter(G)
+        if diam % 2 == 1:
+            maxIter = (diam-1)/2
+
+
+    #### Identification of non-root nodes, parents, children, branches, and depths  
+    Gs.append( G.copy() )
+    # list of nodes for each level
+    Vs = []
+    Vs.append( sorted(G.nodes()) )
+    # list of leaves for each level
+    Ds = []
+    # initialize branch, children and depth dicts
+    B = { x:[] for x in Vs[0] }
+    C = { x:[] for x in Vs[0] }
+    # iteration step    
+    lvl = 0
+    while True:
+        # find all leaves of current level
+        Ds.append( sorted( [ x for x, deg in nx.degree(Gs[lvl], Vs[lvl]).items() if deg==1] ) )
+        # check if leaves not empty
+        if (len(Ds[lvl]) == 0) or (lvl == maxIter):
+            break
+        # continue if leaves present
+        else:
+            # define nodes and graphs of next level
+            Vs.append( sorted( [ x for x in Vs[lvl] if x not in Ds[lvl] ] ) )
+            Gs.append( nx.subgraph(Gs[lvl], nbunch=Vs[lvl+1]))
+            # calculate further measures
+            for x in Ds[lvl]:
+                # add leaves to non-root nodes
+                N.append(x)
+                # add leaves to parent's list of children
+                p[x] = nx.neighbors(Gs[lvl], x)[0]
+                C[p[x]].append(x)
+                # determine branch for each removed leave
+                B[x] = branch(x,C)
+                # save depth of each leave
+                delta[x] = lvl
+            # increase level counter
+            lvl+=1
+    
+    #### Identification of root nodes
+    #determine all root nodes
+    R = list( set( [ p[x] for x in N if (p[x] not in N) ] ) )
+    #determine branches and depth levels of roots
+    for r in R:
+        # determine branch for roots
+        B[r] = branch(r,C)
+        # calculate depth of root
+        delta[r] = 1 + max( [ delta[x] for x in C[r] ] )
+             
+    #calculate branch sizes for all nodes
+    beta = { x:len(branch) for x,branch in B.items() }
+    
+    #### Identification of heights (this implementation is still a bit clumsy)
+    Hs = []
+    Hs.append(R)   
+    Ws = []
+    Ws.append(R)    
+    i = 0 
+    while len(Hs[i]) != 0:
+        Hn = []
+        for x in Hs[i]:
+            for c in C[x]:
+                if c not in Ws[i]:
+                    Hn.append(c)
+            # save the height value
+            eta[x]=i
+        Hs.append( Hn )
+        
+        Wn = []
+        for l in Ws:
+            for x in l:
+                Wn.append(x)
+        Ws.append( Wn )
+    
+        i+=1
+        
+    #### Identification of non-Tree nodes
+    for v in Vs[0]:
+        if (v not in R) and (v not in N):
+            X.append(v)
+            
+    
+    # return all determined values
+    return Gs, N, R, X, B, C, p, delta, eta, beta
+
+# procedure which returns a dict of node categories indexed by the node number
+def node_categories(G, denseThres=6.):
+    Gs, N, R, X, B, C, p, delta, eta, beta = full_node_classification(G)
+    
+    avndeg = nx.average_neighbor_degree(G)
+    cat = {}
+    
+    for x in G.nodes():
+         if x in X: cat[x]="bulk"
+         if x in R: cat[x]="root"
+         if x in N: 
+             cat[x]="inner"
+             #leaves
+             if delta[x]==0:
+                 cat[x]="leaf"
+                 #sprouts
+                 if eta[x]==1:
+                     cat[x]="sparseSprout"
+                     if avndeg[x]>=denseThres:
+                         cat[x]="denseSprout"
+    
+    return cat
+    
+#%%
+if __name__ == '__main__':    
+    
+    # generate connected random geometric graph Gtest
+    while True:
+        Gtest = nx.random_geometric_graph(50, 0.2)
+        if nx.is_connected(Gtest):
+            break
+    pos = nx.get_node_attributes(Gtest, 'pos')
+    
+    # generate minimum spanning tree Gmst
+    Gmst = nx.minimum_spanning_tree(Gtest)
+    nx.draw_networkx(Gmst, pos)
+    
+    #%%  do the full root classification of Gmst
+    # add some arbitrary edge such that we have a cylce
+    Gmst.add_edge(35,41)
+    nx.draw_networkx(Gmst, pos)
+    #Gs, N, R, X, B, C, p, delta, eta, beta = full_node_classification(Gmst)
+    
+    cats=node_categories(Gmst, denseThres=6)
+    
+
+
+
+
+