diff --git a/inc/lizfcm.h b/inc/lizfcm.h index 85a29bb..4cdba8d 100644 --- a/inc/lizfcm.h +++ b/inc/lizfcm.h @@ -76,4 +76,6 @@ extern double fixed_point_secant_bisection_method(double (*f)(double), extern double dominant_eigenvalue(Matrix_double *m, Array_double *v, double tolerance, size_t max_iterations); +extern Matrix_double *leslie_matrix(Array_double *age_class_surivor_ratio, + Array_double *age_class_offspring); #endif // LIZFCM_H diff --git a/src/eigen.c b/src/eigen.c index aaacedc..6ea26c3 100644 --- a/src/eigen.c +++ b/src/eigen.c @@ -4,6 +4,23 @@ #include #include +Matrix_double *leslie_matrix(Array_double *age_class_surivor_ratio, + Array_double *age_class_offspring) { + assert(age_class_surivor_ratio->size + 1 == age_class_offspring->size); + + Matrix_double *leslie = InitMatrixWithSize(double, age_class_offspring->size, + age_class_offspring->size, 0.0); + + free_vector(leslie->data[0]); + leslie->data[0] = age_class_offspring; + + for (size_t i = 0; i < age_class_surivor_ratio->size; i++) { + leslie->data[i + 1]->data[i] = age_class_surivor_ratio->data[i]; + } + + return leslie; +} + double dominant_eigenvalue(Matrix_double *m, Array_double *v, double tolerance, size_t max_iterations) { assert(m->rows == m->cols); diff --git a/test/eigen.t.c b/test/eigen.t.c index 985a304..9b2f6f0 100644 --- a/test/eigen.t.c +++ b/test/eigen.t.c @@ -1,5 +1,24 @@ #include "lizfcm.test.h" +UTEST(eigen, leslie_matrix) { + Array_double *felicity = InitArray(double, {0.0, 1.5, 0.8}); + Array_double *survivor_ratios = InitArray(double, {0.8, 0.55}); + + Matrix_double *m = InitMatrixWithSize(double, 3, 3, 0.0); + m->data[0]->data[0] = 0.0; + m->data[0]->data[1] = 1.5; + m->data[0]->data[2] = 0.8; + m->data[1]->data[0] = 0.8; + m->data[2]->data[1] = 0.55; + + Matrix_double *leslie = leslie_matrix(survivor_ratios, felicity); + + EXPECT_TRUE(matrix_equal(leslie, m)); + + free_matrix(leslie); + free_matrix(m); +} + UTEST(eigen, dominant_eigenvalue) { Matrix_double *m = InitMatrixWithSize(double, 2, 2, 0.0); m->data[0]->data[0] = 2.0;