Repository: StochSS/SpatialPy (language: python)

Issue:
Problematic results for a simple reversible autocatalytic reaction
This issue was originally reported [here](https://github.com/StochSS/GillesPy2/issues/869)

## Fix
- Remove the `0.5` term when creating homo-bimolecular propensities
- Update unit tests
Fix 325



Discussion / hints:

Fix patch:
diff --git a/spatialpy/__version__.py b/spatialpy/__version__.py
index 3d08e052..11f878a1 100644
--- a/spatialpy/__version__.py
+++ b/spatialpy/__version__.py
@@ -21,7 +21,7 @@
 # @website https://github.com/StochSS/SpatialPy
 # =============================================================================
 
-__version__      = '1.1.0'
+__version__      = '1.1.1'
 __title__        = 'SpatialPy'
 __description__  = 'Python Interface for Spatial Stochastic Biochemical Simulations'
 __url__          = 'https://spatialpy.github.io/SpatialPy/'
diff --git a/spatialpy/core/reaction.py b/spatialpy/core/reaction.py
index cad800c9..e8591195 100644
--- a/spatialpy/core/reaction.py
+++ b/spatialpy/core/reaction.py
@@ -293,7 +293,7 @@ def _create_mass_action(self):
                 reactant = reactant.name
             # Case 1: 2X -> Y
             if stoichiometry == 2:
-                propensity_function = f"0.5 * {propensity_function} * {reactant} * ({reactant} - 1) / vol"
+                propensity_function = f"{propensity_function} * {reactant} * ({reactant} - 1) / vol"
                 ode_propensity_function += f" * {reactant} * {reactant}"
             else:
                 # Case 3: X1, X2 -> Y;
diff --git a/spatialpy/solvers/c_base/ssa_sdpd-c-simulation-engine/src/model.c b/spatialpy/solvers/c_base/ssa_sdpd-c-simulation-engine/src/model.c
index 279bfb4a..346b505e 100644
--- a/spatialpy/solvers/c_base/ssa_sdpd-c-simulation-engine/src/model.c
+++ b/spatialpy/solvers/c_base/ssa_sdpd-c-simulation-engine/src/model.c
@@ -25,7 +25,10 @@ Based on a Matlab program by Bruno Jacob (UCSB)
 This program is distributed under the terms of the GNU GENERAL PUBLIC LICENSE Version 3.
 See the file LICENSE.txt for details.
 ***************************************************************************************** */
-#include <math.h>
+#if defined(WIN32) || defined(_WIN32) || defined(__MINGW32__)
+#define _USE_MATH_DEFINES
+#endif
+#include <cmath>
 #include <stdio.h>
 #include <stdlib.h>
 
diff --git a/spatialpy/solvers/c_base/ssa_sdpd-c-simulation-engine/src/particle.cpp b/spatialpy/solvers/c_base/ssa_sdpd-c-simulation-engine/src/particle.cpp
index c2825347..4efd1b28 100644
--- a/spatialpy/solvers/c_base/ssa_sdpd-c-simulation-engine/src/particle.cpp
+++ b/spatialpy/solvers/c_base/ssa_sdpd-c-simulation-engine/src/particle.cpp
@@ -26,7 +26,10 @@ See the file LICENSE.txt for details.
 
 #include <stdlib.h>
 #include <stdio.h>
-#include <math.h>
+#if defined(WIN32) || defined(_WIN32) || defined(__MINGW32__)
+#define _USE_MATH_DEFINES
+#endif
+#include <cmath>
 #include <cstdlib>
 #include <vector>
 #include <queue>
@@ -84,13 +87,13 @@ namespace Spatialpy{
 
     void Particle::check_particle_nan(){
         if(
-            isnan(x[0]) || !isfinite(x[0]) ||
-            isnan(x[1]) || !isfinite(x[1]) ||
-            isnan(x[2]) || !isfinite(x[2]) ||
-            isnan(v[0]) || !isfinite(v[0]) ||
-            isnan(v[1]) || !isfinite(v[1]) ||
-            isnan(v[2]) || !isfinite(v[2]) ||
-            isnan(rho)  || !isfinite(rho) ){
+            std::isnan(x[0]) || !std::isfinite(x[0]) ||
+            std::isnan(x[1]) || !std::isfinite(x[1]) ||
+            std::isnan(x[2]) || !std::isfinite(x[2]) ||
+            std::isnan(v[0]) || !std::isfinite(v[0]) ||
+            std::isnan(v[1]) || !std::isfinite(v[1]) ||
+            std::isnan(v[2]) || !std::isfinite(v[2]) ||
+            std::isnan(rho)  || !std::isfinite(rho) ){
             printf("ERROR: nan/inf detected!!!\n");
             printf("number of neighbors: %li\n", neighbors.size()) ;
             printf("id=%i\n",id);
@@ -183,7 +186,7 @@ namespace Spatialpy{
     	// Eq 28 of Drawert et al 2019, Tartakovsky et. al., 2007, JCP
     	double D_i_j = -2.0*(mass*neighbor->mass)/(mass+neighbor->mass)*(rho+neighbor->rho)/(rho*neighbor->rho) * r2 * wfd / (r2+0.01*h*h);
 
-    	if(isnan(D_i_j)){
+    	if(std::isnan(D_i_j)){
     	    printf("Got NaN calculating D_i_j for me=%i, neighbor=%i\n",id, neighbor->id);
     	    printf("system->dimension=%i ",system->dimension);
     	    printf("r2=%e ",r2);
diff --git a/spatialpy/solvers/solver.py b/spatialpy/solvers/solver.py
index b1788c3a..0b4cc337 100644
--- a/spatialpy/solvers/solver.py
+++ b/spatialpy/solvers/solver.py
@@ -470,41 +470,50 @@ def compile(self, debug=False, profile=False):
         self.__create_propensity_file(stoich_matrix, dep_graph, file_name=self.prop_file_name)
 
         # Build the solver
+        for make_exe_location in ["make", "mingw32-make"]:
+            make_exe = shutil.which("make")
+            if make_exe is not None:
+                break
+        if make_exe is None:
+            raise SimulationError("Make executable could not be found")
         makefile = self.spatialpy_rootdir+'/build/Makefile'
         makefile_ann = self.spatialpy_rootdir+'/external/ANN/src/Makefile.spatialpy'
-        cmd_list = ['cd', self.core_dir, '&&',
-            'make', '-f' , makefile_ann, 
-            'ROOTINC="' + self.spatialpy_rootinc+'"',
-            '&&',
-            'make', 'CORE', '-f',  makefile,
-            'ROOT="' + self.spatialpy_rootparam+'"',
-            'ROOTINC="' + self.spatialpy_rootinc+'"',
-            'BUILD='+self.core_dir, '&&'
-            'cd', self.build_dir, '&&', 'make', '-I', self.core_dir, '-f', makefile,
-            'ROOT="' + self.spatialpy_rootparam+'"',
-            'ROOTINC="' + self.spatialpy_rootinc+'"',
-            'COREDIR="' + self.core_dir + '"',
-            'MODEL=' + self.prop_file_name, 'BUILD='+self.build_dir]
+
+        cmd_ann = [
+            make_exe, '-d', '-C', self.core_dir, '-f', makefile_ann,
+            f'ROOTINC={self.spatialpy_rootinc}',
+        ]
+        cmd_core = [
+            make_exe, '-d', '-C', self.core_dir, 'CORE', '-f',  makefile,
+            f'ROOT={self.spatialpy_rootparam}',
+            f'ROOTINC={self.spatialpy_rootinc}',
+            f'BUILD={self.core_dir}',
+        ]
+        cmd_build = [
+            make_exe, '-d', '-C', self.build_dir, '-I', self.core_dir, '-f', makefile,
+            'ROOT=' + self.spatialpy_rootparam,
+            'ROOTINC=' + self.spatialpy_rootinc,
+            'COREDIR=' + self.core_dir,
+            'MODEL=' + self.prop_file_name, 'BUILD='+self.build_dir
+        ]
         if profile:
-            cmd_list.append('GPROFFLAG=-pg')
+            cmd_build.append('GPROFFLAG=-pg')
         if profile or debug:
-            cmd_list.append('GDB_FLAG=-g')
-        cmd = " ".join(cmd_list)
+            cmd_build.append('GDB_FLAG=-g')
         if self.debug_level > 1:
+            cmd = " && ".join([*cmd_ann, *cmd_core, *cmd_build])
             print(f"cmd: {cmd}\n")
         try:
-            with subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, shell=True) as handle:
-                stdout, _ = handle.communicate()
-                return_code = handle.wait()
-                if return_code != 0:
-                    try:
-                        print(stdout.decode("utf-8"))
-                    except Exception:
-                        pass
-                    raise SimulationError(f"Compilation of solver failed, return_code={return_code}")
-
+            for cmd_target in [cmd_ann, cmd_core, cmd_build]:
+                result = subprocess.check_output(cmd_target, stderr=subprocess.STDOUT)
                 if self.debug_level > 1:
-                    print(stdout.decode("utf-8"))
+                    print(result.stdout.decode("utf-8"))
+        except subprocess.CalledProcessError as err:
+            try:
+                print(err.stdout.decode("utf-8"))
+            except Exception:
+                pass
+            raise SimulationError(f"Compilation of solver failed, return_code={result.return_code}")
         except OSError as err:
             print(f"Error, execution of compilation raised an exception: {err}")
             print(f"cmd = {cmd}")
@@ -559,7 +568,7 @@ def run(self, number_of_trajectories=1, seed=None, timeout=None,
             result = Result(self.model, outfile)
             if self.debug_level >= 1:
                 print(f"Running simulation. Result dir: {outfile}")
-            solver_cmd = f'cd {outfile};{os.path.join(self.build_dir, self.executable_name)}'
+            solver_cmd = os.path.join(self.build_dir, self.executable_name)
 
             if number_of_threads is not None:
                 solver_cmd += " -t " + str(number_of_threads)
@@ -573,7 +582,7 @@ def run(self, number_of_trajectories=1, seed=None, timeout=None,
             start = time.monotonic()
             return_code = None
             try:
-                with subprocess.Popen(solver_cmd, shell=True,
+                with subprocess.Popen(solver_cmd, cwd=outfile, shell=True,
                         stdout=subprocess.PIPE, stderr=subprocess.STDOUT,
                         start_new_session=True) as process:
                     try:


Test patch:
diff --git a/test/unit_tests/test_reaction.py b/test/unit_tests/test_reaction.py
index d32e51b5..3220fb22 100644
--- a/test/unit_tests/test_reaction.py
+++ b/test/unit_tests/test_reaction.py
@@ -450,7 +450,7 @@ def test__create_mass_action__2X_to_Y(self):
         self.valid_ma_reaction.reactants = {"X": 2}
         self.valid_ma_reaction.products = {"Y": 1}
         self.valid_ma_reaction._create_mass_action()
-        self.assertEqual(self.valid_ma_reaction.propensity_function, "((((0.5*k1)*X)*(X-1))/vol)")
+        self.assertEqual(self.valid_ma_reaction.propensity_function, "(((k1*X)*(X-1))/vol)")
         self.assertEqual(self.valid_ma_reaction.ode_propensity_function, "((k1*X)*X)")
 
     def test__create_mass_action__species_obj_reactant(self):
@@ -581,7 +581,7 @@ def test_add_reactant__massaction_reaction_stoich_0_to_2(self):
         self.valid_ma_reaction.add_reactant("X", 2)
         self.assertIn("X", self.valid_ma_reaction.reactants)
         self.assertEqual(self.valid_ma_reaction.reactants["X"], 2)
-        self.assertEqual(self.valid_ma_reaction.propensity_function, "((((0.5*k1)*X)*(X-1))/vol)")
+        self.assertEqual(self.valid_ma_reaction.propensity_function, "(((k1*X)*(X-1))/vol)")
         self.assertEqual(self.valid_ma_reaction.ode_propensity_function, "((k1*X)*X)")
 
     def test_add_reactant__invalid_species(self):