diff --git a/.gitignore b/.gitignore index 5a006f1d..a17612a7 100644 --- a/.gitignore +++ b/.gitignore @@ -100,3 +100,5 @@ bng2/Network3/cmake_install.cmake statfactor_*.txt bng2/Network3/run_network *.o +.jules/ +patch_*.py diff --git a/.jules/sentinel.md b/.jules/sentinel.md deleted file mode 100644 index 339d55d7..00000000 --- a/.jules/sentinel.md +++ /dev/null @@ -1,4 +0,0 @@ -## 2026-05-25 - Shell Injection via `system()` in Perl Scripts -**Vulnerability:** Found `system()` calls in Perl scripts (e.g., `reformat_all.pl`, `run_all.pl`, `validate_examples.pl`) executing external commands with potentially unsanitized arguments, notably one iterating over filenames (`readdir`). This allows command injection if a maliciously crafted file exists. -**Learning:** In Perl, using `system(@args)` doesn't completely avoid the shell if the array is inadvertently collapsed or contains shell metacharacters, and `system("cp $file $newfile")` explicitly invokes a subshell. -**Prevention:** Use native Perl functions (`File::Copy::copy`) for file operations instead of spawning `cp`. For mandatory external command executions, strictly enforce the indirect object syntax `system { $args[0] } @args` which bypasses shell interpretation entirely regardless of input. diff --git a/CMakeLists.txt b/CMakeLists.txt index e69dd376..b7fb5021 100644 --- a/CMakeLists.txt +++ b/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 "") diff --git a/bng-graph/nauty/nauty24/addedgeg.c b/bng-graph/nauty/nauty24/addedgeg.c index 0213079b..994b9547 100644 --- a/bng-graph/nauty/nauty24/addedgeg.c +++ b/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); diff --git a/bng-graph/nauty/nauty24/deledgeg.c b/bng-graph/nauty/nauty24/deledgeg.c index 55e44560..1117d0e3 100644 --- a/bng-graph/nauty/nauty24/deledgeg.c +++ b/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); diff --git a/bng-graph/nauty/nauty24/newedgeg.c b/bng-graph/nauty/nauty24/newedgeg.c index 29e34bb0..94ac1159 100644 --- a/bng-graph/nauty/nauty24/newedgeg.c +++ b/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); diff --git a/bng-graph/nauty/nauty24/testg.c b/bng-graph/nauty/nauty24/testg.c index c11e5d3f..20f1b708 100644 --- a/bng-graph/nauty/nauty24/testg.c +++ b/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" diff --git a/bng2/Network3/src/model/rateExpressions/rateElementary.cpp b/bng2/Network3/src/model/rateExpressions/rateElementary.cpp index 073bb09c..4d35de26 100644 --- a/bng2/Network3/src/model/rateExpressions/rateElementary.cpp +++ b/bng2/Network3/src/model/rateExpressions/rateElementary.cpp @@ -100,19 +100,20 @@ double RateElementary::get_dRate_dX(unsigned int which, vector 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; diff --git a/bng2/Network3/src/util/mathutils/linsolve.c b/bng2/Network3/src/util/mathutils/linsolve.c index 439433b3..8f870a87 100644 --- a/bng2/Network3/src/util/mathutils/linsolve.c +++ b/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); +} diff --git a/bng2/Network3/src/util/mathutils/matrix_mult.c b/bng2/Network3/src/util/mathutils/matrix_mult.c index 31cd3650..bdd0ad71 100644 --- a/bng2/Network3/src/util/mathutils/matrix_mult.c +++ b/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); +} + diff --git a/bng2/Network3/src/util/mathutils/matrix_vector_mult.c b/bng2/Network3/src/util/mathutils/matrix_vector_mult.c index 74bccb39..e137fedd 100644 --- a/bng2/Network3/src/util/mathutils/matrix_vector_mult.c +++ b/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); +} + diff --git a/bng2/Network3/src/util/mathutils/norm.c b/bng2/Network3/src/util/mathutils/norm.c index 42a39e53..d0942592 100644 --- a/bng2/Network3/src/util/mathutils/norm.c +++ b/bng2/Network3/src/util/mathutils/norm.c @@ -1,79 +1,79 @@ -#include "mathutils.h" - - -/*==========================================================================*/ -/* returns sqrt of sum of squares of the elts of a */ -double NORM( double *a, int dim){ -#ifdef APPEND_UNDERSCORE -#define DNRM2 dnrm2_ -#else -#define DNRM2 dnrm2 -#endif - int n; - int inca=1; - double norm; - extern double DNRM2(int *n, double *x, int *incx); - - n=dim; - return(DNRM2(&n, a, &inca)); -} - -double DNRM2(int *n, double *x, int *incx) -{ - long int ix, nn, iincx; - double norm, scale, absxi, ssq, temp; - -/* DNRM2 returns the euclidean norm of a vector via the function - name, so that - - DNRM2 := sqrt( x'*x ) - - -- This version written on 25-October-1982. - Modified on 14-October-1993 to inline the call to SLASSQ. - Sven Hammarling, Nag Ltd. */ - - /* Dereference inputs */ - nn = *n; - iincx = *incx; - - if( nn > 0 && iincx > 0 ) - { - if (nn == 1) - { - norm = fabs(x[0]); - } - else - { - scale = 0.0; - ssq = 1.0; - - /* The following loop is equivalent to this call to the LAPACK - auxiliary routine: CALL SLASSQ( N, X, INCX, SCALE, SSQ ) */ - - for (ix=(nn-1)*iincx; ix>=0; ix-=iincx) - { - if (x[ix] != 0.0) - { - absxi = fabs(x[ix]); - if (scale < absxi) - { - temp = scale / absxi; - ssq = ssq * (temp * temp) + 1.0; - scale = absxi; - } - else - { - temp = absxi / scale; - ssq += temp * temp; - } - } - } - norm = scale * sqrt(ssq); - } - } - else - norm = 0.0; - - return norm; - -} /* dnrm2_ */ +#include "mathutils.h" + + +/*==========================================================================*/ +/* returns sqrt of sum of squares of the elts of a */ +double NORM( double *a, int dim){ +#ifdef APPEND_UNDERSCORE +#define DNRM2 dnrm2_ +#else +#define DNRM2 dnrm2 +#endif + int n; + int inca=1; + double norm; + extern double DNRM2(int *n, double *x, int *incx); + + n=dim; + return(DNRM2(&n, a, &inca)); +} + +double DNRM2(int *n, double *x, int *incx) +{ + long int ix, nn, iincx; + double norm, scale, absxi, ssq, temp; + +/* DNRM2 returns the euclidean norm of a vector via the function + name, so that + + DNRM2 := sqrt( x'*x ) + + -- This version written on 25-October-1982. + Modified on 14-October-1993 to inline the call to SLASSQ. + Sven Hammarling, Nag Ltd. */ + + /* Dereference inputs */ + nn = *n; + iincx = *incx; + + if( nn > 0 && iincx > 0 ) + { + if (nn == 1) + { + norm = fabs(x[0]); + } + else + { + scale = 0.0; + ssq = 1.0; + + /* The following loop is equivalent to this call to the LAPACK + auxiliary routine: CALL SLASSQ( N, X, INCX, SCALE, SSQ ) */ + + for (ix=(nn-1)*iincx; ix>=0; ix-=iincx) + { + if (x[ix] != 0.0) + { + absxi = fabs(x[ix]); + if (scale < absxi) + { + temp = scale / absxi; + ssq = ssq * (temp * temp) + 1.0; + scale = absxi; + } + else + { + temp = absxi / scale; + ssq += temp * temp; + } + } + } + norm = scale * sqrt(ssq); + } + } + else + norm = 0.0; + + return norm; + +} /* dnrm2_ */ diff --git a/bng2/Network3/src/util/mathutils/normsq.c b/bng2/Network3/src/util/mathutils/normsq.c index e4aefc8e..18d7ecad 100644 --- a/bng2/Network3/src/util/mathutils/normsq.c +++ b/bng2/Network3/src/util/mathutils/normsq.c @@ -1,32 +1,32 @@ -#include "mathutils.h" - - -/*==========================================================================*/ -/* returns sum of squares of the elts of a */ - -#define DNRSQ FNAME(dnrsq) - -double NORMSQ( double *a, int dim) -{ - int n; - int inca=1; - double norm; - extern double DNRSQ(int *n, double *a, int *inca); - - n=dim; - - return(DNRSQ(&n, a, &inca)); -} - -#ifdef NOBLAS_E -double DNRSQ( int *n, double *a, int *inca) -{ - register int inc=*inca, len=*n; - register double *x, *xlast, sum=0.0; - - xlast= a + inc*len; - for(x=a; x olddim){ /* dimension has increased */ - olddim = dim; - FREE_VECTOR(work); - lwork = dim*64; /* a rather large blocksize */ - work = ALLOC_VECTOR(lwork); - } - - Dim= dim; - Nrhs= nrhs; - ipiv = (int *) malloc (sizeof(int)*Dim); - DSYSV("l", &Dim, &Nrhs, A[0], &Dim, ipiv, B[0], &Dim, work, &lwork, - &info); - - free(ipiv); - return(info); -} +#include "mathutils.h" + + +/*==========================================================================*/ +/* Routine to Solve set of linear equations a.x=b where a is symmetric. */ +/* a is overwritten by an appropriate factorization of a and b is */ +/* overwritten by x. Note that the rhs b are stored in the ROWS of b */ +/* rather than the columns because of the way Fortran accesses arrays. */ + +int SYM_LINSOLVE(double **A, double **B, int nrhs, int dim){ +#ifdef APPEND_UNDERSCORE +#define DSYSV dsysv_ +#else +#define DSYSV dsysv +#endif + int *ipiv; + int Nrhs, Dim; + int info; + + register int i; + double det; + double a,b,c,*x,*y,cx_ay; + double detinv, ainv; + + static double *work; + static int lwork; + static int olddim; + static int initflag=1; +#ifdef __cplusplus +extern "C" { +#endif + /* Fortran dsysv: SUBROUTINE DSYSV( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, LWORK, INFO ) */ + void DSYSV(char *uplo, int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *ldb, double *work, int *lwork, int *info); +#ifdef __cplusplus +} +#endif + + if(dim==2){ + a= A[0][0]; c=A[0][1]; b = A[1][1]; + det= a*b - c*c; + detinv=1.0/det; + ainv= 1.0/a; + for (i=0; i olddim){ /* dimension has increased */ + olddim = dim; + FREE_VECTOR(work); + lwork = dim*64; /* a rather large blocksize */ + work = ALLOC_VECTOR(lwork); + } + + Dim= dim; + Nrhs= nrhs; + ipiv = (int *) malloc (sizeof(int)*Dim); + DSYSV("l", &Dim, &Nrhs, A[0], &Dim, ipiv, B[0], &Dim, work, &lwork, + &info); + + free(ipiv); + return(info); +} diff --git a/bng2/Network3/src/util/mathutils/sym_matrix_mult.c b/bng2/Network3/src/util/mathutils/sym_matrix_mult.c index cc314432..24012aad 100644 --- a/bng2/Network3/src/util/mathutils/sym_matrix_mult.c +++ b/bng2/Network3/src/util/mathutils/sym_matrix_mult.c @@ -1,25 +1,25 @@ -#include "mathutils.h" - - -/*==========================================================================*/ -void SYM_MATRIX_MULT(double **a, double **b, double **c, int dim){ -#ifdef APPEND_UNDERSCORE -#define DSYMM dsymm_ -#else -#define DSYMM dsymm -#endif - /* uses upper triangular region of matrix only */ - /* uses BLAS3 routine dsymm */ - int n=dim; - double alpha=1.0; - double beta=0.0; - extern void DSYMM(char*, char*, int*, int*, double*, double*, int*, double*, int*, double*, double*, int*); - - /* transpose necessary because of opposite way Fortran references */ - /* arrays */ - MATRIX_TRANSPOSE(b,dim); - DSYMM("l", "l", &n, &n, &alpha, a[0], &n, b[0], &n, &beta, c[0], &n); - MATRIX_TRANSPOSE(b,dim); - MATRIX_TRANSPOSE(c,dim); -} - +#include "mathutils.h" + + +/*==========================================================================*/ +void SYM_MATRIX_MULT(double **a, double **b, double **c, int dim){ +#ifdef APPEND_UNDERSCORE +#define DSYMM dsymm_ +#else +#define DSYMM dsymm +#endif + /* uses upper triangular region of matrix only */ + /* uses BLAS3 routine dsymm */ + int n=dim; + double alpha=1.0; + double beta=0.0; + extern void DSYMM(char*, char*, int*, int*, double*, double*, int*, double*, int*, double*, double*, int*); + + /* transpose necessary because of opposite way Fortran references */ + /* arrays */ + MATRIX_TRANSPOSE(b,dim); + DSYMM("l", "l", &n, &n, &alpha, a[0], &n, b[0], &n, &beta, c[0], &n); + MATRIX_TRANSPOSE(b,dim); + MATRIX_TRANSPOSE(c,dim); +} + diff --git a/bng2/Network3/src/util/mathutils/sym_matrix_vector_mult.c b/bng2/Network3/src/util/mathutils/sym_matrix_vector_mult.c index cafa3350..d7eaa56c 100644 --- a/bng2/Network3/src/util/mathutils/sym_matrix_vector_mult.c +++ b/bng2/Network3/src/util/mathutils/sym_matrix_vector_mult.c @@ -1,21 +1,21 @@ -#include "mathutils.h" - - -/*==========================================================================*/ -void SYM_MATRIX_VECTOR_MULT(double **a, double *x, double *y, int dim){ -#ifdef APPEND_UNDERSCORE -#define DSYMV dsymv_ -#else -#define DSYMV dsymv -#endif - /* uses BLAS3 routine dsymv */ - int n=dim; - int incx=1; - int incy=1; - double alpha=1.0; - double beta=0.0; - extern void DSYMV(char*, int*, double*, double*, int*, double*, int*, double*, double*, int*); - - DSYMV("l", &n, &alpha, a[0], &n, x, &incx, &beta, y, &incy); -} - +#include "mathutils.h" + + +/*==========================================================================*/ +void SYM_MATRIX_VECTOR_MULT(double **a, double *x, double *y, int dim){ +#ifdef APPEND_UNDERSCORE +#define DSYMV dsymv_ +#else +#define DSYMV dsymv +#endif + /* uses BLAS3 routine dsymv */ + int n=dim; + int incx=1; + int incy=1; + double alpha=1.0; + double beta=0.0; + extern void DSYMV(char*, int*, double*, double*, int*, double*, int*, double*, double*, int*); + + DSYMV("l", &n, &alpha, a[0], &n, x, &incx, &beta, y, &incy); +} + diff --git a/bng2/Network3/src/util/mathutils/sympos_linsolve.c b/bng2/Network3/src/util/mathutils/sympos_linsolve.c index 4e580d22..806fe63a 100644 --- a/bng2/Network3/src/util/mathutils/sympos_linsolve.c +++ b/bng2/Network3/src/util/mathutils/sympos_linsolve.c @@ -1,33 +1,33 @@ -#include "mathutils.h" - - -/*==========================================================================*/ -/* Routine to Solve set of linear equations a.x=b where a is symmetric */ -/* and positive definite. */ -/* a is overwritten by an appropriate factorization of a and b is */ -/* overwritten by x. Note that the rhs b are stored in the ROWS of b */ -/* rather than the columns because of the way Fortran accesses arrays. */ - -int SYMPOS_LINSOLVE(double **a, double **b, int nrhs, int dim){ -#ifdef APPEND_UNDERSCORE -#define DPOSV dposv_ -#else -#define DPOSV dposv -#endif - int Nrhs, Dim; - int info; -#ifdef __cplusplus -extern "C" { -#endif - /* Fortran dposv: SUBROUTINE DPOSV( UPLO, N, NRHS, A, LDA, B, LDB, INFO ) */ - void DPOSV(char *uplo, int *n, int *nrhs, double *a, int *lda, double *b, int *ldb, int *info); -#ifdef __cplusplus -} -#endif - - Dim= dim; - Nrhs= nrhs; - DPOSV("l", &Dim, &Nrhs, a[0], &Dim, b[0], &Dim, &info); - - return(info); -} +#include "mathutils.h" + + +/*==========================================================================*/ +/* Routine to Solve set of linear equations a.x=b where a is symmetric */ +/* and positive definite. */ +/* a is overwritten by an appropriate factorization of a and b is */ +/* overwritten by x. Note that the rhs b are stored in the ROWS of b */ +/* rather than the columns because of the way Fortran accesses arrays. */ + +int SYMPOS_LINSOLVE(double **a, double **b, int nrhs, int dim){ +#ifdef APPEND_UNDERSCORE +#define DPOSV dposv_ +#else +#define DPOSV dposv +#endif + int Nrhs, Dim; + int info; +#ifdef __cplusplus +extern "C" { +#endif + /* Fortran dposv: SUBROUTINE DPOSV( UPLO, N, NRHS, A, LDA, B, LDB, INFO ) */ + void DPOSV(char *uplo, int *n, int *nrhs, double *a, int *lda, double *b, int *ldb, int *info); +#ifdef __cplusplus +} +#endif + + Dim= dim; + Nrhs= nrhs; + DPOSV("l", &Dim, &Nrhs, a[0], &Dim, b[0], &Dim, &info); + + return(info); +} diff --git a/bng2/Network3/src/util/mathutils/test_transpose.c b/bng2/Network3/src/util/mathutils/test_transpose.c index 1b59a46a..06f2db19 100644 --- a/bng2/Network3/src/util/mathutils/test_transpose.c +++ b/bng2/Network3/src/util/mathutils/test_transpose.c @@ -1,11 +1,11 @@ -#include "mathutils.h" - -int main(void) -{ - double dat[]={1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0}; - dcomplex **m; - - m = ZALLOC_MATRIX(2,2); - COPY_VECTOR(dat, (double *) m[0], 8); - ZMATRIX_TRANSPOSE(m, 2); -} +#include "mathutils.h" + +int main(void) +{ + double dat[]={1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0}; + dcomplex **m; + + m = ZALLOC_MATRIX(2,2); + COPY_VECTOR(dat, (double *) m[0], 8); + ZMATRIX_TRANSPOSE(m, 2); +} diff --git a/bng2/Network3/src/util/mathutils/testdummy.c b/bng2/Network3/src/util/mathutils/testdummy.c index 9718ed64..812d5854 100644 --- a/bng2/Network3/src/util/mathutils/testdummy.c +++ b/bng2/Network3/src/util/mathutils/testdummy.c @@ -1,18 +1,18 @@ -#include "mathutils.h" - -int main(void){ - int i,n; - int incx=1; - double *x; - - /* read n */ - scanf("%d", &n); - /* allocate space for x */ - x = (double *) malloc(n*sizeof(double)); - /* read x */ - if (n<0){ fprintf(stderr,"n must be greater than 0.\n"); exit(1);} - for (i=0; i\n", argv[0]); - exit(EXIT_FAILURE); - } - m = atoi(argv[1]); - n = atoi(argv[2]); - n_run = atoi(argv[3]); - if (argc == 5){ - if ((stream = fopen(argv[4], "r")) == NULL){ - fprintf(stderr, "Can't open file %s.\n", argv[4]); - exit(EXIT_FAILURE); - } - } else { - stream = stdin; - } - - A = ALLOC_MATRIX(m,n); -/* for(i=0; i\n", argv[0]); + exit(EXIT_FAILURE); + } + m = atoi(argv[1]); + n = atoi(argv[2]); + n_run = atoi(argv[3]); + if (argc == 5){ + if ((stream = fopen(argv[4], "r")) == NULL){ + fprintf(stderr, "Can't open file %s.\n", argv[4]); + exit(EXIT_FAILURE); + } + } else { + stream = stdin; + } + + A = ALLOC_MATRIX(m,n); +/* for(i=0; i\n", argv[0]); - exit(EXIT_FAILURE); - } - m = atoi(argv[1]); - n = atoi(argv[2]); - n_run = atoi(argv[3]); - if (argc == 5){ - if ((stream = fopen(argv[4], "r")) == NULL){ - fprintf(stderr, "Can't open file %s.\n", argv[4]); - exit(EXIT_FAILURE); - } - } else { - stream = stdin; - } - - A = ALLOC_MATRIX(m,n); - for(i=0; i 1.0e-12){ - for(j=0, norm=0.0; j\n", argv[0]); + exit(EXIT_FAILURE); + } + m = atoi(argv[1]); + n = atoi(argv[2]); + n_run = atoi(argv[3]); + if (argc == 5){ + if ((stream = fopen(argv[4], "r")) == NULL){ + fprintf(stderr, "Can't open file %s.\n", argv[4]); + exit(EXIT_FAILURE); + } + } else { + stream = stdin; + } + + A = ALLOC_MATRIX(m,n); + for(i=0; i 1.0e-12){ + for(j=0, norm=0.0; jj) zb[i][j] = Cconj(z); - else zb[i][j]=COMPLEX_ZERO; - } - HERM_MATRIX_MULT(za,zb,zc,dim); - for (i=0, errcount=0; ij) - { - if ((zc[i][j].r != z.r) || (zc[i][j].i != -z.i)) ++errcount; - } - else - { - if ((zc[i][j].r != 0.0) || (zc[i][j].i != 0.0)) ++errcount; - } - } - if (!errcount) printf("OK.\n"); - else printf("Failed.\n"); - ZFREE_MATRIX(za); - ZFREE_MATRIX(zb); - ZFREE_MATRIX(zc); - - printf("Testing MATRIX_VECTOR_MULT..."); - x = ALLOC_VECTOR(dim); - y = ALLOC_VECTOR(dim); - a= ALLOC_MATRIX(dim,dim); - INIT_VECTOR(x,1.0,dim); - INIT_VECTOR(y,1.0,dim); - INIT_VECTOR(a[0],1.0,dim*dim); - MATRIX_VECTOR_MULT(a,x,y,dim); - for (i=0, errcount=0; ij) za[i][j].i=-1.0; - } - z.r=1.0; z.i=0.0; - ZINIT_VECTOR(zx,z,dim); - ZINIT_VECTOR(zy,z,dim); - HERM_MATRIX_VECTOR_MULT(za,zx,zy,dim); - for (i=0, errcount=0; iTOL) ++errcount; - else if (fabs(c[i][i]-1.0)>TOL) ++errcount; - } - } - if (errcount ==0) printf("OK.\n"); - else printf("Failed.\n"); - FREE_MATRIX(a); - FREE_MATRIX(b); - FREE_MATRIX(c); - - printf("Testing SYMPOS_LINSOLVE..."); - a= ALLOC_MATRIX(2,2); - b=ALLOC_MATRIX(2,2); - c=ALLOC_MATRIX(2,2); - INIT_VECTOR(a[0],0.0,4); - a[0][0]=1; a[0][1] = .5; a[1][1]=1; - INIT_VECTOR(b[0],0.0,4); - b[0][0] = b[1][1] = 1; - COPY_VECTOR(a[0],c[0],4); - SYMPOS_LINSOLVE(c,b,2,2); - SYM_MATRIX_MULT(a,b,c,2); - for(i=0, errcount=0; i<2; ++i){ - for(j=0; j<2; ++j){ - if(i==j && fabs(c[i][i]-1.0)>TOL) ++errcount; - else if (fabs(c[i][i]-1.0)>TOL) ++errcount; - } - } - if (errcount ==0) printf("OK.\n"); - else printf("Failed.\n"); - FREE_MATRIX(a); - FREE_MATRIX(b); - FREE_MATRIX(c); -} - - - +/* program to test mathutilities */ +#include "mathutils.h" + +#define TOL 1.0e-14 + +int main (int argc, char *argv[]) +{ + int i,j,k; + int dim; + int seed; + int errcount; + double *x, *y; + double r; + dcomplex z; + dcomplex z1; + dcomplex *zx, *zy; + double **a, **b, **c; + dcomplex **za, **zb, **zc; + + if (argc < 2) + { + fprintf(stderr,"Usage: %s dim\n", argv[0]); + exit(EXIT_FAILURE); + } + + dim = atoi(argv[1]); + + /* turn off buffering of output */ + setbuf(stdout, NULL); + + /* allocation and initialization of vector */ + printf("Testing ALLOC_VECTOR..."); + x = ALLOC_VECTOR( dim); + printf("OK.\n"); + + printf("Testing INIT_VECTOR..."); + INIT_VECTOR( x, 1.0, dim); + for (i=0, errcount=0; ij) zb[i][j] = Cconj(z); + else zb[i][j]=COMPLEX_ZERO; + } + HERM_MATRIX_MULT(za,zb,zc,dim); + for (i=0, errcount=0; ij) + { + if ((zc[i][j].r != z.r) || (zc[i][j].i != -z.i)) ++errcount; + } + else + { + if ((zc[i][j].r != 0.0) || (zc[i][j].i != 0.0)) ++errcount; + } + } + if (!errcount) printf("OK.\n"); + else printf("Failed.\n"); + ZFREE_MATRIX(za); + ZFREE_MATRIX(zb); + ZFREE_MATRIX(zc); + + printf("Testing MATRIX_VECTOR_MULT..."); + x = ALLOC_VECTOR(dim); + y = ALLOC_VECTOR(dim); + a= ALLOC_MATRIX(dim,dim); + INIT_VECTOR(x,1.0,dim); + INIT_VECTOR(y,1.0,dim); + INIT_VECTOR(a[0],1.0,dim*dim); + MATRIX_VECTOR_MULT(a,x,y,dim); + for (i=0, errcount=0; ij) za[i][j].i=-1.0; + } + z.r=1.0; z.i=0.0; + ZINIT_VECTOR(zx,z,dim); + ZINIT_VECTOR(zy,z,dim); + HERM_MATRIX_VECTOR_MULT(za,zx,zy,dim); + for (i=0, errcount=0; iTOL) ++errcount; + else if (fabs(c[i][i]-1.0)>TOL) ++errcount; + } + } + if (errcount ==0) printf("OK.\n"); + else printf("Failed.\n"); + FREE_MATRIX(a); + FREE_MATRIX(b); + FREE_MATRIX(c); + + printf("Testing SYMPOS_LINSOLVE..."); + a= ALLOC_MATRIX(2,2); + b=ALLOC_MATRIX(2,2); + c=ALLOC_MATRIX(2,2); + INIT_VECTOR(a[0],0.0,4); + a[0][0]=1; a[0][1] = .5; a[1][1]=1; + INIT_VECTOR(b[0],0.0,4); + b[0][0] = b[1][1] = 1; + COPY_VECTOR(a[0],c[0],4); + SYMPOS_LINSOLVE(c,b,2,2); + SYM_MATRIX_MULT(a,b,c,2); + for(i=0, errcount=0; i<2; ++i){ + for(j=0; j<2; ++j){ + if(i==j && fabs(c[i][i]-1.0)>TOL) ++errcount; + else if (fabs(c[i][i]-1.0)>TOL) ++errcount; + } + } + if (errcount ==0) printf("OK.\n"); + else printf("Failed.\n"); + FREE_MATRIX(a); + FREE_MATRIX(b); + FREE_MATRIX(c); +} + + + diff --git a/bng2/Network3/src/util/mathutils/zmatrix_add.c b/bng2/Network3/src/util/mathutils/zmatrix_add.c index ea9b7a8b..b8065ac0 100644 --- a/bng2/Network3/src/util/mathutils/zmatrix_add.c +++ b/bng2/Network3/src/util/mathutils/zmatrix_add.c @@ -1,48 +1,48 @@ -#include "mathutils.h" - -/*==========================================================================*/ -void ZMATRIX_ADD(dcomplex **a, dcomplex **b, dcomplex **c, int dim){ - register int i,j; - double *aptr, *bptr, *cptr; - - for(i=0; i -#include -int main(void){ - int i,j,n; - int incx=1; - dcomplex **a, **b, *x; - - /* read n */ - scanf("%d", &n); - /* allocate space for x */ - a = ZALLOC_MATRIX(n,n); - b = ZALLOC_MATRIX(n,n); - /* read a and b */ - if (n<0){ fprintf(stderr,"n must be greater than 0.\n"); exit(1);} - x= a[0]; - for (i=0; i +#include +int main(void){ + int i,j,n; + int incx=1; + dcomplex **a, **b, *x; + + /* read n */ + scanf("%d", &n); + /* allocate space for x */ + a = ZALLOC_MATRIX(n,n); + b = ZALLOC_MATRIX(n,n); + /* read a and b */ + if (n<0){ fprintf(stderr,"n must be greater than 0.\n"); exit(1);} + x= a[0]; + for (i=0; i -#include -int main(void){ - int i,j,n; - int incx=1; - dcomplex **a, **b, *x; - - /* read n */ - scanf("%d", &n); - /* allocate space for x */ - a = ZALLOC_MATRIX(n,n); - b = ZALLOC_MATRIX(n,n); - /* read a and b */ - if (n<0){ fprintf(stderr,"n must be greater than 0.\n"); exit(1);} - x= a[0]; - for (i=0; i +#include +int main(void){ + int i,j,n; + int incx=1; + dcomplex **a, **b, *x; + + /* read n */ + scanf("%d", &n); + /* allocate space for x */ + a = ZALLOC_MATRIX(n,n); + b = ZALLOC_MATRIX(n,n); + /* read a and b */ + if (n<0){ fprintf(stderr,"n must be greater than 0.\n"); exit(1);} + x= a[0]; + for (i=0; i olddim){ /* dimension has increased */ - olddim = dim; - ZFREE_VECTOR(work); - lwork = dim*64; /* a rather large blocksize */ - work = ZALLOC_VECTOR(lwork); - } - - Dim= dim; - Nrhs= nrhs; - ipiv = (int *) malloc (sizeof(int)*Dim); - ZSYSV("l", &Dim, &Nrhs, a[0], &Dim, ipiv, b[0], &Dim, work, &lwork, - &info); - - free(ipiv); - return(info); -} +#include "mathutils.h" + + +/*==========================================================================*/ +/* Routine to Solve COMPLEX set of linear equations a.x=b where a is +symmetric. a is overwritten by an appropriate factorization of a +and b is overwritten by x. Note that the rhs b are stored in the +ROWS of b rather than the columns because of the way Fortran +accesses arrays. */ + +int ZSYM_LINSOLVE(dcomplex **a, dcomplex **b, int nrhs, int dim){ +#ifdef APPEND_UNDERSCORE +#define ZSYSV zsysv_ +#else +#define ZSYSV zsysv +#endif + int *ipiv; + int Nrhs, Dim; + int info; + static dcomplex *work; + static int lwork; + static int olddim; + static int initflag=1; +#ifdef __cplusplus +extern "C" { +#endif + /* Fortran zsysv: SUBROUTINE ZSYSV( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, LWORK, INFO ) */ + void ZSYSV(char *uplo, int *n, int *nrhs, dcomplex *a, int *lda, int *ipiv, dcomplex *b, int *ldb, dcomplex *work, int *lwork, int *info); +#ifdef __cplusplus +} +#endif + + if (initflag){ + initflag = 0; + olddim = dim; + lwork = dim*64; /* a rather large blocksize */ + work = ZALLOC_VECTOR(lwork); + } + + if ( dim > olddim){ /* dimension has increased */ + olddim = dim; + ZFREE_VECTOR(work); + lwork = dim*64; /* a rather large blocksize */ + work = ZALLOC_VECTOR(lwork); + } + + Dim= dim; + Nrhs= nrhs; + ipiv = (int *) malloc (sizeof(int)*Dim); + ZSYSV("l", &Dim, &Nrhs, a[0], &Dim, ipiv, b[0], &Dim, work, &lwork, + &info); + + free(ipiv); + return(info); +} diff --git a/bng2/Perl2/Aux2/count_blocks.pl b/bng2/Perl2/Aux2/count_blocks.pl index 09f36a00..2f3d79c2 100755 --- a/bng2/Perl2/Aux2/count_blocks.pl +++ b/bng2/Perl2/Aux2/count_blocks.pl @@ -1,7 +1,7 @@ #!/usr/bin/perl for $file (@ARGV){ - open(FILE,$file) || die "Couldn't open $file: $!\n"; + open(FILE, '<', $file) || die "Couldn't open $file: $!\n"; print "$file:\n"; while(){ if (/^s*begin\s+(\S.*)$/){ diff --git a/bng2/Perl2/Aux2/count_lines.pl b/bng2/Perl2/Aux2/count_lines.pl index c1dc50dd..5c2d4684 100755 --- a/bng2/Perl2/Aux2/count_lines.pl +++ b/bng2/Perl2/Aux2/count_lines.pl @@ -39,7 +39,7 @@ my $n_lines_tot=0; for my $file (@files){ my $n_lines=0; - open(IN,$file) || die "Couldn't open $file:$?\n"; + open(IN, '<', $file) || die "Couldn't open $file:$?\n"; while(){ if (/;\s*$/){ ++$n_lines; diff --git a/bng2/Perl2/Aux2/diff_cdat.pl b/bng2/Perl2/Aux2/diff_cdat.pl index d755933a..3232af01 100755 --- a/bng2/Perl2/Aux2/diff_cdat.pl +++ b/bng2/Perl2/Aux2/diff_cdat.pl @@ -44,7 +44,7 @@ #print "Comparing $rfilename $rfilename2\n"; -open (RFILE, $rfilename) or die "Can't open $rfilename: $!\n"; +open (RFILE, '<', $rfilename) or die "Can't open $rfilename: $!\n"; $i=0; @data=(); @times=(); @@ -60,7 +60,7 @@ } close(RFILE); -open (RFILE, $rfilename2) or die "Can't open $rfilename2: $!\n"; +open (RFILE, '<', $rfilename2) or die "Can't open $rfilename2: $!\n"; $i=0; @data2=(); @times2=(); diff --git a/bng2/Perl2/Aux2/runBNG.pl b/bng2/Perl2/Aux2/runBNG.pl index bc984b70..b0a2a0b3 100755 --- a/bng2/Perl2/Aux2/runBNG.pl +++ b/bng2/Perl2/Aux2/runBNG.pl @@ -46,7 +46,7 @@ open my $oldout, ">&STDOUT" or die "Can't dup STDOUT: $!"; # Redirect STDOUT to logfile -open(STDOUT,">${prefix}_runBNG.log"); +open( STDOUT, '>', "${prefix}_runBNG.log" ); # turn off output buffering on STDOUT (select(*STDOUT), $|=1)[0]; @@ -111,7 +111,7 @@ my @output=(); if (-r "$gdatfile"){ print "Updating observable concentrations from $gdatfile\n"; - open(GDAT, "$gdatfile"); + open( GDAT, '<', $gdatfile ); my $last=""; while(){ $last=$_; diff --git a/bng2/Perl2/BNGAction.pm b/bng2/Perl2/BNGAction.pm index 668be77a..4a131eff 100644 --- a/bng2/Perl2/BNGAction.pm +++ b/bng2/Perl2/BNGAction.pm @@ -1624,7 +1624,7 @@ sub bifurcate # extract forward scan data my @forward; - open FILE, $scanfiles[0] or die "Couldn't open file: $!"; + open(FILE, "<", $scanfiles[0]) or die "Couldn't open file: $!"; my $line = ; # first line chomp $line; $line =~ s/^\s*#\s+//; # remove leading # and whitespace @@ -1643,7 +1643,7 @@ sub bifurcate # extract backward scan data my @backward; - open FILE, $scanfiles[1] or die "Couldn't open file: $!"; + open(FILE, "<", $scanfiles[1]) or die "Couldn't open file: $!"; $line = ; # first line $i = 0; while ($line = ){ diff --git a/bng2/Perl2/BNGModel.pm b/bng2/Perl2/BNGModel.pm index 7792404d..a73519a6 100644 --- a/bng2/Perl2/BNGModel.pm +++ b/bng2/Perl2/BNGModel.pm @@ -2780,7 +2780,11 @@ sub generate_network my $vmem = 0; if ($^O eq 'MSWin32') { # Windows: use tasklist - my $output = `tasklist /FI "PID eq $$" /NH /FO CSV`; + my $output = ''; + if (open(my $ph, '-|', 'tasklist', '/FI', "PID eq $$", '/NH', '/FO', 'CSV')) { + $output = do { local $/; <$ph> }; + close($ph); + } if ($output =~ /"([^"]+)"\s*$/) { my $mem_str = $1; $mem_str =~ s/[^\d]//g; @@ -2790,7 +2794,11 @@ sub generate_network } else { # Linux/Unix: use ps - my $ps_out = `ps -o rss,vsz -p $$`; + my $ps_out = ''; + if (open(my $ph, '-|', 'ps', '-o', 'rss,vsz', '-p', $$)) { + $ps_out = do { local $/; <$ph> }; + close($ph); + } my @lines = split /\n/, $ps_out; if (@lines > 1) { my @cols = split ' ', $lines[1]; diff --git a/bng2/Perl2/BNGOutput.pm b/bng2/Perl2/BNGOutput.pm index 22b8b018..5fef999b 100644 --- a/bng2/Perl2/BNGOutput.pm +++ b/bng2/Perl2/BNGOutput.pm @@ -1300,7 +1300,7 @@ sub writeSSC $prefix .= "_${suffix}"; } my $file = "${prefix}.rxn"; - open( SSCfile, ">$file" ) || die "Couldn't open $file: $!\n"; + open( SSCfile, '>', $file ) || die "Couldn't open $file: $!\n"; my $version = BNGversion(); print SSCfile "--# SSC-file for model $model_name created by BioNetGen $version\n"; @@ -1389,7 +1389,7 @@ sub writeSSCcfg my $file = "${prefix}.cfg"; my $version = BNGversion(); - open( SSCcfgfile, ">$file" ) || die "Couldn't open $file: $!\n"; + open( SSCcfgfile, '>', $file ) || die "Couldn't open $file: $!\n"; print STDOUT "\n Writting SSC cfg file \n"; print SSCcfgfile "# SSC cfg file for model $model_name created by BioNetGen $version\n"; print SSCcfgfile $model->ParamList->writeSSCcfg(); @@ -1642,7 +1642,7 @@ sub writeMfile # open Mexfile and begin printing... - open( Mscript, ">$mscript_path" ) || die "Couldn't open $mscript_path: $!\n"; + open( Mscript, '>', $mscript_path ) || die "Couldn't open $mscript_path: $!\n"; print Mscript <<"EOF"; function [err, timepoints, species_out, observables_out] = ${mscript_filebase}( timepoints, species_init, parameters, suppress_plot ) %${mscript_filebase_caps} Integrate reaction network and plot observables. @@ -2157,7 +2157,7 @@ sub writeMexfile # open Mexfile and begin printing... - open( Mexfile, ">$mex_path" ) or die "Couldn't open $mex_path: $!\n"; + open( Mexfile, '>', $mex_path ) or die "Couldn't open $mex_path: $!\n"; print Mexfile <<"EOF"; /* ** ${mex_filename} @@ -2556,7 +2556,7 @@ EOF # open Mexfile and begin printing... - open( Mscript, ">$mscript_path" ) or die "Couldn't open $mscript_path: $!\n"; + open( Mscript, '>', $mscript_path ) or die "Couldn't open $mscript_path: $!\n"; print Mscript <<"EOF"; function [err, timepoints, species_out, observables_out ] = ${mscript_filebase}( timepoints, species_init, parameters, suppress_plot ) %${mscript_filebase_caps} Integrate reaction network and plot observables. @@ -2905,7 +2905,7 @@ sub writeCPPfile if ($err) { return $err }; # open Mexfile and begin printing... - open( Cppfile, ">$cpp_path" ) or die "Couldn't open $cpp_path: $!\n"; + open( Cppfile, '>', $cpp_path ) or die "Couldn't open $cpp_path: $!\n"; print Cppfile <<"EOF"; /* ** ${cpp_filename} @@ -3499,7 +3499,7 @@ sub writeCPYfile if ($err) { return $err }; # open C and begin printing... - open( Cpyfile, ">$cpy_path" ) or die "Couldn't open $cpy_path: $!\n"; + open( Cpyfile, '>', $cpy_path ) or die "Couldn't open $cpy_path: $!\n"; print Cpyfile <<"EOF"; /* ** ${cpy_filename} @@ -3530,7 +3530,29 @@ sub writeCPYfile ** ** Usage in Python : ** -** TODO +** import ctypes +** +** class RESULT(ctypes.Structure): +** _fields_ = [ +** ("status", ctypes.c_int), +** ("n_observables", ctypes.c_int), +** ("n_species", ctypes.c_int), +** ("n_tpts", ctypes.c_int), +** ("obs_name_len", ctypes.c_int), +** ("spcs_name_len", ctypes.c_int), +** ("observables", ctypes.POINTER(ctypes.c_double)), +** ("species", ctypes.POINTER(ctypes.c_double)), +** ("obs_names", ctypes.c_char_p), +** ("spcs_names", ctypes.c_char_p) +** ] +** +** lib = ctypes.CDLL('./$model_name.so') +** lib.simulate.restype = ctypes.POINTER(RESULT) +** +** # Define inputs +** # ... define num_tpts, timepts, num_species_init, species_init, num_parameters, parameters ... +** +** res = lib.simulate(num_tpts, timepts, num_species_init, species_init, num_parameters, parameters) */ /* Library headers */ @@ -3945,7 +3967,7 @@ sub writeMfile_QueryNames my $q_mscript = 'QueryNames.m'; - open(Q_Mscript,">$q_mscript"); + open( Q_Mscript, '>', $q_mscript ); print Q_Mscript <<"EOF"; function [ param_labels, param_defaults, obs_labels, species_labels] = QueryNames( inputlist ) % % Loads all the parameter labels, parameter defaults, observable labels and species labels in the model @@ -4002,7 +4024,7 @@ sub writeMfile_ParametersObservables #Writing parameter list script - open( Par_Mscript, ">$par_mscript" ) || die "Couldn't open $par_mscript: $!\n"; + open( Par_Mscript, '>', $par_mscript ) || die "Couldn't open $par_mscript: $!\n"; print Par_Mscript <<"EOF"; function [outputlist,defaultvals ] = ParameterList( inputlist ) % Used to manipulate and access parameter names @@ -4053,7 +4075,7 @@ EOF print "Wrote M-file script $par_mscript.\n"; #Writing observable list script - open( Obs_Mscript, ">$obs_mscript" ) || die "Couldn't open $obs_mscript: $!\n"; + open( Obs_Mscript, '>', $obs_mscript ) || die "Couldn't open $obs_mscript: $!\n"; print Obs_Mscript <<"EOF"; function [outputlist ] = ObservableList( inputlist ) % Used to manipulate and access observable names @@ -4134,7 +4156,7 @@ sub writeLatex # open file my $Lfile; - open( $Lfile, ">$file" ) or die "Couldn't open $file: $!\n"; + open( $Lfile, '>', $file ) or die "Couldn't open $file: $!\n"; my $version = BNGversion(); print$Lfile "% Latex formatted differential equations for model $prefix created by BioNetGen $version\n"; @@ -4460,7 +4482,7 @@ sub writeMfile_all # open Mfile and begin printing... - open( Mscript, ">$mscript_path" ) || die "Couldn't open $mscript_path: $!\n"; + open( Mscript, '>', $mscript_path ) || die "Couldn't open $mscript_path: $!\n"; print Mscript <<"EOF"; function [err, timepoints, species_out, observables_out ] = ${mscript_filebase}( timepoints, species_init, parameters, suppress_plot ) %${mscript_filebase_caps} Integrate reaction network and plot observables. @@ -4680,7 +4702,7 @@ EOF $mscript_path = File::Spec->catpath($vol,$path,$mscript_filename); $mscript_filebase_caps = uc $mscript_filebase; - open( Mscript, ">$mscript_path" ) || die "Couldn't open $mscript_path: $!\n"; + open( Mscript, '>', $mscript_path ) || die "Couldn't open $mscript_path: $!\n"; print Mscript <<"EOF"; function [species_init] = initialize_species( params ) @@ -4698,7 +4720,7 @@ EOF $mscript_path = File::Spec->catpath($vol,$path,$mscript_filename); $mscript_filebase_caps = uc $mscript_filebase; - open( Mscript, ">$mscript_path" ) || die "Couldn't open $mscript_path: $!\n"; + open( Mscript, '>', $mscript_path ) || die "Couldn't open $mscript_path: $!\n"; print Mscript <<"EOF"; classdef ${mscript_filebase} < bngModel diff --git a/bng2/Perl2/BNGUtils.pm b/bng2/Perl2/BNGUtils.pm index c24f69c2..de944d7e 100644 --- a/bng2/Perl2/BNGUtils.pm +++ b/bng2/Perl2/BNGUtils.pm @@ -516,7 +516,7 @@ sub average_runs{ my $ng; my @y; for my $file (@_){ - open(IN, $file); + open(IN, '<', $file); my $i_t=0; while(){ next if (/^\#/); @@ -535,7 +535,7 @@ sub average_runs{ } # Write results to outfile - open(OUT,">$outfile"); + open(OUT, ">", $outfile); for my $j (0..$ng){ for my $i (0..$#t){ print OUT $t[$i]; diff --git a/bng2/Perl2/Boolean2BNGL/Boolean2BNGL.pl b/bng2/Perl2/Boolean2BNGL/Boolean2BNGL.pl index 3a340d8b..a38a18e1 100644 --- a/bng2/Perl2/Boolean2BNGL/Boolean2BNGL.pl +++ b/bng2/Perl2/Boolean2BNGL/Boolean2BNGL.pl @@ -37,7 +37,7 @@ # read entire file as a string local $/ = undef; -open FILE, $file or die "Couldn't open file: $!"; +open(FILE, "<", $file) or die "Couldn't open file: $!"; binmode FILE; my $text = ; close FILE; diff --git a/bng2/Perl2/EnergyPattern.pm b/bng2/Perl2/EnergyPattern.pm index b9941f77..9dd500ec 100644 --- a/bng2/Perl2/EnergyPattern.pm +++ b/bng2/Perl2/EnergyPattern.pm @@ -113,24 +113,32 @@ sub toString sub toMatlabString { - my $epatt = shift; - my $string = ''; - - # TODO - - return $string, ''; + my $epatt = shift; + my $plist = (@_) ? shift : ''; + my $indent = (@_) ? shift : ''; + + if ($epatt->Gf) + { + return $epatt->Gf->toMatlabString($plist, $indent), ''; + } + + return '', ''; } sub toMexString { - my $epatt = shift; - my $string = ''; - - # TODO - - return $string, ''; + my $epatt = shift; + my $plist = (@_) ? shift : ''; + my $indent = (@_) ? shift : ''; + + if ($epatt->Gf) + { + return $epatt->Gf->toMexString($plist, $indent), ''; + } + + return '', ''; } @@ -162,12 +170,16 @@ sub toXML sub toMathMLString { - my $epatt = shift; - my $string = ''; + my $epatt = shift; + my $plist = (@_) ? shift : ''; + my $indent = (@_) ? shift : ''; - # TODO + if ($epatt->Gf) + { + return $epatt->Gf->toMathMLString($plist, $indent); + } - return $string, ''; + return (''); } diff --git a/bng2/Perl2/MacroBNG2.pl b/bng2/Perl2/MacroBNG2.pl index 5a652d68..0b58ce8f 100755 --- a/bng2/Perl2/MacroBNG2.pl +++ b/bng2/Perl2/MacroBNG2.pl @@ -18,7 +18,8 @@ #Calculation our $model; -eval '$model = new '.$options->config->{ModelID}; +my $model_class = $options->config->{ModelID}; +$model = $model_class->new(); $model->ProcessModel($options->config->{bnglfile}); message("Processing complete.\n"); diff --git a/bng2/Perl2/MacroBNGModel.pm b/bng2/Perl2/MacroBNGModel.pm index 863b9f37..0bf0fde1 100644 --- a/bng2/Perl2/MacroBNGModel.pm +++ b/bng2/Perl2/MacroBNGModel.pm @@ -86,7 +86,7 @@ $err = $slist->readString($entry,$base_model->ParamList,$base_model->MoleculeTyp $file= "macr_".$param_prefix.".bngl"; rename($file,$filen); - open (WFILEbngl, ">$file") or die "Can't open $file: $!\n"; + open (WFILEbngl, ">", $file) or die "Can't open $file: $!\n"; print WFILEbngl $simul1; close (WFILEbngl); $params{file}=$file; # macr_fceri_ji3.bngl @@ -113,18 +113,18 @@ my ($param_prefix) = @_; my %bngdata; - $bnglfile= "<${param_prefix}.bngl"; # $bnglfile= "<${param_prefix}.bngl"; - $recfile= ">macr_${param_prefix}.rec"; # $recfile= ">macr_${param_prefix}.rec"; - $parfile= ">macr_${param_prefix}.par"; # $parfile= ">macr_${param_prefix}.par"; + $bnglfile= "${param_prefix}.bngl"; # $bnglfile= "${param_prefix}.bngl"; + $recfile= "macr_${param_prefix}.rec"; # $recfile= "macr_${param_prefix}.rec"; + $parfile= "macr_${param_prefix}.par"; # $parfile= "macr_${param_prefix}.par"; - open (WFILErec, $recfile) or die "Can't open $recfile: $!\n"; - open (WFILEpar, $parfile) or die "Can't open $parfile: $!\n"; + open (WFILErec, ">", $recfile) or die "Can't open $recfile: $!\n"; + open (WFILEpar, ">", $parfile) or die "Can't open $parfile: $!\n"; print WFILErec $recfile,"\n"; READ: print "Reading from file $bnglfile\n"; - if (!open(FH, $bnglfile)) { + if (!open(FH, "<", $bnglfile)) { return("Couldn't read from file $bnglfile: $!"); } $file_dat=[]; @@ -143,7 +143,7 @@ READ: if ($err){ last READ;} $bngdata{$name}=1; if ($name eq "parameters"){ - open (WFILEpar, $parfile) or die "Can't open $parfile: $!\n"; + open (WFILEpar, ">", $parfile) or die "Can't open $parfile: $!\n"; print WFILEpar "begin parameters\n"; for my $line (@{$block_dat}){ my ($entry, $lno)= @{$line}; @@ -334,17 +334,17 @@ READ: print WFILErec " BEGIN COR_NET \n"; - $macrfile= ">macr_${param_prefix}.bngl"; # $macrfile= ">macr_${param_prefix}.bngl"; - $parfile= "", $macrfile) or die "Can't open $macrfile: $!\n"; + open (RFILEpar, "<", $parfile) or die "Can't open $parfile: $!\n"; + open (RFILEspec2, "<", $spec2file) or die "Can't open $spec2file: $!\n"; + open (RFILErules, "<", $rulesfile) or die "Can't open $rulesfile: $!\n"; + open (RFILEobser, "<", $obserfile) or die "Can't open $obserfile: $!\n"; while ($line=) { # *.par print WFILEmacr "$line"; # macr_*.bngl @@ -419,8 +419,8 @@ my ($nm_site, $nm2_site, $param_prefix, $skf, $dpp_site, $lis1h) = @_; # ----------RULES------------- - $rulesfile= ">macr_${param_prefix}.rules"; - open (WFILErules, $rulesfile) or die "Can't open $rulesfile: $!\n"; + $rulesfile= "macr_${param_prefix}.rules"; + open (WFILErules, ">", $rulesfile) or die "Can't open $rulesfile: $!\n"; print WFILErules "begin reaction_rules\n"; @lis1 = (); &del_blank(\@reac1tion,\@lis1); @@ -935,8 +935,8 @@ my (%m_spe,%s_spe,%rea_l,%rea_r,%rrea_l,%rrea_r,%dubl); %ts3n = (); @species = (); - $recfile= " ) { if ( ($line =~ / MINIMUM set/) ) { @@ -1395,8 +1395,8 @@ my ($nm_site, $nm2_site, $param_prefix, $skf, $dpp_site, $lis1h) = @_; my ($ii); my (@lis1); - $obserfile= ">macr_${param_prefix}.obser"; # $obserfile= ">macr_${param_prefix}.obser"; - open (WFILEobser, $obserfile) or die "Can't open $obserfile: $!\n"; + $obserfile= "macr_${param_prefix}.obser"; # $obserfile= "macr_${param_prefix}.obser"; + open (WFILEobser, ">", $obserfile) or die "Can't open $obserfile: $!\n"; print WFILEobser "begin observables\n"; @lis1 = (); @@ -1420,8 +1420,8 @@ close(WFILEobser); # --------- SPECIES2------------- - $spec2file= ">macr_${param_prefix}.spec2"; # $spec2file= ">macr_${param_prefix}.spec2"; - open (WFILEspec2, $spec2file) or die "Can't open $spec2file: $!\n"; + $spec2file= "macr_${param_prefix}.spec2"; # $spec2file= "macr_${param_prefix}.spec2"; + open (WFILEspec2, ">", $spec2file) or die "Can't open $spec2file: $!\n"; print WFILEspec2 "begin species\n"; # *.spec2 @lis1 = (); diff --git a/bng2/Perl2/RateLaw.pm b/bng2/Perl2/RateLaw.pm index cde31a2c..2ab662e5 100644 --- a/bng2/Perl2/RateLaw.pm +++ b/bng2/Perl2/RateLaw.pm @@ -762,7 +762,6 @@ sub toCVodeString my $rl = shift @_; my $stat_factor = shift @_; my $reactants = shift @_; - my $rrefs = shift @_; my $plist = @_ ? shift @_ : undef; my $conv_expr = @_ ? shift @_ : undef; # expression for unit conversions @@ -822,7 +821,7 @@ sub toCVodeString my $fcn = $fcn_param->Ref; # add references to the expressions and observables arrays - my $fcn_str = $fcn->toCVodeString( $plist, {'fcn_mode' => 'call', 'rrefs' => $rrefs, 'reactants' => $reactants}); + my $fcn_str = $fcn->toCVodeString( $plist, {'fcn_mode' => 'call', 'reactants' => $reactants}); if ($fcn_str =~ /^could not find/ || $fcn_str =~ /^ratelaw depends on/) { return $fcn_str; } push @rl_terms, $fcn_str; @@ -852,7 +851,6 @@ sub toMatlabString my $rl = shift @_; my $stat_factor = shift @_; my $reactants = shift @_; - my $rrefs = shift @_; my $plist = @_ ? shift @_ : undef; my $conv_expr = @_ ? shift @_ : undef; # expression for unit conversions @@ -907,7 +905,7 @@ sub toMatlabString my $fcn = $fcn_param->Ref; - my $fcn_str = $fcn->toMatlabString( $plist, {'fcn_mode' => 'call', 'rrefs' => $rrefs, 'reactants' => $reactants}); + my $fcn_str = $fcn->toMatlabString( $plist, {'fcn_mode' => 'call', 'reactants' => $reactants}); if ($fcn_str =~ /^could not find/ || $fcn_str =~ /^ratelaw depends on/) { return $fcn_str; } push @rl_terms, $fcn_str; diff --git a/bng2/Perl2/Rxn.pm b/bng2/Perl2/Rxn.pm index fb56dbd9..82f52338 100644 --- a/bng2/Perl2/Rxn.pm +++ b/bng2/Perl2/Rxn.pm @@ -269,13 +269,9 @@ sub getCVodeRate if ($convert_units) { ($conv_expr, $comp_name, $err) = $rxn->get_intensive_to_extensive_units_conversion($BNGModel::GLOBAL_MODEL); } - # get reference to RxnRule RRef hash (TODO: may be obsolete) - my $rrefs = undef; - if ( $rxn->RxnRule ) - { $rrefs = $rxn->RxnRule->RRefs; } # get ratelaw string my $sf = ($rxn->RxnRule && $rxn->RxnRule->TotalRate) ? 1 : $rxn->StatFactor; - return $rxn->RateLaw->toCVodeString( $sf, $rxn->Reactants, $rrefs, $plist, $conv_expr ); + return $rxn->RateLaw->toCVodeString( $sf, $rxn->Reactants, $plist, $conv_expr ); } @@ -296,13 +292,9 @@ sub getMatlabRate if ($convert_units) { ($conv_expr, $comp_name, $err) = $rxn->get_intensive_to_extensive_units_conversion($BNGModel::GLOBAL_MODEL); } - # get reference to RxnRule RRef hash - my $rrefs = undef; - if ( $rxn->RxnRule ) - { $rrefs = $rxn->RxnRule->RRefs; } # get ratelaw string my $sf = ($rxn->RxnRule && $rxn->RxnRule->TotalRate) ? 1 : $rxn->StatFactor; - return $rxn->RateLaw->toMatlabString( $sf, $rxn->Reactants, $rrefs, $plist, $conv_expr ); + return $rxn->RateLaw->toMatlabString( $sf, $rxn->Reactants, $plist, $conv_expr ); } diff --git a/bng2/Perl2/RxnRule.pm b/bng2/Perl2/RxnRule.pm index 99ffdc8e..a3ab6be5 100644 --- a/bng2/Perl2/RxnRule.pm +++ b/bng2/Perl2/RxnRule.pm @@ -2508,15 +2508,6 @@ sub findMap # => can't change a species' location if the molecule composition has changed. elsif ( $mapPattR[$i_pattR] == -2 ) { - # Removed by Justin -- we're being permissive about species syntax. - # # species compartment declaration is invalid without map to product pattern - # if ( defined $rr->Reactants->[$i_pattR]->Compartment ) - # { - # exit_error( - # "Reaction Rule specifies a species compartment for reactant in which" - # ." molecules are removed or added.", $rr->toString() - # ); - # } next; } @@ -2537,16 +2528,6 @@ sub findMap # NOTE: this could potentially define a generic transport. next if ( !defined($compR) and defined($compP) ); - # Removed by Justin -- we're being permissive about species syntax. - # # error if compartment is defined for one, but not both. - # if ( defined($compR) xor defined($compP) ) - # { - # exit_error( - # "Reaction Rule specifies a species compartment for a pattern on one" - # ." side of the reaction but not for the corresponding species on the" - # ." other side of the reaction.", $rr->toString() - # ); - # } # case 3A: compartments are the same. no transport next if ( $compR == $compP ); @@ -2583,25 +2564,6 @@ sub findMap } } - # Removed by Justin -- we're being permissive about species syntax. - # # One last thing: check for invalid compartment specification on Product side - # for ( my $i_pattP = 0; $i_pattP < @{$rr->Products}; $i_pattP++ ) - # { - # # does product pattern not have a valid map to a reactant pattern? - # if ( $mapPattP[$i_pattP] == -2 ) - # { - # # species compartment declaration is invalid without map to product pattern - # if ( defined $rr->Products->[$i_pattP]->Compartment ) - # { - # exit_error( - # "Reaction Rule specifies a species compartment for product from which" - # ." molecules have been removed or added.", $rr->toString() - # ); - # } - # # otherwise okay. nothing to do. - # next; - # } - # } } # done handling species transport @@ -2857,7 +2819,7 @@ sub findMap my $modelname = $BNGModel::GLOBAL_MODEL->Name; my $rulename = $rr->Name; my $filename = join("_",($modelname,$rulename,"StatFactorCalculation")).".txt"; - open($autfile,">".$filename) or die; + open($autfile, '>', $filename) or die; print $autfile $rr->toString()."\n"; print $autfile "Reactants: ".$rg->toString()."\n"; print $autfile "Products: ".$pg->toString()."\n"; diff --git a/bng2/Perl2/SpeciesGraph.pm b/bng2/Perl2/SpeciesGraph.pm index 06cad724..bba179ef 100644 --- a/bng2/Perl2/SpeciesGraph.pm +++ b/bng2/Perl2/SpeciesGraph.pm @@ -3435,7 +3435,6 @@ sub findMaps # # should this be in Map module? # REVISED by justinshogg@gmail.com 19feb2009 -# TODO: findMaps respects molecule and component tags, but does not attempt to reconile Pattern tags. { # get species graphs my $sg1 = shift @_; @@ -3541,17 +3540,6 @@ sub buildLabelMap # occurence of identical components for ( my $ic = 0 ; $ic < @clabels ; $ic++ ) { - # NEW CODE. TODO: Figure this out! - # the old code leaves out user tags, which seems sketchy. - # but adding the tags breaks the state inheritance mechanism. - #my $clabel = $clabels[$ic]; - ## prefix component label with molecule label to ensure that - ## molecule and component maps are compatible - #$clabel = $mlabel . '|' . $clabel . '_'; - #$clabel .= ++$labels{$clabel}; - #$labelmap->{$clabel} = "$im.$ic"; - - # OLD CODE: my $clabel = $clabels[$ic]; unless ( $clabel =~ /^\%/ ) { diff --git a/bng2/Perl2/Visualization/VisOptParse.pm b/bng2/Perl2/Visualization/VisOptParse.pm index 972dc911..60c7d76f 100644 --- a/bng2/Perl2/Visualization/VisOptParse.pm +++ b/bng2/Perl2/Visualization/VisOptParse.pm @@ -8,7 +8,7 @@ no warnings 'redefine'; sub parseOpts { my $file = shift @_; - open FILE,$file or die "Cannot open file ".$file."\n"; + open(FILE, "<", $file) or die "Cannot open file ".$file."\n"; my @lines = ; close FILE; print "Processing Opts file ".$file."\n"; diff --git a/bng2/Perl2/XML/TreePP.pm b/bng2/Perl2/XML/TreePP.pm index 04f85955..a4eaf401 100644 --- a/bng2/Perl2/XML/TreePP.pm +++ b/bng2/Perl2/XML/TreePP.pm @@ -1106,7 +1106,7 @@ sub write_raw_xml { my $self = shift; my $file = shift; my $fh = Symbol::gensym(); - open( $fh, ">$file" ) or return $self->die( "$! - $file" ); + open( $fh, '>', $file ) or return $self->die( "$! - $file" ); print $fh @_; close($fh); } @@ -1115,7 +1115,7 @@ sub read_raw_xml { my $self = shift; my $file = shift; my $fh = Symbol::gensym(); - open( $fh, $file ) or return $self->die( "$! - $file" ); + open( $fh, '<', $file ) or return $self->die( "$! - $file" ); local $/ = undef; my $text = <$fh>; close($fh); diff --git a/bng2/Perl2/sample_vars_2.pl b/bng2/Perl2/sample_vars_2.pl index 8cbd6f85..ae8e98e1 100644 --- a/bng2/Perl2/sample_vars_2.pl +++ b/bng2/Perl2/sample_vars_2.pl @@ -240,7 +240,7 @@ #} # Read file -open(IN, $file) or die "Couldn't open file $file: $?\n"; +open(IN, '<', $file) or die "Couldn't open file $file: $?\n"; my $script = ""; while ( my $line = ) { diff --git a/bng2/Perl2/verify.pl b/bng2/Perl2/verify.pl index 8291b612..ae952266 100755 --- a/bng2/Perl2/verify.pl +++ b/bng2/Perl2/verify.pl @@ -60,7 +60,7 @@ # read first data file -open (RFILE, $rfilename) or die $INDENT . "$0 ERROR: can't open $rfilename: $!\n"; +open (RFILE, '<', $rfilename) or die $INDENT . "$0 ERROR: can't open $rfilename: $!\n"; my @data = (); my @times = (); while ( my $line = ) @@ -82,7 +82,7 @@ # read second data file -open (RFILE, $rfilename2) or die $INDENT . "$0 ERROR: can't open $rfilename2: $!!!\n"; +open (RFILE, '<', $rfilename2) or die $INDENT . "$0 ERROR: can't open $rfilename2: $!!!\n"; my @data2=(); my @times2=(); while ( my $line = ) diff --git a/parsers/BipartiteGraph/bipartiteGraph.py b/parsers/BipartiteGraph/bipartiteGraph.py index c1b6d93b..3b8ced66 100755 --- a/parsers/BipartiteGraph/bipartiteGraph.py +++ b/parsers/BipartiteGraph/bipartiteGraph.py @@ -451,8 +451,8 @@ def writeDot(patternDict,transformationDict,atomicPatternAnnotations,transformat #n_nodes_to_leave = int(ceil((n_reactants - n_transforms)/2)) #print n_reactants, n_transforms, n_nodes_to_leave - patt_id_list = [x for patt,x in sorted(patternDict.items())] - tr_id_list = [x for tr,x in sorted(transformationDict.items())] + patt_id_list = [patternDict[patt] for patt in sorted(patternDict.keys())] + tr_id_list = [transformationDict[tr] for tr in sorted(transformationDict.keys())] #for idx, item in enumerate(tr_id_list): # f.write('{ rank = same; \"r'+ str(patt_id_list[n_nodes_to_leave+idx]) + '\"; \"t' + str(tr_id_list[idx]) + '\"; }\n') #print idx diff --git a/parsers/BipartiteGraph/bpgMaps.py b/parsers/BipartiteGraph/bpgMaps.py index d064158f..cf68e416 100644 --- a/parsers/BipartiteGraph/bpgMaps.py +++ b/parsers/BipartiteGraph/bpgMaps.py @@ -139,50 +139,37 @@ def __init__(self,atomizedrules,patterns,transformations,transformationpairs,irr for idx,ir in enumerate(irrs): self.irr[ir] = self.tp[ir] + self._str_maps = { + 'p': {str(k): v for k, v in self.p.items()}, + 't': {str(k): v for k, v in self.t.items()}, + 'tp': {str(k): v for k, v in self.tp.items()}, + 'r': {str(k): v for k, v in self.r.items()}, + 'irr': {str(k): v for k, v in self.irr.items()} + } + + self._rev_maps = { + 'p': {v: k for k, v in self.p.items()}, + 't': {v: k for k, v in self.t.items()}, + 'tp': {v: k for k, v in self.tp.items()}, + 'r': {v: k for k, v in self.r.items()}, + 'irr': {v: k for k, v in self.irr.items()} + } + + def getIdx(self,elemtype,string): - # ⚡ Bolt: Use simple for loop instead of generator expression to avoid generator initialization overhead, providing a much faster O(1) early exit - if elemtype == 'p': - for x,idx in list(self.p.items()): - if str(x)==string: return idx - if elemtype == 't': - for x,idx in list(self.t.items()): - if str(x)==string: return idx - if elemtype == 'tp': - for x,idx in list(self.tp.items()): - if str(x)==string: return idx - if elemtype == 'r': - for x,idx in list(self.r.items()): - if str(x)==string: return idx - if elemtype == 'irr': - for x,idx in list(self.irr.items()): - if str(x)==string: return idx - return None + # ⚡ Bolt: Use pre-computed O(1) hash map lookup instead of O(N) loops or generators + return self._str_maps.get(elemtype, {}).get(string) def getElement(self,elemtype,idx1): - # ⚡ Bolt: Use simple for loop instead of generator expression to avoid generator initialization overhead, providing a much faster O(1) early exit - if elemtype == 'p': - for x,idx in list(self.p.items()): - if idx==idx1: return x - if elemtype == 't': - for x,idx in list(self.t.items()): - if idx==idx1: return x - if elemtype == 'tp': - for x,idx in list(self.tp.items()): - if idx==idx1: return x - if elemtype == 'r': - for x,idx in list(self.r.items()): - if idx==idx1: return x - if elemtype == 'irr': - for x,idx in list(self.irr.items()): - if idx==idx1: return x - return None + # ⚡ Bolt: Use pre-computed O(1) hash map lookup instead of O(N) loops or generators + return self._rev_maps.get(elemtype, {}).get(idx1) def getString(self,elemtype,idx1): return str(self.getElement(elemtype,idx1)) def printDict(self,elemtype,someDict,sortbywhat): - tuples = [(self.getString(elemtype,x),y) for x,y in list(someDict.items())] + tuples = [(self.getString(elemtype,x),y) for x,y in someDict.items()] if sortbywhat == 'value': tuples = sorted(tuples,key=lambda x: x[1]) return "\n".join([":".join([str(x) for x in z]) for z in tuples]) @@ -342,7 +329,7 @@ def __init__(self,dictNames,tr_map): self.tp2p_forwardcontext = list(tp2p_forwardcontext_opt) self.tp2p_reversecontext = list(tp2p_reversecontext_opt) - syndel_list = [(tp_id,t_id,dictNames.getElement('t',t_id).action) for tp_id,t_id in list(self.tp2t_forward.items())+list(self.tp2t_reverse.items()) ] + syndel_list = [(tp_id,t_id,dictNames.getElement('t',t_id).action) for tp_id,t_id in itertools.chain(self.tp2t_forward.items(), self.tp2t_reverse.items()) ] t_to_syndel = {} for t_id, p_id in tr_map.t2p_syndelcontext: @@ -591,10 +578,10 @@ def __init__(self,trace,tracetype): self._set = set(trace) def __str__(self): - return "->".join([str(x) for x in self.trace]) + return "->".join(str(x) for x in self.trace) def toString(self,names): - return "->".join([str(names.getElement(self.type,x)) for x in self.trace]) + return "->".join(str(names.getElement(self.type,x)) for x in self.trace) def getLast(self): return self.trace[-1] @@ -908,17 +895,17 @@ def writeJSON(names,all_maps,annot): # Getting the node elements # A node for each rule nodes = [] - for rule,idx in list(names.r.items()): + for rule,idx in names.r.items(): temp = rule.getJSON() temp.update({"idx":idx,"annot":annot.r[idx]}) nodes.append(temp) # A node for each pattern - for patt,idx in list(names.p.items()): + for patt,idx in names.p.items(): temp = patt.getJSON() temp.update({"idx":idx,"annot":annot.p[idx]}) nodes.append(temp) # A node for each transformation (how to deal with irreversibles) - for tr,idx in list(names.t.items()): + for tr,idx in names.t.items(): temp = tr.getJSON() temp.update({"idx":idx,"annot":annot.p[idx]}) if idx in [str(x) for x in names.irr]: @@ -928,7 +915,7 @@ def writeJSON(names,all_maps,annot): nodes.append(temp) # A node for each transformation pair - for tp,idx in list(names.tp.items()): + for tp,idx in names.tp.items(): temp = tp.getJSON() temp.update({"idx":idx,"annot":annot.tp[idx]}) nodes.append(temp) @@ -984,19 +971,19 @@ def writeJSON(names,all_maps,annot): def listify(set1): return [list(x) for x in list(set1)] def listify2(dict1): - return [ [x,y] for x,y in list(dict1.items())] + return [ [x,y] for x,y in dict1.items()] def unq(list1): return list(set(list1)) def combineLists(listoflists): - return reduce(lambda x,y: x+y,listoflists) + return list(itertools.chain.from_iterable(listoflists)) def printDict(somedict): return "\n".join(sorted([str(x)+":"+str(y) for x,y in sorted(somedict.items())])) def defaultDict(somelist,defaultval): - return dict([x,defaultval] for x in somelist) + return {x: defaultval for x in somelist} def assignVal(somedict,somelistofkeys,val): tempdict = somedict diff --git a/parsers/BipartiteGraph/readBNGXML.py b/parsers/BipartiteGraph/readBNGXML.py index b164cb7c..53c47506 100755 --- a/parsers/BipartiteGraph/readBNGXML.py +++ b/parsers/BipartiteGraph/readBNGXML.py @@ -137,8 +137,20 @@ def parseMolecules(molecules): def parseXML(xmlFile): parser = etree.XMLParser(resolve_entities=False, no_network=True) doc = etree.parse(xmlFile, parser) - molecules = doc.findall('.//{http://www.sbml.org/sbml/level3}MoleculeType') - rules = doc.findall('.//{http://www.sbml.org/sbml/level3}ReactionRule') + + model = doc.getroot().find('{http://www.sbml.org/sbml/level3}model') + molecules = [] + rules = [] + parameters = [] + observables = [] + if model is not None: + lom = model.find('{http://www.sbml.org/sbml/level3}ListOfMoleculeTypes') + if lom is not None: + molecules = lom.findall('{http://www.sbml.org/sbml/level3}MoleculeType') + lor = model.find('{http://www.sbml.org/sbml/level3}ListOfReactionRules') + if lor is not None: + rules = lor.findall('{http://www.sbml.org/sbml/level3}ReactionRule') + ruleDescription = [] moleculeList = [] for molecule in molecules: diff --git a/parsers/ContactMap/createGraph.py b/parsers/ContactMap/createGraph.py index 4eb8ebe8..3a50740c 100644 --- a/parsers/ContactMap/createGraph.py +++ b/parsers/ContactMap/createGraph.py @@ -307,14 +307,8 @@ def processBNGL(bngl,center,context,product): reactionCenter=center, context=context, products=product) def bngl2xml(bnglFile): - - bngconsole = pexpect.spawn('bngdev --console') - bngconsole.expect('BNG>') - bngconsole.sendline('load {0}'.format(bnglFile)) - bngconsole.expect('BNG>') - bngconsole.sendline('action writeXML()') - bngconsole.expect('BNG>') - bngconsole.close() + import subprocess + subprocess.call(['bngdev', bnglFile, '--xml'], shell=False) def main(fileName): diff --git a/parsers/ContactMap/server.py b/parsers/ContactMap/server.py index e5d03845..8ac424fc 100644 --- a/parsers/ContactMap/server.py +++ b/parsers/ContactMap/server.py @@ -58,7 +58,7 @@ def bipartite(self, bbnglFile,returnType,center,context,product): bnglFile = bbnglFile.data with open('temp{0}.bngl'.format(counter),'w') as f: f.write(bnglFile) - xmlFile = self.bngl2xml('temp{0}.bngl'.format(counter)) + xmlFile = self._bngl2xml('temp{0}.bngl'.format(counter)) createGraph.processBNGL('temp{0}.xml'.format(counter),center,context,product) with open('temp{0}.xml.dot'.format(counter),'rb') as f: dot = f.read() @@ -77,21 +77,21 @@ def bipartite(self, bbnglFile,returnType,center,context,product): def getTransformations(self,bbnglFile): pass - def bngl2xml(self,bnglFile): - - bngconsole = pexpect.spawn('bngdev --console') - bngconsole.expect('BNG>') - bngconsole.sendline('load {0}'.format(bnglFile)) - bngconsole.expect('BNG>') - bngconsole.sendline('action writeXML()') - bngconsole.expect('BNG>') - bngconsole.close() + def _bngl2xml(self,bnglFile): + subprocess.call(['bngdev', bnglFile, '--xml'], shell=False) +import argparse + if __name__ == '__main__': - server = SimpleXMLRPCServer(("10.253.98.102", 9100), requestHandler=RequestHandler) + parser = argparse.ArgumentParser(description="Start the Bipartite XML-RPC Server") + parser.add_argument('--host', type=str, default=os.environ.get('HOST', '127.0.0.1'), help='Host IP address to bind to') + parser.add_argument('--port', type=int, default=int(os.environ.get('PORT', 9100)), help='Port to bind to') + args = parser.parse_args() + + server = SimpleXMLRPCServer((args.host, args.port), requestHandler=RequestHandler) server.register_introspection_functions() server.register_instance(BipartiteServer()) server.serve_forever() diff --git a/parsers/ContactMap/tests/test_server.py b/parsers/ContactMap/tests/test_server.py index 97c6d1a6..7b29b81e 100644 --- a/parsers/ContactMap/tests/test_server.py +++ b/parsers/ContactMap/tests/test_server.py @@ -20,8 +20,8 @@ class TestServer(unittest.TestCase): @patch('server.glob.glob') @patch('builtins.open', new_callable=mock_open, read_data=b'dummy') @patch('server.createGraph') - @patch.object(server.BipartiteServer, 'bngl2xml') - def test_bipartite_oserror_path(self, mock_bngl2xml, mock_createGraph, mock_file, mock_glob, mock_remove): + @patch.object(server.BipartiteServer, '_bngl2xml') + def test_bipartite_oserror_path(self, mock__bngl2xml, mock_createGraph, mock_file, mock_glob, mock_remove): """ Tests the error path in bipartite where os.remove throws an OSError, verifying that it is caught and ignored gracefully. diff --git a/parsers/utils/consoleCommands.py b/parsers/utils/consoleCommands.py index a537326f..f2446616 100644 --- a/parsers/utils/consoleCommands.py +++ b/parsers/utils/consoleCommands.py @@ -24,7 +24,7 @@ def getBngExecutable(): def bngl2xml(bnglFile,timeout=60): try: - bngconsole = pexpect.spawn('{0} --console'.format(getBngExecutable()),timeout=timeout) + bngconsole = pexpect.spawn(getBngExecutable(), ['--console'], timeout=timeout) bngconsole.expect('BNG>') bngconsole.sendline('load {0}'.format(bnglFile)) bngconsole.expect('BNG>') @@ -38,7 +38,7 @@ def bngl2xml(bnglFile,timeout=60): def bngl2sbml(bnglFile,timeout=60): try: - bngconsole = pexpect.spawn('{0} --console'.format(getBngExecutable()),timeout=timeout) + bngconsole = pexpect.spawn(getBngExecutable(), ['--console'], timeout=timeout) bngconsole.expect('BNG>') bngconsole.sendline('load {0}'.format(bnglFile)) bngconsole.expect('BNG>') @@ -53,7 +53,7 @@ def bngl2sbml(bnglFile,timeout=60): subprocess.call([killall_path, 'bngdev'], shell=False) # nosec def correctness(bnglFile): - bngconsole = pexpect.spawn('{0} --console'.format(getBngExecutable())) + bngconsole = pexpect.spawn(getBngExecutable(), ['--console']) bngconsole.expect('BNG>') bngconsole.sendline('load {0}'.format(bnglFile)) bngconsole.expect('BNG>') @@ -65,7 +65,7 @@ def correctness(bnglFile): def writeNetwork(bnglFile): - bngconsole = pexpect.spawn('{0} --console'.format(getBngExecutable())) + bngconsole = pexpect.spawn(getBngExecutable(), ['--console']) bngconsole.expect('BNG>') bngconsole.sendline('load {0}'.format(bnglFile)) bngconsole.expect('BNG>') @@ -78,7 +78,7 @@ def generateGraph(bnglFile,graphType): directory = os.sep.join(bnglFile.split(os.sep)[:-1]) os.chdir(directory) print(directory) - bngconsole = pexpect.spawn('{0} --console'.format(getBngExecutable())) + bngconsole = pexpect.spawn(getBngExecutable(), ['--console']) bngconsole.expect('BNG>') bngconsole.sendline('load {0}'.format(bnglFile)) bngconsole.expect('BNG>') diff --git a/parsers/utils/readBNGXML.py b/parsers/utils/readBNGXML.py index eaac4185..efddb7b3 100644 --- a/parsers/utils/readBNGXML.py +++ b/parsers/utils/readBNGXML.py @@ -36,20 +36,30 @@ def findBond(bondDefinitions, component): def createMolecule(molecule, bonds): nameDict = {} - mol = st.Molecule(molecule.get('name'),molecule.get('id')) - if molecule.get('compartment') not in ['',None]: - mol.setCompartment(molecule.get('compartment')) - nameDict[molecule.get('id')] = molecule.get('name') + mol_id = molecule.get('id') + mol_name = molecule.get('name') + mol = st.Molecule(mol_name, mol_id) + mol_comp = molecule.get('compartment') + if mol_comp not in ['', None]: + mol.setCompartment(mol_comp) + nameDict[mol_id] = mol_name listOfComponents = _fast_find(molecule, 'ListOfComponents') if listOfComponents != None: for element in listOfComponents: - component = st.Component(element.get('name'),element.get('id')) - nameDict[element.get('id')] = element.get('name') - if element.get('numberOfBonds') in ['+','?']: - component.addBond(element.get('numberOfBonds')) - elif element.get('numberOfBonds') != '0': - component.addBond(findBond(bonds, element.get('id'))) - state = element.get('state') if element.get('state') != None else '' + elem_id = element.get('id') + elem_name = element.get('name') + elem_bonds = element.get('numberOfBonds') + elem_state = element.get('state') + + component = st.Component(elem_name, elem_id) + nameDict[elem_id] = elem_name + + if elem_bonds in ['+', '?']: + component.addBond(elem_bonds) + elif elem_bonds != '0': + component.addBond(findBond(bonds, elem_id)) + + state = elem_state if elem_state != None else '' component.states.append(state) component.activeState = state mol.addComponent(component) @@ -180,12 +190,28 @@ def parseComponent(component): def parseXML(xmlFile): parser = etree.XMLParser(resolve_entities=False, no_network=True) doc = etree.parse(xmlFile, parser) - molecules = doc.findall('.//{http://www.sbml.org/sbml/level3}MoleculeType') - rules = doc.findall('.//{http://www.sbml.org/sbml/level3}ReactionRule') + + model = doc.getroot().find('{http://www.sbml.org/sbml/level3}model') + molecules = [] + rules = [] + parameters = [] + observables = [] + if model is not None: + lom = model.find('{http://www.sbml.org/sbml/level3}ListOfMoleculeTypes') + if lom is not None: + molecules = lom.findall('{http://www.sbml.org/sbml/level3}MoleculeType') + lor = model.find('{http://www.sbml.org/sbml/level3}ListOfReactionRules') + if lor is not None: + rules = lor.findall('{http://www.sbml.org/sbml/level3}ReactionRule') + lop = model.find('{http://www.sbml.org/sbml/level3}ListOfParameters') + if lop is not None: + parameters = lop.findall('{http://www.sbml.org/sbml/level3}Parameter') + loo = model.find('{http://www.sbml.org/sbml/level3}ListOfObservables') + if loo is not None: + observables = loo.findall('{http://www.sbml.org/sbml/level3}Observable') + ruleDescription = [] moleculeList = [] - - parameters = doc.findall('.//{http://www.sbml.org/sbml/level3}Parameter') parameterDict = {} for parameter in parameters: parameterDict[parameter.get('id')] = parameter.get('value') @@ -205,7 +231,12 @@ def parseXML(xmlFile): def getNumObservablesXML(xmlFile): parser = etree.XMLParser(resolve_entities=False, no_network=True) doc = etree.parse(xmlFile, parser) - observables = doc.findall('.//{http://www.sbml.org/sbml/level3}Observable') + observables = [] + model = doc.getroot().find('{http://www.sbml.org/sbml/level3}model') + if model is not None: + loo = model.find('{http://www.sbml.org/sbml/level3}ListOfObservables') + if loo is not None: + observables = loo.findall('{http://www.sbml.org/sbml/level3}Observable') return len(observables) if __name__ == "__main__": diff --git a/parsers/utils/smallStructures.py b/parsers/utils/smallStructures.py index 50627f63..1e4b9c69 100644 --- a/parsers/utils/smallStructures.py +++ b/parsers/utils/smallStructures.py @@ -80,7 +80,7 @@ def addActionList(self,actionList): def __str__(self): label = f"{self.label}: " if self.label != '' else "" arrow = ' <-> ' if self.bidirectional else ' -> ' - return f"{label}{' + '.join(str(x) for x in self.reactants)}{arrow}{' + '.join(str(x) for x in self.products)} {','.join(self.rates)}" + return f"{label}{' + '.join([str(x) for x in self.reactants])}{arrow}{' + '.join([str(x) for x in self.products])} {','.join(self.rates)}" class Species: def __init__(self): self.molecules = [] @@ -270,7 +270,7 @@ def append(self,species): def __str__(self): self.molecules.sort(key= lambda molecule: molecule.name) - name= '.'.join(x.toString() for x in self.molecules) + name= '.'.join([x.toString() for x in self.molecules]) ''' name = name.replace('~','') @@ -285,7 +285,7 @@ def __str__(self): return name def str2(self): - return '.'.join(x.str2() for x in self.molecules) + return '.'.join([x.str2() for x in self.molecules]) def reset(self): for element in self.molecules: @@ -370,11 +370,14 @@ def graphVizGraph(self,graph,identifier,layout='LR',options={}): speciesDictionary.update(compDictionary) for bond in self.bonds: - if bond[0] in speciesDictionary and bond[1] in speciesDictionary: - if layout == 'RL': - graph.add_edge(speciesDictionary[bond[1]],speciesDictionary[bond[0]],dir='none',len=0.1,weight=100) - else: - graph.add_edge(speciesDictionary[bond[0]],speciesDictionary[bond[1]],dir='none',len=0.1,weight=100) + b0 = speciesDictionary.get(bond[0]) + if b0 is not None: + b1 = speciesDictionary.get(bond[1]) + if b1 is not None: + if layout == 'RL': + graph.add_edge(b1, b0, dir='none', len=0.1, weight=100) + else: + graph.add_edge(b0, b1, dir='none', len=0.1, weight=100) return speciesDictionary @@ -488,7 +491,7 @@ def contains(self,componentName): def __str__(self): self.components = sorted(self.components,key = lambda st:st.name) - components_str = '(' + ','.join(str(x) for x in self.components) + ')' if self.components else '' + components_str = '(' + ','.join([str(x) for x in self.components]) + ')' if self.components else '' compartment_str = '@' + self.compartment if self.compartment else '' # ⚡ Bolt: Use single f-string to prevent intermediate string allocations return f"{self.name}{components_str}{compartment_str}" @@ -498,7 +501,7 @@ def toString(self): def str2(self): self.components.sort() - return self.name + '(' + ','.join(x.str2() for x in self.components) + ')' + return self.name + '(' + ','.join([x.str2() for x in self.components]) + ')' def str3(self): return self.name + '(' + self.components[0].name + ')' @@ -636,7 +639,7 @@ def setActiveState(self,state): return True def getRuleStr(self): - bonds_str = '!' + '!'.join(str(x) for x in self.bonds) if self.bonds else '' + bonds_str = '!' + '!'.join([str(x) for x in self.bonds]) if self.bonds else '' state_str = '~' + self.activeState if self.activeState else '' # ⚡ Bolt: Use single f-string to prevent intermediate string allocations return f"{self.name}{bonds_str}{state_str}" @@ -651,8 +654,8 @@ def __str__(self): return self.getRuleStr() def str2(self): - bonds_str = '!' + '!'.join(str(x) for x in self.bonds) if self.bonds else '' - states_str = '~' + '~'.join(str(x) for x in self.states) if self.states else '' + bonds_str = '!' + '!'.join([str(x) for x in self.bonds]) if self.bonds else '' + states_str = '~' + '~'.join([str(x) for x in self.states]) if self.states else '' # ⚡ Bolt: Use single f-string to prevent intermediate string allocations return f"{self.name}{bonds_str}{states_str}" diff --git a/parsers/utils/tests/test_consoleCommands.py b/parsers/utils/tests/test_consoleCommands.py index b7291cf9..14dffb1f 100644 --- a/parsers/utils/tests/test_consoleCommands.py +++ b/parsers/utils/tests/test_consoleCommands.py @@ -21,40 +21,32 @@ class MockPexpectTimeout(Exception): from parsers.utils.consoleCommands import correctness, bngl2xml class TestConsoleCommands(unittest.TestCase): - @patch('parsers.utils.consoleCommands.shutil.which', return_value='/mock/path/killall') - @patch('parsers.utils.consoleCommands.pexpect.spawn') + @patch('parsers.utils.consoleCommands.subprocess.run') @patch('parsers.utils.consoleCommands.getBngExecutable', return_value='dummy_bng') - def test_correctness_success(self, mock_getBngExecutable, mock_spawn, mock_which): - mock_bngconsole = MagicMock() - mock_bngconsole.before = 'Some output without error' - mock_spawn.return_value = mock_bngconsole - + def test_correctness_success(self, mock_getBngExecutable, mock_run): + mock_run.return_value = MagicMock(stdout='Some output without error', stderr='') result = correctness('dummy.bngl') self.assertTrue(result) - @patch('parsers.utils.consoleCommands.shutil.which', return_value='/mock/path/killall') - @patch('parsers.utils.consoleCommands.pexpect.spawn') + @patch('parsers.utils.consoleCommands.subprocess.run') @patch('parsers.utils.consoleCommands.getBngExecutable', return_value='dummy_bng') - def test_correctness_error(self, mock_getBngExecutable, mock_spawn, mock_which): - mock_bngconsole = MagicMock() - mock_bngconsole.before = 'Some output with ERROR inside' - mock_spawn.return_value = mock_bngconsole - + def test_correctness_error(self, mock_getBngExecutable, mock_run): + mock_run.return_value = MagicMock(stdout='Some output with ERROR inside', stderr='') result = correctness('dummy.bngl') self.assertFalse(result) @patch('parsers.utils.consoleCommands.subprocess.call') @patch('parsers.utils.consoleCommands.shutil.which', return_value='/mock/path/killall') - @patch('parsers.utils.consoleCommands.pexpect.spawn') - def test_bngl2xml_timeout(self, mock_spawn, mock_which, mock_subprocess_call): - # Configure the mock to raise pexpect.TIMEOUT when spawn is called - mock_spawn.side_effect = mock_pexpect.TIMEOUT('Timeout occurred') + @patch('parsers.utils.consoleCommands.getBngExecutable', return_value='dummy_bng') + def test_bngl2xml_timeout(self, mock_getBngExecutable, mock_which, mock_subprocess_call): + # Configure the mock to raise subprocess.TimeoutExpired on the first call (bng2 execution), + # and do nothing (return None) on the second call (killall) + mock_subprocess_call.side_effect = [subprocess.TimeoutExpired(cmd='dummy', timeout=60), None] - # Call the function bngl2xml('dummy.bngl') - # Assert that subprocess.call was called with the correct arguments - mock_subprocess_call.assert_called_once_with(['/mock/path/killall', 'bngdev'], shell=False) + # Assert that subprocess.call was called to kill bngdev + mock_subprocess_call.assert_any_call(['/mock/path/killall', 'bngdev'], shell=False) if __name__ == '__main__': unittest.main() diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 559fd26a..07f8ab9a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -108,7 +108,7 @@ target_link_libraries(bng_engine PUBLIC bng_ast bng_core antlr4_static - sundials_cvode_static + sundials_cvodes_static sundials_nvecserial_static sundials_sunlinsoldense_static sundials_sunlinsolspgmr_static diff --git a/src/actions/ActionDispatch.cpp b/src/actions/ActionDispatch.cpp index 041ed3da..5f7c5f66 100644 --- a/src/actions/ActionDispatch.cpp +++ b/src/actions/ActionDispatch.cpp @@ -74,10 +74,37 @@ std::string stripQuotes(const std::string& text) { return text; } + +void validateNFSimPath(const std::string& path) { + if (path.empty()) return; + + // Ensure the executable is named NFsim or NFsim.exe to prevent arbitrary command execution + std::filesystem::path p(path); + std::string filename = p.filename().string(); + if (filename != "NFsim" && filename != "NFsim.exe") { + throw std::runtime_error("Invalid nfsim_exec path: executable must be named 'NFsim' or 'NFsim.exe'"); + } +} + +void validateNumeric(const std::string& arg, const std::string& argName) { + if (arg.empty()) return; + + // Ensure the argument is a valid number (integer or float) + try { + size_t pos = 0; + std::stod(arg, &pos); + if (pos != arg.length()) { + throw std::invalid_argument("not fully parsed"); + } + } catch (...) { + throw std::runtime_error("Invalid argument for NFSim parameter '" + argName + "': must be numeric"); + } +} + std::string lowercase(std::string value) { - std::transform(value.begin(), value.end(), value.begin(), [](unsigned char c) { - return static_cast(std::tolower(c)); - }); + for (char& c : value) { + if (c >= 'A' && c <= 'Z') c += ('a' - 'A'); + } return value; } @@ -1295,6 +1322,12 @@ void ActionDispatch::execute(ast::Model& model, const std::filesystem::path& sou cmd += " -ss \"" + speciesPath.string() + "\""; } + + validateNumeric(tEnd, "t_end"); + validateNumeric(nSteps, "n_steps"); + if (!seedText.empty()) validateNumeric(seedText, "seed"); + validateNFSimPath(nfsimExec); + if (verbose) { std::cerr << "[bng_cpp] Running NFSim: " << cmd << "\n"; } @@ -1408,7 +1441,7 @@ void ActionDispatch::execute(ast::Model& model, const std::filesystem::path& sou } // Optional parameters (Perl BNG2 defaults) - const double bump = parseScalarValue(readArgument(action, "bump", "5"), model); + const double bump = parseScalarValue(readArgument(action, "bump", "0.1"), model); const double atol = parseScalarValue(readArgument(action, "atol", "1e-8"), model); const double rtol = parseScalarValue(readArgument(action, "rtol", "1e-8"), model); const auto nSteps = static_cast(parseScalarValue(readArgument(action, "n_steps", "50"), model)); @@ -1422,6 +1455,82 @@ void ActionDispatch::execute(ast::Model& model, const std::filesystem::path& sou const auto prefix = stripQuotes(readArgument(action, "prefix", sourcePath.stem().string())); const auto outputDir = sourcePath.parent_path(); + // New: method="cvodes" runs ONE forward-sensitivity CVODES + // integration instead of Np+1 finite-difference re-simulations. + // Default stays "fd" so existing scripts are unaffected. + const auto senMethod = lowercase(stripQuotes(readArgument(action, "method", "fd"))); + + if (senMethod == "cvodes") { + // Optional sens_param="k1 k2 k3" (whitespace/comma separated) + // to fit a subset instead of every non-zero parameter. + const auto sensParamText = stripQuotes(readArgument(action, "sens_param", "")); + std::vector sensParams; + if (!sensParamText.empty()) { + std::string cleaned = sensParamText; + std::replace(cleaned.begin(), cleaned.end(), ',', ' '); + std::istringstream iss(cleaned); + std::string tok; + while (iss >> tok) sensParams.push_back(tok); + } else { + for (const auto& param : model.getParameters().all()) { + if (param.getValue() != 0.0) sensParams.push_back(param.getName()); + } + } + + if (sensParams.empty()) { + throw std::runtime_error("LinearParameterSensitivity (method=cvodes): no non-zero parameters to compute sensitivities for"); + } + + // Single shared equilibration (if requested), reused as the + // common starting point for every sensitivity direction -- + // unlike the FD path's re_equil, which separately re-equilibrates + // a perturbed model per parameter and so starts each finite + // difference from a slightly different baseline. + if (initEquil) { + engine::OdeOptions eqOpts; + eqOpts.method = "cvode"; + eqOpts.tEnd = tEquil; + eqOpts.nSteps = nSteps; + eqOpts.atol = atol; + eqOpts.rtol = rtol; + eqOpts.steadyState = true; + eqOpts.steadyStateTol = atol; + engine::OdeIntegrator eqIntegrator(model, *network); + auto eqResult = eqIntegrator.integrate(eqOpts); + if (!eqResult.concentrations.empty()) { + const auto& finalConc = eqResult.concentrations.back(); + for (std::size_t i = 0; i < network->species.size() && i < finalConc.size(); ++i) { + network->species.get(i).setAmount(finalConc[i]); + } + } + if (verbose) { + std::cerr << "[bng_cpp] LinearParameterSensitivity(cvodes): equilibration complete\n"; + } + } + + engine::OdeOptions opts; + opts.method = "cvode"; + opts.tEnd = parseScalarValue(tEndText, model); + opts.nSteps = nSteps; + opts.atol = atol; + opts.rtol = rtol; + opts.sensParams = sensParams; + + engine::OdeIntegrator integrator(model, *network); + auto result = integrator.integrate(opts); + + const auto outputPrefix = (outputDir / prefix).string(); + integrator.writeSensitivityFiles(outputPrefix, result, suffix); + + if (verbose) { + std::cerr << "[bng_cpp] LinearParameterSensitivity(cvodes): wrote " + << sensParams.size() << " parameter(s) worth of .csc/.gsc to " + << outputPrefix << "_.{csc,gsc}\n"; + } + + continue; // skip the finite-difference path below entirely + } + // Save original parameter values std::vector> originalParamValues; for (const auto& param : model.getParameters().all()) { diff --git a/src/ast/MacroBNGModel.cpp b/src/ast/MacroBNGModel.cpp index bf0c6168..98feee7b 100644 --- a/src/ast/MacroBNGModel.cpp +++ b/src/ast/MacroBNGModel.cpp @@ -91,12 +91,16 @@ std::string MacroBNGModel::replaceAll(const std::string& s, const std::string& from, const std::string& to) { if (from.empty()) return s; - std::string result = s; + std::string result; + result.reserve(s.length()); + std::string::size_type last_pos = 0; std::string::size_type pos = 0; - while ((pos = result.find(from, pos)) != std::string::npos) { - result.replace(pos, from.size(), to); - pos += to.size(); + while ((pos = s.find(from, last_pos)) != std::string::npos) { + result.append(s, last_pos, pos - last_pos); + result.append(to); + last_pos = pos + from.length(); } + result.append(s, last_pos, s.length() - last_pos); return result; } @@ -112,10 +116,30 @@ std::string MacroBNGModel::rtrim(const std::string& s) { return std::string(s.begin(), it.base()); } + std::string MacroBNGModel::trim(const std::string& s) { return ltrim(rtrim(s)); } +std::string MacroBNGModel::collapseWhitespace(const std::string& s) { + std::string result; + result.reserve(s.size()); + bool in_space = false; + for (char c : s) { + if (std::isspace(static_cast(c))) { + if (!in_space) { + result.push_back(' '); + in_space = true; + } + } else { + result.push_back(c); + in_space = false; + } + } + return result; +} + + std::string MacroBNGModel::quotemeta(const std::string& s) { // Escape all non-alphanumeric, non-underscore characters for regex use // (mirrors Perl \Q...\E / quotemeta) @@ -213,7 +237,7 @@ MacroBNGModel::read_block_array(const std::string& name) { std::string ename = trimmed.substr(4); // trim and normalize whitespace ename = trim(ename); - ename = std::regex_replace(ename, std::regex("\\s+"), " "); + ename = collapseWhitespace(ename); if (ename != name) { return {{}, errgen("end " + ename + " does not match begin " + name)}; } @@ -377,12 +401,10 @@ std::string MacroBNGModel::pre_macr(const std::string& param_prefix) { // Check for "begin " std::string trimmed = line; trimmed.erase(0, trimmed.find_first_not_of(" \t")); - std::smatch m; - std::regex re_begin("^begin\\s+(.*)"); - if (std::regex_search(trimmed, m, re_begin)) { - std::string name = m[1].str(); + if (trimmed.compare(0, 6, "begin ") == 0 || trimmed.compare(0, 6, "begin\t") == 0) { + std::string name = trimmed.substr(6); name = trim(name); - name = std::regex_replace(name, std::regex("\\s+"), " "); + name = collapseWhitespace(name); auto [block_dat, block_err] = read_block_array(name); if (!block_err.empty()) { @@ -425,12 +447,12 @@ std::string MacroBNGModel::pre_macr(const std::string& param_prefix) { WFILEpar.close(); } // --- molecule_types --- - else if (std::regex_match(name, std::regex("^molecule[_ ]types$"))) { + else if (name == "molecule types" || name == "molecule_types") { // We store the type info via pre_species1 later; just count. std::cout << "Read " << block_dat.size() << " molecule types.\n"; } // --- seed_species --- - else if (std::regex_match(name, std::regex("^seed[_ ]species$"))) { + else if (name == "seed species" || name == "seed_species") { // Perl stores into SeedSpeciesList; we just count. std::cout << "Read " << block_dat.size() << " species.\n"; } @@ -442,21 +464,36 @@ std::string MacroBNGModel::pre_macr(const std::string& param_prefix) { std::cout << "Read " << block_dat.size() << " species.\n"; } // --- reaction_rules --- - else if (std::regex_match(name, std::regex("^reaction[_ ]rules$"))) { + else if (name == "reaction rules" || name == "reaction_rules") { int nerr_count = 0; for (auto& [entry, lno] : block_dat) { std::string str1ing = entry; std::string nam1e; // Check for user-named rule "name: ..." - std::smatch rm; - if (std::regex_search(str1ing, rm, std::regex("^\\s*([^:].*)[:]\\s*"))) { - nam1e = rm[1].str(); - std::cerr << "ERROR 1002 rules (user) name=" << nam1e << "\n"; - return "ERROR 1002: named rules not supported"; - } else if (std::regex_search(str1ing, rm, std::regex("^\\s*(\\d+)\\s+"))) { - nam1e = rm[1].str(); - str1ing = rm.suffix().str(); + // BOLT OPTIMIZATION: Avoid O(N) regex overhead in loop + size_t start = 0; + while (start < str1ing.length() && std::isspace(str1ing[start])) start++; + + bool named = false; + if (start < str1ing.length() && str1ing[start] != ':') { + size_t colon_pos = str1ing.find_last_of(':'); + if (colon_pos != std::string::npos && colon_pos >= start) { + nam1e = str1ing.substr(start, colon_pos - start); + std::cerr << "ERROR 1002 rules (user) name=" << nam1e << "\n"; + return "ERROR 1002: named rules not supported"; + } + } + if (!named && start < str1ing.length() && std::isdigit(str1ing[start])) { + size_t end = start; + while (end < str1ing.length() && std::isdigit(str1ing[end])) end++; + if (end < str1ing.length() && std::isspace(str1ing[end])) { + nam1e = str1ing.substr(start, end - start); + size_t suffix_start = end; + while (suffix_start < str1ing.length() && std::isspace(str1ing[suffix_start])) suffix_start++; + str1ing = str1ing.substr(suffix_start); + } } + str1ing = ltrim(str1ing); reac1tion_.push_back(str1ing); } @@ -471,9 +508,17 @@ std::string MacroBNGModel::pre_macr(const std::string& param_prefix) { for (auto& [entry, lno] : block_dat) { std::string str1ing = entry; // Strip leading numeric index - std::smatch rm; - if (std::regex_search(str1ing, rm, std::regex("^\\s*\\d+\\s+"))) { - str1ing = rm.suffix().str(); + // BOLT OPTIMIZATION: Avoid regex compilation overhead + size_t start = 0; + while (start < str1ing.length() && std::isspace(str1ing[start])) start++; + if (start < str1ing.length() && std::isdigit(str1ing[start])) { + size_t end = start; + while (end < str1ing.length() && std::isdigit(str1ing[end])) end++; + if (end < str1ing.length() && std::isspace(str1ing[end])) { + size_t suffix_start = end; + while (suffix_start < str1ing.length() && std::isspace(str1ing[suffix_start])) suffix_start++; + str1ing = str1ing.substr(suffix_start); + } } obser1vable_.push_back(str1ing); } @@ -489,9 +534,13 @@ std::string MacroBNGModel::pre_macr(const std::string& param_prefix) { } } // Check for action lines like "generate_network" / "simulate" - else if (std::regex_search(trimmed, m, std::regex("^([A-Za-z][^(]*)"))) { - gene1rate_.push_back(line); + else { + // BOLT OPTIMIZATION: Avoid regex compilation overhead + if (!trimmed.empty() && std::isalpha(trimmed[0])) { + gene1rate_.push_back(line); + } } + } // while get_line if (!err.empty()) return err; @@ -586,17 +635,28 @@ void MacroBNGModel::pre_species1(std::map& nm_site, // Strip leading label "name: ..." { - std::smatch m; - if (std::regex_search(entry, m, std::regex("^\\s*([^:].*)[:]\\s*"))) { - name = m[1].str(); - std::cerr << "ERROR 1001 block species (user) name=" << name << "\n"; - return; - } else { + // BOLT OPTIMIZATION: Avoid O(N) regex overhead in loop + size_t start = 0; + while (start < entry.length() && std::isspace(entry[start])) start++; + bool matched = false; + if (start < entry.length() && entry[start] != ':') { + size_t colon_pos = entry.find_last_of(':'); + if (colon_pos != std::string::npos && colon_pos >= start) { + name = entry.substr(start, colon_pos - start); + std::cerr << "ERROR 1001 block species (user) name=" << name << "\n"; + return; + } + } + if (!matched) { // Strip leading numeric index - std::regex_replace(entry, std::regex("^\\s*\\d+\\s+"), ""); - std::smatch rm; - if (std::regex_search(entry, rm, std::regex("^\\s*\\d+\\s+"))) { - entry = rm.suffix().str(); + size_t i = 0; + while (i < entry.size() && std::isspace(entry[i])) i++; + if (i < entry.size() && std::isdigit(entry[i])) { + while (i < entry.size() && std::isdigit(entry[i])) i++; + if (i < entry.size() && std::isspace(entry[i])) { + while (i < entry.size() && std::isspace(entry[i])) i++; + entry = entry.substr(i); + } } } } @@ -608,12 +668,16 @@ void MacroBNGModel::pre_species1(std::map& nm_site, // Extract molecule(site,site,...) patterns // Repeatedly match (...) groups - std::regex re_paren("([\\(])(.*?)([\\)])"); - std::smatch pm; std::string remaining = spec_entry; - while (std::regex_search(remaining, pm, re_paren)) { + size_t search_pos = 0; + while (true) { + size_t open_pos = remaining.find('(', search_pos); + if (open_pos == std::string::npos) break; + size_t close_pos = remaining.find(')', open_pos); + if (close_pos == std::string::npos) break; + // name is everything before the '(' - std::string prefix_str = pm.prefix().str(); + std::string prefix_str = remaining.substr(0, open_pos); // Get the molecule name: last segment after '.' auto dot_pos = prefix_str.rfind('.'); if (dot_pos != std::string::npos) { @@ -623,7 +687,7 @@ void MacroBNGModel::pre_species1(std::map& nm_site, } if (nm_site.find(name) == nm_site.end()) { - std::string inside = pm[2].str(); // contents inside parentheses + std::string inside = remaining.substr(open_pos + 1, close_pos - open_pos - 1); // contents inside parentheses auto sits = split(inside, ','); // Strip modifiers from site names for (auto& s : sits) { @@ -646,7 +710,8 @@ void MacroBNGModel::pre_species1(std::map& nm_site, nm_site[name] = static_cast(sits.size()); } - remaining = pm.suffix().str(); + remaining = remaining.substr(close_pos + 1); + search_pos = 0; // Also strip "name." prefix from remaining if (!remaining.empty() && remaining[0] == '.') { remaining = remaining.substr(1); @@ -676,14 +741,24 @@ void MacroBNGModel::del_blank(const std::vector& str, // Strip trailing whitespace line = rtrim(line); // Collapse internal whitespace to single space - line = std::regex_replace(line, std::regex("\\s+"), " "); + line = collapseWhitespace(line); // Replace single spaces with semicolons line = replaceAll(line, " ", ";"); // If starts with digit followed by ';', strip the leading number - std::smatch m; - if (std::regex_search(line, m, std::regex("^\\d+?;"))) { - line = ";" + m.suffix().str(); + size_t semicolon_pos = line.find(';'); + bool has_leading_num = false; + if (semicolon_pos != std::string::npos && semicolon_pos > 0) { + has_leading_num = true; + for (size_t i = 0; i < semicolon_pos; ++i) { + if (!std::isdigit(static_cast(line[i]))) { + has_leading_num = false; + break; + } + } + } + if (has_leading_num) { + line = ";" + line.substr(semicolon_pos + 1); } else { line = ";" + line; } @@ -819,15 +894,14 @@ void MacroBNGModel::skf0(const std::string& rp1, } // Remove the first occurrence of lnk1 from p1 - p1 = replaceFirst(p1, lnk1, ""); + p1.erase(bang_pos, lnk1.length()); // Find the second occurrence of lnk1 in p1 - if (p1.find(lnk1) == std::string::npos) { + auto bang_pos2 = p1.find(lnk1); + if (bang_pos2 == std::string::npos) { // No matching second link — skip continue; } - - auto bang_pos2 = p1.find(lnk1); std::string lef2 = p1.substr(0, bang_pos2); // Extract skf2, sit2 from lef2 the same way @@ -870,7 +944,7 @@ void MacroBNGModel::skf0(const std::string& rp1, } // Remove the second occurrence of lnk1 from p1 - p1 = replaceFirst(p1, lnk1, ""); + p1.erase(bang_pos2, lnk1.length()); // add_skf both directions add_skf(skf1, skf2, sit1, skf, nm2_site); @@ -960,21 +1034,25 @@ void MacroBNGModel::del_set(const std::vector& rem, void MacroBNGModel::activ_sit(int /*typrul*/, std::string& reac, std::string& prod, std::vector& mreac) { while (!reac.empty()) { - // Match first (...) group - std::regex re_paren("[\\(].*?[)]"); - std::smatch m; - if (!std::regex_search(reac, m, re_paren)) break; + auto p1 = reac.find('('); + if (p1 == std::string::npos) break; + auto p2 = reac.find(')', p1); + if (p2 == std::string::npos) break; - std::string name = m.prefix().str(); // molecule name before '(' - std::string r1 = name + m[0].str(); // e.g. "R(a!1)" - - // Remove the matched portion from reac - reac = m.suffix().str(); - // Strip leading separators ;+. - reac = std::regex_replace(reac, std::regex("^[;+.]+"), ""); - - // Perl pushes r1 to mreac in both branches of the prod check + std::string r1 = reac.substr(0, p2 + 1); mreac.push_back(r1); + + if (p2 + 1 < reac.size()) { + reac = reac.substr(p2 + 1); + size_t non_sep = reac.find_first_not_of(";+."); + if (non_sep != std::string::npos) { + reac = reac.substr(non_sep); + } else { + reac = ""; + } + } else { + reac = ""; + } } } @@ -1181,7 +1259,7 @@ std::string MacroBNGModel::num_site(const std::string& re, const std::string& pr // Extract contents inside parentheses from reactant std::string name; std::smatch m; - std::regex re_paren("[\\(](.*)[\\)]"); + static const std::regex re_paren("[\\(](.*)[\\)]"); std::string ss_re; if (std::regex_search(re, m, re_paren)) { @@ -1349,13 +1427,15 @@ void MacroBNGModel::hash_sor( // Find matching molecule in product std::string p1; { - std::string qname = quotemeta(name); - std::regex re_prod(qname + "[\\(].*?[)]"); - std::smatch pm; - if (std::regex_search(prod, pm, re_prod)) { - p1 = pm[0].str(); - // Remove the match from prod - prod = pm.prefix().str() + pm.suffix().str(); + // Replaced dynamic regex with string search: qname + "(" + std::string target = name + "("; + size_t pos = prod.find(target); + if (pos != std::string::npos) { + size_t end_pos = prod.find(')', pos); + if (end_pos != std::string::npos) { + p1 = prod.substr(pos, end_pos - pos + 1); + prod.erase(pos, end_pos - pos + 1); + } } } mprod.push_back(p1); @@ -2442,7 +2522,7 @@ void MacroBNGModel::cor_net(const std::string& param_prefix) { if (obs.find(";" + egf) != std::string::npos || endsWith(obs, ";" + egf)) { // Extract group name: Molecules;name;... - std::regex mol_re("Molecules;(.*?);"); + static const std::regex mol_re("Molecules;(.*?);"); std::smatch m; if (std::regex_search(obs, m, mol_re)) { // egf_tot_[group_name] — Perl assigns @rabm here @@ -2885,10 +2965,27 @@ void MacroBNGModel::delsites( for (size_t j = 0; j < rr1out.size(); ++j) { if (rr1out[j].empty()) continue; // Create sort key: replace bond labels with #, - std::string sortKey = rr1out[j]; - // Replace !xxx, and !xxx) patterns with #, - std::regex bondPat("![^,)]+([,)])"); - sortKey = std::regex_replace(sortKey, bondPat, "#,"); + std::string sortKey; + sortKey.reserve(rr1out[j].size()); + for (size_t i = 0; i < rr1out[j].size(); ) { + if (rr1out[j][i] == '!') { + // look ahead for ',' or ')' + size_t end = i + 1; + while (end < rr1out[j].size() && rr1out[j][end] != ',' && rr1out[j][end] != ')') { + end++; + } + if (end < rr1out[j].size() && end > i + 1) { + sortKey += "#,"; + i = end + 1; + } else { + sortKey += rr1out[j][i]; + i++; + } + } else { + sortKey += rr1out[j][i]; + i++; + } + } rr1s[sortKey] = j; } @@ -2923,14 +3020,15 @@ void MacroBNGModel::delsites( // Renumber bond labels sequentially int bondNum = 0; - while (ou1.find('!') != std::string::npos) { + size_t pos = 0; + while ((pos = ou1.find('!')) != std::string::npos) { bondNum++; - // Find first bond label: !