Source code for graphs

#!/usr/bin/env python3
#  Digraph3 graphs.py module
#  Python3.3+ computing resources
#  Copyright (C)  2011-2015 Raymond Bisdorff
#
#    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, write to the Free Software Foundation, Inc.,
#    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
#############################################
from decimal import Decimal

[docs]class Graph(object): """ In the `graphs` module, the root :py:class:`graphs.Graph` class provides a generic graph model. A given object consists in: 1. a vertices dictionary 2. a characteristic valuation domain, {-1,0,+1} by default 3. an edges dictionary, characterising each edge in the given valuation domain 4. a gamma function dictionary, holding the neighborhood vertices of each vertex General structure:: vertices = {'v1': {'name': ...,'shortName': ...}, 'v2': {'name': ...,'shortName': ...}, 'v3': {'name': ...,'shortName': ...}, ... } valuationDomain = {'min': -1, 'med': 0, 'max': 1} edges = {frozenset({'v1','v2'}): 1, frozenset({'v1','v3'}): 1, frozenset({'v2','v3'}): -1, ...} ## links from each vertex to its neighbors gamma = {'v1': {'v2',v3'}, 'v2': {'v1'}, 'v3': {'v1'}, ... } Example python3 session: >>> from graphs import Graph >>> g = Graph(numberOfVertices=5,edgeProbability=0.5) # random instance >>> g.showShort() *----- show short --------------* *---- short description of the graph ----* Name : 'random' Vertices : ['v1', 'v2', 'v3', 'v4', 'v5'] Valuation domain : {'med': 0, 'max': 1, 'min': -1} Gamma function : v1 -> ['v4'] v2 -> [] v3 -> ['v4'] v4 -> ['v1', 'v3'] v5 -> [] """ def __init__(self, fileName=None, Empty=False, numberOfVertices=7, edgeProbability=0.5): """ Constructor for Graph objects. """ from decimal import Decimal if Empty: self.name = 'emptyInstance' self.vertices = dict() self.order = len(self.vertices) self.edges = dict() self.valuationDomain = {'min': Decimal('-1'), 'med': Decimal('0'), 'max': Decimal('1')} self.gamma = dict() self.size = 0 elif fileName==None: g = RandomGraph(order=numberOfVertices,\ edgeProbability=edgeProbability) self.name = g.name self.vertices = g.vertices self.order = len(self.vertices) self.edges = g.edges self.valuationDomain = g.valuationDomain self.size = self.computeSize() self.gamma = self.gammaSets() else: fileNameExt = fileName+'.py' argDict = {} exec(compile(open(fileNameExt).read(), fileNameExt, 'exec'),argDict) self.name = fileName self.vertices = argDict['vertices'] self.order = len(self.vertices) self.valuationDomain = argDict['valuationDomain'] self.edges = argDict['edges'] self.size = self.computeSize() self.gamma = self.gammaSets() def __neg__(self): """ Make the negation operator -self available for Graph instances. Returns a DualGraph instance of self. """ new = DualGraph(self) new.__class__ = self.__class__ return new #-----------Dias/Castonguay/Longo/Jradi--------* def _degreeLabelling(self): """ Inspired from Dias, Castonguay, Longo & Jradi, Algorithmica 2015, p 14 """ degree = {} color = {} for v in self.vertices: degree[v] = len(self.gamma[v]) color[v] = 'white' labelling = {} for i in range(1,self.order+1): minDegree = self.order for x in self.vertices: if color[x] == 'white' and degree[x] < minDegree: v = x minDegree = degree[x] labelling[v] = i color[v] = 'black' #print(v,i,minDegree) for u in self.gamma[v]: if color[u] == 'white': degree[u] -= 1 self.labelling = labelling return labelling def _triplets(self,Comments=False): """ p.15 Inspired from Dias, Castonguay, Longo & Jradi, Algorithmica 2015. """ from itertools import product labelling = self.labelling tG = [] cycles = set() for u in self.vertices: for x,y in product(self.gamma[u],repeat=2): if x != y: ## print(u,self.labelling[u], ## x,self.labelling[x], ## y,self.labelling[y]) if labelling[u] < labelling[x] and \ labelling[x] < labelling[y]: ## if u < x and x < y: if self.edges[frozenset([x,y])] < self.valuationDomain['med']: if Comments: print('inital triple:',x,u,y) tG.append((x,u,y)) else: if Comments: print('3-cycle:',x,u,y) cycles.add((x,u,y)) #self.tG = tG #self.cycles = cycles return tG,cycles ## def _computeChordlessCyclesMP(self,Odd=False,\ ## Threading=False,nbrOfCPUs=None,\ ## Comments=False,Debug=False): ## """ ## Multiprocessing version of computeChordlessCycles(). ## ## Renders the set of all chordless odd cycless detected in an undirrected graph. ## Result (possible empty list) stored in <self.cyclesList> ## holding a possibly empty list tuples with at position 0 the ## list of adjacent actions of the circuit and at position 1 ## the set of actions in the stored circuit. ## Inspired by Dias, Castonguay, Longo, Jradi, Algorithmica (2015). ## ## Returns a possibly empty list of tuples (cycle,frozenset(cycle)). ## ## If Odd == True, only cycless of odd length are retained in the result. ## """ ## ## tG,self.cyclesList = self._triplets(Comments=Comments) ## if Comments: ## print('There are %d starting triplets !' % len(tG) ) ## self.blocked = {} ## if Threading: ## self.Odd = Odd ## self.Comments = Comments ## from multiprocessing import Pool ## from os import cpu_count ## if nbrOfCPUs == None: ## nbrOfCPUs= cpu_count() ## with Pool(nbrOfCPUs) as proc: ## cycless = proc.map(self._computeChordlessPathsFromInitialTriplet,tG) ## #print(circuits) ## for i in range(len(tG)): ## if Debug: ## print(i,cycles[i]) ## if cycles[i] != []: ## for c in cycles[i]: ## #print(circ) ## self.cyclesList.append(c) ## else: ## for p in tG: ## u = p[1] #### if Debug: #### print('===>>>',p,u) ## gammaU = (self.gamma[u][1] | self.gamma[u][0]) ## for x in gammaU: ## #print(x) ## self.blocked[x] += 1 ## self.cyclesList.append(self._ccVisit(p,u, ## Odd=Odd, ## Comments=Comments)) ## for x in gammaU: ## if self.blocked[x] > 0: ## self.blocked[x] -= 1 ## if Debug: ## print(self.cyclessList) ## return self.cyclessList ## ## def _computeChordlessPathsFromInitialTriplet(self,p,Debug=False): ## if self.Comments: ## print('===>> thread : ',p) ## u = p[1] ## blocked = {} ## for x in self.actions: ## blocked[x] = 0 ## circuits = [] ## gammaU = (self.gamma[u][1] | self.gamma[u][0]) ## for x in gammaU: ## blocked[x] += 1 ## cycless,blocked = self._ccVisitMP(cycless,blocked,p,u,Odd=self.Odd,Comments=self.Comments) ## for x in gammaU: ## if blocked[x] > 0: ## blocked[x] -= 1 ## if self.Comments: ## print(p,cycless) ## for x in self.actions: ## if blocked[x] > 1: ## blocked[x] = 0 ## if Debug: ## print(p,'return',c<cless) ## return cycless ## ## def _ccVisitMP(self,cyless,blocked,p,u, ## Odd=False,Comments=False,Debug=False): ## """ p.15 """ ## Med = self.valuationdomain['med'] ## ut = p[-1] ## u1 = p[0] ## inAsymGammaUt = self.gamma[ut][1] - self.gamma[ut][0] ## gammaUt = self.gamma[ut][0] | self.gamma[ut][1] ## if Debug: ## print(p,self.gamma[ut][1],ut,self.gamma[ut][0]) ## for x in gammaUt: ## blocked[x] += 1 ## for v in inAsymGammaUt: ## if str(v) > str(u) and blocked[v] == 1: ## p1 = p + tuple([v]) ## if Debug: ## print(p,p1) ## if self.relation[u1][v] > Med and\ ## self.relation[v][u1] <= Med: ## if Odd: ## if (len(p1) % 2) != 1: ## OddFlag=False ## else: ## OddFlag = True ## else: ## OddFlag = True ## if OddFlag: ## circ = list(reversed(p1)) ## if Comments: ## print(p,'cycle certificate: ',circ) ## circuits.append((circ,frozenset(circ))) ## ## elif self.relation[u1][v] <= Med and\ ## self.relation[v][u1] <= Med : ## if Debug: ## print(p,'continue with ', p1) ## cycles,blocked = self._ccVisitMP(cycles,blocked, ## p1,u,Odd=Odd, ## Comments=Comments) #### circuits.append(circuits1) ## if Debug: ## print(p,cycles) ## for x in (gammaUt): ## if blocked[x] > 0: ## blocked[x] -= 1 ## ## return cycles,blocked
[docs] def computeChordlessCycles(self,Cycle3=False,Comments=False,Debug=False): """ Renders the set of all chordless cycles observed in a Graph intance. Inspired from Dias, Castonguay, Longo & Jradi, Algorithmica 2015. .. note:: By default, a chordless cycle must have at least length 4.If the Cycle3 flag is set to True, the cyclicly closed triplets will be inserted as 3-cycles in the result. """ if Debug: Comments=True #self.visitedChordlessPathsNew = [] self._degreeLabelling() triplets,cycles3 = self._triplets(Comments=Comments) if Comments: print('# of initial triplets:',len(triplets)) print('# of 3-cycles :',len(cycles3)) if Cycle3: cycles = cycles3 else: cycles = set() self.blocked = {} for u in self.vertices: self.blocked[u] = 0 for p in triplets: if Comments: print(p,self.blocked) if Comments: print('===>>>', p) u = p[1] for x in self.gamma[u]: self.blocked[x] += 1 cycles = self._ccVisit(p,cycles,u,Comments=Comments) for x in self.gamma[u]: if self.blocked[x] > 0: self.blocked[x] -= 1 return cycles
def _ccVisit(self,p,cycles,u,Comments=False): """ p.15 """ #labelling = self.labelling #self.visitedChordlessPathsNew.append(p) ut = p[-1] u1 = p[0] ## if Comments: ## print(p,u1,u,ut) for x in self.gamma[ut]: self.blocked[x] += 1 for v in self.gamma[ut]: if self.labelling[v] > self.labelling[u] and self.blocked[v] == 1: p1 = p + tuple([v]) if self.edges[frozenset([u1,v])] > self.valuationDomain['med']: if Comments: print('Cycle certificate: ', p1) cycles.add(p1) else: if Comments: print('continue ...',p1) cycles = self._ccVisit(p1,cycles,u,Comments=Comments) for x in self.gamma[ut]: if self.blocked[x] > 0: self.blocked[x] -= 1 return cycles #---------------------------------------- def _chordlessPaths(self,Pk,v0,Cycle3=False,Comments=False,Debug=False): """ recursice chordless precycle (len > 3) construction: Pk is the current pre chordless cycle v0 is the initial vertex of the precycle vn is the last vertex of the precycle """ vn = Pk[-1] detectedChordlessCycle = False self.visitedChordlessPaths.add(frozenset(Pk)) Med = self.valuationDomain['med'] if Cycle3: minPreCycleLength = 2 else: # only cycles of length 4 and more are holes in fact minPreCycleLength = 3 #if len(e) > 1: if v0 != vn: # not a reflexive link e = frozenset([v0,vn]) if self.edges[e] > Med and len(Pk) > 2: # we close the chordless pre cycle here detectedChordlessCycle = True ## if Debug: ## print('Pk, len(Pk)', Pk, len(Pk)) if len(Pk) > minPreCycleLength: # only cycles of length 4 and more are holes in fact self.xCC.append(Pk) Pk.append(v0) if Comments: print('Chordless cycle certificate -->>> ', Pk) return detectedChordlessCycle if detectedChordlessCycle == False: NBvn = set(self.gamma[vn]) - set(Pk[1:len(Pk)]) # exterior neighborhood of vn ## if Debug: ## print('vn, NBvn, Pk, Pk[1:len(Pk)] = ', vn, NBvn, Pk, Pk[1:len(Pk)]) #while NBvn != set(): for v in NBvn: # we try in turn all neighbours of vn #v = NBvn.pop() vCP = set(Pk) vCP.add(v) vCP = frozenset(vCP) ## if Debug: ## print('v,vCP =', v,vCP) if vCP not in self.visitedChordlessPaths: # test history of paths P = list(Pk) ## if Debug: ## print('P,P[:-1] = ', P,P[:-1]) noChord = True for x in P[:-1]: ## if Debug: ## print('x = ', x) if x != v0: # we avoid the initial vertex # to stay with a chordless precycle ex = frozenset([x,v]) ## if Debug: ## print('x, v, ex = ',x,v,ex) if self.edges[ex] > Med: # there is a chord noChord = False break if noChord: P.append(v) ## if Debug: ## print('P,v0',P,v0) if self._chordlessPaths(P,v0,Cycle3=Cycle3,Comments=Comments,Debug=Debug): # we continue with the current chordless precycle detectedChordlessCycle=True ## if Debug: ## print('No further chordless precycles from ',vn,' to ',v0) return detectedChordlessCycle def _MISgen(self,S,I): """ generator of maximal independent choices (voir Byskov 2004): * S ::= remaining nodes; * I ::= current independent choice .. note:: - Initialize self.misset = set() before using ! - Inititalize S with a set of vertex keys, and I with an empty set. - See self.showMIS() for usage instructions. """ if S == set(): add = 1 self.missetit = self.misset.copy() for mis in self.missetit: if mis < I: self.misset.remove(mis) else: if I <= mis: add = 0 break if add == 1: self.misset = self.misset | frozenset([I]) yield I else: v = S.pop() Sv = S - (self.gamma[v]) Iv = I | set([v]) for choice in self._MISgen(Sv,Iv): yield choice for choice in self._MISgen(S,I): yield choice def _saveEdges(self,fileName='graphEdges',Agrum=False,Decimal=True): """ Saving graph instances as list of edges, ie node node on each line for enumChordlessCycles C++/agrum progam. """ print('*--- Saving graph edges in file: <' + fileName + '.text> ---*') verticesKeys = [x for x in self.vertices] verticesKeys.sort() Med = self.valuationDomain['med'] fileNameExt = str(fileName)+str('.text') fo = open(fileNameExt, 'w') for i in range( len(verticesKeys) ): for j in range(i+1, len(verticesKeys)): if self.edges[frozenset([verticesKeys[i],verticesKeys[j]])] > Med: if Agrum: fo.write('%d %d\n' % ( i+1, j+1 ) ) else: fo.write('%s %s\n' % ( str(verticesKeys[i]),str(verticesKeys[j]) ) ) fo.close() def _singletons(self): """ List of singletons in frozenset format with neighborhood and not neighbourhood sets. """ s = [] vertices = set([x for x in self.vertices]) for x in vertices: indep = vertices - (self.gamma[x]) s = s + [(frozenset([x]),self.gamma[x],indep)] return s def _computeChordlessCycles(self,Cycle3=False,Comments=True,Debug=False): """ Renders the set of all chordless cycles observed in a Graph intance. Obsolete home brewed version. """ verticesKeys = [x for x in self.vertices] self.visitedChordlessPaths = set() chordlessCycles = [] for v in verticesKeys: P = [v] self.xCC = [] if self._chordlessPaths(P,v,Cycle3=Cycle3,Comments=Comments,Debug=Debug): chordlessCycles += self.xCC self.chordlessCycles = chordlessCycles chordlessCyclesList = [ (x,frozenset(x)) for x in chordlessCycles] if Debug: print('Return list of chordless cycles as a tuple (path, set)') for cc in chordlessCyclesList: print(cc) return chordlessCyclesList
[docs] def computeCliques(self,Comments=False): """ Computes all cliques, ie maximal complete subgraphs in self: .. Note:: - Computes the maximal independent vertex sets in the dual of self. - Result is stored in self.cliques. """ import time,copy if Comments: print('*--- Maximal Cliques ---*') t0 = time.time() dualSelf = -self dualSelf.misset = set() vertices = set([x for x in self.vertices]) self.cliques = [m[0] for m in dualSelf.generateIndependent(dualSelf._singletons()) if m[0] == m[2]] t1 = time.time() if Comments: n = len(vertices) v = [0 for i in range(n+1)] cliqueList = [(len(clique),clique) for clique in self.cliques] cliqueList.sort() for clq in cliqueList: # clq = (len(clique),clique) clique = list(clq[1]) clique.sort() print(clique) v[clq[0]] += 1 print('number of solutions: ', len(self.cliques)) print('cardinality distribution') print('card.: ', list(range(n+1))) print('freq.: ', v) print('execution time: %.5f sec.' % (t1-t0)) print('Results in self.cliques')
[docs] def computeComponents(self): """ Computes the connected components of a graph instance. Returns a partition of the vertices as a list """ components = [] dfs = self.depthFirstSearch() for tree in dfs: components.append(set(tree)) return components
[docs] def computeDegreeDistribution(self,Comments=False): """ Renders the distribution of vertex degrees. """ degreeDistribution = [0 for i in range(self.order)] for v in self.vertices: dv = len(self.gamma[v]) degreeDistribution[dv] += 1 if Comments: print('degrees : ', list(range(self.order))) print('distribution : ', degreeDistribution) return degreeDistribution
[docs] def computeDiameter(self, Oriented = False): """ Renders the diameter (maximal neighbourhood depth) of the digraph instance. .. Note:: The diameter of a disconnected graph is considered to be *infinite* (results in a value -1) ! """ order = self.order nbDepths = self.computeNeighbourhoodDepthDistribution() nbDepths.reverse() if nbDepths[0] != 0: diameter = -1 else: diameter = 0 for i in range(len(nbDepths)): if nbDepths[i+1] != 0: diameter = order - (i+1) break return diameter
[docs] def computeMIS(self,Comments=False): """ Prints all maximal independent vertex sets: .. Note:: - Result is stored in self.misset ! """ import time if Comments: print('*--- Maximal Independent Sets ---*') t0 = time.time() self.misset = set() vertices = set([x for x in self.vertices]) self.misset = [m[0] for m in self.generateIndependent(self._singletons()) if m[0] == m[2]] t1 = time.time() if Comments: n = len(vertices) v = [0 for i in range(n+1)] misList = [(len(mis),mis) for mis in self.misset] misList.sort() for m in misList: # m = (len(mis),mis)) mis = list(m[1]) mis.sort() print(mis) v[m[0]] += 1 print('number of solutions: ', len(misList)) print('cardinality distribution') print('card.: ', list(range(n+1))) print('freq.: ', v) print('execution time: %.5f sec.' % (t1-t0)) print('Results in self.misset')
[docs] def computeNeighbourhoodDepth(self,vertex,Debug=False): """ Renders the distribtion of neighbourhood depths. """ import copy order = self.order vertices = set([x for x in self.vertices]) if Debug: print('-->',vertex) nbx = 0 neighbx = set([vertex]) restVertices = vertices - neighbx while restVertices != set() and nbx < order: if Debug: print('nbx,restVertices', nbx,restVertices) nbx += 1 iterneighbx = copy.copy(neighbx) for y in iterneighbx: neighbx = neighbx | self.gamma[y] if Debug: print('y,self.gamma[y],neighbx', y,self.gamma[y],neighbx) restVertices = vertices - neighbx if Debug: print('nbx,restVertices',nbx,restVertices) if restVertices != set(): return order else: return nbx
[docs] def computeNeighbourhoodDepthDistribution(self,Comments=False,Debug=False): """ Renders the distribtion of neighbourhood depths. """ import copy vertices = set([x for x in self.vertices]) order = self.order vecNeighbourhoodDepth = [0 for i in range(order+1)] for x in vertices: if Debug: print('-->',x) nbx = 0 neighbx = set([x]) restVertices = vertices - neighbx while restVertices != set() and nbx < order: if Debug: print('nbx,restVertices', nbx,restVertices) nbx += 1 iterneighbx = copy.copy(neighbx) for y in iterneighbx: neighbx = neighbx | self.gamma[y] if Debug: print('y,self.gamma[y],neighbx', y,self.gamma[y],neighbx) restVertices = vertices - neighbx if Debug: print('nbx,restVertices',nbx,restVertices) if restVertices != set(): vecNeighbourhoodDepth[order] += 1 else: vecNeighbourhoodDepth[nbx] += 1 if Comments: depths = list(range(self.order)) depths.append('inf.') print('nbh depths : ', depths) print('distribution : ', vecNeighbourhoodDepth) return vecNeighbourhoodDepth
[docs] def computeSize(self): """ Renders the number of positively characterised edges of this graph instance (result is stored in self.size). """ size = 0 Med = self.valuationDomain['med'] for edge in self.edges: if self.edges[edge] > Med: size += 1 self.size = size return size
[docs] def depthFirstSearch(self,Debug=False): """ Depth first search through a graph in lexicographical order of the vertex keys. """ def visitVertex(self, x, Debug = False): """ Visits all followers of vertex x. """ self.vertices[x]['color'] = 1 ## self.date += 1 self.vertices[x]['startDate'] = self.date self.dfsx.append(x) if Debug: print(' dfs %s, date = %d' % (str(self.dfs), self.vertices[x]['startDate'])) nextVertices = [y for y in self.gamma[x]] nextVertices.sort() if Debug: print(' next ', nextVertices) for y in nextVertices: if self.vertices[y]['color'] == 0: self.date += 1 visitVertex(self,y, Debug = Debug) if self.vertices[x]['color'] == 1: self.dfsx.append(x) self.vertices[x]['color'] = 2 self.vertices[x]['endDate'] = self.date self.date += 1 def visitAllVertices(self, Debug=False): """ Mark the starting date for all vertices and controls the progress of the search with vertices colors: White (0), Grey (1), Black (2) """ self.dfs = [] for x in self.vertices: self.vertices[x]['color'] = 0 self.date = 0 verticesList = [x for x in self.vertices] verticesList.sort() for x in verticesList: self.dfsx = [] if self.vertices[x]['color'] == 0: if Debug: print('==>> Starting from %s ' % x) visitVertex(self, x, Debug = Debug) self.dfs.append(self.dfsx) #self.vertices[x]['color'] = 2 #self.vertices[x]['endDate'] = self.date # ---- main ----- visitAllVertices(self, Debug=Debug) return self.dfs
[docs] def exportGraphViz(self,fileName=None,verticesSubset=None, noSilent=True, graphType='png',graphSize='7,7', withSpanningTree=False, layout=None, arcColor='black', lineWidth=1): """ Exports GraphViz dot file for graph drawing filtering. Example: >>> g = Graph(numberOfVertices=5,edgeProbability=0.3) >>> g.exportGraphViz('randomGraph')) .. image:: randomGraph.png :alt: Random graph :width: 300 px :align: center """ import os if noSilent: print('*---- exporting a dot file for GraphViz tools ---------*') if verticesSubset == None: vertexkeys = [x for x in self.vertices] else: vertexkeys = [x for x in verticesSubset] n = len(vertexkeys) edges = self.edges Med = self.valuationDomain['med'] i = 0 if fileName == None: name = self.name else: name = fileName dotName = name+'.dot' if noSilent: print('Exporting to '+dotName) ## if bestChoice != set(): ## rankBestString = '{rank=max; ' ## if worstChoice != set(): ## rankWorstString = '{rank=min; ' fo = open(dotName,'w') fo.write('strict graph G {\n') fo.write('graph [ bgcolor = cornsilk, fontname = "Helvetica-Oblique",\n fontsize = 12,\n label = "') fo.write('\\nGraphs Python module (graphviz), R. Bisdorff, 2015", size="') fo.write(graphSize),fo.write('"];\n') for i in range(n): try: nodeName = str(self.vertices[vertexkeys[i]]['shortName']) except: try: nodeName = self.vertices[vertexkeys[i]]['name'] except: nodeName = str(vertexkeys[i]) node = 'n'+str(i+1)+' [shape = "circle", label = "' +nodeName+'"' try: if self.vertices[vertexkeys[i]]['spin'] == 1: node += ', style = "filled", color = %s' % spinColor except: pass node += '];\n' fo.write(node) if withSpanningTree: try: dfs = self.dfs except: print('no spanning tree yet computed. Run self.randomDepthFirstSearch() !') edgesColored = set() print(dfs) for tree in dfs: for i in range((len(tree)-1)): #print(i,tree[i],tree[i+1]) edgesColored.add(frozenset([tree[i],tree[i+1]])) #print('Spanning tree: ', edgesColored) for i in range(n): for j in range(i+1, n): if i != j: edge = 'n'+str(i+1) if edges[frozenset( [vertexkeys[i], vertexkeys[j]])] > Med: if withSpanningTree and \ frozenset( [vertexkeys[i], vertexkeys[j]]) in edgesColored: arrowFormat = \ ' [dir=both,style="setlinewidth(3)",color=red, arrowhead=none, arrowtail=none] ;\n' else: arrowFormat = \ ' [dir=both,style="setlinewidth(%d)",color=%s, arrowhead=none, arrowtail=none] ;\n' % (lineWidth,arcColor) edge0 = edge+'-- n'+str(j+1)+arrowFormat fo.write(edge0) elif edges[frozenset([vertexkeys[i],vertexkeys[j]])] == Med: edge0 = edge+'-- n'+str(j+1)+\ ' [dir=both, color=grey, arrowhead=none, arrowtail=none] ;\n' fo.write(edge0) fo.write('}\n') fo.close() # choose layout model if isinstance(self,(GridGraph,TriangulatedGrid,RandomTree)): if layout == None: layout = 'neato' elif isinstance(self,(CycleGraph)): if layout == None: layout = 'circo' else: if layout == None: layout = 'fdp' commandString = layout+' -T'+graphType+' '+dotName+' -o '+name+'.'+graphType if noSilent: print(commandString) try: os.system(commandString) except: if noSilent: print('graphViz tools not avalaible! Please check installation.') print('On Ubuntu: ..$ sudo apt-get install graphviz')
[docs] def gammaSets(self,Debug=False): """ renders the gamma function as dictionary """ vkeys = [x for x in self.vertices] if Debug: print('vkeys', vkeys) edges = self.edges gamma = dict() for v in vkeys: gamma[v] = set() for e in edges: if edges[e] > 0: ## if Debug: ## print('e', e) pair = set(e) e1 = pair.pop() e2 = pair.pop() gamma[e1].add(e2) gamma[e2].add(e1) return gamma
[docs] def generateIndependent(self,U): """ Generator for all independent vertices sets with neighborhoods of a graph instance: .. note:: * Initiate with U = self._singletons(). * Yields [independent set, covered set, all vertices - covered set)]. * If independent set == (all vertices - covered set), the given independent set is maximal ! """ if U == []: vertices = set([x for x in self.vertices]) yield [frozenset(),set(),vertices] else: x = list(U.pop()) for S in self.generateIndependent(U): yield S if x[0] <= S[2]: Sxgam = S[1] | x[1] Sxindep = S[2] & x[2] Sxchoice = S[0] | x[0] Sx = [Sxchoice,Sxgam,Sxindep] yield Sx
[docs] def graph2Digraph(self): """ Converts a Graph object into a symmetric Digraph object. """ from copy import deepcopy from digraphs import EmptyDigraph dg = EmptyDigraph(order=self.order) dg.name = deepcopy(self.name) dg.actions = deepcopy(self.vertices) dg.order = len(dg.actions) dg.valuationdomain = deepcopy(self.valuationDomain) dg.convertValuationToDecimal() dg.relations = {} actionsKeys = [x for x in dg.actions] for x in actionsKeys: dg.relation[x] = {} for y in actionsKeys: if x == y: dg.relation[x][y] = dg.valuationdomain['med'] else: edgeKey = frozenset([x,y]) dg.relation[x][y] = self.edges[edgeKey] dg.convertRelationToDecimal() dg.gamma = dg.gammaSets() dg.notGamma = dg.notGammaSets() return dg
[docs] def isConnected(self): """ Cheks if self is a connected graph instance. """ dfs = self.depthFirstSearch() if len(dfs) == 1: return True else: return False
[docs] def isTree(self): """ Checks if self is a tree by verifing the required number of edges: order-1; and the existence of leaves. """ n = self.order m = self.size if n == 1: return True elif m != n-1: return False else: nbrOfLeaves = 0 for x in self.vertices: degreex = len(self.gamma[x]) if degreex == 0: # isolated vertex return False elif degreex == 1: nbrOfLeaves += 1 if nbrOfLeaves < 2: # a cycle graph return False else: return True
[docs] def randomDepthFirstSearch(self,seed=None,Debug=False): """ Depth first search through a graph in random order of the vertex keys. .. Note:: The resulting spanning tree (or forest) is by far not uniformly selected among all possible trees. Spanning stars will indeed be much less probably selected then streight walks ! """ import random random.seed(seed) def visitVertex(self,x,Debug=False): """ Visits all followers of vertex x. """ self.vertices[x]['color'] = 1 ## self.date += 1 self.vertices[x]['startDate'] = self.date self.dfsx.append(x) if Debug: print(' dfs %s, date = %d' % (str(self.dfs), self.vertices[x]['startDate'])) nextVertices = [y for y in self.gamma[x]] nextVertices.sort() if Debug: print(' next ', nextVertices) while nextVertices != []: y = random.choice(nextVertices) if self.vertices[y]['color'] == 0: self.date += 1 visitVertex(self,y,Debug=Debug) if self.vertices[x]['color'] == 1: self.dfsx.append(x) nextVertices.remove(y) self.vertices[x]['color'] = 2 self.vertices[x]['endDate'] = self.date self.date += 1 def visitAllVertices(self,Debug=False): """ Mark the starting date for all vertices and controls the progress of the search with vertices colors: White (0), Grey (1), Black (2) """ self.dfs = [] for x in self.vertices: self.vertices[x]['color'] = 0 self.date = 0 verticesList = [x for x in self.vertices] verticesList.sort() while verticesList != []: x = random.choice(verticesList) self.dfsx = [] if self.vertices[x]['color'] == 0: if Debug: print('==>> Starting from %s ' % x) visitVertex(self,x,Debug=Debug) self.dfs.append(self.dfsx) verticesList.remove(x) # ---- main ----- visitAllVertices(self,Debug=Debug) return self.dfs
[docs] def recodeValuation(self,newMin=-1,newMax=1,Debug=False): """ Recodes the characteristic valuation domain according to the parameters given. .. note:: Default values gives a normalized valuation domain """ from copy import deepcopy oldMax = self.valuationDomain['max'] oldMin = self.valuationDomain['min'] oldMed = self.valuationDomain['med'] oldAmplitude = oldMax - oldMin if Debug: print(oldMin, oldMed, oldMax, oldAmplitude) newMin = Decimal(str(newMin)) newMax = Decimal(str(newMax)) newMed = Decimal('%.3f' % ((newMax + newMin)/Decimal('2.0'))) newAmplitude = newMax - newMin if Debug: print(newMin, newMed, newMax, newAmplitude) verticesList = [x for x in self.vertices] oldEdges = self.edges newEdges = {} for i in range(self.order): x = verticesList[i] for j in range(i+1,self.order): y = verticesList[j] edge = frozenset([x,y]) if oldEdges[edge] == oldMax: newEdges[edge] = newMax elif oldEdges[edge] == oldMin: newEdges[edge] = newMin elif oldEdges[edge] == oldMed: newEdges[edge] = newMed else: newEdges[edge] = newMin + ((oldEdges[edge] - oldMin)/oldAmplitude)*newAmplitude if Debug: print(edge,oldEdges[edge],newEdges[edge]) # install new values in self self.valuationDomain['max'] = newMax self.valuationDomain['min'] = newMin self.valuationDomain['med'] = newMed self.valuationDomain['hasIntegerValuation'] = False self.edges = deepcopy(newEdges)
[docs] def save(self,fileName='tempGraph',Debug=False): """ Persistent storage of a Graph class instance in the form of a python source code file. """ print('*--- Saving graph in file: <' + fileName + '.py> ---*') verticesKeys = [x for x in self.vertices] verticesKeys.sort() #order = len(self.vertices) edges = self.edges Min = self.valuationDomain['min'] Med = self.valuationDomain['med'] Max = self.valuationDomain['max'] fileNameExt = str(fileName)+str('.py') fo = open(fileNameExt, 'w') fo.write('# Graph instance saved in Python format\n') fo.write('from decimal import Decimal\n') fo.write('vertices = {\n') for x in verticesKeys: fo.write('\'%s\': %s,\n' % (x,self.vertices[x])) fo.write('}\n') fo.write('valuationDomain = {\'min\':'+ str(Min)+',\'med\':'+str(Med)+',\'max\':'+str(Max)+'}\n') fo.write('edges = {\n') for i in range(self.order): for j in range(i+1,self.order): fo.write('frozenset([\'%s\',\'%s\']) : %d, \n' % (verticesKeys[i],verticesKeys[j],edges[frozenset([verticesKeys[i],verticesKeys[j]])])) fo.write( '}\n') fo.close()
[docs] def setEdgeValue(self,edge,value,Comments=False): """ Wrapper for updating the charactreistic valuation of a Graph instance. The egde parameter consists in a pair of vertices; edge = ('v1','v2') for instance. The new value must be in the limits of the valuation domain. """ from decimal import Decimal # check vertices' existance verticesIds = [x for x in self.vertices] if (edge[0] not in verticesIds) or (edge[1] not in verticesIds): #self.showShort() print('!!! Error: edge %s not found !!!' % str(edge)) return # check new edge value Min = self.valuationDomain['min'] Max = self.valuationDomain['max'] newValue = Decimal(str(value)) if newValue > Max or value < Min: print('!!! Error: edge value %s out of range !!!' % str(value)) print(self.valuationDomain) return # set new value self.edges[frozenset(edge)] = Decimal(str(value)) self.gamma = self.gammaSets() if Comments: print('edge %s put to value %s.' % (edge,str(value)))
[docs] def showCliques(self): self.computeCliques(Comments=True)
[docs] def showMIS(self): self.computeMIS(Comments=True)
[docs] def showMore(self): """ Generic show method for Graph instances. """ print('*---- Properties of the graph ----*') print('Name : \'%s\'' % (self.name) ) vKeys = [x for x in self.vertices] vKeys.sort() print('Vertices : ', vKeys) print('Order : ', self.order) self.showMIS() self.showCliques()
[docs] def showShort(self): """ Generic show method for Graph instances. """ print('*---- short description of the graph ----*') print('Name : \'%s\'' % (self.name) ) vKeys = [x for x in self.vertices] vKeys.sort() print('Vertices : ', vKeys) print('Valuation domain : ', self.valuationDomain) print('Gamma function : ') for v in vKeys: print('%s -> %s' % (v, list(self.gamma[v]))) self.computeDegreeDistribution(Comments=True) self.computeNeighbourhoodDepthDistribution(Comments=True)
#----------------
[docs]class EmptyGraph(Graph): """ Intantiates graph of given order without any positively valued edge. *Parameter*: * order (positive integer) """ def __init__(self,order=5,seed=None): self.name = 'emptyGraph' self.order = order nd = len(str(order)) vertices = dict() for i in range(order): vertexKey = ('v%%0%dd' % nd) % (i+1) vertices[vertexKey] = {'shortName':vertexKey, 'name': 'random vertex'} self.vertices = vertices self.valuationDomain = {'min':Decimal('-1'),'med':Decimal('0'),'max':Decimal('1')} Min = self.valuationDomain['min'] edges = dict() verticesList = [v for v in vertices] verticesList.sort() for x in verticesList: for y in verticesList: if x != y: edgeKey = frozenset([x,y]) edges[edgeKey] = Min self.edges = edges self.size = self.computeSize() self.gamma = self.gammaSets()
[docs]class CompleteGraph(Graph): """ Instances of complete graphs bipolarly valuated in {-1,0,+1}. Each vertex x is positively linked to all the other vertices (edges[{x,y}] = +1) *Parameter*: * order (positive integer) """ def __init__(self,order=5,seed=None): self.name = 'completeGraph' self.order = order nd = len(str(order)) vertices = dict() for i in range(order): vertexKey = ('v%%0%dd' % nd) % (i+1) vertices[vertexKey] = {'shortName':vertexKey, 'name': 'random vertex'} self.vertices = vertices self.valuationDomain = {'min':Decimal('-1'),'med':Decimal('0'),'max':Decimal('1')} Max = self.valuationDomain['max'] edges = dict() verticesList = [v for v in vertices] verticesList.sort() for x in verticesList: for y in verticesList: if x != y: edgeKey = frozenset([x,y]) edges[edgeKey] = Max self.edges = edges self.size = self.computeSize() self.gamma = self.gammaSets()
[docs]class DualGraph(Graph): """ Instantiates the dual Graph object of a given other Graph instance. The relation constructor returns the dual of self.relation with formula: relationOut[a][b] = Max - self.relation[a][b] + Min where Max (resp. Min) equals valuation maximum (resp. minimum). """ def __init__(self,other): from copy import deepcopy self.name = 'dual_' + str(other.name) try: self.description = deepcopy(other.description) except AttributeError: pass self.valuationDomain = deepcopy(other.valuationDomain) Max = self.valuationDomain['max'] Min = self.valuationDomain['min'] self.vertices = deepcopy(other.vertices) self.order = len(self.vertices) self.edges = {} for e in other.edges: self.edges[e] = Max - other.edges[e] + Min self.__class__ = other.__class__ self.gamma = self.gammaSets()
[docs]class CycleGraph(Graph): """ Instances of cycle graph characterized in [-1,1]. *Parameter*: * order (positive integer) Example of 7-cycle graph instance: .. image:: 7cycle.png :alt: 7-cycle instance :width: 300 px :align: center """ def __init__(self,order=5,seed=None,Debug=False): from collections import OrderedDict self.name = 'cycleGraph' self.order = order nd = len(str(order)) vertices = OrderedDict() for i in range(order): vertexKey = ('v%%0%dd' % nd) % (i+1) vertices[vertexKey] = {'shortName':vertexKey, 'name': 'random vertex'} self.vertices = vertices verticesList = [key for key in vertices] self.valuationDomain = {'min':Decimal('-1'),'med':Decimal('0'),'max':Decimal('1')} Min = self.valuationDomain['min'] Max = self.valuationDomain['max'] edges = OrderedDict() #verticesList = [v for v in vertices] #verticesList.sort() for x in vertices: for y in vertices: if x != y: edgeKey = frozenset([x,y]) edges[edgeKey] = Min for i in range(order-1): edgeKey = frozenset(verticesList[i:i+2]) edges[edgeKey] = Max if Debug: print(edgeKey) x = verticesList[-1] y = verticesList[0] edgeKey = frozenset([x,y]) edges[edgeKey] = Max if Debug: print(edgeKey) self.edges = edges self.size = self.computeSize() self.gamma = self.gammaSets()
[docs]class RandomGraph(Graph): """ Random instances of the Graph class *Parameters*: * order (positive integer) * edgeProbability (in [0,1]) """ def __init__(self,order=5,edgeProbability=0.4,seed=None): import random random.seed(seed) self.name = 'randomGraph' self.order = order nd = len(str(order)) vertices = dict() for i in range(order): vertexKey = ('v%%0%dd' % nd) % (i+1) vertices[vertexKey] = {'shortName':vertexKey, 'name': 'random vertex'} self.vertices = vertices self.valuationDomain = {'min':Decimal('-1'),'med':Decimal('0'),'max':Decimal('1')} Min = self.valuationDomain['min'] Max = self.valuationDomain['max'] edges = dict() verticesList = [v for v in vertices] verticesList.sort() for i in range(order): x = verticesList[i] for j in range(i+1,order): y = verticesList[j] edgeKey = frozenset([x,y]) if random.random() > 1.0 - edgeProbability: edges[edgeKey] = Max else: edges[edgeKey] = Min self.edges = edges self.size = self.computeSize() self.gamma = self.gammaSets()
[docs]class RandomValuationGraph(Graph): """ Specialization of the genuine Graph class for generating temporary randomly valuated graphs in the range [-1.0;1.0]. *Parameter*: * order (positive integer) * ndigits (decimal precision) """ def __init__(self,order=5,ndigits=2,seed=None): import random random.seed(seed) self.name = 'randomGraph' self.order = order nd = len(str(order)) vertices = dict() for i in range(order): vertexKey = ('v%%0%dd' % nd) % (i+1) vertices[vertexKey] = {'shortName':vertexKey, 'name': 'random vertex'} self.vertices = vertices self.valuationDomain = {'min':Decimal('-1'),'med':Decimal('0'),'max':Decimal('1')} Min = float(self.valuationDomain['min']) Max = float(self.valuationDomain['max']) edges = dict() verticesList = [v for v in vertices] verticesList.sort() for i in range(order): x = verticesList[i] for j in range(i+1,order): y = verticesList[j] edgeKey = frozenset([x,y]) weightString = ('%%.%df' % ndigits) % random.uniform(Min,Max) edges[edgeKey] = Decimal(weightString) self.edges = edges self.size = self.computeSize() self.gamma = self.gammaSets()
[docs]class RandomRegularGraph(Graph): """ Specialization of the general Graph class for generating temporary random regular graphs of fixed degrees. """ def __init__(self,order=7,degree=2,seed=None): from randomDigraphs import RandomRegularDigraph rdg = RandomRegularDigraph(order=order, degree=degree, seed=seed) rg = rdg.digraph2Graph() self.vertices = rg.vertices self.valuationDomain = rg.valuationDomain self.edges = rg.edges self.name = 'randomRegularGraph' self.order = len(self.vertices) self.size = self.computeSize() self.gamma = self.gammaSets()
[docs]class RandomFixedSizeGraph(Graph): """ Generates a random graph with a fixed size (number of edges), by instantiating a fixed numbers of arcs from random choices in the set of potential pairs of vertices numbered from 1 to order. """ def __init__(self,order=7,size=14,seed=None,Debug=False): import random random.seed(seed) # check feasability r = ((order * order) - order)//2 if size > r : print('Graph not feasable: size exceeds number of potential edges = %d !!' % r) return print(order,size,r) self.name = 'randomFixedSize' self.order = order nd = len(str(order)) vertices = dict() for i in range(order): vertexKey = ('v%%0%dd' % nd) % (i+1) vertices[vertexKey] = {'shortName':vertexKey, 'name': 'random vertex'} self.vertices = vertices if Debug: print(self.vertices) self.valuationDomain = {'min':Decimal('-1'),'med':Decimal('0'),'max':Decimal('1')} edges = dict() Min = self.valuationDomain['min'] verticesList = [v for v in vertices] verticesList.sort() for x in verticesList: for y in verticesList: if x != y: edgeKey = frozenset([x,y]) edges[edgeKey] = Min if Debug: print(edges) edgesKeys = [key for key in edges] edgesKeys.sort() if Debug: print(edgesKeys) Max = self.valuationDomain['max'] for i in range(size): edgeKey = random.choice(edgesKeys) edges[edgeKey] = Max edgesKeys.remove(edgeKey) if Debug: print(i,edgeKey,edgesKeys) self.edges = edges self.size = size self.gamma = self.gammaSets()
[docs]class RandomFixedDegreeSequenceGraph(Graph): """ Specialization of the general Graph class for generating temporary random graphs with a fixed sequence of degrees. .. warning:: The implementation is not guaranteeing a uniform choice among all potential valid graph instances. """ def __init__(self,order=7,degreeSequence=[3,3,2,2,1,1,0],seed=None): from randomDigraphs import RandomFixedDegreeSequenceDigraph rdg = RandomFixedDegreeSequenceDigraph(order=order, degreeSequence=degreeSequence, seed=seed) rg = rdg.digraph2Graph() self.vertices = rg.vertices self.order = order self.valuationDomain = rg.valuationDomain self.edges = rg.edges self.size = self.computeSize() self.name = 'randomFixedDegreeSequenceGraph' self.gamma = self.gammaSets()
[docs]class GridGraph(Graph): """ Specialization of the general Graph class for generating temporary Grid graphs of dimension n times m. *Parameters*: * n,m > 0 * valuationDomain ={'min':-1, 'med':0, 'max':+1} Default instantiation (5 times 5 Grid Digraph): * n = 5, * m=5, * valuationDomain = {'min':-1.0,'max':1.0}. Example of 5x5 GridGraph instance: .. image:: grid-5-5.png :alt: 5x5 grid instance :width: 300 px :align: center """ def __init__(self,n=5,m=5,valuationMin=-1,valuationMax=1): self.name = 'grid-'+str(n)+'-'+str(m) self.n = n self.m = m na = list(range(n+1)) na.remove(0) ma = list(range(m+1)) ma.remove(0) vertices = {} gridNodes={} for x in na: for y in ma: vertex = str(x)+'-'+str(y) gridNodes[vertex]=(x,y) vertices[vertex] = {'name': 'gridnode', 'shortName': vertex} order = len(vertices) self.order = order self.vertices = vertices self.gridNodes = gridNodes Min = Decimal(str(valuationMin)) Max = Decimal(str(valuationMax)) Med = Decimal(str((Max + Min)/Decimal('2'))) self.valuationDomain = {'min':Min,'med':Med,'max':Max} edges = {} # instantiate edges verticesKeys = [x for x in vertices] for x in verticesKeys: for y in verticesKeys: if x != y: if gridNodes[x][1] == gridNodes[y][1]: if gridNodes[x][0] == gridNodes[y][0]-1 : edges[frozenset([x,y])] = Max elif gridNodes[x][0] == gridNodes[y][0]+1: edges[frozenset([x,y])] = Max else: edges[frozenset([x,y])] = Min elif gridNodes[x][0] == gridNodes[y][0]: if gridNodes[x][1] == gridNodes[y][1]-1: edges[frozenset([x,y])] = Max elif gridNodes[x][1] == gridNodes[y][1]+1: edges[frozenset([x,y])] = Max else: edges[frozenset([x,y])] = Min else: edges[frozenset([x,y])] = Min self.edges = edges self.size = self.computeSize() self.gamma = self.gammaSets()
[docs] def showShort(self): print('*----- show short --------------*') print('Grid graph : ', self.name) print('n : ', self.n) print('m : ', self.m) print('order : ', self.order) print('size : ', self.size)
#--------------------------
[docs]class SnakeGraph(GridGraph): """ Snake graphs S(p/q) are made up of all the integer grid squares between the lower and upper Christofel paths of the rational number p/q, where p and q are two coprime integers such that 0 <= p <= q, i.e. p/q gives an irreducible ratio between 0 and 1. *Reference*: M. Aigner, Markov's Theorem and 100 Years of the Uniqueness Conjecture, Springer, 2013, p. 141-149 S(4/7) snake graph instance:: >>> from graphs import SnakeGraph >>> s4_7 = SnakeGraph(p=4,q=7) >>> s4_7.showShort() *---- short description of the snake graph ----* Name : 'snakeGraph' Rational p/q : 4/7 Christoffel words: Upper word : BBABABA Lower word : ABABABB >>> s4_7.exportGraphViz('4_7_snake',lineWidth=3,arcColor='red') .. image:: 4_7_snake.png :alt: 4/7 snake graph instance :width: 300 px :align: center """ def __init__(self,p,q): from math import floor, ceil self.name = '%d_%d_snakeGraph' % (p,q) # vertices self.n = q self.m = p vertices = {} gridNodes={} for x in range(q+1): for y in range(p+1): vertex = str(x)+'-'+str(y) gridNodes[vertex]=(x,y) vertices[vertex] = {'name': 'gridnode', 'shortName': vertex} order = len(vertices) self.order = order self.vertices = vertices self.gridNodes = gridNodes # valuation domain Min = Decimal('-1') Med = Decimal('0') Max = Decimal('1') self.valuationDomain = {'min': Min, 'med': Med, 'max': Max} # edges edges = {} # instantiate edges verticesKeys = [x for x in vertices] for x in verticesKeys: for y in verticesKeys: if x != y: if gridNodes[x][1] == gridNodes[y][1]: if gridNodes[x][0] == gridNodes[y][0]-1 : edges[frozenset([x,y])] = Med elif gridNodes[x][0] == gridNodes[y][0]+1: edges[frozenset([x,y])] = Med else: edges[frozenset([x,y])] = Min elif gridNodes[x][0] == gridNodes[y][0]: if gridNodes[x][1] == gridNodes[y][1]-1: edges[frozenset([x,y])] = Med elif gridNodes[x][1] == gridNodes[y][1]+1: edges[frozenset([x,y])] = Med else: edges[frozenset([x,y])] = Min else: edges[frozenset([x,y])] = Min self.edges = edges # snake lower Christoffel path k = [0 for i in range(q+1)] for i in range(q+1): k[i] = floor((p/q)*i) print(k) for i in range(q): if k[i] == k[i+1]: x = '%d-%d' % (i,k[i]) y = '%d-%d' % (i+1,k[i+1]) self.setEdgeValue((x,y),Max) else: x = '%d-%d' % (i,k[i]) y = '%d-%d' % (i+1,k[i]) self.setEdgeValue((x,y),Max) x = '%d-%d' % (i+1,k[i]) y = '%d-%d' % (i+1,k[i+1]) self.setEdgeValue((x,y),Max) # snake upper Christoffel path K = [0 for i in range(q+1)] for i in range(q+1): K[i] = ceil((p/q)*i) print(K) for i in range(q): if K[i] == K[i+1]: x = '%d-%d' % (i,K[i]) y = '%d-%d' % (i+1,K[i]) print(x,y) self.setEdgeValue((x,y),Max) else: x = '%d-%d' % (i,K[i]) y = '%d-%d' % (i,K[i+1]) print(x,y) self.setEdgeValue((x,y),Max) x = '%d-%d' % (i,K[i+1]) y = '%d-%d' % (i+1,K[i+1]) print(x,y) self.setEdgeValue((x,y),Max) # storing graph instance self.size = self.computeSize() self.gamma = self.gammaSets() # ChristoffelWord lcw = '' ucw = '' for i in range(1,q+1): if k[i]-k[i-1] == 0: lcw += 'A' else: lcw += 'B' if K[i]-K[i-1] == 0: ucw += 'A' else: ucw += 'B' self.cw = (lcw,ucw)
[docs] def showShort(self,WithVertices=False): """ Show method for SnakeGraph instances. """ print('*---- short description of the snake graph ----*') print('Name : \'%s\'' % (self.name) ) print('Rational p/q : %d/%d' % (self.m,self.n)) print('Christoffel words:') print('Upper word : ',self.cw[1]) print('Lower word : ',self.cw[0]) if WithVertices: vKeys = [x for x in self.vertices] vKeys.sort() print('Vertices : ', vKeys) print('Valuation domain : ', self.valuationDomain) print('Gamma function : ') for v in vKeys: if self.gamma[v] != set(): print('%s -> %s' % (v, list(self.gamma[v])))
def __repr__(self): """ Show method for SnakeGraph instances. """ print('*---- short description of the snake graph ----*') print('Name : \'%s\'' % (self.name) ) print('Rational p/q : %d/%d' % (self.m,self.n)) print('Christoffel words:') print('Upper word : ',self.cw[1]) print('Lower word : ',self.cw[0]) return 'graphs.SnakeGraph(GridGraph) instance'
#---------------------
[docs]class TriangulatedGrid(Graph): """ Specialization of the general Graph class for generating temporary triangulated grids of dimension n times m. *Parameters*: * n,m > 0 * valuationDomain = {'min':m, 'max':M} Example of 5x5 triangulated grid instance: .. image:: triangular-5-5.png :alt: triangulated 5x5 grid :width: 300 px :align: center """ def __init__(self,n=5,m=5,valuationMin=-1,valuationMax=1): self.name = 'triangulated-'+str(n)+'-'+str(m) self.n = n self.m = m na = list(range(n+1)) na.remove(0) ma = list(range(m+1)) ma.remove(0) vertices = {} gridNodes={} for x in na: for y in ma: vertex = str(x)+'-'+str(y) gridNodes[vertex]=(x,y) vertices[vertex] = {'name': 'gridnode', 'shortName': vertex} order = len(vertices) self.order = order self.vertices = vertices self.gridNodes = gridNodes Min = Decimal(str(valuationMin)) Max = Decimal(str(valuationMax)) Med = Decimal((Max + Min)/Decimal('2')) self.valuationDomain = {'min':Min,'med':Med,'max':Max} edges = {} # instantiate edges verticesKeys = [x for x in vertices] for x in verticesKeys: for y in verticesKeys: if x != y: if gridNodes[x][1] == gridNodes[y][1]: if gridNodes[x][0] == gridNodes[y][0]-1 : edges[frozenset([x,y])] = Max elif gridNodes[x][0] == gridNodes[y][0]+1: edges[frozenset([x,y])] = Max else: edges[frozenset([x,y])] = Min elif gridNodes[x][0] == gridNodes[y][0]: if gridNodes[x][1] == gridNodes[y][1]-1: edges[frozenset([x,y])] = Max elif gridNodes[x][1] == gridNodes[y][1]+1: edges[frozenset([x,y])] = Max else: edges[frozenset([x,y])] = Min elif gridNodes[x][0] == gridNodes[y][0]-1: if gridNodes[x][1] == gridNodes[y][1]-1: edges[frozenset([x,y])] = Max else: edges[frozenset([x,y])] = Min elif gridNodes[x][0] == gridNodes[y][0]+1: if gridNodes[x][1] == gridNodes[y][1]+1: edges[frozenset([x,y])] = Max else: edges[frozenset([x,y])] = Min else: edges[frozenset([x,y])] = Min self.edges = edges self.size = self.computeSize() self.gamma = self.gammaSets()
[docs] def showShort(self): print('*----- show summary --------------*') print('Triangular graph : ', self.name) print('n x m : ', self.n, 'x', self.m) print('order : ', self.order) print('size : ', self.size)
[docs]class RandomTree(Graph): """ Instance of a tree generated from a random (or a given) Prüfer code. .. image:: randomTree.png :alt: radomTree instance :width: 300 px :align: center """ def __init__(self,order=None, vertices= None, prueferCode = None, seed = None, Debug=False): import random random.seed(seed) self.name='randomTree' if order == None: if prueferCode == None: order = 6 else: order = len(prueferCode) + 2 self.order = order if Debug: print(self.name, self.order) if vertices == None: vertices = dict() for i in range(order): vertexKey = 'v%d' % (i+1) vertices[vertexKey] = {'shortName':str(vertexKey), 'name': 'random vertex'} self.vertices = vertices if Debug: print('vertices = ', self.vertices) self.valuationDomain = {'min':Decimal('-1'), 'med':Decimal('0'), 'max':Decimal('1')} Min = self.valuationDomain['min'] Med = self.valuationDomain['med'] Max = self.valuationDomain['max'] if Debug: print('valuationDomain = ', self.valuationDomain) verticesList = [x for x in self.vertices.keys()] verticesList.sort() if Debug: print(verticesList) edges = dict() for i in range(order): for j in range(i+1,order): edgeKey = frozenset([verticesList[i],verticesList[j]]) edges[edgeKey] = Min self.edges = edges if Debug: print('edges = ',self.edges) if prueferCode == None: prueferCode = [] for k in range(order-2): prueferCode.append( verticesList[random.choice( list(range(order)) )] ) self.prueferCode = prueferCode if Debug: print('prueferCode = ', self.prueferCode) degree = {} for x in verticesList: degree[x] = 1 for i in range(order-2): degree[prueferCode[i]] += 1 if Debug: print('degrees = ', degree) for i in range(order-2): for j in range(order): if degree[verticesList[j]] == 1: edgeKey = frozenset([prueferCode[i],verticesList[j]]) self.edges[edgeKey] = Max degree[verticesList[j]] -= 1 degree[prueferCode[i]] -= 1 break lastEdgeKey = frozenset([verticesList[i] for i in range(order) if degree[verticesList[i]] > 0]) self.edges[lastEdgeKey] = Max if Debug: print('updated edges = ', self.edges) self.size = self.computeSize() self.gamma = self.gammaSets(Debug) if Debug: print('gamma = ', self.gamma)
[docs] def tree2Pruefer(self,vertices=None,Debug=False): """ Renders the Pruefer code of a given tree. """ if vertices == None: vertices = self.vertices verticesList = [x for x in vertices] verticesList.sort() if Debug: print('verticesList = ',verticesList) np = len(verticesList)-2 Med = self.valuationDomain['med'] edges = [e for e in self.edges if self.edges[e] > Med] degree = dict() for v in verticesList: degree[v] = len(self.gamma[v]) if Debug: print('degrees = ', degree) prueferCode = [] for i in range(np): leaves = [v for v in verticesList if degree[v] == 1] leaves.sort() leavesEdges = [e for e in edges if leaves[0] in e] for v in leavesEdges[0]: if degree[v] > 1: prueferCode.append(v) degree[v] -= 1 edges.remove(leavesEdges[0]) return prueferCode
[docs]class RandomSpanningForest(RandomTree): """ Random instance of a spanning forest (one or more trees) generated from a random depth first search graph g traversal. .. image:: spanningForest.png :alt: randomSpanningForest instance :width: 300 px :align: center """ def __init__(self,g,seed=None,Debug=False): from copy import deepcopy import random random.seed(seed) self.name= g.name+'_randomSpanningTree' if Debug: print(self.name) self.vertices = deepcopy(g.vertices) order = len(self.vertices) self.order = order self.valuationDomain = deepcopy(g.valuationDomain) Min = self.valuationDomain['min'] Med = self.valuationDomain['med'] Max = self.valuationDomain['max'] if Debug: print('valuationDomain = ', self.valuationDomain) verticesList = [x for x in self.vertices] verticesList.sort() edges = dict() for i in range(order): for j in range(i+1,order): edgeKey = frozenset([verticesList[i],verticesList[j]]) if g.edges[edgeKey] > Med: edges[edgeKey] = Med else: edges[edgeKey] = g.edges[edgeKey] self.edges = deepcopy(edges) if Debug: print('edges = ',self.edges) self.dfs = g.randomDepthFirstSearch(seed=seed,Debug=Debug) if Debug: print('dfs = ', self.dfs) components = [] for tree in self.dfs: component = set() n = len(tree) for i in range(n-1): component.add(tree[i]) component.add(tree[i+1]) edgeKey = frozenset([tree[i],tree[i+1]]) self.edges[edgeKey] = Max if Debug: print('tree = ',tree) print('component = ',component) components.append(component) if Debug: print('updated edges = ', self.edges) self.size = self.computeSize() self.gamma = self.gammaSets(Debug) if Debug: print('gamma = ', self.gamma) prueferCodes = [] for tree in self.dfs: component = set(tree) if len(component) > 2: prueferCodes.append({'component': component, 'code': self.tree2Pruefer(vertices=component) }) else: prueferCodes.append({'component': component, 'code': []}) self.prueferCodes = prueferCodes
[docs]class RandomSpanningTree(RandomTree): """ Uniform random instance of a spanning tree generated with Wilson's algorithm from a connected Graph g instance. .. Note:: Wilson's algorithm only works for connecte graphs. .. image:: randomSpanningTree.png :alt: randomSpanningTree instance :width: 300 px :align: center """ def __init__(self,g,seed=None,Debug=False): from copy import copy as copy if not g.isConnected(): print('Error !: Wilson\'s algorithm requires a connected graph!') return import random random.seed(seed) self.name= g.name+'_randomSpanningTree' if Debug: print(self.name) self.vertices = copy(g.vertices) order = len(self.vertices) self.order = order self.valuationDomain = copy(g.valuationDomain) Min = self.valuationDomain['min'] Med = self.valuationDomain['med'] Max = self.valuationDomain['max'] if Debug: print('valuationDomain = ', self.valuationDomain) verticesList = [x for x in self.vertices] random.shuffle(verticesList) if Debug: print(verticesList) spannedVertices = set() spanningTree = [] randomWalk = [verticesList[0]] i = 0 while randomWalk[i] != verticesList[-1]: nextVertex = random.choice(list(g.gamma[randomWalk[i]])) randomWalk.append(nextVertex) i += 1 if Debug: print('random walk: ',randomWalk) spanningBranch = self._reduceCycles(randomWalk,Debug=Debug) spanningTree.append(spanningBranch) spannedVertices = spannedVertices | set(spanningBranch) remainingVertices = set(verticesList) - set(spanningBranch) if Debug: print('spanningBranch = ', spanningBranch) print('spannedVertices = ', spannedVertices) print('remainingVertices = ',remainingVertices) print('spanningTree = ', spanningTree) while remainingVertices != set(): remainingVerticesList = list(remainingVertices) random.shuffle(remainingVerticesList) randomWalk = [random.choice(remainingVerticesList)] i = 0 while randomWalk[i] not in spannedVertices: nextVertex = random.choice(list(g.gamma[randomWalk[i]])) randomWalk.append(nextVertex) i += 1 if Debug: print('randomWalk = ',randomWalk) spanningBranch = self._reduceCycles(randomWalk,Debug=Debug) spannedVertices = spannedVertices | set(spanningBranch) remainingVertices = set(remainingVerticesList) - set(spanningBranch) spanningTree.append(spanningBranch) if Debug: print('spanningBranch = ', spanningBranch) print('spannedVertices = ', spannedVertices) print('remainingVertices = ',remainingVertices) print('spanningTree = ', spanningTree) self.edges = copy(g.edges) for edgeKey in self.edges: if self.edges[edgeKey] > Med: self.edges[edgeKey] = Med for branch in spanningTree: for i in range(len(branch)-1): edgeKey = frozenset([branch[i],branch[i+1]]) self.edges[edgeKey] = g.edges[edgeKey] if Debug: print('edges = ', self.edges) self.size = self.computeSize() self.gamma = self.gammaSets() self.dfs = self.depthFirstSearch() self.prueferCode = self.tree2Pruefer() def _reduceCycles(self,randomWalk,Debug=False): if Debug: print('randomWalk', randomWalk) reducedWalk = [randomWalk[0]] n = len(randomWalk) t = 0 while t < n-1 and t < 100: k = t+1 for j in range(t+1,n): if randomWalk[t] == randomWalk[j]: k = j+1 if Debug: print(t, k, j, n, randomWalk[t:k]) reducedWalk.append(randomWalk[k]) if Debug: print('reducedWalk', reducedWalk) t = k return reducedWalk
#----------
[docs]class BestDeterminedSpanningForest(RandomTree): """ Constructing the most determined spanning tree (or forest if not connected) using Kruskal's greedy algorithm on the dual valuation. Example Python session: >>> from graphs import * >>> g = RandomValuationGraph(seed=2) >>> g.showShort() *---- short description of the graph ----* Name : 'randomGraph' Vertices : ['v1', 'v2', 'v3', 'v4', 'v5'] Valuation domain : {'med': Decimal('0'), 'min': Decimal('-1'), 'max': Decimal('1')} Gamma function : v1 -> ['v2', 'v3'] v2 -> ['v4', 'v1', 'v5', 'v3'] v3 -> ['v1', 'v5', 'v2'] v4 -> ['v5', 'v2'] v5 -> ['v4', 'v2', 'v3'] >>> mt = BestDeterminedSpanningForest(g) >>> mt.exportGraphViz('spanningTree',withSpanningTree=True) *---- exporting a dot file for GraphViz tools ---------* Exporting to spanningTree.dot [['v4', 'v2', 'v1', 'v3', 'v1', 'v2', 'v5', 'v2', 'v4']] neato -Tpng spanningTree.dot -o spanningTree.png .. image:: spanningTree.png :alt: 7-cycle instance :width: 300 px :align: center """ def __init__(self,g,seed=None,Debug=False): from copy import deepcopy import random random.seed(seed) self.name= g.name+'_randomSpanningTree' if Debug: print(self.name) self.vertices = deepcopy(g.vertices) verticesList = [x for x in self.vertices] verticesList.sort() order = len(self.vertices) self.order = order self.valuationDomain = deepcopy(g.valuationDomain) Min = self.valuationDomain['min'] Med = self.valuationDomain['med'] Max = self.valuationDomain['max'] if Debug: print('valuationDomain = ', self.valuationDomain) edges = dict() for i in range(order): for j in range(i+1,order): edgeKey = frozenset([verticesList[i],verticesList[j]]) if g.edges[edgeKey] > Med: edges[edgeKey] = Med else: edges[edgeKey] = g.edges[edgeKey] forest = [] weightedEdgesList = [(g.edges[e],e) \ for e in g.edges if g.edges[e] > Med ] weightedEdgesList.sort(reverse=True) if Debug: print(weightedEdgesList) for e in weightedEdgesList: if Debug: print('===>>> ', e) edgeKeys = set([v for v in e[1]]) Included = False connectedSubtrees = [] for tree in forest: if Debug: print(tree) test = tree[1] & edgeKeys if len(test) == 2: # already connected Included = True if Debug: print('already included') break elif len(test) == 1: # extending connectedSubtrees.append(tree) if Debug: print('extending',tree) if not Included: edges[e[1]] = e[0] if Debug: print(e[1],edges[e[1]]) print('connecting subtrees',connectedSubtrees) if connectedSubtrees == []: forest.append([set([e[1]]),edgeKeys]) if Debug: print('Add new subtree') elif len(connectedSubtrees) == 1: connectedSubtrees[0][0].add(e[1]) connectedSubtrees[0][1] = connectedSubtrees[0][1] | edgeKeys if Debug: print('Extending an existing subtree') else: connectedSubtrees[0][0].add(e[1]) connectedSubtrees[0][0] =\ connectedSubtrees[0][0] | connectedSubtrees[1][0] connectedSubtrees[0][1] =\ connectedSubtrees[0][1] | connectedSubtrees[1][1] | edgeKeys forest.remove(connectedSubtrees[1]) if Debug: print('connected two subtrees') print(forest) self.edges = deepcopy(edges) if Debug: print('edges = ',self.edges) self.size = self.computeSize() self.gamma = self.gammaSets(Debug) if Debug: print('gamma = ', self.gamma) self.dfs = self.depthFirstSearch()
[docs]class Q_Coloring(Graph): """ Generate a q-coloring of a Graph instance via a Gibbs MCMC sampler in nSim simulation steps (default = len(graph.edges)). Example 3-coloring of a grid 6x6 : >>> from graphs import * >>> g = GridGraph(n=6,m=6) >>> g.showShort() >>> g.exportGraphViz() *----- show short --------------* Grid graph : grid-6-6 n : 6 m : 6 order : 36 >>> qc = Q_Coloring(g,colors=['gold','lightblue','lightcoral']) Running a Gibbs Sampler for 630 step ! >>> qc.checkFeasibility() The q-coloring with 3 colors is feasible !! >>> qc.exportGraphViz() *---- exporting a dot file for GraphViz tools ---------* Exporting to grid-6-6-qcoloring.dot fdp -Tpng grid-6-6-qcoloring.dot -o grid-6-6-qcoloring.png .. image:: grid-6-6-qcoloring.png :alt: 3 coloring of a 6x6 grid :width: 300 px :align: center """ def __init__(self,g,colors=['gold','lightcoral','lightblue'], nSim=None,maxIter=20,seed=None, Comments=True,Debug=False): from copy import deepcopy self.gClass = g.__class__ self.name = '%s-qcoloring' % g.name if isinstance(g.vertices,dict): self.vertices = deepcopy(g.vertices) else: self.vertices = dict() for v in g.vertices: self.vertices[v] = {'name':v,'shortName':v} self.order = len(self.vertices) self.colors = deepcopy(colors) self.valuationDomain = deepcopy(g.valuationDomain) self.edges = deepcopy(g.edges) self.size = len(self.edges) self.gamma = deepcopy(g.gamma) ## print(self.colors[0]) ## for v in self.vertices: ## self.vertices[v]['color'] = colors[0] if nSim == None: nSim = len(self.edges)*2 self.nSim = nSim infeasibleEdges = set([e for e in self.edges]) _iter = 0 while infeasibleEdges != set() and _iter < maxIter: _iter += 1 print('Iteration:', _iter) self.generateFeasibleConfiguration(seed=seed,Reset=True,Debug=Debug) infeasibleEdges = self.checkFeasibility(Comments=Comments)
[docs] def showConfiguration(self): for v in self.vertices: print(v,self.vertices[v]['color'])
[docs] def generateFeasibleConfiguration(self,Reset=True,nSim=None,seed=None,Debug=False): import random random.seed(seed) if Reset: for v in self.vertices: self.vertices[v]['color'] = self.colors[0] if nSim == None: nSim = self.nSim print('Running a Gibbs Sampler for %d step !' % nSim) for s in range(nSim): verticesKeys = [v for v in self.vertices] v = random.choice(verticesKeys) neighborColors = [self.vertices[x]['color']\ for x in self.gamma[v]] feasibleColors = list(set(self.colors) - set(neighborColors)) try: self.vertices[v]['color'] = random.choice(feasibleColors) except: if Debug: print(s,v,'Warning !! Not feasible coloring') self.vertices[v]['color'] = random.choice(self.colors) if Debug: print(s, v, self.vertices[v]['color'])
[docs] def checkFeasibility(self,Comments=True,Debug=False): infeasibleEdges = set() for e in self.edges: if self.edges[e] > self.valuationDomain['med']: neighborColors = set([self.vertices[x]['color'] for x in e]) if len(neighborColors) < 2: if Debug: print('Infeasible Confirguration !! See edge %s' % e) infeasibleEdges.add(e) else: if Debug: print(e,neighborColors) if Comments: if infeasibleEdges == set(): print('The q-coloring with %d colors is feasible !!'\ % len(self.colors)) else: print('The q-coloring with %d colors is apparently not feasible !!'\ % len(self.colors)) print('Either augment nSim=%d or add one more color'\ % self.nSim) return infeasibleEdges
[docs] def exportGraphViz(self,fileName=None, noSilent=True, graphType='png', graphSize='7,7', layout=None): """ Exports GraphViz dot file for q-coloring drawing filtering. The graph drawing layout is depending on the graph type, but can be forced to either 'fdp', 'circo' or 'neato' with the layout parameter. Example: >>> g = Graph(numberOfVertices=10,edgeProbability=0.4) >>> g.showShort() *---- short description of the graph ----* Name : 'randomGraph' Vertices : ['v1','v10','v2','v3','v4','v5','v6','v7','v8','v9'] Valuation domain : {'max': 1, 'min': -1, 'med': 0} Gamma function : v1 -> ['v7', 'v2', 'v3', 'v5'] v10 -> ['v4'] v2 -> ['v1', 'v7', 'v8'] v3 -> ['v1', 'v7', 'v9'] v4 -> ['v5', 'v10'] v5 -> ['v6', 'v7', 'v1', 'v8', 'v4'] v6 -> ['v5', 'v8'] v7 -> ['v1', 'v5', 'v8', 'v2', 'v3'] v8 -> ['v6', 'v7', 'v2', 'v5'] v9 -> ['v3'] >>> qc = Q_Coloring(g,nSim=1000) Running a Gibbs Sampler for 1000 step ! >>> qc.checkFeasibility() The q-coloring with 3 colors is feasible !! >>> qc.exportGraphViz() *---- exporting a dot file for GraphViz tools ---------* Exporting to randomGraph-qcoloring.dot fdp -Tpng randomGraph-qcoloring.dot -o randomGraph-qcoloring.png .. image:: randomGraph-qcoloring.png :alt: 3-coloring of a random graph :width: 300 px :align: center """ import os if noSilent: print('*---- exporting a dot file for GraphViz tools ---------*') vertexkeys = [x for x in self.vertices] n = len(vertexkeys) edges = self.edges Med = self.valuationDomain['med'] i = 0 if fileName == None: name = self.name else: name = fileName dotName = name+'.dot' if noSilent: print('Exporting to '+dotName) ## if bestChoice != set(): ## rankBestString = '{rank=max; ' ## if worstChoice != set(): ## rankWorstString = '{rank=min; ' fo = open(dotName,'w') fo.write('strict graph G {\n') fo.write('graph [ bgcolor = cornsilk, fontname = "Helvetica-Oblique",\n fontsize = 12,\n label = "') fo.write('\\nGraphs Python module (graphviz), R. Bisdorff, 2015", size="') fo.write(graphSize),fo.write('"];\n') for i in range(n): try: nodeName = str(self.vertices[vertexkeys[i]]['shortName']) except: try: nodeName = self.vertices[vertexkeys[i]]['name'] except: nodeName = str(vertexkeys[i]) node = 'n'+str(i+1)+' [shape = "circle", label = "' +nodeName+'"' node += ', style = "filled", color = %s' \ % self.vertices[vertexkeys[i]]['color'] node += '];\n' fo.write(node) for i in range(n): for j in range(i+1, n): if i != j: edge = 'n'+str(i+1) if edges[frozenset( [vertexkeys[i], vertexkeys[j]])] > Med: edge0 = edge+'-- n'+str(j+1)+' [dir=both,style="setlinewidth(1)",color=black, arrowhead=none, arrowtail=none] ;\n' fo.write(edge0) elif edges[frozenset([vertexkeys[i],vertexkeys[j]])] == Med: edge0 = edge+'-- n'+str(j+1)+' [dir=both, color=grey, arrowhead=none, arrowtail=none] ;\n' fo.write(edge0) fo.write('}\n') fo.close() if layout == None: if self.gClass in (GridGraph,RandomTree): commandString = 'neato -T'+graphType+' ' +dotName+' -o '+name+'.'+graphType elif self.gClass == CycleGraph: commandString = 'circo -T'+graphType+' ' +dotName+' -o '+name+'.'+graphType else: commandString = 'fdp -T'+graphType+' ' +dotName+' -o '+name+'.'+graphType else: commandString = layout+' -T'+graphType+' ' +dotName+' -o '+name+'.'+graphType if noSilent: print(commandString) try: os.system(commandString) except: if noSilent: print('graphViz tools not avalaible! Please check installation.')
[docs]class IsingModel(Graph): """ Specialisation of a Gibbs Sampler for the Ising model Example: >>> from graphs import GridGraph, IsingModel >>> g = GridGraph(n=15,m=15) >>> g.showShort() *----- show short --------------* Grid graph : grid-6-6 n : 6 m : 6 order : 36 >>> im = IsingModel(g,beta=0.3,nSim=100000,Debug=False) Running a Gibbs Sampler for 100000 step ! >>> im.exportGraphViz(colors=['lightblue','lightcoral']) *---- exporting a dot file for GraphViz tools ---------* Exporting to grid-15-15-ising.dot fdp -Tpng grid-15-15-ising.dot -o grid-15-15-ising.png .. image:: grid-15-15-ising.png :alt: ising configuration of a 15x15 grid with beta=0.3 :width: 300 px :align: center """ def __init__(self,g,beta=0, nSim=None, Debug=False): from copy import deepcopy self.gClass = g.__class__ self.name = '%s-ising' % g.name if isinstance(g.vertices,dict): self.vertices = deepcopy(g.vertices) else: self.vertices = dict() for v in g.vertices: self.vertices[v] = {'name':v,'shortName':v} self.order = len(self.vertices) self.valuationDomain = deepcopy(g.valuationDomain) self.edges = deepcopy(g.edges) self.size = g.size self.gamma = deepcopy(g.gamma) for v in self.vertices: self.vertices[v]['spin'] = 0 if nSim == None: nSim = len(self.edges) self.nSim = nSim self.generateSpinConfiguration(beta=beta,Debug=Debug) self.SpinEnergy = self.computeSpinEnergy()/self.size
[docs] def generateSpinConfiguration(self,beta=0,nSim=None,Debug=False): from random import choice, random from math import exp if nSim == None: nSim = self.nSim print('Running a Gibbs Sampler for %d step !' % nSim) for s in range(nSim): verticesKeys = [v for v in self.vertices] v = choice(verticesKeys) plusNeighbors = [x for x in self.gamma[v] if self.vertices[x]['spin'] == 1] nPlus = len(plusNeighbors) minusNeighbors = [x for x in self.gamma[v] if self.vertices[x]['spin'] == -1] nMinus = len(minusNeighbors) numerator = exp(2*beta*(nPlus-nMinus)) threshold = numerator/(numerator+1) U = random() if U < threshold: self.vertices[v]['spin'] = 1 else: self.vertices[v]['spin'] = -1 if Debug: print('Spin energy: %d' % (self.computeSpinEnergy()) ) print('s,v,nPlus,nMinus,numerator,threshold,U,spin\n',\ s,v,nPlus,nMinus,numerator,threshold,U,self.vertices[v]['spin'])
[docs] def computeSpinEnergy(self): """ Spin energy H(c) of a spin configuration is H(c) = -sum_{{x,y} in self.edges}[spin_c(x)*spin_c(y)] """ Hc = 0 for e in self.edges: pair = list(e) x = pair[0] y = pair[1] Hc -= self.vertices[x]['spin'] * self.vertices[y]['spin'] return Hc
[docs] def exportGraphViz(self,fileName=None, noSilent=True, graphType='png', graphSize='7,7', edgeColor='black', colors=['gold','lightblue']): """ Exports GraphViz dot file for Ising models drawing filtering. """ import os if noSilent: print('*---- exporting a dot file for GraphViz tools ---------*') vertexkeys = [x for x in self.vertices] n = len(vertexkeys) edges = self.edges Med = self.valuationDomain['med'] i = 0 if fileName == None: name = self.name else: name = fileName dotName = name+'.dot' if noSilent: print('Exporting to '+dotName) ## if bestChoice != set(): ## rankBestString = '{rank=max; ' ## if worstChoice != set(): ## rankWorstString = '{rank=min; ' fo = open(dotName,'w') fo.write('strict graph G {\n') fo.write('graph [ bgcolor = cornsilk, fontname = "Helvetica-Oblique",\n fontsize = 12,\n label = "') fo.write('\\nGraphs Python module (graphviz), R. Bisdorff, 2015", size="') fo.write(graphSize),fo.write('"];\n') for i in range(n): try: nodeName = str(self.vertices[vertexkeys[i]]['shortName']) except: try: nodeName = self.vertices[vertexkeys[i]]['name'] except: nodeName = str(vertexkeys[i]) node = 'n'+str(i+1)+' [shape = "circle", label = "' +nodeName+'"' if self.vertices[vertexkeys[i]]['spin'] == 1: color=colors[0] elif self.vertices[vertexkeys[i]]['spin'] == -1: color=colors[1] else: color=None if color != None: node += ', style = "filled", color = %s' % color node += '];\n' fo.write(node) for i in range(n): for j in range(i+1, n): if i != j: edge = 'n'+str(i+1) if edges[frozenset( [vertexkeys[i], vertexkeys[j]])] > Med: edge0 = edge+'-- n'+str(j+1)+\ ' [dir=both,style="setlinewidth(1)",color='+edgeColor+\ ' arrowhead=none, arrowtail=none] ;\n' fo.write(edge0) elif edges[frozenset([vertexkeys[i],vertexkeys[j]])] == Med: edge0 = edge+'-- n'+str(j+1)+' [dir=both, color=grey, arrowhead=none, arrowtail=none] ;\n' fo.write(edge0) fo.write('}\n') fo.close() if self.gClass in (GridGraph,TriangulatedGrid,RandomTree): commandString = 'neato -T'+graphType+' ' +dotName+' -o '+name+'.'+graphType elif self.gClass == CycleGraph: commandString = 'circo -T'+graphType+' ' +dotName+' -o '+name+'.'+graphType else: commandString = 'fdp -T'+graphType+' ' +dotName+' -o '+name+'.'+graphType if noSilent: print(commandString) try: os.system(commandString) except: if noSilent: print('graphViz tools not avalaible! Please check installation.')
[docs]class MetropolisChain(Graph): """ Specialisation of the graph class for implementing a generic Metropolis Markov Chain Monte Carlo sampler with a given probability distribution probs = {'v1': x, 'v2': y, ...} Usage example: >>> from graphs import * >>> g = Graph(numberOfVertices=5,edgeProbability=0.5) >>> g.showShort() *---- short description of the graph ----* Name : 'randomGraph' Vertices : ['v1', 'v2', 'v3', 'v4', 'v5'] Valuation domain : {'max': 1, 'med': 0, 'min': -1} Gamma function : v1 -> ['v2', 'v3', 'v4'] v2 -> ['v1', 'v4'] v3 -> ['v5', 'v1'] v4 -> ['v2', 'v5', 'v1'] v5 -> ['v3', 'v4'] >>> probs = {} >>> n = g.order >>> i = 0 >>> verticesList = [x for x in g.vertices] >>> verticesList.sort() >>> for v in verticesList: ... probs[v] = (n - i)/(n*(n+1)/2) ... i += 1 >>> met = MetropolisChain(g,probs) >>> frequency = met.checkSampling(verticesList[0],nSim=30000) >>> for v in verticesList: ... print(v,probs[v],frequency[v]) v1 0.3333 0.3343 v2 0.2666 0.2680 v3 0.2 0.2030 v4 0.1333 0.1311 v5 0.0666 0.0635 >>> met.showTransitionMatrix() * ---- Transition Matrix ----- Pij | 'v1' 'v2' 'v3' 'v4' 'v5' -----|------------------------------------- 'v1' | 0.23 0.33 0.30 0.13 0.00 'v2' | 0.42 0.42 0.00 0.17 0.00 'v3' | 0.50 0.00 0.33 0.00 0.17 'v4' | 0.33 0.33 0.00 0.08 0.25 'v5' | 0.00 0.00 0.50 0.50 0.00 """ def __init__(self,g, probs = None): from copy import deepcopy from random import choice self.name = '%s-metro' % g.name if isinstance(g.vertices,dict): self.vertices = deepcopy(g.vertices) else: self.vertices = dict() for v in g.vertices: self.vertices[v] = {'name':v,'shortName':v} self.order = len(self.vertices) if probs == None: for v in self.vertices: self.vertices[v]['prob'] = 1.0/self.order else: for v in self.vertices: self.vertices[v]['prob'] = probs[v] self.valuationDomain = deepcopy(g.valuationDomain) self.edges = deepcopy(g.edges) self.size = g.size self.gamma = deepcopy(g.gamma) self.transition = self.computeTransitionMatrix()
[docs] def computeTransitionMatrix(self): from decimal import Decimal transition = {} for si in self.vertices: transition[si] = {} di = len(self.gamma[si]) pi = self.vertices[si]['prob'] for sj in self.vertices: if si != sj: if sj in self.gamma[si]: dj = len(self.gamma[sj]) pj = self.vertices[sj]['prob'] transition[si][sj] =\ Decimal( str( min(1.0,(pj*di)/(pi*dj))/di ) ) else: transition[si][sj] = Decimal('0') else: sp = 0.0 for sx in self.gamma[si]: dx = len(self.gamma[sx]) px = self.vertices[sx]['prob'] sp += min(1.0,(px*di)/(pi*dx))/di transition[si][sj] = Decimal(str(1.0 - sp)) return transition
[docs] def MCMCtransition(self,si,Debug=False): from random import random,choice neighborsSi = [x for x in self.gamma[si]] di = len(self.gamma[si]) if di == 0: return si pi = self.vertices[si]['prob'] sj = choice(neighborsSi) dj = len(self.gamma[sj]) pj = self.vertices[sj]['prob'] U = random() threshold = min(1.0,(pj*di)/(pi*dj)) if Debug: print(si,di,pi,sj,dj,pj,U,threshold) if U < threshold: return sj else: return si
[docs] def checkSampling(self,si,nSim): frequency = {} for v in self.vertices: frequency[v] = 0.0 sc = si for i in range(nSim): sc = self.MCMCtransition(sc) frequency[sc] += 1.0 for x in frequency: frequency[x] /= nSim return frequency
[docs] def saveCSVTransition(self,fileName='transition',Debug=False): """Persistent storage of the transition matrix in the form of a csv file. """ import csv from decimal import Decimal if Debug: print('*--- Saving transition matrix P_ij into file: %s.csv> ---*' %\ (fileName)) fileNameExt = str(fileName)+str('.csv') fo = open(fileNameExt, 'w') csvfo = csv.writer(fo,quoting=csv.QUOTE_NONNUMERIC) verticesList = [x for x in self.vertices] verticesList.sort() headerText = ["P_ij"] + verticesList if Debug: print(headerText) csvfo.writerow(headerText) relation = self.transition for x in verticesList: rowText = [x] for y in verticesList: rowText.append( Decimal('%.5f' % float(relation[x][y])) ) if Debug: print(rowText) csvfo.writerow(rowText) fo.close()
[docs] def showTransitionMatrix(self,Sorted=True,\ IntegerValues=False,\ vertices=None,\ relation=None,\ ndigits=2,\ ReflexiveTerms=True): """ Prints on stdout the transition probabilities in vertices X vertices table format. """ if vertices == None: vertices = self.vertices if relation == None: try: relation = self.transition except: relation = self.computeTransitionMatrix() print('* ---- Transition Matrix -----\n', end=' ') print(' S | ', end=' ') verticesList = [] for x in vertices: if isinstance(x,frozenset): try: verticesList += [(vertices[x]['shortName'],x)] except: verticesList += [(verticess[x]['name'],x)] else: verticesList += [(str(x),x)] if Sorted: verticesList.sort() print(verticesList) verticesList.sort() for x in verticesList: print("'"+x[0]+"'\t ", end=' ') print('\n-----|------------------------------------------------------------') for x in verticesList: print("'"+x[0]+"' | ", end=' ') for y in verticesList: if x != y: formatString = '%%2.%df\t' % ndigits print(formatString % (relation[x[1]][y[1]]), end=' ') else: if ReflexiveTerms: formatString = '%%2.%df\t' % ndigits print(formatString % (relation[x[1]][y[1]]), end=' ') else: formatString = ' - \t' print(formatString, end=' ') print() print('\n')
[docs]class MISModel(Graph): """ Specialisation of a Gibbs Sampler for the hard code model, that is a random MIS generator. Example: >>> from graphs import MISModel >>> from digraphs import CirculantDigraph >>> dg = CirculantDigraph(order=15) >>> g = dg.digraph2Graph() >>> g.showShort() *---- short description of the graph ----* Name : 'c15' Vertices : ['1', '10', '11', '12', '13', '14', '15', '2', '3', '4', '5', '6', '7', '8', '9'] Valuation domain : {'med': 0, 'min': -1, 'max': 1} Gamma function : 1 -> ['2', '15'] 10 -> ['11', '9'] 11 -> ['10', '12'] 12 -> ['13', '11'] 13 -> ['12', '14'] 14 -> ['15', '13'] 15 -> ['1', '14'] 2 -> ['1', '3'] 3 -> ['2', '4'] 4 -> ['3', '5'] 5 -> ['6', '4'] 6 -> ['7', '5'] 7 -> ['6', '8'] 8 -> ['7', '9'] 9 -> ['10', '8'] >>> mis = MISModel(g) Running a Gibbs Sampler for 1050 step ! >>> mis.checkMIS() {'2','4','7','9','11','13','15'} is maximal ! >>> mis.exportGraphViz() *---- exporting a dot file for GraphViz tools ---------* Exporting to c15-mis.dot fdp -Tpng c15-mis.dot -o c15-mis.png .. image:: c15-mis.png :alt: 15-cycle with colored MIS :width: 300 px :align: center """ def __init__(self,g, nSim=None, maxIter=20, seed=None, Debug=False): from copy import deepcopy self.gClass = deepcopy(g.__class__) self.name = '%s-mis' % g.name if isinstance(g.vertices,dict): self.vertices = deepcopy(g.vertices) else: self.vertices = dict() for v in g.vertices: self.vertices[v] = {'name':v,'shortName':v} self.order = len(self.vertices) self.valuationDomain = deepcopy(g.valuationDomain) self.edges = deepcopy(g.edges) self.size = g.size self.gamma = deepcopy(g.gamma) if nSim == None: nSim = len(self.edges)*10 self.nSim = nSim unCovered = set([x for x in self.vertices]) _iter = 0 while unCovered != set() and _iter < maxIter: _iter += 1 print('Iteration: ', _iter) self.generateMIS(Reset=True, nSim=nSim,seed=seed,Debug=Debug) mis,misCover,unCovered = self.checkMIS()
[docs] def generateMIS(self,Reset=True,nSim=None,seed=None,Comments=True,Debug=False): import random random.seed(seed) from math import exp if nSim == None: nSim = self.nSim verticesKeys = [v for v in self.vertices] if Reset: for v in verticesKeys: self.vertices[v]['mis'] = 0 if Comments: print('Running a Gibbs Sampler for %d step !' % nSim) for s in range(nSim): v = random.choice(verticesKeys) Potential = True for x in self.gamma[v]: if self.vertices[x]['mis'] == 1: Potential = False break if Potential: self.vertices[v]['mis'] = 1 else: self.vertices[v]['mis'] = -1 if Debug: print('s,v,neighbors,mis\n',\ s,v,self.gamma[v],self.vertices[v]['mis']) self.mis,self.misCover,self.unCovered = self.checkMIS(Comments=Debug)
[docs] def checkMIS(self,Comments=True): """ Verify maximality of independent set. .. note:: Returns three sets: an independent choice, the covered vertices, and the remaining uncovered vertices. When the last set is empty, the independent choice is maximal. """ cover = set() mis = set() for x in self.vertices: if self.vertices[x]['mis'] == 1: cover = cover | self.gamma[x] mis.add(x) misCover = mis | cover unCovered = set(self.vertices.keys()) - misCover if Comments: if unCovered == set(): print(mis,' is maximal !') else: print(mis,' is not maximal, uncovered = ',unCovered) return mis,misCover,unCovered
[docs] def exportGraphViz(self,fileName=None, noSilent=True, graphType='png', graphSize='7,7', misColor='lightblue'): """ Exports GraphViz dot file for MIS models drawing filtering. """ import os if noSilent: print('*---- exporting a dot file for GraphViz tools ---------*') vertexkeys = [x for x in self.vertices] n = len(vertexkeys) edges = self.edges Med = self.valuationDomain['med'] i = 0 if fileName == None: name = self.name else: name = fileName dotName = name+'.dot' if noSilent: print('Exporting to '+dotName) ## if bestChoice != set(): ## rankBestString = '{rank=max; ' ## if worstChoice != set(): ## rankWorstString = '{rank=min; ' fo = open(dotName,'w') fo.write('strict graph G {\n') fo.write('graph [ bgcolor = cornsilk, fontname = "Helvetica-Oblique",\n fontsize = 12,\n label = "') fo.write('\\nGraphs Python module (graphviz), R. Bisdorff, 2015", size="') fo.write(graphSize),fo.write('"];\n') for i in range(n): try: nodeName = str(self.vertices[vertexkeys[i]]['shortName']) except: try: nodeName = self.vertices[vertexkeys[i]]['name'] except: nodeName = str(vertexkeys[i]) node = 'n'+str(i+1)+' [shape = "circle", label = "' +nodeName+'"' if self.vertices[vertexkeys[i]]['mis'] == 1: color=misColor else: color=None if color != None: node += ', style = "filled", color = %s' % color node += '];\n' fo.write(node) for i in range(n): for j in range(i+1, n): if i != j: edge = 'n'+str(i+1) if edges[frozenset( [vertexkeys[i], vertexkeys[j]])] > Med: edge0 = edge+'-- n'+str(j+1)+' [dir=both,style="setlinewidth(1)",color=black, arrowhead=none, arrowtail=none] ;\n' fo.write(edge0) elif edges[frozenset([vertexkeys[i],vertexkeys[j]])] == Med: edge0 = edge+'-- n'+str(j+1)+' [dir=both, color=grey, arrowhead=none, arrowtail=none] ;\n' fo.write(edge0) fo.write('}\n') fo.close() if self.gClass in (GridGraph,RandomTree): commandString = 'neato -T'+graphType+' ' +dotName+' -o '+name+'.'+graphType elif self.gClass == CycleGraph: commandString = 'circo -T'+graphType+' ' +dotName+' -o '+name+'.'+graphType else: commandString = 'fdp -T'+graphType+' ' +dotName+' -o '+name+'.'+graphType if noSilent: print(commandString) try: os.system(commandString) except: if noSilent: print('graphViz tools not avalaible! Please check installation.')
# --------------testing the module ---- if __name__ == '__main__': from graphs import SnakeGraph S = SnakeGraph(p=3,q=7) S.showShort() S.exportGraphViz('4_7_snake',lineWidth=3,arcColor="red") ## from time import time ## #g = GridGraph(4,4) ## g = RandomGraph(order=30,seed=10) ## #g.exportGraphViz() ## #print(g._degreeLabelling()) ## #print(g._triplets(Comments=True)) ## t0 = time();print(len(g._computeChordlessCycles(Cycle3=True,Comments=False)));print(time()-t0) ## t0 = time();print(len(g.computeChordlessCycles(Cycle3=True,Comments=False)));print(time()-t0) #g.save('test') #g = Graph('test') ## ust = RandomSpanningTree(g,seed=200,Debug=False) ## ust.showShort() ## ust.showMore() ## ust.dfs = ust.randomDepthFirstSearch() ## ust.exportGraphViz(withSpanningTree=True) ## print(ust.prueferCode)