.packageName <- "MppMLE"
dista=function(par,parm){
  return(sqrt(sum((par-parm)^2)))
}

weights.mpp=function(inhompar0,interactpar0,inhompar,interactpar,samples){
  w=exp(c(inhompar-inhompar0,interactpar-interactpar0)%*%samples)
  sumw=sum(w)
  return(list(w=w/sumw,sumw=sumw))
}

newtonraphson.mpp=function(inhompar,interactpar,interactradii,smpsize,space,data,border=0,model,trustfactor=100,torus=F,stopcrit=10e-8,L=c(),constraint=T,fix=c()){

  
  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])



  obspoints=cbind(data$x,data$y)
  if (length(data$m)>0)
    obspoints=cbind(obspoints,data$m)
  else
    obspoints=cbind(obspoints,rep(0,length(data$x)))


  window=c(data$window$xrange[2],data$window$yrange[2])

  
  maxpar=10
  itrno=1
  noinhompar=length(inhompar)
  nointeractpar=length(interactpar)
  nopar=noinhompar+nointeractpar
  inhompar0=inhompar
  interactpar0=interactpar
  side=c(window,0,0)
  condint=0
  condext=0
  if (border>0){
    condext=1
    obspoints[,1]=obspoints[,1]-border
    obspoints[,2]=obspoints[,2]-border
    side[1]=window[1]-2*border
    side[2]=window[2]-2*border
    side[3]=border
    side[4]=border
  }
  #compute sufficient statistic
  print(c("torus",torus))
  temp=calcsuffobs(interactradii,condint,condext,obspoints,model,side,torus)
  suffobs=c(temp[1:noinhompar],temp[maxpar+c(1:nointeractpar)])
  print(c("sufficient obs",suffobs))

  if (length(fix)==0)
    fix=rep(F,nopar)
  
  samples=sim.mpp(c(inhompar0,rep(0,maxpar-noinhompar)),c(interactpar0,rep(0,maxpar-nointeractpar)),interactradii,smpsize,space,side,condint,condext,obspoints=obspoints,initpoints=matrix(0,0,0),model,torus,F)
  samples=rbind(samples[1:noinhompar,],samples[(maxpar+1):(maxpar+nointeractpar),])
  nosamp=1
  repeat {
    #beregner NR-update
    w=weights.mpp(inhompar0,interactpar0,inhompar,interactpar,samples)
    if (min(w$w)<(1/smpsize)/trustfactor){
      inhompar0=inhompar
      interactpar0=interactpar
      samples=sim.mpp(c(inhompar0,rep(0,maxpar-noinhompar)),c(interactpar0,rep(0,maxpar-nointeractpar)),interactradii,smpsize,space,side,condint,condext,obspoints=obspoints,initpoints=matrix(0,0,0),model,torus,F)
      samples=rbind(samples[1:noinhompar,],samples[(maxpar+1):(maxpar+nointeractpar),])
      nosamp=nosamp+1
      w=weights.mpp(inhompar0,interactpar0,inhompar,interactpar,samples)
    }

    u=sweep(samples,2,w$w,FUN="*")
    u=apply(u,1,sum)
    j=matrix(0,nopar,nopar)
    for (l in 1:nopar)
      for (k in 1:nopar)
        j[l,k]=sum(w$w*(samples[l,]-u[l])*(samples[k,]-u[k]))
    u=suffobs-u
    par=c(inhompar,interactpar)
    parm=rep(0,nopar)
    parm[!fix]=par[!fix]+u[!fix]%*%solve(j[!fix,!fix])
    parm[fix]=par[fix]

    if (constraint){#for certain models, positive interaction parameters are not allowed...
      interact=rep(T,nopar)
      interact[1:noinhompar]=F
      parm[parm>0 & interact]=0
    }
    
    print(c("itrno",itrno))
    print("u")
    print(u)
    #print("j")
    #print(j)
    print("parameters")
    print(par)
    print(parm)
    #checker for convergence
    if ((dista(par,parm)<=stopcrit) & (constraint | (sqrt(sum(u[!fix]^2)))<=stopcrit))
      break
    inhompar=parm[1:noinhompar]
    interactpar=parm[(noinhompar+1):(noinhompar+nointeractpar)]
    itrno=itrno+1
  }

  if (is.matrix(L)){
    L=L[,!fix]
    if (!is.matrix(L))
      L=matrix(L,nrow=1)
    wald=(parm[!fix]%*%t(L))%*%solve(L%*%solve(j[!fix,!fix])%*%t(L))%*%(L%*%parm[!fix])
  }
  else
    wald=NULL
  
  return(list(suffobs=suffobs,u=u,j=j,par=parm,nosamples=nosamp,wald=wald))  
}

newtonraphson.missing.mpp=function(smpsize,space,inhompar,interactpar,interactradii,stopcrit,obspoints,window,extension,model,trustfactor,L=c(),constraints=F,hardcore=c(),torus=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])

  
  obspoints=cbind(data$x,data$y)
  window=c(data$window$xrange[2],data$window$yrange[2])
    
  maxpar=10
  itrno=1
  noinhompar=length(inhompar)
  nointeractpar=length(interactpar)
  nopar=noinhompar+nointeractpar
  inhompar0=inhompar
  interactpar0=interactpar
  #compute sufficient statistic for points in observation window
  temp=calcsuffobs(interactradii,0,0,obspoints,model,side=c(window,0,0),0)
  suffobs=c(temp[1:noinhompar],temp[10+c(1:nointeractpar)])
  
  if (length(fixpar)>0)
    fix[hardcore]=T
  print("parameters fixed at -1e10")
  print(fix)
  
  #uncond sample
  samples=sim.mpp(c(inhompar0,rep(0,maxpar-noinhompar)),c(interactpar0,rep(0,maxpar-nointeractpar)),interactradii,smpsize,space,c(window,extension,extension),F,F,obspoints=obspoints,initpoints=matrix(0,0,0),model,torus,F)
  samplesunco=rbind(samples[1:noinhompar,],samples[(maxpar+1):(maxpar+nointeractpar),])
  #cond sample
  samples=sim.mpp(c(inhompar0,rep(0,maxpar-noinhompar)),c(interactpar0,rep(0,maxpar-nointeractpar)),interactradii,smpsize,space,c(window,extension,extension),T,T,obspoints=obspoints,initpoints=matrix(0,0,0),model,torus,F)
  samplescond=rbind(samples[1:noinhompar,],samples[(maxpar+1):(maxpar+nointeractpar),])    
  nosamp=1
  repeat {
    #beregner NR-update
    wunco=weights.mpp(inhompar0,interactpar0,inhompar,interactpar,samplesunco)
    wcond=weights.mpp(inhompar0,interactpar0,inhompar,interactpar,samplescond)
    if ((min(wunco$w)<(1/smpsize)/trustfactor) | (min(wcond$w)<(1/smpsize)/trustfactor)){
      inhompar0=inhompar
      interactpar0=interactpar
      #uncond sample
      samples=sim.mpp(c(inhompar0,rep(0,maxpar-noinhompar)),c(interactpar0,rep(0,maxpar-nointeractpar)),interactradii,smpsize,space,c(window,extension,extension),F,F,obspoints=obspoints,initpoints=matrix(0,0,0),model,torus,F)
     samplesunco=rbind(samples[1:noinhompar,],samples[(maxpar+1):(maxpar+nointeractpar),])
      #cond sample
      samples=sim.mpp(c(inhompar0,rep(0,maxpar-noinhompar)),c(interactpar0,rep(0,maxpar-nointeractpar)),interactradii,smpsize,space,c(window,extension,extension),T,F,obspoints=obspoints,initpoints=matrix(0,0,0),model,torus,F)
      samplescond=rbind(samples[1:noinhompar,],samples[(maxpar+1):(maxpar+nointeractpar),])    
      nosamp=nosamp+1
      wunco=weights.mpp(inhompar0,interactpar0,inhompar,interactpar,samplesunco)
      wcond=weights.mpp(inhompar0,interactpar0,inhompar,interactpar,samplescond)
    }

    tunco=sweep(samplesunco,2,wunco$w,FUN="*")
    tunco=apply(tunco,1,sum)
    junco=matrix(0,nopar,nopar)
    for (l in 1:nopar)
      for (k in 1:nopar)
        junco[l,k]=sum(wunco$w*(samplesunco[l,]-tunco[l])*(samplesunco[k,]-tunco[k]))
    tcond=sweep(samplescond,2,wcond$w,FUN="*")
    tcond=apply(tcond,1,sum)
    jcond=matrix(0,nopar,nopar)
    for (l in 1:nopar)
      for (k in 1:nopar)
        jcond[l,k]=sum(wcond$w*(samplescond[l,]-tcond[l])*(samplescond[k,]-tcond[k]))
    tcond=tcond+suffobs
    u=tcond-tunco
    j=junco-jcond
    par=c(inhompar,interactpar)
    parm=rep(0,nopar)
    parm[!fix]=par[!fix]+u[!fix]%*%solve(j[!fix,!fix])
    parm[fix]=-10^9

    if (constraint){
      interact=rep(T,nopar)
      interact[1:noinhompar]=F
      parm[parm>0 & interact]=0
    }
    print(c("itrno",itrno))
    print("u")
    print(u)
    print("j")
    print(j)
    print("parameters")
    print(par)
    print(parm)
    #checker for convergence
    if ((dista(par,parm)<=stopcrit) & (constraint | sqrt(sum(u^2))<=stopcrit))
      break
    inhompar=parm[1:noinhompar]
    interactpar=parm[(noinhompar+1):(noinhompar+nointeractpar)]
    itrno=itrno+1
  }

  if (is.matrix(L)){
    L=L[,!fix]
    if (!is.matrix(L))
      L=matrix(L,nrow=1)
    wald=(parm[!fix]%*%t(L))%*%solve(L%*%solve(j[!fix,!fix])%*%t(L))%*%(L%*%parm[!fix])
  }
  else
    wald=NULL
  
  return(list(suffobs=suffobs,u=u,j=j,par=parm,nosamples=nosamp,wald=wald))  
}

pathsampling.mpp=function(inhompar0,interactpar0,inhompar1,interactpar1,interactradii,nogrid,smpsize,space,data,border,model,torus=F){

  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])

  
  obspoints=cbind(data$x,data$y)
  if (length(data$m)>0)
    obspoints=cbind(obspoints,data$m)
  else
    obspoints=cbind(obspoints,rep(0,length(data$x)))

  window=c(data$window$xrange[2],data$window$yrange[2])
  
  maxpar=10
  noinhompar=length(inhompar0)
  nointeractpar=length(interactpar0)
  nopar=noinhompar+nointeractpar
  side=c(window,0,0)
  condint=0
  condext=0
  if (border>0){
    condext=1
    obspoints[,1]=obspoints[,1]-border
    obspoints[,2]=obspoints[,2]-border
    side[1]=window[1]-2*border
    side[2]=window[2]-2*border
    side[3]=border
    side[4]=border
  }

  tss=c(0:nogrid)/nogrid
  par0=c(inhompar0,interactpar0)
  par1=c(inhompar1,interactpar1)  
  pardiff=par1-par0
  part=tss%*%t(pardiff)
  part=sweep(part,2,FUN="+",par0)
  us=rep(0,nogrid+1)
  impsamp=rep(0,nogrid+1)
  
  temp=calcsuffobs(interactradii,condint,condext,obspoints,model,side,torus)
  suffobs=c(temp[1:noinhompar],temp[maxpar+c(1:nointeractpar)])
  
  
  for (i in 1:(nogrid+1)){
    print(c("sample",i))
    samples=sim.mpp(c(part[i,1:noinhompar],rep(0,maxpar-noinhompar)),c(part[i,(noinhompar+1):(noinhompar+nointeractpar)],rep(0,maxpar-nointeractpar)),interactradii,smpsize,space,side,condint,condext,obspoints=obspoints,initpoints=matrix(0,0,0),model,torus,F)
    samples=rbind(samples[1:noinhompar,],samples[(maxpar+1):(maxpar+nointeractpar),]) 
    us[i]=apply(samples,1,mean)%*%pardiff
    if (i<=nogrid)
      impsamp[i+1]=log(mean(exp((part[i+1,]-part[i,])%*%samples)))
  }
  loglike=rep(0,nogrid+1)
  for (i in 1:nogrid)
    loglike[i+1]=(us[i]+us[i+1])/(2*nogrid)
  for (i in 1:nogrid)
    loglike[i+1]=-loglike[i+1]+suffobs%*%(part[i+1,]-part[i,])
  for (i in 1:nogrid)
    impsamp[i+1]=-impsamp[i+1]+suffobs%*%(part[i+1,]-part[i,])
  
  return(list(loglikeratio=cumsum(loglike),loglikeratioimp=cumsum(impsamp),part=part,integrand=us))
}



pathsamp.miss=function(inhompar0,interactpar0,inhompar1,interactpar1,nogrid,smpsize,space,obspoints,data,extension,torus){

  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])

  
  obspoints=cbind(data$x,data$y)
  window=c(data$window$xrange[2],data$window$yrange[2])
  
  maxpar=10
  noinhompar=length(inhompar0)
  nointeractpar=length(interactpar0)
  nopar=noinhompar+nointeractpar
  side=c(window,0,0)
  condint=0
  condext=0

  tss=c(0:nogrid)/nogrid
  par0=c(inhompar0,interactpar0)
  par1=c(inhompar1,interactpar1)  
  pardiff=par1-par0
  part=tss%*%t(pardiff)
  part=sweep(part,2,FUN="+",par0)
  us=rep(0,nogrid+1)
  impsamp=rep(0,nogrid+1)
  
  temp=calcsuffobs(0,0,obspoints,model,side=c(window,0,0),0)
  suffobs=c(temp[[5]][1:noinhompar],temp[[6]][1:nointeractpar])

  
  for (i in 1:(nogrid+1)){
    print(c("sample",i))
    samples=mksim(c(part[i,1:noinhompar],rep(0,maxpar-noinhompar)),c(part[i,(noinhompar+1):(noinhompar+nointeractpar)],rep(0,maxpar-nointeractpar)),smpsize,space,0,0,obspoints=obspoints,initpoints=matrix(0,0,0),model,c(window,extension,extension),torus,0)[[3]]
  samplesunco=rbind(samples[1:noinhompar,],samples[(maxpar+1):(maxpar+nointeractpar),])
    #cond sample
    samples=mksim(c(part[i,1:noinhompar],rep(0,maxpar-noinhompar)),c(part[i,(noinhompar+1):(noinhompar+nointeractpar)],rep(0,maxpar-nointeractpar)),smpsize,space,1,0,obspoints=obspoints,initpoints=matrix(0,0,0),model,c(window,extension,extension),torus,0)[[3]]
    samplescond=rbind(samples[1:noinhompar,],samples[(maxpar+1):(maxpar+nointeractpar),])

    us[i]=(apply(samplescond,1,mean)+suffobs-apply(samplesunco,1,mean))%*%pardiff
    
    if (i<=nogrid)
      impsamp[i+1]=log(exp((part[i+1,]-part[i,])%*%suffobs))+log(mean(exp((part[i+1,]-part[i,])%*%samplescond)))-log(mean(exp((part[i+1,]-part[i,])%*%samplesunco)))
  }
  loglike=rep(0,nogrid+1)
  for (i in 1:nogrid)
    loglike[i+1]=(us[i]+us[i+1])/(2*nogrid)
  
  return(list(us=us,loglike=cumsum(loglike),loglikeimp=cumsum(impsamp),part=part))
}
sim.mpp=function(inhompar,interactpar,interactradii,smpsize,space,side=c(1,1,0,0),condint=F,condext=F,obspoints=matrix(0,0,0),initpoints=matrix(0,0,0),model=1,torus=F,boot=F){

  timeseries=matrix(as.double(0),20,smpsize)
  antobspkt=dim(obspoints)[1]
  ninitpkt=dim(initpoints)[1]

  if (antobspkt==0){
    obsx=c()
    obsy=c()
    obsm=c()
  } else {
    obsx=obspoints[,1]
    obsy=obspoints[,2]
    obsm=obspoints[,3]
  }
  if (ninitpkt==0){
    initpx=c()
    initpy=c()
    initpm=c()
  } else {
    initpx=initpoints[,1]
    initpy=initpoints[,2]
    initpm=initpoints[,3]
  }

  out=.C("markovsim",as.double(inhompar),as.double(interactpar),as.double(interactradii),timeseries,as.integer(smpsize),as.integer(space),as.integer(condint),as.integer(condext),as.double(obsx),as.double(obsy),as.double(obsm),as.integer(antobspkt),as.double(rep(0,2000)),as.double(rep(0,2000)),as.double(rep(0,2000)),as.integer(0),as.double(initpx),as.double(initpy),as.double(initpm),as.integer(ninitpkt),as.integer(model),as.double(side),as.integer(torus),as.integer(boot))[[4]]

  return(out)

}



calcsuffobs=function(interactradii,condint=F,condext=F,obspoints=matrix(0,0,0),model=1,side=c(1,1,0,0),torus=F){

  antobspkt=dim(obspoints)[1]

  timeseries=matrix(as.double(0),20,1)

  if (antobspkt==0){
    obsx=c()
    obsy=c()
    obsm=c()
  } else {
    obsx=obspoints[,1]
    obsy=obspoints[,2]
    obsm=obspoints[,3]
  }

  initpoints=obspoints[obspoints[,1]<side[1] & obspoints[,2]<side[2] & obspoints[,1]>0 & obspoints[,2]>0,]

  
  out=.C("markovsim",as.double(rep(0,10)),as.double(rep(0,10)),as.double(interactradii),timeseries,as.integer(1),as.integer(0),as.integer(condint),as.integer(condext),as.double(obsx),as.double(obsy),as.double(obsm),as.integer(antobspkt),as.double(rep(0,2000)),as.double(rep(0,2000)),as.double(rep(0,2000)),as.integer(0),as.double(initpoints[,1]),as.double(initpoints[,2]),as.double(initpoints[,3]),as.integer(dim(initpoints)[1]),as.integer(model),as.double(side),as.integer(torus),as.integer(0))[[4]]
}

 .First.lib <- function(lib, pkg) library.dynam("MppMLE", pkg, lib) 
