Commit d6872967 authored by FAGES Francois's avatar FAGES Francois
Browse files

merging

Merge branch 'develop' of gitlab.inria.fr:lifeware/biocham into develop
parents c9d56e2c d8ee360d
......@@ -14,6 +14,7 @@
).
:- use_module(library(clpfd)).
:- use_module(doc).
......@@ -432,23 +433,25 @@ add_merge_reaction(M, Reaction1, Reaction2):-
add_reaction(Kin1*Kin2,Reac,Inhib,Prod,false).
:- dynamic(sur_invariant/1).
% Dynamic predicate for storing structural information of the graph
:- 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)),
(
normalize(B, BB),
assertz(sur_invariant(BB)),
debug(rate, "sur_invariant: ~w", [BB]),
(
B = [Out]
->
......@@ -470,16 +473,14 @@ test_rate_independence_invariants :-
),
reaction_graph,
get_current_graph(Id),
find_emptyable_siphons(Id),
(
output(_)
->
(
forall(
output(Out),
test_rate_independence_invariants(Out, Id)
)
->
writeln('The model is rate independent.')
output(Out),
test_rate_independence_invariants(Out, Id),
fail
;
true
)
......@@ -488,88 +489,148 @@ test_rate_independence_invariants :-
).
test_rate_independence_invariants(Output, Id) :-
%! 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]),
(
is_covered_by_p_invariant_with_single_surinv(Output, Output, _)
is_persistent(Output),
p_invariant(P),
member(Output, P)
->
debug(rate, "~w is the single output of a P-invariant", [Output])
format("The output ~w is rate independent.", [Output])
;
format("Undecided, ~w not isolated in any P-invariant.", [Output]),
format("Undecided, ~w not persistent or not covered by any P-invariant.", [Output]),
fail
),
retractall(seen(_)),
retractall(visiting(_)),
iterate_on_parents(Output, Id).
is_covered_by_p_invariant_with_single_surinv(Species, Output, P) :-
p_invariant(P),
member(Species, P),
output(Output),
\+ (
sur_invariant(S),
S \= [Output],
ord_subset(S, 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) :-
assertz(seen(Species)),
debug(rate, "Searching for ~w parents", [Species]),
forall(
%! 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: edge, id: EdgeId1, item: (Reaction -> Species)]),
\+ (
item([parent: GraphId, kind: edge, id: EdgeId2, item: (Species -> Reaction)]),
debug(rate, "parent reaction ~w also has ~w as input", [Reaction, Species]),
reaction_graphs:get_stoichiometry(EdgeId1, Stoich1),
reaction_graphs:get_stoichiometry(EdgeId2, Stoich2),
debug(rate, "stoichiometries ~w -> ~w", [Stoich2, Stoich1]),
Stoich1 = Stoich2
)
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,
(
debug(rate, "parent reaction: ~w", [Reaction]),
retractall(visiting(Species-_)),
assertz(visiting(Species-[])),
forall(
item([parent: GraphId, kind: edge, item: (Reactant -> Reaction)]),
(
debug(rate, "parent: ~w", [Reactant]),
visiting(Species-ParentPinvs),
(
seen(Reactant)
->
debug(rate, "parent ~w already seen", [Reactant])
;
is_covered_by_sur_invariant(Reactant)
->
debug(rate, "parent ~w non-decreasing", [Reactant])
;
is_covered_by_p_invariant_with_single_surinv(Reactant, _, P),
\+ member(P, ParentPinvs)
->
retract(visiting(Species-ParentPinvs)),
assertz(visiting(Species-[P | ParentPinvs])),
debug(rate, "parent ~w in P-invariant, iterating", [Reactant]),
iterate_on_parents(Reactant, GraphId)
;
format("Undecided, ~w needed for an output and non-absorbed.", [Reactant]),
fail
)
)
)
% 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
).
normalize(Invar, Normalized) :-
maplist(nusmv:reactant_to_name, Invar, NoStoichiometry),
sort(NoStoichiometry, Normalized).
%! 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).
......@@ -29,28 +29,15 @@ test('Max', [setup(max_model), cleanup(clear_model)]) :-
with_output_to(atom(Result), test_rate_independence),
atom_concat('The model is rate independent', _, Result).
test('input_sinks', [
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))),
true(File-Output == File-Expected)
]) :-
with_output_to(atom(OutputFull), command(test_rate_independence_inputs_sinks)),
split_string(OutputFull, "\n", "\r\t ", OutputLines),
reverse(OutputLines, [_, OutputLong | _]),
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))),
setup((with_output_to(atom(ExpectedLong), command(load(File))), sub_string(ExpectedLong, 0, 3, _, Expected))),
true(File-Output == File-Expected)
]) :-
with_output_to(atom(OutputFull), command(test_rate_independence_invariants)),
split_string(OutputFull, "\n", "\r\t ", [OutputLong | _]),
sub_string(OutputLong, 0, 9, _, Output).
sub_string(OutputLong, 0, 3, _, Output).
:- end_tests(rate_independence).
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment