jacobians logic and config implemented but failing to converge

Author: ZocoLini

README.md

@@ -159,6 +159,7 @@ There are three types of input files that must be provided:
 
   - **tolerances:** rtol and atol used by the cvodes library.
   - **ode_expr:** ODE model expression or identifier.
+  - **jacobian_provider:** Jacobian provider defined in Rust 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
@@ -191,6 +192,7 @@ There are three types of input files that must be provided:
             p1 * y11 * y7 - p24 * p22 * y14
             p28 + p27 * y7 - p29 * y15
         </ode_expr>
+        <jacobian_provider>nfkb</jacobian_provider> <!-- 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
@@ -273,3 +275,6 @@ used indicating the identifier in the XML config file like this:
 ```
 
 Note that y_count and p_count are still present.
+
+The software also provides a way to define a Jacobian provider and a states transformer inside the 
+Rust module. See the respective Rust modules `jacobians.rs` and `states_transformers.rs`.

modules/cuqdyn-c/include/config.h

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

modules/cuqdyn-c/include/functions.h

@@ -1,11 +1,14 @@
 #ifndef LOTKA_VOLTERRA_H
 #define LOTKA_VOLTERRA_H
-#include <sundials/sundials_direct.h>
 #include <sundials/sundials_nvector.h>
+#include "sundials/sundials_matrix.h"
 
 /// 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
+/// Default obj function used by the sacess library
 void *obj_func(double *x, void *data);
+/// Jacobian function used to solve the ODE using cvodes
+int get_jac(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2,
+            N_Vector tmp3);
 
 #endif // LOTKA_VOLTERRA_H

modules/cuqdyn-c/src/functions.c

@@ -9,8 +9,11 @@
 #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);
+extern void jacobian(sunrealtype t, sunrealtype *y, sunrealtype *fy, sunrealtype *J, sunrealtype *params,
+                     sunrealtype *tmp1, sunrealtype *tmp2, sunrealtype *tmp3, CuqDynContext context);
 
 int ode_model_fun(sunrealtype t, N_Vector y, N_Vector ydot, void *user_data)
 {
@@ -94,3 +97,24 @@ void *obj_func(double *x, void *data)
 
     return res;
 }
+
+int get_jac(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2,
+            N_Vector tmp3)
+{
+#ifndef JE
+#define JE(J, i, j) SM_ELEMENT_D((J), (i), (j))
+#endif
+
+    for (int j = 0; j < SM_COLUMNS_D(J); ++j)
+    {
+        for (int i = 0; i < SM_ROWS_D(J); ++i)
+        {
+            SM_ELEMENT_D(J, i, j) = 0.0;
+        }
+    }
+
+    jacobian(t, NV_DATA_S(y), NV_DATA_S(fy), NV_DATA_S(J), (sunrealtype *) user_data, NV_DATA_S(tmp1), NV_DATA_S(tmp2),
+             NV_DATA_S(tmp3), get_cuqdyn_context());
+
+    return 0; // éxito
+}

modules/cuqdyn-c/src/ode_solver.c

@@ -12,121 +12,8 @@
 #include "cvodes/cvodes_ls.h"
 #include "sundials/sundials_types.h"
 
-#ifndef JE
-#define JE(J,i,j) SM_ELEMENT_D((J),(i),(j))
-#endif
-
 static int check_retval(void *, const char *, int);
 
-int nfk_jac(sunrealtype t, N_Vector y, N_Vector fy,
-                       SUNMatrix J, void *user_data,
-                       N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
-{
-    (void)t; (void)fy; (void)tmp1; (void)tmp2; (void)tmp3;
-    
-        const sunrealtype *Y = NV_DATA_S(y);
-        const sunrealtype y1=Y[0], y2=Y[1], y3=Y[2], y4=Y[3], y5=Y[4],
-                       y6=Y[5], y7=Y[6], y8=Y[7], y9=Y[8], y10=Y[9],
-                       y11=Y[10], y12=Y[11], y13=Y[12], y14=Y[13], y15=Y[14];
-    
-        const sunrealtype *p = (const sunrealtype*) user_data;
-        const sunrealtype p1 = p[0],  p2 = p[1],  p3 = p[2],  p4 = p[3],  p5 = p[4],
-                       p6 = p[5],  p7 = p[6],  p8 = p[7],  p9 = p[8],  p10= p[9],
-                       p11= p[10], p12= p[11], p13= p[12], p14= p[13], p15= p[14],
-                       p16= p[15], p17= p[16], p18= p[17], p19= p[18], p20= p[19],
-                       p21= p[20], p22= p[21], p23= p[22], p24= p[23], p25= p[24],
-                       p26= p[25], p27= p[26], p28= p[27], p29= p[28];
-    
-        // Limpia la matriz (si tu versión de SUNDIALS no garantiza cero)
-        for (int j = 0; j < SM_COLUMNS_D(J); ++j) 
-        {
-            for (int i = 0; i < SM_ROWS_D(J); ++i) 
-            {
-                SM_ELEMENT_D(J,i,j) = 0.0;
-            }
-        }
-    
-        // f1
-        JE(J, 0, 0) = -p21 -p17;
-        
-        // f2
-        JE(J, 1, 0) = p17;
-        JE(J, 1, 1) = -p19 - p18*y8 - p21 - p2*y10 - p4*y13;
-        JE(J, 1, 3) = p3;
-        JE(J, 1, 4) = p5;
-        JE(J, 1, 7) = -p18*y2;
-        JE(J, 1, 9) = -p2*y2;
-        JE(J, 1, 12) = -p4*y2;
-        
-        // f3
-        JE(J, 2, 1) = p19 + p18*y8;
-        JE(J, 2, 2) = -p21;
-        JE(J, 2, 7) = p18*y2;
-    
-        // f4
-        JE(J, 3, 1) = p2*y10;
-        JE(J, 3, 3) = -p3;
-        JE(J, 3, 9) = p2*y2;
-        
-        // f5
-        JE(J, 4, 1) = p4*y13;
-        JE(J, 4, 4) = -p5;
-        JE(J, 4, 12) = p4*y2;
-        
-        // f6
-        JE(J, 5, 4) = p5;
-        JE(J, 5, 5) = -p1*y10 -p23;
-        JE(J, 5, 9) = -p1*y6;
-        JE(J, 5, 12) = p11;
-        
-        // f7
-        JE(J, 6, 5) = p23*p22;
-        JE(J, 6, 6) = -p1*y11;
-        JE(J, 6, 10) = -p1*y7;
-        
-        // f8
-        JE(J, 7, 7) = -p16;
-        JE(J, 7, 8) = p15;
-        
-        // f9
-        JE(J, 8, 6) = p12; 
-        JE(J, 8, 8) = -p14;
-        
-        // f10
-        JE(J, 9, 1) = -p2*y10;
-        JE(J, 9, 5) = -p1*y10;
-        JE(J, 9, 9) = -p2*y2 - p1*y6 - p10 - p25;
-        JE(J, 9, 10) = p26;
-        JE(J, 9, 11) = p9;
-        
-        // f11
-        JE(J, 10, 6) = -p1*y11; 
-        JE(J, 10, 9) = p25*p22; 
-        JE(J, 10, 10) = -p1*y7 - p26*p22; 
-        
-        // f12
-        JE(J, 11, 6) = p6;
-        JE(J, 11, 11) = -p8;
-        
-        // f13
-        JE(J, 12, 1) = -p4*y13; 
-        JE(J, 12, 5) = p1*y10;
-        JE(J, 12, 9) = p1*y6;
-        JE(J, 12, 12) = -p11 -p4*y2;
-        JE(J, 12, 13) = p24;
-        
-        // f14
-        JE(J, 13, 6) = p1*y11; 
-        JE(J, 13, 10) = p1*y7;
-        JE(J, 13, 13) = -p24*p22;
-        
-        // f15
-        JE(J, 14, 6) = p27;
-        JE(J, 14, 14) = -p29;
-        
-        return 0; // éxito
-}
-
 
 TransposedStates solve_ode(N_Vector parameters, N_Vector initial_values, sunrealtype t0, N_Vector times)
 {
@@ -170,10 +57,13 @@ TransposedStates solve_ode(N_Vector parameters, N_Vector initial_values, sunreal
         return NULL;
     }
 
-    retval = CVodeSetJacFn(cvode_mem, nfk_jac);
-    if (check_retval(&retval, "CVodeSetJacFn", 1))
+    if (cuqdyn_conf->jacobian_provider_present) 
     {
-        return NULL;
+        retval = CVodeSetJacFn(cvode_mem, get_jac);
+        if (check_retval(&retval, "CVodeSetJacFn", 1))
+        {
+            return NULL;
+        }
     }
 
     char *min_step_s = getenv("SACESS_CVODES_MIN_STEP");