Add new file

evaluation.py

+from geo2net import *
+from matplotlib.pyplot import figure,imshow,plot,hist,show as pltshow, draw as pltdraw
+from skimage.future.graph import show_rag
+from skimage.exposure import rescale_intensity
+sys.path.append('C:\\Program Files\\QGIS 2.18\\apps\\Python27\\Lib\\site-packages')
+from osgeo import gdal
+import numpy as np
+import networkx as nx
+import os
+from utilities import multilayer_format_geo2net
+import matplotlib.pyplot as plt
+import glob
+import network_analysis
+from os.path import basename
+
+
+fld = './testdata'
+
+# cropmask
+ds = gdal.Open(os.path.join(fld,'koumbia_cropland.tif'))
+crop = ds.ReadAsArray()
+ds = None
+
+# seg singola immagine
+segfile = os.path.join(fld,'koumbia_seg.tif')
+ds = gdal.Open(segfile)
+seg = ds.ReadAsArray()
+ds = None
+
+
+communities_file = ".\\output\\infomap\\mlnet_K20_cmeans.clu"
+
+
+
+
+def readCommunities(clu_file):
+    coms = {}
+    f = open(clu_file, 'r')
+    for line in f:
+        if not line.startswith("#"):
+            vals = line.split(' ')
+            com = int(vals[1])
+            coms[int(vals[0])]=com
+    return coms
+
+def readCommunities_list(clu_file):
+    coms = {}
+    f = open(clu_file, 'r')
+    for line in f:
+        if not line.startswith("#"):
+            vals = line.split(' ')
+            com = int(vals[1])
+            if com not in coms:
+                coms[com]=[]
+            coms[com].append(int(vals[0]))
+    return coms
+
+def histo(coms,seg,crop):
+    #print(len(coms),seg.shape,crop.shape)
+    classes = {}
+    comset=set()
+    for i, j in np.ndindex(seg.shape): #for each pixel
+        if seg[i,j]!=2246:
+            curr_com = coms[seg[i,j]]
+            comset.add(curr_com)
+            if curr_com not in classes:
+                classes[curr_com]=np.zeros(2)
+            classes[curr_com][crop[i,j]-1]+=1
+
+    N = len(comset)
+
+    crop = []
+    nat = []
+    for c in classes:
+        crop.append(classes[c][0])
+        nat.append(classes[c][1])
+
+    ind = np.arange(N)  # the x locations for the groups
+    width = 0.35  # the width of the bars: can also be len(x) sequence
+
+    p1 = plt.bar(ind, crop, width,)
+    p2 = plt.bar(ind, nat, width,
+                 bottom=crop)
+
+    #plt.ylabel('Scores')
+    #plt.title('Scores by group and gender')
+    #plt.xticks(ind, ('C1', 'C2', 'C3', 'C4', 'C5'))
+    plt.yticks(np.arange(0, 81, 10))
+    plt.legend((p1[0], p2[0]), ('Cropland', 'Natural Areas'))
+    plt.savefig("hist"+str(len(classes))+".png")
+
+    #plt.show()
+
+def percs(coms, seg, crop):
+    classes = {}
+    comset = set()
+    for i, j in np.ndindex(seg.shape):  # for each pixel
+        if seg[i, j] != 2246:
+            curr_com = coms[seg[i, j]]
+            comset.add(curr_com)
+            if curr_com not in classes:
+                classes[curr_com] = np.zeros(2)
+            classes[curr_com][crop[i, j] - 1] += 1
+
+    for c in classes:
+        tot = np.sum(classes[c])
+        classes[c]/=tot
+        #print(np.max(classes[c]),c)
+    return classes
+
+
+# Implementazione multislice modularity definita in :
+#8. Mucha, P. J., Richardson, T., Macon, K., Porter, M. A. & Onnela, J.-P.
+# Community structure in time - dependent, multiscale, and multiplex networks.
+# Science 328, 876-8(2010)
+
+def multislice_modularity(layer_graphs,interlayer,communities):
+    gamma = 1.0 #teoricamente se ne puo definire uno per layer
+    w = 1.0 #interlayer weight: 1 per i coupling, 0 altrimenti [CHECK]
+    assignments = {}
+    c = 0
+    for D_key in communities:
+        D = communities[D_key]
+        for n in D:
+            assignments[n]=c
+        c+=1
+    mu = 0.0 #total edge weight
+    Q = 0.0
+    for l in layer_graphs:
+        for i in layer_graphs[l].nodes():
+            sum_i = 0.0
+            ni = layer_graphs[l].neighbors(i)
+            deg_i = 0
+            if ni is not None:
+                deg_i = len(list(ni))
+            for j in layer_graphs[l].neighbors(i):
+                if i in assignments and j in assignments:
+                    if assignments[i]==assignments[j]: #teoricamente gli assegn. possono essere diversi per ogni layer
+                        deg_j = len(list(layer_graphs[l].neighbors(j)))
+                        P_ij = float(deg_i*deg_j)/float(2*len(layer_graphs[l].edges()))
+                        #P_ij = float(deg_i-deg_j)/float(2*len(layer_graphs[l].edges()))
+                        sum_i+= 1.0-gamma*P_ij
+                        if i in interlayer:
+                            if j in interlayer[i]:
+                                sum_i += w # contribution interlayer edge
+            Q+=sum_i
+            mu+=deg_i+(len(layer_graphs.keys())-1)*w    # grado su ogni layer + peso coupling edges (uno per layer)
+    Q=float(Q)/float(2*mu)
+    return Q
+
+
+if __name__=='__main__':
+
+    coms = readCommunities(".\\output\\infomap\\mlnet_K10_cmeans_pc50.clu")
+    histo(coms,seg,crop)
+    exit()
+    pr = percs(coms, seg, crop)
+    maxes = []
+    for c in pr:
+        maxes.append(np.max(pr[c]))
+    print maxes
+    print len(maxes),np.max(maxes),np.min(maxes),np.mean(maxes),np.std(maxes)
+    coms = readCommunities_list(".\\output\\infomap\\mlnet_K10_cmeans_pc50.clu")
+    layer_graphs, node_layers, interlayer = network_analysis.readNet("D:\\Mes Donnees\\Amoris\\geo2net-master\\geo2net-master\\ml_nets\\mlnet_K10_cmeans_pc50.ncol")
+    m = multislice_modularity(layer_graphs, interlayer, coms)
+    print(m)
+    exit()
+
+    os.chdir("D:\\Mes Donnees\\Amoris\\geo2net-master\\geo2net-master\\ml_nets\\")
+
+    files = glob.glob("*.ncol")
+    for f in files:
+        name = basename(f).split('.')[0]
+        comfile = "D:\\Mes Donnees\\Amoris\\geo2net-master\\geo2net-master\\output\\infomap\\"+name+".clu"
+        coms = readCommunities_list(comfile)
+        layer_graphs, node_layers, interlayer = network_analysis.readNet(f)
+        m = multislice_modularity(layer_graphs,interlayer,coms)
+        print name,m
+
+    #files = glob.glob(".\\output\\infomap\\*.clu")
+
+    #for f in files:
+    #   coms = readCommunities(f)
+
+
+    #coms = readCommunities(communities_file)
+    #histo(coms,seg,crop)
+    #percs(coms,seg,crop)
+
+
+    files = glob.glob(".\\output\\infomap\\*.clu")
+
+    for f in files:
+        coms = readCommunities(f)
+        histo(coms,seg,crop)
+        pr = percs(coms, seg, crop)
+        maxes = []
+        for c in pr:
+            maxes.append(np.max(pr[c]))
+        print f,maxes
+        #print f,len(maxes),np.max(maxes),np.min(maxes),np.mean(maxes),np.std(maxes)