LDsynLikSpat.R~

library(sl)
library(LaplacesDemon)
library(testthat)
source("RanalysisFunctions.R")
source("logLik.r") # override some sl functions and add synLik

source("param.r")

## load data 
## use simulated data as in base_regression.R
source("maps_basic_regression.R") # gives maps
source("functions_migration.R")

#==================
#  Params, if modifying stats: look for multiGilStat
#==================
Nrep=500;
genIntervals <- seq(0, 250, 15) # distance classes for the general variogram

#===================
# Prep for simulations
# declare Data for LD
#===================
threshold <- 2000
dist_mat <- nearest.dist(x=sp, y=NULL, method="euclidian", delta=threshold, upper=NULL);          
dist_mat <- as.matrix(dist_mat)
probMat <- generate_prob_mat(halfDistJ, halfDistH, 1, useDelta, delta, rateHopInMove, rateSkipInMove, rateJumpInMove, threshold, sp, dist_mat)

#redefine probMat as the cumulative sum on each line
cumulProbMat<-probMat

for(line in 1:L){
		cumulProbMat[line,]<-cumsum(probMat[line,])
	}

# transpose so that redeable by lines in C
cumulProbMat<-t(cumulProbMat)

### the vector of stats for the data
infestH<-which(maps$infest3==1)
outData <- multiGilStat(cumulProbMat=cumulProbMat, infestH, timeH=maps$ages[infestH], endTime = nbit, rateMove, Nrep, coords = maps[, c("X", "Y")], breaks = genIntervals,simul=FALSE)

statsData<-outData$statsTable[,1]

### starting point for simulations
infestH <- which(maps$infest2 > 0)
timeH <- maps$ages[infestH]

### Priors (also the place to change the parameters)
paramNames<-c("rateMove")
priorMeans<-0.03

### LD formalism for Data (no setting should be made past this declaration)
PGF<-function(Data){ # parameters generating functions (for init etc...)
	priorMeans<-Data$priorMeans
  values<-rlnorm(length(priorMeans),meanlog=log(priorMeans),sdlog=0.15)
  return(values)
}


MyData <- list(y=statsData,
	     trans=NULL,
	     probMat=cumulProbMat,
	     infestH=infestH,
	     timeH=timeH,
	     endTime=nbit,
	     maps=maps,
	     Nrep=Nrep,
	     priorMeans=priorMeans,
	     genIntervals=genIntervals,
	     PGF=PGF,
	     mon.names=c("ll"), # monitored variables (like in Model)
	     parm.names=paramNames # parameters names (like in Model and Initial.Values)
	     )

#==================================
## declare Model for LD
#==================================
Model<-function(theta,Data){
	# coerce theta, a priori all positive
	theta<-interval(theta,a=0)

	# simulations
	start <- Sys.time()
	out <- multiGilStat(Data$cumulProbMat, 
			    Data$infestH, Data$timeH, endTime = nbit, 
			    theta$rateMove, 
			    Data$Nrep, 
			    coords = Data$maps[, c("X", "Y")], 
			    breaks = Data$genIntervals)
	yhat<-out$statsTable[,1]
	end <- Sys.time()
	cat("t multiGil:",end-start,"\n")

	# synthetic likelihood
	ll<-synLik(out$statsTable,Data$y,Data$trans)

	# get likelihood with priors
	LL<-ll
	attributes(LL)<-NULL
	LP<-LL+sum(dlnorm(theta,meanlog=log(Data$priorMeans),sdlog=1))

	# return
	Modelout <- list(LP=LP, # joint posterior
			 Dev=-2*LL, # deviance, probably not to be changed
			 Monitor=c(LL), # to be monitored/ploted
			 yhat=yhat, # data generated for that set of parameter
			 # will be used for posterior check
			 parm=theta # the parameters, possibly constrained by the model
			 )
} 

#===========================
# Init values and Testing Model/Data
#===========================

Initial.Values <- GIV(Model, MyData, PGF=TRUE) #GIV: generate initial values

ModelOut<-Model(Initial.Values,MyData)
stop()

#===========================
## launch LD
#===========================

#===========================
## results
#===========================