diff --git a/evaluation.py b/evaluation.py new file mode 100644 index 0000000000000000000000000000000000000000..12f59fc237b6954d5c7be0cae4a69ae45448edab --- /dev/null +++ b/evaluation.py @@ -0,0 +1,201 @@ +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)