butterflies: implement and debug dispersal (random walks)

src/modules/butterflies/butterflies_dispersal.c

@@ -12,5 +12,6 @@ void add_dispersed_eggs(struct cats_configuration *conf, struct cats_grid *grid,
         //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;
+        assert(eggs >= 0);
         data->eggs[row][col] += eggs;
 }
\ No newline at end of file

src/modules/butterflies/butterflies_main.c

@@ -51,16 +51,25 @@ void load_butterflies_species_params(struct cats_configuration *conf, struct cat
         load_conf_vital_rate(&data->butterfly_egg_to_adult_survival, conf, ini, section_name, param);
 
         //bool _get_int32_config_value(const struct cats_ini *ini, const char *section, const char *key, bool required, int32_t *value);
-        int32_t rw_max;
+
         cats_dt_rates female_percentage;
-        load_conf_value(true, ini, section_name, "butterflies random walk steps maximum", &rw_max);
+        load_conf_value(true, ini, section_name, "butterflies random walk steps maximum", &data->animal_dispersal_max_radius);
+
+        if (data->animal_dispersal_max_radius <= 0) {
+                log_message(LOG_ERROR, "butterflies random walk steps maximum must be > 0, is %d", data->animal_dispersal_max_radius);
+                exit_cats(EXIT_FAILURE);
+        }
+
         load_conf_value(true, ini, section_name, "butterflies female fraction", &female_percentage);
         data->generations_max = (int32_t) ceill(data->butterfly_generations.max_rate);
+        param->plant_dispersal_max_radius = data->animal_dispersal_max_radius;
         if (! data->actions_added) {
                 add_actions(conf);
                 data->actions_added = true;
         }
 
+
+
 }
 
 

src/modules/butterflies/butterflies_main.h

@@ -13,8 +13,12 @@ struct grid_data_butterflies {
 
 };
 
+#define BF_DEBUG_ROW  4319
+#define BF_DEBUG_COL 10502
+
 struct conf_data_butterflies {
 
+        // fixme -> move to grid data;
         int32_t generations_max;
         int32_t current_generation;
         int32_t animal_dispersal_max_radius; ///< maximal flight/dispersal distance

src/modules/butterflies/butterflies_populations.c

@@ -23,22 +23,42 @@ void butterflies_cell_maturation(struct cats_grid *grid, struct cats_thread_info
         const int module_id = CATS_MODULE_ID;
         struct grid_data_butterflies *data = grid->grid_modules[module_id].module_data;
         int32_t local_max_generation = (int32_t) ceilf(data->generations[row][col]);
-        if (data->generation_current  > local_max_generation ) return;
-        if (data->eggs[row][col] == 0) return;
+        if (data->eggs[row][col] < 0) {
+                log_message(LOG_ERROR, "Number of eggs < 0: row %d col %d: %f", row, col, data->eggs[row][col]);
+                exit(EXIT_FAILURE);
+        }
+
 
         float proportion = 1.0f;
         if (data->generation_current == local_max_generation && (float) local_max_generation > proportion) {
-                proportion = data->generations[row][col] - (float) local_max_generation;
-        }
+                proportion =  (float) local_max_generation-  data->generations[row][col];
 
+        }
+        assert(proportion >= 0);
+        assert(proportion <= 1.0);
+        if (data->eggs[row][col] == 0) return; // FIXME MOVE UP
 
         assert (conf->grid_count == 1);
         float eggs = proportion * data->eggs[row][col];
+
+        if (eggs > data->eggs[row][col]) {
+                log_message(LOG_ERROR, "Removing more eggs than present: %d %d: %f/%f", row, col, data->eggs[row][col], eggs);
+                exit(EXIT_FAILURE);
+        }
+
         data->eggs[row][col] -= eggs;
+
+        assert(data->eggs[row][col] >= 0);
+
+
         cats_dt_rates maturation_rate = calculate_rate(&module_conf->butterfly_egg_to_adult_survival, NAN, &grid->param, grid, row, col, NULL);
         cats_dt_population K = get_carrying_capacity(grid, row, col);
         cats_dt_population adults = poisson_population_capped(ts->rng, eggs * proportion * maturation_rate, K);
-        set_population_ignore_cc(grid, row, col, adults);
+        /*
+        if (adults >0 ) {
+                printf("adults %d %d: %d\n", row, col, adults);
+        }*/
+        set_population(grid, row, col, adults);
 }
 
 

src/modules/butterflies/butterfly_actions.c

@@ -7,9 +7,14 @@
 #include "butterflies_main.h"
 #include "inline_overlays.h"
 #include "butterflies_populations.h"
-#include "inline_population.h"
 #include "random_walk.h"
 #include "butterflies_inline.h"
+#include "grids/grid_wrapper.h"
+#include "grids/gdal_save.h"
+#include "paths/output_paths.h"
+#include "paths/paths.h"
+#include "populations/population.h"
+#include "inline_population.h"
 
 
 enum action_status action_butterfly_stats_reset(struct cats_grid *grid, struct cats_configuration *conf)
@@ -28,7 +33,7 @@ void grid_butterflies_maturation(struct cats_grid *grid, struct cats_thread_info
         const cats_dt_coord end_col = ts->area.end_col;
 
 
-
+        ts->temp = 0;
         for (cats_dt_coord row = start_row; row < end_row; row++) {
                 for (cats_dt_coord col = start_col; col < end_col; col++) {
 
@@ -36,10 +41,11 @@ void grid_butterflies_maturation(struct cats_grid *grid, struct cats_thread_info
                             || cell_excluded_by_generation(grid, row, col)) {
                                 continue;
                         }
-
                         butterflies_cell_maturation(grid, ts, row, col, false);
+
                 }
         }
+
 }
 
 
@@ -52,22 +58,25 @@ void grid_butterflies_dispersal(struct cats_grid *grid, struct cats_thread_info
         const cats_dt_coord start_col = ts->area.start_col;
         const cats_dt_coord end_col = ts->area.end_col;
 
-        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;
-
+        // 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;
 
+        ts->temp = 0;
+        ts->temp1 = 0;
+        ts->temp2 = 0;
         for (cats_dt_coord row = start_row; row < end_row; row++) {
                 for (cats_dt_coord col = start_col; col < end_col; col++) {
 
                         if (cell_excluded_by_overlay(conf, row, col)
-                            || data->generation_current > (int32_t) ceilf(data->generations[row][col])) {
+                            || cell_excluded_by_generation(grid, row, col)) {
                                 continue;
                         }
 
-                        butterflies_random_walk()
+                        butterflies_random_walk(grid, ts, row, col, false);
                 }
         }
+        // if (ts->temp)   printf("thread %d: %ld random walks (%ld cells) avg %f\n", ts->id, ts->temp, ts->temp1, (float) ts->temp / (float) ts->temp1);
 }
 
 
@@ -94,6 +103,46 @@ enum action_status action_butterflies_maturation(struct cats_grid *grid, struct
         return ACTION_RUN;
 }
 
+char *get_butterfly_population_filename(struct cats_configuration *conf, struct cats_grid *grid)
+{
+
+        assert(grid != NULL);
+        assert(conf != NULL);
+        assert(conf->grid_count == 1);
+
+        int module_id = CATS_MODULE_ID;
+        struct grid_data_butterflies *data = grid->grid_modules[module_id].module_data;
+
+
+        struct string_array *path = get_output_directory(conf, "butterfly-adults"); // FIXME MAKE DIRECTORY
+
+        char *extension = get_extension(conf, "adults");
+        struct string_array *name = standard_output_file_name(conf, NULL, NULL, NULL);
+        string_array_add_int(name, data->generation_current, "g%03d");
+
+        char *filename = assemble_filename(path, name, "_", extension);
+
+        free_string_array(&path);
+        free_string_array(&name);
+        free(extension);
+
+        return filename;
+
+
+
+}
+
+
+enum action_status action_butterflies_save_grid(struct cats_grid *grid, struct cats_configuration *conf)
+{
+        int32_t id = grid->id;
+        char *filename = get_butterfly_population_filename(conf, grid);
+        struct grid_wrapper data = gridwrapper(grid->population, grid->dimension);
+        save_grid_to_gdal(&data, GDT_Int32, conf, filename, conf->param[id].species_name);
+        free(filename);
+        return ACTION_RUN;
+}
+
 
 enum action_status action_butterflies_kill_adults(struct cats_grid *grid, struct cats_configuration *conf)
 {
@@ -105,6 +154,8 @@ enum action_status action_butterflies_kill_adults(struct cats_grid *grid, struct
 
 enum action_status action_butterflies_dispersal(struct cats_grid *grid, struct cats_configuration *conf)
 {
+        threaded_action(&grid_butterflies_dispersal, grid, conf, TS_DISPERSAL);
+
         return ACTION_RUN;
 }
 
@@ -151,6 +202,7 @@ void grid_update_generations(struct cats_grid *grid, struct cats_thread_info *ts
                         }
                         cats_dt_rates gen = calculate_rate(rate, 0, conf->param, grid, row, col, NULL);
                         data->generations[row][col] = (float) gen;
+
                         //printf("GENERATIONS::thread %03d:: %d %d %f\n", ts->id ,row, col, (float) gen);
                 }
         }
@@ -159,15 +211,47 @@ void grid_update_generations(struct cats_grid *grid, struct cats_thread_info *ts
 
 enum action_status action_butterfly_update_generations(struct cats_grid *grid, struct cats_configuration *conf)
 {
-        struct conf_data_butterflies *data = CATS_MODULE_DATA;
+        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;
+
+        data->generation_current = module_conf->generations_max;
+        log_message(LOG_IMPORTANT, "resetting generation to %d", module_conf->generations_max);
+        if (grid->param.initial_population.adjusted == false) {
+                if (grid->param.initial_population.set_to_cc == true) increase_initial_population_to_cc(grid, conf);
+                if (grid->param.initial_population.suitability_threshold > 0.0) {
+                        prune_initial_population_under_threshold(conf, grid);
+                }
+                grid->param.initial_population.adjusted = true;
+        }
+
         threaded_action(&grid_update_generations, grid, conf, TS_DEFAULT);
         return ACTION_RUN;
 }
 
 
-enum action_status action_update_generation(struct cats_grid *grid, struct cats_configuration *conf)
+enum action_status action_finish_generation(struct cats_grid *grid, struct cats_configuration *conf)
 {
-        struct conf_data_butterflies *data = CATS_MODULE_DATA;
+        int module_id = CATS_MODULE_ID;
+
+        struct grid_data_butterflies *data = grid->grid_modules[module_id].module_data;
+
+        data->generation_current--;
+        log_message(LOG_IMPORTANT, "Setting generation to %d\n", data->generation_current);
+        assert(data->generation_current >= 0);
+        // struct conf_data_butterflies *data = CATS_MODULE_DATA;
+        return ACTION_RUN;
+}
+
+enum action_status action_start_generation(struct cats_grid *grid, struct cats_configuration *conf)
+{
+        int module_id = CATS_MODULE_ID;
+
+        struct grid_data_butterflies *data = grid->grid_modules[module_id].module_data;
+
+        log_message(LOG_IMPORTANT, "Starting generation to %d\n", data->generation_current);
+        assert(data->generation_current >= 0);
+        // struct conf_data_butterflies *data = CATS_MODULE_DATA;
         return ACTION_RUN;
 }
 
@@ -186,17 +270,27 @@ void add_butterfly_generation_action(struct cats_configuration *conf, action_fun
 void add_actions(struct cats_configuration *conf)
 {
         struct conf_data_butterflies *data = CATS_MODULE_DATA;
+
         printf("MAXIMUM GENERATIONS %d\n", data->generations_max);
 
         register_action_function(conf, action_butterfly_stats_reset, "butterfly_action_reset_stats",
                                  "resetting butterfly statistics");
         register_action_function(conf, action_butterfly_update_generations, "butterfly_action_update_generations",
                                  "update generations");
+
+        register_action_function(conf, action_finish_generation, "action_finish_generation",
+                                 "update generation");
+        register_action_function(conf, action_start_generation, "action_start_generation",
+                                 "start generation");
+
+
         register_action_function(conf, action_butterflies_maturation, "butterfly_action_egg_to_adult",
                                  "transition eggs to adults");
 
         register_action_function(conf, action_butterflies_dispersal, "butterfly_action_dispersal", "egg dispersal");
         register_action_function(conf, action_butterflies_kill_adults, "butterfly_kill_adults", "kill adults");
+        register_action_function(conf, action_butterflies_save_grid, "action_butterflies_save_grid", "output");
+
         list_actions_full(conf);
 
 
@@ -204,12 +298,16 @@ void add_actions(struct cats_configuration *conf)
         append_action_by_name(conf, "action_load_environments", ALL_STAGES, "environment update", module_name);
         append_action_by_name(conf, "action_overlay_update", ALL_STAGES, "overlay update", module_name);
         append_action(conf, action_butterfly_update_generations, ALL_STAGES, "update generations", module_name);
+
         for (int32_t generation = data->generations_max; generation > 0; generation--) {
-                add_butterfly_generation_action(conf, action_update_generation, "update generation", generation);
+
+                add_butterfly_generation_action(conf, action_start_generation, "start generation", generation);
                 add_butterfly_generation_action(conf, action_butterflies_maturation, "transition eggs to adults",
                                                 generation);
                 add_butterfly_generation_action(conf, action_butterflies_dispersal, "dispersal", generation);
+                add_butterfly_generation_action(conf, action_butterflies_save_grid, "output", generation);
                 add_butterfly_generation_action(conf, action_butterflies_kill_adults, "kill adults", generation);
+                add_butterfly_generation_action(conf, action_finish_generation, "finish generation", generation);
         }
 
 }
\ No newline at end of file

src/modules/butterflies/random_walk.c

@@ -39,44 +39,52 @@ static inline bool suitable_for_egg_laying(const struct cats_configuration *conf
 
         if (conf->overlays.overlay[OL_HABITAT_TYPE_CC].enabled) {
                 double multiplier = conf->overlays.habitat_cc->data[row][col];
-                if (multiplier < 1.0) { return false; }
+                if (multiplier < 0.0) { return false; }
         }
 
 
         return true;
 }
 
-void 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)
+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;
 
-        double max_distance = 0;
-        cats_dt_coord max_steps = 10;
-        double distance = 0;
+        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;
+                //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;
@@ -85,17 +93,34 @@ void single_random_walk(struct cats_thread_info *ts, struct cats_grid *grid, cat
                         accepted = false;
                         done = true;
                 }
+                  if (!accepted) return;
+*/
 
+                int32_t eggs_to_deposit = 0;
                 if (suitable_for_egg_laying(conf, grid, row, col)) {
-                        int32_t eggs_to_deposit= (int32_t) ceilf((float) eggs/(float) deposition_count);
-                        data->eggs[row][col] += (float) eggs_to_deposit;
                         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) done = true;
-                steps += 1;
+                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);
 }
 
 
@@ -105,14 +130,27 @@ void butterflies_random_walk(struct cats_grid *grid, struct cats_thread_info *ts
         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;
+        // 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;
-        int32_t eggs = 0;
+        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, eggs);
+
+                single_random_walk(ts, grid, row, col, e, p);
+                ts->temp++;
+
+
         }
+        //exit_cats(EXIT_SUCCESS);
 }