butterflies: updates to scale factor, dispersal, demographics

src/modules/butterflies/butterflies_dispersal.c

@@ -45,7 +45,7 @@ static void inline single_random_walk(struct cats_thread_info *ts, struct cats_g
         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;
 
@@ -53,7 +53,7 @@ static void inline single_random_walk(struct cats_thread_info *ts, struct cats_g
         const cats_dt_coord cols = grid->dimension.cols;
 
 
-        while (steps < max_steps) {
+        for (cats_dt_coord step = 0; step < max_steps; step++){
                 const unsigned long int direction = gsl_rng_uniform_int(ts->rng, N_DIRECTIONS);
 
 
@@ -74,8 +74,8 @@ static void inline single_random_walk(struct cats_thread_info *ts, struct cats_g
                 // is the cell a valid dispersal target location?
                 if (!(data->info_layer[row][col] & BF_CELL_VALID_DISPERSAL_TARGET)) {
                         if (debug_rw) {
-                                fprintf(module_conf->debug_rw_file, "%d,%d,%d,%d,%d,%d,%d\n", rw_num, row, col, steps + 1,
-                                        module_conf->animal_dispersal_max_radius - steps - 1, 0, eggs_left);
+                                fprintf(module_conf->debug_rw_file, "%d,%d,%d,%d,%d,%d,%d\n", rw_num, row, col, step + 1,
+                                        module_conf->animal_dispersal_max_radius - step - 1, 0, eggs_left);
                         }
 
                         continue;
@@ -93,13 +93,13 @@ static void inline single_random_walk(struct cats_thread_info *ts, struct cats_g
                 data->eggs[row][col] += (float) eggs_to_deposit;
                 ts->temp2++;
                 if (debug_rw) {
-                        fprintf(module_conf->debug_rw_file, "%d,%d,%d,%d,%d,%d,%d\n", rw_num, row, col, steps + 1,
-                                module_conf->animal_dispersal_max_radius - steps - 1, eggs_to_deposit, eggs_left);
+                        fprintf(module_conf->debug_rw_file, "%d,%d,%d,%d,%d,%d,%d\n", rw_num, row, col, step + 1,
+                                module_conf->animal_dispersal_max_radius - step - 1, eggs_to_deposit, eggs_left);
                 }
 
 
                 if (eggs_left == 0) break;
-                steps++;
+
         }
 
         assert(eggs_left >= 0);
@@ -206,6 +206,12 @@ butterflies_cell_dispersal(struct cats_grid *grid, struct cats_thread_info *ts,
         //printf("thread %d: row %d col %d: doing %d rws\n", ts->id, row, col, wandering_females);
 
         if (debug_rw) {
+                fprintf(module_conf->debug_rw_file, "# total source cell eggs = %d\n", (int)ceill(source_cell_eggs));
+                fprintf(module_conf->debug_rw_file, "# stationary females = %d\n", stationary_females);
+                fprintf(module_conf->debug_rw_file, "# wandering females = %d\n", wandering_females);
+                fprintf(module_conf->debug_rw_file, "# total females = %d\n", wandering_females + stationary_females);
+                fprintf(module_conf->debug_rw_file, "# wandering females source cell eggs = %d\n", (int) (wandering_females * (eggs_per_female - eggs_to_disperse_per_female)));
+                fprintf(module_conf->debug_rw_file, "# stationary females source cell eggs = %d\n", (int) (stationary_females * eggs_per_female));
 
                 fprintf(module_conf->debug_rw_file, "id,row,col,step,steps left,eggs deposited,eggs left\n");
 
@@ -215,7 +221,7 @@ butterflies_cell_dispersal(struct cats_grid *grid, struct cats_thread_info *ts,
         for (int32_t rw_number = 0; rw_number < wandering_females; rw_number++) {
                 if (debug_rw) {
 
-                        fprintf(module_conf->debug_rw_file, "%d,%d,%d,%d,0,%d,%d\n", rw_number, row, col,
+                        fprintf(module_conf->debug_rw_file, "%d,%d,%d,0,%d,%d,%d\n", rw_number, row, col,
                                 module_conf->animal_dispersal_max_radius,
                                 (int) ceill(eggs_per_female * module_conf->egg_fraction_source),
                                 (int) ceill(eggs_to_disperse_per_female));

src/modules/butterflies/butterflies_main.c

@@ -10,10 +10,39 @@
 #include "butterflies_action_helpers.h"
 #include "paths/output_paths.h"
 #include "paths/directory_helper.h"
-
+#include "butterflies_populations.h"
 struct cats_global global;
 struct cats_debug_options cats_debug;
 
+#include "lambda/leslie_matrix.h"
+
+
+double *butterflies_matrix(struct cats_configuration *conf, struct lambda_parameters *l_param,
+                                             bool silent, int32_t *N_out){
+        bool print_rate = !silent;
+
+        struct conf_data_butterflies *module_conf = CATS_MODULE_DATA;
+        cats_dt_rates eggs_per_female = calculate_rate_for_matrix(&module_conf->eggs_per_female, l_param, print_rate);
+        cats_dt_rates adults_per_female = calculate_rate_for_matrix(&module_conf->adults_per_female, l_param, print_rate);
+        cats_dt_rates female_fraction = module_conf->female_fraction;
+        cats_dt_rates stationary = module_conf->probability_to_stay;
+        cats_dt_rates mobile = 1.0 - stationary;
+        cats_dt_rates egg_fraction_source = module_conf->egg_fraction_source;
+        cats_dt_rates max_eggs = module_conf->eggs_per_female.max_rate;
+
+        cats_dt_rates eggs_to_adults_rate = bf_egg_to_adult_survival_rate(adults_per_female, max_eggs);
+        cats_dt_rates local_eggs =  (stationary + mobile * egg_fraction_source ) * eggs_per_female;
+        // female -> female;
+        cats_dt_rates result =  local_eggs * eggs_to_adults_rate * female_fraction;
+
+        double *matrix = calloc_or_die(1, sizeof(double));
+
+        *matrix = (double) result;
+        *N_out = 1;
+        return matrix;
+
+}
+
 
 void *butterfly_grid_init(__attribute__((unused)) struct cats_configuration *conf, struct cats_grid *grid,
                           __attribute__((unused)) void *ignored)
@@ -65,9 +94,10 @@ void load_butterflies_species_params(struct cats_configuration *conf, struct cat
 
         load_conf_vital_rate(&data->eggs_per_female, conf, ini, section_name, param);
         load_conf_vital_rate(&data->butterfly_generations, conf, ini, section_name, param);
+        load_conf_vital_rate(&data->adults_per_female, conf, ini, section_name, param);
 
 
-        load_conf_vital_rate(&data->butterfly_egg_to_adult_survival, conf, ini, section_name, param);
+        //load_conf_vital_rate(&data->butterfly_egg_to_adult_survival, conf, ini, section_name, param);
 
         load_conf_value(true, ini, section_name, "butterflies random walk steps maximum",
                         &data->animal_dispersal_max_radius);
@@ -151,6 +181,7 @@ void cats_module_init(struct cats_configuration *conf)
         int32_t id = register_module(conf, module_name, data, flags);
         register_cats_grid_init_function(conf, butterfly_grid_init, butterfly_grid_cleanup);
         register_load_species_param_config_func(conf, load_butterflies_species_params);
+        register_create_leslie_matrix_func(conf, butterflies_matrix);
         bf_add_vital_rates(conf, data);
 
         log_message(LOG_INFO, "Hello from %s (id: %d)\n", module_name, id);

src/modules/butterflies/butterflies_main.h

@@ -47,6 +47,7 @@ struct conf_data_butterflies {
     struct cats_vital_rate eggs_per_female;
     struct cats_vital_rate butterfly_egg_to_adult_survival;
     struct cats_vital_rate butterfly_generations;
+    struct cats_vital_rate adults_per_female;
     cats_dt_rates female_fraction;
 
     bool actions_added;

src/modules/butterflies/butterflies_populations.c

@@ -13,8 +13,17 @@
 #include "butterflies_main.h"
 #include "butterflies_inline.h"
 #include "inline_population.h"
+#include "inline.h"
 
 
+cats_dt_rates bf_egg_to_adult_survival_rate(cats_dt_rates adults_per_female, cats_dt_rates max_eggs)
+{
+
+
+        return adults_per_female/max_eggs;
+
+}
+
 void bf_cell_maturation(struct cats_grid *grid, struct cats_thread_info *ts, cats_dt_coord row, cats_dt_coord col,
                         bool check_exclusion)
 {
@@ -37,7 +46,7 @@ void bf_cell_maturation(struct cats_grid *grid, struct cats_thread_info *ts, cat
                 exit(EXIT_FAILURE);
         }
 
-
+        int orig_eggs = data->eggs[row][col];
         // the number of generations per cell is usually not an integer value
         // the non-integer part of the number of generations is used in the generation which is
         // one greater than the integer part of the number of generations
@@ -62,7 +71,7 @@ void bf_cell_maturation(struct cats_grid *grid, struct cats_thread_info *ts, cat
                 assert(this_generation_fraction <= 1.0);
         }
 
-
+        cats_dt_population max_cc = (cats_dt_population) grid->param.carrying_capacity.max_rate;
         float eggs = this_generation_fraction * data->eggs[row][col];
 
         if (eggs > data->eggs[row][col]) {
@@ -75,14 +84,26 @@ void bf_cell_maturation(struct cats_grid *grid, struct cats_thread_info *ts, cat
 
         assert(data->eggs[row][col] >= 0);
 
-        struct conf_data_butterflies *module_conf = CATS_MODULE_DATA;
-        cats_dt_rates maturation_rate = calculate_rate(&module_conf->butterfly_egg_to_adult_survival, NAN, &grid->param,
-                                                       grid, row, col, NULL);
 
         // not capped, we can have more adults than CC
-        cats_dt_population adults = poisson(ts->rng, eggs * maturation_rate);
+        struct conf_data_butterflies *module_conf = CATS_MODULE_DATA;
+        cats_dt_rates adults_per_female = calculate_rate(&module_conf->adults_per_female, NAN, &grid->param,
+                                                         grid, row, col, NULL);
+
+        cats_dt_rates max_eggs = module_conf->eggs_per_female.max_rate;
+        cats_dt_rates survival = bf_egg_to_adult_survival_rate(adults_per_female, max_eggs);
+        cats_dt_population adults = poisson(ts->rng, eggs * survival);
         assert(adults >= 0);
+        cats_dt_environment suit = get_suitability(grid, row, col);
+        printf("ABCD,%d,%d,%0.3f\n", (int) orig_eggs, adults, suit);
         set_population_ignore_cc(grid, row, col, adults);
+
+        if (adults > max_cc * 10) {
+
+                log_message(LOG_ERROR, "row %d col %d: number of adults %d exceeds %d time smaximum carrying capacity %d (maturation %Lf, eggs %f, max adults per f %Lf), suit %f, OT %Lf, survival %Lf", row, col, adults, 10, max_cc,
+                            survival, eggs, module_conf->adults_per_female.max_rate, suit, grid->param.OT, survival);
+                exit_cats(EXIT_FAILURE);
+        }
 }
 
 

src/modules/butterflies/butterflies_populations.h

@@ -7,4 +7,6 @@
 #include "data/cats_grid.h"
 void bf_cell_maturation(struct cats_grid *grid, struct cats_thread_info *ts, cats_dt_coord row, cats_dt_coord col, bool check_exclusion);
 void bf_area_kill_adults(struct cats_grid *grid, struct cats_thread_info *ts);
+cats_dt_rates bf_egg_to_adult_survival_rate(cats_dt_rates adults_per_female, cats_dt_rates max_eggs);
+
 #endif //CATS_BUTTERFLIES_POPULATIONS_H

src/modules/butterflies/butterflies_vital_rates.c

@@ -16,11 +16,18 @@ void bf_add_vital_rates(struct cats_configuration *conf, struct conf_data_butter
         set_vital_rate_suitability_cutoff_hint(&data->eggs_per_female, HYBRID_SUIT_TS_ZT);
         set_vital_rate_link_hybrid_function(&data->eggs_per_female, conf, LINK_SUITABILITY_SIGMOID);
 
+        // eggs to adults
+        register_module_vital_rate(conf, &data->adults_per_female, "butterflies adults per female");
+        set_vital_rate_name(&data->adults_per_female, "butterflies adults per female");
+        set_vital_rate_suitability_cutoff_hint(&data->adults_per_female, HYBRID_SUIT_TS_ZT);
+        set_vital_rate_link_hybrid_function(&data->adults_per_female, conf, LINK_SUITABILITY_SIGMOID);
 
         // egg to adult survival
+        /*
         register_module_vital_rate(conf, &data->butterfly_egg_to_adult_survival,"butterflies egg to adult survival rate");
         set_vital_rate_name(&data->butterfly_egg_to_adult_survival, "butterflies eggs to adult survival rate");
         set_vital_rate_suitability_cutoff_hint(&data->butterfly_egg_to_adult_survival, HYBRID_SUIT_TS_ZT);
         set_vital_rate_link_hybrid_function(&data->butterfly_egg_to_adult_survival, conf, LINK_SUITABILITY_SIGMOID);
         set_vital_density(&data->butterfly_egg_to_adult_survival, NO_DENSITY_DEP);
+         */
 }