#' Adding a population to the model
#' @param modelObject the scheduler object for which the populations has to be added
#' @param nbIndividuals Number of individuals of this population that has to be created
#' @param evtCtctRate Contact rate with the environment
#' @param intraSpecCtctRate Intra-specific contact rate within this population
#' @param lifespan Average lifespan of indivuduals within this population
#' @param name (optional) Name of the population
#' @return an object of class scheduler for update purpose
#' @export
addPopulation<-function(modelObject,nbIndividuals,evtCtctRate,intraSpecCtctRate,lifespan,name=NULL){
temp=attr(modelObject,"populations");
if(length(name)==0){
name=sprintf("pop%d",length(temp)+1);
}
pop=new("population",evtCtctRate=evtCtctRate,intraSpecCtctRate=intraSpecCtctRate,name=name,id=length(temp)+1);
temp=c(temp,pop);
attr(modelObject,"populations")=temp;
for(i in 1:nbIndividuals){
ind=new("individual",popId=pop,lifespan=lifespan)
temp=attr(modelObject,"individuals");
temp=c(temp,ind);
attr(modelObject,"individuals")=temp
}
print(sprintf("%d individuals of population '%s' created",nbIndividuals,attr(pop,"name")))
return(modelObject)
}
#' Adding an environment to the model
#' @param modelObject the scheduler object for which the populations has to be added
#' @param seasonality time series of fluctuating environment
#' @param name (optional) Name of the population
#' @return an object of class scheduler for update purpose
#' @export
addEnvironment<-function(modelObject,seasonality,name=NULL){
temp=attr(modelObject,"evt");
if(length(name)==0){
name=sprintf("evt%d",length(temp)+1);
}
evt=new("evt",seasonality=data.frame(seasonality),name=name,pathogens=list(),pathogenLoad=list(),id=length(temp)+1);
temp=c(temp,evt);
attr(modelObject,"environment")=temp;
print(sprintf("New environment created: '%s' created",attr(evt,"name")))
return(modelObject)
}
#' Adding a pathogen to the model
#' @param modelObject the scheduler object for which the populations has to be added
#' @param mutationRate Mutation rate of the pathogen for each population
#' @param crossProtection Factor for intra-specific pathogen cross protection
#' @param probInfection Probability of infection upon contact for each population
#' @param latencyPeriod Latency period of the pathogen for each population
#' @param infectiousPeriod Infectious period of the pathogen for each population
#' @param recoveryPeriod Immunity period of the pathogen for each population
#' @param name (optional) Name of the pathogen
#' @return an object of class scheduler for update purpose
#' @export
addPathogen<-function(modelObject,mutationRate,crossProtection,probaInfection,latencyPeriod,infectiousPeriod,recoveryPeriod,nbInfectedtoAdd=1,name=NULL){
temp=attr(modelObject,"pathogens");
if(length(name)==0){
name=sprintf("path%d",length(temp)+1);
}
path=new("pathogen",mutationRate=mutationRate,crossProtection=crossProtection,probaInfection=probaInfection,latencyPeriod=latencyPeriod,infectiousPeriod=infectiousPeriod,recoveryPeriod=recoveryPeriod,name=name);
temp=c(temp,path);
attr(modelObject,"pathogens")=temp;
attr(modelObject,"nbInfected")=c(attr(modelObject,"nbInfected"),nbInfectedtoAdd);
print(sprintf("New pathogen created: '%s' created",attr(path,"name")))
return(modelObject)
}
#' Setting the interspecific protection matrix
#' @param modelObject the scheduler object hosting the populations
#' @param matrixContact Matrix of protection among pathogens
#' @return an object of class scheduler for update purpose
#' @export
setMatrixContact<-function(modelObject,matrixContact){
if(dim(matrixContact)[1]!=length(attr(modelObject,"populations"))){
stop('Number of pathogens defined and matrix size don\'t match. All populations have to be defined first');
}
else{
attr(modelObject,"matrixContact")=matrixContact;
}
return(modelObject);
attr(modelObject,"interSpecificProtection")=interSpecificProtection;
return(modelObject);
}
#' Simulate the Individual-Based Model
#' @param modelObject the scheduler object for which the populations has to be added
#' @param numberInfected Number of infected individuals within each population
#' @param timeMax Maximal number of time steps for the simulation
#' @return an object of class scheduler for update purpose
#' @export
simulateModel<-function(modelObject,times){
for(t in times){
print(t)
arrayPop=attr(modelObject,"populations");
nbInfected=attr(modelObject,"numberInfected");
arrayInd=attr(modelObject,"individuals");
arrayEvt=attr(modelObject,"environment");
for(i in 1:length(arrayInd)){
arrayInd[i]=step(modelObject,arrayPop,arrayInd[i]);
}
for(i in 1:length(environment)){
#arrayEvt[i]=stepEvt(arrayEvt[i])
}
for(i in 1:length(arrayInd)){
#arrayInd[i]=update(modelObject,arrayPop,arrayInd[i]);
}
}
}
#####################################################################
# INTERNAL FUNCTIONS #
#####################################################################
###Infection function. Not accessible from outside
infection<-function(model,pPathogen){
}
###Mutation function. Not accessible from outside
mutation<-function(pPathogen){
}
###Step function. Not accessible from outside
stepEvt<-function(modelObject,arrayPop,indObject,evtObj){
# Get the ID of the individual population
idPop=attr(attr(indObject,"population"),"id")
# Get the pathogens that have already infected the individual
pathInd=attr(indObject,"pathogens");
# Get the number of infectious individuals for each pathogens
nbInfected=attr(modelObject,"nbInfected");
matrixContact=attr(modelObject,"matrixContact");
interSpecificProtection=attr(modelObject,"interSpecificProtection");
# Get the pathogen objects
arrayPath=attr(modelObject,"pathogens");
indexPat=sample(1:length(arrayPath))
for(k in length(arrayPath)){
#Infection from Environment
loadEvt=attr(evtObj,"pathogenLoad");
probaInfection=evtCtctRate*(loadEvt/(loadEvt+patArray[k]$loadNeeded))
proba=modulationProbaInfectionHistory(proba,patArray[k]);
if(rand()<proba){
#Infection
infection(arrayPath[[indexPat[k]]]);
print("Infection Evt");
}
}
}
###Step function. Not accessible from outside
step<-function(modelObject,arrayPop,indObject,evtObj){
# Get the ID of the individual population
idPop=attr(attr(indObject,"population"),"id")
# Get the pathogens that have already infected the individual
pathInd=attr(indObject,"pathogens");
# Get the number of infectious individuals for each pathogens
nbInfected=attr(modelObject,"nbInfected");
matrixContact=attr(modelObject,"matrixContact");
interSpecificProtection=attr(modelObject,"interSpecificProtection");
# Get the pathogen objects
arrayPath=attr(modelObject,"pathogens");
indexPat=sample(1:length(arrayPath))
for(k in length(arrayPath)){
#Determine force of infection
proba=attributes(arrayPath[[indexPat[k]]])$probaInfection*sum(matrixContact*as.numeric(nbInfected));
#proba=modulationProbaInfectionHistory(proba,patArray[k]);
if( runif(1)<proba){
#Infection
infection(arrayPath[[indexPat[k]]]);
print("Infection Direct");
}
}
}
#Update function. Not accessible from outside
update<-function(){
if(runif(1)<mutationRate){
#mutation();
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.