IsingOccu-core:

## misc helper functions
getintralayerGraph <- function(distM,link_map,eta,d,int_range = "exp",spp_mat) #it can be used multiple times for interislan and intra-island
{
  # pass all graphs as sparse matrix in package Matrix
  nspp <- nrow(spp_mat) # which is the interspecific neighborhood matrix
  A <- list() # intralayer graphs are passed using lists
  link_map <- as(as.matrix(link_map),"dgCMatrix")
  if(int_range=="arth"){
    A <- lapply(1:nspp,function(i,eta,distM,d){
      eta[i]*as.matrix(1/((distM)^(2+d[i])))
    },eta,distM,d)
  }
  else{
    if(int_range=="exp"){
	  A <- lapply(1:nspp,function(i,eta,d,distM,link_map){
		At <- eta[i]*as.matrix(exp(-exp(d[i])*distM)) * (link_map)
	    diag(At) <- 0
	    return(At)
	  },eta,d,distM,link_map)
    }
    else{
      if(int_range=="nn"){
	  A <- lapply(1:nspp, function(i,eta,link_map){
		    eta[i]*as.matrix((link_map))
	    },eta,link_map)
      }
      else{
        #print("int_range must be exp or arth, will assume exp")
		A <- lapply(1:nspp,function(i,eta,d,distM,link_map){
		  At <- eta[i]*as.matrix(exp(-exp(d[i])*distM)) * (link_map)
	      diag(At) <- 0
	      return(At)
	    },eta,d,distM,link_map)
      }
    }
  }
  return(A) # if link map is sparse, then A is sparse
} 
## TODO: Try to reduce number of if...else...
  # passed 2019/3/18

getfullGraph <- function(A_ex,A_in,spp_mat){
  nspp <- nrow(spp_mat)
  nsite <- nrow(A_ex[[1]])
  A <- Matrix(0,nspp*nsite,nspp*nsite,sparse = T)
  for(i in 2:nspp-1){
    A[1:nsite + (i-1)*nsite,1:nsite + (i-1)*nsite] <- A_ex[[i]] + A_in[[i]] # diagonal part
    for(j in (i+1):nspp){
      
      diag(A[1:nsite + (i-1)*nsite,1:nsite + (j-1)*nsite]) <- spp_mat[i,j]
      diag(A[1:nsite + (j-1)*nsite,1:nsite + (i-1)*nsite]) <- spp_mat[j,i]
      
    }
  }
  i <- nspp
  A[1:nsite + (i-1)*nsite,1:nsite + (i-1)*nsite] <- A_ex[[i]] + A_in[[i]]
  A <- as(A,'symmetricMatrix') # turn sparse
  return(A)
} 
  # passed 2019/3/18

mainland_thr <- function(dist_mainland,link_mainland,eta,d,int_range_inter="exp"){
	A <- 0*dist_mainland
	link_mainland <- (as.matrix(link_mainland))
	if(int_range_inter=="arth"){
			A <- eta*as.matrix(1/((dist_mainland)^(2+d)))
	}
	else{
		if(int_range_inter=="exp"){
			A <- eta*as.matrix(exp(-exp(d)*dist_mainland)) * (link_mainland)
		}
	  else{
	    if(int_range_inter=="nn")
	    A <- eta * (link_mainland)
	  }
	}
	return(A)
	# test for 2spp passed 3/18/2019
}


getMRF <- function(theta,envX,distM,link_map,dist_mainland,link_mainland,int_range_intra="nn",int_range_inter="exp"){
  	nsite <- nrow(envX)
	beta_occu <- theta$beta_occu
	eta_intra <- theta$eta_intra # intra spp, intra island if apply
	d_intra <- theta$d_intra
	spp_mat <- theta$spp_mat
	nspp <- nrow(spp_mat)
	A_in <- getintralayerGraph(distM,link_map$intra,eta_intra,d_intra,int_range = int_range_intra,spp_mat)
	eta_inter <- theta$eta_inter # assume there is one eta_inter 
	d_inter <- theta$d_inter
	A_ex <- getintralayerGraph(distM,link_map$inter,eta_inter,d_inter,int_range = int_range_inter,spp_mat) # graph among islands, if apply, distM should only contain graph 
	A <- getfullGraph(A_ex,A_in,spp_mat)
    rm(A_ex,A_in)
	thr <- lapply(1:nspp,
	  function(i,envX,beta_occu,dist_mainland,link_mainland,eta_inter,d_inter,int_range_inter){
	    envX %*% beta_occu[1:ncol(envX)+(i-1)*ncol(envX)] + 
		      mainland_thr(dist_mainland,link_mainland,eta_inter[i],d_inter[i],int_range_inter)
	    },envX,beta_occu,dist_mainland,link_mainland,eta_inter,d_inter,int_range_inter)
	thr <- Reduce(rbind,thr)
    return(list(A = A,thr = thr)) # return the Markov random field
}

IsingStateProb <- function (s, graph, thresholds, beta, responses = c(-1L, 1L)) 
{
  if (!is.list(s)) 
    s <- list(s)
    N <- length(s[[1]])
    Z <- PartitionCpp(graph,thresholds,beta,responses)
    sapply(s, function(x) exp(-beta * H(graph, x, ( thresholds)))/Z)
}

Hamiltonian <- function(MRF,Z_vec){
	nrep <- ncol(Z_vec)
	Ham <- lapply(1:nrep,function(i,Z,J,h){H(J,Z[,i],h)},Z=Z_vec,J=MRF$A,h=( MRF$thr))
  
	Ham <- Reduce(rbind,Ham)
	return(Ham) # if we have repeat, just make Z_vec has two cols 
	
}


rIsingOccu_multi <- function(MRF,n=1,method = "CFTP",nIter = 100){
	Z <- IsingSamplerCpp(n=n,graph = MRF$A,thresholds=MRF$thr, responses = matrix( c(-1L, 1L),2,1),beta = 1,nIter=nIter,exact = (method=="CFTP"),constrain = NA + MRF$thr)
  return(t(Z))
	# test for 2spp case, passed 3/18/2019
}
  
Pdet_Ising_single_site <- function(thr, Z, dethis, sppmat_det){
	spp_exist <- c(Z)==1
  sppmat_det <- as(sppmat_det,"dgCMatrix")
	if(sum(spp_exist)==0 | sum(!is.na(dethis))==0){return(0)} # no species there, probability one to be no detection, or no observation here
	if(prod(spp_exist)==0){
	  thr_exis <- as.matrix( thr[,spp_exist])
	  # thr_abs <- - apply(matrix(sppmat_det[!spp_exist,spp_exist],sum(!spp_exist),sum(spp_exist)),2,sum) # condition on some species not exist here thus never be detected 
	  # do not include thr_abs since if two species never coexist we cannot infer "what if they coexist", i.e. thr_exis will be total colinear with thr_exis
	  thr <- thr_exis
	  #thr <- apply(matrix(1:ncol(thr_exis)),1,function(k,ww,kk){ww[,k]+kk[k]},thr_exis,( thr_abs))
	}
	graph <- sppmat_det[spp_exist,spp_exist]
	has_obs <- !(is.na(rowSums(dethis)))
	dethis <- dethis[has_obs,spp_exist]# convert it to nrow = nperiod, ncol = nspp for single site, single repeat
	thr <- thr[has_obs,]
	
	Pdet_site <- apply(matrix(1:sum(has_obs)),1,function(k,dethis,thr,graph){
		IsingStateProb(dethis[k,], graph, thr[k,], beta=1, responses = c(-1L, 1L))
	} ,matrix( dethis,sum(has_obs),sum(spp_exist)), matrix( thr,sum(has_obs),sum(spp_exist)), as( as.matrix(graph),'dsCMatrix'))
	
	return(sum(log(Pdet_site + 1e-15)))
	
}


extract_thr <- function(i,thr_list){
	nspp <- length(thr_list)
	thr <- sapply(thr_list,function(thr1,i){t(thr1[i,])},i=i) # thr at site i for all spps, will return a matrix with ncol = nspp, nrow = nperiod
	return(thr)
}

Pdet_Ising <- function(nperiod,envX,detX,beta_det,sppmat_det,Z,detmat){
	#require(IsingSamplerCpp)
	 # This is the full design matrix list of detection probability p at time
	if(is.null(detX)) {
	  detDesign <- lapply(1:nperiod,function(dummy,envX){envX},envX)
	  
	}
	else detDesign <- lapply(detX,function(x,y){ as.matrix( cbind(y,x))},y = envX)
	npardet <- ncol(detDesign[[1]])
	nsite <- nrow(envX)
	nspp <- nrow(sppmat_det)
	thr_list <- lapply( 1:nspp, function(i,detDesign,beta_det,naprdet,n_row,nperiod){ 
		temp <- lapply(detDesign,function(w,beta1,i){w%*%beta1},beta1 = matrix( beta_det[1:npardet + (i-1) * npardet]),i=i)
		thr1 <- (matrix(unlist(temp),nrow = n_row,ncol = nperiod))
		return(thr1) # now here is basically a matrix, for each species at site and period
		},detDesign,beta_det,npardet,nrow(envX),nperiod) # this is gonna be  a list for all species, 
	
	Pdet <- lapply(1:nsite,function(i,thr_list,detmat,Z,sppmat_det,nsite,nspp,nperiod){
		thr1 <- extract_thrCpp(i-1,thr_list,nspp,nperiod,nsite)
		rows1 <- i + (1:nspp-1)*nsite
		dethis <- t(detmat[rows1,])
		Z_site <- Z[rows1,]
		Pdet_Ising_single_siteCpp(thr1, Z_site, dethis, sppmat_det,c(-1L,1L)  )
	},thr_list,detmat,as.matrix( Z), as.matrix( sppmat_det),nsite,nspp,nperiod)# loop over sites
	return(Reduce(rbind,Pdet)) # change 25/8/2019
}


Pdet_Ising_rep <- function(nrep,nperiod,envX,detX,beta_det,sppmat_det,Z,detmat){
  Pdets <- lapply(1:nrep,function(k,nperiod,envX,detX,beta_det,sppmat_det,Z,detmat){
    sum(Pdet_Ising(nperiod,envX,detX[[k]],beta_det,sppmat_det,Z[,k],detmat[[k]])) # change 25/8/2019
  },nperiod,envX,detX,beta_det,sppmat_det,Z,detmat)
  return(Reduce('+',Pdets))
}



IsingOccu_Ising_det_multi_logL_innorm <- function(MRF, beta_det,sppmat_det, Z ,detmat, envX,detX){ # the in-normalized log likelihood of IsingOccu Model beta is matrix here detX should be a list of list detmat should be a list, they should have the same length
  nspp <- nrow(sppmat_det)
  nperiod <- ncol(detmat[[1]])
  negPot <- Hamiltonian(MRF,Z)
  negPot <- -sum(negPot)
  nrep <- ncol(Z)
  logLdata <- Pdet_Ising_rep(nrep,nperiod,envX,detX,beta_det,sppmat_det,Z,detmat)
  return(negPot+logLdata)
}


         
write_json.IsingOccu_samples <- function(x,path){
  n_sample <- nrow(x$Z.mcmc)
  x$theta.mcmc <- lapply(x$theta.mcmc,matrix,nrow = n_sample)
  x$Z.mcmc <- matrix(x$Z.mcmc,nrow = n_sample)
  class(x) <- 'list'
  jsonlite::write_json(x,path)
}
		   
		   
adjacency.matrix <- function(m, n = NULL)
{
    if (missing(n))
    {
        A <- Matrix(0, m^2, m^2,sparse = T)
        for (i in 1:m^2)
        {
            up <- i - m
            down <- i + m
            left <- i - 1
            right <- i + 1
            if (up > 0)
                A[i, up] <- 1
            if (down <= m^2)
                A[i, down] <- 1
            if (left %% m != 0)
                A[i, left] <- 1
            if (i %% m != 0)
                A[i, right] <- 1
        }
    }
    else
    {
        A <- Matrix(0, m * n, m * n,sparse = T)
        for (i in 1:(m * n))
        {
            up <- i - n
            down <- i + n
            left <- i - 1
            right <- i + 1
            if (up > 0)
                A[i, up] <- 1
            if (down <= (m * n))
                A[i, down] <- 1
            if (left %% n != 0)
                A[i, left] <- 1
            if (i %% n != 0)
                A[i, right] <- 1
        }
    }
    A
}