dispersal_kernel.hpp

// Copyright 2021 Arnaud Becheler    <abechele@umich.edu>
 
/*********************************************************************** * This program is free software; you can
 *redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free
 *Software Foundation; either version 2 of the License, or    * (at your option) any later version. *
 *                                                                      *
 ***************************************************************************/
 
#pragma once
 
#include <cassert>
#include <cmath> // tgamma pow
#include <numbers> // pi 
 
#include <mp-units/math.h>
#include <mp-units/systems/si.h>
 
namespace quetzal
{
namespace demography
{
namespace dispersal_kernel
{
 
using namespace mp_units;
using namespace mp_units::si::unit_symbols; // shortened units notation: km, m, s
using std::numbers::pi;
 
template <QuantityOf<isq::length> Distance = mp_units::quantity<mp_units::si::metre>> 
struct gaussian
{
    using pdf_result_type = quantity<inverse(pow<2>(Distance::reference)), typename Distance::rep>;
 
    struct param_type
    {
        Distance a;
    };
 
    static constexpr pdf_result_type pdf(Distance r, param_type const &p)
    {
        assert(p.a > 0. * m && r >= 0. * m);
        return (1. / (pi * p.a * p.a)) * exp(-(r * r) / (p.a * p.a));
    }
 
    static constexpr Distance mean_dispersal_distance(param_type const &p)
    {
        using std::pow;
        assert(p.a > 0. * m);
        return (p.a * pow(pi, 0.5)) / 2.;
    }
};
 
template <QuantityOf<isq::length> Distance = mp_units::quantity<mp_units::si::metre>> struct logistic
{
    using pdf_result_type = quantity<inverse(pow<2>(Distance::reference)), typename Distance::rep>;
 
    struct param_type
    {
        Distance a;
 
        double b;
    };
 
    static constexpr pdf_result_type pdf(Distance r, param_type const &p)
    {
        using std::exp;
        using std::log;
        using std::pow;
        using std::tgamma;
        Distance a = p.a;
        double b = p.b;
        assert(a > 0. * m && b > 2. && r >= 0. * m);
        return (b / (2. * pi * (a * a) * tgamma(2. / b) * tgamma(1. - 2. / b))) *
               (1. / (1. + (pow(r, b) / (pow(a, b)))));
    }
 
    static constexpr Distance mean_dispersal_distance(param_type const &p)
    {
        using std::tgamma;
        Distance a = p.a;
        double b = p.b;
        assert(a > 0. * m && b > 2);
        return (a * tgamma(3. / b) * tgamma(1. - 3. / b)) / (tgamma(2. / b) * tgamma(1. - 2. / b));
    }
};
 
template <QuantityOf<isq::length> Distance = mp_units::quantity<mp_units::si::metre>> 
struct negative_exponential
{
    using pdf_result_type = quantity<inverse(pow<2>(Distance::reference)), typename Distance::rep>;
 
    struct param_type
    {
        Distance a;
    };
 
    static constexpr pdf_result_type pdf(Distance r, param_type const &p)
    {
        using std::exp;
        Distance a = p.a;
        assert(a > 0. * m && r >= 0. * m);
        return 1. / (2. * pi * a * a) * exp(-r / a);
    }
 
    static constexpr Distance mean_dispersal_distance(param_type const &p)
    {
        Distance a = p.a;
        assert(a > 0. * m);
        return 2. * a;
    }
};
 
template <QuantityOf<isq::length> Distance = mp_units::quantity<mp_units::si::metre>> struct exponential_power
{
    using pdf_result_type = quantity<inverse(pow<2>(Distance::reference)), typename Distance::rep>;
 
    struct param_type
    {
        Distance a;
 
        double b;
    };
 
    static constexpr pdf_result_type pdf(Distance r, param_type const &p)
    {
        //using std::pow;
        //using std::tgamma;
        Distance a = p.a;
        double b = p.b;
        assert(a > 0. * m && b > 0. && r >= 0. * m);
        return b / (2. * pi * a * a * tgamma(2. / b)) * exp(- pow(r.numerical_value_in(m), b) / pow(a.numerical_value_in(m),b));
    }
 
    static constexpr Distance mean_dispersal_distance(param_type const &p)
    {
        using std::tgamma;
        Distance a = p.a;
        double b = p.b;
        assert(a > 0. * m && b > 0.);
        return a * tgamma(3. / b) / tgamma(2. / b);
    }
};
 
template <QuantityOf<isq::length> Distance = mp_units::quantity<mp_units::si::metre>> struct two_dt
{
    using pdf_result_type = quantity<inverse(pow<2>(Distance::reference)), typename Distance::rep>;
 
    struct param_type
    {
        Distance a;
 
        double b;
    };
 
    static constexpr pdf_result_type pdf(Distance r, param_type const &p)
    {
        using std::pow;
        Distance a = p.a;
        double b = p.b;
        assert(a > 0. * m && b > 1. && r >= 0. * m);
        return (b - 1.) / (pi * a * a) * pow((1. + (r * r) / (a * a)), -b);
    }
 
    static constexpr Distance mean_dispersal_distance(param_type const &p)
    {
        using std::pow;
        using std::tgamma;
        Distance a = p.a;
        double b = p.b;
        assert(a > 0. * m && b > 1.);
        if (b < 3. / 2.)
        {
            return std::numeric_limits<double>::infinity();
        }
        return (a * pow(pi, 0.5) * tgamma(b - 3. / 2.)) / (2. * tgamma(b - 1.));
    }
};
 
template <QuantityOf<isq::length> Distance = mp_units::quantity<mp_units::si::metre>> struct inverse_power_law
{
    using pdf_result_type = quantity<inverse(pow<2>(Distance::reference)), typename Distance::rep>;
 
    struct param_type
    {
        Distance a;
 
        double b;
    };
 
    static constexpr pdf_result_type pdf(Distance r, param_type const &p)
    {
        using std::pow;
        Distance a = p.a;
        double b = p.b;
        assert(a > 0. * m && b > 2. && r >= 0. * m);
        return pow(1. + r / a, -b) * (b - 2.) * (b - 1.) / (2. * pi * a * a);
    }
 
    static constexpr Distance mean_dispersal_distance(param_type const &p)
    {
        Distance a = p.a;
        double b = p.b;
        assert(a > 0. * m && b > 2.);
        if (b < 3.)
        {
            return std::numeric_limits<double>::infinity();
        }
        return 2. * a / (b - 3.);
    }
};
 
template <QuantityOf<isq::length> Distance = mp_units::quantity<mp_units::si::metre>> struct lognormal
{
    using pdf_result_type = quantity<inverse(pow<2>(Distance::reference)), typename Distance::rep>;
 
    struct param_type
    {
        Distance a;
 
        double b;
    };
 
    // static constexpr pdf_result_type pdf(Distance r, param_type const& p)
    // {
    //   Distance a = p.a;
    //   double b = p.b;
    //   assert(a > 0. * m && b > 0. && r >= 0. * m);
    //   return (1./(pow(2.*pi, 3./.2)*b*r*r)) * exp(- (mp_units::log(r/a)*mp_units::log(r/a) / (2.*b*b) ) );
    // }
 
    static constexpr Distance mean_dispersal_distance(param_type const &p)
    {
        using std::exp;
        Distance a = p.a;
        double b = p.b;
        assert(a > 0. * m && b > 0.);
        return a * exp(b * b / 2.);
    }
};
 
template <QuantityOf<isq::length> Distance = mp_units::quantity<mp_units::si::metre>> struct gaussian_mixture
{
    using pdf_result_type = quantity<inverse(pow<2>(Distance::reference)), typename Distance::rep>;
 
    struct param_type
    {
        Distance a1;
        Distance a2;
        double p;
 
      public:
    };
 
    static constexpr pdf_result_type pdf(Distance r, param_type const &param)
    {
        using std::exp;
        Distance a1 = param.a1;
        Distance a2 = param.a2;
        double p = param.p;
        assert(a1 > 0. * m && a2 > 0. * m && p > 0. && p < 1. && r >= 0. * m);
        return p / (pi * a1 * a1) * exp(-(r * r) / (a1 * a1)) +
               (1. - p) / (pi * a2 * a2) * exp(-(r * r) / (a2 * a2));
    }
 
    static constexpr Distance mean_dispersal_distance(param_type const &param)
    {
        using std::pow;
        Distance a1 = param.a1;
        Distance a2 = param.a2;
        double p = param.p;
        assert(a1 > 0. * m && a2 > 0. * m && p > 0. && p < 1.);
        return pow(pi, 0.5) * (p * a1 + (1. - p) * a2) / 2.;
    }
};
} // namespace dispersal_kernel
} // namespace demography
} // namespace quetzal