conservation_laws.pl 4.37 KB
Newer Older
Thierry Martinez's avatar
Thierry Martinez committed
1 2 3
:- module(
  conservation_laws,
  [
4 5 6 7 8 9 10
    add_conservation/1,
    delete_conservation/1,
    delete_conservations/0,
    list_conservations/0,
    check_conservations/0,
    search_conservations/0,
    search_conservations/1
Thierry Martinez's avatar
Thierry Martinez committed
11 12
  ]
).
13 14 15 16 17


add_conservation(Conservation) :-
  biocham_command,
  type(Conservation, solution),
SOLIMAN Sylvain's avatar
SOLIMAN Sylvain committed
18 19 20 21 22
  % Avoid a SWI-Prolog warning about long strings…
  doc(
    'declares a new mass conservation law for all molecules given with the
    corresponding weight in \\argument{Conservation}. During a numerical
    simulation, one of those variables will be
23 24
    eliminated thanks to this conservation law. Be careful if you declare
    conservation laws and then plot the result of a previous simulation, the
Thierry Martinez's avatar
FOLTL  
Thierry Martinez committed
25 26
    caption might not be correct.
When added, the conservation law will be checked against the reactions
27 28 29
    (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.'
SOLIMAN Sylvain's avatar
SOLIMAN Sylvain committed
30
  ),
31
  solution_to_conservation(Conservation, C),
Thierry Martinez's avatar
Thierry Martinez committed
32
  add_item([kind: conservation, item: Conservation, id: Id]),
33
  set_annotation(Id, conservation_list, C).
34 35 36 37 38 39


delete_conservation(Conservation) :-
  biocham_command,
  type(Conservation, solution),
  doc('removes the given mass conservation law.'),
40 41
  % model: current_model?
  find_item([id: Id, type: conservation, item: Conservation]),
42 43 44 45 46 47
  delete_item(Id).


delete_conservations :-
  biocham_command,
  doc('removes all mass conservation laws.'),
48 49
  % model: current_model?
  delete_items([kind: conservation]).
50 51 52 53 54


list_conservations :-
  biocham_command,
  doc('prints out all the mass conservation laws.'),
55 56
  % model: current_model?
  list_items([kind: conservation]).
57 58 59 60 61


check_conservations :-
  biocham_command,
  doc(
62
    'checks all conservation laws against reactions, and if necessary kinetics
SOLIMAN Sylvain's avatar
SOLIMAN Sylvain committed
63 64
    (see also \\command{add_conservation/1}).'
  ),
Thierry Martinez's avatar
Thierry Martinez committed
65
  ode_system,
66
  forall(
67
    (
68
      % model: current_model? compute_ode already supposes a single model
69 70
      % otherwise for each selected model, get all conservations and check
      item([kind: conservation, item: Conserv, id: Id]),
71 72
      get_annotation(Id, conservation_list, C)
    ),
73
    check_conservation(C, Conserv)
74
  ).
75 76 77 78


search_conservations :-
  biocham_command,
SOLIMAN Sylvain's avatar
SOLIMAN Sylvain committed
79
  doc('calls \\command{search_conservations/1} with the default size of 4.'),
80
  search_conservations(4).
81 82


SOLIMAN Sylvain's avatar
SOLIMAN Sylvain committed
83
search_conservations(Integer) :-
84
  biocham_command,
SOLIMAN Sylvain's avatar
SOLIMAN Sylvain committed
85
  type(Integer, integer),
86
  doc(
SOLIMAN Sylvain's avatar
SOLIMAN Sylvain committed
87 88 89 90
    'computes and displays the P-invariants of the system up to the maximal
    size \\argument{Integer}. Such P-invariants are particular mass
    conservation laws that are independent from the kinetics.'
  ),
91
  invariants:find_pinvar(Integer).
92 93


94 95 96 97 98 99 100 101 102
%%% 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).


103
check_conservation(C, Conserv) :-
104 105 106
  (
    check_conserv_reactions(C)
  ->
107
    print_message(informational, conservation_checked(Conserv, reactions, true))
108
  ;
109 110 111 112
    print_message(
      informational,
      conservation_checked(C, reactions, false)
    ),
113 114
    check_conserv_num(C)
  ->
115
    print_message(informational, conservation_checked(Conserv, kinetics, true))
116
  ;
117 118
    print_message(informational, conservation_checked(Conserv, kinetics, false)),
    print_message(warning, conservation_trusted(Conserv))
119 120 121 122 123 124
  ).


check_conserv_reactions(C) :-
  forall(
    (
125 126
      % model: current_model?
      item([kind: reaction, item: Item]),
127 128 129 130 131 132 133 134 135
      reaction(Item, _Kinetics, Reactants, Products)
    ),
    (
      scalar_conservation(Reactants, C, Moiety),
      scalar_conservation(Products, C, Moiety)
    )
  ).


136 137 138 139 140 141 142 143 144 145 146
check_conserv_num(C) :-
  get_weighted_ode(C, WeightedSum),
  simplify(WeightedSum, Sum),
  Sum == 0.


get_weighted_ode([], 0).

get_weighted_ode([K*M | L], K*O + S) :-
  ode(M, O),
  get_weighted_ode(L, S).
147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162


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.


163 164 165 166 167 168 169 170 171 172 173 174
prolog:message(conservation_trusted(C)) -->
  ['Conservation ~w trusted without proof'-[C]].

prolog:message(conservation_checked(C, From, Truth)) -->
  {
    Truth == true
  ->
    Not = ''
  ;
    Not = 'not '
  },
  ['Conservation ~w ~wchecked from ~w'-[C, Not, From]].
175

176
% vi: set ts=2 sts=2 sw=2