.packageName <- "InhomCluster"
#matrices J and G
J=function(z,beta,rule){
# z is N x (p+1) where N is the number of grid points
  p=dim(z)[2]-1
  ztz=t(apply(z,1,ztz.row))
  if (!is.matrix(ztz))
    ztz=matrix(ztz,ncol=1)
  
  linpred=z%*%matrix(beta,ncol=1)
  w=rule$w*exp(linpred) #rule: x y w hvor x y skal svare til punkter hvor z er evalueret 
  integrand=sweep(ztz,MARGIN=1,STATS=w,FUN="*")
  
  integral=apply(integrand,2,sum)
  
  return(matrix(integral,ncol=p+1,byrow=T))
}

ztz.row=function(z){
  z=matrix(z,ncol=1)
  return(c(t(z%*%t(z))))
}

G=function(z,beta,srule,grule,neighbhd,omega){ # x y koordinater hvor z er observeret
  p=dim(z)[2]-1
  xy=cbind(srule$x,srule$y)
  cdim=length(grule$x) # grule: x y w
  inner=matrix(0,cdim[1],dim(z)[2])
  linpred=z%*%matrix(beta,ncol=1)
  
  dimx=length(unique(xy[,1]))
  dimy=length(unique(xy[,2]))
  discrx=xy[dimy+1,1]-xy[1,1]
  discry=xy[2,2]-xy[1,2]

  inner=.C("Ginner",as.double(rep(0,cdim*dim(z)[2])),as.double(t(xy)),as.double(t(z)),as.double(grule$x),as.double(grule$y),as.double(srule$w),as.double(linpred),as.integer(dim(z)[2]),as.integer(cdim),as.integer(dimx),as.integer(dimy),as.double(discrx),as.double(discry),as.double(neighbhd),as.double(omega))[[1]]

  inner=matrix(inner,ncol=p+1,byrow=T)
  inner=t(apply(inner,1,ztz.row))
  inner=apply(sweep(inner,MARGIN=1,STATS=grule$w,FUN="*" ),2,sum)
  return(matrix(inner,ncol=p+1))
}

inhom.thomas.asympcov=function(z,beta,srule,grule,neighbhd,kappa,omega){
  
  lrx=min(srule$x)
  lry=min(srule$y)
  srule$x = srule$x - lrx
  srule$y = srule$y - lry
  grule$x = grule$x - lrx
  grule$y = grule$y - lry

  Jmat=J(z,beta,srule)
  Jinv=solve(Jmat)
  Gmat=G(z,beta,srule,grule,neighbhd,omega)
  asycov=Jinv+Jinv%*%Gmat%*%Jinv/kappa

  return(list(asycov=asycov,asycovpois=Jinv))
}



myKest<-function(data,lambda,maxt){
  
  if (data$window$type!="rectangle"){
    print("stop only rectangular observation windows allowed")
    stop
  }

  data$x=data$x-data$window$xrange[1]
  data$y=data$y-data$window$yrange[1]
  data$window$xrange=c(0,data$window$xrange[2]-data$window$xrange[1])
  data$window$yrange=c(0,data$window$yrange[2]-data$window$yrange[1])


  x=as.double(data$x)
  y=as.double(data$y)
  
  sideEW=as.double(data$window$xrange[2])
  sideNS=as.double(data$window$yrange[2])
  
  par=as.double(lambda)
  antpktdat=as.integer(length(x))
  maxt=as.double(maxt)
  Kestimate=as.double(rep(0,101))

  parametric=as.integer(F)
  correctype=as.integer(3)

  #print(c(length(x),length(y),sideEW,sideNS,length(par),antpktdat,maxt,parametric,correctype))
  Kestimate<-.C("Ktrans",antpktdat,x,y,par,maxt,Kestimate,correctype,sideEW,sideNS,parametric)[[6]]
  return(list(t=c(0:100)*maxt/100,K=Kestimate[1:101]))
}

myKest.cross<-function(datai,dataj,lambdai,lambdaj,maxt){


  if (datai$window$type!="rectangle" & dataj$window$type!="rectangle"){
    print("stop only rectangular observation windows allowed")
    stop
  }

  datai$x=datai$x-datai$window$xrange[1]
  datai$y=datai$y-datai$window$yrange[1]
  datai$window$xrange=c(0,datai$window$xrange[2]-datai$window$xrange[1])
  datai$window$yrange=c(0,datai$window$yrange[2]-datai$window$yrange[1])
  
  xi=as.double(datai$x)
  yi=as.double(datai$y)
  pari=as.double(lambdai)
  antpktdati=as.integer(length(xi))

  if (datai$window$type!="rectangle"){
    print("stop only rectangular observation windows allowed")
    stop
  }

  dataj$x=dataj$x-dataj$window$xrange[1]
  dataj$y=dataj$y-dataj$window$yrange[1]
  dataj$window$xrange=c(0,dataj$window$xrange[2]-dataj$window$xrange[1])
  dataj$window$yrange=c(0,dataj$window$yrange[2]-dataj$window$yrange[1])
  
  xj=as.double(dataj$x)
  yj=as.double(dataj$y)
  parj=as.double(lambdaj)
  antpktdatj=as.integer(length(xj))  

  maxt=as.double(maxt)
  
  sideEW<-as.double(datai$window$xrange[2])
  sideNS<-as.double(datai$window$yrange[2])

  if (sideEW!=dataj$window$xrange[2] | sideNS!=dataj$window$yrange[2]){
    print("different windows for two types of points")
    stop
  }
  
  parametric=as.integer(F)
  correctype=as.integer(3)
  
  Kestimate1=as.double(rep(0,101))
  Kestimate1<-.C("Ktrans_cross",antpktdati,xi,yi,antpktdatj,xj,yj,pari,parj,maxt,Kestimate1,correctype,sideEW,sideNS,parametric)[[10]]
  Kestimate2=as.double(rep(0,101))
  Kestimate2<-.C("Ktrans_cross",antpktdatj,xj,yj,antpktdati,xi,yi,parj,pari,maxt,Kestimate2,correctype,sideEW,sideNS,parametric)[[10]]
  return(list(t=c(0:100)*maxt/100,Kij=Kestimate1[1:101],Kji=Kestimate2[1:101]))
}

mygest<-function(data,lambda,maxt,bw,adaptive=1,kerneltype=1)
{
  if (data$window$type!="rectangle"){
    print("stop only rectangular observation windows allowed")
    stop
  }

  data$x=data$x-data$window$xrange[1]
  data$y=data$y-data$window$yrange[1]
  data$window$xrange=c(0,data$window$xrange[2]-data$window$xrange[1])
  data$window$yrange=c(0,data$window$yrange[2]-data$window$yrange[1])


  x=as.double(data$x)
  y=as.double(data$y)
  
  sideEW<-as.double(data$window$xrange[2])
  sideNS<-as.double(data$window$yrange[2])
  
  par = as.double(lambda)
  antpktdat = as.integer(length(x))
  maxt = as.double(maxt)
  bw = as.double(bw)
  gestimate = as.double(rep(0, 100))

  parametric=as.integer(F)
  correctype=as.integer(3)
  adaptive=as.integer(adaptive)
  kerneltype=as.integer(kerneltype)
  gestimate <- .C("gtrans",
                  antpktdat,
                  x,
                  y,
                  par,
                  maxt,
                  gestimate,
                  bw,
                  correctype, sideEW,sideNS,parametric,adaptive,kerneltype)[[6]]
  return(list(t=c(1:100)*maxt/100,g=gestimate))
}
thomasestK=function(t,K,startpar){
  theoret=function(t,par){
    pi*t^2+(1-exp(-t^2/(4*par[2]^2)))/par[1] 
  }
  object=function(par,tK){
    t=tK$t
    K=tK$K
    return(sum((K^0.25-theoret(t,par)^0.25)^2))
  }
  tK=list(t=t,K=K)
  minimum=optim(startpar,object,tK=tK)
  return(list(minimum=minimum,Ktheoret=theoret(t,minimum$par)))
}


lgcpestK=function(t,K,startpar){
  integrand=function(r,par){
    2*pi*r*exp(par[1]*exp(-r/par[2]))
  }
  object=function(par,tK){
    t=tK$t
    K=tK$K 
    Ktheoret=K
    if (t[1]!=0)
      K[1]=integrate(integrand,lower=0,upper=t[1],par=par)$value
    for (i in 1:(length(K)-1))
      Ktheoret[i+1]=Ktheoret[i]+integrate(integrand,lower=t[i],upper=t[i+1],par=par)$value
    return(sum((K^0.25-Ktheoret^0.25)^2))
  }
  tK=list(t=t,K=K)
  minimum=optim(startpar,object,tK=tK)

  Ktheoret=K
  if (t[1]!=0)
    K[1]=integrate(integrand,lower=0,upper=t[1],par=minimum$par)$value
  for (i in 1:(length(K)-1))
  Ktheoret[i+1]=Ktheoret[i]+integrate(integrand,lower=t[i],upper=t[i+1],par=minimum$par)$value
  return(list(minimum=minimum,Ktheoret=Ktheoret))
}
 .First.lib <- function(lib, pkg) library.dynam("InhomCluster", pkg, lib) 
