workiing with the input data matrix transposed

Author: ZocoLini

modules/cuqdyn-c/src/cuqdyn.c

@@ -12,6 +12,7 @@
 #include "matlab.h"
 #include "ode_solver.h"
 #include "states_transformer.h"
+#include "sunmatrix/sunmatrix_dense.h"
 #ifdef MPI
 #include <mpi.h>
 #endif
@@ -65,7 +66,7 @@ CuqdynResult *cuqdyn_algo(const char *data_file, const char *sacess_conf_file, c
     CuqdynConf *config = get_cuqdyn_conf(get_cuqdyn_context());
 
     N_Vector times = NULL;
-    ObservedData observed_data = NULL;
+    TransposedObservedData observed_data = NULL;
 
     if (read_data_file(data_file, &times, &observed_data) != 0)
     {
@@ -78,16 +79,16 @@ CuqdynResult *cuqdyn_algo(const char *data_file, const char *sacess_conf_file, c
 
     if (config->y0.len == 0)
     {
-        initial_condition = copy_matrix_row(observed_data, 0, 0, SM_COLUMNS_D(observed_data));
+        initial_condition = copy_matrix_column(observed_data, 0, 0, SM_ROWS_D(observed_data));
     }
     else
     {
         initial_condition = New_Serial(config->y0.len);
         memcpy(NV_DATA_S(initial_condition), config->y0.array, config->y0.len * sizeof(sunrealtype));
     }
 
-    const long m = SM_ROWS_D(observed_data);
-    const long n = SM_COLUMNS_D(observed_data);
+    const long n = SM_ROWS_D(observed_data);
+    const long m = SM_COLUMNS_D(observed_data);
 
     SUNMatrix resid_loo = NULL;
     MatrixArray media_matrix = create_matrix_array(m - 1);
@@ -156,7 +157,7 @@ CuqdynResult *cuqdyn_algo(const char *data_file, const char *sacess_conf_file, c
         LongArray indices_to_remove = create_array((long[]) {i + 1}, 1);
 
         N_Vector texp = copy_vector_remove_indices(times, indices_to_remove);
-        SUNMatrix yexp = copy_matrix_remove_rows(observed_data, indices_to_remove);
+        SUNMatrix yexp = copy_matrix_remove_columns(observed_data, indices_to_remove);
 
         N_Vector tmp_initial_condition = copy_vector_remove_indices(initial_condition, create_array((long[]) {}, 0));
 
@@ -169,7 +170,7 @@ CuqdynResult *cuqdyn_algo(const char *data_file, const char *sacess_conf_file, c
 
         for (int j = 0; j < n; ++j)
         {
-            sunrealtype observed = SM_ELEMENT_D(observed_data, i, j);
+            sunrealtype observed = SM_ELEMENT_D(observed_data, j, i);
             sunrealtype predicted = SM_ELEMENT_D(predicted_obs_states, j, i);
 
             NV_Ith_S(residuals, j) = fabs(observed - predicted);

modules/cuqdyn-c/src/data_reader.c

@@ -7,6 +7,8 @@
 #include <string.h>
 
 #include "../include/cuqdyn.h"
+#include "sundials/sundials_types.h"
+#include "sunmatrix/sunmatrix_dense.h"
 
 int read_data_file(const char *data_file, N_Vector *t, ObservedData *y)
 {
@@ -25,7 +27,7 @@ int read_data_file(const char *data_file, N_Vector *t, ObservedData *y)
     return 1;
 }
 
-int read_txt_data_file(const char *data_file, N_Vector *t, ObservedData *y)
+int read_txt_data_file(const char *data_file, N_Vector *t, TransposedObservedData *y)
 {
     const char *ext = strrchr(data_file, '.');
     if (ext && strcmp(ext, ".txt") != 0)
@@ -46,19 +48,19 @@ int read_txt_data_file(const char *data_file, N_Vector *t, ObservedData *y)
     *t = New_Serial(rows);
     sunrealtype *data_t = N_VGetArrayPointer(*t);
 
-    *y = NewDenseMatrix(rows, cols - 1);
-    sunrealtype *data_y = ((SUNMatrixContent_Dense)(*y)->content)->data;
+    *y = NewDenseMatrix(cols - 1, rows);
 
     double tmp;
 
     for (int i = 0; i < rows; ++i)
     {
         fscanf(f, "%lf", &tmp);
         data_t[i] = tmp;
+        sunrealtype *data_yi = ((SUNMatrixContent_Dense)(*y)->content)->cols[i];
         for (int j = 0; j < cols - 1; ++j)
         {
             fscanf(f, "%lf", &tmp);
-            data_y[j * rows + i] = tmp;
+            data_yi[j] = tmp;
         }
     }
 

modules/cuqdyn-c/src/functions.c

@@ -44,15 +44,13 @@ void *obj_func(double *x, void *data)
     const long cols = SM_COLUMNS_D(result);
     sunrealtype J = 0.0;
 
-    // We compare rows with cols because the result matrix is transposed
-    if (SM_ROWS_D(exptotal->yexp) != cols)
+    if (SM_ROWS_D(exptotal->yexp) != rows)
     {
         fprintf(stderr, "ERROR: The yexp rows don't match the ode result rows: %ld vs %ld\n", SM_ROWS_D(exptotal->yexp), rows);
         exit(-1);
     }
 
-    // We compare rows with cols because the result matrix is transposed
-    if (SM_COLUMNS_D(exptotal->yexp) != rows)
+    if (SM_COLUMNS_D(exptotal->yexp) != cols)
     {
         fprintf(stderr, "ERROR: The yexp cols don't match the ode result cols: %ld vs %ld\n", SM_COLUMNS_D(exptotal->yexp), cols - 1);
         exit(-1);
@@ -62,7 +60,7 @@ void *obj_func(double *x, void *data)
     {
         for (long j = 0; j < cols; ++j)
         {
-            const sunrealtype diff = SM_ELEMENT_D(result, i, j) - SM_ELEMENT_D(exptotal->yexp, j, i);
+            const sunrealtype diff = SM_ELEMENT_D(result, i, j) - SM_ELEMENT_D(exptotal->yexp, i, j);
             J += diff * diff;
         }
     }

tests/data/lotka_volterra_data_homoc_noise_0.10_size_30_data_1.mat

@@ -1,4 +1,3 @@
-#define LOTKA_VOLTERRA_MAT "data/lotka_volterra_data_homoc_noise_0.10_size_30_data_1.mat"
 #define DATA_TXT "data/test_data.txt"
 
 #include <nvector/nvector_serial.h>
@@ -7,7 +6,6 @@
 
 
 #include "data_reader.h"
-#include "cuqdyn.h"
 
 void test_read_data_txt();
 
@@ -26,10 +24,13 @@ void test_read_data_txt()
 
     assert(read_txt_data_file(DATA_TXT, &t, &y) == 0);
 
+    assert(y != NULL);
+    assert(t != NULL);
+    
     for (int i = 0; i < SM_ROWS_D(y); ++i)
     {
         assert(NV_Ith_S(t, i) == i);
-        assert(SM_ELEMENT_D(y, i, 0) == i);
-        assert(SM_ELEMENT_D(y, i, 1) == i);
+        assert(SM_ELEMENT_D(y, 0, i) == i);
+        assert(SM_ELEMENT_D(y, 1, i) == i);
     }
 }

tests/test_data_reader.c

@@ -3,6 +3,7 @@
 #include "matlab.h"
 #include "config.h"
 #include "method_module/structure_paralleltestbed.h"
+#include "sunmatrix/sunmatrix_dense.h"
 
 #define CUQDYN_CONF "data/obj_function_cuqdyn_config.xml"
 #define SACESS_CONF "data/obj_function_ess_config_nl2sol.dn2fb.xml"
@@ -16,7 +17,7 @@ int main()
     SUNMatrix yexp = NULL;
 
     read_data_file(DATA, &texp, &yexp);
-    N_Vector initial_values = copy_matrix_row(yexp, 0, 0, SM_COLUMNS_D(yexp));
+    N_Vector initial_values = copy_matrix_row(yexp, 0, 0, SM_ROWS_D(yexp));
 
     experiment_total *exp_total = malloc(sizeof(experiment_total));
     create_expetiment_struct(SACESS_CONF, &exp_total[0], 1, 0, "output", 1, texp, yexp, initial_values);

tests/test_ess_obj_function.c

@@ -15,6 +15,7 @@
 #include <cuqdyn.h>
 #include <stdio.h>
 #include <sunmatrix/sunmatrix_dense.h>
+#include <math.h>
 #include "nvector/nvector_serial.h"
 
 #include "data_reader.h"
@@ -64,7 +65,7 @@ void lotka_volterra_ess(char *conf_file)
 
     read_data_file(LOTKA_VOLTERRA_DATA, &texp, &yexp);
 
-    N_Vector initial_values = copy_matrix_row(yexp, 0, 0, SM_COLUMNS_D(yexp));
+    N_Vector initial_values = copy_matrix_column(yexp, 0, 0, SM_ROWS_D(yexp));
 
     N_Vector xbest = execute_ess_solver(conf_file, OUPUT_PATH, texp, yexp, initial_values, NULL);
 
@@ -73,8 +74,7 @@ void lotka_volterra_ess(char *conf_file)
         sunrealtype expected = expected_values[i];
         sunrealtype result = NV_Ith_S(xbest, i);
 
-        sunrealtype a = 6;
-        // assert(fabs(result - expected) < 0.1);
+        assert(fabs(result - expected) < 1);
     }
 
     destroy_cuqdyn_context(context);
@@ -91,7 +91,7 @@ void alpha_pinene_ess(char *conf_file)
 
     read_data_file(ALPHA_PINENE_DATA, &texp, &yexp);
 
-    N_Vector initial_values = copy_matrix_row(yexp, 0, 0, SM_COLUMNS_D(yexp));
+    N_Vector initial_values = copy_matrix_column(yexp, 0, 0, SM_ROWS_D(yexp));
 
     N_Vector xbest = execute_ess_solver(conf_file, OUPUT_PATH, texp, yexp, initial_values, NULL);
 
@@ -100,8 +100,7 @@ void alpha_pinene_ess(char *conf_file)
         sunrealtype expected = expected_values[i];
         sunrealtype result = NV_Ith_S(xbest, i);
 
-        sunrealtype a = 6;
-        // assert(fabs(result - expected) < 0.1);
+        assert(fabs(result - expected) < 1);
     }
 
     destroy_cuqdyn_context(context);
@@ -111,24 +110,23 @@ void logistic_model_ess(char *conf_file)
 {
     CuqDynContext context = init_cuqdyn_context_from_file("data/logistic_model_cuqdyn_config.xml");
 
-    sunrealtype expected_values[2] = {0.1, 100};
+    sunrealtype expected_values[2] = {0.1, 102};
 
     N_Vector texp = NULL;
     SUNMatrix yexp = NULL;
 
     read_data_file(LOGISTIC_MODEL_DATA, &texp, &yexp);
 
-    N_Vector initial_values = copy_matrix_row(yexp, 0, 0, SM_COLUMNS_D(yexp));
+    N_Vector initial_values = copy_matrix_column(yexp, 0, 0, SM_ROWS_D(yexp));
 
     N_Vector xbest = execute_ess_solver(conf_file, OUPUT_PATH, texp, yexp, initial_values, NULL);
 
-    for (int i = 0; i < sizeof(expected_values); ++i)
+    for (int i = 0; i < 2; ++i)
     {
         sunrealtype expected = expected_values[i];
         sunrealtype result = NV_Ith_S(xbest, i);
 
-        sunrealtype a = 6;
-        // assert(fabs(result - expected) < 0.1);
+        assert(fabs(result - expected) < 1);
     }
 
     destroy_cuqdyn_context(context);