assigning J R and g

Author: ZocoLini

modules/cuqdyn-c/include/functions.h

@@ -5,7 +5,8 @@
 
 /// Function used to solve the ODE using cvodes
 int ode_model_fun(sunrealtype t, N_Vector y, N_Vector ydot, void *data);
-/// Objetive function for the lotka_volterra problem used by the sacess library
+/// Objetive functions used by the sacess library
 void *obj_func(double *x, void *data);
+void *obj_func2(double *x, void *data);
 
 #endif // LOTKA_VOLTERRA_H

modules/cuqdyn-c/src/ess_solver.c

@@ -105,7 +105,7 @@ N_Vector execute_ess_solver(const char *file, const char *path, N_Vector texp, S
         }
     }
 
-    execute_Solver(exptotal, &result, obj_func);
+    execute_Solver(exptotal, &result, obj_func2);
 
     N_Vector predicted_params = New_Serial(exptotal[0].test.bench.dim);
     destroyexp(exptotal);

modules/cuqdyn-c/src/functions.c

@@ -2,13 +2,15 @@
 #include <method_module/structure_paralleltestbed.h>
 #include <nvector/nvector_serial.h>
 #include <stdio.h>
+#include <stdlib.h>
 #include <string.h>
 #include <sunmatrix/sunmatrix_dense.h>
 
 #include "config.h"
 #include "cuqdyn.h"
 #include "ode_solver.h"
 #include "states_transformer.h"
+#include "sundials/sundials_types.h"
 
 extern void eval_f_exprs(sunrealtype t, sunrealtype *y, sunrealtype *ydot, sunrealtype *params, CuqDynContext context);
 
@@ -62,7 +64,6 @@ void *obj_func(double *x, void *data)
 
     const long rows = SM_ROWS_D(result);
     const long cols = SM_COLUMNS_D(result);
-    sunrealtype J = 0.0;
 
     if (SM_ROWS_D(exptotal->yexp) != rows)
     {
@@ -78,19 +79,111 @@ void *obj_func(double *x, void *data)
         exit(-1);
     }
 
-    for (long i = 0; i < rows; ++i)
+    double *J = malloc(sizeof(double));
+    *J = 0.0;
+    double *g = malloc(sizeof(double));
+    *g = 0.0;
+    double *R = malloc(sizeof(double) * rows * cols);
+
+    for (long j = 0; j < cols; ++j)
     {
-        for (long j = 0; j < cols; ++j)
+        for (long i = 0; i < rows; ++i)
         {
             const sunrealtype diff = SM_ELEMENT_D(result, i, j) - SM_ELEMENT_D(exptotal->yexp, i, j);
-            J += diff * diff;
+            R[j * rows + i] = diff;
+            *J += diff * diff;
         }
     }
 
-    res->value = J;
+    res->size_j = 1;
+    res->J = J;
+    res->value = *J;
+    res->size_r = rows * cols;
+    res->R = R;
+    res->g = g;
 
     N_VDestroy(parameters);
     SUNMatDestroy(result);
 
     return res;
 }
+
+void *obj_func2(double *x, void *data)
+{
+    CuqdynConf *conf = get_cuqdyn_conf(get_cuqdyn_context());
+    experiment_total *exptotal = data;
+    output_function *res = calloc(1, sizeof(output_function));
+
+    N_Vector parameters = New_Serial(conf->ode_expr.p_count);
+    memcpy(NV_DATA_S(parameters), x, NV_LENGTH_S(parameters) * sizeof(sunrealtype));
+
+    N_Vector texp = New_Serial(NV_LENGTH_S(exptotal->texp));
+    memcpy(NV_DATA_S(texp), NV_DATA_S(exptotal->texp), NV_LENGTH_S(exptotal->texp) * sizeof(sunrealtype));
+
+#ifdef DEBUG
+    fprintf(stdout, "[DEBUG] [OBJ FUNC] ");
+    fprintf(stdout, "Params: ");
+    for (int i = 0; i < NV_LENGTH_S(parameters); i++)
+    {
+        fprintf(stdout, "%f ", NV_Ith_S(parameters, i));
+    }
+    fprintf(stdout, "\n");
+#endif
+
+    const sunrealtype t0 = NV_Ith_S(texp, 0);
+
+    TransposedStates ode_solution = solve_ode(parameters, exptotal->initial_values, t0, texp);
+    ObservablesTransposedStates result = transform_states(ode_solution);
+
+    const long rows = SM_ROWS_D(result);
+    const long cols = SM_COLUMNS_D(result);
+
+    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);
+    }
+
+    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);
+    }
+
+    double *J = malloc(sizeof(double));
+    *J = 0.0;
+    double *g = malloc(sizeof(double));
+    *g = 0.0;
+    double *R = malloc(sizeof(double) * rows * cols);
+
+    for (long j = 0; j < cols; ++j)
+    {
+        double col_mean = 0.0;
+        for (long i = 0; i < rows; ++i)
+        {
+            col_mean += SM_ELEMENT_D(exptotal->yexp, i, j);
+        }
+        col_mean /= rows;
+
+        for (long i = 0; i < rows; ++i)
+        {
+            const sunrealtype diff = SM_ELEMENT_D(result, i, j) - SM_ELEMENT_D(exptotal->yexp, i, j);
+            R[j * rows + i] = diff / col_mean;
+            *J += diff * diff;
+        }
+    }
+
+    res->size_j = 1;
+    res->J = J;
+    res->value = *J*100;
+    res->size_r = rows * cols;
+    res->R = R;
+    res->g = g;
+
+    N_VDestroy(parameters);
+    SUNMatDestroy(result);
+
+    return res;
+}