The Source File for the Action Class Virus

At the beginning of Virus.cpp we must put the header files we need:

#include <omp.h>

#include "MessLoggerT.h"

#include "clsutils.cpp"
#include "ParamProvider2.h"
#include "SPopulation.h"
#include "SCell.h"
#include "WELL512.h"
#include "QDFUtils.h"
#include "QDFUtilsT.h"
#include "Virus.h"

We also need to define attribute names to read and write the action’s attributes.

template<typename T>
const char *Virus<T>::asNames[] = {
    ATTR_VIRUS_INFECTION_PROB_NAME,
    ATTR_VIRUS_INITIAL_LOAD_NAME,
    ATTR_VIRUS_GROWTH_RATE_NAME,
    ATTR_VIRUS_CONTAGION_LEVEL_NAME,
    ATTR_VIRUS_LETHALITY_LEVEL_NAME,
};

constructor

There is not much to do in the constructor: intialisation of the variables, and saving the attribute names.

template<typename T>
Virus<T>::Virus(SPopulation<T> *pPop, SCellGrid *pCG, std::string sID, WELL512 **apWELL)
    : Action<T>(pPop,pCG,ATTR_VIRUS_NAME,sID),
    m_apWELL(apWELL),
    m_fInfectionProb(0),
    m_fInitialLoad(0),
    m_fGrowthRate(0),
    m_fContagionLevel(0),
    m_fLethalityLevel(0),
    m_iNumThreads(omp_get_max_threads()) {

    this->m_vNames.insert(this->m_vNames.end(), asNames, asNames+sizeof(asNames)/sizeof(char*));

    // create array with one int map for each thread
    m_amCellAgents = new cellagentmap[m_iNumThreads];

    // create array with one map for each thread
    m_amAgentInfects = new  agentloadmap[m_iNumThreads];

}

destructor

The destructor deletes the arrays that were allocated in the constructor

template<typename T>
Virus<T>::~Virus() {
    // delete the arrays
    delete[] m_amCellAgents;
    delete[] m_amAgentInfects;
}

initialize

The method initialize is called before execute(). Here the array of maps relating cell indexes to agent indexes (m_amCellAGents) is filled, and the array of maps relating agent IDs to virus loads (m_amAgentInfects) is resaet.

Note that since initialize is not called in a parllelized context, we can use parallelized loops with the full set of available threads.

template<typename T>
int Virus<T>::initialize(float fT) {
    int iResult = 0;

    // clear the agent arrays and the infection registerr
    for (int iT = 0; iT < m_iNumThreads; iT++) {
        m_amCellAgents[iT].clear();
        m_amAgentInfects[iT].clear();
    }

    int iFirstAgent = this->m_pPop->getFirstAgentIndex();
    int iLastAgent  = this->m_pPop->getLastAgentIndex();

#pragma omp parallel for
    for (int iA = iFirstAgent; iA <= iLastAgent; iA++) {
        int iT = omp_get_thread_num();
        T* pA = &(this->m_pPop->m_aAgents[iA]);
        m_amCellAgents[iT][pA->m_iCellIndex]->push_back(iA);
    }

    return iResult;
}

execute

The execute() method is applied to every agent (in parallel). This method is usually the one most of the action’s work is done. Here is what we need the virus to do for each agent: - let the the agent’s virus load grow logistically - if the agent’s virus load exceeds the lethality level, register it for death - if the agent’s virus load exceeds the contagion level, it infects all agents in the cell

template<typename T>
int Virus<T>::execute(int iAgentIndex, float fT) {
    int iResult = 0;

    int iT = omp_get_thread_num();

    T *pa = &(this->m_pPop->m_aAgents[iAgentIndex]);
    // we only handle agents not already marked as dead
    if (pa->m_iLifeState > 0) {

        // perform the logistic growth of the agent's viral load
        float fM = pa->m_fViralLoad;
        pa->m_fViralLoad += m_fGrowthRate*fM*(1-fM);

        int iCellIndex = pa->m_iCellIndex;

        // check if the agent must die
        if (pa->m_fViralLoad > m_fLethalityLevel) {
            this->m_pPop->registerDeath(iCellIndex, iAgentIndex);

            // is the agent contagious?
        } else if  (pa->m_fViralLoad > m_fContagionLevel) {
            // loop through agents in cell
            for (uint i = 0; i < m_amCellAgents[iT][iCellIndex].size(); i++) {
                int iOtherAgent = m_amCellAgents[iCellIndex][i];
                // we only infect other agents
                if (iOtherAgent != iAgentIndex) {
                    T *pao = &(this->m_pPop->m_aAgents[iOtherAgent]);
                    // infect!
                    double dR =  this->m_apWELL[omp_get_thread_num()]->wrandd();
                    // we reduce the effective infection probability by multiplying it with (1-immunity)
                    if (dR < m_fInfectionProb * (1-pao->m_fImmunity)) {
                        // register infection load for agent
                        m_amAgentInfects[iT][iOtherAgent] += m_fInitialLoad;
                    }
                }
            }
        }
    }
    return iResult;
}

finalize

The method finalize is called after all calls to execute()`. Here we finally add to each agent’s virus load the colllected virus load passed from its neighbors. During execute() several maps may have entries for the same agent.

template<typename T>
int Virus<T>::finalize(float fT) {
    int iResult = 0;

    for (int iT = 0; iT < m_iNumThreads; iT++) {
        for (agentloadmap::const_iterator it = m_amAgentInfects[iT].begin(); it !=  m_amAgentInfects[iT].end(); ++it) {
            // first:  agent index
            // second: infecttload

            T *pa = &(this->m_pPop->m_aAgents[it->first]);
            pa->m_fViralLoad += it->second;
        }
    }
    return iResult;
}

The remaining methods deal with reading and writing the virus’ attributes.

extractAttributesQDF

This method reads virus attributes from a QDF file.

template<typename T>
int Virus<T>::extractAttributesQDF(hid_t hSpeciesGroup) {

    int iResult = 0;

    if (iResult == 0) {
        iResult = qdf_extractAttribute(hSpeciesGroup, ATTR_VIRUS_INFECTION_PROB_NAME, 1, &m_fInfectionProb);
        if (iResult != 0) {
            LOG_ERROR("[Virus] couldn't read attribute [%s]", ATTR_VIRUS_INFECTION_PROB_NAME);
        }
    }

    if (iResult == 0) {
        iResult = qdf_extractAttribute(hSpeciesGroup, ATTR_VIRUS_INITIAL_LOAD_NAME, 1, &m_fInitialLoad);
        if (iResult != 0) {
            LOG_ERROR("[Virus] couldn't read attribute [%s]", ATTR_VIRUS_INITIAL_LOAD_NAME);
        }
    }

    if (iResult == 0) {
        iResult = qdf_extractAttribute(hSpeciesGroup, ATTR_VIRUS_GROWTH_RATE_NAME, 1, &m_fGrowthRate);
        if (iResult != 0) {
            LOG_ERROR("[Virus] couldn't read attribute [%s]", ATTR_VIRUS_GROWTH_RATE_NAME);
        }
    }

    if (iResult == 0) {
        iResult = qdf_extractAttribute(hSpeciesGroup, ATTR_VIRUS_CONTAGION_LEVEL_NAME, 1, &m_fContagionLevel);
        if (iResult != 0) {
            LOG_ERROR("[Virus] couldn't read attribute [%s]", ATTR_VIRUS_CONTAGION_LEVEL_NAME);
        }
    }

    if (iResult == 0) {
        iResult = qdf_extractAttribute(hSpeciesGroup, ATTR_VIRUS_LETHALITY_LEVEL_NAME, 1, &m_fLethalityLevel);
        if (iResult != 0) {
            LOG_ERROR("[Virus] couldn't read attribute [%s]", ATTR_VIRUS_LETHALITY_LEVEL_NAME);
        }
    }
    return iResult;
}

writeAttributesQDF

This method writes the virus’ attributes to a QDF file.

template<typename T>
int Virus<T>:: writeAttributesQDF(hid_t hSpeciesGroup) {
    int iResult = 0;

    iResult += qdf_insertAttribute(hSpeciesGroup, ATTR_VIRUS_INFECTION_PROB_NAME, 1, &m_fInfectionProb);
    iResult += qdf_insertAttribute(hSpeciesGroup, ATTR_VIRUS_INITIAL_LOAD_NAME,  1, &m_fInitialLoad);
    iResult += qdf_insertAttribute(hSpeciesGroup, ATTR_VIRUS_GROWTH_RATE_NAME,  1, &m_fGrowthRate);
    iResult += qdf_insertAttribute(hSpeciesGroup, ATTR_VIRUS_CONTAGION_LEVEL_NAME,  1, &m_fContagionLevel);
    iResult += qdf_insertAttribute(hSpeciesGroup, ATTR_VIRUS_LETHALITY_LEVEL_NAME,  1, &m_fLethalityLevel);

    return iResult;
}

tryGetAttributes

This method reads the virus attributes passed via the xml population file<pop_xml_ref>.

template<typename T>
int Virus<T>::tryGetAttributes(const ModuleComplex *pMC) {
    int iResult = -1;

    const stringmap &mParams = pMC->getAttributes();
    if (this->checkAttributes(mParams) == 0) {
        iResult = 0;
        iResult += getAttributeVal(mParams,  ATTR_VIRUS_INFECTION_PROB_NAME, &m_fInfectionProb);
        iResult += getAttributeVal(mParams,  ATTR_VIRUS_INITIAL_LOAD_NAME,  &m_fInitialLoad);
        iResult += getAttributeVal(mParams,  ATTR_VIRUS_GROWTH_RATE_NAME,  &m_fGrowthRate);
        iResult += getAttributeVal(mParams,  ATTR_VIRUS_CONTAGION_LEVEL_NAME,  &m_fContagionLevel);
        iResult += getAttributeVal(mParams,  ATTR_VIRUS_LETHALITY_LEVEL_NAME,  &m_fLethalityLevel);
    }
    return iResult;
}

The entire file can be found here: Virus.cpp

Implementation of the population class VirusHostPop

Since our Virus` action makes certain assumptions about the agents, we can’t use an arbitary population class. In particular, it needs fields for the virus load and immunity.

The header file

First we have the includes for the header files of all required actions (among them Virus.h), and the header file for the population’s base class SPopulation

#include "GetOld.h"
#include "ATanDeath.h"
#include "RandomMove.h"
#include "Fertility.h"
#include "Verhulst.h"
#include "RandomPair.h"
#include "Virus.h"
#include "SPopulation.h"

Some constants needed to extract mutation rate and inheritance types

const std::string VAR_VIRUSHOST_MUT_RATE_NAME    = "MutationRate";
const std::string VAR_VIRUSHOST_IMM_INHERIT_NAME = "ImmunityInheritance";
const std::string INH_TYPES[] = {
    "mix",
    "mat",
    "pat",
    "min",
    "max",
    };

const int INH_MIX     = 0;
const int INH_MAT     = 1;
const int INH_PAT     = 2;
const int INH_MIN     = 3;
const int INH_MAX     = 4;

Now the declaration of the agent structure: we have members for agent fertility, and the members for the viral infections virus load, immunity and incoming virus load.

struct VirusHostAgent : Agent {

    float m_fAge;
    float m_fLastBirth;
    int m_iMateIndex;

    float m_fViralLoad;
    float m_fImmunity;
};

The declaration of the VirusHostPop class derived from SPopulation:

class VirusHostPop : public  SPopulation<VirusHostAgent> {

The public methods are the usual: constructor, destructor and methods to add agent data from a populatiom dat file<pop_dat_ref> and to write agent data to a QDF file. In addition to that there is a method for setting the inital state of the baby (which we’ve also encountered before), and a method to get a parameter for the population from the populaition xml file<pop_xml_ref> (in this case we need the mutation rate, and the inheritance type).

public:
    VirusHostPop(SCellGrid *pCG, PopFinder *pPopFinder, int iLayerSize, IDGen **apIDG, uint32_t *aulState, uint *aiSeeds);
    virtual ~VirusHostPop_SexualPop();

    int addPopSpecificAgentData(int iAgentIndex, char **ppData);
    void addPopSpecificAgentDataTypeQDF(hid_t *hAgentDataType);

    int makePopSpecificOffspring(int iAgent, int iMother, int iFather);
    int getPopParams(const stringmap &mVarDefs);

In the protected section there are pointers to the actions of this population and the variable for the mutation rate.

protected:
    GetOld<VirusHostAgent>      *m_pGO;
    ATanDeath<VirusHostAgent>   *m_pAD;
    RandomMove<VirusHostAgent>  *m_pRM;
    Fertility<VirusHostAgen>    *m_pFert;
    Verhulst<VirusHostAgen>     *m_pVerhulst;
    RandomPair<VirusHostAgen>   *m_pPair;
    Virus<VirusHostAgen>        *m_pVirus;

    float m_fMutationRate;

constructor

Here we instantiate all actions and initialize the fields.

VirusHostPop::VirusHostPop(SCellGrid *pCG, PopFinder *pPopFinder, int iLayerSize, IDGen **apIDG, uint32_t *aulState, uint *aiSeeds)
    : SPopulation<VirusHostAgent>(pCG, pPopFinder, iLayerSize, apIDG, aulState, aiSeeds) {


    m_pGO  = new GetOld<VirusHostAgent>(this, m_pCG, "");
    m_pAD  = new ATanDeath<VirusHostAgent>(this, m_pCG, "", m_apWELL);
    m_pRM  = new RandomMove<VirusHostAgent>(this, m_pCG, "", m_apWELL);
    m_pFert     = new Fertility<VirusHostAgent>(this, m_pCG, "");
    m_pVerhulst = new Verhulst<VirusHostAgent>(this, m_pCG, "", m_apWELL);

    m_pPair = new RandomPair<VirusHostAgent>(this, m_pCG, "", m_apWELL);
    m_pVirus = new Virus<VirusHostAgent>(this, m_pCG, "", m_apWELL);

    m_prio.addAction(m_pGO);
    m_prio.addAction(m_pAD);
    m_prio.addAction(m_pRM);
    m_prio.addAction(m_pFert);
    m_prio.addAction(m_pVerhulst);

    m_prio.addAction(m_pPair);
    m_prio.addAction(m_pVirus);

    m_fMutationRate = 0.0;
    m_sInheritType = "mix";
}

destructor

The destructor deletes all action created in the constructor.

VirusHostPop::~VirusHostPop() {

    if (m_pGO != NULL) {
        delete m_pGO;
    }
    if (m_pAD != NULL) {
        delete m_pAD;
    }
    if (m_pRM != NULL) {
        delete m_pRM;
    }
    if (m_pFert != NULL) {
        delete m_pFert;
    }
    if (m_pVerhulst != NULL) {
        delete m_pVerhulst;
    }
    if (m_pPair != NULL) {
        delete m_pPair;
    }

    if (m_pVirus != NULL) {
        delete m_pVirus;
    }

}

getPopParams

This method tries to extract mutation rate and inheritance type from the stringmap created by the population XML reader.

int  VirusHostPop::getPopParams(const stringmap &mVarDefs) {
    int iResult = -1;

    stringmap::const_iterator it = mVarDefs.find(VAR_VIRUSHOST_MUT_RATE_NAME);
    if (it != mVarDefs.end()) {
        float fMutTemp = 0;
        if (strToNum(it->second, &fMutTemp)) {
            m_fMutationRate = fMutTemp;
            printf("[VirusHostPop::getPopParams] have mutation rate %f\n", m_fMutationRate);
            iResult = 0;
        } else {
            // couldn't convert string to float
        }
    } else {
        // there is no key VAR_VIRUSHOST_MUT_RATE_NAME
    }

    it = mVarDefs.find(VAR_VIRUSHOST_IMM_INHERIT_NAME);
    if (it != mVarDefs.end()) {
        std::string sInheritType = it->second;

        // handle the "+"
        if (!sInheritType.ends_with("+")) {
            m_fMutationRate = 0;
        } else {
            // drop the plus
            sInheritType = sInheritType.substr(0, sInheritType.size()-1);
        }
        // check if valid inherit type
        bool bSearching = true;
        for (int i  = 0; bSearching && (i < 5); i++) {
            if (sInheritType == INH_TYPES[i]) {
                bSearching = false;
            }
        }

        if (bSearching) {
            printf("[VirusHostPop::getPopParams] unknown inherit type [%s]\n", sInheritType.c_str());
        } else {
            m_sInheritType = sInheritType;
            printf("[VirusHostPop::getPopParams] have inherit type [%s]\n", m_sInheritType.c_str());
            iResult = 0;
        }

    } else {
        // there is no key VAR_VIRUSHOST_IMM_INHERIT_NAME
    }
    return iResult;
}

makePopSpecificOffspring

This method is called when a baby is born in order to initialize it’s field. In particular, the baby’s immunity level is determined by the inhertiance type.

int VirusHostPop::makePopSpecificOffspring(int iAgent, int iMother, int iFather) {
    int iResult = 0;

    m_aAgents[iAgent].m_fAge = 0.0;
    m_aAgents[iAgent].m_fLastBirth = 0.0;
    m_aAgents[iAgent].m_iMateIndex = -3;
    // SPopulation assigns random genders to a baby, ehich is ok here
    m_aAgents[iAgent].m_fViralLoad = 0;


    // register this date as the last birth of the mother
    m_aAgents[iAgent].m_fLastBirth = m_fCurTime;

    // the offspring's immunity is the average of the parents' immunities plus a little fuzz factor.
    float fImm = 0.0;
    if (m_sInheritType == INH_TYPES[INH_MIX]) {
        fImm =  (m_aAgents[iMother].m_fImmunity + m_aAgents[iFather].m_fImmunity)/2;
    } else if (m_sInheritType == INH_TYPES[INH_MAT]) {
        fImm =  m_aAgents[iMother].m_fImmunity;
    } else if (m_sInheritType == INH_TYPES[INH_PAT]) {
        fImm =  m_aAgents[iFather].m_fImmunity;
    } else if (m_sInheritType == INH_TYPES[INH_MIN]) {
        fImm = (m_aAgents[iMother].m_fImmunity < m_aAgents[iFather].m_fImmunity)?m_aAgents[iMother].m_fImmunity:m_aAgents[iFather].m_fImmunity;
    } else if (m_sInheritType == INH_TYPES[INH_MAX]) {
        fImm = (m_aAgents[iMother].m_fImmunity > m_aAgents[iFather].m_fImmunity)?m_aAgents[iMother].m_fImmunity:m_aAgents[iFather].m_fImmunity;
    } else {
        // shouldn't happen
        iResult = -1;
    }

    //mutate
    fImm += m_apWELL[omp_get_thread_num()]->wgauss(m_fMutationRate);
    if (fImm < 0) {
        fImm = 0;
    }else if (fImm > 1) {
        fImm = 1;
    }
    m_aAgents[iAgent].m_fImmunity = fImm;


    if (m_aAgents[iMother].m_fViralLoad > 0) {
        m_aAgents[iAgent].m_fViralLoad = 0.2;
        printf("Babyinfect: %d -> %d\nn", iMother, iAgent);
    }


    //printf("[VirusHostPop::makePopSpecificOffsprin] baby %d vl %f, imm %f\n",  iAgent, m_aAgents[iAgent].m_fViralLoad, m_aAgents[iAgent].m_fImmunity);
    return iResult;
}

addPopSpecificAgentData

This method reads agent specific data from the population dat file<pop_dat_ref>.

int VirusHostPop::addPopSpecificAgentData(int iAgentIndex, char **ppData) {

    int iResult = 0;

    if (iResult == 0) {
        iResult = addAgentDataSingle(ppData, &m_aAgents[iAgentIndex].m_fAge);
    }

    if (iResult == 0) {
        iResult = addAgentDataSingle(ppData, &m_aAgents[iAgentIndex].m_fLastBirth);
    }

    if (iResult == 0) {
        iResult = addAgentDataSingle(ppData, &m_aAgents[iAgentIndex].m_fViralLoad);
    }

    if (iResult == 0) {
        iResult = addAgentDataSingle(ppData, &m_aAgents[iAgentIndex].m_fImmunity);
    }

    return iResult;
}

addPopSpecificAgentDataTypeQDF

This method creates a HDF data type to handlle I/O of agents.

void VirusHostPop::addPopSpecificAgentDataTypeQDF(hid_t *hAgentDataType) {

    VirusHostAgent va;
    H5Tinsert(*hAgentDataType, "Age",       qoffsetof(va, m_fAge), H5T_NATIVE_FLOAT);
    H5Tinsert(*hAgentDataType, "LastBirth", qoffsetof(va, m_fLastBirth), H5T_NATIVE_FLOAT);
    H5Tinsert(*hAgentDataType, "ViralLoad", qoffsetof(va, m_fViralLoad), H5T_NATIVE_FLOAT);
    H5Tinsert(*hAgentDataType, "Immunity",  qoffsetof(va, m_fImmunity), H5T_NATIVE_FLOAT);
}