TheoreticalEcology / EcoPhyloSim

R package for simulation of biogeographical and phylogenetic data
http://theoreticalecology.github.io/EcoPhyloSim/
0 stars 0 forks source link

Species m_Mean updated but not m_Variance #79

Closed TankredO closed 7 years ago

TankredO commented 7 years ago
void Species::addIndividual(double env, double comp, double neutral) {
    m_Count += 1;
    m_MeanSum += env;
    m_CompetitionSum += comp;
    m_NeutralSum += neutral;
    updateMean(); // why is variance not updated?
}

void Species::removeIndividual(double env, double comp, double neutral, int generation) {

    if (m_Count < 2) {
        m_Date_of_Extinction = generation;
        m_Count -= 1;
    } else {
        m_Count -= 1;
        m_MeanSum -= env;
        m_CompetitionSum -= comp;
        m_NeutralSum -= neutral;
        updateMean(); // why is variance not updated?
    }
}

void Species::updateMean() {
    m_Mean = m_MeanSum / (double(m_Count));
    m_CompetitionMean = m_CompetitionSum / (double(m_Count));
    m_NeutralMean = m_NeutralSum / (double(m_Count));
}

When new individuals are added to or removed from a species, the trait means are updated, but the variance is not. Is there a reason for this? It seems to me that the the variance values are never changed, but are used in the fitness calculation in Individuals.cpp:

double Individual::getFitness(double temp, bool env, bool dd, int generation, double redQueenStrength, double redQueen) {
    double out = (DBL_MIN * 100.0); // TODO: Why this?

    if (env) out += m_envStrength * exp(-0.5 * pow((temp - m_Mean) / m_Variance, 2.0)) + 1 - m_envStrength; // environmental niche
    ...
}
florianhartig commented 7 years ago

I think m_Variance is simply a parameter that controls the attraction to the trait mean? If you set it wider, it the trait distribution within a species will become wider.

The whole procedure was introduced to make sure infraspecific variability remains smaller than interspecific trait variability.

TankredO commented 7 years ago

The m_Variance variable controls only the "width" of the environmental niche niche. Maybe this should be renamed?

Right now the fitness can take values above 1. According to the mortality based reproduction, an individual can't die when its fitness is greater than or equal to 1. However, we are staying in this while loop till we reach numerOfRuns deaths. This seems to me like a huge ressource sink. The m_mortalityStrength variable tries to solve this problem. Every m_mortalityStrength'th step the randomly chosen individual dies, no matter its fitness. Still all steps in between can potentially "get stuck" for a long time.

florianhartig commented 7 years ago

about m_Variance - no, that's not the niche, it's the within species trait attraction during evolution, something very different. I'll open another threat for the mortality problem.

florianhartig commented 7 years ago

About the variance - we could still rename it, I guess it's not a particularly good name. I also wonder if it is properly document - can this parameter be changed from R? Might be nice.

TankredO commented 7 years ago

The variable is hardcoded in c++ at the moment.

TankredO commented 7 years ago

I am not sure about the trait attraction. The variable (variance) is only used for the calculation of the fitness. The evolution method has its own width variable.

void Individual::evolve() {

    //if (m_X_coordinate == 0 && m_Y_coordinate == 0) printInfo();

    double width = 0.01;

    double upperBound = 1.0;
    double lowerBound = 0.0;

    double weightSpecies = 0.2;

    // Environment

    m_Mean = (1.0 - weightSpecies) * m_Mean + weightSpecies * m_Species->m_Mean +
             m_RandomGenerator.randomDouble(-width, width);
    if (m_Mean > upperBound) m_Mean = upperBound - (m_Mean - upperBound);
    else if (m_Mean < lowerBound) m_Mean = lowerBound + std::abs(m_Mean);

    // Competition

    m_CompetitionMarker = (1.0 - weightSpecies) * m_CompetitionMarker + weightSpecies * m_Species->m_CompetitionMean +
                          m_RandomGenerator.randomDouble(-width, width);
    if (m_CompetitionMarker > upperBound) m_CompetitionMarker = upperBound - (m_CompetitionMarker - upperBound);
    else if (m_CompetitionMarker < lowerBound) m_CompetitionMarker = lowerBound + std::abs(m_CompetitionMarker);

    //Neutral

    m_NeutralMarker = (1.0 - weightSpecies) * m_NeutralMarker + weightSpecies * m_Species->m_NeutralMean +
                      m_RandomGenerator.randomDouble(-width, width);
    if (m_NeutralMarker > upperBound) m_NeutralMarker = upperBound - (m_NeutralMarker - upperBound);
    else if (m_NeutralMarker < lowerBound) m_NeutralMarker = lowerBound + std::abs(m_NeutralMarker);

    reportBirth(); // ATTENTION: reportBirth is called here!
}
florianhartig commented 7 years ago

I would multiply instead of dividing in the Gauss function - then you have a parameter that is 1/m_variance, and it would make sense to call this parameter "intraspecificTraitAttraction".

florianhartig commented 7 years ago

Ah, shit, sorry ... ok, then this is the niche width after all. Then maybe better call it niche width

TankredO commented 7 years ago

efb4cffb24d571f4058aa6429462021915c891ad