# schelling example
# introduction to R programminig
# modificato luglio 2005

# try move to free location
movetofree <- function(i,j){
    z1 <- sample(1:dim,1)
    z2 <- sample(1:dim,1)
    if(m[z1,z2]==0){m[z1,z2] <<- m[i,j];m[i,j] <<- 0
	if(m[z1,z2]==1) incrementap(z1,z2,i,j)
	if(m[z1,z2]==-1)incrementam(z1,z2,i,j)
    }
}

# upgrade mp
incrementap <- function(i,j,u,v){
    zz <- matrix(0,dim,dim)
    zz[(i-1):min(i+1,dim),(j-1):min(j+1,dim)] <- 1
    zz[i,j] <- 0
    mp <<- mp+zz
    zz <- matrix(0,dim,dim)
    zz[(u-1):min(u+1,dim),(v-1):min(v+1,dim)] <- 1
    zz[u,v] <- 0
    mp <<- mp-zz
}

# upgrade mm
incrementam <- function(i,j,u,v){
    zz <- matrix(0,dim,dim)
    zz[(i-1):min(i+1,dim),(j-1):min(j+1,dim)] <- 1
    zz[i,j] <- 0
    mm <<- mm+zz
    zz <- matrix(0,dim,dim)
    zz[(u-1):min(u+1,dim),(v-1):min(v+1,dim)] <- 1
    zz[u,v] <- 0
    mm <<- mm-zz
}

# how many "plus" in a neighborhood
hmanyp <- function(i,j){
    z <- m[(i-1):min(i+1,dim),(j-1):min(j+1,dim)]
    sum(z[z==1])-(m[i,j]==1)
}

# how many "minus" in a neighborhood
hmanym <- function(i,j){
    z <- m[(i-1):min(i+1,dim),(j-1):min(j+1,dim)]
    sum(z[z==-1])+(m[i,j]==-1)
}

happyQ <- function(i,j){
    z <- F
    tot <- mp[i,j]+mm[i,j]
    if(m[i,j]==1 & tot!=0 & mp[i,j]/tot > perc) z <- T
    if(m[i,j]==-1 & tot!=0 & mm[i,j]/tot > perc) z <- T
    z
}

hmuchperc <- function(i,j){
    tot <- mm[i,j]+mp[i,j]
    z <- 0
    if(m[i,j]==1 & tot!=0) z<- mp[i,j]/tot
    if(m[i,j]==-1 & tot!=0) z<- mm[i,j]/tot
    z
}

    

# main

# colori <- c("grey","white","green")
colori <- c("red","black","green")


happyts <- rep(0,100)
dim <- 25
perc <- 0.6
m <- matrix(trunc(runif(dim**2)*3)-1,dim,dim)
numagent <- sum(m!=0)
mp <- matrix(0,dim,dim)
mm <- mp
mh <- mp

for(i in 1:dim)for(j in 1:dim)mp[i,j] <- hmanyp(i,j)
for(i in 1:dim)for(j in 1:dim)mm[i,j] <- -hmanym(i,j)
for(i in 1:dim)for(j in 1:dim)mh[i,j] <- happyQ(i,j)

happy <- 0
image(m,col=colori,main=paste("empty spaces in",colori[2]))
volte <- 1;while(happy!=1){
    for(i in 1:dim)for(j in 1:dim)if(mh[i,j]==F) movetofree(i,j)
    for(i in 1:dim)for(j in 1:dim)mh[i,j] <- happyQ(i,j)
    happy <- sum(mh)/numagent
    happyts[volte] <- happy
    volte <- volte+1
#     image(m,col=c(2,1,3)) # -1 is 2=red, 1 is 1=black, 1 is 3=green    
    image(m,col=colori,add=T)
#     if(volte==2)readline()
    cat("iterazione ", volte, ": happy ",happy,"\n")
}
readline("premi invio")
plot(happyts[happyts>0],t="l",xlab="Iterazioni",ylab="Perc. contenti")

mperc <- matrix(0,dim,dim)
for(i in 1:dim)for(j in 1:dim)mperc[i,j] <- hmuchperc(i,j)
cat("percentuale simili ",sum(mperc)/numagent,"\n")