Commit c40732e6 authored by jfuser's avatar jfuser
Browse files

rename to cpp src files

plug R and c++ in demo
+ gen doc
parent bee9834e
Package: briskaRepi
Package: landsepi
Type: Package
Encoding: UTF-8
Title: Biological Risk Assessment for epidemology
Title: Landscape Epidemio
Version: 0.0.1
Date: 2016-01-04
Date: 2017-11-23
Authors@R: c(person("Loup", "Rimbaud", role = "aut", email = "loup.rimbaud@csiro.au"),
person("Jean-Francois", "Rey", role = "cre", email = "jean-francois.rey@inra.fr"),
person("Julien", "Papaix", role = c("aut","cre"), email = "julien.papaix@inra.fr"))
......@@ -31,20 +31,19 @@ Imports:
rgeos,
maptools,
fields,
splancs,
sf
Suggests:
knitr(>= 1.11),
rmarkdown(>= 0.8.1),
testthat(>= 1.0.0)
Collate:
'AgriLand.R'
'RcppExports.R'
'allocation.R'
'graphLand.R'
'periodic_cov.R'
'multiN.R'
'landscapeGenerator.R'
'graphLand.R'
'AgriLand.R'
'RcppExports.R'
'demo_landsepi.R'
'landsepi.R'
'output.R'
LinkingTo: Rcpp
RoxygenNote: 6.0.1
File mode changed from 100644 to 100755
# Generated by roxygen2: do not edit by hand
export(AgriLand)
export(demo_landepi)
export(modeleLandsLPI)
import(Rcpp)
import(fields)
import(maptools)
import(methods)
import(raster)
import(rgdal)
import(rgeos)
import(sp)
import(splancs)
import(stats)
useDynLib(landsepi)
#' @title AgriLand
#' @incluce multiN.R periodic_cov.R landscapeGenerator.R graphLand.R
#' @include multiN.R periodic_cov.R graphLand.R
#' @export
AgriLand <- function(landscape,filename="landscapeTEST",propSR,isolSR,propRR,isolRR,strat,Nhote,nYears,Cmax0,Cmax1){
......
......@@ -9,10 +9,10 @@
#' @param landscape list of landscape file name and layer name
#' @param inits list inits variables
#' @param val_seed seed value
#' @param host
#' @param host host param
#' @return nothing yet
#' @export
modeleLandsLPI <- function(times, landscape, dispersion, inits, val_seed, host) {
invisible(.Call('_briskaRepi_modeleLandsLPI', PACKAGE = 'briskaRepi', times, landscape, dispersion, inits, val_seed, host))
invisible(.Call('_landsepi_modeleLandsLPI', PACKAGE = 'landsepi', times, landscape, dispersion, inits, val_seed, host))
}
#################################################
#### allocation.R ####
#################################################
## Allocation of cultivar using multivariate normal distribution
mainPath <- "MAINPATH"
pathRES <- "PATHRES"
graphOn <- GRAPHON
NpolyTarg <- NPOLYTARG
Nhote <- NHOTE
nYears <- NYEARS
idLAN <- IDLAN
seed <- SEED
prop <- c(PROPSR, PROPRR)
isol <- c(ISOLSR, ISOLRR)
strat <- "STRAT"
# mainPath <- "/home/loup/Documents/csiro/model_resistance"
# pathRES <- "/home/loup/Documents/csiro/model_resistance/RES/test/"
# graphOn <- 1
# NpolyTarg <- 500
# Nhote <- 2
# nYears <- 50
# idLAN <- 1
# seed <- 1
# prop <- c(3/6, 0.5)
# isol <- c(3,1)
# strat <- "PY"
# library(rgdal)
library(RCALI) ## function readpoly2
library(fields) ## function Exponential
library(MASS) ## function mvnorm
library(Matrix) ## function nearPD to convert a covariance matrix in a positive-definite matrix
set.seed(seed)
## Parameters of cultivar allocation
nAlloc <- (Nhote>1) ## basic number of allocations to perform
if (strat=="MO" & Nhote>2)
nAlloc <- 2
## Isolation/aggregation parameter
aggreg <- c(-180, 200, -2000, 0)
## selection of the appropriate aggregation parameter with isol
## isol = 1 --> low
## isol = 2 --> moderate
## isol = 3 --> high
## isol = 4 --> random
## Landscape structure
landscape <- readpoly2(paste(mainPath,"/landscapes/p",NpolyTarg,"/p", NpolyTarg, "-",idLAN,".txt",sep=""), delim="\t")
area <- scan(paste(mainPath,"/landscapes/p",NpolyTarg,"/p",NpolyTarg,"-",idLAN,"_area.txt",sep=""),sep=",", dec=".")
nPoly <- length(landscape)
## Host dispersal matrix
# write(as.vector(diag(nPoly)),file=paste(mainPath,"/landscapes/p",NpolyTarg,"-",idLAN,"_dispH.txt",sep=""),sep=",")
## Graphic
# plot(landscape, color=F, border="darkgrey", lwd=1.5)
# title(paste("Landscape structure",idLAN))
# mtext("m", side=2, line=2.5, padj=0, cex=1, las=1)
# mtext("m", side=1, line=0, padj=3.5, cex=1, las=1)
# summary(unlist(landscape)) ## Bounds of the landscape
## Centroid of the paddocks
centroid <- data.frame(x=NA, y=NA)
for (j in 1:nPoly)
centroid[j,] <- apply(landscape[[j]],2,mean)
d <- as.matrix(dist(centroid)) ## 2-by-2 distance between centroid of each paddock
neigh <- (d <= 399) ## matrix of neighborood
## periodic covariance function
periodic_cov <- function(d,range,phi=1) { return( exp(-2 * sin(d*pi/(2*range))^2 / phi^2) ) }
## take the matrix d to compute the allocation of cultivars
multiN <- function(d, area, aggreg, prop) {
nPoly <- nrow(area)
## Multivariate normal distribution
if (aggreg != 0) {
if (aggreg > 0) ## aggregation
covMat <- Exponential(d, aggreg) ## function from 'fields' package
if (aggreg < 0) { ## repulsion
covMat_tmp <- periodic_cov(d, aggreg)
covMatPD <- nearPD(covMat_tmp, keepDiag=TRUE, ensureSymmetry=TRUE, maxit=1000) ## conversion in a positive-definite matrix
if (covMatPD$converge==TRUE) {
print("OK: Periodic covariance could be converted in a positive-definite matrix")
covMat <- covMatPD$mat
} else {
print(paste("WARNING: covariance cannot be converted in a positive-definite matrix -- Aggreg =", aggreg))
print("Exponential kernel used instead")
aggreg2 <- 1*(aggreg > -200) + 2000*(aggreg < -200)
covMat <- Exponential(d, abs(aggreg2))
}
} ## if aggreg < 0
s <- mvrnorm(1,mu=rep(0,nPoly),Sigma=covMat) ## function from 'MASS' package
} else {
s <- rnorm(nPoly,mean=0,sd=1) ## random allocation
}
area$s <- s
## Ordered gaussian beam
ord_S <- order(s)
area.ord_S <- area[ord_S,]
area.ord_S$cumsum <- cumsum(area.ord_S$area)
## Threshold for allocation of the R cultivar
prop.areaTot <- prop * area.ord_S$cumsum[nPoly]
area.ord_S$dif <- abs(area.ord_S$cumsum - prop.areaTot)
th_index <- which.min(area.ord_S$dif)
area.ord_S$cultivar <- as.numeric( (1:nPoly) < th_index )
## Final landscape
habitat <- area.ord_S[order(area.ord_S$num), c("num", "cultivar")]
return(habitat)
# plot(centroid, col=grey((s-min(s))/diff(range(s))), pch=16, cex=2)
# plot(centroid, col=(s<0)+1, pch=16, cex=2)
}
## Multivariate distribution
area.df <- data.frame(num=1:nPoly, area)
habitat <- data.frame(num=1:nPoly, cultivar=rep(0,nPoly))
i=0
num_index=1
while( sum(num_index)>0 & i<nAlloc ) {
## update area, d, qnd i for allocation of the cultivar
area.tmp <- area.df[habitat$cultivar==i,]
d.tmp <- d[habitat$cultivar==i, habitat$cultivar==i]
i <- i+1
## Compute the multivariate normal distribution
habitat.tmp <- multiN(d.tmp, area.tmp, aggreg[isol[i]], prop[i])
## update habitat by allocating cultivar i
num_index <- habitat.tmp$cultivar==1
habitat[habitat.tmp$num[num_index],"cultivar"] <- i
}
## Writing the habitat files for C function
habitat1 <- habitat$cultivar
habitat2 <- habitat1
if (strat=="RO" | strat=="TO"| strat=="MI")
habitat2[habitat2==1] <- 2
write(as.vector(habitat1),file=paste(pathRES,"/habitat1.txt",sep=""),sep=",")
write(as.vector(habitat2),file=paste(pathRES,"/habitat2.txt",sep=""),sep=",")
if (graphOn==1) {
source(paste(mainPath,"/scripts/graphLand.R",sep="")) ## function plot.land + invlogit
title.hab <- "Simulated landscape"
isol.name <- c("low", "medium", "high", "random")
alloc.name <- c("R/(S+R)","R2/(R1+R2)")
subtitle.hab <- paste("Cropping ratio",alloc.name,"=",round(prop,2), " Aggregation",alloc.name,"=",isol.name[isol])
if (strat!="MO" | Nhote==2)
subtitle.hab <- subtitle.hab[1]
col.hab <- c("white", "gray80","gray45") ## habitat 0 = S ; habitat 1 = R1 ; habitat 2 = R2
dens.hab <- c(0,0,0)
angle.hab <- c(0,0,0)
legend.hab <- c("Susceptible","Resistant")
if ((strat=="MO" | strat=="RO" | strat=="TO") & Nhote>2)
legend.hab <- c("Susceptible","Resistant 1","Resistant 2")
if (strat=="MI")
legend.hab <- c("Susceptible","Resistant 1 + Resistant 2")
if (strat=="PY")
legend.hab <- c("Susceptible","Resistant 1+2")
png(file=paste(pathRES,"/landscape",idLAN, "_propSR",round(prop[1],2),"_propRR",round(prop[2],2),"_isolSR",isol[1],"_isolRR",isol[2],"_",seed,".png",sep=""),width=1000,height=1000)
plot.land(landscape, col.hab[habitat1+1], dens.hab[habitat1+1], angle.hab[habitat1+1], col.hab, dens.hab, angle.hab, title.hab, subtitle.hab, legend.hab)
dev.off()
if (strat=="RO" | strat=="TO") {
png(file=paste(pathRES,"/landscape",idLAN, "_propSR",round(prop[1],2),"_propRR",round(prop[2],2),"_isolSR",isol[1],"_isolRR",isol[2],"_",seed,"bis.png",sep=""),width=1000,height=1000)
plot.land(landscape, col.hab[habitat2+1], dens.hab[habitat2+1], angle.hab[habitat2+1], col.hab, dens.hab, angle.hab, title.hab, subtitle.hab, legend.hab)
dev.off()
}
} ## if graphOn
## Crop rotations & Turn-over: calculation of time series
rotation <- data.frame(y=1:(nYears+1), habitat=rep(0,nYears+1)) ## need to have nYears+1 because of C algorithm for host plantation
if (strat=="RO" | strat=="TO") {
## isolRR gives the number of years for each habitat
rotation$habitat <- rep(c(1,0), each=isol[2], length=nYears+1)
## Graphic
if (graphOn==1) {
landscape.rot <- crlistpoly()
for (y in 1:nYears) {
x.poly=c(y-1,y-1,y,y)
y.poly=c(0,1,1,0)
# polygon(x.poly,y.poly,col=rotation$habitat[y]+1)
landscape.rot[[y]] <- as.poly(x.poly,y.poly)
}
habitat.rot <- c(max(habitat1),max(habitat2))[rotation$habitat+1] ## each year, index of the cultivar used as resistant
col.rot <- col.hab[2:3]
dens.rot <- dens.hab[2:3]
angle.rot <- angle.hab[2:3]
title.rot <- "Simulated crop rotation"
subtitle.rot <- paste("Cropping ratio",alloc.name,"=",round(prop,2), " Aggregation",alloc.name,"=",isol.name[isol])
subtitle.rot <- subtitle.rot[2]
legend.rot <- c("Resistant 1","Resistant 2")
png(file=paste(pathRES,"/rotation_propRR",round(prop[2],2),"_isolRR",isol[2],"_",seed,".png",sep=""),width=500,height=270)
par (mar=c(3,0,4.5,0), cex=0.5)
plot.land(landscape.rot, col.rot[habitat.rot], dens.rot[habitat.rot], angle.rot[habitat.rot]
, col.rot, dens.rot, angle.rot, title.rot, subtitle.rot, legend.rot, XMAX=nYears, YMAX=1)
for (y in 1:nYears)
text(y-0.5,0.5,labels=y, cex=0.9)
dev.off()
} ## if graphOn
}
## Not useful any more
# if (strat=="RO" | strat=="TO") {
# rho1=c(-1, NA, 1, 0, -1, 0.1, 1) ## weight of previous year
# rho2=c( 0, 0, 0, 1, 1, 1, 1) ## weight of stochasticity
# ## index 1 = deterministic rotation 1 year each
# ## index 2 = deterministic rotation 2 year each
# ## index 3 = deterministic turn-over
# ## index 4 = stochastic random
# ## index 5 = stochastic low
# ## index 6 = stochastic moderate
# ## index 7 = stochastic high
# id_rho <- isol[2]
# if (strat=="TO")
# id_rho <- 3
#
# s <- numeric()
# s[1] <- -1
# for (i in 2:nYears)
# s[i] <- rho1[id_rho] * s[i-1] + rho2[id_rho] * rnorm(1,0,1)
#
# thres <- round(nYears * prop[2]) ## threshold for habitat allocation
# id_hab1 <- order(s)[1:thres]
#
# rotation$habitat[id_hab1] <- 1
#
# if (is.na(rho1[id_rho])) {
# ## 2 years for each habitat --> complete cycle in 4 years
# for (i in 1:nYears)
# rotation$habitat[i] <- ((i+1) %% 4 > 1)
# }
#
# # plot(rotation$s, type="l")
# # points(rotation$s, col=rotation$habitat+1, pch=16, cex=2)
# }
write(as.vector(rotation$habitat),file=paste(pathRES,"/rotation.txt",sep=""),sep=",")
# source("graphpoly.R")
# plot.epi(landscape, habitat, 2, hab=NULL)
## Read landscape (fonction classique de RCALI + récupère nom du polygone)
# readpoly2.custom <- function (filename, delim = "\t") {
# lu = readLines(filename)
# npoly = as.integer(lu[1])
# retour = list()
# k = 1
# nomP <- NULL
# for (i in 1:npoly) {
# k = k + 1
# id = strsplit(lu[k], delim)[[1]][1]
# id = id[id != ""]
# nomP = c(nomP,as.integer(strsplit(lu[k], delim)[[1]][2]))
# k = k + 1
# xlu = as.double(strsplit(lu[k], delim)[[1]])
# k = k + 1
# ylu = as.double(strsplit(lu[k], delim)[[1]])
# coo = matrix(c(xlu, ylu), ncol = 2)
# dimnames(coo) = list(NULL, c("xcoord", "ycoord"))
# retour[[id]] = coo
# }
# names(retour) <- nomP
# class(retour) <- "listpoly"
# attr(retour, "call") <- match.call()
# return(retour)
# }
#
# paysage <- readpoly2.custom("p150-1.txt")
#-----------------------------------------------------------------------
......@@ -14,7 +14,7 @@
#' Run demo landepi
#' @include RcppExports.R AgriLand.R graphLand.R landscapeGenerator.R multiN.R periodic_cov.R
#' @include RcppExports.R AgriLand.R graphLand.R multiN.R periodic_cov.R
#' @export
demo_landepi <- function(){
......@@ -38,5 +38,8 @@ landscape <- AgriLand(landscapeTEST,filename="paysageTEST",propSR,isolSR,propRR,
#plot(habitat)
#plot(resultat)
nTSpY = 120
modeleLandsLPI(times = list(nYears=nYears,nTSpY=nTSpY) , landscape, dispersion=list(dispP="fichier"), inits=list(C_0=0.1, PI0=5e-4), val_seed=12345, host=list(Nhote=3))
}
......@@ -4,7 +4,7 @@ logit <- function(p) log(p/(1-p))
invlogit <- function(n) exp(n) / (1+exp(n))
#### FUNCTION TO DRAW THE LANDSCAPE ####
library("splancs") ## function polymap
#library("splancs") ## function polymap
plot.land <- function(landscape, COL=rep(0,length(landscape)), DENS=rep(0,length(landscape)), ANGLE=rep(30,length(landscape))
, COL.LEG=unique(COL), DENS.LEG=unique(DENS), ANGLE.LEG=unique(ANGLE)
......@@ -25,4 +25,4 @@ plot.land <- function(landscape, COL=rep(0,length(landscape)), DENS=rep(0,length
legend(XMAX/2.66, -YMAX/40, legend=LEGEND1, col="black", density=2*DENS.LEG, angle=ANGLE.LEG, bty="n")
legend(-XMAX/5, YMAX, legend=LEGEND2, fill=COL.LEG, bty="n", title=TITLE.LEG2)
}
}
\ No newline at end of file
}
###########################################################
#### Landscape structure & Dispersal kernel ####
###########################################################
## Tesselation objects are in class S4
library(RLiTe)
## parameters
NpolyTarg <- 150
setwd(paste("/home/loup/Documents/csiro/model_resistance/landscapes/p", NpolyTarg, sep=""))
Niter <- 50000 ## number of MCMC iterations
idLAN <- 6 ## index of the landscapes
itau <- c(3) #47
ialpha <- 0
ibeta <- 0 #3
plan <- expand.grid(idLAN=idLAN, itau=itau, ialpha=ialpha, ibeta=ibeta)
# plan <- data.frame(cbind(itau, ialpha, ibeta))
N <- nrow(plan)
nPoly <- numeric()
plan
# par(mfrow=c(4,4))
for (i in 1:N) {
print(paste(i,"/",N))
## T-tesselation
# seed <- sample(1:152384,1)
# print(paste("seed =",seed))
# seed <- plan$idLAN[i]
seed <- 14
setLiTeSeed(seed) ## seed
myTes <- new(TTessel) ## start an empty tesselation (object of class "TTessel")
myTes$setDomain(1,1) ## dimensions of the landscape
## Model for T-tesselation (object of class "model" to be applied to an existing tesselation)
# myMod <- ModelCRTT(tau=1.9) ## Completely Random T-Tesselation
# myMod <- ModelSquaredAreas(tau=0.0043,alpha=10000) ## Penalty on the area
if (plan$ialpha[i] == 0) {
myMod <- ModelAcuteAngles(tau=plan$itau[i],beta=plan$ibeta[i]) ## Penalty on the angles
}else{
myMod <- ModelAreasAngles(tau=plan$itau[i],alpha=plan$ialpha[i],beta=plan$ibeta[i]) ## Penalty on both areas and angles
}
myMod$set_ttessel(myTes) ## application of the model to the tesselation
mySim <- new(SMFChain,myMod,0.33,0.33) ## generation of the landscape by Markov Chain Monte Carlo (0.33: prob of Split/Merge/Flip events)
mySim$step(Niter) ## number of iterations of the Markov Chains
areas <- myTes$getCellAreas()
# areas <- areas * landHeight * landWidth / 1E4 ## area in ha
nPoly[i] <- length(areas)
## Graphics
# hist(areas, col="grey")
# plot(myTes, main=paste("tau =", plan$itau[i], "\nalpha =", plan$ialpha[i], "\nbeta =", plan$ibeta[i])
# , sub=paste("nPoly =", nPoly[i], "\nMedian area =", round(median(areas),1), "[", paste(round(quantile(areas),1)[c(2,4)], collapse="-"), "]"))
## Exportation of the landscape
sink(paste("p",NpolyTarg,"-",plan$idLAN[i],"-tmp.txt",sep="")) ## open
myTes$printRCALI() ## write with format RCALI
sink() ## close
}
####################################################
#### Generation of a regular lattice ####
####################################################
regular_lattice <- function(fileOut=getwd(), Nx=2, Ny=2) {
## Coordinates of the polygons
dim <- c(Nx, Ny)
nPoly <- prod(dim)
coord.x <- list()
coord.y <- list()
for (i in 1:Nx)
coord.x[[i]] <- (i-1)*(1/Nx) + c(0, 1/Nx, 1/Nx, 0)
for (j in 1:Ny)
coord.y[[j]] <- (j-1)*(1/Ny) + c(0, 0, 1/Ny, 1/Ny)
landscape <- list()
l=0
for (i in 1:Nx) {
for (j in 1:Ny) {
l <- l+1
landscape[[l]] <- cbind(coord.x[[i]], coord.y[[j]])
}
}
## Write the landscape
cat(nPoly,file=fileOut)
for (l in 1:nPoly) {
n <- 4 ## number of vertices
cat(NULL,file=fileOut,append=T,sep="\n")
cat(c(l,l,n),file=fileOut,append=T,sep=" ")
cat(NULL,file=fileOut,append=T,sep="\n")
cat(landscape[[l]][,1],file=fileOut,append=T,sep=" ")
cat(NULL,file=fileOut,append=T,sep="\n")
cat(landscape[[l]][,2],file=fileOut,append=T,sep=" ")
}
cat(NULL,file=fileOut,append=T,sep="\n")
# landscape <- readpoly2(fileOut, delim=" ")
# plot(landscape, color=F, border="darkgrey", lwd=1.5)
}
regular_lattice(fileOut=paste("p",NpolyTarg,"-",7,"-tmp.txt",sep=""), Nx=12, Ny=12)
regular_lattice(fileOut=paste("p",NpolyTarg,"-",8,"-tmp.txt",sep=""), Nx=4, Ny=40)
####################################
#### Dispersal Kernel ####
####################################
# graphics.off()
# library(RCALI)
# ?califlopp
library(splancs) ## function polymap
source("/home/loup/Documents/csiro/model_resistance/RCALI/R/sourceDir.R")
sourceDir("/home/loup/Documents/csiro/model_resistance/RCALI/R")
dyn.load("/home/loup/Documents/csiro/model_resistance/RCALI/libs/RCALI.so")
## power-law distribution with parameters a = 40 and b = 7 (f1), 3.4 (f2), 3.04 (f3), 2.5 (f5)
## CAREFUL!! BUG with califlopp function: dispersal matrix must be generated one at a time
## (i.e. not all the dispersal kernels in the same time)
NpolyTarg <- 150
landHeight <- 2000 ## landscape length (m)
landWidth <- 2000 ## landscape width (m)
dispF <- 1:4
lan1 <- 1
lanE <- 8
for (idLAN in lan1:lanE) {
print(idLAN)
# fileIn.tmp <- paste("p",NpolyTarg,"-",idLAN,"-tmp.txt",sep="")
fileIn <- paste("p",NpolyTarg,"-",idLAN,".txt",sep="")
fileOut <- paste("p", NpolyTarg, "-", idLAN, "_f.res",sep="")
## Conversion of the landscape for califlopp function
## replace " " by "\t"
## re-order the vertices in clockwise order
## multiply by landscape length and width
## add 1 because califlopp doesn't like coordinates at 0
landscape <- readpoly2(fileIn.tmp, delim=" ")
nPoly <- length(landscape)
cat(nPoly,file=fileIn)
for (l in 1:nPoly) {
n <- nrow(landscape[[l]])
cat(NULL,file=fileIn,append=T,sep="\n")
cat(c(l,l,n),file=fileIn,append=T,sep="\t")
cat(NULL,file=fileIn,append=T,sep="\n")
cat(landscape[[l]][n:1,1]*landWidth+1,file=fileIn,append=T,sep="\t")
cat(NULL,file=fileIn,append=T,sep="\n")
cat(landscape[[l]][n:1,2]*landHeight+1,file=fileIn,append=T,sep="\t")
}
cat(NULL,file=fileIn,append=T,sep="\n")
## Graphic
landscape <- readpoly2(fileIn, delim="\t")
nPoly <- length(landscape)
write(nPoly, "nPoly.txt", append=TRUE)
png(file=paste("graphics/p",NpolyTarg,"-",idLAN,".png",sep=""),width=1000,height=1000)
par(cex=1, cex.axis=1.5, cex.main=2)
plot(landscape, color=F, border="darkgrey", lwd=1.5)
title(paste("Landscape structure",idLAN))
mtext("m", side=2, line=2.5, padj=0, cex=1.5, las=1)
mtext("m", side=1, line=0, padj=3.5, cex=1.5, las=1)
dev.off()
## Computation of the dispersal kernel: power-law distribution
for (f in dispF) {
param <- list(input=2, output=0, method="cub", dp=rep(6000, 1), dz=rep(6000, 1), warn.poly=FALSE, warn.conv=FALSE, verbose=FALSE)
disp.res <- califlopp(file=fileIn, dispf=f, param=param, resfile=fileOut)
## Vectors of polygons and areas
disp.res <- getRes(fileOut)
#distpoly <- read.table(paste(,sep=""),header=T)
emitter <- c(disp.res$poly1,disp.res$poly2)
receiver <- c(disp.res$poly2,disp.res$poly1)
area_e <- c(disp.res$area1,disp.res$area2)
area_r <- c(disp.res$area2,disp.res$area1)
area <- as.vector(by(area_e,emitter,mean))
write(area,file=paste("p",NpolyTarg,"-",idLAN,"_area.txt",sep=""),sep=",")
## Generation of the dispersal matrix
# name.f <- paste("mean.flow.f", f, sep="")
# if (length(dispF)==1)
name.f <- "mean.flow"
flow.mean <- c(disp.res[,name.f],disp.res[,name.f])
flow.f <- cbind(emitter,receiver,flow.mean,area_e,area_r)
## remove the doublons (i.e. half the lines where emitter == receiver)
flow.f[1:nrow(disp.res),][(disp.res$poly2-disp.res$poly1)==0,] <- NA
flow.f <- flow.f[is.na(apply(flow.f,1,sum))==F,]
flow.f <- as.data.frame(flow.f)
colnames(flow.f) <- c("emitter","receiver","flow","area_e","area_r")
flow.f <- flow.f[order(flow.f$emitter),]