Closes DiODeProject#170

mumot/models.py

@@ -22,6 +22,7 @@
     solve,
     Symbol,
     symbols,
+    Function
 )
 
 from . import (
@@ -468,7 +469,8 @@ def showMasterEquation(self):
 
         """
 
-        P, t = symbols('P t')
+        t = symbols('t')
+        P = Function('P')
         out_rhs = ""
         stoich = self._stoichiometry
         nvec = []
@@ -2208,7 +2210,8 @@ def _get_orderedLists_vKE(stoich):
 
 def _getFokkerPlanckEquation(_get_orderedLists_vKE, stoich):
     """Return the Fokker-Planck equation."""
-    P, t = symbols('P t')
+    t = symbols('t')
+    P = Function('P')
     V = Symbol(r'\overline{V}', real=True, constant=True)
     Vlist_lhs, Vlist_rhs, substring = _get_orderedLists_vKE(stoich)
     rhsFPE = 0
@@ -2261,7 +2264,10 @@ def _getNoiseEOM(_getFokkerPlanckEquation, _get_orderedLists_vKE, stoich):
     Returns equations of motion for noise.
 
     """
-    P, M_1, M_2, t = symbols('P M_1 M_2 t')
+    t = symbols('t')
+    P = Function('P')
+    M_1 = Function('M_1')
+    M_2 = Function('M_2')
 
     # A,B, alpha, beta, gamma = symbols('A B alpha beta gamma')
     # custom_stoich= {'reaction1': {'rate': alpha, A: [0,1]}, 'reaction2': {'rate': gamma, A: [2,0], B: [0,1]},
@@ -2440,7 +2446,10 @@ def _getNoiseStationarySol(_getNoiseEOM, _getFokkerPlanckEquation, _get_orderedL
     Returns analytical solution for stationary noise.
 
     """
-    P, M_1, M_2, t = symbols('P M_1 M_2 t')
+    t = symbols('t')
+    P = Function('P')
+    M_1 = Function('M_1')
+    M_2 = Function('M_2')
 
     EQsys1stOrdMom, EOM_1stOrderMom, NoiseSubs1stOrder, EQsys2ndOrdMom, EOM_2ndOrderMom, NoiseSubs2ndOrder = _getNoiseEOM(_getFokkerPlanckEquation, _get_orderedLists_vKE, stoich)
 
@@ -2562,7 +2571,8 @@ def _getNoiseStationarySol(_getNoiseEOM, _getFokkerPlanckEquation, _get_orderedL
 
 def _getODEs_vKE(_get_orderedLists_vKE, stoich):
     """Return the ODE system derived from Master equation."""
-    P, t = symbols('P t')
+    t = symbols('t')
+    P = Function('P')
     V = Symbol(r'\overline{V}', real=True, constant=True)
     Vlist_lhs, Vlist_rhs, substring = _get_orderedLists_vKE(stoich)
     rhsODE = 0

mumot/views.py

@@ -1189,7 +1189,8 @@ def _plot_NumSolODE(self, _=None):
             for key in EOM_1stOrdMomDict:
                 noiseCorrEOMdict[sym * key] = sympy.expand(sym * EOM_1stOrdMomDict[key])
 
-        M_1, M_2 = sympy.symbols('M_1 M_2')
+        M_1 = Function('M_1')
+        M_2 = Function('M_2')
         eta_SV1 = sympy.Symbol(f"eta_{self._stateVariable1}")
         cVar1 = sympy.symbols('cVar1')
         if self._stateVariable2:
@@ -1333,7 +1334,8 @@ def _numericSol2ndOrdMoment(self, EOM_2ndOrdMomDict, steadyStateDict, argDict):
         # if len(self._stateVarList) == 3:
         #     eta_SV3 = sympy.Symbol('eta_' + str(self._stateVarList[2]))
 
-        M_1, M_2 = sympy.symbols('M_1 M_2')
+        M_1 = Function('M_1')
+        M_2 = Function('M_2')
         eta_SV1 = sympy.Symbol('eta_' + str(self._stateVariable1))
         if self._stateVariable2:
             eta_SV2 = sympy.Symbol('eta_' + str(self._stateVariable2))
@@ -5644,8 +5646,10 @@ def _deriveMasterEquation(stoichiometry):
     """
     substring = None
 
-    P, E_op, x, y, v, w, t, m = symbols('P E_op x y v w t m')
-    x, y, v, w = symbols('x y v w', cls=Function)
+    x, y, v, w, t, m = symbols('x y v w t m')
+    E_op = Function('E_op')
+    z = Function('z')
+    P = Function('P')
     V = Symbol(r'\overline{V}', real=True, constant=True)
 
     stoich = stoichiometry
@@ -5668,7 +5672,7 @@ def _deriveMasterEquation(stoichiometry):
 
     rhs = 0
     sol_dict_rhs = {}
-    f = lambdify((x(y, v - w)), x(y, v - w), modules='sympy')
+    f = lambdify(z(y, v - w), z(y, v - w), modules='sympy')
     g = lambdify((x, y, v), (factorial(x) / factorial(x - y)) / v**y, modules='sympy')
     for key1 in stoich:
         prod1 = 1
@@ -5695,7 +5699,9 @@ def _deriveMasterEquation(stoichiometry):
 
 def _doVanKampenExpansion(rhs, stoich):
     """Return the left-hand side and right-hand side of van Kampen expansion."""
-    P, E_op, x, y, v, w, t, m = symbols('P E_op x y v w t m')
+    x, y, v, w, t, m = symbols('x y v w t m')
+    E_op = Function('E_op')
+    P = Function('P')
     V = Symbol(r'\overline{V}', real=True, constant=True)
     nvec = []
     nconstvec = []

setup.py

@@ -39,7 +39,7 @@
         'pydstool>=0.90.3',  # min version that allows scipy >= 1.0.0 to be used
         'pyzmq<17',  # needed if using tornado < 5
         'scipy',
-        'sympy',
+        'sympy>=1.4',
         'tornado<5'  # needed to avoid errors with older ipykernel
     ],
     extras_require={