Commit 22cbdf64 authored by Mathieu Hemery's avatar Mathieu Hemery
Browse files

Merge branch 'develop' into feature/remove_fraction

parents 51084982 9cd83751
......@@ -6,7 +6,8 @@
additive_normal_form/2,
always_negative/1,
always_positive/1,
normalize_number/2
normalize_number/2,
is_null/1
]).
% Insert here for separate compilation and linting
......@@ -767,3 +768,23 @@ normalize_number(N, Norm) :-
;
Norm = N
).
%! is_null(+Expr)
%
% Check if an expression is uniformly null
is_null(Expr) :-
once(is_null_sr(Expr)).
is_null_sr(0).
is_null_sr(0.0).
is_null_sr(A*B) :- is_null_sr(A); is_null_sr(B).
is_null_sr(A/_B) :- is_null_sr(A).
is_null_sr(A+B) :- is_null_sr(A), is_null_sr(B).
is_null_sr(A-B) :- is_null_sr(A), is_null_sr(B).
is_null_sr(_A^(-_N)) :- !, false.
is_null_sr(A^_N) :- is_null_sr(A).
......@@ -5,4 +5,8 @@
test('distribute', [true(Out == a * a + a * b + c)]) :-
distribute(a * (a + b) + c, Out).
test('is_null', []) :-
is_null(a*0*1/2),
\+(is_null(a*b/0)).
:- end_tests(arithmetic_rules).
......@@ -695,13 +695,6 @@ binomial_bounded_search(First_Var, Varlist, Currentlist, Sumlist, PODE, Bound, R
% clean_writing(1*input, input).
% clean_writing(input^2*0, 0).
is_null(Zero) :-
catch(
Zero =:= 0,
error(_A,_B),
fail
).
is_one(One) :-
catch(
One =:= 1,
......
......@@ -30,6 +30,7 @@ commands = [
"check_conservations",
"check_ctl",
"check_ltl",
"check_model",
"check_multistability",
"check_oscillations",
"cleanup_ctl",
......
......@@ -35,6 +35,7 @@ Coder: T. Martinez
set_model_name/1,
inherits/1,
parametrize/0,
check_model/0,
% Public API
get_model_name/1,
get_model_name/2,
......@@ -72,13 +73,19 @@ Coder: T. Martinez
get_model/2,
load_all/2,
parametrize_initial_concentration/0,
parametrize_reaction/0
parametrize_reaction/0,
test_wf_reactant/2,
test_wf_inhibitor/2,
test_strictness/2
]
).
% Only for separate compilation/linting
:- use_module(arithmetic_rules).
:- use_module(doc).
:- use_module(formal_derivation).
:- use_module(reaction_rules).
:- use_module(util).
:- devdoc('This secton needs documentation for the developper, in particular about the invariants concerning item, parent, id, key, kind.
......@@ -297,6 +304,42 @@ parametrize :-
parametrize_initial_concentration,
parametrize_reaction.
check_model :-
biocham_command,
doc('Check if the current model is well formed and strict, (see TCS15_FGS).'),
doc('If it is not the case, proposes a witness.'),
single_model(ModelId),
findall(
Molecule,
identifier_kind(ModelId, Molecule, object),
List_Molecule
),
(
List_Molecule = []
->
format("There is currently no model to work on!~n", [])
;
reactions_with_species(List_Molecule, List_Reaction),
(
forall(
member(Reaction, List_Reaction),
test_wellformed(List_Molecule, Reaction)
)
->
format("The model is wellformed.~n", [])
;
format("The model is NOT wellformed.~n", [])
),
(
maplist(test_all_reactions(catalyst, test_strictness), List_Molecule)
->
format("The model is strict.~n", [])
;
format("The model is NOT strict.~n", [])
)
).
:- devdoc('\\section{Public API}').
......@@ -1261,3 +1304,91 @@ introduce_param(Func, Func2) :-
introduce_param(O,O).
%! test_all_reactions(+Type, +Test, +Molecule)
%
% apply test to all the reactions where Molecule is present as Type (catalyst, etc.)
% Type is concatenate with 'reactions_with_', cf reaction_rules.pl
test_all_reactions(Type, Test, Molecule) :-
atomic_concat('reactions_with_', Type, Apply),
Constructor =.. [Apply, [Molecule], List_Reactions],
call(Constructor),
forall(
member(Expression for _Reaction, List_Reactions),
(The_test =.. [Test, Molecule, Expression],
call(The_test))
).
%! test_wellformed(+List_Molecule, +Reaction)
%
% loop on molecule to determine if they follow the condition of wellformedness
test_wellformed([], _Reaction).
test_wellformed([Molecule|Tail], Expression for Reaction) :-
(
is_reactant(Molecule, Expression for Reaction)
->
test_wf_reactant(Molecule, Expression)
;
is_inhibitor(Molecule, Expression for Reaction)
->
test_wf_inhibitor(Molecule, Expression)
;
test_wf_other(Molecule, Expression)
),
test_wellformed(Tail, Expression for Reaction).
%! test_wf_reactant(+Molecule, +Expression)
%
% Check if the derivative of Expression wrt Molecule is positive somewhere
test_wf_reactant(Molecule, Expression) :-
formal_derivation:derivate(Expression, Molecule, Derivative),
(
\+(arithmetic_rules:always_negative(Derivative))
;
format("Fail for wellformedness of ~w with ~w as reactant.~n", [Expression, Molecule]),
fail
).
%! test_wf_inhibitor(+Molecule, +Expression)
%
% Check if the derivative of Expression wrt Molecule is negative somewhere
test_wf_inhibitor(Molecule, Expression) :-
formal_derivation:derivate(Expression, Molecule, Derivative),
(
\+(arithmetic_rules:always_positive(Derivative))
;
format("Fail for wellformedness of ~w with ~w as inhibitor.~n", [Expression, Molecule]),
fail
).
%! test_wf_other(+Molecule, +Expression)
test_wf_other(Molecule, Expression) :-
formal_derivation:derivate(Expression, Molecule, Derivative),
(
is_null(Derivative)
;
format("Fail for wellformedness of ~w with ~w as not present.~n", [Expression, Molecule]),
fail
).
%! test_strictness(+Molecule, +Expression)
%
% Check if Expression is null when Molecule is absent
test_strictness(Molecule, Expression) :-
substitute([Molecule], [0], Expression, Output),
(
is_null(Output)
;
format("Due to: ~w while checking ~w;~n", [Expression, Molecule]),
fail
).
......@@ -39,5 +39,16 @@ test(
list_model
).
test('test_wf_reactant', []):-
with_output_to(atom(_), (
test_wf_reactant(a, a*2),
\+(test_wf_reactant(m, -m))
)).
test('test_strictness', []):-
with_output_to(atom(_), (
test_strictness(a, a*b+a^2.5/c),
\+(test_strictness(mapk, j0*mk/(1+(mapk/j1)^j2)/(j1+mk)))
)).
:- end_tests(models).
......@@ -3,13 +3,7 @@
[
% Commands
test_rate_independence/0,
test_rate_independence_inputs_sinks/0,
test_rate_independence/2,
rate_independence_reduction_inputs_sinks/0,
rate_independence_reduction/2,
test_rate_independence_invariants/0,
%API
add_merge_reaction/3
test_rate_independence_invariants/0
]
).
......@@ -44,49 +38,7 @@ test_rate_independence :-
format("Time: ~p s~n", [TTime]).
% to generalize to influence models too
% Test rate independence for given inputs and outputs by performing rate-independence preserving reduction
test_rate_independence(Inputs, Outputs):-
biocham_command,
type(Inputs, {object}),
type(Outputs, {object}),
doc('Test graphical sufficient conditions for rate independence of the current model (assuming well-formed kinetics)
for the computation of molecular species in \\argument{Outputs} from species in \\argument{Inputs}.
Warning: this command destroys the current model by applying the \\command{rate_independence_reduction}.'),
rate_independence_reduction(Inputs, Outputs),
test_rate_independence.
test_rate_independence_inputs_sinks :-
biocham_command,
doc('Test graphical sufficient conditions for rate independence of the current model for the computation of graph output species from graph input species (assuming well-formed kinetics). The input species are the species that are not a reaction product or a strict catalyst (i.e. with same stoichiometry). The output species are the species that are not reactant.'),
input_species(Inputs),
sink_species(Outputs),
test_rate_independence(Inputs, Outputs).
% Performs a model reduction which preserves the rate independence property for given inputs and outputs
rate_independence_reduction(Inputs, Outputs):-
biocham_command,
type(Inputs, {object}),
type(Outputs, {object}),
doc('Reduces the current reaction model while preserving rate independence for computing \\argument{Outputs} species from \\argument{Inputs} species.'),
write('inputs: '), writeln(Inputs),
write('outputs: '), writeln(Outputs),
simplify(Inputs, Outputs),
!.
rate_independence_reduction_inputs_sinks:-
biocham_command,
doc('Reduces the current reaction model while preserving the rate independence property for computing the graph sink species from the input species. Warning: the reduction does not preserve the computation, only the rate independence property which may be true or false.'),
input_species(Inputs),
sink_species(Outputs),
write('sources: '), writeln(Inputs),
write('sinks: '), writeln(Outputs),
rate_independence_reduction(Inputs, Outputs).
% API
%%% API
%! is_synthesis_free is det.
%
......@@ -172,271 +124,6 @@ dfs_no_loop(GraphId, Vertex) :-
assertz(seen(Vertex)).
:- dynamic(removableInputs/0).
:- dynamic(restrictionOnTrivialReaction/0).
restrictionOnTrivialReaction.
:- dynamic(noIsolatedReactant/0).
noIsolatedReactant.
:- dynamic(noBetaEmpty/0).
noBetaEmpty. % incorrect otherwise on eliminating d in loop a+b<=>c+d=>e=>c+b
% Elimination of trivial reactions and isolated species with no outgoing fork
simplify(Inputs, Outputs):-
enumerate_molecules(Molecules),
(removableInputs -> Mol=Molecules ; subtract(Molecules, Inputs, Mol)),
subtract(Mol, Outputs, Removable),
reduce(Removable).
reduce(Molecules):-
list_model,
\+ trivial_reaction(Molecules),
\+ trivial_loop(Molecules),
\+ same_reaction,
\+ (member(M, Molecules), eliminate(M)).
reduce(Molecules):-
reduce(Molecules).
% Isolated reactant in one trivial reaction loop
trivial_loop(Molecules):-
member(M, Molecules),
reactions_with_product([M], [Reaction1]),
reactions_with_reactant([M], [Reaction2]),
reaction(Reaction1, _, Reactants, _, [C*M]),
reaction(Reaction2, _, [C*M], _, Reactants),
forall(member(_*X, Reactants), member(X, Molecules)), % not on output
write('--> eliminating '), write(M), write(' isolated in trivial loop: '), write(Reaction1), write(' <--> '), writeln(Reaction2),
delete_reaction(Reaction1),
delete_reaction(Reaction2).
% Elimination of trivial reactions with products=reactants apart from inputs and outputs
trivial_reaction(Molecules):-
item([kind:reaction, id: Id, item: Reaction]),
reaction(Reaction, _Kinetics, Reactants, _Inhibitors, Reactants),
(restrictionOnTrivialReaction -> forall(member((_*X), Reactants), member(X, Molecules)); true),
write('--> removing trivial reaction '), writeln(Reaction),
delete_item(Id).
% Elimination of doublon reactions with same reactants and products
% Todo: the case of equivalent stiochiometry?
same_reaction:-
item([kind:reaction, id: Id, item: Reaction]),
item([kind:reaction, id: Id2, item: Reaction2]),
Id\=Id2,
reaction(Reaction, _Kinetics, Reactants, _Inhibitors, Products),
reaction(Reaction2, _, Reactants, _, Products),
write('--> removing double reaction '), writeln(Reaction2),
delete_item(Id2).
% Isolated reactant in one reaction
% to do? eliminates d in c+d=>e=>c+output and does not conclude to rate independence
% or not to do? does not reduce c+d=>e=>c e=>output and does not conclude to rate independence
eliminate(M):-
\+ noIsolatedReactant,
reactions_with_product([M], []),
reactions_with_reactant([M], [Reaction1]),
reaction(Reaction1, Kinetics1, Reactants1, Inhibitors1, Products1),
delete(Reactants1, (_*M), Reac1),
Reac1\=[],
write('--> eliminating '), write(M), write(' isolated reactant in: '), writeln(Reaction1),
delete_reaction(Reaction1),
add_reaction(Kinetics1, Reac1, Inhibitors1, Products1, false).
% Isolated non increasing catalyst in one reaction
eliminate(M):-
reactions_with_reactant([M], [Reaction1]),
reactions_with_product([M], [Reaction1]),
reaction(Reaction1, Kinetics1, Reactants1, Inhibitors1, Products1),
member(C*M, Reactants1),
member(D*M, Products1),
C>=D,
write('--> eliminating '), write(M), write(' isolated catalyst in: '), writeln(Reaction1),
delete(Reactants1, (_*M), Reac1),
delete(Products1, (_*M), Prod1),
delete_reaction(Reaction1),
add_reaction(Kinetics1, Reac1, Inhibitors1, Prod1, false).
% Elimination of species M isolated in reactions alpha_i =>2 beta_i + M and M =>1 gamma
% if forall i alpha_i\cap\gamma=\emptyset
% by merging reactions resulting in alpha_i => beta_i + gamma
eliminate(M):-
reactions_with_reactant([M], [Reaction1]),
reaction(Reaction1, _, [C*M], _, _Products1), % gamma empty
reactions_with_product([M], Reactions),
\+ member(Reaction1, Reactions),
forall(
member(Reaction2, Reactions),
(
reaction(Reaction2, _, _Reactants2, _, Products2),
member(C*M, Products2)
%forall(member(D*N, Reactants2), \+ member(D*N, Products1)) % no creation of trivial loop (harmful on output)
)
),
write('--> eliminating '), write(M), write(' reactant only in '), write(Reaction1), write(' product of '), writeln(Reactions),
forall(
member(Reaction2, Reactions),
(
add_merge_reaction(M, Reaction1, Reaction2),
delete_reaction(Reaction2)
)
),
delete_reaction(Reaction1).
% Todo: accept reversible reactions with syntax <=>
% Elimination of M in alpha 2<=>1 M =>3 gamma rewritten to alpha=>gamma
% Elimination of M in alpha 2<=>1 beta+M =>3 gamma rewritten to alpha<=>beta =>gamma
eliminate(M):-
reactions_with_product([M], [Reaction1]),
reactions_with_reactant([M], [ReactionA, ReactionB]),
reaction(Reaction1, Kinetics1, Reactants1, Inhibitors1, Products1),
member((C*M), Products1),
(
reaction(ReactionA, _, _, _, Reactants1)
->
Reaction2= ReactionA,
Reaction3= ReactionB
;
reaction(ReactionB, _, _, _, Reactants1),
Reaction2= ReactionB,
Reaction3= ReactionA
),
Reaction1\=Reaction2,
reaction(Reaction2, Kinetics2, Reactants2, Inhibitors2, Products2),
member((C*M), Reactants2),
reaction(Reaction3, Kinetics3, Reactants3, Inhibitors3, Products3),
member((C*M), Reactants3),
delete(Reactants2, (C*M), Beta),
delete(Reactants3, (C*M), Beta),
delete(Products1, (C*M), Beta),
(noBetaEmpty -> Beta=[]; true),
write('--> eliminating '), write(M), write(' in loop '), write(Reaction1),
write(' <--> '),write(Reaction2), write(' --> '), writeln(Reaction3),
(
Beta=[]
-> % merge reactions
add_merge_reaction(M, Reaction1, Reaction3)
; % eliminate M from reactions
add_reaction(Kinetics2, Beta, Inhibitors2, Products2, false),
add_reaction(Kinetics3, Beta, Inhibitors3, Products3, false),
add_reaction(Kinetics1, Reactants1, Inhibitors1, Beta, false)
),
delete_reaction(Reaction2),
delete_reaction(Reaction3),
delete_reaction(Reaction1).
% Elimination of M in alpha 2<=>1 M 4<=>3 gamma rewritten to alpha<=>gamma
% Elimination of M in alpha 2<=>1 beta+M 4<=>3 gamma rewritten to alpha<=>beta <=>gamma
eliminate(M):-
reactions_with_product([M], [Reaction1, Reaction4]),
reactions_with_reactant([M], [ReactionA, ReactionB]),
reaction(Reaction1, Kinetics1, Reactants1, Inhibitors1, Products1),
member((C*M), Products1),
(
reaction(ReactionA, _, _, _, Products1)
->
Reaction2= ReactionA,
Reaction3= ReactionB
;
reaction(ReactionB, _, _, _, Products1),
Reaction2= ReactionB,
Reaction3= ReactionA
),
reaction(Reaction2, Kinetics2, Reactants2, Inhibitors2, Products2),
member((C*M), Reactants2),
reaction(Reaction3, Kinetics3, Reactants3, Inhibitors3, Products3),
member((C*M), Reactants3),
reaction(Reaction4, Kinetics4, Reactants4, Inhibitors4, Products4),
member((C*M), Products4),
delete(Reactants2, (C*M), Beta),
delete(Reactants3, (C*M), Beta),
delete(Products1, (C*M), Beta),
delete(Products4, (C*M), Beta),
(noBetaEmpty -> Beta=[]; true),
write('--> eliminating '), write(M), write(' in double loop '), write(Reaction1),
write(' <--> '),write(Reaction2), write(Reaction3), write(' <--> '), writeln(Reaction4),
(
Beta=[]
-> % merge reactions
add_merge_reaction(M, Reaction1, Reaction3),
add_merge_reaction(M, Reaction4, Reaction2)
; % eliminate M from reactions
add_reaction(Kinetics2, Beta, Inhibitors2, Products2, false),
add_reaction(Kinetics3, Beta, Inhibitors3, Products3, false),
add_reaction(Kinetics1, Reactants1, Inhibitors1, Beta, false),
add_reaction(Kinetics4, Reactants4, Inhibitors4, Beta, false)
),
delete_reaction(Reaction2),
delete_reaction(Reaction3),
delete_reaction(Reaction1),
delete_reaction(Reaction4).
:- devdoc(' eliminate M in two reactions by merging them and not removing them ').
add_merge_reaction(M, Reaction1, Reaction2):-
reaction(Reaction1, Kinetics1, Reactants1, Inhibitors1, Products1),
delete(Reactants1, (_*M), Reac1),
delete(Products1, (_*M), Prod1),
delete(Inhibitors1, M, Inhib1),
reaction(Reaction2, Kinetics2, Reactants2, Inhibitors2, Products2),
delete(Reactants2, (_*M), Reac2),
delete(Products2, (_*M), Prod2),
delete(Inhibitors2, M, Inhib2),
kinetics(Reac1, Inhib1, Kinetics1, Kin1), % M will remain in kinetics
kinetics(Reac2, Inhib2, Kinetics2, Kin2),
append(Reac1, Reac2, Reac),
append(Prod1, Prod2, Prod),
append(Inhib1, Inhib2, Inhib),
add_reaction(Kin1*Kin2,Reac,Inhib,Prod,false).
% Dynamic predicate for storing structural information of the graph
:- dynamic(p_invariant/1).
%% :- dynamic(sur_invariant/1).
......
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