Commit b0a9576c authored by SOLIMAN Sylvain's avatar SOLIMAN Sylvain
Browse files

Invariant-based rate-indep first version

parent 3fedc0ef
......@@ -7,6 +7,7 @@
test_rate_independence/2,
rate_independence_reduction_inputs_sinks/0,
rate_independence_reduction/2,
test_rate_independence_invariants/0,
%API
add_merge_reaction/3
]
......@@ -415,4 +416,139 @@ add_merge_reaction(M, Reaction1, Reaction2):-
append(Inhib1, Inhib2, Inhib),
add_reaction(Kin1*Kin2,Reac,Inhib,Prod,false).
\ No newline at end of file
:- dynamic(sur_invariant/1).
:- dynamic(p_invariant/1).
:- dynamic(output/1).
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(_)),
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]
->
debug(rate, "output: ~w", [B]),
assertz(output(Out))
;
true
)
)
),
with_output_to(atom(_), invariants:find_invar_aux(4, is_place, is_transition, '#=')),
forall(
retract(invariants:base_mol(B)),
(
normalize(B, BB),
debug(rate, "p_invariant: ~w", [BB]),
assertz(p_invariant(BB))
)
),
reaction_graph,
get_current_graph(Id),
(
output(_)
->
(
forall(
output(Out),
test_rate_independence_invariants(Out, Id)
)
->
writeln('The model is rate independent.')
;
true
)
;
write('Undecided: no outputs!')
).
test_rate_independence_invariants(Output, Id) :-
debug(rate, "Trying for output ~w", [Output]),
(
is_covered_by_p_invariant_with_single_surinv(Output, Output)
->
debug(rate, "~w is the single output of a P-invariant", [Output])
;
format("Undecided, ~w not isolated in any P-invariant.", [Output]),
fail
),
retractall(seen(_)),
iterate_on_parents(Output, Id).
is_covered_by_p_invariant_with_single_surinv(Species, Output) :-
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),
!.
iterate_on_parents(Species, GraphId) :-
assertz(seen(Species)),
debug(rate, "Searching for ~w parents", [Species]),
forall(
(
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
)
),
(
debug(rate, "parent reaction: ~w", [Reaction]),
forall(
item([parent: GraphId, kind: edge, item: (Reactant -> Reaction)]),
(
debug(rate, "parent: ~w", [Reactant]),
(
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, _)
->
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
)
)
)
)
).
normalize(Invar, Normalized) :-
maplist(nusmv:reactant_to_name, Invar, NoStoichiometry),
sort(NoStoichiometry, Normalized).
......@@ -29,19 +29,29 @@ 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('everything', [
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(_), command(search_conservations)),
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))),
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).
:- end_tests(rate_independence).
max_model :-
......
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