From 03f83534a200b52c5d1130e5dd38f286855981be Mon Sep 17 00:00:00 2001 From: ASinanSaglam Date: Fri, 8 Apr 2022 11:07:38 -0400 Subject: [PATCH 1/6] first working verion of writeCppfile --- .github/workflows/validation.yml | 2 +- bng2/Perl2/BNGOutput.pm | 581 +++++++++++++++++++++++++++++++ 2 files changed, 582 insertions(+), 1 deletion(-) diff --git a/.github/workflows/validation.yml b/.github/workflows/validation.yml index 4f7b8db9..6cc98801 100644 --- a/.github/workflows/validation.yml +++ b/.github/workflows/validation.yml @@ -7,7 +7,7 @@ on: push: branches: - master - # - gh_actions_ci + - writeCppfile jobs: build_and_validate: diff --git a/bng2/Perl2/BNGOutput.pm b/bng2/Perl2/BNGOutput.pm index e892f6d4..7b07c353 100644 --- a/bng2/Perl2/BNGOutput.pm +++ b/bng2/Perl2/BNGOutput.pm @@ -2669,6 +2669,587 @@ EOF ### +sub writeCppfile +{ + my $model = shift; + my $params = (@_) ? shift : {}; + + # a place to hold errors + my $err; + + # nothing to do if NO_EXEC is true + return '' if $BNGModel::NO_EXEC; + + # nothing to do if there are no reactions + if ( @{$model->RxnList->Array}==0 ) + { + return "writeMexfile() has nothing to do: no reactions in current model. " + ."Did you remember to call generate_network() before attempting to " + ."write network output?"; + } + + # get reference to parameter list + my $plist = $model->ParamList; + + # get model name + my $model_name = $model->Name; + + # Strip prefixed path + my $prefix = defined $params->{prefix} ? $model->getOutputPrefix( $params->{prefix} ) : $model->getOutputPrefix(); + my $suffix = ( defined $params->{suffix} ) ? $params->{suffix} : undef; + if ( $suffix ) + { $prefix .= "_${suffix}"; } + + # split prefix into volume, path and filebase + my ($vol, $path, $filebase) = File::Spec->splitpath($prefix); + + # define mexfile name + my $cpp_filebase = "${filebase}_cvode"; + my $cpp_filename = "${cpp_filebase}.h"; + my $cpp_path = File::Spec->catpath($vol,$path,$cpp_filename); + + # configure options + my $cvode_abstol = 1e-6; + if ( exists $params->{'atol'} ) + { $cvode_abstol = $params->{'atol'}; } + + my $cvode_reltol = 1e-8; + if ( exists $params->{'rtol'} ) + { $cvode_reltol = $params->{'rtol'}; } + + my $cvode_max_num_steps = 2000; + if ( exists $params->{'max_num_steps'} ) + { $cvode_max_num_steps = $params->{'max_num_steps'}; } + + my $cvode_max_err_test_fails = 7; + if ( exists $params->{'max_err_test_fails'} ) + { $cvode_max_err_test_fails = $params->{'max_err_test_fails'}; } + + my $cvode_max_conv_fails = 10; + if ( exists $params->{'max_conv_fails'} ) + { $cvode_max_conv_fails = $params->{'max_conv_fails'}; } + + my $cvode_max_step = '0.0'; + if ( exists $params->{'max_step'} ) + { $cvode_max_step = $params->{'max_step'}; } + + # Stiff = CV_BDF,CV_NEWTON (Default); Non-stiff = CV_ADAMS,CV_FUNCTIONAL + my $cvode_linear_multistep = 'CV_BDF'; + my $cvode_nonlinear_solver = 'CV_NEWTON'; + if ( exists $params->{'stiff'} ) + { + # if stiff is FALSE, then change to CV_ADAMS and CV_FUNCTIONAL + unless ( $params->{'stiff'} ) + { + $cvode_linear_multistep = 'CV_ADAMS'; + $cvode_nonlinear_solver = 'CV_FUNCTIONAL'; + } + } + + # set sparse option (only permitted with CV_NEWTON) + my $cvode_linear_solver; + if ( ($cvode_nonlinear_solver eq 'CV_NEWTON') and ($params->{'sparse'}) ) + { + $cvode_linear_solver = "flag = CVSpgmr(cvode_mem, PREC_NONE, 0);\n" + ." if (check_flag(&flag, \"CVSpgmr\", 1))"; + } + else + { + $cvode_linear_solver = "flag = CVDense(cvode_mem, __N_SPECIES__);\n" + ." if (check_flag(&flag, \"CVDense\", 1))"; + } + + # time options for mscript + my $t_start = 0; + if ( exists $params->{'t_start'} ) + { $t_start = $params->{'t_start'}; } + + my $t_end = 10; + if ( exists $params->{'t_end'} ) + { $t_end = $params->{'t_end'}; } + + my $n_steps = 20; + if ( exists $params->{'n_steps'} ) + { $n_steps = $params->{'n_steps'}; } + + # code snippet for cleaning up dynamic memory before exiting CVODE-MEX + my $cvode_cleanup_memory = "{ \n" + ." N_VDestroy_Serial(expressions);\n" + ." N_VDestroy_Serial(observables);\n" + ." N_VDestroy_Serial(ratelaws); \n" + ." N_VDestroy_Serial(species); \n" + ." CVodeFree(&cvode_mem); \n" + ." return_status = 1; \n" + ." return Result {return_status, species_out, observables_out}; \n" + ." } "; + + # Index parameters associated with Constants, ConstantExpressions and Observables + ($err) = $plist->indexParams(); + if ($err) { return $err }; + + # and retrieve a string of expression definitions + my $n_parameters = $plist->countType( 'Constant' ); + my $n_expressions = $plist->countType( 'ConstantExpression' ) + $n_parameters; + (my $calc_expressions_string, $err) = $plist->getCVodeExpressionDefs(); + if ($err) { return $err }; + + # generate CVode references for species + # (Do this now, because we need references to CVode species for Observable definitions and Rxn Rates) + my $n_species = scalar @{$model->SpeciesList->Array}; + + + # retrieve a string of observable definitions + my $n_observables = scalar @{$model->Observables}; + my $calc_observables_string; + ($calc_observables_string, $err) = $plist->getCVodeObservableDefs(); + if ($err) { return $err }; + + # Construct user-defined functions + my $user_fcn_declarations = ''; + my $user_fcn_definitions = ''; + foreach my $param ( @{ $model->ParamList->Array } ) + { + if ( $param->Type eq 'Function' ) + { + # get reference to the actual Function + my $fcn = $param->Ref; + + # don't write function if it depends on a local observable evaluation (this is useless + # since CVode can't do local evaluations) + next if ( $fcn->checkLocalDependency($plist) ); + + # get function declaration, add it to the user_fcn_declarations string + $user_fcn_declarations .= $fcn->toCVodeString( $plist, {fcn_mode=>'declare',indent=>''} ); + + # get function definition + my $fcn_defn = $fcn->toCVodeString( $plist, {fcn_mode=>'define', indent=>''} ); + + # add definition to the user_fcn_definitions string + $user_fcn_definitions .= $fcn_defn . "\n"; + } + } + + # index reactions + ($err) = $model->RxnList->updateIndex( $plist ); + if ($err) { return $err }; + + # retrieve a string of reaction rate definitions + my $n_reactions = scalar @{$model->RxnList->Array}; + my $calc_ratelaws_string; + ($calc_ratelaws_string, $err) = $model->RxnList->getCVodeRateDefs( $plist ); + if ($err) { return $err }; + + + # get stoichiometry matrix (sparse encoding in a hashmap) + my $stoich_hash = {}; + ($err) = $model->RxnList->calcStoichMatrix( $stoich_hash ); + + # retrieve a string of species deriv definitions + my $calc_derivs_string; + ($calc_derivs_string, $err) = $model->SpeciesList->toCVodeString( $model->RxnList, $stoich_hash, $plist ); + if ($err) { return $err }; + + # open Mexfile and begin printing... + open( Cppfile, ">$cpp_path" ) or die "Couldn't open $cpp_path: $!\n"; + print Cppfile <<"EOF"; +/* +** ${cpp_filename} +** +** Cvode-C++ implementation of BioNetGen model '$model_name'. +** +** Code Adapted from templates provided by Mathworks and Sundials. +** +** Requires the CVODE libraries: sundials_cvode and sundials_nvecserial. +** https://computation.llnl.gov/casc/sundials/main.html +** +**----------------------------------------------------------------------------- +** +** Compilation notes: +** +** include the model in your C++ file with +** +** #include <$model_name.h> +** +** and compile with +** +** g++ -I/path/to/$model_name.h -c your_script.cpp +** g++ your_script.o -L/path/to/sundials/lib -lsundials_cvode -lsundials_nvecserial -o run_model +** ./run_model +** +** note1: if cvode is in your library path, you can omit path specifications. +** +**----------------------------------------------------------------------------- +** +** Usage in C++ : +** +** Result result; +** vector time_points; +** vector species_initial_values; +** vector parameters; +** result = simulate(timepoints, species_initial_values, parameters); +** +** result.status: is the simulation status +** result.species: is the time series for each species, result.species[N][t] +** is the value of species N at time t +** result.observables: is the time series for each observable, result.observables[N][t] +** is the value of observable N at time t +*/ + +/* Library headers */ +#include +#include +#include +#include + +#include /* prototypes for CVODE */ +#include /* serial N_Vector */ +#include /* prototype for CVDense */ +#include /* prototype for CVSpgmr */ + +using namespace std; + +/* Problem Dimensions */ +#define __N_PARAMETERS__ $n_parameters +#define __N_EXPRESSIONS__ $n_expressions +#define __N_OBSERVABLES__ $n_observables +#define __N_RATELAWS__ $n_reactions +#define __N_SPECIES__ $n_species + +/* user-defined function declarations */ +$user_fcn_declarations + +/* user-defined function definitions */ +$user_fcn_definitions + +/* Calculate expressions */ +void +calc_expressions ( N_Vector expressions, double * parameters ) +{ +$calc_expressions_string +} + +/* Calculate observables */ +void +calc_observables ( N_Vector observables, N_Vector species, N_Vector expressions ) +{ +$calc_observables_string +} + +/* Calculate ratelaws */ +void +calc_ratelaws ( N_Vector ratelaws, N_Vector species, N_Vector expressions, N_Vector observables ) +{ +$calc_ratelaws_string +} + + +/* Calculate species derivatives */ +int +calc_species_deriv ( realtype time, N_Vector species, N_Vector Dspecies, void * f_data ) +{ + int return_val; + N_Vector * temp_data; + + N_Vector expressions; + N_Vector observables; + N_Vector ratelaws; + + /* cast temp_data */ + temp_data = (N_Vector*)f_data; + + /* sget ratelaws Vector */ + expressions = temp_data[0]; + observables = temp_data[1]; + ratelaws = temp_data[2]; + + /* calculate observables */ + calc_observables( observables, species, expressions ); + + /* calculate ratelaws */ + calc_ratelaws( ratelaws, species, expressions, observables ); + + /* calculate derivatives */ +$calc_derivs_string + + return(0); +} + +/* Check function return value... + * opt == 0 means SUNDIALS function allocates memory so check if + * returned NULL pointer + * opt == 1 means SUNDIALS function returns a flag so check if + * flag >= 0 + * opt == 2 means function allocates memory so check if returned + * NULL pointer + */ +int check_flag(void *flagvalue, char *funcname, int opt) +{ + int *errflag; + + /* Check if SUNDIALS function returned NULL pointer - no memory allocated */ + if (opt == 0 && flagvalue == NULL) + { + cout << "\\nSUNDIALS_ERROR: " << funcname << "() failed - returned NULL pointer" << endl; + return(1); + } + + /* Check if flag < 0 */ + else if (opt == 1) + { + errflag = (int *) flagvalue; + if (*errflag < 0) + { + cout << "\\nSUNDIALS_ERROR: " << funcname << "() failed with flag = " << errflag << endl; + return(1); + } + } + + /* Check if function returned NULL pointer - no memory allocated */ + else if (opt == 2 && flagvalue == NULL) + { + cout << "\\nMEMORY_ERROR: " << funcname << "() failed - returned NULL pointer" << endl; + return(1); + } + + return(0); +} + +/* +** ======== +** result structure +** ======== +*/ + +struct Result { + int status; + vector> species; + vector> observables; +}; + +/* +** ======== +** main simulate command +** ======== +*/ +Result simulate( vector timepoints, vector species_init, vector parameters) +{ + int n_timepoints; + size_t i; + size_t j; + + /* intermediate data vectors */ + N_Vector expressions; + N_Vector observables; + N_Vector ratelaws; + + /* array to hold pointers to data vectors */ + N_Vector temp_data[3]; + + /* CVODE specific variables */ + realtype reltol; + realtype abstol; + realtype time; + N_Vector species; + void * cvode_mem; + int flag; + + /* make sure timepoints has correct dimensions */ + if ( (timepoints.size() < 2) ) + { cout << "TIMEPOINTS must be a vector with 2 or more elements." <> species_out( __N_SPECIES__ , vector (n_timepoints, 0)); + vector> observables_out( __N_OBSERVABLES__ , vector (n_timepoints, 0)); + + double * ptr_parameters = ¶meters[0]; + double * ptr_timepoints = &timepoints[0]; + double * ptr_species_init = &species_init[0]; + + /* initialize intermediate data vectors */ + expressions = NULL; + expressions = N_VNew_Serial(__N_EXPRESSIONS__); + if (check_flag((void *)expressions, "N_VNew_Serial", 0)) + { + return_status = 1; + return Result {return_status, species_out, observables_out}; + } + + observables = NULL; + observables = N_VNew_Serial(__N_OBSERVABLES__); + if (check_flag((void *)observables, "N_VNew_Serial", 0)) + { + N_VDestroy_Serial(expressions); + return_status = 1; + return Result {return_status, species_out, observables_out}; + } + + ratelaws = NULL; + ratelaws = N_VNew_Serial(__N_RATELAWS__); + if (check_flag((void *)ratelaws, "N_VNew_Serial", 0)) + { + N_VDestroy_Serial(expressions); + N_VDestroy_Serial(observables); + return_status = 1; + return Result {return_status, species_out, observables_out}; + } + + /* set up pointers to intermediate data vectors */ + temp_data[0] = expressions; + temp_data[1] = observables; + temp_data[2] = ratelaws; + + /* calculate expressions (expressions are constant, so only do this once!) */ + calc_expressions( expressions, ptr_parameters ); + + /* SOLVE model equations! */ + species = NULL; + cvode_mem = NULL; + + /* Set the scalar relative tolerance */ + reltol = $cvode_reltol; + abstol = $cvode_abstol; + + /* Create serial vector for Species */ + species = N_VNew_Serial(__N_SPECIES__); + if (check_flag((void *)species, "N_VNew_Serial", 0)) + { + N_VDestroy_Serial(expressions); + N_VDestroy_Serial(observables); + N_VDestroy_Serial(ratelaws); + return_status = 1; + return Result {return_status, species_out, observables_out}; + } + for ( i = 0; i < __N_SPECIES__; i++ ) + { NV_Ith_S(species,i) = species_init[i]; } + + /* write initial species populations into species_out */ + for ( i = 0; i < __N_SPECIES__; i++ ) + { species_out[i][0] = species_init[i]; } + + /* write initial observables populations into species_out */ + calc_observables( observables, species, expressions ); + for ( i = 0; i < __N_OBSERVABLES__; i++ ) + { observables_out[i][0] = NV_Ith_S(observables,i); } + + /* Call CVodeCreate to create the solver memory: + * CV_ADAMS or CV_BDF is the linear multistep method + * CV_FUNCTIONAL or CV_NEWTON is the nonlinear solver iteration + * A pointer to the integrator problem memory is returned and stored in cvode_mem. + */ + cvode_mem = CVodeCreate($cvode_linear_multistep, $cvode_nonlinear_solver); + if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) + $cvode_cleanup_memory + + + + /* Call CVodeInit to initialize the integrator memory: + * cvode_mem is the pointer to the integrator memory returned by CVodeCreate + * rhs_func is the user's right hand side function in y'=f(t,y) + * T0 is the initial time + * y is the initial dependent variable vector + */ + flag = CVodeInit(cvode_mem, calc_species_deriv, timepoints[0], species); + if (check_flag(&flag, "CVodeInit", 1)) + $cvode_cleanup_memory + + /* Set scalar relative and absolute tolerances */ + flag = CVodeSStolerances(cvode_mem, reltol, abstol); + if (check_flag(&flag, "CVodeSStolerances", 1)) + $cvode_cleanup_memory + + /* pass params to rhs_func */ + flag = CVodeSetUserData(cvode_mem, &temp_data); + if (check_flag(&flag, "CVodeSetFdata", 1)) + $cvode_cleanup_memory + + /* select linear solver */ + $cvode_linear_solver + $cvode_cleanup_memory + + flag = CVodeSetMaxNumSteps(cvode_mem, $cvode_max_num_steps); + if (check_flag(&flag, "CVodeSetMaxNumSteps", 1)) + $cvode_cleanup_memory + + flag = CVodeSetMaxErrTestFails(cvode_mem, $cvode_max_err_test_fails); + if (check_flag(&flag, "CVodeSetMaxErrTestFails", 1)) + $cvode_cleanup_memory + + flag = CVodeSetMaxConvFails(cvode_mem, $cvode_max_conv_fails); + if (check_flag(&flag, "CVodeSetMaxConvFails", 1)) + $cvode_cleanup_memory + + flag = CVodeSetMaxStep(cvode_mem, $cvode_max_step); + if (check_flag(&flag, "CVodeSetMaxStep", 1)) + $cvode_cleanup_memory + + /* integrate to each timepoint */ + for ( i=1; i < n_timepoints; i++ ) + { + flag = CVode(cvode_mem, timepoints[i], species, &time, CV_NORMAL); + if (check_flag(&flag, "CVode", 1)) + { + N_VDestroy_Serial(expressions); + N_VDestroy_Serial(observables); + N_VDestroy_Serial(ratelaws); + N_VDestroy_Serial(species); + CVodeFree(&cvode_mem); + return_status = 1; + return Result {return_status, species_out, observables_out} ; + } + + /* copy species output from nvector to matlab array */ + for ( j = 0; j < __N_SPECIES__; j++ ) + { species_out[j][i] = NV_Ith_S(species,j); } + + /* copy observables output from nvector to matlab array */ + calc_observables( observables, species, expressions ); + for ( j = 0; j < __N_OBSERVABLES__; j++ ) + { observables_out[j][i] = NV_Ith_S(observables,j); } + } + + /* Free vectors */ + N_VDestroy_Serial(expressions); + N_VDestroy_Serial(observables); + N_VDestroy_Serial(ratelaws); + N_VDestroy_Serial(species); + + /* Free integrator memory */ + CVodeFree(&cvode_mem); + + /* return status, species and obserables */ + return Result {return_status, species_out, observables_out}; +} +EOF + close(Cppfile); + + + + # open Mexfile and begin printing... + + print "Wrote Cppfile $cpp_path.\n"; + return (); +} + + + +### +### +### + + + sub writeMfile_QueryNames { From a1bcfdd65b99b7e8e3b88cce76def609edd27aa8 Mon Sep 17 00:00:00 2001 From: ASinanSaglam Date: Wed, 20 Apr 2022 13:55:35 -0400 Subject: [PATCH 2/6] functional version of writeCPYfile --- bng2/Perl2/BNGOutput.pm | 623 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 619 insertions(+), 4 deletions(-) diff --git a/bng2/Perl2/BNGOutput.pm b/bng2/Perl2/BNGOutput.pm index 7b07c353..5d336bfb 100644 --- a/bng2/Perl2/BNGOutput.pm +++ b/bng2/Perl2/BNGOutput.pm @@ -3022,8 +3022,8 @@ int check_flag(void *flagvalue, char *funcname, int opt) struct Result { int status; - vector> species; - vector> observables; + vector< vector > species; + vector< vector > observables; }; /* @@ -3070,8 +3070,8 @@ Result simulate( vector timepoints, vector species_init, vector< /* Create output arrays */ int return_status; - vector> species_out( __N_SPECIES__ , vector (n_timepoints, 0)); - vector> observables_out( __N_OBSERVABLES__ , vector (n_timepoints, 0)); + vector< vector > species_out( __N_SPECIES__ , vector (n_timepoints, 0)); + vector< vector > observables_out( __N_OBSERVABLES__ , vector (n_timepoints, 0)); double * ptr_parameters = ¶meters[0]; double * ptr_timepoints = &timepoints[0]; @@ -3243,6 +3243,621 @@ EOF } +### +### +### + + +sub writeCPYfile +{ + my $model = shift; + my $params = (@_) ? shift : {}; + + # a place to hold errors + my $err; + + # nothing to do if NO_EXEC is true + return '' if $BNGModel::NO_EXEC; + + # nothing to do if there are no reactions + if ( @{$model->RxnList->Array}==0 ) + { + return "writeCPYfile() has nothing to do: no reactions in current model. " + ."Did you remember to call generate_network() before attempting to " + ."write network output?"; + } + + # get reference to parameter list + my $plist = $model->ParamList; + + # get model name + my $model_name = $model->Name; + + # Strip prefixed path + my $prefix = defined $params->{prefix} ? $model->getOutputPrefix( $params->{prefix} ) : $model->getOutputPrefix(); + my $suffix = ( defined $params->{suffix} ) ? $params->{suffix} : undef; + if ( $suffix ) + { $prefix .= "_${suffix}"; } + + # split prefix into volume, path and filebase + my ($vol, $path, $filebase) = File::Spec->splitpath($prefix); + + # define c file name + my $cpy_filebase = "${filebase}_cvode_py"; + my $cpy_filename = "${cpy_filebase}.c"; + my $cpy_path = File::Spec->catpath($vol,$path,$cpy_filename); + + # configure options + my $cvode_abstol = 1e-6; + if ( exists $params->{'atol'} ) + { $cvode_abstol = $params->{'atol'}; } + + my $cvode_reltol = 1e-8; + if ( exists $params->{'rtol'} ) + { $cvode_reltol = $params->{'rtol'}; } + + my $cvode_max_num_steps = 2000; + if ( exists $params->{'max_num_steps'} ) + { $cvode_max_num_steps = $params->{'max_num_steps'}; } + + my $cvode_max_err_test_fails = 7; + if ( exists $params->{'max_err_test_fails'} ) + { $cvode_max_err_test_fails = $params->{'max_err_test_fails'}; } + + my $cvode_max_conv_fails = 10; + if ( exists $params->{'max_conv_fails'} ) + { $cvode_max_conv_fails = $params->{'max_conv_fails'}; } + + my $cvode_max_step = '0.0'; + if ( exists $params->{'max_step'} ) + { $cvode_max_step = $params->{'max_step'}; } + + # Stiff = CV_BDF,CV_NEWTON (Default); Non-stiff = CV_ADAMS,CV_FUNCTIONAL + my $cvode_linear_multistep = 'CV_BDF'; + my $cvode_nonlinear_solver = 'CV_NEWTON'; + if ( exists $params->{'stiff'} ) + { + # if stiff is FALSE, then change to CV_ADAMS and CV_FUNCTIONAL + unless ( $params->{'stiff'} ) + { + $cvode_linear_multistep = 'CV_ADAMS'; + $cvode_nonlinear_solver = 'CV_FUNCTIONAL'; + } + } + + # set sparse option (only permitted with CV_NEWTON) + my $cvode_linear_solver; + if ( ($cvode_nonlinear_solver eq 'CV_NEWTON') and ($params->{'sparse'}) ) + { + $cvode_linear_solver = "flag = CVSpgmr(cvode_mem, PREC_NONE, 0);\n" + ." if (check_flag(&flag, \"CVSpgmr\", 1))"; + } + else + { + $cvode_linear_solver = "flag = CVDense(cvode_mem, __N_SPECIES__);\n" + ." if (check_flag(&flag, \"CVDense\", 1))"; + } + + # time options for mscript + my $t_start = 0; + if ( exists $params->{'t_start'} ) + { $t_start = $params->{'t_start'}; } + + my $t_end = 10; + if ( exists $params->{'t_end'} ) + { $t_end = $params->{'t_end'}; } + + my $n_steps = 20; + if ( exists $params->{'n_steps'} ) + { $n_steps = $params->{'n_steps'}; } + + # code snippet for cleaning up dynamic memory before exiting CVODE + my $cvode_cleanup_memory = "{ \n" + ." N_VDestroy_Serial(expressions);\n" + ." N_VDestroy_Serial(observables);\n" + ." N_VDestroy_Serial(ratelaws); \n" + ." N_VDestroy_Serial(species); \n" + ." CVodeFree(&cvode_mem); \n" + ." result->status = 1; \n" + ." return result; \n" + ." } "; + + # Index parameters associated with Constants, ConstantExpressions and Observables + ($err) = $plist->indexParams(); + if ($err) { return $err }; + + # and retrieve a string of expression definitions + my $n_parameters = $plist->countType( 'Constant' ); + my $n_expressions = $plist->countType( 'ConstantExpression' ) + $n_parameters; + (my $calc_expressions_string, $err) = $plist->getCVodeExpressionDefs(); + if ($err) { return $err }; + + # generate CVode references for species + # (Do this now, because we need references to CVode species for Observable definitions and Rxn Rates) + my $n_species = scalar @{$model->SpeciesList->Array}; + + # generate species string + my $species_string = ""; + foreach my $spec ( @{$model->SpeciesList->Array} ) + { + $species_string .= $spec->SpeciesGraph->toString() . ":"; + } + + # retrieve a string of observable definitions + my $n_observables = scalar @{$model->Observables}; + my $calc_observables_string; + ($calc_observables_string, $err) = $plist->getCVodeObservableDefs(); + if ($err) { return $err }; + + # generate observables string + my $observables_string = ""; + foreach my $obs ( @{$model->Observables} ) + { + $observables_string .= $obs->Name . ":"; + } + + # Construct user-defined functions + my $user_fcn_declarations = ''; + my $user_fcn_definitions = ''; + foreach my $param ( @{ $model->ParamList->Array } ) + { + if ( $param->Type eq 'Function' ) + { + # get reference to the actual Function + my $fcn = $param->Ref; + + # don't write function if it depends on a local observable evaluation (this is useless + # since CVode can't do local evaluations) + next if ( $fcn->checkLocalDependency($plist) ); + + # get function declaration, add it to the user_fcn_declarations string + $user_fcn_declarations .= $fcn->toCVodeString( $plist, {fcn_mode=>'declare',indent=>''} ); + + # get function definition + my $fcn_defn = $fcn->toCVodeString( $plist, {fcn_mode=>'define', indent=>''} ); + + # add definition to the user_fcn_definitions string + $user_fcn_definitions .= $fcn_defn . "\n"; + } + } + + # index reactions + ($err) = $model->RxnList->updateIndex( $plist ); + if ($err) { return $err }; + + # retrieve a string of reaction rate definitions + my $n_reactions = scalar @{$model->RxnList->Array}; + my $calc_ratelaws_string; + ($calc_ratelaws_string, $err) = $model->RxnList->getCVodeRateDefs( $plist ); + if ($err) { return $err }; + + + # get stoichiometry matrix (sparse encoding in a hashmap) + my $stoich_hash = {}; + ($err) = $model->RxnList->calcStoichMatrix( $stoich_hash ); + + # retrieve a string of species deriv definitions + my $calc_derivs_string; + ($calc_derivs_string, $err) = $model->SpeciesList->toCVodeString( $model->RxnList, $stoich_hash, $plist ); + if ($err) { return $err }; + + # open C and begin printing... + open( Cpyfile, ">$cpy_path" ) or die "Couldn't open $cpy_path: $!\n"; + print Cpyfile <<"EOF"; +/* +** ${cpy_filename} +** +** Cvode-C for Python Ctypes implementation of BioNetGen model '$model_name'. +** +** Code Adapted from templates provided by Mathworks and Sundials. +** +** Requires the CVODE libraries: sundials_cvode and sundials_nvecserial. +** https://computation.llnl.gov/casc/sundials/main.html +** +**----------------------------------------------------------------------------- +** +** Compilation notes: +** +** include the model in your C++ file with +** +** #include <$model_name.h> +** +** and compile with +** +** gcc -fPIC -I/path/to/cvode_lib/include/ -c '$model_name'.c +** gcc '$model_name'.o --shared -fPIC -L/path/to/cvode_lib/lib/ -lsundials_cvode -lsundials_nvecserial -o '$model_name'.so +** +** note1: if cvode is in your library path, you can omit path specifications. +** +**----------------------------------------------------------------------------- +** +** Usage in Python : +** +** TODO +*/ + +/* Library headers */ +#include +#include +#include +#include +#include /* prototypes for CVODE */ +#include /* serial N_Vector */ +#include /* prototype for CVDense */ +#include /* prototype for CVSpgmr */ + +/* Problem Dimensions */ +#define __N_PARAMETERS__ $n_parameters +#define __N_EXPRESSIONS__ $n_expressions +#define __N_OBSERVABLES__ $n_observables +#define __N_RATELAWS__ $n_reactions +#define __N_SPECIES__ $n_species + +// return struct +typedef struct result { + int status; + int n_observables; + int n_species; + int n_tpts; + int obs_name_len; + int spcs_name_len; + double *observables; + double *species; + char *obs_names; + char *spcs_names; +} RESULT; + + +/* core function declarations */ +RESULT *simulate ( int num_tpts, double *timepts, int num_species_init, double *species_init, int num_parameters, double *parameters ); +int check_flag ( void *flagvalue, char *funcname, int opt ); +void calc_expressions ( N_Vector expressions, double * parameters ); +void calc_observables ( N_Vector observables, N_Vector species, N_Vector expressions ); +void calc_ratelaws ( N_Vector ratelaws, N_Vector species, N_Vector expressions, N_Vector observables ); +int calc_species_deriv ( realtype time, N_Vector species, N_Vector Dspecies, void * f_data ); + +/* user-defined function declarations */ +$user_fcn_declarations + +/* user-defined function definitions */ +$user_fcn_definitions + +/* Calculate expressions */ +void +calc_expressions ( N_Vector expressions, double * parameters ) +{ +$calc_expressions_string +} + +/* Calculate observables */ +void +calc_observables ( N_Vector observables, N_Vector species, N_Vector expressions ) +{ +$calc_observables_string +} + +/* Calculate ratelaws */ +void +calc_ratelaws ( N_Vector ratelaws, N_Vector species, N_Vector expressions, N_Vector observables ) +{ +$calc_ratelaws_string +} + + +/* Calculate species derivatives */ +int +calc_species_deriv ( realtype time, N_Vector species, N_Vector Dspecies, void * f_data ) +{ + int return_val; + N_Vector * temp_data; + + N_Vector expressions; + N_Vector observables; + N_Vector ratelaws; + + /* cast temp_data */ + temp_data = (N_Vector*)f_data; + + /* sget ratelaws Vector */ + expressions = temp_data[0]; + observables = temp_data[1]; + ratelaws = temp_data[2]; + + /* calculate observables */ + calc_observables( observables, species, expressions ); + + /* calculate ratelaws */ + calc_ratelaws( ratelaws, species, expressions, observables ); + + /* calculate derivatives */ +$calc_derivs_string + + return(0); +} + +/* Check function return value... + * opt == 0 means SUNDIALS function allocates memory so check if + * returned NULL pointer + * opt == 1 means SUNDIALS function returns a flag so check if + * flag >= 0 + * opt == 2 means function allocates memory so check if returned + * NULL pointer + */ +int check_flag(void *flagvalue, char *funcname, int opt) +{ + int *errflag; + + /* Check if SUNDIALS function returned NULL pointer - no memory allocated */ + if (opt == 0 && flagvalue == NULL) + { + printf( "\\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\\n", funcname ); + return(1); + } + + /* Check if flag < 0 */ + else if (opt == 1) + { + errflag = (int *) flagvalue; + if (*errflag < 0) + { + printf( "\\nSUNDIALS_ERROR: %s() failed with flag = %d\\n", funcname, *errflag ); + return(1); + } + } + + /* Check if function returned NULL pointer - no memory allocated */ + else if (opt == 2 && flagvalue == NULL) + { + printf( "\\nMEMORY_ERROR: %s() failed - returned NULL pointer\\n", funcname ); + return(1); + } + + return(0); +} + +/* +** ======================== +** main simulate command +** ======================== +*/ +// we have an input array of (a) time points, (b) species init and (c) parameters +// we have an output array of (a) status, (b) species and (c) observables we need to edit +// plhs = inputs, prhs = outputs, nrhs/nrhs are gone +RESULT *simulate( int num_tpts, double *timepoints, int num_species_init, double *species_init, int num_parameters, double *parameters ) +{ + /* variables */ + size_t i; + size_t j; + + /* intermediate data vectors */ + N_Vector expressions; + N_Vector observables; + N_Vector ratelaws; + + /* array to hold pointers to data vectors */ + N_Vector temp_data[3]; + + /* CVODE specific variables */ + realtype reltol; + realtype abstol; + realtype time; + N_Vector species; + void * cvode_mem; + int flag; + + /* make sure timepoints has correct dimensions */ + if ( num_tpts <= 1 ) + { printf("TIMEPOINTS must be a column vector with 2 or more elements."); } + + /* make sure species_init has correct dimensions */ + if ( num_species_init != __N_SPECIES__ ) + { printf("SPECIES_INIT must be a row vector with 4 elements."); } + + /* make sure params has correct dimensions */ + if ( num_parameters != __N_PARAMETERS__ ) + { printf("PARAMS must be a column vector with 4 elements."); } + + // set output result object + RESULT res_obj; + RESULT *result; + result = malloc(sizeof(RESULT)); + res_obj.n_observables = __N_OBSERVABLES__; + res_obj.n_species = __N_SPECIES__; + res_obj.n_tpts = num_tpts; + double *res_species_ptr; + double *res_observables_ptr; + res_species_ptr = malloc(num_tpts*__N_SPECIES__*sizeof(double)); + res_observables_ptr = malloc(num_tpts*__N_OBSERVABLES__*sizeof(double)); + res_obj.species = res_species_ptr; + res_obj.observables = res_observables_ptr; + res_obj.status = 0; + res_obj.obs_name_len = 0; + res_obj.spcs_name_len = 0; + // store the observable and species names + char onames[] = "$observables_string"; + char snames[] = "$species_string"; + char *ptr_obs_names = strdup(onames); + char *ptr_spcs_names = strdup(snames); + size_t olen = sizeof(onames); + size_t slen = sizeof(snames); + res_obj.obs_names = ptr_obs_names; + res_obj.spcs_names = ptr_spcs_names; + res_obj.obs_name_len = olen; + res_obj.spcs_name_len = slen; + + *result = res_obj; + + /* initialize intermediate data vectors */ + expressions = NULL; + expressions = N_VNew_Serial(__N_EXPRESSIONS__); + if (check_flag((void *)expressions, "N_VNew_Serial", 0)) + { + result->status = 1; + return result; + } + + observables = NULL; + observables = N_VNew_Serial(__N_OBSERVABLES__); + if (check_flag((void *)observables, "N_VNew_Serial", 0)) + { + N_VDestroy_Serial(expressions); + result->status = 1; + return result; + } + + ratelaws = NULL; + ratelaws = N_VNew_Serial(__N_RATELAWS__); + if (check_flag((void *)ratelaws, "N_VNew_Serial", 0)) + { + N_VDestroy_Serial(expressions); + N_VDestroy_Serial(observables); + result->status = 1; + return result; + } + + /* set up pointers to intermediate data vectors */ + temp_data[0] = expressions; + temp_data[1] = observables; + temp_data[2] = ratelaws; + + /* calculate expressions (expressions are constant, so only do this once!) */ + calc_expressions( expressions, parameters ); + + /* SOLVE model equations! */ + species = NULL; + cvode_mem = NULL; + + /* Set the scalar relative tolerance */ + reltol = $cvode_reltol; + abstol = $cvode_abstol; + + /* Create serial vector for Species */ + species = N_VNew_Serial(__N_SPECIES__); + if (check_flag((void *)species, "N_VNew_Serial", 0)) + { + N_VDestroy_Serial(expressions); + N_VDestroy_Serial(observables); + N_VDestroy_Serial(ratelaws); + result->status = 1; + return result; + } + for ( i = 0; i < __N_SPECIES__; i++ ) + { NV_Ith_S(species,i) = species_init[i]; } + + /* write initial species populations into species_out */ + for ( i = 0; i < __N_SPECIES__; i++ ) + { res_species_ptr[i*num_tpts] = species_init[i]; } + + /* write initial observables populations into species_out */ + calc_observables( observables, species, expressions ); + for ( i = 0; i < __N_OBSERVABLES__; i++ ) + { res_observables_ptr[i*num_tpts] = NV_Ith_S(observables,i); } + + /* Call CVodeCreate to create the solver memory: + * CV_ADAMS or CV_BDF is the linear multistep method + * CV_FUNCTIONAL or CV_NEWTON is the nonlinear solver iteration + * A pointer to the integrator problem memory is returned and stored in cvode_mem. + */ + cvode_mem = CVodeCreate($cvode_linear_multistep, $cvode_nonlinear_solver); + if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) + $cvode_cleanup_memory + + + + /* Call CVodeInit to initialize the integrator memory: + * cvode_mem is the pointer to the integrator memory returned by CVodeCreate + * rhs_func is the user's right hand side function in y'=f(t,y) + * T0 is the initial time + * y is the initial dependent variable vector + */ + flag = CVodeInit(cvode_mem, calc_species_deriv, timepoints[0], species); + if (check_flag(&flag, "CVodeInit", 1)) + $cvode_cleanup_memory + + /* Set scalar relative and absolute tolerances */ + flag = CVodeSStolerances(cvode_mem, reltol, abstol); + if (check_flag(&flag, "CVodeSStolerances", 1)) + $cvode_cleanup_memory + + /* pass params to rhs_func */ + flag = CVodeSetUserData(cvode_mem, &temp_data); + if (check_flag(&flag, "CVodeSetFdata", 1)) + $cvode_cleanup_memory + + /* select linear solver */ + $cvode_linear_solver + $cvode_cleanup_memory + + flag = CVodeSetMaxNumSteps(cvode_mem, $cvode_max_num_steps); + if (check_flag(&flag, "CVodeSetMaxNumSteps", 1)) + $cvode_cleanup_memory + + flag = CVodeSetMaxErrTestFails(cvode_mem, $cvode_max_err_test_fails); + if (check_flag(&flag, "CVodeSetMaxErrTestFails", 1)) + $cvode_cleanup_memory + + flag = CVodeSetMaxConvFails(cvode_mem, $cvode_max_conv_fails); + if (check_flag(&flag, "CVodeSetMaxConvFails", 1)) + $cvode_cleanup_memory + + flag = CVodeSetMaxStep(cvode_mem, $cvode_max_step); + if (check_flag(&flag, "CVodeSetMaxStep", 1)) + $cvode_cleanup_memory + + /* integrate to each timepoint */ + for ( i=1; i < num_tpts; i++ ) + { + flag = CVode(cvode_mem, timepoints[i], species, &time, CV_NORMAL); + if (check_flag(&flag, "CVode", 1)) + { + N_VDestroy_Serial(expressions); + N_VDestroy_Serial(observables); + N_VDestroy_Serial(ratelaws); + N_VDestroy_Serial(species); + CVodeFree(&cvode_mem); + result->status = 1; + return result; + } + + /* copy species output from nvector to matlab array */ + for ( j = 0; j < __N_SPECIES__; j++ ) + { res_species_ptr[j*num_tpts + i] = NV_Ith_S(species,j); } + + /* copy observables output from nvector to matlab array */ + calc_observables( observables, species, expressions ); + for ( j = 0; j < __N_OBSERVABLES__; j++ ) + { res_observables_ptr[j*num_tpts + i] = NV_Ith_S(observables,j); } + } + + /* Free vectors */ + N_VDestroy_Serial(expressions); + N_VDestroy_Serial(observables); + N_VDestroy_Serial(ratelaws); + N_VDestroy_Serial(species); + + /* Free integrator memory */ + CVodeFree(&cvode_mem); + + return result; +} + +void free_result(RESULT *r) { + free(r->obs_names); + free(r->spcs_names); + free(r->observables); + free(r->species); + free(r); +} +EOF + close(Cpyfile); + + + + # open Mexfile and begin printing... + + print "Wrote Cpyfile $cpy_path.\n"; + return (); +} ### ### From 1efe521993af801e883230d47af3f6469915e451 Mon Sep 17 00:00:00 2001 From: ASinanSaglam Date: Mon, 25 Apr 2022 16:05:30 -0400 Subject: [PATCH 3/6] fixed CPY max parameters/species counts being hardcoded --- bng2/Perl2/BNGOutput.pm | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/bng2/Perl2/BNGOutput.pm b/bng2/Perl2/BNGOutput.pm index 5d336bfb..5c1a884a 100644 --- a/bng2/Perl2/BNGOutput.pm +++ b/bng2/Perl2/BNGOutput.pm @@ -3651,11 +3651,11 @@ RESULT *simulate( int num_tpts, double *timepoints, int num_species_init, double /* make sure species_init has correct dimensions */ if ( num_species_init != __N_SPECIES__ ) - { printf("SPECIES_INIT must be a row vector with 4 elements."); } + { printf("SPECIES_INIT must be a row vector with $n_species elements."); } /* make sure params has correct dimensions */ if ( num_parameters != __N_PARAMETERS__ ) - { printf("PARAMS must be a column vector with 4 elements."); } + { printf("PARAMS must be a column vector with $n_parameters elements."); } // set output result object RESULT res_obj; From 33938a7d402a8f1ffa9430ddf687a5de8772dca1 Mon Sep 17 00:00:00 2001 From: ASinanSaglam Date: Tue, 26 Apr 2022 14:09:35 -0400 Subject: [PATCH 4/6] fixed max param/species in CPP file --- bng2/Perl2/BNGOutput.pm | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/bng2/Perl2/BNGOutput.pm b/bng2/Perl2/BNGOutput.pm index 5c1a884a..73d1186d 100644 --- a/bng2/Perl2/BNGOutput.pm +++ b/bng2/Perl2/BNGOutput.pm @@ -3059,11 +3059,11 @@ Result simulate( vector timepoints, vector species_init, vector< /* make sure species_init has correct dimensions */ if ( (species_init.size() != __N_SPECIES__) ) - { cout << "SPECIES_INIT must be a vector with 4 elements." << endl; } + { cout << "SPECIES_INIT must be a vector with $n_species elements." << endl; } /* make sure params has correct dimensions */ if ( (parameters.size() != __N_PARAMETERS__) ) - { cout << "PARAMS must be a column vector with 4 elements." << endl; } + { cout << "PARAMS must be a column vector with $n_parameters elements." << endl; } /* get number of timepoints */ n_timepoints = timepoints.size(); From 6388b3292eca1785f830faa67ebf2004269716b2 Mon Sep 17 00:00:00 2001 From: ASinanSaglam Date: Tue, 26 Apr 2022 14:28:51 -0400 Subject: [PATCH 5/6] removing branch from GH actions --- .github/workflows/validation.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/validation.yml b/.github/workflows/validation.yml index 6cc98801..25bc52fe 100644 --- a/.github/workflows/validation.yml +++ b/.github/workflows/validation.yml @@ -7,7 +7,6 @@ on: push: branches: - master - - writeCppfile jobs: build_and_validate: From 68cd56ae974d6f6f1dde19696e40d534457aa476 Mon Sep 17 00:00:00 2001 From: ASinanSaglam Date: Tue, 26 Apr 2022 14:35:08 -0400 Subject: [PATCH 6/6] changed cpp command to be all capital CPP --- bng2/Perl2/BNGOutput.pm | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bng2/Perl2/BNGOutput.pm b/bng2/Perl2/BNGOutput.pm index 73d1186d..1a2fdca4 100644 --- a/bng2/Perl2/BNGOutput.pm +++ b/bng2/Perl2/BNGOutput.pm @@ -2669,7 +2669,7 @@ EOF ### -sub writeCppfile +sub writeCPPfile { my $model = shift; my $params = (@_) ? shift : {};