generate leslie matrix

3 files changed, 38 insertions, 0 deletions

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 <stdio.h>
 #include <string.h>
 
+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;