## REQUIREMENT INSTALL THE (GREAT) R PACKAGE spatstat

require(spatstat)


#########################################################################
### 1. Estimation of the intensity of a homogeneous spatial point process
#########################################################################

## 3 trees datasets. Which one has the highest intensity?

jap<-japanesepines
swe<-unmark(swedishpines)
fin<-unmark(finpines)

jap
swe
fin

par(mfrow=c(1,3),mar=c(2,2,2,2))
plot(jap,main="japanespines")
plot(swe,main="swedishpines")
plot(fin,main="finpines")

## Do the confidence interval overlap?

icRho<-function(X,alpha=.05){
	uni<-X$window$units[3]$mul
	v<-area.owin(X)* uni^2 
	rhoEst<-npoints(X)/v
	(sqrt(rhoEst)+ c(1,-1)* qnorm(alpha/2)/2/sqrt(v))^2->ic
	list(rhoHat=rhoEst,ic=ic)

}

icRho(jap)
icRho(swe)
icRho(fin)

## Quadrat based test to test H1: X does not follow a staionary Poisson point process.

par(mfrow=c(1,1))
quadrat.test(jap,nx=3)->q;plot(q);points(jap);q
quadrat.test(swe,nx=3)->q;plot(q);points(swe);q
quadrat.test(fin,nx=3)->q;plot(q);points(fin);q

####################################################################################
## 2. Intensity estimation of inhomogeneous spatial point processes
###################################################################################

## dataset bei: is this point process stationary or Poisson?

bei
par(mfrow=c(1,2))
plot(bei,cex=.4)
quadrat.test(bei,nx=3)->q;plot(q,lwd=2,col='red');points(bei,cex=.1);q

## non-parametric estimation of the intensity

density(bei)->dbei
par(mfrow=c(1,1));plot(dbei);points(bei,cex=.1);q

## Plot with extra covariates: elevation field and slope of elevation.

image(bei.extra$elev);points(bei,pch=3,cex=.4)
image(bei.extra$grad);points(bei,pch=3,cex=.4)

## Parametric estimation of the intensity: log-linear model for the intensity
## Estimation using the "Poisson likelihood" contrast function (by default)

ppm(bei,~elev+grad,covariates=list(elev=bei.extra$elev,grad=bei.extra$grad),nd=200)->ppmbei
ppmbei

par(mfrow=c(1,2))
plot(ppmbei,ngrid=c(128,128),se=FALSE,superimpose=FALSE,pause=FALSE)->trend
points(bei,pch=3,cex=.4)

plot(ppmbei,ngrid=c(128,128),se=TRUE,trend=FALSE,superimpose=FALSE,pause=FALSE)->se
points(bei,pch=3,cex=.4)

par(mfrow=c(1,1))
Mdbei<-as.matrix(dbei)
Mtrendbei<-as.matrix(trend$trend[[1]])
image(as.im(Mdbei-Mtrendbei),main='difference between non parametric and parametric estimations')


## Do simulated replications of Poisson point processes using the estimated intensity look like the dataset? 

rhoEst<-function(x,y,coef=ppmbei$coef){
	r<-coef[1]+interp.im(bei.extra$elev,x,y)*coef[2] + interp.im(bei.extra$grad,x,y)*coef[3]
	exp(r)
}


par(mfrow=c(3,3))
for (i in 1:9){
sim<-rpoispp(rhoEst,lmax=0.03,win=as.owin(bei))
if (i==5) plot(bei,pch=3,cex=.4,main="") else plot(sim,pch=3,cex=.4,main="")
}
par(mfrow=c(1,1))

## Chorley dataset

cases<-unmark(chorley)
incin<-as.ppp(chorley.extra$incin,W=as.owin(chorley))
# dIncin<-crossdist(rpoispp(lam=1e5/area.owin(chorley),win=as.owin(chorley) ))
dIncin<-function(x,y) sqrt( (incin$x-x)^2+(incin$y-y)^2 )

ppmcases<- ppm(cases,~1+Z,covariates=list(Z=dIncin))
ppmcases
plot(ppmcases,super=FALSE,se=FALSE,pause=FALSE,ngrid=c(128,128))
points(cases,pch=3,cex=.5);points(incin,col="red",cex=2,pch=19)
plot(density(cases))
points(cases,pch=3,cex=.5) ;points(incin,col="red",cex=2,pch=19)



########################################
########## Summary statistics
########################################


## dataset cells:

X<-cells
plot(X)
cells
intensity(cells)
plot(density(cells))

rr<-seq(0,.15,l=100)
plot(pcf(cells,r=rr))

plot(envelope(cells,pcf,r=rr,nsim=500,nrank=trunc(.025*500)))
mad.test(cells,rinterval=range(rr))

plot(Kest(cells,r=rr))
plot(envelope(cells,Kest,r=rr,nsim=500,nrank=trunc(.025*500),savepatterns=TRUE)->eX)
mad.test(eX,rinterval=range(rr))

plot(Lest(cells,r=rr))
plot(envelope(cells,Lest,r=rr,nsim=500,nrank=trunc(.025*500)))

plot(Fest(cells,r=rr))
plot(envelope(cells,Fest,r=rr,nsim=500,nrank=trunc(.025*500)))

plot(Gest(cells,r=rr))
plot(envelope(cells,Gest,r=rr,nsim=500,nrank=trunc(.025*500)))

plot(Jest(cells,r=rr))
plot(envelope(cells,Jest,r=rr,nsim=500,nrank=trunc(.025*500),savepatterns=TRUE)->eX)
mad.test(eX,rinterval=range(rr))



##### dataset bei:

plot(bei)

ns<-100

rr<-seq(0,100,l=100)
plot(Lest(bei,r=rr))
plot(envelope(bei,Lest,r=rr,nsim=ns,nrank=trunc(.025*ns),savepatterns=TRUE)->eX)
# mad.test(eX,rinterval=range(rr))

rhoEstbei<-trend$trend[[1]]
plot(Linhom(bei,lambda=rhoEstbei,r=rr))
# mad.test(eX,rinterval=range(rr))

rhoEstbei2<-density(bei)
plot(Linhom(bei,lambda=rhoEstbei2,r=rr))
# mad.test(eX,rinterval=range(rr))


##########################################################
####  Gibbs point processes estimation
##########################################################


## dataset cells modeled by a stationary Strauss point process

plot(cells)
R<-.12
ppm.cells<-ppm(cells,inter=Strauss(r=R),method="mpl",nd=200)
best<-exp(ppm.cells$coef)[1]
gest<-exp(ppm.cells$coef)[2]
mod <- list(cif="strauss",par=list(beta=best,gamma=gest,r=R),w=c(0,1,0,1))
par(mfrow=c(3,3))
for (i in 1:9){
	X<- rmh(model=mod,start=list(n.start=400),control=list(nrep=1e6,nverb=1e5))
	if (i==5) plot(cells) else plot(X)
}

par(mfrow=c(1,1))
rr<-seq(0,.2,l=100)
plot(Lest(cells,r=rr)->L.cells)
Lest.X<-NULL
for (i in 1:100){
	cat('\r i:',i)
	X<- rmh(model=mod,start=list(n.start=400),control=list(nrep=1e6,nverb=1e6))
	Lest.X<-cbind(Lest.X,Lest(X,r=rr)$border)
}
q025<-apply(Lest.X,1,quantile,0.025)
q975<-apply(Lest.X,1,quantile,0.975)
points(rr,q025,col="pink",lwd=2,type="l")
points(rr,q975,col="pink",lwd=2,type="l")


#### Is the Strauss Hard-core model better for this dataset?

ppm.cells<-ppm(cells,inter=StraussHard(hc=.08,r=R),method="mpl",nd=200)
best<-exp(ppm.cells$coef)[1]
gest<-exp(ppm.cells$coef)[2]
mod <- list(cif="straush",par=list(beta=best,gamma=gest,r=R,hc=.08),w=c(0,1,0,1))
par(mfrow=c(3,3))
for (i in 1:9){
	X<- rmh(model=mod,start=list(n.start=400),control=list(nrep=1e6,nverb=1e5))
	if (i==5) plot(cells) else plot(X)
}

par(mfrow=c(1,1))
rr<-seq(0,.2,l=100)
plot(Lest(cells,r=rr)->L.cells)
Lest.X<-NULL
for (i in 1:100){
	cat('\r i:',i)
	X<- rmh(model=mod,start=list(n.start=400),control=list(nrep=1e6,nverb=1e6))
	Lest.X<-cbind(Lest.X,Lest(X,r=rr)$border)
}
q025<-apply(Lest.X,1,quantile,0.025)
q975<-apply(Lest.X,1,quantile,0.975)
points(rr,q025,col="pink",lwd=2,type='l')
points(rr,q975,col="pink",lwd=2,type='l')