Local Parameters

Table of Contents

• Environmental Niche Models
  • Niche models as mathematical expressions
  • Mathematical composition in C++
  • Mathematical composition with quetzal::expressive
• Preparing a suitability landscape for the simulation
• Adjusting a suitability map using a Digital Elevation Model
  • Why separating Elevation from Suitability ?
• Defining a friction model for trans-oceanic dispersal events

Environmental Niche Models

Niche models as mathematical expressions

In the field of ecology, the niche refers to the compatibility between a species and a specific environmental condition. It explains how an organism or a population adapts to the availability of resources and competition.

For instance, when resources are plentiful and the presence of predators, parasites, and pathogens is low, the organism or population tends to grow. In locations and times of resources scarcity and high predation, the populations tend to shrink.

The composition and significance of environmental variables that constitute the dimensions of a niche can vary among species, and the importance of specific environmental factors for a species may change depending on geographical and ecological contexts.

Note

In the broader context of this library, the niche model lacks a specific definition as it relies on the preferences and requirements of individual users. What we can understand is that it involves some combination of spatial and temporal factors that are closely intertwined with the demographic process. This justifies the need for a highly adaptable approach for representing the concept of niche, as it is bound to vary depending on the specific species, population, or geographic area under study.

Given this understanding, the primary goal of Quetzal is to provide users with quetzal::expressive, a user-friendly means to effectively combine geospatial and temporal elements. This enables users to establish meaningful mathematical connections between a diverse and heterogeneous landscape and the underlying demographic processes, based on their own perspectives and preferences. Generally speaking, users will want to use quetzal::expressive to mathematically compose spatio-temporal expressions with a signature double (location_type, time_type) until they reach a satisfying definition of the local demographic parameter of interest.

Mathematical composition in C++

Mathematical composition in C++ is not trivial as unlike languages, C++ does not understand mathematical expressions and their composition. In this strongly-typed language, you can not simply multiply the integer 42 and the function \( f(x,y) = x^y \) and expect it to work. Indeed, if you try to rescale a lambda function, e.g.:

int a = 42;
auto f = [](auto x, auto y){return std::pow(x,y); };
my_expression = a * f; // error

it will not compile. As the compiler does not know natively what the rules are to multiply a number with whatever type f is, it will (verbosely) insult you saying the operator*() is not defined for the type int and an anonymous type (the lambda function).

Mathematical composition with quetzal::expressive

You need a library to tell the compiler what to do. This is where quetzal::expressive comes handy: it allows you to add, substract, multiply, or compose mathematical expressions by giving them a uniform interface.

In the context of demographic models, my_expression in the example below could represent for example r(x,t) the growth rate at location x at time t. Being able to compose several smaller expressions into bigger expressions gives much freedom to users who can for example define r in one part on the program and compose it at a later time with a different expression e.g. k(x,t) the carrying capacity at location x at time t. By treating these local demographic parameters as independent expressions that can be mixed and matched, quetzal::expressive allows great freedom in the implementation details of a spatio-temporal process.

Input

#include "quetzal/quetzal.hpp"
#include <cassert>
#include <cmath> // std::pow
#include <iostream>
 
int main()
{
    // This is used to transform literals into callables
    using quetzal::expressive::literal_factory;
    // This is used to make lambdas usable by expressive
    using quetzal::expressive::use;
 
    // This object will create little functions that can be called with (double, double) parameters
    literal_factory<double, double> lit;
 
    // a is no longer a literal: it's a callable always returning 42
    auto a = lit(42.);
    // check a is indeeed callable with (double, double) parameters
    assert(a(1.0, 3.14) == 42.0);
 
    // This enriches the lambda interface by giving it access to mathematical operators
    auto f = use([](double x, double y) { return std::pow(x, y); });
    // check a is indeeed callable with (double, double) parameters
    assert(f(1.0, 3.14) == 1.0);
 
    // Both a and f can now be invoked with the same interface (double, double) -> double
 
    // We can now compose them easily ...
    auto my_expression = a * f;
 
    // ... and call them possibly (much) later eg in the library's internals
    std::cout << my_expression(1.0, 3.14) << std::endl;
}

Output

42

---

Preparing a suitability landscape for the simulation

Note

The objective of this section is to load a suitability map with quetzal::geography::raster and prepare its integration into the simulation framework with quetzal::expressive.

In the context of ecological niche modeling (ENM), suitability refers to the degree of appropriateness or favorability of a particular environmental condition for a species to thrive and persist. ENM is a technique used in ecology to predict the distribution of species across landscapes based on their known occurrences and environmental variables.

Suitability can be thought of as a spectrum, ranging from conditions that are highly suitable for a species' survival and reproduction to conditions that are unsuitable or even detrimental. ENM attempts to map this spectrum across geographical space by analyzing the relationships between species occurrences and various environmental factors, such as temperature, precipitation, elevation, soil type, and vegetation cover.

It's important to note that suitability does not solely depend on a single environmental variable. Instead, it's a complex interplay of multiple factors that determine whether a species can survive, reproduce, and establish a viable population in a given area. Suitability GIS maps derived from previous analysis steps (see e.g. the niche models in Quetzal-crumbs python package) can then be integrated into the coalescence simulation framework.

The following example code

• loads a suitability map using quetzal::geography::raster
• prepares its integration into the simulation framework using quetzal::expressive
• it gives to the suitability the semantic of a function of space and time f(x,t)
• it gives it access to mathematical operators with quetzal::expressive::use
• then it builds upon it to define an (arbitrary) expression of e.g. the growth rate (but it could be friction, or carrying capacity...)

Input

#include "quetzal/quetzal.hpp"
#include <iostream>
 
using namespace quetzal;
 
int main()
{
 
    // defining some useful type aliases
    using time_type = int;
    using raster_type = geography::raster<time_type>;
    using location_type = raster_type::location_descriptor;
 
    // Transforms simple numbers into functions of space and time
    expressive::literal_factory<location_type, time_type> lit;
 
    // The raster has 10 bands that we will assign to 2001 ... 2011.
    std::vector<time_type> times(10);
    std::iota(times.begin(), times.end(), 2001);
 
    // Load the suitability map
    auto file = std::filesystem::current_path() / "data/suitability.tif";
    auto raster = raster_type::from_file(file, times);
 
    // We need to make choices here concerning how NA are handled
    auto f1 = raster.to_view(); // lightweight functor returning empty optionals for NA
    auto f2 = [&](location_type x, time_type t) {
        return f1(x, t).value_or(0.0);
    };                             // remap empty optionals (NA) to 0.0 suitability value
    auto f3 = expressive::use(f2); // f3 has now access to math operators
    auto r = lit(2) * f3 + lit(1); // the growth rate r is a linear function of the suitability
 
    // r can be later invoked with a location_type and a time_type argument
    auto x = *raster.locations().begin();
    auto t = *raster.times().begin();
    std::cout << "r( x = " << x << ", t = " << t << " ) = " << r(x, t) << std::endl;
}

Output

r( x = 0, t = 0 ) = 1.67259

---

Adjusting a suitability map using a Digital Elevation Model

Note

The objective of this section is to load both a suitability map and a Digital Elevation Model (DEM) with quetzal::geography::landscape and prepare their integration into the simulation framework with quetzal::expressive by dynamically modulating the suitability value as a function of the local elevation.

Why separating Elevation from Suitability ?

It is a common approach to pack a wealth of information into a sole raster file. For instance, this can involve using a single friction file and designating ocean cells with NA values, and obstacles for dispersal with values of 1. Another way is to try to inject elevational data into the suitability map during the Ecological Niche Modeling step.

While these approaches are feasible and occasionally advantageous, they do not consistently represent the optimal or most adaptable solution.

In a broader sense, it's advisable to consider each Geographic Information System (GIS) input file as a means of representing a distinct attribute of the landscape such as suitability, friction, or elevation, utilizing quetzal::geography::landscape to read and align these input files then utilizing quetzal::expressive to adjust or combine these quantities according to the user's modeling intentions.

As an illustration we will show in this example how to adjust the stochasticity of a suitability-driven process as a function of the elevation model. This case and its variations can be applicable in various contexts:

• Conveniently enforcing a threshold value beyond which the suitability undergoes a significant change (e.g., becoming 0). This approach is valuable for examining sky-island systems or for easily estimating an isoline by dynamically adjusting the cutoff parameter value at runtime, as opposed to altering input files for every simulation run.
• Depicting the likelihood of a propagule reaching a nunatak, that is a small-scale suitable shelter protruding from the snow sheet (or ice cap). Typically nunataks are the only places where plants can survive in these environments.

Input

#include "quetzal/quetzal.hpp"
#include <iostream>
#include <random>
 
using namespace quetzal;
 
int main()
{
 
    // defining some useful type aliases
    using key_type = std::string;
    using time_type = int;
    using landscape_type = geography::landscape<key_type, time_type>;
    using location_type = landscape_type::location_descriptor;
 
    // Load the suitability map
    auto file1 = std::filesystem::current_path() / "data/suitability.tif";
    auto file2 = std::filesystem::current_path() / "data/elevation.tif";
 
    // The raster has 10 bands that we will assign to 2001 ... 2011.
    std::vector<time_type> times(10);
    std::iota(times.begin(), times.end(), 2001);
 
    // Makes sure the rasters are aligned
    auto landscape = landscape_type::from_file({{"suit", file1}, {"DEM", file2}}, times);
 
    // lightweight functor returning empty optionals where NA
    auto s_view = landscape["suit"].to_view();
    auto e_view = landscape["DEM"].to_view();
 
    // We need to make choices here concerning how NA are handled
    auto suit = expressive::use([&](location_type x, time_type t) { return s_view(x, t).value_or(0.0); });
    auto elev = expressive::use([&](location_type x, time_type t) { return e_view(x, t).value_or(0.0); });
 
    std::random_device rd;  // a seed source for the random number engine
    std::mt19937 gen(rd()); // mersenne_twister_engine seeded with rd()
 
    // 1) hostile environment above an isoline: particularly useful if 123.0 is a randomly sampled parameter to be
    // estimated
    auto isoline_suitability = [&](location_type x, time_type t) {
        return (elev(x, t) >= 123.0) ? 0.0 : suit(x, t);
    };
 
    // 2) small-scale suitable ice-free shelters randomly popping above the snow cover at high-altitude
    auto nunatak_suitability = [&](location_type x, time_type t) {
        std::bernoulli_distribution d(0.1); // give "false" 9/10 of the time
        bool is_nunatak = d(gen);
        return (elev(x, t) >= 123.0) ? static_cast<double>(is_nunatak) * suit(x, t) : suit(x, t);
    };
 
    // expressions can be passed to a simulator core and later invoked with a location_type and a time_type argument
    auto x = *landscape.locations().begin();
    auto t = *landscape.times().begin();
 
    std::cout << "Suitability w/ isoline ( x = " << x << ", t = " << t << " ) = " << isoline_suitability(x, t)
              << std::endl;
    std::cout << "Suitability w/ shelter ( x = " << x << ", t = " << t << " ) = " << nunatak_suitability(x, t)
              << std::endl;
}

Output

Suitability w/ isoline ( x = 0, t = 0 ) = 0
Suitability w/ shelter ( x = 0, t = 0 ) = 0

---

Defining a friction model for trans-oceanic dispersal events

Note

The objective of this section is to load both a suitability map and a Digital Elevation Model (DEM) with quetzal::geography::landscape and prepare their integration into the simulation framework with quetzal::expressive and account for trans-oceanic dispersal events by dynamically defining the friction and the carrying capacity of a cell as functions of the suitability and elevation value.

Drafting events, in the context of biological dispersal, refer to a phenomenon where organisms are carried across large distances by clinging to or being transported by other organisms, objects, or air currents. This can be thought of as a form of passive dispersal where an organism takes advantage of external forces to move beyond its usual range.

The term "drafting" is borrowed from concepts in physics. In biology, this concept is applied to scenarios where organisms, often small and lightweight, catch a ride on other organisms (like birds, insects, or larger animals), objects (like debris), or wind currents. Drafting events can play a significant role in the dispersal of organisms, especially those with limited mobility.

These drafting events provide an opportunity for organisms to reach new areas that they might not be able to access through their own locomotion. It's an interesting ecological phenomenon that highlights the intricate ways in which different species interact and influence each other's distribution and survival.

The general idea is to define lambda expressions that embed stochastic aspects of the model, while preserving the callable interface (space, time) -> double.

Input

#include "quetzal/quetzal.hpp"
#include <iostream>
#include <random>
 
using namespace quetzal;
 
int main()
{
 
    // defining some useful type aliases
    using key_type = std::string;
    using time_type = int;
    using landscape_type = geography::landscape<key_type, time_type>;
    using location_type = landscape_type::location_descriptor;
 
    // Load the suitability map
    auto file1 = std::filesystem::current_path() / "data/suitability.tif";
    auto file2 = std::filesystem::current_path() / "data/elevation.tif";
 
    // The raster has 10 bands that we will assign to 2001 ... 2011.
    std::vector<time_type> times(10);
    std::iota(times.begin(), times.end(), 2001);
 
    // Makes sure the rasters are aligned
    auto landscape = landscape_type::from_file({{"suit", file1}, {"DEM", file2}}, times);
 
    // lightweight functor returning empty optionals where NA
    auto s_view = landscape["suit"].to_view();
    auto e_view = landscape["DEM"].to_view();
 
    // We need to make choices here concerning how NA are handled
    auto suit = expressive::use([&](location_type x, time_type t) { return s_view(x, t).value_or(0.0); });
    auto elev = expressive::use([&](location_type x, time_type t) { return e_view(x, t).value_or(0.0); });
 
    std::random_device rd;  // a seed source for the random number engine
    std::mt19937 gen(rd()); // mersenne_twister_engine seeded with rd()
 
    auto rafting_capacity = [&](location_type x, time_type t) {
        std::bernoulli_distribution d(0.1);          // give "false" 9/10 of the time
        if (suit(x, t) == 0.0 and elev(x, t) == 0.0) // ocean cell case
            return static_cast<double>(d(gen)) * 2;  // will (rarely) allow 2 individuals to survive in the ocean cell
        else if (suit(x, 0) == 0.0 and elev(x, t) > 0.0) // unhabitable continental cell case
            return 0.0;                                  // suitability is minimal, so should be the capacity
        else                                             // habitable continental cells
            return 100. * suit(x, 0);                    // simple rescaling by the max number of individuals in a cell
    };
 
    auto rafting_friction = [&](location_type x, time_type t) {
        if (suit(x, t) == 0.0 and elev(x, t) == 0.0)     // ocean cell case
            return 0.0;                                  // the raft should move freely
        else if (suit(x, 0) == 0.0 and elev(x, t) > 0.0) // hostile continental cell case
            return 1.00;                                 // max friction as the cell is not attractive
        else                                             // favorable continental cell case
            return 1.0 - suit(x, 0);                     // higher the suitability, easier the travel
    };
 
    // expressions can be passed to a simulator core and later invoked with a location_type and a time_type argument
    auto x = *landscape.locations().begin();
    auto t = *landscape.times().begin();
 
    std::cout << "Friction f( x = " << x << ", t = " << t << " ) = " << rafting_friction(x, t) << std::endl;
    std::cout << "Carrying capacity K( x = " << x << ", t = " << t << " ) = " << rafting_capacity(x, t) << std::endl;
}

Output

Friction f( x = 0, t = 0 ) = 0.663704
Carrying capacity K( x = 0, t = 0 ) = 33.6296

---