############################
# R 2.4.0.                 #
# spatstat spatstat 1.10-3 # 
############################


rm(list=ls(all=TRUE))

################################
# R Code : Firemen emergencies #
################################

# Load packages #
library(spatstat) # rlabel function
#library(splancs) # rlabel function / Just used for computing the optimal h by the Diggle's method
library(KernSmooth)
library(ks)
library(RandomFields)
library(grDevices)
library(geoR) # boxcox.fit
library(car) # box.cox
library(evd)


# Load data #
load("CSBIGS.Rdata")

ls()

# Emergencies : data frame of emergencies - locations (X,Y), workload mark (M), time of occurence (T), month (Month).
# Region : Window of observation.
# Pop : Population in the window observation for each administrative unit. 


# Initialisation #
perturb <- TRUE
disp <- 50
titleplot <- ""

if (perturb) {
# Perturbation of locations # 
x <- Emergencies[,1] + rnorm(length(Emergencies[,1]), mean=0, sd=disp)
y <- Emergencies[,2] + rnorm(length(Emergencies[,1]), mean=0, sd=disp)
ok <- inside.owin(x,y,Region)
Emergencies <- Emergencies[ok,]
Emergencies[,1] <- x[ok]
Emergencies[,2] <- y[ok] 
}


# Marked Point Patterns #

Emergencies01 <- Emergencies[(Emergencies[,5]==1),] 
Emergencies02 <- Emergencies[(Emergencies[,5]==2),]
Emergencies03 <- Emergencies[(Emergencies[,5]==3),] 
Emergencies04 <- Emergencies[(Emergencies[,5]==4),]
Emergencies05 <- Emergencies[(Emergencies[,5]==5),] 
Emergencies06 <- Emergencies[(Emergencies[,5]==6),]
Emergencies07 <- Emergencies[(Emergencies[,5]==7),] 
Emergencies08 <- Emergencies[(Emergencies[,5]==8),]
Emergencies09 <- Emergencies[(Emergencies[,5]==9),] 
Emergencies10 <- Emergencies[(Emergencies[,5]==10),]
Emergencies11 <- Emergencies[(Emergencies[,5]==11),] 
Emergencies12 <- Emergencies[(Emergencies[,5]==12),]

pworkload <- ppp(Emergencies[,1],Emergencies[,2],window=Region,marks=Emergencies[,3])
ptime <- ppp(Emergencies[,1],Emergencies[,2],window=Region,marks=Emergencies[,4])
p01 <- ppp(Emergencies01[,1],Emergencies01[,2],window=Region,marks=Emergencies01[,3])
p02 <- ppp(Emergencies02[,1],Emergencies02[,2],window=Region,marks=Emergencies02[,3])
p03 <- ppp(Emergencies03[,1],Emergencies03[,2],window=Region,marks=Emergencies03[,3])
p04 <- ppp(Emergencies04[,1],Emergencies04[,2],window=Region,marks=Emergencies04[,3])
p05 <- ppp(Emergencies05[,1],Emergencies05[,2],window=Region,marks=Emergencies05[,3])
p06 <- ppp(Emergencies06[,1],Emergencies06[,2],window=Region,marks=Emergencies06[,3])
p07 <- ppp(Emergencies07[,1],Emergencies07[,2],window=Region,marks=Emergencies07[,3])
p08 <- ppp(Emergencies08[,1],Emergencies08[,2],window=Region,marks=Emergencies08[,3])
p09 <- ppp(Emergencies09[,1],Emergencies09[,2],window=Region,marks=Emergencies09[,3])
p10 <- ppp(Emergencies10[,1],Emergencies10[,2],window=Region,marks=Emergencies10[,3])
p11 <- ppp(Emergencies11[,1],Emergencies11[,2],window=Region,marks=Emergencies11[,3])
p12 <- ppp(Emergencies12[,1],Emergencies12[,2],window=Region,marks=Emergencies12[,3])

rm(Emergencies01,Emergencies02,Emergencies03,Emergencies04,Emergencies05,Emergencies06,Emergencies07,Emergencies08,Emergencies09,Emergencies10,Emergencies11,Emergencies12)


# Intensities #

# Choice of h #
# Diggle method - Splancs Package #
Nbpas <- 1000
poly <- data.frame(x=p06$window$bdry[[1]]$x,y=p06$window$bdry[[1]]$y)
Mse2d <- mse2d(as.points(unmark(p06)), poly, nsmse=Nbpas, range=6000)
plot(Mse2d$h[1:Nbpas],Mse2d$mse[1:Nbpas], type="l")
# Problem : close to zero #

# Plug-in #
matvarcov06 <- Hpi.diag(cbind(p06$x,p06$y),binned=TRUE)

h <- 800
epsgrid <- 200

intens <- density.ppp(pworkload,sigma=h,eps=epsgrid)
intens01 <- density.ppp(p01,sigma=h,eps=epsgrid)
intens02 <- density.ppp(p02,sigma=h,eps=epsgrid)
intens03 <- density.ppp(p03,sigma=h,eps=epsgrid)
intens04 <- density.ppp(p04,sigma=h,eps=epsgrid)
intens05 <- density.ppp(p05,sigma=h,eps=epsgrid)
intens06 <- density.ppp(p06,sigma=h,eps=epsgrid)
intens07 <- density.ppp(p07,sigma=h,eps=epsgrid)
intens08 <- density.ppp(p08,sigma=h,eps=epsgrid)
intens09 <- density.ppp(p09,sigma=h,eps=epsgrid)
intens10 <- density.ppp(p10,sigma=h,eps=epsgrid)
intens11 <- density.ppp(p11,sigma=h,eps=epsgrid)
intens12 <- density.ppp(p12,sigma=h,eps=epsgrid)

# Function in order to avoid points outside the pixel image of the intensity #

no.na <- function(p,intens){
p <- p[!(is.na(intens[p,drop=FALSE]))]
}

pworkload <- no.na(pworkload,intens)
ptime <- no.na(ptime,intens)
p01 <- no.na(p01,intens01)
p02 <- no.na(p02,intens02)
p03 <- no.na(p03,intens03)
p04 <- no.na(p04,intens04)
p05 <- no.na(p05,intens05)
p06 <- no.na(p06,intens06)
p07 <- no.na(p07,intens07)
p08 <- no.na(p08,intens08)
p09 <- no.na(p09,intens09)
p10 <- no.na(p10,intens10)
p11 <- no.na(p11,intens11)
p12 <- no.na(p12,intens12)


# 2.1. Estimation of the intensity lambda #
###########################################

# Choice of a month for exploratory analysis #
pmonth <- p06

x11()
#postscript("Graphiques/PPJune.eps", horizontal=FALSE)
#titleplot <- "Locations of emergencies in June"
plot(unmark(pmonth),main=titleplot)
#dev.off()

# Transformation function for the intensity estimation #
# 2 possible choices : choice=0 : log
#                      choice=p (p>0) : ^(1/p)

Transfo <- function(intens,choice){
value <- intens
if (choice==0) {
	ok <- value$v<=0
	if (sum(ok,na.rm=TRUE)!=0) {value$v[ok] <- min(value$v[(ok==0)],na.rm=TRUE)}
	value$v <- log(value$v)}
else {value$v <- (value$v)^(1/choice)}
value
}

x11()
#postscript("Graphiques/IntensJune.eps", horizontal=FALSE)
#titleplot <- "Intensity estimation in June"
plot(Transfo(intens06,0),main=titleplot)
#dev.off()



###########################################
# 3. Time variation and spatial variation #
###########################################

## Comparison of month intensities ##
choice <- 0

# Minimum and maximum values of the ribbon #
Zlimmin <- min(Transfo(intens01,choice)$v,Transfo(intens02,choice)$v,Transfo(intens03,choice)$v,Transfo(intens04,choice)$v,Transfo(intens05,choice)$v,Transfo(intens06,choice)$v,Transfo(intens07,choice)$v,
Transfo(intens08,choice)$v,Transfo(intens09,choice)$v,Transfo(intens10,choice)$v,Transfo(intens11,choice)$v,Transfo(intens12,choice)$v,na.rm=TRUE)

Zlimmax <- max(Transfo(intens01,choice)$v,Transfo(intens02,choice)$v,Transfo(intens03,choice)$v,Transfo(intens04,choice)$v,Transfo(intens05,choice)$v,Transfo(intens06,choice)$v,Transfo(intens07,choice)$v,
Transfo(intens08,choice)$v,Transfo(intens09,choice)$v,Transfo(intens10,choice)$v,Transfo(intens11,choice)$v,Transfo(intens12,choice)$v,na.rm=TRUE)

Zlim <- c(Zlimmin,Zlimmax)


x11()
#postscript("Graphiques/IntensMonthTrans1.eps", horizontal=FALSE)
par(mfrow=c(3,2))
plot(Transfo(intens01,choice),main="January",sub=substitute(paste("N = ", n), list(n=p01$n)),zlim=Zlim)
plot(Transfo(intens02,choice),main="February",sub=substitute(paste("N = ", n), list(n=p02$n)),zlim=Zlim)
plot(Transfo(intens03,choice),main="March",sub=substitute(paste("N = ", n), list(n=p03$n)),zlim=Zlim)
plot(Transfo(intens04,choice),main="April",sub=substitute(paste("N = ", n), list(n=p04$n)),zlim=Zlim)
plot(Transfo(intens05,choice),main="May",sub=substitute(paste("N = ", n), list(n=p05$n)),zlim=Zlim)
plot(Transfo(intens06,choice),main="June",sub=substitute(paste("N = ", n), list(n=p06$n)),zlim=Zlim)
#dev.off()
x11()
#postscript("Graphiques/IntensMonthTrans2.eps", horizontal=FALSE)
par(mfrow=c(3,2))
plot(Transfo(intens07,choice),main="July",sub=substitute(paste("N = ", n), list(n=p07$n)),zlim=Zlim)
plot(Transfo(intens08,choice),main="August",sub=substitute(paste("N = ", n), list(n=p08$n)),zlim=Zlim)
plot(Transfo(intens09,choice),main="September",sub=substitute(paste("N = ", n), list(n=p09$n)),zlim=Zlim)
plot(Transfo(intens10,choice),main="October",sub=substitute(paste("N = ", n), list(n=p10$n)),zlim=Zlim)
plot(Transfo(intens11,choice),main="November",sub=substitute(paste("N = ", n), list(n=p11$n)),zlim=Zlim)
plot(Transfo(intens12,choice),main="December",sub=substitute(paste("N = ", n), list(n=p12$n)),zlim=Zlim)
#dev.off()



## Ratio between time variation and spatial variation ##

# Variability in time #
# Mean intensity on time for each point of the grid #
MeanIntens <- intens
MeanIntens$v <- (intens01$v + intens02$v + intens03$v + intens04$v + intens05$v + intens06$v + intens07$v + intens08$v + intens09$v + intens10$v + intens11$v + intens12$v)/12

# Mean square error between intensity by month and mean intensity #
MSTime <- intens
MSTime$v <- ( (intens01$v-MeanIntens$v)^2 + (intens01$v-MeanIntens$v)^2 + (intens01$v-MeanIntens$v)^2
+ (intens02$v-MeanIntens$v)^2 + (intens03$v-MeanIntens$v)^2 + (intens04$v-MeanIntens$v)^2
+ (intens05$v-MeanIntens$v)^2 + (intens06$v-MeanIntens$v)^2 + (intens07$v-MeanIntens$v)^2
+ (intens08$v-MeanIntens$v)^2 + (intens09$v-MeanIntens$v)^2 + (intens10$v-MeanIntens$v)^2
+ (intens11$v-MeanIntens$v)^2 + (intens12$v-MeanIntens$v)^2 )/12


# Variability in space #
# Mean square error between year intensity and constant intensity (homogeneous case) for each point #
WindowArea <- area.owin(Region)
MSSpace <- intens
MSSpace$v <- (intens$v-(pworkload$n/WindowArea))^2 

# Average mean square error #
AMSSpace <- mean(MSSpace$v,na.rm=TRUE)


# Ratio : Variability in time (MSTime) / Variability in space (AMSSpace) #
Ratio <- intens
Ratio$v <- MSTime$v / AMSSpace
x11()
#postscript("Graphiques/LogImageRatio.eps", horizontal=FALSE)
#titleplot <- "Logarithm of the ratio of measures of variabilities \n between space and time"
plot(Transfo(Ratio,0),main=titleplot)
#dev.off()



## Separability test : positions * time ##

# Data : Subample #
Nech <- 1000 # 1000 : 2628.23 secs
sptime <- ptime[sample(ptime$n,Nech)]
data <- cbind(sptime$x,sptime$y,sptime$marks)

# Regular grid of points #
lam <- density.ppp(sptime,sigma=h,eps=epsgrid)
nx <- lam$dim[1]; ny <- lam$dim[2]; nt <- 6
sxy <- gridcentres(Region,ny,nx)
sxyrep <- cbind(rep(sxy$x,nt),rep(sxy$y,nt))
st <- seq(min(sptime$marks), max(sptime$marks), length = 2 * nt + 1)[2 * (1:nt)]
strep <- rep(st, each=(nx*ny))
s <- cbind(sxyrep,strep)

# Bandwidth #
hx <- hy <- h
ht <- dpik(sptime$marks)      # (max(sptime$marks) - min(sptime$marks))/nt
hxyt <- c(hx,hy,ht)

# 3D-Density estimation #
system.time(dens3d <- kde(data,H=diag(hxyt)^2,eval.points=s))[3]

# dens3d Array #
val <- dens3d$estimate
ok <- !inside.owin(sxyrep[,1],sxyrep[,2],Region)
val[ok] <- NA
dens3dArray <- aperm(array(val,c(ny,nx,nt)),c(2,1,3))

# Positions #
# Density estimation of positions #
sxyVec <- cbind(sxy$x,sxy$y)
densPos <- kde(data[,1:2],H=diag(c(hx^2,hy^2)),eval.points=sxyVec)

# densPos Array #
val <- densPos$estimate
ok <- !inside.owin(sxy$x,sxy$y,Region)
val[ok] <- NA
densPosArray <- aperm(array(val,c(ny,nx)),c(2,1))

# Time #
densT <- kde(data[,3],h=ht,eval.points=st)

# lambdaHat #
lambdaHat <- Nech*dens3dArray

# lambdaTilde #
valPos <- rep(densPos$estimate,nt)
valT <- rep(densT$estimate,nx*ny)
val <- valPos*valT
ok <- !inside.owin(sxyrep[,1],sxyrep[,2],Region)
val[ok] <- NA
densTildeArray <- aperm(array(val,c(ny,nx,nt)),c(2,1,3))

lambdaTilde <- Nech*densTildeArray

# Separability Test #
S1PosTime <- max(abs(lambdaHat-lambdaTilde)/sqrt(lambdaTilde),na.rm=TRUE)
S1PosTime
S2PosTime <- min(abs(lambdaHat-lambdaTilde)/sqrt(lambdaTilde),na.rm=TRUE)
S2PosTime
S5PosTime <- mean((lambdaHat-lambdaTilde)^2,na.rm=TRUE)
S5PosTime
S6PosTime <- max((lambdaHat-lambdaTilde)^2,na.rm=TRUE)
S6PosTime


# Monte-Carlo Test with a Poisson point process model assumption #

# lambdaPos #
lambdaPos <- lam
lambdaPos$v <- Nech*aperm(array(densPos$estimate,c(ny,nx)),c(2,1))

# A pixel image #
ppt <- runifpoint(100,win=owin(c(min(densT$x),max(densT$x)),c(min(densT$x),max(densT$x))) )
lamt <- density.ppp(ppt,sigma=ht,dimyx=nt)

lambdat <- lamt
lambdat$v <- rep(1/nt,nt)%*%t(densT$estimate)  

# lambdaTilde #
lambdaTilde <- array(dim=c(nx,ny,nt))
for (j in 1:nt) {
lambdaTilde[,,j] <- lambdaPos$v * densT$estimate[j]
}

nsim <- 19
S1PosTimeSimul <- rep(0,nsim)
S2PosTimeSimul <- rep(0,nsim)
S5PosTimeSimul <- rep(0,nsim)
S6PosTimeSimul <- rep(0,nsim)

system.time(for (i in 1:nsim) {

cat(paste(i, if (i < nsim) ", "
            else ".", if (i%%10 == 0) 
                "\n", sep = ""))

# Simulation of a point pattern of positions #
simulxy <- rpoispp(lambdaPos)

# Simulation of occurences #
simult <- rpoint(simulxy$n,lambdat)$x

# Data #
data <- cbind(simulxy$x,simulxy$y,simult)

# 3D-Density estimation #
dens3d <- kde(data,H=diag(hxyt)^2,eval.points=s)
val <- dens3d$estimate
ok <- !inside.owin(sxyrep[,1],sxyrep[,2],Region)
val[ok] <- NA
dens3dArray <- aperm(array(val,c(ny,nx,nt)),c(2,1,3))

# lambdaHat #
lambdaHat <- Nech*dens3dArray

S1PosTimeSimul[i] <- max(abs(lambdaHat-lambdaTilde)/sqrt(lambdaTilde),na.rm=TRUE)
S2PosTimeSimul[i] <- min(abs(lambdaHat-lambdaTilde)/sqrt(lambdaTilde),na.rm=TRUE)
S5PosTimeSimul[i] <- mean((lambdaHat-lambdaTilde)^2,na.rm=TRUE)
S6PosTimeSimul[i] <- max((lambdaHat-lambdaTilde)^2,na.rm=TRUE)

})[3]


if ( S1PosTime < max(S1PosTimeSimul,na.rm=TRUE) ) { print("S1PosTime : Accept the separability assumption")
} else { print("S1PosTime : Reject the separability assumption") }
if ( S2PosTime < max(S2PosTimeSimul,na.rm=TRUE) ) { print("S2PosTime : Accept the separability assumption")
} else { print("S2PosTime : Reject the separability assumption") }
if ( S5PosTime < max(S5PosTimeSimul,na.rm=TRUE) ) { print("S5PosTime : Accept the separability assumption")
} else { print("S5PosTime : Reject the separability assumption") }
if ( S6PosTime < max(S6PosTimeSimul,na.rm=TRUE) ) { print("S6PosTime : Accept the separability assumption")
} else { print("S6PosTime : Reject the separability assumption") }




####################################
# 4. Analysis of the workload mark #
####################################

logworkload <- log(pworkload$marks)

# 4.1. Dependence between marks and position #
##############################################

## Discretization by minimising the variance inside each group with the kmeans method ##
cl <- kmeans(logworkload,3,iter.max=100,nstart=500)
tri <- order(c(min(logworkload[cl$cluster==1]),min(logworkload[cl$cluster==2]),min(logworkload[cl$cluster==3])))
bornes <- quantile(logworkload,c(0,1/3,2/3,1))
names(bornes) <- NULL
bornes[2:3] <- c(min(logworkload[cl$cluster==tri[2]]),min(logworkload[cl$cluster==tri[3]]))
bornes <- bornes + c(-0.1,0,0,0.1)
pp <- ppp(x=pworkload$x,y=pworkload$y,marks=logworkload,window=pworkload$w)
Multitypepp <- cut(pp,bornes,labels=c('Low','Medium','High'))
rm(cl,tri)


## Comparison of the 3 intensities corresponding to each level of the created multitype point pattern ##
Low <- split(Multitypepp)$Low
Medium <- split(Multitypepp)$Medium
High <- split(Multitypepp)$High

intensLow <- density.ppp(Low,sigma=h,eps=epsgrid)
intensMedium <- density.ppp(Medium,sigma=h,eps=epsgrid)
intensHigh <- density.ppp(High,sigma=h,eps=epsgrid)

choice <- 0

# Minimum and maximum values of the ribbon #
Zlimmin <- min(Transfo(intensLow,choice)$v,Transfo(intensMedium,choice)$v,Transfo(intensHigh,choice)$v,na.rm=TRUE)

Zlimmax <- max(Transfo(intensLow,choice)$v,Transfo(intensMedium,choice)$v,Transfo(intensHigh,choice)$v,na.rm=TRUE)

Zlim <- c(Zlimmin,Zlimmax)

x11()
par(mfrow=c(1,3))

#postscript("Graphiques/intensMulti1.eps", horizontal=FALSE)
#titleplot <- "Intensity estimate of emergencies with low mark"
plot(Transfo(intensLow,choice),main=titleplot,zlim=Zlim)
#dev.off()

#postscript("Graphiques/intensMulti2.eps", horizontal=FALSE)
#titleplot <- "Intensity estimate of emergencies with medium mark"
plot(Transfo(intensMedium,choice),main=titleplot,zlim=Zlim)
#dev.off()

#postscript("Graphiques/intensMulti3.eps", horizontal=FALSE)
#titleplot <- "Intensity estimate of emergencies with high mark"
plot(Transfo(intensHigh,choice),main=titleplot,zlim=Zlim)
#dev.off()


## separability test ##
# Data : June #
pwork <- pmonth

# Data : Subsample #
Nech <- 100  #pmonth$n   #  pmonth$n : 2503.5 secs  100
pwork <- pwork[sample(pwork$n,Nech)]


data <- cbind(pwork$x,pwork$y,log(pwork$marks))

# Regular grid of points #
lam <- density.ppp(pwork,sigma=h,eps=epsgrid)
nx <- lam$dim[1]; ny <- lam$dim[2]; nm <- 3
sxy <- gridcentres(Region,ny,nx)
sxyrep <- cbind(rep(sxy$x,nm),rep(sxy$y,nm))
sm <- seq(min(data[,3]), max(data[,3]), length = 2 * nm + 1)[2 * (1:nm)]
smrep <- rep(sm, each=(nx*ny))
s <- cbind(sxyrep,smrep)

# Bandwidth #
hx <- hy <- h
hm <- dpik(data[,3])
hxym <- c(hx,hy,hm)

# 3D-Density estimation #
system.time(dens3d <- kde(data,H=diag(hxym)^2,eval.points=s))[3]

# dens3d Array #
val <- dens3d$estimate
ok <- !inside.owin(sxyrep[,1],sxyrep[,2],Region)
val[ok] <- NA
dens3dArray <- aperm(array(val,c(ny,nx,nm)),c(2,1,3))

# Positions #
# Density estimation of positions #
sxyVec <- cbind(sxy$x,sxy$y)
densPos <- kde(data[,1:2],H=diag(c(hx^2,hy^2)),eval.points=sxyVec)

# densPos Array #
val <- densPos$estimate
ok <- !inside.owin(sxy$x,sxy$y,Region)
val[ok] <- NA
densPosArray <- aperm(array(val,c(ny,nx)),c(2,1))

# Marks #
densM <- kde(data[,3],h=hm,eval.points=sm)

# lambdaHat #
lambdaHat <- Nech*dens3dArray

# lambdaTilde #
valPos <- rep(densPos$estimate,nm)
valM <- rep(densM$estimate,nx*ny)
val <- valPos*valM
ok <- !inside.owin(sxyrep[,1],sxyrep[,2],Region)
val[ok] <- NA
densTildeArray <- aperm(array(val,c(ny,nx,nm)),c(2,1,3))

lambdaTilde <- Nech*densTildeArray

# Separability Test #
S1PosMark <- max(abs(lambdaHat-lambdaTilde)/sqrt(lambdaTilde),na.rm=TRUE)
S1PosMark
S2PosMark <- min(abs(lambdaHat-lambdaTilde)/sqrt(lambdaTilde),na.rm=TRUE)
S2PosMark
S5PosMark <- mean((lambdaHat-lambdaTilde)^2,na.rm=TRUE)
S5PosMark
S6PosMark <- max((lambdaHat-lambdaTilde)^2,na.rm=TRUE)
S6PosMark


# Monte-Carlo Test with a Poisson point process model assumption #
# lambdaPos #
lambdaPos <- lam
lambdaPos$v <- Nech*aperm(array(densPos$estimate,c(ny,nx)),c(2,1))

# A pixel image #
ppm <- runifpoint(100,win=owin(c(min(densM$x),max(densM$x)),c(min(densM$x),max(densM$x))) )
lamm <- density.ppp(ppm,sigma=hm,dimyx=nm)

lambdam <- lamm
lambdam$v <- rep(1/nm,nm)%*%t(densM$estimate)  

# lambdaTilde #
lambdaTilde <- array(dim=c(nx,ny,nm))
for (j in 1:nm) {
lambdaTilde[,,j] <- lambdaPos$v * densM$estimate[j]
}

nsim <- 19
S1PosMarkSimul <- rep(0,nsim)
S2PosMarkSimul <- rep(0,nsim)
S5PosMarkSimul <- rep(0,nsim)
S6PosMarkSimul <- rep(0,nsim)

system.time(for (i in 1:nsim) {

cat(paste(i, if (i < nsim) ", "
            else ".", if (i%%10 == 0) 
                "\n", sep = ""))

# Simulation of a point pattern of positions #
simulxy <- rpoispp(lambdaPos)

# Simulation of occurences #
simulm <- rpoint(simulxy$n,lambdam)$x

# Data #
data <- cbind(simulxy$x,simulxy$y,simulm)

# 3D-Density estimation #
dens3d <- kde(data,H=diag(hxym)^2,eval.points=s)
val <- dens3d$estimate
ok <- !inside.owin(sxyrep[,1],sxyrep[,2],Region)
val[ok] <- NA
dens3dArray <- aperm(array(val,c(ny,nx,nm)),c(2,1,3))

# lambdaHat #
lambdaHat <- Nech*dens3dArray

S1PosMarkSimul[i] <- max(abs(lambdaHat-lambdaTilde)/sqrt(lambdaTilde),na.rm=TRUE)
S2PosMarkSimul[i] <- min(abs(lambdaHat-lambdaTilde)/sqrt(lambdaTilde),na.rm=TRUE)
S5PosMarkSimul[i] <- mean((lambdaHat-lambdaTilde)^2,na.rm=TRUE)
S6PosMarkSimul[i] <- max((lambdaHat-lambdaTilde)^2,na.rm=TRUE)

})[3]


if ( S1PosMark < max(S1PosMarkSimul,na.rm=TRUE) ) { print("S1PosMark : Accept the separability assumption")
} else { print("S1PosMark : Reject the separability assumption") }
if ( S2PosMark < max(S2PosMarkSimul,na.rm=TRUE) ) { print("S2PosMark : Accept the separability assumption")
} else { print("S2PosMark : Reject the separability assumption") }
if ( S5PosMark < max(S5PosMarkSimul,na.rm=TRUE) ) { print("S5PosMark : Accept the separability assumption")
} else { print("S5PosMark : Reject the separability assumption") }
if ( S6PosMark < max(S6PosMarkSimul,na.rm=TRUE) ) { print("S6PosMark : Accept the separability assumption")
} else { print("S6PosMark : Reject the separability assumption") }



# 4.2. Dependence between marks and time #
##########################################

## Separability test ##

# Point pattern #
time <- ptime$marks
t <- (time-min(time))/(max(time)-min(time))
logworkload <- log(pworkload$marks)
m <- (logworkload-min(logworkload))/(max(logworkload)-min(logworkload))
pmarkstime <- ppp(x=t,y=m,c(0,1),c(0,1))

# Sub-sample #
Nech <- 2000   # 2000 : 372.64 secs
spmarkstime <- pmarkstime[sample(pmarkstime$n,Nech)]
data <- cbind(spmarkstime$x,spmarkstime$y)

# Regular grid of points #
matvarcov <- Hpi.diag(cbind(spmarkstime$x,spmarkstime$y),binned=TRUE)
lam <- density.ppp(spmarkstime,varcov=matvarcov,eps=0.01)
nt <- lam$dim[1]; nm <- lam$dim[2]
sxy <- gridcentres(spmarkstime$w,nm,nt)
sxyVec <- cbind(sxy$x,sxy$y) 

# Bandwidth #
ht <- sqrt(diag(matvarcov)[1])  
hm <- sqrt(diag(matvarcov)[2])  

# dens2D Array #
system.time(dens2d <- kde(data,H=diag(c(ht^2,hm^2)),eval.points=sxyVec))[3]
# lambda2d <- density.ppp(spmarkstime,varcov=matvarcov,eps=epsgrid,edge=FALSE)
val <- dens2d$estimate
dens2dArray <- aperm(array(val,c(nm,nt)),c(2,1))

# Time #
system.time(densT <- kde(data[,1],h=ht,eval.points=sxyVec[1:nt,1]))[3]

# Marks #
system.time(densM <- kde(data[,2],h=hm,eval.points=unique(sxyVec[,2])))[3]

# lambdaHat #
lambdaHat <- Nech*dens2dArray

# lambdaTilde #
valT <- rep(densT$estimate,nt)
valM <- rep(densM$estimate,nm)
val <- valT*valM
lambdaTilde <- Nech*array(val,c(nt,nm))

# LambdaHat and LambdaTilde #
lambdaHatImage <- lam
lambdaHatImage$v <- lambdaHat

lambdaTildeImage <- lam
lambdaTildeImage$v <- lambdaTilde

x11();
par(mfrow=c(1,3))

#postscript("Graphiques/lambdaHat.eps", horizontal=FALSE)
#titleplot <- "$\hat{\lambda}$"
plot(Transfo(lambdaHatImage,4),main=titleplot)
#dev.off()

#postscript("Graphiques/lambdaTilde.eps", horizontal=FALSE)
#titleplot <- "$\tilde{\lambda}$"
plot(Transfo(lambdaTildeImage,4),main=titleplot)
#dev.off()

# Separability Test #
S1MarkTime <- max(abs(lambdaHat-lambdaTilde)/sqrt(lambdaTilde),na.rm=TRUE)
S1MarkTime
S2MarkTime <- min(abs(lambdaHat-lambdaTilde)/sqrt(lambdaTilde),na.rm=TRUE)
S2MarkTime
S5MarkTime <- mean((lambdaHat-lambdaTilde)^2,na.rm=TRUE)
S5MarkTime
S6MarkTime <- max((lambdaHat-lambdaTilde)^2,na.rm=TRUE)
S6MarkTime


## Monte-Carlo Test with a Poisson point process model assumption ##
lambdaTildeImage <- lam
lambdaTildeImage$v <- lambdaTilde

nsim <- 19
S1MarkTimeSimul <- rep(0,nsim)
S2MarkTimeSimul <- rep(0,nsim)
S5MarkTimeSimul <- rep(0,nsim)
S6MarkTimeSimul <- rep(0,nsim)

for (i in 1:nsim) {

cat(paste(i, if (i < nsim) ", "
            else ".", if (i%%10 == 0) 
                "\n", sep = ""))

# Simulation of a point pattern #
simultm <- rpoispp(lambdaTildeImage)

# Data #
data <- cbind(simultm$x,simultm$y)

# dens2D Array #
system.time(dens2d <- kde(data,H=diag(c(ht^2,hm^2)),eval.points=sxyVec))[3]
# lambda2d <- density.ppp(spmarkstime,varcov=matvarcov,eps=epsgrid,edge=FALSE)
val <- dens2d$estimate
dens2dArray <- aperm(array(val,c(nm,nt)),c(2,1))

# lambdaHat #
lambdaHat <- Nech*dens2dArray

S1MarkTimeSimul[i] <- max(abs(lambdaHat-lambdaTilde)/sqrt(lambdaTilde),na.rm=TRUE)
S2MarkTimeSimul[i] <- min(abs(lambdaHat-lambdaTilde)/sqrt(lambdaTilde),na.rm=TRUE)
S5MarkTimeSimul[i] <- mean((lambdaHat-lambdaTilde)^2,na.rm=TRUE)
S6MarkTimeSimul[i] <- max((lambdaHat-lambdaTilde)^2,na.rm=TRUE)

}


if ( S1MarkTime < max(S1MarkTimeSimul,na.rm=TRUE) ) { print("S1MarkTime : Accept the separability assumption")
} else { print("S1MarkTime : Reject the separability assumption") }
if ( S2MarkTime < max(S2MarkTimeSimul,na.rm=TRUE) ) { print("S2MarkTime : Accept the separability assumption")
} else { print("S2MarkTime : Reject the separability assumption") }
if ( S5MarkTime < max(S5MarkTimeSimul,na.rm=TRUE) ) { print("S5MarkTime : Accept the separability assumption")
} else { print("S5MarkTime : Reject the separability assumption") }
if ( S6MarkTime < max(S6MarkTimeSimul,na.rm=TRUE) ) { print("S6MarkTime : Accept the separability assumption")
} else { print("S6MarkTime : Reject the separability assumption") }



# 4.3. Dependence between mark categories #
###########################################

## Kcross.inhom functions ##

# Data : June #
SampleMultitypepp <- ppp(x=pmonth$x,y=pmonth$y,marks=log(pmonth$marks),window=pmonth$window)
SampleMultitypepp <- cut(SampleMultitypepp,bornes,labels=c('Low','Medium','High'))

# Data : Subsample #
#Nech <- 4000
#SampleMultitypepp <- Multitypepp[sample(Multitypepp$n,Nech)]

# Split #
Low <- split(SampleMultitypepp)$Low
Medium <- split(SampleMultitypepp)$Medium
High <- split(SampleMultitypepp)$High

LowMedium <- SampleMultitypepp[SampleMultitypepp$marks=="Low" | SampleMultitypepp$marks=="Medium"]
levels(LowMedium) <- c("Low","Medium")
LowHigh <- SampleMultitypepp[SampleMultitypepp$marks=="Low" | SampleMultitypepp$marks=="High"]
levels(LowHigh) <- c("Low","High")
MediumHigh <- SampleMultitypepp[SampleMultitypepp$marks=="Medium" | SampleMultitypepp$marks=="High"]
levels(MediumHigh) <- c("Medium","High")

# Intensity estimation #
denspop <- density.ppp(Pop,sigma=900,weights=Pop$marks,eps=epsgrid)

# Semi-parametric #
XSamp <- unmark(SampleMultitypepp)
fit <- ppm(XSamp, ~ log(POP), covariates=list(POP=denspop))
intensSampleM <- denspop
intensSampleM$v <- exp (as.numeric(fit$coef[1]) + as.numeric(fit$coef[2]) * log(denspop$v))

intensLow <- denspop
intensLow$v <- (Low$n/SampleMultitypepp$n) * intensSampleM$v
intensMedium <- denspop
intensMedium$v <- (Medium$n/SampleMultitypepp$n) * intensSampleM$v
intensHigh <- denspop
intensHigh$v <- (High$n/SampleMultitypepp$n) * intensSampleM$v


# Kcross functions #
correc <- "border"
rmax <- min(diff(Region$xrange),diff(Region$yrange))/4
rvalue <- seq(0,rmax,length.out=512)

x11()
par(mfrow=c(1,3))
KcrossInhomEnv <- envelope(LowMedium, Kcross.inhom, nsim=39, i="Low", j="Medium", lambdaI=intensLow, lambdaJ=intensMedium, r=rvalue, correction=correc, simulate=expression(rlabel(LowMedium)))
#postscript("Graphiques/KLowMedium.eps", horizontal=FALSE)
#titleplot <- "Kcross Low/Medium"
plot(KcrossInhomEnv,cbind(sqrt(./pi)-r,0) ~ r,main=titleplot)
#dev.off()

KcrossInhomEnv <- envelope(LowHigh, Kcross.inhom, nsim=39, i="Low", j="High", lambdaI=intensLow, lambdaJ=intensHigh, r=rvalue, correction=correc, simulate=expression(rlabel(LowHigh)))
#postscript("Graphiques/KLowHigh.eps", horizontal=FALSE)
#titleplot <- "Kcross Medium/High"
plot(KcrossInhomEnv,cbind(sqrt(./pi)-r,0) ~ r,main=titleplot)
#dev.off()

KcrossInhomEnv <- envelope(MediumHigh, Kcross.inhom, nsim=39, i="Medium", j="High", lambdaI=intensMedium, lambdaJ=intensHigh, r=rvalue, correction=correc, simulate=expression(rlabel(MediumHigh)))
#postscript("Graphiques/KMediumHigh.eps", horizontal=FALSE)
#titleplot <- "Kcross Medium/High"
plot(KcrossInhomEnv,cbind(sqrt(./pi)-r,0) ~ r,main=titleplot)
#dev.off()



# 4.4. Marginal distribution #
##############################

logworkload <- log(pworkload$marks)

# Histogram #
x11()
#postscript("Graphiques/Histo.eps", horizontal=FALSE)
#titleplot <- "Histogram of the logarithm of the workload marks"
hist(logworkload,main=titleplot)
#dev.off()

# Q-Q plot #
y <- logworkload
x11()
#postscript("Graphiques/QQplot.eps", horizontal=FALSE)
#titleplot <- "QQplot of the logarithm of the workload marks"
qqnorm(y)
qqline(y)
#dev.off()

# Normality test #
# Shapiro.test : number of values between 3 and 5000 #
ysub <- sample(y,5000)
shapiro.test(ysub)

# Box-Cox #
ml <- boxcox.fit(pworkload$marks)
y <- box.cox(pworkload$marks,ml$lambda)    # ml$lambda = 0.126
x11()
hist(y)
x11()
qqnorm(y)
qqline(y)


## Exponential ##
lambda <- 1/mean(pworkload$marks)
x11()
hist(rexp(pworkload$n,rate=lambda))

u <- seq((1/pworkload$n),(1-(1/pworkload$n)),len=pworkload$n)
qe <- qexp(u, rate = lambda)
qy <- sort(pworkload$marks)

x11()
plot(qe,qy)



## Pareto ##
execPareto <- function(){ 

dpareto <- function(x,a=0.5,b=1) {
a*b^a/x^(a+1)
}

ppareto <- function(x,a=0.5,b=1) {
(x > b)*(1-(b/x)^a)
}

qpareto <- function(u,a=0.5,b=1) {
b/(1-u)^(1/a)
}

rpareto <- function(n,a=0.5,b=1) {
qpareto(runif(n),a,b)
}

# When b=1 #

# Five estimators : http://astrostatistics.psu.edu/datasets/2006tutorial/2006mle.html #
five <- function(d) {
     lsd <- log(sort(d))
     n <- length(d)
     lseq <- log((n:1)/n)
     m1=lm(lseq ~ lsd)$coef
     hd <- hist(d,nclass=n/5,plot=F)
     counts <- hd$counts
     ldens <- log(hd$density[counts>0])
     lmidpts <- log(hd$mids[counts>0])
     m2 <- lm(ldens~lmidpts)$coef
     m3 <- lm(ldens~lmidpts, weights=counts[counts>0])$coef
     out <- c(max.lik=1/mean(log(d)), meth.mom=1/(1-1/mean(d)),
          EDF=-as.numeric(m1[2]), unwt.hist=-1-as.numeric(m2[2]),
          wt.hist=-1-as.numeric(m3[2]))
     return(out)
   }


estim.a <- five(pworkload$marks)

for (i in 1:5) {
 
a <- estim.a[i]

u <- seq((1/pworkload$n),(1-(1/pworkload$n)),len=pworkload$n)
qp <- qpareto(u,a=a,b=1)
qy <- sort(pworkload$marks)

x11()
plot(qp,qy)

x11()
hist(log(rpareto(pworkload$n,a=a,b=1)))
}

# When b = 2 #
b <- 2
max.lik.a <- pworkload$n/(sum(log(pworkload$marks)) - b*pworkload$n)
meth.mom.a <- 1/(1-b/mean(pworkload$marks))
# 1 < meth.mom.a < 2 : Pb with the estimation of b by the method of moments # 
estim.a <- c(max.lik.a,meth.mom.a)

for (i in 1:2) {

a <- estim.a[i]

u <- seq((1/pworkload$n),(1-(1/pworkload$n)),len=pworkload$n)
qp <- qpareto(u,a=a,b=2)
qy <- sort(pworkload$marks)

x11()
plot(qp,qy)

x11()
hist(log(rpareto(pworkload$n,a=a,b=2)))
} 

# Sequence of b and moment estimator for a #
b <- seq(0.1,2,0.3)
estim.a <- 1/(1-b/mean(pworkload$marks))

for (i in 1:length(estim.a)) {

a <- estim.a[i]

u <- seq((1/pworkload$n),(1-(1/pworkload$n)),len=pworkload$n)
qp <- qpareto(u,a=a,b=2)
qy <- sort(pworkload$marks)

x11()
plot(qp,qy)

x11()
hist(log(rpareto(pworkload$n,a=a,b=b[i])))
} 
}

#execPareto()


## Gev ##

# with logworkload #
Mgev <- fgev(logworkload)
simulgev <- rgev(pworkload$n, loc=Mgev$estimate[1], scale=Mgev$estimate[2], shape=Mgev$estimate[3])
x11()
hist(simulgev,nclass=15)

u <- seq((1/pworkload$n),(1-(1/pworkload$n)),len=pworkload$n)
qgev <- qgev(u,  loc=Mgev$estimate[1], scale=Mgev$estimate[2], shape=Mgev$estimate[3])
qy <- sort(logworkload)

x11()
plot(qgev,qy)


# with pworkload$marks #
Mgev <- fgev(pworkload$marks)
simulgev <- rgev(pworkload$n, loc=Mgev$estimate[1], scale=Mgev$estimate[2], shape=Mgev$estimate[3])
x11()
hist(simulgev,nclass=15)

u <- seq((1/pworkload$n),(1-(1/pworkload$n)),len=pworkload$n)
qgev <- qgev(u,  loc=Mgev$estimate[1], scale=Mgev$estimate[2], shape=Mgev$estimate[3])
qy <- sort(pworkload$marks)

x11()
plot(qgev,qy)




###################
# 5. Modelisation #
###################

pmonth <- p06
correc <- "border"
rmax <- min(diff(Region$xrange),diff(Region$yrange))/4
rvalue <- seq(0,rmax,length.out=512)

# Point pattern #
X <- unmark(pmonth)



# 5.1. Parametric and Nonparametric estimation #
################################################

## Parametric ##

# Polynomial trend #
fit <- ppm(X, ~ polynom(x,y,4))
# degree > 4 -> Warning message

lam <- predict(fit,type="lambda",ngrid=intens$dim)

x11();
par(mfrow=c(1,3))

# Background intensity #
#postscript("Graphiques/ParamLamKen1.eps", horizontal=FALSE)
#titleplot <- "Logarithm transformation of a parametric \n estimate of intensity for emergencies in June"
plot(Transfo(lam,0),main=titleplot)
#dev.off()

# Simulation #
simul <- rpoispp(lam)
simul$window <- Region
#postscript("Graphiques/ParamLamKen2.eps", horizontal=FALSE)
#titleplot <- "Simulation of a Poisson point process \n with the estimated intensity"
plot(simul,main=titleplot)
#dev.off()

# Envelope #
Ken <- envelope(X, Kinhom, nsim=39, simulate=expression(rpoispp(lam)), lambda=lam, r=rvalue, correction=correc)
#postscript("Graphiques/ParamLamKen3.eps", horizontal=FALSE)
#titleplot <- "Test of no interaction with L(r)-r \n Parametric"
plot(Ken, cbind(sqrt(./pi)-r,0) ~ r,main=titleplot)
#dev.off()



## Nonparametric estimation ##

## "sigma" choose by the user ##
lam <- density.ppp(X,sigma=h,eps=epsgrid)
x11();
par(mfrow=c(1,3))

# Background intensity #
#postscript("Graphiques/NonParamLamKenh1.eps", horizontal=FALSE)
#titleplot <- "Logarithm transformation of a nonparametric \n estimate of intensity with \n h=800 for emergencies in June"
plot(Transfo(lam,0),main=titleplot)
#dev.off()

# Simulation #
simul <- rpoispp(lam)
simul$window <- Region
#postscript("Graphiques/NonParamLamKenh2.eps", horizontal=FALSE)
#titleplot <- "Simulation of a Poisson point process \n with the estimated intensity"
plot(simul,main=titleplot)
#dev.off()

# Envelope #
Ken <- envelope(X, Kinhom, nsim=39, simulate=expression(rpoispp(lam)), lambda=lam, r=rvalue, correction=correc)
#postscript("Graphiques/NonParamLamKenh3.eps", horizontal=FALSE)
#titleplot <- "Test of no interaction with L(r)-r \n Nonparametric, h=800, lambda Hat"
plot(Ken, cbind(sqrt(./pi)-r,0) ~ r,main=titleplot)
#dev.off()



## case of over-fitting => use the lambda tilde in Baddeley(2000) ##
Ken <- envelope(X, Kinhom, nsim=39, nrank=1, simulate=expression(rpoispp(lam)), sigma=800, r=rvalue, correction=correc)
x11();
#postscript("Graphiques/NonParamLOVKenh.eps", horizontal=FALSE)
#titleplot <- "Test of no interaction with L(r)-r \n Nonparametric, h=800, lambda Tilde"
plot(Ken, cbind(sqrt(./pi)-r,0) ~ r,main=titleplot)
#dev.off()
# Pb with Leave-one-out estimation : nonconventional form #



## With intensity estimated on the emergencies of an other month ##
pmonthbef <- p05   # p05 : relatively "good" result -- p07 : "artefacts"
Xbef <- unmark(pmonthbef)

lam <- density.ppp(Xbef,sigma=h,eps=epsgrid)
lam$v <- (X$n/Xbef$n)*lam$v

x11();
par(mfrow=c(1,3))

# Intensity #
#postscript("Graphiques/NonParamMonthBefLamKenh1.eps", horizontal=FALSE)
#titleplot <- "Logarithm transformation of a nonparametric \n estimate of intensity with \n h=800 for emergencies in May"
plot(Transfo(lam,0),main=titleplot)
#dev.off()

# Simulation #
simul <- rpoispp(lam)
simul$window <- Region
#postscript("Graphiques/NonParamMonthBefLamKenh2.eps", horizontal=FALSE)
#titleplot <- "Simulation of a Poisson point process \n with the estimated intensity"
plot(simul,main=titleplot)
#dev.off()

# Envelope #
Ken <- envelope(X, Kinhom, nsim=39, simulate=expression(rpoispp(lam)), lambda=lam, r=rvalue, correction=correc)
#postscript("Graphiques/NonParamMonthBefLamKenh3.eps", horizontal=FALSE)
#titleplot <- "Test of no interaction with L(r)-r \n Nonparametric, h=800, lambda Hat"
plot(Ken, cbind(sqrt(./pi)-r,0) ~ r,main=titleplot)
#dev.off()



# 5.2. Semiparametric estimation with one covariable #
######################################################

# 5.2.1. Two density estimates #
################################

x11();
par(mfrow=c(1,2))

# IRIS population #
#postscript("Graphiques/Pop1.eps", horizontal=FALSE)
#titleplot <- "Population represented by disc centered on IRIS with radii proportional to the number of inhabitants"
plot(Pop,main=titleplot)
#dev.off()


## Density estimation of the population ##
 
# Global bandwith #
intenspop1 <- density.ppp(Pop,sigma=900,weights=Pop$marks,eps=epsgrid)
denspop1 <- intenspop1
denspop1$v <- denspop1$v/(sum(Pop$marks))


# Local bandwidth : adaptive kernel estimation by k-nearest neighbours #
XPop <- unmark(Pop)
si <- cbind(XPop$x,XPop$y)
wi <- Pop$marks
nx <- denspop1$dim[1]
ny <- denspop1$dim[2]
w <- denspop1$window

s <- gridcentres(Region,ny,nx)
s <- cbind(s$x,s$y)

# Number of neighbours #
k <- 5

spp <- ppp(s[,1],s[,2],window=Region)
d <- crossdist(spp,XPop)
dsort <- t(apply(d,1,sort))
vectlocalh <- dsort[,k] 
ok <- !inside.owin(s[,1],s[,2],Region)
s[ok] <- NA
vectlocalh[ok] <- NA

# Density estimation by an Epanechnikov kernel #

DensEpanech <- function(s,si,h,wi) {
h <- ceiling(h)
m <- dim(s)[1]
n <- dim(si)[1]
N <- sum(wi)

Xs <- s[,1]
matXs <- matrix(Xs,m,n)
Xsi <- si[,1]
matXsi <- matrix(Xsi,m,n,byrow=TRUE)
Ys <- s[,2]
matYs <- matrix(Ys,m,n)
Ysi <- si[,2]
matYsi <- matrix(Ysi,m,n,byrow=TRUE)
matW <- matrix(wi,m,n,byrow=TRUE)
norm2 <- (1/h^2)*( (matXs-matXsi)^2 + (matYs-matYsi)^2 )
dummy <- norm2 < 1
sumnorm2 <- rowSums(matW*norm2*dummy)

val <- (1/n)*(3/4)*(1/h^2)*(1 - (1/N)*sumnorm2)
matval <- matrix(val,nx,ny,byrow=TRUE)

}

system.time(Dens <- DensEpanech(s,si,vectlocalh,wi))[3]
denspop2 <- denspop1
denspop2$v <- Dens


# Compare 2 estimations #

# Global bandwidth #
#postscript("Graphiques/Pop2.eps", horizontal=FALSE)
#titleplot <- "Logarithm transformation \n of the population density estimated \n by a nonparametric kernel method"
plot(Transfo(denspop1,0),main=titleplot)
#dev.off()

# Local bandwidth #
#postscript("Graphiques/Pop3.eps", horizontal=FALSE)
#titleplot <- "Logarithm transformation \n of the population density estimated \n by an adaptive nonparametric kernel method"
plot(Transfo(denspop2,0),main=titleplot)
#dev.off()



# 5.2.2. and 5.2.3. Models #
###########################

## Choice of the intensity estimation ##
#denspop <- denspop1
denspop <- denspop2

# Fitting #
fit <- ppm(X, ~ log(POP), covariates=list(POP=denspop))
lam <- denspop
lam$v <- exp (as.numeric(fit$coef[1]) + as.numeric(fit$coef[2]) * log(denspop$v))



# Poisson #
###########

# Envelope #
x11();
par(mfrow=c(1,2))

# Simulation #
simul <- rpoispp(lam)
simul$window <- Region
#postscript("Graphiques/PoissonK1.eps", horizontal=FALSE)
#titleplot <- "Simulation of a Poisson point process \n with the estimated intensity"
plot(simul,main=titleplot)
#dev.off()

# Envelope #
Ken <- envelope(X, Kinhom, nsim=39, simulate=expression(rpoispp(lam)), lambda=lam , r=rvalue, correction=correc)
#postscript("Graphiques/PoissonK2.eps", horizontal=FALSE)
#titleplot <- "Test of no interaction with L(r)-r \n Semiparametric"
plot(Ken, cbind(sqrt(./pi)-r,0) ~ r,main=titleplot)
#dev.off()



# Inhomogeneous Thomas Cluster by thinning #
############################################

Ki <- Kinhom(X, lambda=lam, r=rvalue, correction=correc)
rmin <- 0
rmax <- 4000

v <- thomas.estK(Ki, startpar=c(kappa=1e-4,sigma2=200),rmin=rmin, rmax=rmax)
v

x11();
#postscript("Graphiques/MinContrastThomas2.eps", horizontal=FALSE)
#titleplot <- "Minimum contrast estimation \n in the case of a Thomas Cluster process"
plot(v,sqrt(./pi)-r ~ r,main=titleplot)
#dev.off()

# Simulation Inhomogeneous Thomas Cluster #
kappa <- v$modelpar[1]
sigma <- v$modelpar[2]
expBeta0 <- exp( as.numeric(fit$coef[1]) )
alpha <- expBeta0 / kappa
expCov <- lam
expCov$v <- lam$v/expBeta0
M <- max(expCov$v,na.rm=TRUE)
Malpha <- M*alpha

simulY <- rThomas(kappa, sigma, mu=Malpha, win = Region)
simulY <- no.na(simulY,expCov)

RetentProb <- expCov
RetentProb$v <- expCov$v/M

simulX <- rthin(simulY,RetentProb)


x11();
par(mfrow=c(1,2))

# Simulation #
#postscript("Graphiques/ThomasK1.eps", horizontal=FALSE)
#titleplot <- "Simulation of an inhomogeneous \n Thomas cluster point process \n obtained by thinning"
plot(simulX,main=titleplot)
#dev.off()

# Envelope #
Ken <- envelope(X, Kinhom, nsim=39, simulate=expression(rthin(no.na(rThomas(kappa, sigma, mu=Malpha, win = Region),expCov),RetentProb)), lambda=lam , r=rvalue, correction=correc)
#postscript("Graphiques/ThomasK2.eps", horizontal=FALSE)
#titleplot <- "Test of goodness-of-fit \n for an inhomogeneous Thomas cluster point process model"
plot(Ken, cbind(sqrt(./pi)-r,0) ~ r,main=titleplot)
#dev.off()



# Inhomogeneous Matern Cluster by thinning #
############################################

Ki <- Kinhom(X, lambda=lam, r=rvalue, correction=correc)
rmin <- 0
rmax <- 4000

v <- matclust.estK(Ki, startpar=c(kappa=1e-4,R=200), rmin=rmin, rmax=rmax)
v

x11();
#postscript("Graphiques/MinContrastMatern.eps", horizontal=FALSE)
#titleplot <- "Minimum contrast estimation \n in the case of a Matern Cluster process"
plot(v,sqrt(./pi)-r ~ r,main=titleplot)
#dev.off()

# Simulation #
kappa <- v$modelpar[1]
r <- v$modelpar[2]
expBeta0 <- exp( as.numeric(fit$coef[1]) )
alpha <- expBeta0 / kappa
expCov <- lam
expCov$v <- lam$v/expBeta0
M <- max(expCov$v,na.rm=TRUE)
Malpha <- M*alpha

simulY <- rMatClust(kappa, r, mu=Malpha, win = Region)
simulY <- no.na(simulY,expCov)

RetentProb <- expCov
RetentProb$v <- expCov$v/M

simulX <- rthin(simulY,RetentProb)

x11();
par(mfrow=c(1,2))

# Simulation #
#postscript("Graphiques/MaternK1.eps", horizontal=FALSE)
#titleplot <- "Simulation of an inhomogeneous \n Matern cluster point process \n obtained by thinning"
plot(simulX,main=titleplot)
#dev.off()

# Envelope #
Ken <- envelope(X, Kinhom, nsim=39, simulate=expression(rthin(no.na(rMatClust(kappa, r, mu=Malpha, win = Region),expCov),RetentProb)), lambda=lam , r=rvalue, correction=correc)
#titleplot <- "Test of goodness-of-fit \n for an inhomogeneous Matern cluster point process model"
plot(Ken, cbind(sqrt(./pi)-r,0) ~ r,main=titleplot) 
#dev.off()



# Inhomogeneous LGCP by thinning #
##################################

# with RandomFields Package #

Ki <- Kinhom(X, lambda=lam, r=rvalue, correction=correc)
rmin <- 0
rmax <- 4000

v <- lgcp.estK(Ki, c(sigma2=1, alpha=300), rmin=rmin, rmax=rmax)
v

x11();
#postscript("Graphiques/MinContrastLGCP.eps", horizontal=FALSE)
#titleplot <- "Minimum contrast estimation \n in the case of a LGCP"
plot(v,sqrt(./pi)-r ~ r,main=titleplot)
#dev.off()

# Parameters #
variance <- v$modelpar[1]
nugget <- 0
scale <- v$modelpar[2]

BetaTilde0 <- as.numeric(fit$coef[1])
BetaHat0 <- BetaTilde0 - variance/2
lambda <- lam
lambda$v <- BetaHat0 + as.numeric(fit$coef[2]) * log(denspop$v)
M <- max(lambda$v,na.rm=TRUE)

mean <- M

(p <- c(mean,variance,nugget,scale))

RetentProb <- lambda
RetentProb$v <- exp(lambda$v-M)

xstep <- ystep <- epsgrid
nxy <- denspop$dim
nx <- nxy[1]
ny <- nxy[2]
xr <- Region$xrange
yr <- Region$yrange
x <- seq(xr[1], xr[2], length = 2 * nx + 1)[2 * (1:nx)]
y <- seq(yr[1], yr[2], length = 2 * ny + 1)[2 * (1:ny)]
x <- rep(x, ny)
y <- rep(y, rep(nx, ny))


# Espectation of the number of points for Y #
EspNbPointsY <- as.numeric(exp(M+variance/2)*area.owin(Region))
EspNbPointsY


LGCPmodel <- function(p,xstep,ystep,RetentProb,Region,nx,ny) {

f <- GaussRF(x=c(Region$yrange[1],Region$yrange[2],ystep),y=c(Region$xrange[1],Region$xrange[2],xstep),grid=TRUE,model="exponen", param=p, gridtriple=TRUE)

Y <- denspop
Y$v <- f

v <- gridcenters(Region,nx,ny)
ok <- !inside.owin(v$x,v$y,Region)

aux <- Y
aux$v[ok] <- NA
Ydomain <- RetentProb
Ydomain$v <- exp(aux$v)

# for avoiding Warning message : low memory #
lmax <- max(Ydomain$v,na.rm=TRUE)
limitNb <- 1e6
limitInt <- limitNb/area.owin(Region)

if (lmax > limitInt) {
numPP <- trunc(lmax/limitInt)
Nbrestant <- round(lmax*area.owin(Region) - numPP*limitNb)

Ydomain2 <- Ydomain
Ydomain2$v <- (limitInt/lmax)*Ydomain$v
Ydomain3 <- Ydomain
Ydomain3$v <- ((Nbrestant/area.owin(Region))/lmax)*Ydomain$v

LGCP <- rpoispp(Ydomain2)

for (i in 1:numPP) {
LGCPaux <- rpoispp(Ydomain2)
LGCP <- superimpose(LGCP,LGCPaux)
}

LGCPaux <- rpoispp(Ydomain3)
LGCP <- superimpose(LGCP,LGCPaux)

} else { LGCP <- rpoispp(Ydomain) }


LGCP$window <- Region

# Thinning #
simulX <- rthin(no.na(LGCP,RetentProb),RetentProb)
simulX

}

simulX <- LGCPmodel(p,xstep,ystep,RetentProb,Region,nx,ny)


x11();
par(mfrow=c(1,2))

# Simulation #
#postscript("Graphiques/LGCPthinningK1.eps", horizontal=FALSE)
#titleplot <- "Simulation of an inhomogeneous LGCP \n obtained by thinning"
plot(simulX,main=titleplot)
#dev.off()

# Envelope #
Ken <- envelope(X, Kinhom, nsim=39, simulate=expression(LGCPmodel(p,xstep,ystep,RetentProb,Region,nx,ny)), lambda=lam , r=rvalue, correction=correc)
#postscript("Graphiques/LGCPthinningK2.eps", horizontal=FALSE)
#titleplot <- "Test of goodness-of-fit \n for an inhomogeneous LGCP model"
plot(Ken, cbind(sqrt(./pi)-r,0) ~ r,main=titleplot)
#dev.off()
