R/methodDefinitions.R

In rocheben/ibmR: Individual-Based Model in R for phylodynamics studies

#' 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();
  }
}