# rm(list=ls())
# poisson clusters


# generate n random points on 0:100 x 0:100
n=400
x=runif(n)*100
y=runif(n)*100
dat=data.frame(x,y)

# assign row and column and position for every dot
dat$row=0
dat$col=0
dat$pos=0
for (j in 1:nrow(dat)){
  dat$row=ceiling(dat$y/10)
  dat$col=ceiling(dat$x/10)
  dat$pos=10*(dat$row-1)+dat$col
}

# plot the points
plot(dat$x, dat$y, ylim=c(0,100), xlim=c(0,100), 
     asp=1, frame.plot=FALSE, axes=FALSE,
     xlab=' ', ylab=' ', col='aquamarine', pch=16,
     las=1, xaxs="i", yaxs="i",
)

#draw gridlines
axis(1, at=0:10*10)
axis(2, at=0:10*10, las=1, pos=0)
axis(3, at=0:10*10,  labels=FALSE)
axis(4, at=0:10*10,   labels=FALSE, pos=100)
for (j in 1:9){
  lines(c(0,100), j*c(10,10))
  lines(j*c(10,10), c(0,100))
}

# calculate number of points for each position
dat$count=1
ppp=aggregate(count~pos, dat, sum)
dat=dat[,-6]


# add positions with zero counts, if any
npp=nrow(ppp)
if(npp<100){
  w=which(!(1:100 %in% ppp$pos))
  addrows=(npp+1):(npp+length(w))
  ppp[addrows,1]=w
  ppp[addrows,2]=0
  ppp=ppp[order(ppp$pos),]
}

tb=table(ppp$count)
print(tb)

# identify lasrgest cluster/s
mx=max(ppp$count)
loc=which(ppp$count==mx)
clusters=dat[dat$pos %in% loc,]
points(clusters$x, clusters$y, col='red', pch=16)

# identify second lasrgest cluster/s
loc=which(ppp$count==mx-1)
clusters=dat[dat$pos %in% loc,]
points(clusters$x, clusters$y, col='blue', pch=16)

# identify third lasrgest cluster/s
loc=which(ppp$count==mx-2)
clusters=dat[dat$pos %in% loc,]
points(clusters$x, clusters$y, col='darkcyan', pch=16)

# identify 0/1 cluster/s
loc=which(ppp$count %in% 0:1)
lows=data.frame(loc)
lows$row=ceiling(lows$loc/10)
lows$col=lows$loc %% 10
w=which(lows$col==0)
if(length(w)>0){
  lows$col[w]=10
}
for (j in 1:nrow(lows)){
  h1=(lows$col[j]-1)*10
  h2=h1+10
  v1=(lows$row[j]-1)*10
  v2=v1+10
  lines(c(h1,h2), c(v1,v1), lwd=2, col='purple')
  lines(c(h1,h2), c(v2,v2), lwd=2, col='purple')
  lines(c(h1,h1), c(v1,v2), lwd=2, col='purple')
  lines(c(h2,h2), c(v1,v2), lwd=2, col='purple')
}

# test goodness of fit for poisson distribution
observed=as.numeric(tb)
lambda=nrow(dat)/100
df=length(observed)-1
expected=100*dpois(0:df, lambda)
g2=sum(observed*log(observed/expected))
pv=1-pchisq(g2,9)