cm.cpp

Go to the documentation of this file.

#include <algorithm>
#include <cassert>

#include "solver.h"
//#include "vector_insert.h"
using namespace std;

// state must be copied to cohorts before calling this function
void Solver::calcRates_CM(double t, vector<double>::iterator S, vector<double>::iterator dSdt){

    //auto ss = species_vec[0]; 
    //cout << "svec: " << t << " " << S[0] << " " << S[1] << " " << S[2] << " " << S[3] << "\n";
    //cout << "coho: " << t << " " << ss->getX(0) << " " << ss->getU(0) << " " << ss->getX(1) << " " << ss->getU(1) << "\n";
    
    vector<double>::iterator its = S    + n_statevars_system; // Skip system variables
    vector<double>::iterator itr = dSdt + n_statevars_system;

    for (int s = 0; s<species_vec.size(); ++s){
        Species_Base* spp = species_vec[s];
        cout << "calcRates(species = " << s << ")\n";
        
        double pe = spp->establishmentProbability(t, env);
        double gb = spp->growthRate(-1, spp->xb, t, env);
        if (spp->birth_flux_in < 0){    
            double birthFlux = calcSpeciesBirthFlux(s,t) * pe;
            spp->set_ub(birthFlux/(gb+1e-20));
        }
        else{
            spp->calc_boundary_u(gb, pe); // this will set u in boundary cohort
        }

        
        for (int i=0; i<spp->J; ++i){
            double x = spp->getX(i);
            vector<double> g_gx = spp->growthRateGradient(i, x, t, env, control.cm_grad_dx);  // FIXME: x can go
            double dx =  g_gx[0];
            double du = -(spp->mortalityRate(i, x, t, env) + g_gx[1]); 
            if (!use_log_densities) du *= spp->getU(i);
        
            *itr++ = dx;
            *itr++ = du;    
            
            its+=2;
        }

        
        if (spp->n_extra_statevars > 0){
            auto itr_prev = itr;
            spp->getExtraRates(itr);
            assert(distance(itr_prev, itr) == spp->n_extra_statevars*spp->J); 
            its += spp->n_extra_statevars*spp->J;   
        }
        //cout << "---\n";
    }
}


void Solver::addCohort_CM(){
    cout << ".......... add cohorts ...............\n";
    updateEnv(current_time, state.begin(), rates.begin());
    
    for (int s = 0; s< species_vec.size(); ++s){
        auto spp = species_vec[s];
        
        // calculate u for boundary cohort
        double gb = spp->growthRate(-1, spp->xb, current_time, env);
        double pe = spp->establishmentProbability(current_time, env);
        if (spp->birth_flux_in < 0){    
            double birthFlux = calcSpeciesBirthFlux(s,current_time) * pe;
            spp->set_ub(birthFlux/(gb+1e-20));
        }
        else{
            spp->calc_boundary_u(gb, pe);  // init density of boundary cohort
        }

        spp->initBoundaryCohort(current_time, env); // init extra state variables and birth time of the boundary cohort
        spp->addCohort();   // introduce copy of boundary cohort in system
    }
    
    resizeStateFromSpecies();
    copyCohortsToState();
}


void Solver::removeCohort_CM(){

    for (auto& spp : species_vec){
        if (spp->J > control.max_cohorts) // remove a cohort if number of cohorts in the species exceeds a threshold
            spp->removeDensestCohort();
    }
    
    resizeStateFromSpecies();
    copyCohortsToState();
    cout << "........................................\n";
    
}



//template<class Model, class Environment>
//double Solver<Model,Environment>::calc_u0_CM(){
        

    //return 0;
//}