import numpy as np
import scipy as sp
import networkx as nx
import matplotlib.pyplot as plt
%matplotlib inline
Before we start, look at here
On this seminar we are implementing Schelling's Segregation Model for 2 classes.
To recall the model setting:
We will run it no Knight Moves Network. However for the sake of demonstration it would be better to run it on Regular Lattice with 8 neighbours
def GenKnightNetwork(boardSize):
G = nx.Graph()
pos = {}
for row in xrange(boardSize):
for col in xrange(boardSize):
nodeId = row + col*boardSize
pos[nodeId] = np.array([1.0*row/boardSize, 1.0*col/boardSize])
newPos = GetLegalMoves(row, col, boardSize)
for e in newPos:
nid = e[0] + e[1]*boardSize
G.add_edge(nodeId, nid)
return G, pos
def GetLegalMoves(x,y,boardSize):
newMoves = []
moveOffsets = [(-1,-2),(-1,2),(-2,-1),(-2,1),
( 1,-2),( 1,2),( 2,-1),( 2,1)]
for i in moveOffsets:
newX = x + i[0]
newY = y + i[1]
if isLegalCoord(newX,boardSize) and \
isLegalCoord(newY,boardSize):
newMoves.append((newX,newY))
return newMoves
def isLegalCoord(x,boardSize):
if x >= 0 and x < boardSize:
return True
else:
return False
def SegregationModel(A, modelParams):
# Grabbing model parameters
dens = modelParams.get('density', 0.5)
propA = modelParams.get('propA', 0.5)
thresh = modelParams.get('threshold', 0.5)
maxIters = modelParams.get('maxIters', 100)
utilFunc = modelParams.get('utilFunc', 1)
# Some initializations
# Vertices Number
V = A.shape[0]
nAB = np.round(V*dens)
nA = np.round(nAB*propA)
nB = nAB - nA
nE = V - nAB
# Masks and Indices
vAB,vA,vB,vE = np.zeros((V,1), dtype=bool),np.zeros((V,1), dtype=bool),np.zeros((V,1), dtype=bool),np.zeros((V,1), dtype=bool)
idx = np.random.permutation(V)
idxA = idx[0:nA]
vA[idxA] = 1
idxB = idx[nA:(nAB)]
vB[idxB] = 1
vE = ~(vA | vB)
idxE = idx[(nAB):]
idxAB = np.concatenate((idxA, idxB))
# Calculate neighourhood types for the 1st step
P = np.empty((V, maxIters))
P[idxA,1] = 1.0*A[idxA,].dot(vB).flatten()/np.sum(A[idxA,], axis=1)
P[idxB,1] = 1.0*A[idxB,].dot(vA).flatten()/np.sum(A[idxB,], axis=1)
for t in xrange(1, maxIters):
# Calculate all utilities
utilAll = utilFunc(P[idxAB,t-1], thresh)
# Check utility state
if np.all(utilAll):
break
else: # Choose node to move
utilCur = True
while utilCur:
idxFrom = np.random.choice(idxAB)
utilCur = utilFunc(P[idxFrom, t-1], thresh)
# Choose candidate node
candP = np.empty((V,)) * np.nan
if idxFrom in idxA:
candP[idxE] = 1.0*A[idxE,].dot(vB).flatten()/np.sum(A[idxE,], axis=1)
else:
candP[idxE] = 1.0*A[idxE,].dot(vA).flatten()/np.sum(A[idxE,], axis=1)
# Calculate candidate utility
utilCand = utilFunc(candP, thresh)
# Move node
if np.any(utilCand > utilCur):
idxTo = np.random.choice(np.where(utilCand > utilCur)[0])
if idxFrom in idxA:
vA[idxFrom], vA[idxTo] = vA[idxTo], vA[idxFrom]
idxA[idxA == idxFrom] = idxTo
else:
vB[idxFrom], vB[idxTo] = vB[idxTo], vB[idxFrom]
idxB[idxB == idxFrom] = idxTo
idxAB[idxAB == idxFrom], idxE[idxE == idxTo] = idxTo, idxFrom
# Update neighbourhood state
P[idxA,t] = 1.0*A[idxA,].dot(vB).flatten()/np.sum(A[idxA,], axis=1)
P[idxB,t] = 1.0*A[idxB,].dot(vA).flatten()/np.sum(A[idxB,], axis=1)
print 'Finished on {0} iteration'.format(t)
return vA.astype('double')+vB.astype('double')*2
# Getting the graph
g,pos = GenKnightNetwork(10)
A = np.asarray(nx.adj_matrix(g).todense())
utilFunc = lambda x,p: x<p
modelParams = {}
modelParams['density'] = 0.5
modelParams['propA'] = 0.5
modelParams['threshold'] = 0.5
modelParams['maxIters'] = 100
modelParams['utilFunc'] = utilFunc
col = SegregationModel(A, modelParams)
nx.draw(g, pos=pos, node_color=col)
Finished on 10 iteration
/home/shestakoff/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:1: RuntimeWarning: invalid value encountered in less if __name__ == '__main__':
Implement Freeman's segregation, Schilling’s mixing coefficients and Assortative mixing coefficient. Calculate values results for each iteration of the model.