Program Listing for File SimulateDispersal.cpp

Return to documentation for file (necsim/SimulateDispersal.cpp)

// 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.

#include <utility>
#include <thread>
#include <functional>
#include <limits>

#include "SimulateDispersal.h"
#include "Logging.h"
#include "custom_exceptions.h"
#include "file_system.h"

namespace necsim
{
    void SimulateDispersal::setSequential(bool bSequential)
    {
        is_sequential = bSequential;
    }

    void SimulateDispersal::setSimulationParameters(shared_ptr<SimParameters> sim_parameters, bool print)
    {
        if(print)
        {
            writeInfo("********************************\nSetting simulation current_metacommunity_parameters...\n");
        }
        simParameters = std::move(sim_parameters);
        if(simParameters == nullptr)
        {
            throw FatalException("Simulation parameters are nullptr. Please report this bug.");
        }
        if(print)
        {
            simParameters->printSpatialVars();
        }
        simParameters->deme = 1.0;
        simParameters->deme_sample = 1.0;
    }

    void SimulateDispersal::importMaps()
    {
        writeInfo("Starting map import...\n");
        if(simParameters == nullptr)
        {
            throw FatalException("Simulation current_metacommunity_parameters have not been set.");
        }
        density_landscape->setDims(simParameters);
        dispersal_coordinator.setMaps(density_landscape);
        setDispersalParameters();
        density_landscape->calcFineMap();
        density_landscape->calcCoarseMap();
        density_landscape->calcOffset();
        density_landscape->calcHistoricalFineMap();
        density_landscape->calcHistoricalCoarseMap();
        density_landscape->setLandscape(simParameters->landscape_type);
        density_landscape->recalculateHabitatMax();
        data_mask.importSampleMask(simParameters);
    }

    void SimulateDispersal::setSizes()
    {
        checkMaxParameterReference();
        unsigned long i = max_parameter_reference;
        if(num_steps.empty())
        {
            num_steps.insert(1);
        }
        for(const auto &item : num_steps)
        {
            parameter_references.insert(make_pair(item, i));
            i++;
        }
        distances.clear();
        distances.resize(num_repeats * num_steps.size());

    }

    void SimulateDispersal::setDispersalParameters()
    {
        dispersal_coordinator.setRandomNumber(random);
        dispersal_coordinator.setGenerationPtr(&generation);
        dispersal_coordinator.setDispersal(simParameters);

    }

    void SimulateDispersal::setOutputDatabase(string out_database)
    {
        // Check the file is a database
        if(out_database.substr(out_database.length() - 3) != ".db")
        {
            throw FatalException("Output database is not a .db file, check file name.");
        }
        // Open our SQL connection to the database
        database.open(out_database);
        checkMaxParameterReference();
    }

    void SimulateDispersal::setNumberRepeats(unsigned long n)
    {
        num_repeats = n;
    }

    void SimulateDispersal::setNumberSteps(const vector<unsigned long> &s)
    {
        for(const auto item : s)
        {
            num_steps.insert(item);
        }
    }

    void SimulateDispersal::setNumberWorkers(unsigned long n)
    {
        num_workers = max(n, 1UL);
    }

    unsigned long SimulateDispersal::getMaxNumberSteps()
    {
        unsigned long max_number_steps = 0;
        if(!num_steps.empty())
        {
            max_number_steps = *num_steps.rbegin();
        }
        else
        {
            throw FatalException("No steps have been set.");
        }
        return max_number_steps;
    }

    void SimulateDispersal::storeCellList()
    {
        unsigned long total = 0;
        unsigned long cell_total = 0;
        // First count the number of density cells and pick a cell size
        for(unsigned long i = 0; i < simParameters->sample_y_size; i++)
        {
            for(unsigned long j = 0; j < simParameters->sample_x_size; j++)
            {
                if(data_mask.getVal(j, i, 0, 0))
                {
                    cell_total++;
                    total += density_landscape->getVal(j, i, 0, 0, 0.0);
                }
            }
        }
        writeInfo("Choosing from " + to_string(cell_total) + " cells of " + to_string(total) + " individuals.\n");
        cells.resize(total);
        unsigned long ref = 0;
        for(unsigned long i = 0; i < simParameters->sample_y_size; i++)
        {
            for(unsigned long j = 0; j < simParameters->sample_x_size; j++)
            {
                for(unsigned long k = 0; k < density_landscape->getVal(j, i, 0, 0, 0.0); k++)
                {
                    cells[ref].x = j;
                    cells[ref].y = i;
                    ref++;
                }
            }
        }
    }

    const Cell &SimulateDispersal::getRandomCell()
    {
        if(cells.size() == 0)
        {
            throw FatalException("No cells in landscape to simulate.");
        }
        auto index = static_cast<unsigned long>(floor(random->d01() * cells.size()));
        return cells[index];
    }

    void SimulateDispersal::getEndPoint(Cell &this_cell)
    {
        getEndPoint(this_cell, dispersal_coordinator);
    }

    void SimulateDispersal::getEndPoint(Cell &this_cell, DispersalCoordinator &dispersal_coordinator)
    {
        Step tmp_step(this_cell);
        dispersal_coordinator.disperse(tmp_step);
        this_cell.x = tmp_step.x + tmp_step.xwrap * simParameters->sample_x_size;
        this_cell.y = tmp_step.y + tmp_step.ywrap * simParameters->sample_y_size;
        //  return (this->*getValFptr)(dist, angle, this_cell, end_cell);
    }

    void SimulateDispersal::runMeanDispersalDistance()
    {
        writeInfo("Simulating dispersal " + to_string(num_repeats) + " times.\n");
        storeCellList();
        Cell this_cell{};
        this_cell = getRandomCell();
        // Set up the parameter reference
        setSizes();
        for(unsigned long i = 0; i < num_repeats; i++)
        {
            Cell start_cell{};
            if(!is_sequential)
            {
                // This takes into account rejection sampling based on density due to
                // setup process for the cell list
                this_cell = getRandomCell();
            }
            start_cell = this_cell;
            // Check the end point
            getEndPoint(this_cell);
            // Now store the output location
            distances[i] = make_tuple(1, start_cell, distanceBetweenCells(this_cell, start_cell));
        }
        writeInfo("Dispersal simulation complete.\n");
    }

    template <bool chooseRandomCells=false>
    void SimulateDispersal::runDistanceLoop(const unsigned long bidx,
                                            const unsigned long eidx,
                                            const unsigned long num_repeats,
                                            std::mutex &mutex,
                                            unsigned long &finished,
                                            DispersalCoordinator &dispersal_coordinator,
                                            double &generation)
    {
        Cell this_cell{}, start_cell{};

        vector<double> distance_accumulator;
        distance_accumulator.resize(num_steps.size());

        const unsigned long max_number_steps = getMaxNumberSteps();

        for(unsigned long i = bidx; i < eidx; i++)
        {
            mutex.lock();
            writeRepeatInfo(finished);
            finished++;
            mutex.unlock();

            std::fill(distance_accumulator.begin(), distance_accumulator.end(), 0.0);

            if (chooseRandomCells)
            {
                start_cell = getRandomCell();
            } else {
                start_cell = cells[i];
            }

            for(unsigned long k = 0; k < num_repeats; k++)
            {
                // iterator for elements in the set.
                auto step_iterator = num_steps.begin();
                unsigned long step_index = 0;
                this_cell = start_cell;
                generation = 0.0;

                // Keep looping until we get a valid end point
                for(unsigned long j = 1; j <= max_number_steps; j++)
                {
                    getEndPoint(this_cell, dispersal_coordinator);
                    generation += 0.5;

                    if(j == *step_iterator)
                    {
                        distance_accumulator[step_index] += distanceBetweenCells(start_cell, this_cell);

                        step_iterator++;
                        step_index++;
                    }
                }
            }

            unsigned long step_index = 0;

            for(auto step_iterator = num_steps.begin(); step_iterator != num_steps.end(); step_iterator++)
            {
                distances[i * num_steps.size() + step_index] = make_tuple(*step_iterator,
                                                                          start_cell,
                                                                          distance_accumulator[step_index]
                                                                          / static_cast<double>(num_repeats));

                step_index++;
            }
        }
    }

    template <bool chooseRandomCells=false>
    void SimulateDispersal::runDistanceWorker(const unsigned long seed,
                                              const unsigned long bidx,
                                              const unsigned long eidx,
                                              const unsigned long num_repeats,
                                              std::mutex &mutex,
                                              unsigned long &finished)
    {
        mutex.lock();

        shared_ptr<RNGController> thread_random = make_shared<RNGController>();
        thread_random->wipeSeed();
        thread_random->setSeed(seed);

        double thread_generation = 0.0;

        DispersalCoordinator thread_dispersal_coordinator;
        thread_dispersal_coordinator.setMaps(density_landscape);
        thread_dispersal_coordinator.setRandomNumber(thread_random);
        thread_dispersal_coordinator.setGenerationPtr(&thread_generation);
        thread_dispersal_coordinator.setDispersal(simParameters);

        mutex.unlock();

        runDistanceLoop<chooseRandomCells>(bidx, eidx, num_repeats, mutex, finished, thread_dispersal_coordinator, thread_generation);
    }

    void SimulateDispersal::runMeanDistanceTravelled()
    {
        stringstream ss;
        ss << "Simulating dispersal in " << num_repeats << " random habitable cells once for (";
        // The maximum number of steps
        setSizes();
        unsigned long max_number_steps = getMaxNumberSteps();
        for(const auto &item : num_steps)
        {
            ss << item;
            if(item != max_number_steps)
            {
                ss << ", ";
            }
        }
        ss << ") generations ";
        if (num_workers > 1)
        {
            ss << "using " << num_workers << " threads.\n";
        } else {
            ss << "sequentially.\n";
        }
        writeInfo(ss.str());
        storeCellList();

        std::mutex mutex;
        unsigned long finished = 0;

        if (num_workers <= 1)
        {
            runDistanceLoop<true>(0, num_repeats, 1, std::ref(mutex), std::ref(finished), std::ref(dispersal_coordinator), std::ref(generation));
        } else {
            vector<std::thread> threads;
            threads.resize(num_workers);

            for(unsigned long i = 0; i < num_workers; i++)
            {
                threads[i] = std::thread(&SimulateDispersal::runDistanceWorker<true>,
                                         this,
                                         random->i0(std::numeric_limits<unsigned long>::max() - 1),
                                         num_repeats * i / num_workers,
                                         num_repeats * (i + 1) / num_workers,
                                         1,
                                         std::ref(mutex),
                                         std::ref(finished));
            }

            for(unsigned long i = 0; i < num_workers; i++)
            {
                threads[i].join();
            }
        }

        writeRepeatInfo(num_repeats);
        writeInfo("\nDispersal simulation complete.\n");
    }

    void SimulateDispersal::runAllDistanceTravelled()
    {
        unsigned long old_num_repeats = this->num_repeats;

        storeCellList();
        setNumberRepeats(cells.size());

        stringstream ss;
        ss << "Simulating dispersal in all " << num_repeats << " habitable cells " << old_num_repeats << " times for (";
        // The maximum number of steps
        setSizes();
        unsigned long max_number_steps = getMaxNumberSteps();
        for(const auto &item : num_steps)
        {
            ss << item;
            if(item != max_number_steps)
            {
                ss << ", ";
            }
            writeInfo("Dispersal simulation complete.\n");
        }
        ss << ") generations ";
        if (num_workers > 1)
        {
            ss << "using " << num_workers << " threads.\n";
        } else {
            ss << "sequentially.\n";
        }
        writeInfo(ss.str());

        std::mutex mutex;
        unsigned long finished = 0;

        if (num_workers <= 1)
        {
            runDistanceLoop(0, num_repeats, old_num_repeats, std::ref(mutex), std::ref(finished), std::ref(dispersal_coordinator), std::ref(generation));
        } else {
            vector<std::thread> threads;
            threads.resize(num_workers);

            for(unsigned long i = 0; i < num_workers; i++)
            {
                threads[i] = std::thread(&SimulateDispersal::runDistanceWorker,
                                         this,
                                         random->i0(std::numeric_limits<unsigned long>::max() - 1),
                                         num_repeats * i / num_workers,
                                         num_repeats * (i + 1) / num_workers,
                                         old_num_repeats,
                                         std::ref(mutex),
                                         std::ref(finished));
            }

            for(unsigned long i = 0; i < num_workers; i++)
            {
                threads[i].join();
            }
        }

        writeRepeatInfo(num_repeats);
        writeInfo("\nDispersal simulation complete.\n");

        setNumberRepeats(old_num_repeats);
    }

    void SimulateDispersal::runSampleDistanceTravelled(const vector<Cell> &samples)
    {
        cells = samples;

        unsigned long old_num_repeats = this->num_repeats;
        setNumberRepeats(cells.size());

        stringstream ss;
        ss << "Simulating dispersal in " << num_repeats << " sampled habitable cells " << old_num_repeats
           << " times for (";
        // The maximum number of steps
        setSizes();
        unsigned long max_number_steps = getMaxNumberSteps();
        for(const auto &item : num_steps)
        {
            ss << item;
            if(item != max_number_steps)
            {
                ss << ", ";
            }
        }
        ss << ") generations ";
        if (num_workers > 1)
        {
            ss << "using " << num_workers << " threads.\n";
        } else {
            ss << "sequentially.\n";
        }
        writeInfo(ss.str());

        std::mutex mutex;
        unsigned long finished = 0;

        if (num_workers <= 1)
        {
            runDistanceLoop(0, num_repeats, old_num_repeats, std::ref(mutex), std::ref(finished), std::ref(dispersal_coordinator), std::ref(generation));
        } else {
            vector<std::thread> threads;
            threads.resize(num_workers);

            for(unsigned long i = 0; i < num_workers; i++)
            {
                threads[i] = std::thread(&SimulateDispersal::runDistanceWorker,
                                         this,
                                         random->i0(std::numeric_limits<unsigned long>::max() - 1),
                                         num_repeats * i / num_workers,
                                         num_repeats * (i + 1) / num_workers,
                                         old_num_repeats,
                                         std::ref(mutex),
                                         std::ref(finished));
            }

            for(unsigned long i = 0; i < num_workers; i++)
            {
                threads[i].join();
            }
        }

        writeRepeatInfo(num_repeats);
        writeInfo("\nDispersal simulation complete.\n");

        setNumberRepeats(old_num_repeats);
    }

    void SimulateDispersal::writeRepeatInfo(unsigned long i)
    {
        stringstream os;
        os << "\rSimulating dispersal " << i << "/" << num_repeats;
        writeInfo(os.str());
    }

    void SimulateDispersal::writeDatabase(string table_name)
    {
        if(database.isOpen())
        {
            if(table_name != "DISTANCES_TRAVELLED" && table_name != "DISPERSAL_DISTANCES")
            {
                string message = "Table name " + table_name;
                message += "  is not one of 'DISTANCES_TRAVELLED' or 'DISPERSAL_DISTANCES'.";
                throw FatalException(message);
            }
            // Write out the current_metacommunity_parameters
            checkMaxParameterReference();
            writeParameters(table_name);
            // Do the sql output
            // First create the table
            string create_table = "CREATE TABLE IF NOT EXISTS " + table_name + " (id INT PRIMARY KEY not null, ";
            create_table += " x INT NOT NULL, y INT NOT NULL, distance DOUBLE not null, parameter_reference INT NOT NULL);";
            database.execute(create_table);
            string insert_table =
                    "INSERT INTO " + table_name + " (id, x, y, distance, parameter_reference) VALUES (?, ?, ?, ?, ?);";
            auto stmt = database.prepare(insert_table);
            database.beginTransaction();
            unsigned long max_id = checkMaxIdNumber(table_name);
            database.useStatement(stmt); // this could be cleaned up if checkMaxIDNumber comes before the insert statement.
            for(unsigned long i = 0; i < distances.size(); i++)
            {
                sqlite3_bind_int(stmt->stmt, 1, static_cast<int>(max_id + i));
                auto iter = parameter_references.find(get<0>(distances[i]));

#ifdef DEBUG
                if(iter == parameter_references.end())
                {
                    throw FatalException("Cannot find parameter reference. Please report this bug.");
                }
#endif // DEBUG
                unsigned long reference = iter->second;
                if(reference > max_parameter_reference)
                {
                    max_parameter_reference = reference;
                }
                sqlite3_bind_int(stmt->stmt, 2, get<1>(distances[i]).x);
                sqlite3_bind_int(stmt->stmt, 3, get<1>(distances[i]).y);
                sqlite3_bind_double(stmt->stmt, 4, get<2>(distances[i]));
                sqlite3_bind_int(stmt->stmt, 5, static_cast<int>(reference));
                int step = stmt->step();
                if(step != SQLITE_DONE)
                {
                    stringstream ss;
                    ss << "Could not insert into database." << endl;
                    ss << database.getErrorMsg(step);
                    throw FatalException(ss.str());
                }
                stmt->clearAndReset();
            }
            database.endTransaction();
            database.finalise();
        }
        else
        {
            throw FatalException("Database connection has not been opened, check programming.");
        }
        clearParameters();
    }

    void SimulateDispersal::writeParameters(string table_name)
    {
        // Now add the current_metacommunity_parameters
        string create_table = "CREATE TABLE IF NOT EXISTS PARAMETERS (ref INT PRIMARY KEY not null,";
        create_table += "simulation_type TEXT not null, ";
        create_table += " sigma DOUBLE not null, tau DOUBLE not null, m_prob DOUBLE not null, cutoff DOUBLE NOT NULL,";
        create_table += "dispersal_method TEXT not null, map_file TEXT not null, seed INT NOT NULL, number_steps ";
        create_table += "INT NOT NULL, number_repeats INT NOT NULL);";
        database.execute(create_table);
        for(const auto &item : parameter_references)
        {
            string insert_table = "INSERT INTO PARAMETERS VALUES(" + to_string(item.second) + ", '" + table_name + "',";
            insert_table += to_string((long double) simParameters->sigma) + ",";
            insert_table +=
                    to_string((long double) simParameters->tau) + ", " + to_string((long double) simParameters->m_prob);
            insert_table +=
                    ", " + to_string((long double) simParameters->cutoff) + ", '" + simParameters->dispersal_method
                    + "','";
            insert_table +=
                    simParameters->fine_map_file + "', " + to_string(seed) + ", " + to_string(item.first) + ", ";
            insert_table += to_string(num_repeats) + ");";
            database.execute(insert_table);
        }
    }

    void SimulateDispersal::clearParameters()
    {
        distances.clear();
        parameter_references.clear();
        num_steps.clear();
    }

    void SimulateDispersal::checkMaxParameterReference()
    {
        if(database.hasTable("PARAMETERS"))
        {
            string to_exec = "SELECT CASE WHEN COUNT(1) > 0 THEN MAX(ref) ELSE 0 END AS [Value] FROM PARAMETERS;";
            auto stmt = database.prepare(to_exec);
            database.step();
            max_parameter_reference = static_cast<unsigned long>(sqlite3_column_int(stmt->stmt, 0) + 1);
            // close the old statement
            database.finalise();
        }
        else
        {
            max_parameter_reference = 1;
        }

    }

    unsigned long SimulateDispersal::checkMaxIdNumber(string table_name)
    {
        string to_exec = "SELECT CASE WHEN COUNT(1) > 0 THEN MAX(id) ELSE 0 END AS [Value] FROM " + table_name + ";";
        auto stmt = database.prepare(to_exec);
        database.step();
        auto max_id = static_cast<unsigned long>(sqlite3_column_int(stmt->stmt, 0) + 1);
        database.finalise();
        return max_id;
    }

}