Skip to content
Snippets Groups Projects
Commit cb25c082 authored by Andreas Gattringer's avatar Andreas Gattringer
Browse files

backport:

hybrid mode: vital rate fixes and testing updates

# Conflicts:
#	src/cats/command_line/command_line_options.c
#	src/cats/command_line/command_line_options.h
#	src/cats/debug/debug_vital_rates.c
#	src/cats/debug/debug_vital_rates.h
#	src/cats/vital_rates/hybrid_functions.c
parent 8cb73d2c
No related branches found
No related tags found
No related merge requests found
......@@ -51,16 +51,83 @@ static struct option longopts[] = {
{OPTION_DEBUG_CONFIG_UNSPECIFIED_NAME, no_argument, NULL, OPTION_DEBUG_CONFIG_UNSPECIFIED},
{OPTION_DEBUG_LAMBDA_NAME, no_argument, NULL, OPTION_DEBUG_LAMBDA},
{OPTION_DEBUG_VR_NAME, no_argument, NULL, OPTION_DEBUG_VR},
{OPTION_DEBUG_VR_SCALE_NAME, required_argument, NULL, OPTION_DEBUG_VR_SCALE},
{OPTION_DEBUG_VR_SUBDIVISIONS_NAME, required_argument, NULL, OPTION_DEBUG_VR_SUBDIVISIONS},
{OPTION_DEBUG_VR_OT_NAME, required_argument, NULL, OPTION_DEBUG_VR_OT},
{OPTION_POISSON_DAMPENING_NAME, required_argument, NULL, OPTION_POISSON_DAMPENING},
{OPTION_POISSON_DAMPENING_MIN_NAME, required_argument, NULL, OPTION_POISSON_DAMPENING_MIN},
{OPTION_POISSON_MAXIMUM_NAME, required_argument, NULL, OPTION_POISSON_MAXIMUM},
{OPTION_POISSON_MAX_DRAW_DIFF_NAME, required_argument, NULL, OPTION_POISSON_MAX_DRAW_DIFF},
{NULL, 0, NULL, 0}
};
int32_t convert_long_double(char *string, long double *value, const char *name)
{
bool conversion_success = string_to_long_double(string, value);
if (!conversion_success) {
fprintf(stderr,
"Option --%s: invalid argument '%s'. Must be a floating point number\n",
name, string);
return 1;
}
return 0;
}
int32_t convert_integer(char *string, int32_t *value, const char *name)
{
bool conversion_success = string_to_integer(string, value);
if (!conversion_success) {
fprintf(stderr,
"Option --%s: invalid argument '%s'. Must be an integer\n",
name, string);
return 1;
}
return 0;
}
int32_t check_integer_greater_than(int32_t value, int32_t valid_min_excl, const char *name, const char *msg)
{
if (value > valid_min_excl) {
return 0;
}
if (msg != NULL) {
fprintf(stderr, "Option --%s: %s must be > %d, is %d\n", msg, name, valid_min_excl, value);
} else {
fprintf(stderr, "Option --%s: must be > %d, is %d\n", name, valid_min_excl, value);
}
return 1;
}
int32_t expect_integer_greater_0(char *string, const char *name, int32_t *value)
{
int32_t error = convert_integer(string, value, name);
if (error == 0) {
error += check_integer_greater_than(*value, 0, name, NULL);
}
return error;
}
int32_t expect_long_double_0_1_excl(char *string, const char *name, long double *value)
{
int32_t error = convert_long_double(string, value, name);
if (error == 0
&& (*value <= 0 || *value >= 1)) {
fprintf(stderr, "Option --%s: must be in range (0, 1), is %Lf\n", name, *value);
error += 1;
}
return error;
}
/// @brief assures that CATS was called with at least one command line argument (configuration file) and prints usage otherwise
struct program_options check_cats_main_arguments(int argc, char **argv)
{
......@@ -87,7 +154,11 @@ struct program_options check_cats_main_arguments(int argc, char **argv)
.have_poisson_maximum_draw_diff=false,
.no_input_rasters_required = false,
.debug_vr = false,
.need_conf = true};
.need_conf = true,
.debug_vr_subdivisions = 200,
.debug_vr_scale = -1.0,
.debug_vr_ot = -1.0,
.debug_vr_zt = -1.0};
int opt;
......@@ -202,6 +273,27 @@ struct program_options check_cats_main_arguments(int argc, char **argv)
options.debug_vr = true;
options.need_conf = false;
break;
case OPTION_DEBUG_VR_SUBDIVISIONS:
error += expect_integer_greater_0(optarg, OPTION_DEBUG_VR_SUBDIVISIONS_NAME,
&options.debug_vr_subdivisions);
break;
case OPTION_DEBUG_VR_SCALE:
error += expect_long_double_0_1_excl(optarg, OPTION_DEBUG_VR_SCALE_NAME,
&options.debug_vr_scale);
break;
case OPTION_DEBUG_VR_ZT:
error += expect_long_double_0_1_excl(optarg, OPTION_DEBUG_VR_ZT_NAME,
&options.debug_vr_zt);
break;
case OPTION_DEBUG_VR_OT:
error += expect_long_double_0_1_excl(optarg, OPTION_DEBUG_VR_OT_NAME,
&options.debug_vr_ot);
break;
case OPTION_POISSON_MAXIMUM:
conversion_success = string_to_float(optarg, &value);
if (conversion_success) {
......@@ -280,6 +372,7 @@ struct program_options check_cats_main_arguments(int argc, char **argv)
}
}
if (options.lambda_gradient || options.lambda_test || options.calculate_lambda_quit) {
options.no_input_rasters_required = true;
}
......
......@@ -28,6 +28,13 @@
#include <stdlib.h>
#include <stdbool.h>
#include "logging/logging.h"
#include "cats_global.h"
#define OPTION_LOG_FILE 1111
#define OPTION_LOG_FILE_NAME "log-file"
#define OPTION_SUMMARY_FILE 1112
#define OPTION_SUMMARY_FILE_NAME "summary-file"
#define OPTION_DEBUG_MPI 4001
#define OPTION_DEBUG_MPI_NAME "debug-mpi"
......@@ -47,14 +54,29 @@
#define OPTION_DEBUG_VR 5040
#define OPTION_DEBUG_VR_NAME "debug-vital-rates"
#define OPTION_LOG_FILE 1111
#define OPTION_LOG_FILE_NAME "log-file"
#define OPTION_DEBUG_VR_SUBDIVISIONS 5041
#define OPTION_DEBUG_VR_SUBDIVISIONS_NAME "debug-vital-rates-subdivisions"
#define OPTION_SUMMARY_FILE 1112
#define OPTION_SUMMARY_FILE_NAME "summary-file"
#define OPTION_DEBUG_VR_SCALE 5042
#define OPTION_DEBUG_VR_SCALE_NAME "debug-vital-rates-scale"
#define OPTION_JSON 8000
#define OPTION_JSON_NAME "json"
#define OPTION_DEBUG_VR_OT 5043
#define OPTION_DEBUG_VR_OT_NAME "debug-vital-rates-ot"
#define OPTION_DEBUG_VR_ZT 5044
#define OPTION_DEBUG_VR_ZT_NAME "debug-vital-rates-zt"
#define OPTION_POISSON_DAMPENING 6001
#define OPTION_POISSON_DAMPENING_NAME "poisson-dampening"
#define OPTION_POISSON_DAMPENING_MIN 6002
#define OPTION_POISSON_DAMPENING_MIN_NAME "poisson-dampening-min"
#define OPTION_POISSON_MAXIMUM 6003
#define OPTION_POISSON_MAXIMUM_NAME "poisson-maximum"
#define OPTION_POISSON_MAX_DRAW_DIFF 6004
#define OPTION_POISSON_MAX_DRAW_DIFF_NAME "poisson-maximum-draw-difference"
#define OPTION_SCALE_GRADIENT 7777
#define OPTION_SCALE_GRADIENT_NAME "scale-gradient"
......@@ -71,18 +93,8 @@
#define OPTION_SCALE_LAMBDA_YEARS 7773
#define OPTION_SCALE_LAMBDA_YEARS_NAME "scale-test-years"
#define OPTION_POISSON_DAMPENING 6001
#define OPTION_POISSON_DAMPENING_NAME "poisson-dampening"
#define OPTION_POISSON_DAMPENING_MIN 6002
#define OPTION_POISSON_DAMPENING_MIN_NAME "poisson-dampening-min"
#define OPTION_POISSON_MAXIMUM 6003
#define OPTION_POISSON_MAXIMUM_NAME "poisson-maximum"
#define OPTION_POISSON_MAX_DRAW_DIFF 6004
#define OPTION_POISSON_MAX_DRAW_DIFF_NAME "poisson-maximum-draw-difference"
#include "cats_global.h"
#define OPTION_JSON 8000
#define OPTION_JSON_NAME "json"
struct program_options {
......@@ -112,6 +124,10 @@ struct program_options {
bool need_conf;
bool have_poisson_dampening_minimum;
bool have_poisson_maximum_draw_diff;
cats_dt_rates debug_vr_scale;
cats_dt_rates debug_vr_zt;
cats_dt_rates debug_vr_ot;
int32_t debug_vr_subdivisions;
};
......
......@@ -90,7 +90,7 @@ load_configuration_from_file(const char *filename, const struct program_options
conf->link_count = setup_default_links(conf->vital_dependency_registry);
if (command_line_options->debug_vr) {
debug_vital_rates(conf);
debug_vital_rates(conf, command_line_options);
exit_cats(EXIT_SUCCESS);
}
......
......@@ -30,4 +30,3 @@
#include "cats_datatypes.h"
#include "logging.h"
......@@ -27,6 +27,7 @@
#include "memory/cats_memory.h"
#include "inline_carrying_capacity.h"
#include "grids/grid_setup.h"
#include "paths/directory_helper.h"
struct cats_environment *minimal_suitability_environment(void)
......@@ -95,107 +96,182 @@ void set_suitability(struct cats_environment *env, cats_dt_environment suit)
}
void debug_vital_rates(struct cats_configuration *conf)
FILE *init_debug_vr_file(const struct cats_vital_rate *vr, const struct cats_vital_rate *cc,
const struct cats_species_param *param, int32_t dens_idx)
{
struct cats_environment *env = minimal_suitability_environment();
set_suitability(env, 0.5f);
struct cats_grid *grid = minimal_grid(conf, env);
struct cats_species_param *param = &grid->param;
set_param_values(param, 0.5, 0.25, 0.5);
struct cats_vital_rate *cc = &param->carrying_capacity;
log_message(LOG_RAW, "K %Lf\n", param->carrying_capacity.func->func(cc, param, 0.5, 0, 0));
enum cats_density_dependence density_dep[] = {NO_DENSITY_DEP, DENSITY_DEP_NEGATIVE};
int32_t dens_count = sizeof(density_dep) / sizeof(density_dep[0]);
//cats_dt_rates suit_step = 0.01;
struct cats_vital_rate vr = {.max_rate = 1.0, .environment_set=env};
cats_dt_rates result;
cats_dt_rates N;
cats_dt_rates K_max = param->max_adult_cc_fraction * param->carrying_capacity.max_rate;
cats_dt_rates density;
cats_dt_rates K_suit;
cats_dt_environment suit;
cats_dt_environment suit_step = 0.001f;
cats_dt_rates scale_step = 0.1;
cats_dt_rates OT_step = 0.1;
for (int32_t i = LINK_MIN + 1; i < conf->link_count - 1; i++) {
vr.func = &conf->vital_dependency_registry[i];
for (int32_t dens_idx = 0; dens_idx < dens_count; dens_idx++) {
vr.density = density_dep[dens_idx];
param->scale_factor = scale_step;
do {
param->OT = OT_step;
do {
char *filename = NULL;
FILE *f = NULL;
char *filename = NULL;
// assemble filename and open
int rc = asprintf(&filename,
"debug/rate-%s_density-%d_ZT-%Lf_OT-%Lf_scale-%Lf.csv",
vr.func->short_name, dens_idx, param->ZT, param->OT,
param->scale_factor);
"debug/vital-rates-test/rate-%s_density-%d_ZT-%Lf_OT-%Lf_scale-%Lf_cc-%s.csv",
vr->func->short_name, dens_idx, param->ZT, param->OT,
param->scale_factor,
cc->func->short_name);
asprintf_check(rc);
f = fopen(filename, "w");
ENSURE_FILE_OPENED(f, filename)
free(filename);
log_message(LOG_RAW, "%s - density %d - ZT %Lf - OT %Lf - scale %Lf\n",
vr.func->short_name,
// write header
printf("%s - density %d - ZT %Lf - OT %Lf - scale %Lf - cc %s\n",
vr->func->short_name,
dens_idx,
param->ZT,
param->OT,
param->scale_factor);
fprintf(f, "value,scale,suit,N,K,density,K_suit,ZT,OT,density_dependence\n");
N = 0;
do {
//if (vr.density == NO_DENSITY_DEP && N > 0 ) break;
suit = 0;
density = N / K_max;
param->scale_factor,
cc->func->short_name);
do {
set_suitability(env, suit);
K_suit = get_carrying_capacity(grid, 0, 0);
fprintf(f, "value,scale,suit,N,K,density,K_suit,ZT,OT,density_dependence,demographic_slope\n");
return f;
}
static inline cats_dt_rates attempt_rate_calculation(cats_dt_environment suit, cats_dt_rates N, cats_dt_rates *K_suit,
struct cats_grid *grid, struct cats_vital_rate *vr)
{
if (suit < 0.0 || suit > 1.0) {
*K_suit = NAN;
return NAN;
}
set_suitability(grid->suitability, suit);
*K_suit = get_adult_carrying_capacity(grid, 0, 0);
if (N > *K_suit) return NAN;
if (N > K_suit) {
result = NAN;
} else {
result = calculate_rate(&vr, N, param, grid, 0, 0,
NULL);
return calculate_rate(vr, N, &grid->param, grid, 0, 0, NULL);
}
fprintf(f,
"%0.4Lf,%0.4Lf,%0.4f,%0.4Lf,%0.4Lf,%0.4Lf,%0.4Lf,%0.4Lf,%0.4Lf,%d,%0.4Lf\n",
static inline void
debug_vr_write_line(FILE *f, cats_dt_rates result, cats_dt_environment suit, cats_dt_rates N, cats_dt_rates K_suit,
cats_dt_rates K_max, cats_dt_rates density, int32_t dens_idx, const struct cats_grid *grid)
{
fprintf(f, "%0.6Lf,%0.6Lf,%0.6f,%0.6Lf,%0.6Lf,%0.6Lf,%0.6Lf,%0.6Lf,%0.6Lf,%d,%0.6Lf\n",
result,
param->scale_factor,
env->environments[0]->current.values[0][0],
grid->param.scale_factor,
suit,
N,
K_max,
density,
K_suit,
grid->param.ZT,
grid->param.OT,
dens_idx,
grid->param.demographic_slope);
}
void debug_vital_rate(struct cats_vital_rate *vr, struct cats_vital_rate *cc, struct cats_grid *grid,
int32_t steps, int32_t dens_idx)
{
struct cats_species_param *param = &grid->param;
FILE *f = init_debug_vr_file(vr, cc, param, dens_idx);
const cats_dt_rates K_max = param->max_adult_cc_fraction * param->carrying_capacity.max_rate;
const cats_dt_population N_step = (cats_dt_population) ((float) K_max / (float) steps);
const cats_dt_environment suit_step = (cats_dt_environment) 1.0 / (cats_dt_environment) steps;
cats_dt_rates N = 0;
int32_t N_count = 0;
do {
cats_dt_environment suit = 0;
cats_dt_rates density = N / K_max;
int32_t suit_count = 0;
cats_dt_rates result = NAN;
cats_dt_rates K_suit = NAN;
do {
result = attempt_rate_calculation(suit, N, &K_suit, grid, vr);
debug_vr_write_line(f, result, suit, N, K_suit, K_max, density, dens_idx, grid);
suit_count += 1;
suit += suit_step;
} while (suit <= 1.0);
N += 100.0;
} while (N <= K_max);
} while (suit_count <= steps);
N += N_step;
N_count += 1;
} while (N_count <= steps);
fflush(f);
fclose(f);
}
void debug_vital_rates(struct cats_configuration *conf, const struct program_options *command_line_options)
{
struct cats_environment *env = minimal_suitability_environment();
set_suitability(env, 0.5f);
struct cats_grid *grid = minimal_grid(conf, env);
struct cats_species_param *param = &grid->param;
set_param_values(param, 0.5, 0.25, 0.5);
struct cats_vital_rate *cc = &param->carrying_capacity;
enum cats_density_dependence density_dep[] = {NO_DENSITY_DEP, DENSITY_DEP_NEGATIVE};
int32_t dens_count = sizeof(density_dep) / sizeof(density_dep[0]);
struct cats_vital_rate vr = {.max_rate = 1.0, .environment_set=env};
int32_t divisions = 200;
if (command_line_options->debug_vr_subdivisions > 0) {
divisions = command_line_options->debug_vr_subdivisions;
}
cats_dt_rates scale_start = 0.1;
cats_dt_rates scale_step = 0.2;
cats_dt_rates scale_max = 1.0;
if (command_line_options->debug_vr_scale != -1) {
scale_step = scale_start = scale_max = command_line_options->debug_vr_scale;
}
cats_dt_rates OT_step = 0.2;
cats_dt_rates OT_start = 0.1;
cats_dt_rates OT_max = 1.0;
if (command_line_options->debug_vr_ot != -1) {
OT_step = OT_start = OT_max = command_line_options->debug_vr_ot;
}
struct string_array *output_directory = new_string_array();
string_array_add(output_directory, "debug");
string_array_add(output_directory, "vital-rates");
check_and_create_directory_if_needed(output_directory);
free_string_array(&output_directory);
for (int32_t ci = LINK_MIN + 1; ci < conf->link_count - 1; ci++) {
cc->func = &conf->vital_dependency_registry[ci];
for (int32_t i = LINK_MIN + 1; i < conf->link_count - 1; i++) {
vr.func = &conf->vital_dependency_registry[i];
for (int32_t dens_idx = 0; dens_idx < dens_count; dens_idx++) {
vr.density = density_dep[dens_idx];
param->scale_factor = scale_start;
do {
param->OT = OT_start;
do {
debug_vital_rate(&vr, cc, grid, divisions, dens_idx);
param->OT += OT_step;
} while (param->OT < 1.0);
} while (param->OT < OT_max);
param->scale_factor += scale_step;
} while (param->scale_factor <= 1.0);
} while (param->scale_factor <= scale_max);
}
}
}
vr.density = DENSITY_DEP_NEGATIVE;
cleanup_minimum_suitability_environment(&env);
}
\ No newline at end of file
......@@ -25,7 +25,10 @@
#ifndef CATS_DEBUG_VITAL_RATES_H
#define CATS_DEBUG_VITAL_RATES_H
#include "configuration/configuration.h"
void debug_vital_rates(struct cats_configuration *conf);
struct cats_grid *minimal_grid(struct cats_configuration *conf, struct cats_environment *env);
struct cats_environment *minimal_suitability_environment(void);
void debug_vital_rates(struct cats_configuration *conf, const struct program_options *command_line_options);
#endif //CATS_DEBUG_VITAL_RATES_H
......@@ -105,7 +105,7 @@ density_multiplier(enum cats_density_dependence dependence,
assert(N >= 0);
switch (dependence) { // NOLINT(hicpp-multiway-paths-covered)
switch (dependence) {
case NO_DENSITY_DEP:
return 1.0;
case DENSITY_DEP_NEGATIVE:
......
......@@ -39,7 +39,6 @@ static inline cats_dt_rates hybrid_density_multiplier(const struct cats_vital_ra
cats_dt_rates K)
{
return density_multiplier(rate_info->density, N, K, rate_info->density_ts);
}
......@@ -79,7 +78,6 @@ sigmoid(const struct cats_vital_rate *rate_info,
assert(rate_info != NULL);
if (rate_info->environment_set == NULL) {
log_message(LOG_ERROR, "%s: environment is NULL", __func__);
exit_cats(EXIT_FAILURE);
......@@ -93,7 +91,6 @@ sigmoid(const struct cats_vital_rate *rate_info,
const cats_dt_rates max_rate = rate_info->max_rate;
if (max_rate == 0.0) {
return max_rate;
}
......@@ -101,7 +98,6 @@ sigmoid(const struct cats_vital_rate *rate_info,
cats_dt_rates scale = get_rate_scale_factor(rate_info, species_parameter);
const cats_dt_rates rate_at_OT = max_rate * scale;
assert(max_rate >= 0);
......@@ -127,13 +123,6 @@ sigmoid(const struct cats_vital_rate *rate_info,
+ expl(ds * (OT - suitability))
* (max_rate / rate_at_OT - 1.0)
);
/*
if (rate_info->is_carrying_capacity) {
printf("* %Lf scale %Lf, suit %Lf, N %Lf, K %Lf, OT %Lf, ds %Lf\n",rate, scale, suitability, N, K,OT, ds);
}
*/
const cats_dt_rates density = hybrid_density_multiplier(rate_info, species_parameter, suitability, N, K);
......@@ -143,21 +132,17 @@ sigmoid(const struct cats_vital_rate *rate_info,
K = K - K * rate_info->density_ts;
assert(N >= 0);
assert(K >= 0);
cats_dt_rates old_rate = rate;
switch (rate_info->density) {
switch (rate_info->density) {
case DENSITY_DEP_NEGATIVE:
assert(K >= 0);
assert(N >= 0);
if (K == 0) return 0.0;
rate = rate * density + (N / K) * rate_at_OT; // was: max_rate * scale;
if (rate > old_rate) rate = old_rate;
rate = min_rates(rate, rate * density + ((N / K) * rate_at_OT));
break;
case DENSITY_DEP_POSITIVE:
if (K == 0) return 0.0;
rate = rate * density + (1.0 - (N / K)) * rate_at_OT; // was: max_rate * scale;
if (rate < old_rate) rate = old_rate;
rate = max_rates(rate, rate * density + (1.0 - (N / K)) * rate_at_OT);
break;
case NO_DENSITY_DEP:
break;
......@@ -184,8 +169,7 @@ cats_dt_rates constant_value_rate(const struct cats_vital_rate *rate_info,
assert(!isnan(rate_info->max_rate));
cats_dt_rates rate = rate_info->max_rate;
//const cats_dt_rates density = rate_density(rate_info, species_parameter, suitability, N, K);
return rate; // * density;
return rate;
}
......@@ -203,16 +187,20 @@ cats_dt_rates linear_rate(const struct cats_vital_rate *rate_info,
cats_dt_rates max_rate = rate_info->max_rate;
const cats_dt_rates density = hybrid_density_multiplier(rate_info, species_parameter, suitability, N, K);
cats_dt_rates scale = get_rate_scale_factor(rate_info, species_parameter);
cats_dt_rates slope = (1.0 - scale) / (1.0 - species_parameter->OT);
cats_dt_rates intercept = scale - slope * species_parameter->OT; // FIXME
cats_dt_rates y = slope * suitability + intercept;
y = clamp(y, 0.0, 1.0);
return y * max_rate * density;
const cats_dt_rates density = hybrid_density_multiplier(rate_info, species_parameter, suitability, N, K);
cats_dt_rates rate = y * max_rate;
if (density == 1.0) return rate;
cats_dt_rates rate_at_OT = max_rate * scale;
rate = min_rates(rate, rate * density + ((N / K) * rate_at_OT)); // was: max_rate * scale;
return rate;
}
......@@ -138,7 +138,7 @@ bool string_to_integer(char *string, int32_t *value)
long converted = strtol(string, &end_pointer, 10);
if (strlen(end_pointer) != 0) {
log_message(LOG_WARNING, "%s: invalid or unused characters when converting '%s' to integer %d: '%s'",
log_message(LOG_WARNING, "%s: invalid or unused characters when converting '%s' to integer %ld: '%s'",
__func__, string, converted, end_pointer);
return false;
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment