Merge pull request #9 from ROBACON/reversible_rate_update

Reversible rate update

Author: fabriciocravo

debugging_script.py

@@ -33,14 +33,15 @@
              test_assignment_similar_species, test_blocked_names, test_blocked_names_2,
              test_update_parameter_for_multi_model, test_update_parameter_through_str,
              test_update_multiple_parameters_in_expression, test_update_parameter_with_unit,
-             test_species_value_modification, test_all_value_modification]
+             test_species_value_modification, test_all_value_modification, test_new_reversible_reaction_notation,
+             test_2D_reaction_with_units]
 
 # test_no_species_in_asg
 # test_illegal_unit_op_in_assignment
 # temporary_test_removal = [test_parameter_fit_with_units, test_multiple_runs_fit, test_simple_fit]
 test_remov = [test_antimony_compose_model_gen, test_antimony_model]
 sub_test = test_list
-#sub_test = [test_antimony_compose_model_gen]
+# sub_test = [test_parameters_with_sbml]
 def perform_tests():
     any_failed = False
     for test in sub_test:

for_local_use.py

@@ -2,26 +2,20 @@
 from mobspy import modules
 
 if __name__ == '__main__':
+    from mobspy import *
 
-    # TODO: Ask Matthias about removed tests
 
     # Replace parameters using units
-    A = BaseSpecies()
-    A.a1, A.a2
-    B = New(A)
-    B.b1, B.b2
-    k1 = ModelParameters(1)
+    Color, Location = BaseSpecies()
+    Color.red, Color.blue
+    Location.here, Location.there
+    Something = Color*Location
 
-    B >> Zero [k1]
+    2*Something >> 3*Something [lambda r1, r2: 1*u.decimeter**2/u.h if Location(r1) == Location(r2) else 0.5*u.decimeter**2/u.h]
 
-    B(100), B.b2(100)
-    S = Simulation(B)
-    S.level = -1
-    S.compile()
-
-    S.update_model([All[B], 200/u.l])
-
-    print(S.generate_sbml()[0])
+    S = Simulation(Something)
+    S.volume = 1*u.m**2
+    print(S.compile())
 
     # print(A.get_characteristics())
 

mobspy/modules/function_rate_code.py

@@ -55,7 +55,6 @@ def extract_reaction_rate(combination_of_reactant_species, reactant_string_list
                                                      reaction_rate_function, type_of_model, is_count)
 
     elif function_rate_arguments is not None:
-
         # [''] means that it is a function that takes no arguments (empty signature)
         if function_rate_arguments != ['']:
             arguments = prepare_arguments_for_callable(combination_of_reactant_species,
@@ -87,12 +86,14 @@ def extract_reaction_rate(combination_of_reactant_species, reactant_string_list
 
             # Having an expression variable implies it is a constructed expression - not mass action
             if len(rate._expression_variables) > 0:
-                reaction_rate_string, _ = rate.generate_string_operation(skip_check=skip_check)
+                reaction_rate_string, _ = rate.generate_string_operation(skip_check=skip_check, dimension=dimension)
 
             # Having no expression variables implies it is a constant for mass - action kinetics.
             elif len(rate._expression_variables) == 0:
+
                 rate_for_mass_action, is_count = \
-                    rate.generate_string_operation(reaction_order=len(reactant_string_list), skip_check=skip_check)
+                    rate.generate_string_operation(reaction_order=len(reactant_string_list), dimension=dimension,
+                                                   skip_check=skip_check)
                 parameters_in_rate = rate._parameter_set
                 parameters_in_reaction = parameters_in_reaction.union(parameters_in_rate)
 

mobspy/modules/meta_class.py

@@ -4,7 +4,9 @@
     Bear in mind that they are not actually Python meta-classes. The first design utilized this feature but know there
     are just regular classes
 """
-from mobspy.simulation_logging.log_scripts import error as simlog_error, debug as simlog_debug
+from copy import deepcopy
+
+from mobspy.simulation_logging.log_scripts import error as simlog_error, debug as simlog_debug, error
 from mobspy.modules.species_string_generator import construct_all_combinations as ssg_construct_all_combinations
 from mobspy.modules.assignments_implementation import Assign as asgi_Assign, Asg as asgi_Asg
 from pint import Quantity
@@ -43,6 +45,30 @@ def override_get_item(cls, object_to_return, item):
         cls.last_rate = item
         return object_to_return
 
+    @classmethod
+    def process_rate(cls, rate):
+
+        if isinstance(rate, (np_int_, np_float_)):
+            rate = float(rate)
+
+        if isinstance(rate, Species) \
+                or isinstance(rate, Reacting_Species) \
+                or isinstance(rate, Reactions):
+            simlog_error('Reaction rate of type ' + str(type(_Last_rate_storage.last_rate)) + ' not valid',
+                         stack_index=3)
+
+        if not (type(rate) == int or type(rate) == float
+                or callable(rate) or type(rate) == str
+                or isinstance(rate, me_OverrideQuantity)
+                or isinstance(rate, Quantity)
+                or isinstance(rate, mp_Mobspy_Parameter)
+                or rate is None
+                or isinstance(rate, me_ExpressionDefiner)):
+            simlog_error('Reaction rate of type ' + str(type(rate)) + ' not valid',
+                         stack_index=3)
+
+        return rate
+
 
 class Reactions:
     """
@@ -95,7 +121,7 @@ def __getitem__(self, item):
         """
         return _Last_rate_storage.override_get_item(self, item)
 
-    def __init__(self, reactants, products):
+    def __init__(self, reactants, products, rate=None):
         """
             Constructor of the reaction object. For the object construction only the reactants and products are
             necessary - the order and the rate are assigned later by the compiler
@@ -117,7 +143,37 @@ def __init__(self, reactants, products):
                          'A reaction was defined with the assignment context activated.'
                          'The assignment context was deactivated. Please try to redefine the model')
 
-        # Add characteristics in Cts_context to each reactant and product
+        # Setup rate - consider reversible reactions
+        try:
+            flag_len = len(_Last_rate_storage.last_rate) == 2
+        except:
+            flag_len = False
+
+        if flag_len and rate is None:
+
+            store_last_rate = _Last_rate_storage.last_rate[0]
+            Reactions(reactants=products, products=reactants, rate=_Last_rate_storage.last_rate[1])
+            rate = _Last_rate_storage.process_rate(store_last_rate)
+
+        elif rate is not None and not flag_len:
+
+            rate = _Last_rate_storage.process_rate(rate)
+
+        elif rate is None and not flag_len:
+
+            rate = _Last_rate_storage.process_rate(_Last_rate_storage.last_rate)
+
+        elif rate is not None and flag_len:
+
+            rate = _Last_rate_storage.process_rate(rate)
+
+        else:
+            simlog_error('To many rate provided', stack_index=3)
+        self.rate = rate
+        if _Last_rate_storage.last_rate is not None:
+            _Last_rate_storage.last_rate = None
+
+            # Add characteristics in Cts_context to each reactant and product
         if len(Species.meta_specie_named_any_context) != 0:
             for j in Species.meta_specie_named_any_context:
                 for r in reactants:
@@ -143,29 +199,12 @@ def __init__(self, reactants, products):
 
         # Assign default order
         self.order = None
-        # Cast np to float
-        if isinstance(_Last_rate_storage.last_rate, (np_int_, np_float_)):
-            _Last_rate_storage.last_rate = float(_Last_rate_storage.last_rate)
 
-        if isinstance(_Last_rate_storage.last_rate, Species) \
-                or isinstance(_Last_rate_storage.last_rate, Reacting_Species) \
-                or isinstance(_Last_rate_storage.last_rate, Reactions):
-            simlog_error('Reaction rate of type ' + str(type(_Last_rate_storage.last_rate)) + ' not valid',
-                         stack_index=3)
-
-        if not (type(_Last_rate_storage.last_rate) == int or type(_Last_rate_storage.last_rate) == float
-                or callable(_Last_rate_storage.last_rate) or type(_Last_rate_storage.last_rate) == str
-                or isinstance(_Last_rate_storage.last_rate, me_OverrideQuantity)
-                or isinstance(_Last_rate_storage.last_rate, Quantity)
-                or isinstance(_Last_rate_storage.last_rate, mp_Mobspy_Parameter)
-                or _Last_rate_storage.last_rate is None
-                or isinstance(_Last_rate_storage.last_rate, me_ExpressionDefiner)):
-            simlog_error('Reaction rate of type ' + str(type(_Last_rate_storage.last_rate)) + ' not valid',
-                         stack_index=3)
+        # Define reversible reaction with two rates
 
-        self.rate = _Last_rate_storage.last_rate
-        if _Last_rate_storage.last_rate is not None:
-            _Last_rate_storage.last_rate = None
+        # Cast np to float
+        #if _Last_rate_storage.last_rate is not None:
+        #    _Last_rate_storage.last_rate = None
 
         # Here we extract all involved objects to pact them in a set
         # This is done to find the reactions associated with the species when the Compiler is started

mobspy/modules/mobspy_expressions.py

@@ -806,7 +806,7 @@ def __getattr__(self, item):
         return super().__getattr__(item)
 
     # Write string operation return and unit
-    def generate_string_operation(self, skip_check=False, reaction_order=None):
+    def generate_string_operation(self, skip_check=False, reaction_order=None, dimension=None):
         """
             Converts the expression from what has been stored to a string format for the sbml file
 
@@ -815,10 +815,11 @@ def generate_string_operation(self, skip_check=False, reaction_order=None):
         """
         ur = u.unit_registry_object
 
-        if self._dimension is None:
-            dimension = 3
-        else:
-            dimension = self._dimension
+        if dimension is None:
+            if self._dimension is None:
+                dimension = 3
+            else:
+                dimension = self._dimension
 
         operation = str(self._operation)
 
@@ -829,9 +830,9 @@ def generate_string_operation(self, skip_check=False, reaction_order=None):
             self._count_in_expression = True
 
         if self._has_units == 'T' and self._expression_variables == set() and reaction_order is not None:
-            if self._unit_count_op.units == (1 / ur.second):
+            if self._unit_count_op.units == (1 / ur.second).units:
                 return operation, True
-            elif self._unit_conc_op.units == (1*ur.decimeter**(3*reaction_order)/(ur.second*ur.decimeter**dimension)):
+            elif self._unit_conc_op.units == (ur.decimeter**(dimension * (reaction_order - 1)) / ur.second):
                 return operation, False
 
         c1 = self._has_units == 'T'

test_plot_images/read_me.txt

@@ -0,0 +1,3 @@
+This directory is used for test purposes. It received the images generated by the test script.
+Any plotting tests should create images here.
+

test_script.py

@@ -769,7 +769,7 @@ def test_bi_dimensional_rates():
 def test_dimension_in_function_only():
     A = BaseSpecies()
 
-    A + A >> 3 * A[lambda: 1 * u.milliliter / u.second]
+    A + A >> 3 * A [lambda: 1 * u.milliliter / u.second]
 
     A(1)
     S = Simulation(A)
@@ -958,7 +958,7 @@ def test_parameters_with_sbml():
         for data_for_sbml in parameter_sweep:
             model_str += order_model_str(data_for_sbml)
 
-    assert compare_model(model_str, 'test_tools/model_31.txt')
+    assert compare_model_ignore_order(model_str, 'test_tools/model_31.txt')
 
 
 def test_shared_parameter_name():
@@ -1942,7 +1942,7 @@ def test_update_parameter_for_multi_model():
 
     k1 = ModelParameters([1, 2, 3])
 
-    A >> Zero [2*k1]
+    A >> mobspy.Zero [2*k1]
 
     A(100)
     S1 = Simulation(A)
@@ -1952,7 +1952,7 @@ def test_update_parameter_for_multi_model():
     # print(S.generate_sbml()[0])
 
     A.reset_reactions()
-    B >> Zero [1]
+    B >> mobspy.Zero [1]
 
     B(200)
     S2 = Simulation(A | B)
@@ -1974,7 +1974,7 @@ def test_update_parameter_through_str():
     A = BaseSpecies()
     k1 = ModelParameters(0.00000001)
 
-    A >> Zero [k1]
+    A >> mobspy.Zero [k1]
 
     S = Simulation(A)
     S.level = -1
@@ -1989,7 +1989,7 @@ def test_update_multiple_parameters_in_expression():
     A = BaseSpecies()
     k1, k2 = ModelParameters(0.00000001, 10)
 
-    A >> Zero [k1 / (10 + k2 ** 4)]
+    A >> mobspy.Zero [k1 / (10 + k2 ** 4)]
 
     S = Simulation(A)
     S.level = -1
@@ -2005,7 +2005,7 @@ def test_update_parameter_with_unit():
     A = BaseSpecies()
     k1 = ModelParameters(1 / u.h)
 
-    A >> Zero[k1]
+    A >> mobspy.Zero[k1]
 
     S = Simulation(A)
     S.level = -1
@@ -2024,7 +2024,7 @@ def test_species_value_modification():
     B.b1, B.b2
     k1 = ModelParameters(1)
 
-    B >> Zero[k1]
+    B >> mobspy.Zero[k1]
 
     B(100), B.b2(100)
     S = Simulation(B)
@@ -2042,7 +2042,7 @@ def test_all_value_modification():
     B.b1, B.b2
     k1 = ModelParameters(1)
 
-    B >> Zero[k1]
+    B >> mobspy.Zero[k1]
 
     B(100), B.b2(100)
     S = Simulation(B)
@@ -2051,4 +2051,36 @@ def test_all_value_modification():
 
     S.update_model([All[B], 200 / u.l])
 
-    assert compare_model_ignore_order(S.generate_sbml()[0], 'test_tools/model_62.txt')
+    assert compare_model_ignore_order(S.generate_sbml()[0], 'test_tools/model_62.txt')
+
+
+def test_new_reversible_reaction_notation():
+
+    A = BaseSpecies()
+    k1, k2 = ModelParameters(1, 1)
+
+    A >> mobspy.Zero[1, 1]
+    A >> mobspy.Zero[1 / (k1 + k2), k1]
+
+    A >> 2 * A[10]
+
+    S = Simulation(A)
+    S.level = -1
+    assert compare_model(S.compile(), 'test_tools/model_63.txt')
+
+
+def test_2D_reaction_with_units():
+
+    Color, Location = BaseSpecies()
+    Color.red, Color.blue
+    Location.here, Location.there
+    Something = Color * Location
+
+    rate = lambda r1, r2: 1 * u.decimeter ** 2 / u.h if Location(r1) == Location(r2) \
+        else 0.5 * u.decimeter ** 2 / u.h
+    2 * Something >> 3 * Something[rate]
+
+    S = Simulation(Something)
+    S.level = -1
+    S.volume = 1 * u.m ** 2
+    assert compare_model(S.compile(), 'test_tools/model_64.txt')

test_tools/model_63.txt

@@ -0,0 +1,20 @@
+
+Species
+A,0
+
+Mappings
+A :
+A
+
+Parameters
+k1,1
+k2,1
+volume,1
+
+Reactions
+reaction_0,{'re': [(1, 'A')], 'pr': [(2, 'A')], 'kin': 'A * 10'}
+reaction_1,{'re': [(1, 'A')], 'pr': [], 'kin': 'A * (1/(k1+k2))'}
+reaction_2,{'re': [(1, 'A')], 'pr': [], 'kin': 'A * 1'}
+reaction_3,{'re': [], 'pr': [(1, 'A')], 'kin': '1 * volume'}
+reaction_4,{'re': [], 'pr': [(1, 'A')], 'kin': 'k1 * volume'}
+