##### This is for fitting a Hawkes model,  
    ##### on the space S = [0,X1] x [0,Y1] (km), in time [0,T]. 
    ##### lambda(t,x,y) = mu rho(x,y) + K SUM gt(t-t_i)gxy(x-xi,y-yi),  
    ##### with rho(x,y) = 1/(X1Y1), 
    ##### gt(t) = beta e^(-beta t),
    ##### gxy(x,y) = alpha/pi exp(-alpha r^2), with x^2+y^2=r^2. 
    ##### For any theta, the integral of lambda over the space time region is approx. 
    ##### mu T + Kn.
    ##### The parameter vector theta = (mu, K, alpha, beta)

m3 = function(x) signif(x,3)

## If you have: 
## datax 
## datay
## datat 
## datan, 
## then write  
## z = list()
## z$t = datat
## z$lon = datax 
## z$lat = datay 
## z$n = datan 
## library(ggplot2) 
## help(stat_density_2d)
## example(geom_density_2d) 

## Make sure the data are stored in z, and you define T,X1, and Y1 externally. 
## First we will write the loglikelihood function in R. 
loglhawk = function(theta,draw=0){
mu = theta[1]; K = theta[2]; alpha = theta[3]; beta = theta[4] 
cat("\n mu = ",m3(mu),", K = ",m3(K),", alpha = ",m3(alpha),", beta = ",m3(beta),".\n") 
if(min(mu,K,alpha,beta)<0.000000001) return(99999) 
if(K>.99999) return(99999)
if(draw){
r = seq(0,3,length=100)
t = alpha/pi * exp(-alpha * r^2)
lines(r,t,col="orange",lty=2) 
}
sumlog = log(mu/X1/Y1) 
intlam = mu*T + K*z$n
const = K*alpha/pi*beta
for(j in 2:(z$n)){
   gij = 0
   for(i in 1:(j-1)){
    r2 = (z$lon[j]-z$lon[i])^2+(z$lat[j]-z$lat[i])^2
    gij = gij + exp(-beta*(z$t[j]-z$t[i])-alpha*r2)
    }
   lamj = mu / X1 / Y1 + const*gij
   if(lamj < 0){
    cat("lambda ",j," is less than 0.")
    return(99999)
   }
sumlog = sumlog + log(lamj)
}
loglik = sumlog - intlam
cat("loglike is ", loglik, ". sumlog = ", sumlog,". integral = ", intlam,".\n")
if(draw) lines(r,t,col="white",lty=2) 
return(-1.0*loglik)
}