Program Listing for File SpatialTree.h

↰ Return to documentation for file (necsim/SpatialTree.h)

// This file is part of necsim project which is released under MIT license.
// See file **LICENSE.txt** or visit https://opensource.org/licenses/MIT) for full license details.
//
#ifndef SPATIALTREE_H
#define SPATIALTREE_H

#include <cstdio>
#include <fstream>
#include <vector>
#include <iostream>
#include <string>
#include <cstring>
#include <cmath>
#include <iomanip>
#include <cmath>
#include <ctime>
#include <ctime>
#include <sqlite3.h>

#ifndef WIN_INSTALL

#include <unistd.h>

#endif

#include <algorithm>
#include <stdexcept>
#include <memory>
// extra boost include - this requires the installation of boost on the system
// note that this requires compilation with the -lboost_filesystem and -lboost_system linkers.
#include <boost/filesystem.hpp>

// include fast-csv-parser by Ben Strasser (available from https://github.com/ben-strasser/fast-cpp-csv-parser)
// for fast file reading

#ifdef use_csv
#include "fast-cpp-csv-parser/csv.h"
#endif

//#ifndef sql_ram
//#define sql_ram
//#endif

// other includes for required files
#include "Tree.h"
#include "Matrix.h"
#include "RNGController.h"
#include "SimParameters.h"
#include "DataPoint.h"
#include "TreeNode.h"
#include "SpeciesList.h"
#include "Landscape.h"
#include "Community.h"
#include "setup.h"
#include "DispersalCoordinator.h"
#include "ActivityMap.h"
#include "Logging.h"
#include "GillespieCalculator.h"

using namespace std;

namespace necsim
{

    class SpatialTree : public virtual Tree
    {
    protected:
        // Our dispersal coordinator for getting dispersal distances and managing calls from the landscape
        DispersalCoordinator dispersal_coordinator;
        // Death probability values across the landscape
        shared_ptr<ActivityMap> death_map;
        // Reproduction probability values across the landscape
        shared_ptr<ActivityMap> reproduction_map;
        // A lineage_indices of new variables which will contain the relevant information for maps and grids.
        //  strings containing the file names to be imported.
        string fine_map_input, coarse_map_input;
        string historical_fine_map_input, historical_coarse_map_input;
        // Landscape object containing both the coarse and fine maps for checking whether or not there is habitat at a
        // particular location.
        shared_ptr<Landscape> landscape;
        // An indexing spatial positioning of the lineages
        Matrix<SpeciesList> grid;
        unsigned long desired_specnum;
        // contains the DataMask for where we should start lineages from.
        DataMask samplegrid;

        // The gillespie variables
        double gillespie_threshold;
        // Matrix of all the probabilities at every location in the map.
        Matrix<GillespieProbability> probabilities;
        // Vector used for holding the priority queue as a binary heap
        vector<GillespieHeapNode> heap;
        // Index to heap position, or UNUSED if cell is not used.
        Matrix<unsigned long> cellToHeapPositions;
        // matrix of self-dispersal probabilities;
        Matrix<double> self_dispersal_probabilities;

        // Total number of individuals present in the simulated world
        unsigned long global_individuals;
        // Mean death rate across the simulated world
        double summed_death_rate;


        // Defines a cell that is unused
        static const unsigned long UNUSED = static_cast<unsigned long>(-1);
#ifdef DEBUG
        unsigned long gillespie_speciation_events{0};
        pair<EventType, CellEventType> last_event;
#endif // DEBUG
    public:
        SpatialTree() : Tree(), dispersal_coordinator(), death_map(make_shared<ActivityMap>()),
                        reproduction_map(make_shared<ActivityMap>()), fine_map_input("none"), coarse_map_input("none"),
                        historical_fine_map_input("none"), historical_coarse_map_input("none"),
                        landscape(make_shared<Landscape>()), grid(), desired_specnum(1), samplegrid(),
                        gillespie_threshold(0.0), probabilities(), heap(), cellToHeapPositions(),
                        self_dispersal_probabilities(), global_individuals(0), summed_death_rate(1.0)
        {

        }

        ~SpatialTree() override = default;

        void runFileChecks() override;

        void checkFolders();

        void setParameters() override;


        // Imports the maps using the variables stored in the class. This function must be run after the set_mapvars() in
        // order to function correctly.
        void importMaps();

        void importActivityMaps();

        unsigned long getInitialCount() override;

        void setupDispersalCoordinator();

        void setup() override;

        unsigned long fillObjects(const unsigned long &initial_count) override;

        unsigned long getIndividualsSampled(const long &x,
                                            const long &y,
                                            const long &x_wrap,
                                            const long &y_wrap,
                                            const double &current_gen);

        unsigned long getNumberLineagesAtLocation(const MapLocation &location) const;

        unsigned long getNumberIndividualsAtLocation(const MapLocation &location) const;

        void removeOldPosition(const unsigned long &chosen) override;

        void calcMove();

        long double calcMinMax(const unsigned long &current);

        void calcNewPos();

        void calcWrappedCoalescence(const unsigned long &nwrap);

        void switchPositions(const unsigned long &chosen) override;

        void calcNextStep() override;

        unsigned long estSpecnum();

#ifdef historical_mode

        void historicalStepChecks();
#endif

        void incrementGeneration() override;

        void updateStepCoalescenceVariables() override;

        void recordLineagePosition();

        void addLineages(double generation_in) override;

        string simulationParametersSqlInsertion() override;

        void simPause() override;

        void dumpMap(shared_ptr<ofstream> out);

        void dumpGrid(shared_ptr<ofstream> out);

        void simResume() override;

        void loadGridSave(shared_ptr<ifstream> in1);

        void loadMapSave(shared_ptr<ifstream> in1);

        void verifyActivityMaps();

        void addWrappedLineage(unsigned long numstart, long x, long y);

        unsigned long countCellExpansion(const long &x,
                                         const long &y,
                                         const long &xwrap,
                                         const long &ywrap,
                                         const double &generationin,
                                         vector<TreeNode> &data_added);

        void expandCell(long x,
                        long y,
                        long x_wrap,
                        long y_wrap,
                        double generation_in,
                        unsigned long add,
                        vector<TreeNode> &data_added,
                        vector<DataPoint> &active_added);

        void addGillespie(const double &g_threshold) override;

        bool runSimulationGillespie() override;

        void runGillespieLoop();

        void setupGillespie();

        void setupGillespieLineages();

        void setupGillespieMaps();

        Cell getCellOfMapLocation(const MapLocation &location);

        void findLocations();

        void checkMapEvents();

        void checkSampleEvents();

        void gillespieCellEvent(GillespieProbability &origin);

        void gillespieUpdateMap();

        void gillespieSampleIndividuals();

        void gillespieCoalescenceEvent(GillespieProbability &origin);

        void gillespieDispersalEvent(GillespieProbability &origin);

        void gillespieSpeciationEvent(GillespieProbability &origin);

        void gillespieLocationRemainingCheck(GillespieProbability &location);

        template<typename T> double getLocalDeathRate(const T &location) const;

        template<typename T> double getLocalSelfDispersalRate(const T &location) const;

        void clearGillespieObjects();

        void setStepVariable(const necsim::GillespieProbability &origin,
                             const unsigned long &chosen,
                             const unsigned long &coal_chosen);

        //        void updateHeapOrigin(GillespieProbability &origin)
        //        {
        //            // If cell is still inhabited
        //            {
        //                Cell pos = convertMapLocationToCell(origin.getMapLocation());
        //
        //                // Update probability and times of origin
        //                updateInhabitedCellOnHeap(pos);
        //
        //            }
        //            // Else
        //            {
        //                std::pop_heap(heap.begin(), heap.end());
        //                heap.pop_back();
        //            }
        //        }

        void gillespieUpdateGeneration(const unsigned long &lineage);

        //        void updateHeapDestination(GillespieProbability &destination)
        //        {
        //            // If dispersal to other cell
        //
        //            Cell pos = convertMapLocationToCell(destination.getMapLocation());
        //
        //            // If cell is already inhabited
        //            {
        //                // Update probability and times of destination
        //                updateInhabitedCellOnHeap(pos);
        //            }
        //            // Else
        //            {
        //                heap.emplace_back(GillespieHeapNode(pos,
        //                                                    (destination.calcTimeToNextEvent(summed_death_rate,
        //                                                                                     getNumberIndividualsAtLocation(
        //                                                                                             destination
        //                                                                                                     .getMapLocation()),
        //                                                                                     global_individuals) + generation),
        //                                                    &cellToHeapPositions.get(pos.y, pos.x)));
        //                std::push_heap(heap.begin(), heap.end());
        //            }
        //        }

        void updateCellCoalescenceProbability(GillespieProbability &origin, const unsigned long &n);

        void updateInhabitedCellOnHeap(const Cell &pos);

        void updateAllProbabilities();

        void removeHeapTop();

        template<typename T> Cell convertMapLocationToCell(const T &location) const
        {
            unsigned long x = landscape->convertSampleXToFineX(location.x, location.xwrap);
            unsigned long y = landscape->convertSampleXToFineX(location.y, location.ywrap);

            return Cell(x, y);
        }

        void createEventList();

        void sortEvents();

        template<bool restoreHeap = true> void addNewEvent(const unsigned long &x, const unsigned long &y);

        void addLocation(const MapLocation &location);

        void setupGillespieProbability(GillespieProbability &gp, const MapLocation &location);

        void fullSetupGillespieProbability(GillespieProbability &gp, const MapLocation &location);

        double calculateCoalescenceProbability(const MapLocation &location) const;

        template<typename T> double convertGlobalGenerationsToLocalGenerations(const T &location, const double g) const
        {
            return g * getLocalDeathRate(location) / (summed_death_rate / (double) global_individuals);
        }

        unsigned long selectRandomLineage(const MapLocation &location) const;

        pair<unsigned long, unsigned long> selectTwoRandomLineages(const MapLocation &location) const;

        vector<unsigned long> detectLineages(const MapLocation &location) const;

        void assignNonSpeciationProbability(const unsigned long chosen);

#ifdef DEBUG

        void validateHeap();

        void debugDispersal();

        void validateLocations();

        void validateLineages() override;

        void debugAddingLineage(unsigned long numstart, long x, long y);


        void runChecks(const unsigned long &chosen, const unsigned long &coalchosen) override;

        void validateGillespie() const;

        void validateSpeciationEvents() const;

        unsigned long countSpeciationEvents() const;

        void checkNoSpeciation(const unsigned long &chosen) const;
#endif

    };
};

#endif  // SPATIALTREE_H