time sclaing parameter added to the config xml

Author: ZocoLini

README.md

@@ -157,14 +157,14 @@ There are three types of input files that must be provided:
 - **cuqdyn config xml:**
   This file contains the configuration of the cuqdyn solver used in the cuqdyn-c library.
 
-  - The tolerances block defines the rtol and atol used by the cvodes library.
-  - The ode model is defined inside the ode_expr block.
-  - y0 is the initial conditions of the ODE (If no present, the first row of the data matrix will be used) (Optional).
-  - The states_transformer block defines the transformations applied to the different states in case of the observed data
-  being a combination of the different states.
+  - **tolerances:** rtol and atol used by the cvodes library.
+  - **ode_expr:** ODE model expression or identifier.
+  - **time_scaling:** Scaling factor for time. Using negative powers of 10 helps cvodes speed.
+  - **y0:** Initial conditions of the ODE.
+  - **states_transformer:** Transformations expressions or identidier applied to the different states in case of the
+  observed data being a combination of the different states.
 
-  There is an option to accelerate the process of evaluating the ODE by defining it and the states transformer inside the
-  mevalexpr module. We will talk about this later.
+  There is an option to accelerate the process of evaluating the ODE by defining it and the states transformer inside the mevalexpr module. We will talk about this later.
 
     ```xml
     <?xml version="1.0" encoding="UTF-8" ?>
@@ -191,7 +191,8 @@ There are three types of input files that must be provided:
             p1 * y11 * y7 - p24 * p22 * y14
             p28 + p27 * y7 - p29 * y15
         </ode_expr>
-        <y0> <!-- Optional -->
+        <time_scaling>0.001</time_scaling> <!-- Optional (Defaults to 1.0) -->
+        <y0> <!-- Optional (Defaults to the first row of the data file) -->
             0.200000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000316, 0.002296, 0.004783, 0.000003, 0.002507, 0.003436, 0.000003, 0.060000, 0.000079, 0.000003
         </y0>
         <states_transformer count="6"> <!-- Optional -->
@@ -208,7 +209,7 @@ There are three types of input files that must be provided:
 - **Data file:**
   The data file containing a mtrix of observed data and the initial value
   needed to solve the ODE. The data file should be a txt file written with the following format:
-  
+
     ```
     # The first row is gonna be used as the initial condition if y0 is not present in the config file
     31 3 # Matrix dimensions so the parsing is easier

example-files/nfkb_cuqdyn_config.xml

@@ -8,6 +8,7 @@
     <ode_expr y_count="15" p_count="29">
         nfkb
     </ode_expr>
+    <time_scaling>0.001</time_scaling>
     <y0>
         0.200000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000316, 0.002296, 0.004783, 0.000003, 0.002507, 0.003436, 0.000003, 0.060000, 0.000079, 0.000003
     </y0>

modules/cuqdyn-c/include/config.h

@@ -32,6 +32,7 @@ typedef struct
 {
     Tolerances tolerances;
     OdeExpr ode_expr;
+    double time_scaling;
     Y0 y0;
     StatesTransformer states_transformer;
 } CuqdynConf;

modules/cuqdyn-c/src/cuqdyn.c

@@ -74,6 +74,11 @@ CuqdynResult *cuqdyn_algo(const char *data_file, const char *sacess_conf_file, c
         exit(1);
     }
 
+    for (int i = 0; i < NV_LENGTH_S(times); ++i) 
+    {
+        NV_Ith_S(times, i) = NV_Ith_S(times, i) * config->time_scaling;
+    }
+    
     const sunrealtype t0 = NV_Ith_S(times, 0);
     N_Vector initial_condition = NULL;
 
@@ -87,6 +92,8 @@ CuqdynResult *cuqdyn_algo(const char *data_file, const char *sacess_conf_file, c
         memcpy(NV_DATA_S(initial_condition), config->y0.array, config->y0.len * sizeof(sunrealtype));
     }
 
+    // Note: observed_data is a transposed matrix so the rows of the transposed 
+    // matrix are the columns of the original matrix
     const long n = SM_ROWS_D(observed_data);
     const long m = SM_COLUMNS_D(observed_data);
 

modules/cuqdyn-c/src/functions.c

@@ -14,11 +14,15 @@ extern void eval_f_exprs(sunrealtype t, sunrealtype *y, sunrealtype *ydot, sunre
 
 int ode_model_fun(sunrealtype t, N_Vector y, N_Vector ydot, void *user_data)
 {
+    CuqdynConf *conf = get_cuqdyn_conf(get_cuqdyn_context());
+
     N_Vector params_vec = user_data;
     sunrealtype *params = NV_DATA_S(params_vec);
     sunrealtype *ydot_pointer = NV_DATA_S(ydot);
     sunrealtype *y_pointer = NV_DATA_S(y);
 
+    t /= conf->time_scaling;
+
 #ifdef DEBUG
     fprintf(stdout, "[DEBUG] [ODE MODEL FUN] ");
     fprintf(stdout, "t: %f\n", t);

modules/mexpreval/src/context.rs

@@ -39,10 +39,15 @@ pub struct StatesTransformerC {
 pub struct CuqdynConfigC {
     tolerances: TolerancesC,
     ode_expr: OdeExprC,
+    time_scaling: f64,
     y0: Y0C,
     states_transformer: StatesTransformerC,
 }
 
+fn default_time_scaling() -> f64 {
+    1.0
+}
+
 #[derive(Debug, Getters, Deserialize, Clone, PartialEq)]
 #[serde(rename = "cuqdyn-config")]
 pub struct CuqdynConfigRs {
@@ -52,6 +57,9 @@ pub struct CuqdynConfigRs {
     #[serde(rename = "ode_expr")]
     ode_expr: OdeExpr,
     #[get = "pub"]
+    #[serde(rename = "time_scaling", default = "default_time_scaling")]
+    time_scaling: f64,
+    #[get = "pub"]
     #[serde(default)]
     y0: Option<F64Vec>,
     #[get = "pub"]
@@ -103,6 +111,7 @@ impl CuqdynConfigRs {
                 p_count: 4,
                 expr: StringVec(vec!["lotka-volterra".to_string()]),
             },
+            time_scaling: 1.0,
             y0: None,
             states_transformer: None,
         }
@@ -122,6 +131,7 @@ impl CuqdynConfigRs {
                     "-y2 * (p3 - p4 * y1)".to_string(),
                 ]),
             },
+            time_scaling: 1.0,
             y0: None,
             states_transformer: None,
         }
@@ -138,6 +148,7 @@ impl CuqdynConfigRs {
                 p_count: 5,
                 expr: StringVec(vec!["alpha-pinene".to_string()]),
             },
+            time_scaling: 1.0,
             y0: None,
             states_transformer: None,
         }
@@ -160,6 +171,7 @@ impl CuqdynConfigRs {
                     "p4 * y3 - p5 * y5".to_string(),
                 ]),
             },
+            time_scaling: 1.0,
             y0: None,
             states_transformer: None,
         }
@@ -176,6 +188,7 @@ impl CuqdynConfigRs {
                 p_count: 2,
                 expr: StringVec(vec!["p1 * y1 * (1 - y1 / p2)".to_string()]),
             },
+            time_scaling: 1.0,
             y0: None,
             states_transformer: None,
         }
@@ -208,6 +221,7 @@ impl CuqdynConfigRs {
                     "p28 + p27 * y7 - p29 * y15".to_string(),
                 ]),
             },
+            time_scaling: 0.001,
             y0: Some(F64Vec(vec![0.200000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000316, 0.002296, 0.004783, 0.000003, 0.002507, 0.003436, 0.000003, 0.060000, 0.000079, 0.000003])),
             states_transformer: Some(StatesTransformer {count: 6, expr: StringVec(vec![
                 "y7".to_string(),
@@ -234,6 +248,7 @@ impl CuqdynConfigRs {
                 p_count: 29,
                 expr: StringVec(vec!["nfkb".to_string()]),
             },
+            time_scaling: 0.001,
             y0: Some(F64Vec(vec![
                 0.200000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000316, 0.002296, 0.004783,
                 0.000003, 0.002507, 0.003436, 0.000003, 0.060000, 0.000079, 0.000003,
@@ -395,6 +410,7 @@ impl From<CuqdynConfigRs> for CuqdynContext {
         let c_config = CuqdynConfigC {
             tolerances,
             ode_expr,
+            time_scaling: value.time_scaling,
             y0,
             states_transformer: observables,
         };
@@ -485,6 +501,7 @@ mod tests {
     <ode_expr y_count="15" p_count="29">
         nfkb
     </ode_expr>
+    <time_scaling>0.001</time_scaling>
     <y0>
         0.200000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000316, 0.002296, 0.004783, 0.000003, 0.002507, 0.003436, 0.000003, 0.060000, 0.000079, 0.000003
     </y0>