new states_transformer struct added to the cuqdyn config and updated the rust fragment to be able to use it

Author: ZocoLini

example-files/nfkb_cuqdyn_config.xml

@@ -22,4 +22,12 @@
         p1 * y11 * y7 - p24 * p22 * y14
         p28 + p27 * y7 - p29 * y15
     </ode_expr>
+    <states_transformer count="6">
+        y7
+        y10 + y13
+        y9
+        y1 + y2 + y3
+        y2
+        y12
+    </states_transformer>
 </cuqdyn-config>

modules/cuqdyn-c/include/config.h

@@ -5,8 +5,9 @@
 
 typedef struct
 {
-    sunrealtype rtol;
-    N_Vector atol;
+    double rtol;
+    int atol_len;
+    double *atol;
 
 } Tolerances;
 
@@ -23,15 +24,25 @@ typedef struct
 OdeExpr create_ode_expr(int y_count, int p_count, char** exprs);
 void destroy_ode_expr(OdeExpr ode_expr);
 
+typedef struct
+{
+    int count;
+    char** exprs;
+} StatesTransformer;
+
+StatesTransformer create_states_transformer(int count, char** exprs);
+void destroy_states_transformer(StatesTransformer observables);
+
 typedef struct
 {
     Tolerances tolerances;
     OdeExpr ode_expr;
+    StatesTransformer states_transformer;
 } CuqdynConf;
 
 CuqdynConf *init_cuqdyn_conf_from_file(const char *filename);
 int parse_cuqdyn_conf(const char* filename, CuqdynConf* config);
-CuqdynConf *init_cuqdyn_conf(Tolerances tolerances, OdeExpr ode_expr);
+CuqdynConf *init_cuqdyn_conf(Tolerances tolerances, OdeExpr ode_expr, StatesTransformer observables);
 void destroy_cuqdyn_conf();
 CuqdynConf * get_cuqdyn_conf();
 

modules/cuqdyn-c/include/functions.h

@@ -5,7 +5,7 @@
 
 #include "config.h"
 
-void mexpreval_init_wrapper(OdeExpr ode_expr);
+void mexpreval_init_wrapper(CuqdynConf cuqdyn_conf);
 /// 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

modules/cuqdyn-c/src/config.c

@@ -13,11 +13,12 @@
 
 static CuqdynConf *config = NULL;
 
-CuqdynConf *create_cuqdyn_conf(Tolerances tolerances, OdeExpr ode_expr)
+CuqdynConf *create_cuqdyn_conf(Tolerances tolerances, OdeExpr ode_expr, StatesTransformer states_transformer)
 {
     CuqdynConf *cuqdyn_conf = malloc(sizeof(CuqdynConf));
     cuqdyn_conf->tolerances = tolerances;
     cuqdyn_conf->ode_expr = ode_expr;
+    cuqdyn_conf->states_transformer = states_transformer;
     return cuqdyn_conf;
 }
 
@@ -30,6 +31,7 @@ void destroy_cuqdyn_conf()
 
     destroy_tolerances(config->tolerances);
     destroy_ode_expr(config->ode_expr);
+    destroy_states_transformer(config->states_transformer);
     free(config);
     config = NULL;
 }
@@ -39,8 +41,7 @@ CuqdynConf *init_cuqdyn_conf_from_file(const char *filename)
     CuqdynConf *tmp_config = malloc(sizeof(CuqdynConf));
     if (tmp_config == NULL)
     {
-        fprintf(stderr, "ERROR: Memory allocation failed in function "
-                        "init_cuqdyn_conf_from_file()\n");
+        fprintf(stderr, "ERROR: Memory allocation failed in function init_cuqdyn_conf_from_file()\n");
         exit(1);
     }
 
@@ -53,7 +54,7 @@ CuqdynConf *init_cuqdyn_conf_from_file(const char *filename)
     }
 
     config = tmp_config;
-    mexpreval_init_wrapper(config->ode_expr);
+    mexpreval_init_wrapper(*config);
     return config;
 }
 
@@ -71,10 +72,14 @@ int parse_cuqdyn_conf(const char *filename, CuqdynConf *config)
 
     sunrealtype rtol = 1e-6;
     N_Vector atol = NULL;
+
     int y_count = 0;
     char **odes = NULL;
     int p_count = 0;
 
+    int obs_count = 0;
+    char **obs_tranformations = NULL;
+
     for (; cur; cur = cur->next)
     {
         if (cur->type != XML_ELEMENT_NODE)
@@ -185,19 +190,61 @@ int parse_cuqdyn_conf(const char *filename, CuqdynConf *config)
                 token = strtok(NULL, "\n");
             }
         }
+        else if (!xmlStrcmp(cur->name, "states_transformer"))
+        {
+            xmlChar *count_attr = xmlGetProp(cur, "count");
+            if (count_attr == NULL)
+            {
+                fprintf(stderr, "Error: <states_transformer> node does not have a 'count' attribute\n");
+                xmlFreeDoc(doc);
+                return 1;
+            }
+
+            obs_count = atoi((char *) count_attr);
+            xmlFree(count_attr);
+
+            xmlChar *content = xmlNodeGetContent(cur);
+            if (content == NULL)
+            {
+                fprintf(stderr, "Error: no content in <states_transformer> node\n");
+                xmlFreeDoc(doc);
+                return 1;
+            }
+
+            char *str = (char *) content;
+
+            char *token = strtok(str, "\n");
+
+            obs_tranformations = calloc(y_count, sizeof(char *));
+
+            int index = 0;
+            while (token != NULL)
+            {
+                while (*token == ' ')
+                    token++;
+
+                if (*token != '\0')
+                {
+                    obs_tranformations[index] = token;
+                    index++;
+                }
+                token = strtok(NULL, "\n");
+            }
+        }
     }
 
     config->tolerances = create_tolerances(rtol, atol);
     config->ode_expr = create_ode_expr(y_count, p_count, odes);
+    config->states_transformer = create_states_transformer(obs_count, obs_tranformations);
     xmlFreeDoc(doc);
     return 0;
 }
 
-CuqdynConf *init_cuqdyn_conf(Tolerances tolerances, OdeExpr ode_expr)
+CuqdynConf *init_cuqdyn_conf(Tolerances tolerances, OdeExpr ode_expr, StatesTransformer states_transformer)
 {
     destroy_cuqdyn_conf();
 
-    config = create_cuqdyn_conf(tolerances, ode_expr);
+    config = create_cuqdyn_conf(tolerances, ode_expr, states_transformer);
     return config;
 }
 
@@ -217,11 +264,13 @@ Tolerances create_tolerances(sunrealtype scalar_rtol, N_Vector atol)
 {
     Tolerances tolerances;
     tolerances.rtol = scalar_rtol;
-    tolerances.atol = atol;
+    tolerances.atol_len = NV_LENGTH_S(atol);
+    tolerances.atol = malloc(tolerances.atol_len * sizeof(double));
+    memcpy(tolerances.atol, NV_DATA_S(atol), tolerances.atol_len * sizeof(double));
     return tolerances;
 }
 
-void destroy_tolerances(Tolerances tolerances) { N_VDestroy(tolerances.atol); }
+void destroy_tolerances(Tolerances tolerances) { free(tolerances.atol); }
 
 OdeExpr create_ode_expr(int y_count, int p_count, char **exprs)
 {
@@ -233,3 +282,16 @@ OdeExpr create_ode_expr(int y_count, int p_count, char **exprs)
 }
 
 void destroy_ode_expr(OdeExpr ode_expr) {}
+
+StatesTransformer create_states_transformer(int count, char** exprs)
+{
+    StatesTransformer states_transformer;
+    states_transformer.count = count;
+    states_transformer.exprs = exprs;
+    return states_transformer;
+}
+
+void destroy_states_transformer(StatesTransformer states_transformer)
+{
+
+}

modules/cuqdyn-c/src/functions.c

@@ -7,11 +7,11 @@
 #include "config.h"
 #include "cuqdyn.h"
 
-extern void mexpreval_init(OdeExpr ode_expr);
+extern void mexpreval_init(CuqdynConf cuqdyn_conf);
 
-void mexpreval_init_wrapper(OdeExpr ode_expr)
+void mexpreval_init_wrapper(CuqdynConf cuqdyn_conf)
 {
-    mexpreval_init(ode_expr);
+    mexpreval_init(cuqdyn_conf);
 }
 
 extern void eval_f_exprs(sunrealtype t, sunrealtype *y, sunrealtype *ydot, sunrealtype *params);
@@ -60,6 +60,18 @@ void *ode_model_obj_func(double *x, void *data)
 
     sunrealtype J = 0.0;
 
+    if (SM_ROWS_D(exptotal->yexp) != rows)
+    {
+        fprintf(stderr, "ERROR: The yexp rows don't match the ode result rows");
+        exit(-1);
+    }
+
+    if (SM_COLUMNS_D(exptotal->yexp) != cols - 1)
+    {
+        fprintf(stderr, "ERROR: The yexp cols don't match the ode result cols");
+        exit(-1);
+    }
+
     for (long i = 0; i < rows; ++i)
     {
         // Note that the first col of the result matrix is t

modules/cuqdyn-c/src/ode_solver.c

@@ -22,8 +22,8 @@ SUNMatrix solve_ode(N_Vector parameters, N_Vector initial_values, sunrealtype t0
     retval = CVodeInit(cvode_mem, ode_model_fun, t0, initial_values);
     if (check_retval(&retval, "CVodeInit", 1)) { return NULL; }
 
-    N_Vector cloned_abs_tol = New_Serial(NV_LENGTH_S(tolerances.atol));
-    memcpy(NV_DATA_S(cloned_abs_tol), NV_DATA_S(tolerances.atol), NV_LENGTH_S(tolerances.atol) * sizeof(sunrealtype));
+    N_Vector cloned_abs_tol = New_Serial(tolerances.atol_len);
+    memcpy(NV_DATA_S(cloned_abs_tol), tolerances.atol, tolerances.atol_len * sizeof(sunrealtype));
 
     // We clone the tolerances because the CVodeFree function frees the memory allocated for the abs_tol it receives
     retval = CVodeSVtolerances(cvode_mem, tolerances.rtol, cloned_abs_tol);

modules/mexpreval/Cargo.toml

@@ -5,6 +5,7 @@ edition = "2018"
 
 [dependencies]
 meval = { git = "https://github.com/ZocoLini/mexpreval-rs.git", tag = "v0.2.2" }
+getset = { version = "0.1.5" }
 
 [dev-dependencies]
 criterion = { version = "0.5.1" }

modules/mexpreval/benches/eval_f_exprs_bench.rs

@@ -1,23 +1,15 @@
 use criterion::{criterion_group, criterion_main, Criterion};
-use mexpreval::{eval_f_exprs, mexpreval_init, OdeExpr};
-use std::ffi::CString;
-use std::os::raw::c_char;
+use mexpreval::config::CuqdynConfig;
+use mexpreval::{eval_f_exprs, mexpreval_init};
 
 fn lotka_volterra_bench_eval(c: &mut Criterion) {
     let num_exprs = 2;
-    let num_params = 4;
 
     let mut y = vec![1.0, 2.0];
     let mut ydot = vec![0.0; num_exprs];
     let mut params = vec![0.1, 0.2, 0.3, 0.4];
 
-    let expr_strings = [
-        CString::new("y1 * (p1 - p2 * y2)").unwrap(),
-        CString::new("-y2 * (p3 - p4 * y1)").unwrap(),
-    ];
-    let expr_ptrs: Vec<*const c_char> = expr_strings.iter().map(|s| s.as_ptr()).collect();
-    
-    unsafe { mexpreval_init(OdeExpr::new(num_exprs as i32, num_params, expr_ptrs.as_ptr())); }
+    unsafe { mexpreval_init(CuqdynConfig::lotka_volterra_expr().into()); }
 
     c.bench_function("lotka_volterra", |b| {
         b.iter(|| unsafe {
@@ -34,18 +26,12 @@ fn lotka_volterra_bench_eval(c: &mut Criterion) {
 #[allow(dead_code)]
 fn lotka_volterra_predefined_bench_eval(c: &mut Criterion) {
     let num_exprs = 2;
-    let num_params = 4;
 
     let mut y = vec![1.0, 2.0];
     let mut ydot = vec![0.0; num_exprs];
     let mut params = vec![0.1, 0.2, 0.3, 0.4];
 
-    let expr_strings = [
-        CString::new("lotka-volterra").unwrap(),
-    ];
-    let expr_ptrs: Vec<*const c_char> = expr_strings.iter().map(|s| s.as_ptr()).collect();
-
-    unsafe { mexpreval_init(OdeExpr::new(num_exprs as i32, num_params, expr_ptrs.as_ptr())); }
+    unsafe { mexpreval_init(CuqdynConfig::lotka_volterra().into()); }
 
     c.bench_function("lotka_volterra_predefined", |b| {
         b.iter(|| unsafe {
@@ -61,18 +47,12 @@ fn lotka_volterra_predefined_bench_eval(c: &mut Criterion) {
 
 fn logistic_model_bench_eval(c: &mut Criterion) {
     let num_exprs = 1;
-    let num_params = 2;
 
     let mut y = vec![0.0];
     let mut ydot = vec![0.0; num_exprs];
     let mut params = vec![0.1, 100.0];
 
-    let expr_strings = [
-        CString::new("p1 * y1 * (1 - y1 / p2)").unwrap()
-    ];
-    let expr_ptrs: Vec<*const c_char> = expr_strings.iter().map(|s| s.as_ptr()).collect();
-
-    unsafe { mexpreval_init(OdeExpr::new(num_exprs as i32, num_params, expr_ptrs.as_ptr())); }
+    unsafe { mexpreval_init(CuqdynConfig::logistic_growth_expr().into()); }
 
     c.bench_function("logistic_model", |b| {
         b.iter(|| unsafe {
@@ -88,22 +68,12 @@ fn logistic_model_bench_eval(c: &mut Criterion) {
 
 fn alpha_pinene_bench_eval(c: &mut Criterion) {
     let num_exprs = 5;
-    let num_params = 5;
 
     let mut y = vec![1.0, 1.0, 1.0, 1.0, 1.0];
     let mut ydot = vec![0.0; num_exprs];
     let mut params = vec![0.1, 0.2, 0.2, 0.2, 0.2];
 
-    let expr_strings = [
-        CString::new("-(p1 + p2) * y1").unwrap(),
-        CString::new("p1 * y1").unwrap(),
-        CString::new("p2 * y1 - (p3 + p4) * y3 + p5 * y5").unwrap(),
-        CString::new("p3 * y3").unwrap(),
-        CString::new("p4 * y3 - p5 * y5").unwrap(),
-    ];
-    let expr_ptrs: Vec<*const c_char> = expr_strings.iter().map(|s| s.as_ptr()).collect();
-
-    unsafe { mexpreval_init(OdeExpr::new(num_exprs as i32, num_params, expr_ptrs.as_ptr())); }
+    unsafe { mexpreval_init(CuqdynConfig::alpha_pinene_expr().into()); }
 
     c.bench_function("alpha_pinene", |b| {
         b.iter(|| unsafe {