diff --git a/src/modules/butterflies/CMakeLists.txt b/src/modules/butterflies/CMakeLists.txt
index b74426a554c7935df7c46671dbcc7f4ee05a5107..e99e173564e7ec84e560963f8ed394201ff04832 100644
--- a/src/modules/butterflies/CMakeLists.txt
+++ b/src/modules/butterflies/CMakeLists.txt
@@ -1,10 +1,10 @@
-add_library(cats-butterflies SHARED "" butterflies_actions.c butterflies_actions.h butterflies_vital_rates.c butterflies_vital_rates.h module.h butterflies_dispersal.c butterflies_dispersal.h butterflies_populations.c butterflies_populations.h butterflies_inline.h butterflies_generations.c butterflies_generations.h butterflies_stats.c)
+add_library(cats-butterflies SHARED "" butterflies_actions.c butterflies_actions.h butterflies_vital_rates.c butterflies_vital_rates.h module.h butterflies_dispersal.c butterflies_dispersal.h butterflies_populations.c butterflies_populations.h butterflies_inline.h butterflies_generations.c butterflies_generations.h butterflies_stats.c butterflies_overlays.c butterflies_overlays.h)
target_include_directories(cats-butterflies PUBLIC ".")
target_sources(cats-butterflies
PRIVATE
- butterflies_main.c butterflies_main.h random_walk.c random_walk.h
+ butterflies_main.c butterflies_main.h
PUBLIC
diff --git a/src/modules/butterflies/butterflies_dispersal.c b/src/modules/butterflies/butterflies_dispersal.c
index 8cc7b964cb8eb3a833c6654e5f18d90ed2ffac00..dcff01f41c4082ad09c1fe75523b1a04a190ca75 100644
--- a/src/modules/butterflies/butterflies_dispersal.c
+++ b/src/modules/butterflies/butterflies_dispersal.c
@@ -1,12 +1,36 @@
//
// Created by gattringera on 21/11/22.
//
+#include "inline_population.h"
+#include "butterflies_inline.h"
+#include "misc/cats_random.h"
+#include "populations/carrying_capacity.h"
+#include "inline_overlays.h"
+
#include <math.h>
+#include <gsl/gsl_randist.h>
#include "butterflies_main.h"
#include "butterflies_dispersal.h"
#include "configuration/configuration.h"
#include "modules/module_header.h"
+const int N_DIRECTIONS = 9;
+/*
+ * directions
+ * 4 3 2
+ * 5 0 1
+ * 6 7 8
+ */
+
+const cats_dt_coord DIRECTION_OFFSETS[9][2] = {
+ {+0, +0}, {+0, +1}, {-1, +1},
+ {-1, +0}, {-1, -1}, {+0, -1},
+ {+1, -1}, {+1, +0}, {+1, +1}
+};
+
+
+
+
void add_dispersed_eggs(struct cats_configuration *conf, struct cats_grid *grid, cats_dt_coord row, cats_dt_coord col, float eggs)
{
//struct conf_data_butterflies *module_conf = CATS_MODULE_DATA;
@@ -14,4 +38,167 @@ void add_dispersed_eggs(struct cats_configuration *conf, struct cats_grid *grid,
struct grid_data_butterflies *data = grid->grid_modules[module_id].module_data;
assert(eggs >= 0);
data->eggs[row][col] += eggs;
+}
+
+
+
+static void inline single_random_walk(struct cats_thread_info *ts, struct cats_grid *grid, cats_dt_coord source_row, cats_dt_coord source_col, int32_t eggs, cats_dt_rates egg_fraction_step, int32_t rw_num)
+{
+
+ const int module_id = CATS_MODULE_ID;
+ const struct grid_data_butterflies *data = grid->grid_modules[module_id].module_data;
+ const struct conf_data_butterflies *module_conf = CATS_MODULE_DATA;
+
+
+ int32_t eggs_left = eggs;
+ const cats_dt_coord max_steps = module_conf->animal_dispersal_max_radius;
+
+ cats_dt_coord steps = 0;
+ cats_dt_coord row = source_row;
+ cats_dt_coord col = source_col;
+
+ const cats_dt_coord rows = grid->dimension.rows;
+ const cats_dt_coord cols = grid->dimension.cols;
+
+ while (steps < max_steps) {
+ const unsigned long int direction = gsl_rng_uniform_int(ts->rng, N_DIRECTIONS);
+ const cats_dt_coord *offsets = DIRECTION_OFFSETS[direction];
+
+
+ const cats_dt_coord row_offset = offsets[0];
+ const cats_dt_coord col_offset = offsets[1];
+
+ assert(row_offset >= -1 && row_offset <= 1);
+ assert(col_offset >= -1 && col_offset <= 1);
+
+ row += row_offset;
+ col += col_offset;
+
+ if (row >= rows || row < 0 || col >= cols || col < 0) {
+ return; // we escaped the simulation extent and got lost
+ }
+
+
+
+ // is the cell a valid dispersal target location?
+ if (! (data->info_layer[row][col] & BF_CELL_VALID_DISPERSAL_TARGET)) {
+ continue;
+ }
+
+ int32_t eggs_to_deposit = (int32_t) ceilf((float) (eggs_left * egg_fraction_step));
+ assert(eggs_to_deposit >= 0);
+
+ if (eggs_to_deposit > eggs_left) {
+ eggs_to_deposit = eggs_left;
+ }
+
+ eggs_left -= eggs_to_deposit;
+ data->eggs[row][col] += (float) eggs_to_deposit;
+
+
+
+
+
+ if (eggs_left == 0) break;
+ steps++;
+ }
+
+ assert(eggs_left >= 0);
+}
+
+
+void butterflies_cell_dispersal(struct cats_grid *grid, struct cats_thread_info *ts, cats_dt_coord row, cats_dt_coord col, bool check_exclusion)
+{
+ const struct cats_configuration *conf = ts->conf;
+ if (check_exclusion
+ && (cell_excluded_by_overlay(conf, row, col) ||cell_excluded_by_generation(grid, row, col))) return;
+
+
+ // total adults: the number of adults that became adult in this cell, possibly exceeding the carrying capacity (thanks to poisson processes)
+ cats_dt_population total_adults = get_adult_population(grid, row, col);
+ cats_dt_population cc = get_carrying_capacity(grid, row, col);
+
+ // total females: how many of the total adults are female, as drawn from a binomial distribution with p = 0.5
+ // can be safely cast, because gsl_ran_binomial will return a number <= total_adults
+ cats_dt_population total_females = (cats_dt_population) gsl_ran_binomial(ts->rng, 0.5, total_adults);
+ assert(total_females > 0 && total_females <= total_adults);
+
+ // total males: the rest
+ cats_dt_population total_males = total_adults - total_females;
+
+ // we need at least one female and one male
+ if (total_females == 0 || total_males == 0) return;
+
+ // here we calculate the number of eggs per female, so we can return early if the cell is unsuitable
+ struct conf_data_butterflies *module_conf = CATS_MODULE_DATA;
+ cats_dt_rates eggs_per_f = calculate_rate(&module_conf->eggs_per_female, total_adults, &grid->param, grid, row, col, NULL);
+ int32_t eggs_per_female = (int32_t) ceill(eggs_per_f);
+ assert(eggs_per_female >= 0);
+
+ if (eggs_per_female == 0) return;
+
+
+ // supernumerous_adults: the number of adults exceeding the carrying capacity
+ cats_dt_population supernumerous_adults = total_adults - cc;
+
+ // supernumerous_females: the (randomly assigned, drawn from a hypergeometric distribution) number of females
+ // in the number of supernumerous adults
+ cats_dt_population supernumerous_females = 0;
+
+
+ if (supernumerous_adults > 0) {
+ // Too many adults in the cell. These must leave the cell, either to be added to the number of
+ // * dispersing females
+ // * discarded (males)
+ supernumerous_females = (cats_dt_population) gsl_ran_hypergeometric(ts->rng, total_females, total_males, supernumerous_adults);
+ }
+
+ // remaining_females: the number of females within the carrying capacity
+ cats_dt_population remaining_females = total_females - supernumerous_females;
+ // at this point the number of 'remaining' adults in the cell is equal or less than the carrying capacity
+
+
+ // how many females will leave the cell (not counting the number of supernumerous ones)?
+ // wandering_females: the number of females that will leave the cell and do a random walk
+ // we draw from a Poisson distribution with a species-specific parameter, the probability to leave times the
+ // number of females
+ cats_dt_rates probability_to_leave = 1.0 - module_conf->probability_to_stay;
+ cats_dt_population wandering_females = poisson(ts->rng, remaining_females * probability_to_leave);
+
+ // but cap the number of females actually available
+ if (wandering_females > remaining_females) {
+ wandering_females = remaining_females;
+ }
+
+ // we need to add the previously subtracted supernumerous females
+ wandering_females += supernumerous_females;
+
+ // stationary_females: the number of females that will not leave the cell, and leave all their eggs here
+ cats_dt_population stationary_females = total_females - wandering_females;
+
+
+ const int module_id = CATS_MODULE_ID;
+ struct grid_data_butterflies *data = grid->grid_modules[module_id].module_data;
+
+ // how many eggs are deposited in the local cell?
+ // all the eggs of all the stationary females
+ cats_dt_rates local_eggs = (stationary_females * eggs_per_female);
+ // plus the local fraction of the eggs for the wandering females
+ cats_dt_rates wandering_local_eggs = (module_conf->egg_fraction_source * eggs_per_female);
+ local_eggs += (wandering_females * wandering_local_eggs);
+ data->eggs[row][col] += (float) ceill(local_eggs);
+
+ int32_t wandering_eggs = (int32_t) ceill(eggs_per_female - wandering_local_eggs);
+
+ if (wandering_eggs == 0) {
+ log_message(LOG_ERROR, "%s: random walk with no eggs to distribute in cell %d %d", __func__, row, col);
+ return;
+ }
+
+ // fixme maybe define threshold over which the random walks are bundled
+ const cats_dt_rates egg_fraction_step = module_conf->egg_fraction_step;
+ for (int32_t rw_number = 0; rw_number < wandering_females; rw_number++) {
+ single_random_walk(ts, grid, row, col, wandering_eggs, egg_fraction_step, rw_number);
+ }
+
}
\ No newline at end of file
diff --git a/src/modules/butterflies/butterflies_dispersal.h b/src/modules/butterflies/butterflies_dispersal.h
index 4017a409fb645905a05f0a7ab6a89f4dd153cad0..6385792ccd2b7055dbc0ad4fa0dac2bf2d07898c 100644
--- a/src/modules/butterflies/butterflies_dispersal.h
+++ b/src/modules/butterflies/butterflies_dispersal.h
@@ -5,4 +5,8 @@
#ifndef CATS_BUTTERFLIES_DISPERSAL_H
#define CATS_BUTTERFLIES_DISPERSAL_H
+#include "threading/threading.h"
+#include "cats/configuration/configuration.h"
+void butterflies_cell_dispersal(struct cats_grid *grid, struct cats_thread_info *ts, cats_dt_coord row, cats_dt_coord col, bool check_exclusion);
+
#endif //CATS_BUTTERFLIES_DISPERSAL_H
diff --git a/src/modules/butterflies/random_walk.c b/src/modules/butterflies/random_walk.c
deleted file mode 100644
index e4e4c976473a4af8e1fdd7f60e7e1230f6559cd7..0000000000000000000000000000000000000000
--- a/src/modules/butterflies/random_walk.c
+++ /dev/null
@@ -1,156 +0,0 @@
-//
-// Created by gattringera on 22/06/22.
-//
-#include <math.h>
-#include "random_walk.h"
-#include "inline_overlays.h"
-#include "modules/module_header.h"
-#include "populations/carrying_capacity.h"
-#include "misc/cats_random.h"
-#include "populations/population.h"
-#include "butterflies_main.h"
-#include "butterflies_inline.h"
-#include "inline_population.h"
-
-const int N_DIRECTIONS = 9;
-/*
- * directions
- * 4 3 2
- * 5 0 1
- * 6 7 8
- */
-
-cats_dt_coord DIRECTION_OFFSETS[9][2] = {
- {+0, +0}, {+0, +1}, {-1, +1},
- {-1, +0}, {-1, -1}, {+0, -1},
- {+1, -1}, {+1, +0}, {+1, +1}
-};
-
-static inline bool suitable_for_egg_laying(const struct cats_configuration *conf, struct cats_grid *grid, cats_dt_coord row, cats_dt_coord col)
-{
-
- if (cell_excluded_by_overlay(conf, row, col)) return false;
-
- if (conf->overlays.overlay[OL_RESOURCE].enabled) {
- double resources = conf->overlays.resources->data[row][col];
- if (resources <= 0) return false;
-
- }
-
- if (conf->overlays.overlay[OL_HABITAT_TYPE_CC].enabled) {
- double multiplier = conf->overlays.habitat_cc->data[row][col];
- if (multiplier < 0.0) { return false; }
- }
-
-
- return true;
-}
-
-static void inline single_random_walk(struct cats_thread_info *ts, struct cats_grid *grid, cats_dt_coord source_row, cats_dt_coord source_col, int32_t eggs, int32_t rw_num)
-{
- const struct cats_configuration *conf = ts->conf;
-
- const int module_id = CATS_MODULE_ID;
- struct grid_data_butterflies *data = grid->grid_modules[module_id].module_data;
- struct conf_data_butterflies *module_conf = CATS_MODULE_DATA;
-
- int32_t eggs_left = eggs;
- cats_dt_coord max_steps = module_conf->animal_dispersal_max_radius;
- //double max_distance = 0;
- //double distance = 0;
- cats_dt_coord steps = 0;
- bool done = false;
- cats_dt_coord row = source_row;
- cats_dt_coord col = source_col;
- int32_t deposition_count = 0;
- /*
- if (row == BF_DEBUG_ROW && col == BF_DEBUG_COL ) {
- printf("rw,run,row,col,step,deposit,eggs_left,thread_id\n");
- printf("rw,%d,%d,%d,%d,%d,%d,%d\n", rw_num, row, col, 0, 0, eggs_left,ts->id);
- }
- */
- while (! done) {
- //bool accepted = false;
- unsigned long int direction = gsl_rng_uniform_int(ts->rng, N_DIRECTIONS);
- cats_dt_coord *offsets = DIRECTION_OFFSETS[direction];
- cats_dt_coord row_offset = offsets[0];
- cats_dt_coord col_offset = offsets[1];
- assert(row_offset >= -1 && row_offset <= 1);
- assert(col_offset >= -1 && col_offset <= 1);
- row += row_offset;
- col += col_offset;
- if (row >= grid->dimension.rows || row < 0 || col >= grid->dimension.cols || col < 0) {
- //printf("out of bounds %d %d - %d %d\n", row, col, grid->dimension.rows, grid->dimension.cols);
- return; // we escaped the simulation extent and got lost
- }
- steps++;
-/*
- double this_distance = sqrt((row_offset * row_offset) + (col_offset * col_offset));
- if (distance + this_distance <= max_distance) {
- accepted = true;
- distance += this_distance;
- } else {
- accepted = false;
- done = true;
- }
- if (!accepted) return;
-*/
-
- int32_t eggs_to_deposit = 0;
- if (suitable_for_egg_laying(conf, grid, row, col)) {
- deposition_count++;
- eggs_to_deposit = (int32_t) ceilf((float) eggs_left/(2.0f));
- assert(eggs_to_deposit >= 0);
-
- if (eggs_to_deposit > eggs_left) {
- eggs_to_deposit = eggs_left;
- }
- eggs_left -= eggs_to_deposit;
- data->eggs[row][col] += (float) eggs_to_deposit;
-
-
- }
- /*
- if (source_row == BF_DEBUG_ROW && source_col == BF_DEBUG_COL ) {
- printf("rw,%d,%d,%d,%d,%d,%d,%d\n", rw_num, row, col, steps, eggs_to_deposit,eggs_left,ts->id);
- }
- */
-
- if (steps >= max_steps || eggs_left == 0) done = true;
- }
- //printf("rw,%d,%d,%d,%d,%d,%d,%d\n", rw_num, row, col, steps, 0,eggs_left,ts->id);
- fflush(stdout);
- assert(eggs_left >= 0);
-}
-
-
-void butterflies_random_walk(struct cats_grid *grid, struct cats_thread_info *ts, cats_dt_coord row, cats_dt_coord col, bool check_exclusion)
-{
- const struct cats_configuration *conf = ts->conf;
- if (check_exclusion
- && (cell_excluded_by_overlay(conf, row, col) ||cell_excluded_by_generation(grid, row, col))) return;
-
- // struct conf_data_butterflies *module_conf = CATS_MODULE_DATA;
- // const int module_id = CATS_MODULE_ID;
- // struct grid_data_butterflies *data = grid->grid_modules[module_id].module_data;
-
- int32_t population = get_adult_population(grid, row, col);
- int32_t dispersing_population = population;
- struct conf_data_butterflies *module_conf = CATS_MODULE_DATA;
- cats_dt_population N = get_adult_population(grid, row, col);
- cats_dt_rates eggs = calculate_rate(&module_conf->eggs_per_female, N, &grid->param, grid, row, col, NULL);
- int32_t e = (int32_t) ceill(eggs);
- assert(e >= 0);
- if (e == 0) return;
- if (dispersing_population > 0) ts->temp1++;
- if (dispersing_population == 0) return;
-
- for (int32_t p = 0; p < dispersing_population; p++) {
-
- single_random_walk(ts, grid, row, col, e, p);
- ts->temp++;
-
-
- }
- //exit_cats(EXIT_SUCCESS);
-}
diff --git a/src/modules/butterflies/random_walk.h b/src/modules/butterflies/random_walk.h
deleted file mode 100644
index f2cdafb486b3c69c118b73e6680b7e2065656905..0000000000000000000000000000000000000000
--- a/src/modules/butterflies/random_walk.h
+++ /dev/null
@@ -1,7 +0,0 @@
-#ifndef CATS_RANDOM_WALK_H
-#define CATS_RANDOM_WALK_H
-#include "cats/configuration/configuration.h"
-#include "threading/threading.h"
-
-void butterflies_random_walk(struct cats_grid *grid, struct cats_thread_info *ts, cats_dt_coord row, cats_dt_coord col, bool check_exclusion);
-#endif //CATS_RANDOM_WALK_H