Commit 44c5d58b authored by jfuser's avatar jfuser
Browse files

Merge branch 'master' of gitlab.paca.inra.fr:CSIRO-INRA/landsepi

parents 013f2d01 a018cdf2
......@@ -34,18 +34,23 @@ Imports:
maptools,
fields,
splancs,
sf
sf,
RCALI
Suggests:
knitr(>= 1.11),
rmarkdown(>= 0.8.1),
testthat(>= 1.0.0)
Collate:
'invlogit.R'
'logit.R'
'RcppExports.R'
'graphLand.R'
'periodic_cov.R'
'multiN.R'
'AgriLand.R'
'RcppExports.R'
'HLIRdynamics.R'
'demo_landsepi.R'
'landsepi.R'
'plot.evol.QR.R'
LinkingTo: Rcpp
RoxygenNote: 6.0.1
# Generated by roxygen2: do not edit by hand
export(AgriLand)
export(HLIRdynamcis)
export(demo_landsepi)
export(modeleLandsLPI)
export(plot.land)
import(MASS)
import(RCALI)
import(Rcpp)
import(fields)
import(graphics)
......@@ -16,6 +19,7 @@ import(splancs)
import(stats)
importFrom(Matrix,nearPD)
importFrom(sf,st_as_sf)
importFrom(sf,st_read)
importFrom(sf,st_write)
importFrom(utils,data)
useDynLib(landsepi)
This diff is collapsed.
......@@ -12,7 +12,7 @@
#' @param host host param
#' @return nothing yet
#' @export
modeleLandsLPI <- function(times, landscape, dispersion, inits, val_seed, host) {
invisible(.Call('_landsepi_modeleLandsLPI', PACKAGE = 'landsepi', times, landscape, dispersion, inits, val_seed, host))
modeleLandsLPI <- function(times, landscape, dispersion, inits, val_seed, host, epiP, evolP) {
invisible(.Call('landsepi_modeleLandsLPI', PACKAGE = 'landsepi', times, landscape, dispersion, inits, val_seed, host, epiP, evolP))
}
......@@ -19,30 +19,89 @@
#' @export
demo_landsepi <- function(){
propSR=5/6 ## proportion of cultivars >0
isolSR=3 ## Class of aggregation between S and R (0:3 with 1=low, 2=medium, 3=high, 4=random aggregation)
propRR=4/6 ## relative proportion of cultivars >1
isolRR=3 ## Class of aggregation between R1 and R2 (0:3 with 1=low, 2=medium, 3=high, 4=random aggregation)
strat <- "MO"
Nhote <- 3
nYears <- 10
Cmax0 <- 2
Cmax1 <- 2
#landscape parameters
propSR=5/6 ## proportion of cultivars >0
isolSR=3 ## Class of aggregation between S and R (0:3 with 1=low, 2=medium, 3=high, 4=random aggregation)
propRR=4/6 ## relative proportion of cultivars >1
isolRR=3 ## Class of aggregation between R1 and R2 (0:3 with 1=low, 2=medium, 3=high, 4=random aggregation)
strat <- "MO"
Nhote <- 3
nYears <- 30
Cmax0 <- 2
Cmax1 <- 2
data(landscapeTEST)
landscape <- AgriLand(landscapeTEST,filename="paysageTEST",propSR,isolSR,propRR,isolRR,strat,Nhote,nYears,Cmax0,Cmax1)
data(landscapeTEST)
landscape <- AgriLand(landscapeTEST,filename="paysageTEST",propSR,isolSR,propRR,isolRR,strat,Nhote,nYears,Cmax0,Cmax1)
#time parameters
nTSpY = 120
#timesStep=1:120#times step to be saved
paramT <- list(nYears=nYears,nTSpY=nTSpY)#,timesStep=timesStep)
#dispersal
fichierP <- "/home/julien/Documents/Documents/modelisation/CSIRO_INRA/landsepi/data/TEST_dispP.txt"
fichierH <- "/home/julien/Documents/Documents/modelisation/CSIRO_INRA/landsepi/data/TEST_dispH.txt"
dispP <- scan(fichierP,sep=",")
dispH <- scan(fichierH,sep=",")
#host parameters
croisH0 <- 0.10
reproH0 <- 0.0
croisH1 <- 0.10
reproH1 <- 0.0
deathH <- 0.0
RESISTANCE0 <- c(0,0,0,0,0,0,0,0)
RESISTANCE1 <- c(1,0,0,0,0,0,0,0)
RESISTANCE2 <- c(0,1,0,0,0,0,0,0)
RESISTANCE <- as.vector(cbind(RESISTANCE0,RESISTANCE1,RESISTANCE2))
Ntrait <- length(RESISTANCE)/Nhote
khost <- 0.002
sighost <- 1.001
shost <- 1.0
paramH <- list(Nhote=Nhote,croisH0=croisH0,reproH0=reproH0,croisH1=croisH1,reproH1=reproH1,deathH=deathH,Ntrait=Ntrait,resistance=as.integer(RESISTANCE),khost=khost,sighost=sighost,shost=shost)
#resultat <- st_read("paysageTEST.gpkg","results")
#habitat <- st_read("paysageTEST.gpkg","habitat")
#plot(habitat)
#plot(resultat)
#epidemiology
PSURV <- 1e-4
EFF <- 0.4
REPROP <- 3.125
TLATEXP <- 10
TLATVAR <- 9
TSPOEXP <- 24
TSPOVAR <- 104
kpatho <- 5.333
sigpatho <- 3
spatho <- 1.0
paramepi <- list(psurv=PSURV,eff=EFF,reproP=REPROP,TlatEXP=TLATEXP,TlatVAR=TLATVAR,TspoEXP=TSPOEXP,TspoVAR=TSPOVAR,kpatho=kpatho,sigpatho=sigpatho,spatho=spatho)
#evolution
COSTINFECT <- 0.75
COSTAGGR <- 0.5
TAUMUT <- 1e-7
MGEFF <- 1.0
QREFF <- 0.5
BETA <- 1.0
NAGGR <- 6
ADAPTATION <- c(1,1,0,0,0,0,0,0)
paramevol <- list(costinfect=COSTINFECT,costaggr=COSTAGGR,taumut=TAUMUT,MGeff=MGEFF,QReff=QREFF,beta=BETA,Naggr=NAGGR,adaptation=as.integer(ADAPTATION))
#run the model!
modeleLandsLPI(times = paramT,
landscape,
dispersion=list(dispP=dispP,dispH=dispH),
inits=list(C_0=0.1, PI0=5e-4),
val_seed=12345,
host=paramH,
epiP=paramepi,
evolP=paramevol)
nTSpY = 120
fichierP <- "/home/julien/Documents/Documents/modelisation/CSIRO_INRA/landsepi/data/TEST_dispP.txt"
fichierH <- "/home/julien/Documents/Documents/modelisation/CSIRO_INRA/landsepi/data/TEST_dispH.txt"
#outputs
pathRES <- getwd()
graphOn <- 1
times <- paramT
host <- paramH
epiP <- paramepi
evolP <- paramevol
HLIRdynamcis(pathRES, graphOn, times, landscape, host, epiP, evolP)
modeleLandsLPI(times = list(nYears=nYears,nTSpY=nTSpY,timesStep=1:120) , landscape, dispersion=list(dispP=fichierP,dispH=fichierH), inits=list(C_0=0.1, PI0=5e-4), val_seed=12345, host=list(Nhote=3))
}
#### FUNCTIONS TO ESTIMATE INTERVALS WITH DIFFERENT COLORS ####
logit <- function(p) log(p/(1-p))
invlogit <- function(n) exp(n) / (1+exp(n))
#### FUNCTION TO DRAW THE LANDSCAPE ####
#library("splancs") ## function polymap
#' @include RcppExports.R logit.R invlogit.R
#' @export
#'
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)
, TITLE="", SUBTITLE="", LEGEND1=rep("", length(COL.LEG)), LEGEND2=rep("", length(COL.LEG)), TITLE.LEG2="", XMAX=2000, YMAX=2000) {
......
invlogit <- function(n) exp(n) / (1+exp(n))
......@@ -53,5 +53,6 @@
#' @import fields
#' @import sp
#' @import splancs
#' @import RCALI
#' @importFrom utils data
"_PACKAGE"
logit <- function(p) log(p/(1-p))
plot.evol.QR <- function(pathRES, nIncr,trait, I_aggrProp, D, nTS,nYears,nTSpY) {
COL.grey <- gray(0:150/150)
COL.grey <- COL.grey[length(COL.grey):1]
TITLE <- "Aggressiveness trait"
LABELS <- c("\nS specialist",rep(NA,nIncr-2),"\nR specialist")
if (trait < 5) {
TITLE <- "Infectivity gene"
LABELS <- c("\nNot infective","\nInfective")
}
tiff(file=paste(pathRES,"/EVOLpatho_", names(D), ".tiff",sep=""),width=100,height=100,units='mm',compression='lzw',res=300)
par(xpd=F, mar=c(4,4,2,2))
image(x=1:nIncr, y=1:nTS, z=I_aggrProp, col=COL.grey, ylab="", xlab="Phenotype", main=paste(TITLE,names(D)), axes=F, zlim=c(0,1), ylim=c(1,nTS+1))
box()
axis(side=1, at=1:nIncr, labels=LABELS)
if (nYears==1) {
axis(2, at=round(seq(1,nTS,length.out=8)), las=1)
title(ylab="Evolutionnary time (days)")
} else {
axis(2, at=seq(1,nTS+1,nTSpY*((nYears-1)%/%8+1)), labels=seq(0,nYears,((nYears-1)%/%8+1)), las=1)
title(ylab="Evolutionnary time (years)")
}
if (!is.na(D) & D<=nTS)
abline(h=D, col="blue", lty=trait+1, lwd=2.5) ## durability of the trait
dev.off()
}
......@@ -4,7 +4,7 @@
\alias{modeleLandsLPI}
\title{Modele Lands LPI}
\usage{
modeleLandsLPI(times, landscape, dispersion, inits, val_seed, host)
modeleLandsLPI(times, landscape, dispersion, inits, val_seed, host, epiP, evolP)
}
\arguments{
\item{times}{simulation years and step}
......
......@@ -6,8 +6,8 @@
using namespace Rcpp;
// modeleLandsLPI
void modeleLandsLPI(Rcpp::List times, Rcpp::List landscape, Rcpp::List dispersion, Rcpp::List inits, int val_seed, Rcpp::List host);
RcppExport SEXP _landsepi_modeleLandsLPI(SEXP timesSEXP, SEXP landscapeSEXP, SEXP dispersionSEXP, SEXP initsSEXP, SEXP val_seedSEXP, SEXP hostSEXP) {
void modeleLandsLPI(Rcpp::List times, Rcpp::List landscape, Rcpp::List dispersion, Rcpp::List inits, int val_seed, Rcpp::List host, Rcpp::List epiP, Rcpp::List evolP);
RcppExport SEXP landsepi_modeleLandsLPI(SEXP timesSEXP, SEXP landscapeSEXP, SEXP dispersionSEXP, SEXP initsSEXP, SEXP val_seedSEXP, SEXP hostSEXP, SEXP epiPSEXP, SEXP evolPSEXP) {
BEGIN_RCPP
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< Rcpp::List >::type times(timesSEXP);
......@@ -16,17 +16,9 @@ BEGIN_RCPP
Rcpp::traits::input_parameter< Rcpp::List >::type inits(initsSEXP);
Rcpp::traits::input_parameter< int >::type val_seed(val_seedSEXP);
Rcpp::traits::input_parameter< Rcpp::List >::type host(hostSEXP);
modeleLandsLPI(times, landscape, dispersion, inits, val_seed, host);
Rcpp::traits::input_parameter< Rcpp::List >::type epiP(epiPSEXP);
Rcpp::traits::input_parameter< Rcpp::List >::type evolP(evolPSEXP);
modeleLandsLPI(times, landscape, dispersion, inits, val_seed, host, epiP, evolP);
return R_NilValue;
END_RCPP
}
static const R_CallMethodDef CallEntries[] = {
{"_landsepi_modeleLandsLPI", (DL_FUNC) &_landsepi_modeleLandsLPI, 6},
{NULL, NULL, 0}
};
RcppExport void R_init_landsepi(DllInfo *dll) {
R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
R_useDynamicSymbols(dll, FALSE);
}
......@@ -236,7 +236,7 @@ void bottleneck(const gsl_rng *gen, int Npoly, int Nhote, int Npatho, double pSu
/* HOST DYNAMIC */
/* -------------------------- */
/* Compute host reproduction, death and growth and updtate the number of H in the concerned poly */
void host_dynamic(const gsl_rng *gen, int Nhote, int Npatho, int area, double *Cmax, int *H, int *Hjuv, int **L, int **I, int **R, int *N, double *mortH, double *delta){
void host_dynamic(const gsl_rng *gen, int Nhote, int Npatho, int area, double *Cmax, int *H, int *Hjuv, int **L, int **I, int **R, int *N, double *mortH, double *delta, double khote, double sighote, double shote){
/* H, Hjuv, L, I and R are the number of host in a given poly */
int patho, hote;
double f1hote;
......@@ -310,7 +310,7 @@ void host_dynamic(const gsl_rng *gen, int Nhote, int Npatho, int area, double *C
/* CONTAMINATION : spore deposition on healthy sites */
/* -------------------------------------------------------------- */
/* Calculation of the number of contaminated sites */
void contamination(const gsl_rng *gen, int Nhote, int Npatho, int *H, int *S, int *N, int **Hcontaminated){
void contamination(const gsl_rng *gen, int Nhote, int Npatho, int *H, int *S, int *N, int **Hcontaminated, double kpatho, double sigpatho, double spatho){
/* H, S are the number of individuals in a given poly */
int patho, hote;
double f1patho;
......@@ -485,7 +485,7 @@ void infection(const gsl_rng *gen, int t, int nTSpY, int Nhote, int Npatho, int
/* DYNAMIC OF THE EPIDEMIC */
/* --------------------------- */
void dynepi(const gsl_rng *gen, int nYears, int nTSpY, int Npoly, int Nhote, int Npatho, int *area, int **habitat, int *rotation, double *C0, double pI0, double pSurv, double *Cmax, double **mort, double *delta, double *reproH, double **disphote, double **disppatho, double eff, double *Tlat, double repro, double *Tspo, double **mutkernelMG, double **mutkernelQR, int **pathoToAggr, int ********aggrToPatho, int Naggr, double **infect, double **aggr, char *strat, int **resistance, int *adaptation, int printOn, OGRLayer *poLayer, NumericVector timesStep){
void dynepi(const gsl_rng *gen, int nYears, int nTSpY, int Npoly, int Nhote, int Npatho, int *area, int **habitat, int *rotation, double *C0, double pI0, double pSurv, double *Cmax, double **mort, double *delta, double *reproH, double **disphote, double **disppatho, double eff, double *Tlat, double repro, double *Tspo, double **mutkernelMG, double **mutkernelQR, int **pathoToAggr, int ********aggrToPatho, int Naggr, double **infect, double **aggr, char *strat, int **resistance, int *adaptation, int printOn, OGRLayer *poLayer, double khote, double sighote, double shote, double kpatho, double sigpatho, double spatho){//, NumericVector timesStep,
int year, t, poly;
char name_fH[15], name_fHjuv[15], name_fS[15], name_fL[15], name_fI[15], name_fR[15];
......@@ -541,29 +541,29 @@ void dynepi(const gsl_rng *gen, int nYears, int nTSpY, int Npoly, int Nhote, int
fS = fopen(name_fS,"wb");
fR = fopen(name_fR,"wb");
int INDtime=0;
//int INDtime=0;
/* Loop for all the timesteps of the cropping season */
for (t = 0 ; t < nTSpY ; t++){
/* Writing model output for timestep t */
if(t==timesStep[INDtime]){
sortie_layer(poLayer,H,Npoly,Nhote,t,year);
INDtime=INDtime+1;
}
//write_HHjuvSLIR(Npoly, Npatho, Nhote, t, H, Hjuv, S, L, I, R, fH, fHjuv, fS, fL, fI, fR, printOn);
//if(t==timesStep[INDtime]){
//sortie_layer(poLayer,H,Npoly,Nhote,t,year);
write_HHjuvSLIR(Npoly, Npatho, Nhote, t, H, Hjuv, S, L, I, R, fH, fHjuv, fS, fL, fI, fR, printOn);
//INDtime=INDtime+1;
//}
reproDisp(gen, Npoly, Nhote, Npatho, Naggr, H, Hjuv, S, I, reproH, repro, disphote, disppatho, resistance, adaptation, mutkernelMG, mutkernelQR, pathoToAggr, aggrToPatho, aggr);
for (poly = 0 ; poly < Npoly ; poly++){
// printf("poly %d\n", poly);
host_dynamic(gen, Nhote, Npatho, area[poly], Cmax, H[poly], Hjuv[poly], L[poly], I[poly], R[poly], N, mort[poly], delta);
contamination(gen, Nhote, Npatho, H[poly], S[poly], N, Hcontaminated);
host_dynamic(gen, Nhote, Npatho, area[poly], Cmax, H[poly], Hjuv[poly], L[poly], I[poly], R[poly], N, mort[poly], delta, khote, sighote, shote);
contamination(gen, Nhote, Npatho, H[poly], S[poly], N, Hcontaminated, kpatho, sigpatho, spatho);
infection(gen, t, nTSpY, Nhote, Npatho, H[poly], Hcontaminated, L[poly], I[poly], R[poly], L2I[poly], I2R[poly], eff, Tlat, Tspo, resistance, pathoToAggr, infect, aggr);
} /* for poly */
} /* for t */
/* last time-step of the season: bottleneck before starting a new season */
sortie_layer(poLayer,H,Npoly,Nhote,nTSpY,year);
//write_HHjuvSLIR(Npoly, Npatho, Nhote, nTSpY, H, Hjuv, S, L, I, R, fH, fHjuv, fS, fL, fI, fR, printOn); /* Writing model output for last timestep */
//sortie_layer(poLayer,H,Npoly,Nhote,nTSpY,year);
write_HHjuvSLIR(Npoly, Npatho, Nhote, nTSpY, H, Hjuv, S, L, I, R, fH, fHjuv, fS, fL, fI, fR, printOn); /* Writing model output for last timestep */
bottleneck(gen, Npoly, Nhote, Npatho, pSurv, Tspo, L, I, aggr, resistance, pathoToAggr, eqIsurv, printOn); /* Calculation of the equivalent number of I */
init_HHjuvSLIR(Npoly, Nhote, Npatho, H, Hjuv, S, L, I, R); /* Re-initialisation at 0 */
init_L2I2R(Npoly, Npatho, Nhote, nTSpY, L2I, I2R); /* Re-initialisation at 0 */
......@@ -574,7 +574,7 @@ void dynepi(const gsl_rng *gen, int nYears, int nTSpY, int Npoly, int Nhote, int
intro_H(Npoly, Nhote, H, area, habitat, C0, strat, rotation[year]);
/* infection of newly planted hosts to generate the inoculum of the next season */
for (poly = 0 ; poly < Npoly ; poly++){
contamination(gen, Nhote, Npatho, H[poly], S[poly], H[poly], Hcontaminated); /* N = H[poly] in beginning of next season */
contamination(gen, Nhote, Npatho, H[poly], S[poly], H[poly], Hcontaminated, kpatho, sigpatho, spatho); /* N = H[poly] in beginning of next season */
infection(gen, 0, nTSpY, Nhote, Npatho, H[poly], Hcontaminated, L[poly], I[poly], R[poly], L2I[poly], I2R[poly], eff, Tlat, Tspo, resistance, pathoToAggr, infect, aggr);
}
......@@ -627,12 +627,12 @@ void sortie_layer(OGRLayer *poLayer,int **H,int Npoly, int Nhote, int pastemps,
/*------------------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------------------------------------------------------*/
void modeleLandsLPI(Rcpp::List times, Rcpp::List landscape, Rcpp::List dispersion, Rcpp::List inits ,int val_seed, Rcpp::List host) {
void modeleLandsLPI(Rcpp::List times, Rcpp::List landscape, Rcpp::List dispersion, Rcpp::List inits ,int val_seed, Rcpp::List host, Rcpp::List epiP, Rcpp::List evolP) {
int nYears = Rcpp::as<int>(times["nYears"]);
int nTSpY = Rcpp::as<int>(times["nTSpY"]);
NumericVector timesStep=Rcpp::as<NumericVector>(times["timesStep"]);
//NumericVector timesStep=Rcpp::as<NumericVector>(times["timesStep"]);
CharacterVector shapefilename = Rcpp::as<CharacterVector>(landscape["shapefilename"]);
CharacterVector layername_hab = Rcpp::as<CharacterVector>(landscape["layername_hab"]);
......@@ -652,13 +652,58 @@ void modeleLandsLPI(Rcpp::List times, Rcpp::List landscape, Rcpp::List dispersio
rotation[roti]=*it;
}
CharacterVector nomfile_disppatho = Rcpp::as<CharacterVector>(dispersion["dispP"]);
CharacterVector nomfile_disphote = Rcpp::as<CharacterVector>(dispersion["dispH"]);
//CharacterVector nomfile_disppatho = Rcpp::as<CharacterVector>(dispersion["dispP"]);
//CharacterVector nomfile_disphote = Rcpp::as<CharacterVector>(dispersion["dispH"]);
NumericVector dispP_tmp=Rcpp::as<NumericVector>(dispersion["dispP"]);
NumericVector dispH_tmp=Rcpp::as<NumericVector>(dispersion["dispH"]);
double C_0 = Rcpp::as<double>(inits["C_0"]);
double PI0 = Rcpp::as<double>(inits["PI0"]);
int Nhote = Rcpp::as<int>(host["Nhote"]);
double CROISH0 = Rcpp::as<double>(host["croisH0"]);
double CROISH1 = Rcpp::as<double>(host["croisH1"]);
double REPROH0 = Rcpp::as<double>(host["reproH0"]);
double REPROH1 = Rcpp::as<double>(host["reproH1"]);
double MORTH = Rcpp::as<double>(host["deathH"]);
double khote = Rcpp::as<double>(host["khost"]);
double sighote = Rcpp::as<double>(host["sighost"]);
double shote = Rcpp::as<double>(host["shost"]);
int Ntrait = Rcpp::as<int>(host["Ntrait"]);
IntegerVector resistance_tmp = Rcpp::as<IntegerVector>(host["resistance"]);
int **resistance;
resistance = init_i2(Nhote,Ntrait);
for (int ihote=0; ihote<Nhote; ihote++){
for (int itrait=0; itrait<Ntrait; itrait++){
resistance[ihote][itrait] = resistance_tmp[itrait+ihote*Ntrait];
}
}
double PSURV = Rcpp::as<double>(epiP["psurv"]);
double EFF = Rcpp::as<double>(epiP["eff"]);
double REPROP = Rcpp::as<double>(epiP["reproP"]);
double TLATEXP = Rcpp::as<double>(epiP["TlatEXP"]);
double TLATVAR = Rcpp::as<double>(epiP["TlatVAR"]);
double TSPOEXP = Rcpp::as<double>(epiP["TspoEXP"]);
double TSPOVAR = Rcpp::as<double>(epiP["TspoVAR"]);
double kpatho = Rcpp::as<double>(epiP["kpatho"]);
double sigpatho = Rcpp::as<double>(epiP["sigpatho"]);
double spatho = Rcpp::as<double>(epiP["spatho"]);
double COSTINFECT = Rcpp::as<double>(evolP["costinfect"]);
double COSTAGGR = Rcpp::as<double>(evolP["costaggr"]);
double TAUMUT = Rcpp::as<double>(evolP["taumut"]);
double MGEFF = Rcpp::as<double>(evolP["MGeff"]);
double QREFF = Rcpp::as<double>(evolP["QReff"]);
double BETA = Rcpp::as<double>(evolP["beta"]);
int NAGGR = Rcpp::as<int>(evolP["Naggr"]);
IntegerVector adaptation_tmp = Rcpp::as<IntegerVector>(evolP["adaptation"]);
int adaptation[Ntrait];
for (int itrait=0; itrait<Ntrait; itrait++){
adaptation[itrait] = adaptation_tmp[itrait];
}
int NPoly = 0;
GDALAllRegister();
......@@ -718,8 +763,6 @@ void modeleLandsLPI(Rcpp::List times, Rcpp::List landscape, Rcpp::List dispersio
gsl_rng_set(gen, seed);
/* This function initializes (or `seeds') the random number generator. */
int Nhote = Rcpp::as<int>(host["Nhote"]);
/*-----------------------*/
......@@ -787,20 +830,20 @@ void modeleLandsLPI(Rcpp::List times, Rcpp::List landscape, Rcpp::List dispersio
Tspo[1] = (double) TSPOVAR;
/* Resistance formula for each host (resistance genes that are present) */
int **resistance;
int resist_tmp[3];
//int **resistance;
//int resist_tmp[3];
resistance = init_i2(Nhote,8);
resist_tmp[0]=(int) atof(RESISTANCE0);
resist_tmp[1]=(int) atof(RESISTANCE1);
resist_tmp[2]=(int) atof(RESISTANCE2);
for (hote=0; hote<Nhote; hote++)
split_i1(8, resist_tmp[hote], resistance[hote]);
//resistance = init_i2(Nhote,8);
//resist_tmp[0]=(int) atof(RESISTANCE0);
//resist_tmp[1]=(int) atof(RESISTANCE1);
//resist_tmp[2]=(int) atof(RESISTANCE2);
//for (hote=0; hote<Nhote; hote++)
// split_i1(8, resist_tmp[hote], resistance[hote]);
/* Adaptation formula (genes that can evolve) */
int adaptation[8];
int adapt_tmp=(int) atof(ADAPTATION);
split_i1(8, adapt_tmp, adaptation);
//int adaptation[8];
//int adapt_tmp=(int) atof(ADAPTATION);
//split_i1(8, adapt_tmp, adaptation);
int nIG = adaptation[0]+adaptation[1]+adaptation[2]+adaptation[3];
int nAG = adaptation[4]+adaptation[5]+adaptation[6]+adaptation[7];
......@@ -856,28 +899,30 @@ void modeleLandsLPI(Rcpp::List times, Rcpp::List landscape, Rcpp::List dispersio
/*-----------*/
/* Dispersal */
/*-----------*/
int N_elements_dispTot = NPoly*NPoly;
//int N_elements_dispTot = NPoly*NPoly;
//char *nomfile_disppatho = "DISPPATHO";
//char *nomfile_disphote = "DISPHOTE";
double *y_disppathoTot;
//double *y_disppathoTot;
double **disppathoTot;
double **disppatho;
double *y_disphoteTot;
//double *y_disphoteTot;
double **disphoteTot;
double **disphote;
y_disppathoTot = init_d1(N_elements_dispTot);
//y_disppathoTot = init_d1(N_elements_dispTot);
disppathoTot = init_d2(NPoly,NPoly);
disppatho = init_d2(NPoly,NPoly);
y_disphoteTot = init_d1(N_elements_dispTot);
//y_disphoteTot = init_d1(N_elements_dispTot);
disphoteTot = init_d2(NPoly,NPoly);
disphote = init_d2(NPoly,NPoly);
lire(120, N_elements_dispTot, (char*)Rcpp::as<std::string>(nomfile_disppatho).c_str(), pdelim, y_disppathoTot);
lire(120, N_elements_dispTot, (char*)Rcpp::as<std::string>(nomfile_disphote).c_str(), pdelim, y_disphoteTot);
//lire(120, N_elements_dispTot, (char*)Rcpp::as<std::string>(nomfile_disppatho).c_str(), pdelim, y_disppathoTot);
//lire(120, N_elements_dispTot, (char*)Rcpp::as<std::string>(nomfile_disphote).c_str(), pdelim, y_disphoteTot);
for (i=0 ; i<NPoly ; i++){
for (j=0 ; j<NPoly ; j++) {
disppathoTot[j][i] = y_disppathoTot[j+i*NPoly];
disphoteTot[j][i] = y_disphoteTot[j+i*NPoly];
//disppathoTot[j][i] = y_disppathoTot[j+i*NPoly];
//disphoteTot[j][i] = y_disphoteTot[j+i*NPoly];
disppathoTot[j][i] = dispP_tmp[j+i*NPoly];
disphoteTot[j][i] = dispP_tmp[j+i*NPoly];
}
}
......@@ -886,7 +931,7 @@ void modeleLandsLPI(Rcpp::List times, Rcpp::List landscape, Rcpp::List dispersio
/*--------------------------------------*/
FILE *fP;
fP = fopen("parameters.txt","w");
print_param(fP, seed, (char*)Rcpp::as<std::string>(nomfile_disppatho).c_str(), (char*)Rcpp::as<std::string>(nomfile_disphote).c_str(), "nomfile_habitat1", "nomfile_habitat2", nYears, nTSpY, pSurv, eff, Tlat, repro, Tspo, NPoly, NPoly, Nhote, Npatho, Naggr, C0, pI0, Cmax, delta, reproH, area, habitat, rotation, mort, strat, propRR, resistance, adaptation, MGeff, QReff, taumut, costInfect, costAggr, beta);
//print_param(fP, seed, (char*)Rcpp::as<std::string>(nomfile_disppatho).c_str(), (char*)Rcpp::as<std::string>(nomfile_disphote).c_str(), "nomfile_habitat1", "nomfile_habitat2", nYears, nTSpY, pSurv, eff, Tlat, repro, Tspo, NPoly, NPoly, Nhote, Npatho, Naggr, C0, pI0, Cmax, delta, reproH, area, habitat, rotation, mort, strat, propRR, resistance, adaptation, MGeff, QReff, taumut, costInfect, costAggr, beta);
//print_param(stdout, seed, (char*)Rcpp::as<std::string>(nomfile_disppatho).c_str(), (char*)Rcpp::as<std::string>(nomfile_disphote).c_str(), "nomfile_habitat1", "nomfile_habitat2", nYears, nTSpY, pSurv, eff, Tlat, repro, Tspo, NPoly, NPoly, Nhote, Npatho, Naggr, C0, pI0, Cmax, delta, reproH, area, habitat, rotation, mort, strat, propRR, resistance, adaptation, MGeff, QReff, taumut, costInfect, costAggr, beta);
fclose(fP);
......@@ -931,7 +976,7 @@ void modeleLandsLPI(Rcpp::List times, Rcpp::List landscape, Rcpp::List dispersio
/* -------------- */
Rprintf("\n*** SPATIOTEMPORAL MODEL SIMULATING THE SPREAD OF AN EPIDEMIC IN A LANDSCAPE ***\n\n");
poLayer = poDS->GetLayerByName( Rcpp::as<std::string>(layername_res).c_str());
dynepi(gen, nYears, nTSpY, NPoly, Nhote, Npatho, area, habitat, rotation, C0, pI0, pSurv, Cmax, mort, delta, reproH, disphote, disppatho, eff, Tlat, repro, Tspo, mutkernelMG, mutkernelQR, pathoToAggr, aggrToPatho, Naggr, infect, aggr, strat, resistance, adaptation, printOn, poLayer, timesStep);
dynepi(gen, nYears, nTSpY, NPoly, Nhote, Npatho, area, habitat, rotation, C0, pI0, pSurv, Cmax, mort, delta, reproH, disphote, disppatho, eff, Tlat, repro, Tspo, mutkernelMG, mutkernelQR, pathoToAggr, aggrToPatho, Naggr, infect, aggr, strat, resistance, adaptation, printOn, poLayer, khote, sighote, shote, kpatho, sigpatho, spatho);//timesStep
/* ----------- */
/* Free memory */
......@@ -953,10 +998,10 @@ void modeleLandsLPI(Rcpp::List times, Rcpp::List landscape, Rcpp::List dispersio
// free(y_rotation);
free(rotation);
free(y_disphoteTot);
free(y_disppathoTot);
// free(y_disphoteTot);
// free(y_disppathoTot);
// free(y_areaTot);
//free(y_habitat1Tot);
// free(y_habitat1Tot);
// free(y_habitat2Tot);
// free_i2(habitatTot,2);
free_i2(habitat,2);
......
......@@ -23,12 +23,7 @@
#include "printReadWrite.hpp" /* Printing, reading, writing functions */
#include "initialisation.hpp"
#define kpatho 5.333
#define sigpatho 3.0
#define spatho 1.0
#define khote 0.002
#define sighote 1.001
#define shote 1.0
#define ID_E 4 /* index of QR on infection efficiency */
#define ID_L 5 /* index of QR on latent period duration */
#define ID_R 6 /* index of QR on reproduction rate */
......@@ -36,32 +31,40 @@
#define PRINTON 0
#define NAGGR 6
//#define STRAT "RO"
//#define PROPRR 4/6
#define COSTINFECT 0.75
#define COSTAGGR 0.5
#define TAUMUT 1e-7
#define MGEFF 1.0
#define QREFF 0.5
#define BETA 1.0
#define CROISH0 0.10
#define REPROH0 0.0
#define CROISH1 0.10
#define REPROH1 0.0
#define MORTH 0.0
#define PSURV 1e-4
#define EFF 0.4
#define REPROP 3.125
#define TLATEXP 10
#define TLATVAR 9
#define TSPOEXP 24
#define TSPOVAR 104
#define RESISTANCE0 "00000000"
#define RESISTANCE1 "10000000"
#define RESISTANCE2 "01000000"
#define ADAPTATION "11000000"
//#define PSURV 1e-4
//#define EFF 0.4
//#define REPROP 3.125