if(dev.cur() <= 1) {
  dd <- getOption("device")
  if(is.character(dd)) dd <- get(dd)
  dd()
}

oldpar <- par(ask = interactive() &&
              (.Device %in% c("X11", "GTK", "windows", "Macintosh")))
par(mfrow=c(1,1))
oldoptions <- options(warn = -1)

# 
#######################################################
#

X <- rpoispp(function(x,y) { 1000 * exp(- 4 * x)}, 1000)
plot(X, main="Inhomogeneous Poisson pattern")

fit.hom <- ppm(X, ~1, Poisson())
fit.inhom <- ppm(X, ~x, Poisson())

diagnose.ppm(fit.inhom, which="marks", type="Pearson",
             main="Mark plot\nCircles for positive residual mass\nColour for negative residual density")

par(mfrow=c(1,2))
diagnose.ppm(fit.hom, which="marks",
             main=c("Wrong model (homogeneous Poisson)", "raw residuals"))
diagnose.ppm(fit.inhom, which="marks",
             main=c("Right model (inhomogeneous Poisson)", "raw residuals"))

par(mfrow=c(1,1))
diagnose.ppm(fit.inhom, which="smooth", main="Smoothed residual field")

par(mfrow=c(1,2))
diagnose.ppm(fit.hom, which="smooth",
             main=c("Wrong model (homogeneous Poisson)",
                    "Smoothed residual field"))
diagnose.ppm(fit.inhom, which="smooth",
             main=c("Right model (inhomogeneous Poisson)",
                    "Smoothed residual field"))

par(mfrow=c(1,1))
diagnose.ppm(fit.inhom, which="x")

par(mfrow=c(1,2))
diagnose.ppm(fit.hom, which="x",
             main=c("Wrong model (homogeneous Poisson)",
                    "lurking variable plot for x coordinate"))
diagnose.ppm(fit.inhom, which="x",
             main=c("Right model (inhomogeneous Poisson)",
                    "lurking variable plot for x coordinate"))

par(mfrow=c(1,1))
diagnose.ppm(fit.hom, type="Pearson",main="standard diagnostic plots")

par(mfrow=c(1,2))
diagnose.ppm(fit.hom, main="Wrong model (homogeneous Poisson)")
diagnose.ppm(fit.inhom,  main="Right model (inhomogeneous Poisson)")
par(mfrow=c(1,1))

# 
#######################################################
#  Q - Q  PLOTS
#
qqplot.ppm(fit.hom, 40) 
#conclusion: homogeneous Poisson model is not correct
title(main="Q-Q plot of smoothed residuals")

qqplot.ppm(fit.inhom, 40) # TAKES A WHILE...
title(main=c("Right model (inhomogeneous Poisson)",
             "Q-Q plot of smoothed residuals"))
# conclusion: fitted inhomogeneous Poisson model looks OK
# 
#######################################################
#
data(cells)
plot(cells)
fitPoisson <- ppm(cells, ~1, Poisson())
diagnose.ppm(fitPoisson,
             main=c("CSR fitted to cells data",
                    "Raw residuals",
                    "No suggestion of departure from CSR"))
diagnose.ppm(fitPoisson, type="pearson",
             main=c("CSR fitted to cells data",
                    "Pearson residuals",
                    "No suggestion of departure from CSR"))
# These diagnostic plots do NOT show evidence of departure from uniform Poisson

qqplot.ppm(fitPoisson, 40)
title(main=c("CSR fitted to cells data",
        "Q-Q plot of smoothed raw residuals",
        "Strong suggestion of departure from CSR"))
           
# Q-Q plot DOES show strong evidence of departure from uniform Poisson.
#
fitStrauss <- ppm(cells, ~1, Strauss(r=0.15))
diagnose.ppm(fitStrauss,
             main=c("Strauss model fitted to cells data",
                    "Raw residuals"))
diagnose.ppm(fitStrauss, type="pearson",
             main=c("Strauss model fitted to cells data",
                    "Pearson residuals"))
# next line takes a LOOONG time ...
qqplot.ppm(fitStrauss, 40, type="pearson")
title(main=c("Strauss model fitted to cells data",
        "Q-Q plot of smoothed Pearson residuals",
        "Suggests adequate fit")) 
# Conclusion: Strauss model seems OK
# 
#######################################################
#
data(nztrees)
plot(nztrees)
fit <- ppm(nztrees, ~1, Poisson())
diagnose.ppm(fit, type="pearson")
title(main=c("CSR fitted to NZ trees",
             "Pearson residuals"))
diagnose.ppm(fit, type="pearson", cumulative=FALSE)
title(main=c("CSR fitted to NZ trees",
             "Pearson residuals (non-cumulative)"))
lurking(fit, expression(x), type="pearson", cumulative=FALSE,
        splineargs=list(spar=0.3))
# Sharp peak at right is suspicious
qqplot.ppm(fit, 40, type="pearson")
title(main=c("CSR fitted to NZ trees",
        "Q-Q plot of smoothed Pearson residuals"))
# Slight suggestion of departure from Poisson at top right of pattern.
par(oldpar)
options(oldoptions)