Merge branch 'master' of git+ssh://scm.gforge.inria.fr/gitroot/biocham/biocham

conservation_laws.pl

@@ -25,31 +25,31 @@ add_conservation(Conservation) :-
     caption might not be correct.'
   ),
   doc(
-    'When added, the conservation law will be checked against the rules (i.e.
-    purely from stoichiometry), if that fails against the kinetics. Since
-    these checks are not complete, even a failure will be accepted with a
-    warning.'
+    'When added, the conservation law will be checked against the reactions
+    (i.e. purely from stoichiometry), if that fails against the kinetics.
+    Since these checks are not complete, even a failure will be accepted with
+    a warning.'
   ),
-  solution_to_list(Conservation, C),
-  add_item(conservation, C).
+  solution_to_conservation(Conservation, C),
+  add_item(conservation, [], Conservation, Id),
+  set_annotation(Id, conservation_list, C).
 
 
 delete_conservation(Conservation) :-
   biocham_command,
   type(Conservation, solution),
   doc('removes the given mass conservation law.'),
-  solution_to_list(Conservation, C),
-  find_item([model: current_model, id: Id, type: conservation, item: C]),
+  find_item(
+    [model: current_model, id: Id, type: conservation, item: Conservation]),
   delete_item(Id).
 
 
 delete_conservations :-
   biocham_command,
   doc('removes all mass conservation laws.'),
-  \+ (
+  forall(
     item([model: current_model, kind: conservation, id: Id]),
-    delete_item(Id),
-    fail
+    delete_item(Id)
   ).
 
 
@@ -62,10 +62,17 @@ list_conservations :-
 check_conservations :-
   biocham_command,
   doc(
-    'checks all conservation laws against rules, and if necessary kinetics
+    'checks all conservation laws against reactions, and if necessary kinetics
     (see also \\command{add_conservation/1}).'
   ),
-  true.
+  compute_ode,
+  forall(
+    (
+      item([model: current_model, kind: conservation, id: Id]),
+      get_annotation(Id, conservation_list, C)
+    ),
+    check_conservation(C)
+  ).
 
 
 search_conservations :-
@@ -85,4 +92,74 @@ search_conservations(Integer) :-
   true.
 
 
+%%% End of public API
+
+
+solution_to_conservation(Solution, Conservation) :-
+  devdoc('transforms a solution to a factorized list of coeff*object'),
+  reaction_editor:solution_to_list(Solution, List),
+  reaction_editor:simplify_solution(List, Conservation).
+
+
+check_conservation(C) :-
+  (
+    check_conserv_reactions(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 not checked from kinetics either, but trusted anyway', [C]),
+    throw(error(unchecked_conservation_law(C)))
+  ).
+
+
+check_conserv_reactions(C) :-
+  forall(
+    (
+      item([model: current_model, kind: reaction, item: Item]),
+      reaction(Item, _Kinetics, Reactants, Products)
+    ),
+    (
+      scalar_conservation(Reactants, C, Moiety),
+      scalar_conservation(Products, C, Moiety)
+    )
+  ).
+
+
+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).
+
+scalar_conservation([K*M | L], C, Total) :-
+  (
+    member(I*M, C)
+  ->
+    true
+  ;
+    I = 0
+  ),
+  scalar_conservation(L, C, T1),
+  Total is I*K + T1.
+
+
+message(M, O) :-
+  format(M, O),
+  nl.
+
 % vi: set ts=2 sts=2 sw=2

conservation_laws.plt

@@ -9,10 +9,11 @@ test(
    [
       setup(models:clear_model),
       cleanup(models:clear_model),
-      true(Conservations == [[1*'a-a', 2*'a']])
+      forall(member(Input, ['a-a' + 2*'a', 'b', 'a' + 'b'])),
+      true(Conservations == [Input])
    ]
 ) :-
-   command(add_conservation(a-a + 2*a)),
+   command(add_conservation(Input)),
    all_items([model: current_model, kind: conservation], Conservations).
 
 %%% FIXME should test that conservations are properly used for numerical
@@ -21,9 +22,9 @@ test(
 test(
    'delete_conservation deletes the correct conservation',
    [
-      setup((models:clear_model, add_some_conservations)),
+      setup(set_some_conservations),
       cleanup(models:clear_model),
-      true(Conservations == [[1*'a-a', 2*'a'], [1*'c-c', 2*'c']])   
+      true(Conservations == ['a-a'+2*a,'c-c'+2*c])
    ]
 ) :-
    command(delete_conservation(b-b + 2*b)),
@@ -33,7 +34,7 @@ test(
 test(
    'delete_conservations deletes all conservations',
    [
-      setup((models:clear_model, add_some_conservations)),
+      setup(set_some_conservations),
       cleanup(models:clear_model),
       true(Conservations == [])
    ]
@@ -42,18 +43,44 @@ test(
    all_items([model: current_model, kind: conservation], Conservations).
 
 
+test(
+   'list_conservations lists all conservations',
+   [
+      setup(set_some_conservations),
+      cleanup(models:clear_model),
+      true(Conservations ==
+        '[0] a-a+2*a\n[1] b-b+2*b\n[2] c-c+2*c\n')
+   ]
+) :-
+  with_output_to(atom(Conservations), command(list_conservations)).
+
+
 test(
   'check_conservations is ok with trivial P-invariants',
   [
     setup(models:clear_model),
-    cleanup(models:clear_model),
-    blocked(not_implemented)
+    cleanup(models:clear_model)
   ]
 ) :-
   % not written 'a => b' to allow flycheck by separate compilation…
   command(add_reaction('=>'(a, b))),
+  command(add_reaction('=>'(b, a))),
   command(add_conservation(a + b)),
-  command(check_conservations).
+  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(
@@ -61,20 +88,21 @@ test(
   [
     setup(models:clear_model),
     cleanup(models:clear_model),
-    error(unproven_conservation_law([1*'a', 2*'b'])),
-    fixme(not_implemented)
+    throws(error(unchecked_conservation_law([1*'a', 2*'b'])))
   ]
 ) :-
   % not written 'a => b' to allow flycheck by separate compilation…
+  command(add_reaction('=>'(2*a, b))),
   command(add_reaction('=>'(a, b))),
   command(add_conservation(a + 2*b)),
-  command(check_conservations).
+  with_output_to(atom(_), command(check_conservations)).
 
 
-add_some_conservations :-
-   command(add_conservation(a-a + 2*a)),
-   command(add_conservation(b-b + 2*b)),
-   command(add_conservation(c-c + 2*c)).
+set_some_conservations :-
+  models:clear_model,
+  command(add_conservation(a-a + 2*a)),
+  command(add_conservation(b-b + 2*b)),
+  command(add_conservation(c-c + 2*c)).
 
 
 :- end_tests(conservation_laws).

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

@@ -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).