butterflies: moved random walks -> dispersal

src/modules/butterflies/CMakeLists.txt

-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
 

src/modules/butterflies/butterflies_dispersal.c

 //
 // 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

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

src/modules/butterflies/random_walk.c

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

src/modules/butterflies/random_walk.h

-#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