properly handle equivalent/inverse reactions

lemon.pl

@@ -84,12 +84,13 @@ export_lemon_stream(Id, Stream) :-
       Edge = (VertexA -> VertexB),
       species(VertexA, IdA),
       species(VertexB, IdB),
-      count(influences, Count),
       get_attribute(EdgeId, origin = Origin),
       forall(
         member((Sign, Reaction), Origin),
         (
-          reactions(Reaction, RId),
+          count(influences, Count),
+          normalize_reaction(Reaction, Normalized, _Negated),
+          reactions(Normalized, RId),
           sign(Sign, SignNumber),
           format(
             Stream,
@@ -125,16 +126,27 @@ export_lemon_stream(Id, Stream) :-
 %
 % associates a unique number to each reaction
 % such that reaction that are stoichiometrically inverse get opposite numbers
-%
-% @tbd handle reverse reactions
+% and stoichiometrically equivalent get same number
 number_reactions :-
   retractall(reactions(_, _)),
-  set_counter(reactions, 0),
+  set_counter(reactions, 1),
   forall(
     item([kind: reaction, item: Item]),
     (
-      count(reactions, Count),
-      assertz(reactions(Item, Count))
+      normalize_reaction(Item, Normalized, Negated),
+      (
+        reactions(Normalized, _Count)
+      ->
+        true
+      ;
+        reactions(Negated, Index)
+      ->
+        Count is -Index,
+        assertz(reactions(Normalized, Count))
+      ;
+        count(reactions, Count),
+        assertz(reactions(Normalized, Count))
+      )
     )
   ).
 
@@ -144,3 +156,15 @@ number_reactions :-
 % associates number Number to a sign Atom
 sign('+', 1).
 sign('-', -1).
+
+
+%! normalize_reaction(+Item, -NormalizedStoichiometry, -Negated) is det.
+%
+% compute a sorted stoichiometry list for the reaction Item as
+% NormalizedStoichiometry and a sorted list for the inverse reaction as
+% Negated
+normalize_reaction(Item, NormalizedStoichiometry, Negated) :-
+  reaction_editor:get_stoichiometry_and_kinetics(Item, Stoichiometry, _),
+  sort(Stoichiometry, NormalizedStoichiometry),
+  maplist(reaction_editor:negate_coefficient, Stoichiometry, NegatedStoich),
+  sort(NegatedStoich, Negated).

reaction_editor.pl

@@ -600,15 +600,9 @@ make_reaction(
 
 compute_ode_for_reactions :-
   forall(
+    item([kind: reaction, item: Item]),
     (
-      item([kind: reaction, item: Item]),
-      reaction(Item, Kinetics, Reactants, Inhibitors, Products)
-    ),
-    (
-      kinetics(Reactants, Inhibitors, Kinetics, KineticsExpression),
-      maplist(negate_coefficient, Reactants, NegatedReactants),
-      append(NegatedReactants, Products, StoichiometricSolution),
-      simplify_solution(StoichiometricSolution, SimplifiedSolution),
+      get_stoichiometry_and_kinetics(Item, SimplifiedSolution, KineticsExpression),
       add_molecules(SimplifiedSolution, KineticsExpression)
     )
   ).
@@ -627,4 +621,19 @@ add_molecule(Coefficient * Molecule, Kinetics) :-
   ode_add_expression_to_molecule(Coefficient * Kinetics, Molecule).
 
 
-negate_coefficient(C*X, (-C)*X).
+negate_coefficient(C*X, D*X) :-
+  (
+    number(C)
+  ->
+    D is -C
+  ;
+    D = -C
+  ).
+
+
+get_stoichiometry_and_kinetics(Item, SimplifiedSolution, KineticsExpression) :-
+  reaction(Item, Kinetics, Reactants, Inhibitors, Products),
+  kinetics(Reactants, Inhibitors, Kinetics, KineticsExpression),
+  maplist(negate_coefficient, Reactants, NegatedReactants),
+  append(NegatedReactants, Products, StoichiometricSolution),
+  simplify_solution(StoichiometricSolution, SimplifiedSolution).