Proper persistency condition

rate_independence.pl

@@ -14,6 +14,7 @@
 ).
 
 
+:- use_module(library(clpfd)).
 :- use_module(doc).
 
 
@@ -432,16 +433,22 @@ add_merge_reaction(M, Reaction1, Reaction2):-
 
   add_reaction(Kin1*Kin2,Reac,Inhib,Prod,false).
 
+% Dynamic predicate for storing structural information of the graph
 :- dynamic(sur_invariant/1).
 :- dynamic(p_invariant/1).
 :- dynamic(output/1).
+:- dynamic(siphon/1).
 
+%! test_rate_independence_invariants is det.
+%
+% test sufficient condition for rate independence based on Petri net invariants
 test_rate_independence_invariants :-
   biocham_command,
   doc('Tests invariant sufficient conditions for rate independence of the current model for all output species (assuming well-formed kinetics).'),
   retractall(sur_invariant(_)),
   retractall(p_invariant(_)),
   retractall(output(_)),
+  retractall(siphon(_)),
   with_output_to(atom(_), invariants:find_invar_aux(4, is_place, is_transition, '#=<')),
   forall(
     retract(invariants:base_mol(B)),
@@ -470,6 +477,7 @@ test_rate_independence_invariants :-
   ),
   reaction_graph,
   get_current_graph(Id),
+  find_emptyable_siphons(Id),
   (
     output(_)
   ->
@@ -488,6 +496,9 @@ test_rate_independence_invariants :-
   ).
 
 
+%! test_rate_independence_invariants(+Output, -Id) is det.
+%
+% check the rate independence of species Output in reaction graph Id
 test_rate_independence_invariants(Output, Id) :-
   debug(rate, "Trying for output ~w", [Output]),
   (
@@ -503,6 +514,10 @@ test_rate_independence_invariants(Output, Id) :-
   iterate_on_parents(Output, Id).
 
 
+%! is_covered_by_p_invariant_with_single_surinv(+Species, ?Output, -P) is det.
+%
+% check if Species is covered by a P-invariant P, Output is an Output, and P
+% does not contain any other sur-invariant than Output
 is_covered_by_p_invariant_with_single_surinv(Species, Output, P) :-
   p_invariant(P),
   member(Species, P),
@@ -514,12 +529,18 @@ is_covered_by_p_invariant_with_single_surinv(Species, Output, P) :-
   ).
 
 
-is_covered_by_sur_invariant(Species) :-
-  sur_invariant(S),
-  member(Species, S),
-  !.
+%! is_persistent(+Species) is det.
+%
+% check sufficient condition for persistency (i.e. is not in any siphon that does not contain
+% the support of a P-invariant)
+is_persistent(Species) :-
+  siphon(S),  % union of all emptyable siphons
+  \+ member(Species, S).
 
 
+%! iterate_on_parents(+Species, +GraphId) is det.
+%
+% iteratively search if all parents either are presistent or in different P-invariants
 iterate_on_parents(Species, GraphId) :-
   assertz(seen(Species)),
   debug(rate, "Searching for ~w parents", [Species]),
@@ -549,9 +570,9 @@ iterate_on_parents(Species, GraphId) :-
           ->
             debug(rate, "parent ~w already seen", [Reactant])
           ;
-            is_covered_by_sur_invariant(Reactant)
+            is_persistent(Reactant)
           ->
-            debug(rate, "parent ~w non-decreasing", [Reactant])
+            debug(rate, "parent ~w persistent", [Reactant])
           ;
             is_covered_by_p_invariant_with_single_surinv(Reactant, _, P),
             \+ member(P, ParentPinvs)
@@ -570,6 +591,121 @@ iterate_on_parents(Species, GraphId) :-
   ).
 
 
+%! normalize(+Invar, -Normalized) is det.
+%
+% remove stoichiometry and sort
 normalize(Invar, Normalized) :-
   maplist(nusmv:reactant_to_name, Invar, NoStoichiometry),
   sort(NoStoichiometry, Normalized).
+
+
+%! find_emptyable_siphons(+GraphId) is det.
+%
+% find all siphons that do not contain the support of a p-invariant
+find_emptyable_siphons(GraphId) :-
+  findall(
+    Species,
+    (
+      item([parent: GraphId, kind: vertex, id: VertexId, item: Species]),
+      get_attribute(VertexId, kind=place)
+    ),
+    LLSpecies
+  ),
+  sort(LLSpecies, LSpecies),
+  length(LSpecies, N),
+  length(LVars, N),
+  LVars ins 0..1,
+  siphon_constraints(LSpecies, LVars),
+  pinv_constraints(LSpecies, LVars),
+  assertz(siphon([])),
+  (
+    repeat,
+    (
+      % only get maximal such siphons
+      get_max_constr(LSpecies, LVars),
+      labeling([ff, down], LVars)
+    ->
+      invariants:to_mols(LVars, LSpecies, Siphon),
+      retract(siphon(Current)),
+      ord_union(Siphon, Current, Union),
+      assertz(siphon(Union)),
+      fail
+    ;
+      debugging(rate)
+    ->
+      siphon(Union),
+      debug(rate, "Union of emptyable siphons: ~w", [Union])
+    ;
+      true
+    )
+  ->
+    true
+  ).
+
+
+%! pinv_constraints(+LSpecies, ?LVars) is det.
+%
+% add constraints that siphons should not cover any P-invariant
+pinv_constraints(LSpecies, LVars) :-
+  foreach(
+    p_invariant(P),
+    is_not_covered(P, LSpecies, LVars)
+  ).
+
+
+%! is_not_covered(+P, +LSpecies, ?LVars) is det.
+%
+% enforce that P is not covered
+is_not_covered(P, LSpecies, LVars) :-
+  maplist(get_var_from_species(LSpecies, LVars), P, PVars),
+  element(_, PVars, 0).
+
+
+%! get_var_from_species(+LSpecies, ?LVars, +Species, ?Var) is det.
+%
+% Var is the variable at the same position in LVars as Species in LSpecies
+get_var_from_species(LSpecies, LVars, Species, Var) :-
+  nth1(N, LSpecies, Species),
+  !,
+  nth1(N, LVars, Var).
+
+
+%! siphon_constraints(+LSpecies, ?LVars) is det.
+%
+% siphon constraints (or of products => or of reactants) on all reactions
+siphon_constraints(LSpecies, LVars) :-
+  findall(
+    (NR, NP),
+    (
+      item([kind: reaction, item: Reaction]),
+      reaction(Reaction, [reactants: R, products: P]),
+      normalize(R, NR),
+      normalize(P, NP)
+    ),
+    L
+  ),
+  siphon_constraints_aux(L, LSpecies, LVars).
+
+
+%! siphon_constraints_aux(+L, +LSpecies, ?LVars) is det.
+%
+% siphon constraint on a list of reactants/products
+siphon_constraints_aux([], _, _).
+
+siphon_constraints_aux([(NR, NP) | L], LSpecies, LVars) :-
+  maplist(get_var_from_species(LSpecies, LVars), NR, VR),
+  maplist(get_var_from_species(LSpecies, LVars), NP, VP),
+  sum(VR, #=, SR),
+  sum(VP, #=, SP),
+  SP #> 0 #==> SR #> 0,
+  siphon_constraints_aux(L, LSpecies, LVars).
+
+
+%! get_max_constr(+LSpecies, ?LVars) is det.
+%
+% enforce that future siphons have at least one species not yet covered
+get_max_constr(LSpecies, LVars) :-
+  siphon(S),
+  ord_subtract(LSpecies, S, Others),
+  maplist(get_var_from_species(LSpecies, LVars), Others, SV),
+  element(_, SV, 1).

rate_independence.plt

@@ -42,7 +42,6 @@ test('input_sinks', [
   sub_string(OutputLong, 0, 9, _, Output).
 
 test('invariants', [
-  fixme('Work in progress'),
   condition(flag(slow_test, true, true)),
   forall((expand_file_name('./library/examples/tests/rate_independence_*.bc', FileList), member(File, FileList))),
   setup((with_output_to(atom(ExpectedLong), command(load(File))), sub_string(ExpectedLong, 0, 9, _, Expected))),