made Leslie-matrix consistent with start-of-year conditions

+ fixed a bug in juvenile survival rate calculation

src/cats/lambda/leslie_matrix.c

@@ -60,7 +60,7 @@ get_rate_for_matrix(enum cats_vital_rate_id rate_type,
         cats_dt_rates rate = calculate_rate(info, N, param, l_param->grid, l_param->row, l_param->col, &override);
 
         if (print) {
-                log_message(LOG_INFO, "%s: %Lf (maximum %Lf) for population %Lf/%Lf",
+                log_message(LOG_INFO, "%s: %Lf (maximum %Lf) for adult population density %Lf/%Lf",
                             info->name, rate, info->max_rate, N, K);
         }
 
@@ -68,10 +68,22 @@ get_rate_for_matrix(enum cats_vital_rate_id rate_type,
 }
 
 
+
+
+
 double *
 create_leslie_matrix(struct cats_configuration *conf, struct lambda_parameters *l_param, bool silent, int32_t *N_out)
 {
 
+
+#define LM_START
+#ifndef LM_START
+        log_message(LOG_IMPORTANT, "using method LM_MID (old)");
+#else
+        //log_message(LOG_IMPORTANT, "using method LM_START (new)");
+
+#endif
+
         assert (N_out != NULL);
         const int32_t species_id = l_param->species_id;
         const int32_t mat_max = get_vital_age_from_param(&conf->param[species_id], VA_AGE_OF_MATURITY_MAX);
@@ -108,33 +120,35 @@ create_leslie_matrix(struct cats_configuration *conf, struct lambda_parameters *
         cats_dt_rates germination_to_adult_survival = get_rate_for_matrix(VR_GERMINATION_TO_ADULT_SURVIVAL,
                                                                           l_param,
                                                                           print_rate);
-        cats_dt_rates juvenile_survival = get_juvenile_tr_from_germination_to_adult_survival(
+        cats_dt_rates juvenile_transition_rate = get_juvenile_tr_from_germination_to_adult_survival(
                 germination_to_adult_survival, &conf->param[species_id]);
         // printf("%Lf %Lf\n", juvenile_survival, germination_to_adult_survival);
         fflush(stdout);
 
         cats_dt_rates hapaxanthy_survival = 1.0;
         if (l_param->param->hapaxanthy > 0.0) {
-                hapaxanthy_survival = (1.0 - l_param->param->hapaxanthy);
+                hapaxanthy_survival = (1.0 - l_param->param->hapaxanthy) * flowering_frequency;
                 assert(hapaxanthy_survival >= 0);
                 assert(hapaxanthy_survival <= 1.0);
-                log_message(LOG_INFO, "hapaxanthy survival: %Lf", hapaxanthy_survival);
+                log_message(LOG_INFO, "hapaxanthy survival: %Lf, times flowering frequency %Lf", (1.0 - l_param->param->hapaxanthy), hapaxanthy_survival);
         }
 
         if (print_rate) {
-                log_message(LOG_INFO, "juvenile transition: %Lf", juvenile_survival);
+                log_message(LOG_INFO, "juvenile transition: %Lf", juvenile_transition_rate);
                 log_message(LOG_INFO, "local dispersal %Lf", local_dispersal);
         }
 
-        double m0 = (double) age_specific_maturation_rate(0, juvenile_survival, mat_min, mat_max);
-        double d = (double) (flowering_frequency * seed_yield * local_dispersal);
+
+
+        double m0 = (double) age_specific_maturation_rate(0, juvenile_transition_rate, mat_min, mat_max);
+        double dl = (double) (flowering_frequency * seed_yield * local_dispersal);
 
         if (print_rate) {
                 for (int32_t i = 0; i < mat_max + 1; i++) {
                         log_message(LOG_RAW, "m%d = %f\n", i,
-                               (double) age_specific_maturation_rate(i, juvenile_survival, mat_min, mat_max));
+                               (double) age_specific_maturation_rate(i, juvenile_transition_rate, mat_min, mat_max));
                 }
-                log_message(LOG_RAW, "d = %f\n", d);
+                log_message(LOG_RAW, "d = %f\n", dl);
         }
 
         // special case:
@@ -154,34 +168,47 @@ create_leslie_matrix(struct cats_configuration *conf, struct lambda_parameters *
                 //return x;
         }
 
-        // FIRST ROW: transition to seeds1
+        // FIRST ROW: transition to seeds0
 
-        // seeds: Si -> S1
+        // seeds: Si -> S0
         for (int32_t s = 0; s < seed_persistence; s++) {
                 assert(s < N);
-                matrix[0 * N + s] = (double) (germination_rate * m0 * d * seed_survival);
+                matrix[0 * N + s] = (double) (seed_survival * germination_rate * m0 * dl);
         }
 
-        // juveniles: Ji -> S1
+        // juveniles: Ji -> S0
         for (int32_t j = 0; j < mat_max; j++) {
+
                 int32_t col = seed_persistence + j;
                 assert(col < N);
-                cats_dt_rates m = age_specific_maturation_rate(j + 1, juvenile_survival, mat_min, mat_max);
-                matrix[0 * N + col] = (double) (m * d * seed_survival);
+                cats_dt_rates m = age_specific_maturation_rate(j + 1, juvenile_transition_rate, mat_min, mat_max);
+#ifndef LM_START
+                matrix[0 * N + col] = (double) (m * dl * seed_survival);
+
+#else
+                cats_dt_rates juvenile_survival_rate = age_modified_juvenile_survival_rate(
+                        juvenile_transition_rate, j, mat_min, mat_max);
+                matrix[0 * N + col] = (double) (juvenile_survival_rate * m * dl);
+#endif
         }
 
-        // adults: A -> S1
+        // adults: A -> S0
         for (int32_t a = seed_persistence + mat_max; a < N; a++) {
                 assert(a == N - 1);
-                matrix[0 * N + a] = (double) (clonal_growth * d * seed_survival);
+#ifndef LM_START
+                matrix[0 * N + a] = (double) (clonal_growth * dl * seed_survival);
+#else
+                matrix[0 * N + a] = (double) (adult_survival * clonal_growth * dl);
+#endif
         }
 
         // SEED ROWS: S_i -> S_{i+1}
         for (int32_t s = 1; s < seed_persistence; s++) {
-                matrix[s * N + (s - 1)] = (double) ((1.0 - germination_rate) * seed_survival);
+                matrix[s * N + (s - 1)] = (double) (seed_survival * (1.0 - germination_rate));
+
         }
 
-        // JUVENILE ROW 1:
+        // JUVENILE ROW 1: Si -> J0
         for (int32_t s = 0; s < seed_persistence; s++) {
                 int32_t row = seed_persistence;
                 if (row * N + s < N * N) {
@@ -192,17 +219,31 @@ create_leslie_matrix(struct cats_configuration *conf, struct lambda_parameters *
                         log_message(LOG_RAW, "row = %d\n", row);
                         exit(EXIT_FAILURE);
                 }
-                matrix[row * N + s] = (double) (germination_rate * (1.0 - m0) * juvenile_survival);
+#ifndef LM_START
+                cats_dt_rates juvenile_survival_rate = age_modified_juvenile_survival_rate(
+                        juvenile_transition_rate, 0, mat_min, mat_max);
+                matrix[row * N + s] = (double) (germination_rate * (1.0 - m0) * juvenile_survival_rate);
+#else
+                matrix[row * N + s] = (double) ( seed_survival * germination_rate * (1.0 - m0));
+#endif
         }
 
-        // JUVENILE ROWS J2 ...: J_i -> J_{i+1}
+        // JUVENILE ROWS J1+ ...: J_i -> J_{i+1}
 
         for (int32_t j = 0; j < mat_max; j++) {
+                cats_dt_rates juvenile_survival_rate = age_modified_juvenile_survival_rate(
+                        juvenile_transition_rate, j, mat_min, mat_max);
                 int32_t row = j + seed_persistence + 1;
                 int32_t col = j + seed_persistence;
                 assert(row < N);
-                cats_dt_rates m = age_specific_maturation_rate(j + 1, juvenile_survival, mat_min, mat_max);
-                matrix[row * N + col] = (double) ((1.0 - m) * juvenile_survival);
+                cats_dt_rates m = age_specific_maturation_rate(j + 1, juvenile_transition_rate, mat_min, mat_max);
+#ifndef LM_START
+
+                matrix[row * N + col] = (double) ((1.0 - m) * juvenile_survival_rate);
+#else
+
+                matrix[row * N + col] = (double) (juvenile_survival_rate * (1.0 - m));
+#endif
 
         }
 
@@ -210,18 +251,36 @@ create_leslie_matrix(struct cats_configuration *conf, struct lambda_parameters *
         int32_t row = seed_persistence + mat_max;
         assert(row == N - 1);
         for (int32_t s = 0; s < seed_persistence; s++) {
+#ifndef LM_START
                 matrix[row * N + s] = (double) (germination_rate * m0 * adult_survival);
+#else
+                matrix[row * N + s] = (double) (seed_survival *germination_rate * m0 * hapaxanthy_survival);
+#endif
 
         }
         for (int32_t j = 0; j < mat_max; j++) {
+
                 int32_t col = seed_persistence + j;
-                cats_dt_rates m = age_specific_maturation_rate(j + 1, juvenile_survival, mat_min, mat_max);
+                cats_dt_rates m = age_specific_maturation_rate(j + 1, juvenile_transition_rate, mat_min, mat_max);
+#ifndef LM_START
+
                 assert(col < N - 1);
                 matrix[row * N + col] = (double) (m * adult_survival);
+#else
+                cats_dt_rates juvenile_survival_rate = age_modified_juvenile_survival_rate(
+                        juvenile_transition_rate, j, mat_min, mat_max);
+                assert(col < N - 1);
+                matrix[row * N + col] = (double) (juvenile_survival_rate * m * hapaxanthy_survival);
+#endif
         }
 
         // LAST ENTRY
+#ifndef LM_START
         matrix[row * N + N - 1] = (double) (clonal_growth * adult_survival * hapaxanthy_survival);
+#else
+        matrix[row * N + N - 1] = (double) (adult_survival * clonal_growth * hapaxanthy_survival);
+#endif
+
 
         return matrix;
 }