Merge pull request #22 from ayush9pandey/reduce

New features for model reduction

.github/workflows/autoreduce_ci.yml

@@ -29,6 +29,7 @@ jobs:
         python -m pip install --upgrade pip
         python -m pip install flake8 pytest
         if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
+        pip install .
     - name: Lint with flake8
       run: |
         # stop the build if there are Python syntax errors or undefined names
@@ -37,4 +38,6 @@ jobs:
         flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
     - name: Test with pytest
       run: |
+        pip install pytest-cov
+        pip install nbval
         pytest

.gitignore

@@ -102,4 +102,7 @@ dmypy.json
 .pyre/
 
 *.tex
-*.svg
+*.svg
+
+# Package specific
+examples/reduced_gene_expression.xml

Tests/test_time_scale_separation.py

@@ -6,7 +6,7 @@
 from unittest.mock import mock_open, patch
 from test_auto_reduce import TestAutoReduce
 import warnings
-
+import sympy
 
 class TestTimeScaleSeparation(TestAutoReduce):
     def test_reduced_models(self):
@@ -17,8 +17,14 @@ def test_reduced_models(self):
             attempt_states = [self.x[i] for i in attempt]
             answer_AB = [A**2*B*k1*k2/(k2+k3) - A**2*B*k1,
                          A**2*B*k1*k2/(k2+k3) - A**2*B*k1]
-            self.test_solve_timescale_separation(attempt_states = [A, B], 
-                                                mode = 'fail', answer = answer_AB)
+            # On simplification, the answer is:
+            answer_AB_eq = [-A**2*B*k1*k3/(k2 + k3),-A**2*B*k1*k3/(k2 + k3)] 
+            try:
+                self.test_solve_timescale_separation(attempt_states = [A, B], 
+                                                     mode = 'fail', answer = answer_AB)
+            except:
+                self.test_solve_timescale_separation(attempt_states = [A, B], 
+                                                     mode = 'fail', answer = answer_AB_eq)
             self.test_solve_timescale_separation(attempt_states = [A, D], 
                                                 mode = 'fail', answer = [0, 0])
             self.test_solve_timescale_separation(attempt_states = [A, C], 
@@ -40,9 +46,18 @@ def test_reduced_models(self):
             answer_ABD = [k1 * k2 * A**2 * B / (k2 + k3) - k1 * A**2 * B, 
                         k1 * k2 * A**2 * B / (k2 + k3) - k1 * A**2 * B, 
                         k1 * k3 * A**2 * B / (k2 + k3)]
-            self.test_solve_timescale_separation(attempt_states = [A, B, D], 
-                                                mode = 'success',
-                                                answer = answer_ABD)
+            # On simplification, the answer is:
+            answer_ABD_eq = [-A**2*B*k1*k3/(k2 + k3), 
+                             -A**2*B*k1*k3/(k2 + k3),
+                             A**2*B*k1*k3/(k2 + k3)]
+            try:
+                self.test_solve_timescale_separation(attempt_states = [A, B, D], 
+                                                     mode = 'success',
+                                                     answer = answer_ABD)
+            except:
+                self.test_solve_timescale_separation(attempt_states = [A, B, D], 
+                                                     mode = 'success',
+                                                     answer = answer_ABD_eq)
 
     def test_solve_timescale_separation(self, attempt_states = None, mode = None, answer = None):
         if mode == 'fail':
@@ -54,7 +69,7 @@ def test_solve_timescale_separation(self, attempt_states = None, mode = None, an
         elif mode == 'success':
             reduced_system, collapsed_system = self.test_get_reduced_model(x_hat = attempt_states)
             assert answer is not None
-            self.assertEqual(reduced_system.f, answer)
+            self.assertEqual(reduced_system.f , answer)
             # Test slow (retained) states
             self.assertEqual(reduced_system.x, attempt_states)
             # Test fast (collapsed) states

autoreduce/converters.py

@@ -5,8 +5,7 @@
 from .model_reduction import Reduce
 
 def load_ODE_model(n_states, n_params = 0):
-    x, f, P = ode_to_sympy(n_states, n_params)
-    return System(x, f, P)
+    return ode_to_sympy(n_states, n_params)
 
 def ode_to_sympy(odesize, n_params = 0):
     '''
@@ -70,7 +69,7 @@ def load_sbml(filename, **kwargs):
     if doc.getNumErrors(LIBSBML_SEV_FATAL):
         print('Encountered serious errors while reading file')
         print(doc.getErrorLog().toString())
-        sys.exit(1)
+        return
     doc.getErrorLog().clearLog()
     # Convert local params to global params
     props = ConversionProperties()

autoreduce/local_sensitivity.py

@@ -85,8 +85,13 @@ def compute_J(self, x, **kwargs):
         else:
             var = x
         # initialize J
-        J = np.zeros( (self.n, len(var)) )   
+        J = [[0]*self.n]*len(var)
+        # Change this variable name P to params_values or something
         P = self.params_values 
+        if 'set_params_as' in kwargs:
+            set_params_as = kwargs['set_params_as']
+        else:
+            set_params_as = None
         u = self.u
         if 'mode' in kwargs:
             if kwargs.get('mode') == 'accurate':
@@ -106,30 +111,33 @@ def fun_ode(t, x, params):
         X = var 
         for i in range(self.n):
             for j in range(len(var)):
-                F = np.zeros( (4,1) )
-                h = X[j]*0.01
+                # F = np.zeros( (4,1) )
+                F = [None]*4
+                h = 0.01
                 # Gets O(4) central difference on dfi/dvarj
                 if h != 0:
                     var = X
                     var[j] = X[j] + 2*h
-                    f = self.evaluate(fun, var, P, u)
+                    f = self.evaluate(fun, var, P, u, set_params_as = set_params_as)
                     F[0] = f[i]
                     var[j] = X[j] + h
-                    f = self.evaluate(fun, var, P, u)
+                    f = self.evaluate(fun, var, P, u, set_params_as = set_params_as)
                     F[1] = f[i]
                     var[j] = X[j] - h
-                    f = self.evaluate(fun, var, P, u)
+                    f = self.evaluate(fun, var, P, u, set_params_as = set_params_as)
                     F[2] = f[i]
                     var[j] = X[j] - 2*h
-                    f = self.evaluate(fun, var, P, u)
+                    f = self.evaluate(fun, var, P, u, set_params_as = set_params_as)
                     F[3] = f[i]
                     #Store approvar. dfi/dvarj into J
-                    J[i,j]= (-F[0] + 8*F[1] - 8*F[2] + F[3])/(12*h)   
+                    J[i][j]= (-F[0] + 8*F[1] - 8*F[2] + F[3])/(12*h)   
                     # print(J[i,j])
                     # if J[i,j] == np.Inf:
                     #     J[i,j] = 1
                     # elif J[i,j] == np.NaN:
                     #     J[i,j] = 0
+        if set_params_as is None:
+            J = np.array(J, dtype = float)
         return J
 
     def compute_SSM(self, normalize = False, **kwargs):
@@ -231,20 +239,21 @@ def sensitivity_to_parameter(self, x, j, **kwargs):
                     sensitivity of each variable to the specified parameter over
                     time.
         '''
+        pass
         # Build a scipy-integratable derivative-of-sensitivity function.
-        import numdifftools as nd
-        import copy
-        def dS_dt(t, s):
-            xs = ode_sol(t)
-            # Wrapper to let numdifftools calculate df/dp.
-            def ode_as_parameter_call(param):
-                call_params = copy.deepcopy(self.params_values)
-                call_params[j] = self.params_values
-                return ode(t, xs, call_params)
-            df_dp = lambda xs: nd.Jacobian(ode_as_parameter_call)(xs).transpose()[:,0]
-            return df_dp(params[j]) + np.matmul(ode_jac(xs), s)
-        sol = odeint(dS_dt, np.zeros(n_vars), self.timepoints, **kwargs)
-        return sol
+        # import numdifftools as nd
+        # import copy
+        # def dS_dt(t, s):
+        #     xs = ode_sol(t)
+        #     # Wrapper to let numdifftools calculate df/dp.
+        #     def ode_as_parameter_call(param):
+        #         call_params = copy.deepcopy(self.params_values)
+        #         call_params[j] = self.params_values
+        #         return ode(t, xs, call_params)
+        #     df_dp = lambda xs: nd.Jacobian(ode_as_parameter_call)(xs).transpose()[:,0]
+        #     return df_dp(params[j]) + np.matmul(ode_jac(xs), s)
+        # sol = odeint(dS_dt, np.zeros(n_vars), self.timepoints, **kwargs)
+        # return sol
 
 
     ############## Sam Clamons ###########