Add CVODES forward sensitivity integration, fix FD bump default, and merge C prototype fixes

.gitignore

@@ -100,3 +100,5 @@ bng2/Network3/cmake_install.cmake
 statfactor_*.txt
 bng2/Network3/run_network
 *.o
+.jules/
+patch_*.py

CMakeLists.txt

@@ -57,9 +57,12 @@ FetchContent_Declare(
   GIT_SUBMODULES ""  # Skip submodules (not needed for CVODE build)
 )
 
-# Configure SUNDIALS v7 build options - only build CVODE, no tests/examples
+# Configure SUNDIALS v7 build options - only build CVODE(S), no tests/examples
 set(BUILD_CVODE ON CACHE INTERNAL "Build CVODE solver")
-set(BUILD_CVODES OFF CACHE INTERNAL "")
+# CVODES is a superset of CVODE (adds forward/adjoint sensitivity analysis,
+# CVodeSensInit/CVodeGetSens/CVodeAdjInit/etc.) needed for gradient-based
+# parameter estimation. Link against it instead of plain CVODE.
+set(BUILD_CVODES ON CACHE INTERNAL "Build CVODES solver (sensitivity analysis)")
 set(BUILD_IDA OFF CACHE INTERNAL "")
 set(BUILD_IDAS OFF CACHE INTERNAL "")
 set(BUILD_KINSOL OFF CACHE INTERNAL "")

bng-graph/nauty/nauty24/addedgeg.c

@@ -259,7 +259,7 @@ main(int argc, char *argv[])
 #if !MAXN
 		        DYNALLOC2(graph,h,h_sz,n,m,"addedgeg");
 #endif
-	 	        fcanonise(g,m,n,h,NULL,FALSE);  /*FIXME (loops)*/
+		        fcanonise(g,m,n,h,NULL,loopcount(g,m,n)>0);
 		        gq = h;
 	            }
 	            if (outcode == SPARSE6) writes6(outfile,gq,m,n);

bng-graph/nauty/nauty24/deledgeg.c

@@ -14,6 +14,7 @@
 /*************************************************************************/
 
 #include "gtools.h" 
+#include "gutils.h"
 
 /**************************************************************************/
 
@@ -159,7 +160,7 @@ main(int argc, char *argv[])
 #if !MAXN
 		        DYNALLOC2(graph,h,h_sz,n,m,"deledgeg");
 #endif
-	 	        fcanonise(g,m,n,h,NULL,FALSE);  /* FIXME (loops) */
+		        fcanonise(g,m,n,h,NULL,loopcount(g,m,n)>0);
 		        gq = h;
 	            }
 	            if (outcode == SPARSE6) writes6(outfile,gq,m,n);

bng-graph/nauty/nauty24/newedgeg.c

@@ -110,7 +110,7 @@ na_newedge(graph *g1, int m1, int n1, boolean dolabel)
 
                 if (dolabel)
                 {
-                    fcanonise(g2,m2,n2,h,NULL,FALSE);  /* FIXME (loops) */
+                    fcanonise(g2,m2,n2,h,NULL,loopcount(g2,m2,n2)>0);
                     gq = h;
                 }
                 if (outcode == SPARSE6) writes6(outfile,gq,m2,n2);

bng-graph/nauty/nauty24/testg.c

@@ -2,7 +2,7 @@
    both pickg (select by property) and countg (count by property).
    Version of Nov 19, 2003. */
 /* TODO - write a header if input has one */
-/* TODO - USERDEF should be long, not int */
+
 
 #define USAGE \
   "[pickg|countg] [-fp#:#q -V] [--keys] [-constraints -v] [ifile [ofile]]"
@@ -74,11 +74,11 @@ External user-defined parameters:
   case the parameter is selected using the letter 'Q'.  The name of the
   parameter is "userdef" unless USERDEFNAME is defined.  The function
   is called with the parameters (graph *g, int m, int n) and must return
-  an integer value.
+  a long value.
 */
 
 #ifdef USERDEF
-int USERDEF(graph*,int,int);
+long USERDEF(graph*,int,int);
 #endif
 #ifndef USERDEFNAME
 #define USERDEFNAME "userdef"

bng2/Network3/src/model/rateExpressions/rateElementary.cpp

@@ -100,19 +100,20 @@ double RateElementary::get_dRate_dX(unsigned int which, vector<double> X){
 		}
 		else{
 			double dX_which = 0.0;
+			double inv_denom = 1.0;
+			for (int j=0;j < stoich;j++){
+				inv_denom /= ((double)j+1.0);
+			}
 			for (int k=0;k < stoich;k++){ // # of terms in summation
 				double prod = 1.0;
 				for (int j=0;j < stoich;j++){ // # of terms in each product of the summation
-					if (j == k){
-						prod *= 1.0/((double)j+1.0);
-					}
-					else{
-						prod *= (X[i]-(double)j)/((double)j+1.0);
+					if (j != k){
+						prod *= (X[i]-(double)j);
 					}
 				}
 				dX_which += prod;
 			}
-			dRate *= dX_which;
+			dRate *= (dX_which * inv_denom);
 		}
 	}
 	return dRate;

bng2/Network3/src/util/mathutils/linsolve.c

@@ -1,38 +1,38 @@
-#include "mathutils.h"
-
-
-/*==========================================================================*/
-/* Routine to Solve set of linear equations a.x=b where a is symmetric. */
-/* Note that the rhs b are stored in the ROWS of b */
-/* rather than the columns because of the way Fortran accesses arrays. */
-
-
-int LINSOLVE(double **a, double **b, int nrhs, int dim){
-#ifdef APPEND_UNDERSCORE
-#define DGESV dgesv_
-#else
-#define DGESV dgesv
-#endif
-    register int i,j;
-    int *ipiv;
-    int Nrhs, Dim;
-    int info;
-#ifdef __cplusplus
-extern "C" {
-#endif
-    /* Fortran dgesv: SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO ) */
-    void DGESV(int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *ldb, int *info);
-#ifdef __cplusplus
-}
-#endif
-
-    /* calculate transpose of a-- this is done because matrices in Fortran */
-    /* and C have row and column indices interchanged */
-    MATRIX_TRANSPOSE(a,dim);
-    Dim= dim;
-    Nrhs= nrhs;
-    ipiv = (int *) malloc (sizeof(int)*Dim);
-    DGESV(&Dim, &Nrhs, a[0], &Dim, ipiv, b[0], &Dim, &info);
-    free(ipiv);
-    return(info);
-}
+#include "mathutils.h"
+
+		    
+/*==========================================================================*/
+/* Routine to Solve set of linear equations a.x=b where a is symmetric. */
+/* Note that the rhs b are stored in the ROWS of b */
+/* rather than the columns because of the way Fortran accesses arrays. */
+
+
+int LINSOLVE(double **a, double **b, int nrhs, int dim){
+#ifdef APPEND_UNDERSCORE
+#define DGESV dgesv_
+#else
+#define DGESV dgesv
+#endif
+    register int i,j;
+    int *ipiv;
+    int Nrhs, Dim;
+    int info;
+#ifdef __cplusplus
+extern "C" {
+#endif
+    /* Fortran dgesv: SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO ) */
+    void DGESV(int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *ldb, int *info);
+#ifdef __cplusplus
+}
+#endif
+    
+    /* calculate transpose of a-- this is done because matrices in Fortran */
+    /* and C have row and column indices interchanged */
+    MATRIX_TRANSPOSE(a,dim);
+    Dim= dim;
+    Nrhs= nrhs;
+    ipiv = (int *) malloc (sizeof(int)*Dim);
+    DGESV(&Dim, &Nrhs, a[0], &Dim, ipiv, b[0], &Dim, &info);
+    free(ipiv);
+    return(info);
+}

bng2/Network3/src/util/mathutils/matrix_mult.c

@@ -1,23 +1,23 @@
-#include "mathutils.h"
-
-
-/*==========================================================================*/
-void MATRIX_MULT(double **a, double **b, double **c, int dim){
-  /* uses BLAS3 routine dgemm (NAG f06yaf) */
-#ifdef APPEND_UNDERSCORE
-#define DGEMM dgemm_
-#else
-#define DGEMM dgemm
-#endif
-  int n=dim;
-  double alpha=1.0;
-  double beta=0.0;
-  extern void DGEMM(char*, char*, int*, int*, int*, double*, double*, int*, double*, int*, double*, double*, int*);
-
-  /* need to call dgemm_ with transposed matrices because C and Fortran */
-  /* reference matrices in the opposite sense */
-  DGEMM("t", "t", &n, &n, &n, &alpha, a[0], &n, b[0], &n, &beta, c[0], &n);
-  /* tranpose result */
-  MATRIX_TRANSPOSE(c, dim);
-}
-
+#include "mathutils.h"
+
+		    
+/*==========================================================================*/
+void MATRIX_MULT(double **a, double **b, double **c, int dim){
+  /* uses BLAS3 routine dgemm (NAG f06yaf) */
+#ifdef APPEND_UNDERSCORE
+#define DGEMM dgemm_
+#else
+#define DGEMM dgemm
+#endif
+  int n=dim;
+  double alpha=1.0;
+  double beta=0.0;
+  extern void DGEMM(char*, char*, int*, int*, int*, double*, double*, int*, double*, int*, double*, double*, int*);
+
+  /* need to call dgemm_ with transposed matrices because C and Fortran */
+  /* reference matrices in the opposite sense */
+  DGEMM("t", "t", &n, &n, &n, &alpha, a[0], &n, b[0], &n, &beta, c[0], &n);
+  /* tranpose result */
+  MATRIX_TRANSPOSE(c, dim);
+}
+

bng2/Network3/src/util/mathutils/matrix_vector_mult.c

@@ -1,23 +1,23 @@
-#include "mathutils.h"
-
-
-/*==========================================================================*/
-void MATRIX_VECTOR_MULT(double **a, double *x, double *y, int dim){
-    /* uses BLAS3 routine dgemv  */
-#ifdef APPEND_UNDERSCORE
-#define DGEMV dgemv_
-#else
-#define DGEMV dgemv
-#endif
-    int n=dim;
-    int incx=1;
-    int incy=1;
-    double alpha=1.0;
-    double beta=0.0;
-    extern void DGEMV(char*, int*, int*, double*, double*, int*, double*, int*, double*, double*, int*);
-
-    /* take the transpose of the multiplying matrix because matrices */
-    /* referenced opposite in Fortran */
-    DGEMV("t", &n, &n, &alpha, a[0], &n, x, &incx, &beta, y, &incy);
-}
-
+#include "mathutils.h"
+
+		    
+/*==========================================================================*/
+void MATRIX_VECTOR_MULT(double **a, double *x, double *y, int dim){
+    /* uses BLAS3 routine dgemv  */
+#ifdef APPEND_UNDERSCORE
+#define DGEMV dgemv_
+#else
+#define DGEMV dgemv
+#endif
+    int n=dim;
+    int incx=1;
+    int incy=1;
+    double alpha=1.0;
+    double beta=0.0;
+    extern void DGEMV(char*, int*, int*, double*, double*, int*, double*, int*, double*, double*, int*);
+    
+    /* take the transpose of the multiplying matrix because matrices */
+    /* referenced opposite in Fortran */
+    DGEMV("t", &n, &n, &alpha, a[0], &n, x, &incx, &beta, y, &incy);
+}
+