checking conservations from kinetics too, fixes stuff in kinetics on the way

conservation_laws.pl

@@ -65,6 +65,7 @@ check_conservations :-
     'checks all conservation laws against reactions, and if necessary kinetics
     (see also \\command{add_conservation/1}).'
   ),
+  compute_ode,
   forall(
     (
       item([model: current_model, kind: conservation, id: Id]),
@@ -106,11 +107,12 @@ check_conservation(C) :-
   ->
     message('~w checked from reactions', [C])
   ;
+    message('~w not checked from reactions', [C]),
     check_conserv_num(C)
   ->
     message('~w checked from kinetics', [C])
   ;
-    message('~w trusted but not checked', [C]),
+    message('~w not checked from kinetics either, but trusted anyway', [C]),
     throw(error(unchecked_conservation_law(C)))
   ).
 
@@ -128,8 +130,18 @@ check_conserv_reactions(C) :-
   ).
 
 
-check_conserv_num(_) :-
-  fail.
+check_conserv_num(C) :-
+  get_weighted_ode(C, WeightedSum),
+  simplify(WeightedSum, Sum),
+  write(Sum), nl,
+  Sum == 0.
+
+
+get_weighted_ode([], 0).
+
+get_weighted_ode([K*M | L], K*O + S) :-
+  ode(M, O),
+  get_weighted_ode(L, S).
 
 
 scalar_conservation([], _, 0).

conservation_laws.plt

@@ -69,6 +69,20 @@ test(
   with_output_to(atom(_), command(check_conservations)).
 
 
+test(
+  'check_conservations is ok with kinetically conserved moieties',
+  [
+    setup(models:clear_model),
+    cleanup(models:clear_model)
+  ]
+) :-
+  % writing with 'command' produces a parsing error on k…
+  add_reaction('for'(k*[a], '=>'(a, b))),
+  add_reaction('for'(k*[a], '=>'(a, 2*a))),
+  command(add_conservation(a)),
+  with_output_to(atom(_), command(check_conservations)).
+
+
 test(
   'check_conservations is unhappy with non-conserved moieties',
   [

kinetics.pl

@@ -10,6 +10,8 @@ kinetics(Kinetics, Reactants, Expression) :-
     Kinetics = 'MA'(Coefficient)
   ->
     mass_action(Coefficient, Reactants, Expression)
+  ;
+    Expression = Kinetics
   ).
 
 
@@ -21,5 +23,5 @@ mass_action(Coefficient, Reactants, Expression) :-
 concentration_product([], 1).
 
 concentration_product([Coefficient * Reactant | Reactants], Product) :-
-  Product = Coefficient * [Reactant] * ProductTail,
+  Product = [Reactant] ^ Coefficient * ProductTail,
   concentration_product(Reactants, ProductTail).

ode.pl

@@ -20,8 +20,8 @@ list_ODE :-
 compute_ode :-
   retractall(ode(_, _)),
   \+ (
-    item([model: current_model, kind: rule, item: Item]),
-    rule(Item, Kinetics, Reactants, Products),
+    item([model: current_model, kind: reaction, item: Item]),
+    reaction(Item, Kinetics, Reactants, Products),
     \+ (
       kinetics(Kinetics, Reactants, KineticsExpression),
       add_molecules(Reactants, -KineticsExpression),
@@ -58,7 +58,7 @@ add_molecule(Coefficient * Molecule, Kinetics) :-
   ;
     Expression = 0
   ),
-  asserta(ode(Molecule, Expression + Kinetics / Coefficient)).
+  asserta(ode(Molecule, Expression + Coefficient * Kinetics )).
 
 
 print_ode :-

simplify_arithmetic.pl

@@ -118,6 +118,17 @@ simplify_addition(A, Value * B, Result) :-
   ValueOpp is -Value,
   Result = A - ValueOpp * B.
 
+simplify_addition(-A, B, C) :-
+  !,
+  simplify_difference(B, A, C).
+
+simplify_addition(A, -B, C) :-
+  !,
+  simplify_difference(A, B, C).
+
+simplify_addition(A, A, 2 * A) :-
+  !.
+
 simplify_addition(A, B, A + B).
 
 
@@ -148,6 +159,17 @@ simplify_difference(A, Value * B, Result) :-
   ValueOpp is -Value,
   Result = A + ValueOpp * B.
 
+simplify_difference(A, A, 0) :-
+  !.
+
+simplify_difference(A, -B, C) :-
+  !,
+  simplify_addition(A, B, C).
+
+simplify_difference(-A, B, -C) :-
+  !,
+  simplify_addition(A, B, C).
+
 simplify_difference(A, B, A - B).