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;
} |
