1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
|
#include "quetzal/quetzal.h"
#include <boost/math/special_functions/binomial.hpp> // binomial coefficient
#include <random>
#include <algorithm>
#include <cassert>
#include <vector>
int main()
{
/******************************
* Geospatial dataset
*****************************/
using time_type = unsigned int;
using landscape_type = quetzal::geography::DiscreteLandscape<std::string,time_type>;
using coord_type = landscape_type::coord_type;
landscape_type env( {{"rain", "europe_precipitation.tif"}}, {time_type(0)} );
/******************************
* Simulator configuration
*****************************/
using strategy_type = quetzal::demography::strategy::mass_based;
using simulator_type = quetzal::simulators::SpatiallyExplicitCoalescenceSimulator<coord_type, time_type, strategy_type>;
coord_type x_0(40.0, -3.0);
x_0 = env.reproject_to_centroid(x_0);
time_type t_0 = 2000;
strategy_type::value_type N_0 = 100;
simulator_type simulator(x_0, t_0, N_0);
/******************************
* Niche and growth functions
*****************************/
using quetzal::expressive::literal_factory;
using quetzal::expressive::use;
literal_factory<coord_type, time_type> lit;
auto rain = env["rain"];
auto s = [rain](coord_type x, time_type){return rain(x, 0);};
unsigned int factor = 5;
auto K = use(s)*lit(factor);
auto r = lit(5);
auto pop_sizes = simulator.pop_size_history();
auto N = use( [pop_sizes](coord_type x, time_type t){ return pop_sizes(x,t);} );
auto g = N * ( lit(1) + r ) / ( lit(1) + ( (r * N)/K ));
auto sim_children = [g](auto& gen, coord_type const&x, time_type t){
std::poisson_distribution<unsigned int> poisson(g(x,t));
return poisson(gen);
};
std::random_device rd;
std::mt19937 gen(rd());
/******************
* Dispersal
******************/
auto const& space = env.geographic_definition_space();
coord_type::km radius = 60.0;
using quetzal::geography::memoized_great_circle_distance;
auto distance = [](coord_type const& a, coord_type const& b) -> double {return memoized_great_circle_distance(a,b);};
auto weight = [radius](coord_type::km d) -> double { if(d < radius){return 1.0;} else {return 0.0;} ;};
auto f = [weight, distance](coord_type const& a, coord_type const& b){
if(a == b){
return 500.0;
}else{
return weight(distance(a,b));
}
};
auto dispersal = strategy_type::make_sparse_distance_based_dispersal(space, f);
/**************************
* Demographic expansion
**************************/
unsigned int t_sampling = 4000;
simulator.expand_demography(t_sampling, sim_children, dispersal, gen);
/**********************
* Sampling schemes
**********************/
std::vector<coord_type> distribution_area = pop_sizes.get().definition_space(t_sampling);
unsigned int n = 200;
auto last_N = [N, t_sampling](coord_type const& x){return N(x, t_sampling);};
// Uniform sampling
auto unif_sampler = quetzal::sampling_scheme::make_unif_constrained_sampler(distribution_area, last_N, n);
auto global_sample = unif_sampler(gen);
// Clustered sampling
// Radius sampling in Spain
coord_type::km radius_1 = 200.0;
auto weight_1 = [radius_1](coord_type::km d) -> double { if(d < radius_1){return 1.0;} else {return 0.0;} ;};
auto f1 = [weight_1, distance, x_0](coord_type const& x){return weight_1(distance(x_0, x));};
auto spain_sampler = quetzal::sampling_scheme::make_constrained_sampler(distribution_area, f1, last_N, n/2);
auto spain_sample = spain_sampler(gen);
// Radius sampling around Sicilia
coord_type sicilia = env.reproject_to_centroid(coord_type(37.3, 15.0));
coord_type::km radius_2 = 200.0;
auto weight_2 = [radius_2](coord_type::km d) -> double { if(d < radius_2){return 1.0;} else {return 0.0;} ;};
auto f2 = [weight_2, distance, sicilia](coord_type const& x){return weight_2(distance(sicilia, x));};
auto sicilia_sampler = quetzal::sampling_scheme::make_constrained_sampler(distribution_area, f2, last_N, n/2);
auto sicilia_sample = sicilia_sampler(gen);
auto merged_sample = spain_sample;
for(auto const& it: sicilia_sample){
merged_sample[it.first] += it.second;
}
/**********************
* Data visualization
**********************/
std::string N_filename = "/home/me/dev/sampling_project/output/N.tif";
env.export_to_geotiff(N, t_0, t_sampling, [&pop_sizes](time_type const& t){return pop_sizes.get().definition_space(t);}, N_filename);
env.export_to_shapefile(global_sample, "output/uniform.shp");
env.export_to_shapefile(spain_sample, "output/spain.shp" );
env.export_to_shapefile(sicilia_sample, "output/sicilia.shp");
/**********************
* Coalescence
**********************/
unsigned int WF_N = 30;
//std::cout << simulate_newick(simulator, global_sample, t_sampling, t_0, WF_N, gen) << std::endl;
//std::cout << simulate_newick(simulator, spain_sample, t_sampling, t_0, WF_N, gen) << std::endl;
//std::cout << simulate_newick(simulator, sicilia_sample, t_sampling, t_0, WF_N, gen) << std::endl;
auto get_name = [&sicilia_sample, &spain_sample](coord_type const&x, time_type){
if ( sicilia_sample.find(x) == sicilia_sample.end() ) {
return "Spain";
} else {
return "Sicilia";
}
};
std::cout << simulate_newick_leaf_name(simulator, merged_sample, t_sampling, t_0, WF_N, get_name, gen) << std::endl;
return 0;
} |
Partager