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

Ajout de test d'independance des reactions sur entrees/sorties par reecriture...

Ajout de test d'independance des reactions sur entrees/sorties par reecriture de graphe preservant l'independance des reactions et ajout de calcul des sources, puits et pathways dans graphe de reaction
parent 587b8803
......@@ -131,8 +131,15 @@ commands = [
"list_options",
"list_parameters",
"list_reactions",
"list_reactions_with_catalyst",
"list_reactions_with_inhibitor",
"list_reactions_with_product",
"list_reactions_with_reactant",
"list_reactions_with_species",
"list_rows",
"list_rules",
"list_sink_species",
"list_source_species",
"list_stable_states",
"list_tables",
"list_tscc_candidates",
......@@ -163,12 +170,14 @@ commands = [
"pac_learning",
"parameter",
"parametrize",
"pathway",
"pattern_reduction",
"place",
"plot",
"present",
"prolog",
"quit",
"rate_independence_reduction",
"reaction_graph",
"reaction_model",
"reduce_model",
......
:- module(
graph_properties,
[
pathway/2,
list_source_species/0,
list_sink_species/0,
% API
pathway/3,
source_species/1,
sink_species/1,
is_source/2,
is_sink/2
]
).
:- doc('
Analysis of some graph properties of the current model.
').
pathway(Object1, Object2):-
biocham_command,
type(Object1, object),
type(Object2, object),
doc('Gives one reaction pathway from \\argument{Object1} to \\argument{Object2} if one exists in the directed reaction graph of the current model (for more complex queries, see next section on Computation Tree Logic model-checking).'),
pathway(Object1, Object2, Path),
reverse(Path, Reactions),
forall(member(R, Reactions), writeln(R)).
pathway(Source, Target, Path):-
path([(Source, Source)], [], [(Source, Source)], Target, Path).
% Dijsktra algorithm with pathway memory
path([], [], _Forward, _Target, _Path):-
!.
path([], New, Forward, Target, Path):-
path(New, [], Forward, Target, Path).
path([A|Frontier], New, Forward, Target, Path):-
new_successors(A, Forward, List),
(
member((Target, Path), List)
->
true
;
append(List, New, New2),
append(List, Forward, Forward2),
path(Frontier, New2, Forward2, Target, Path)
).
new_successors(A, Forward, List):-
realsetof(B, graph_properties:new_successor(A, B, Forward), List).
new_successor((A, PathA), (B, [Reaction|PathA]), Forward):-
is_reactant(A, Reaction),
is_product(B, Reaction),
\+ member((B,_), Forward).
% Sources and Sinks
list_source_species:-
biocham_command,
doc('Lists the molecular species that are neither a reaction product, nor the target of an influence rule, in the curent model.'),
source_species(Molecules),
writeln(Molecules).
list_sink_species:-
biocham_command,
doc('Lists the molecular species that are neither a reactant, nor a (positive) source of an influence rule, in the curent model.'),
sink_species(Molecules),
writeln(Molecules).
%API
:-devdoc('Prolog lists of graph sources and sinks.').
source_species(L):-
enumerate_molecules(Molecules),
realsetof(M, is_source(M, Molecules), L).
sink_species(L):-
enumerate_molecules(Molecules),
realsetof(M, is_sink(M, Molecules), L).
is_source(M, Molecules):-
member(M, Molecules),
\+ (
is_product(M,_)
),
\+ (
is_target(M,_)
).
is_sink(M, Molecules):-
member(M, Molecules),
\+ (
is_reactant(M,_)
),
\+ (
is_source(M,_)
).
\ No newline at end of file
......@@ -15,7 +15,10 @@
influence/2,
compute_ode_for_influences/0,
print_influence/2,
substract_list/3
substract_list/3,
is_influence_source/2,
is_negative_source/2,
is_target/2
]
).
......@@ -451,3 +454,21 @@ inputs(
list_enumeration(PositiveInputList, PositiveInputEnum),
list_enumeration(NegativeInputList, NegativeInputEnum),
!.
is_influence_source(S, Influence):-
item([kind: influence, item: Influence]),
influence(Influence, [force: _, positive_inputs: Sources, negative_inputs:_, sign;_, output:_]),
member((_*S), Sources).
is_negative_source(S, Influence):-
item([kind: influence, item: Influence]),
influence(Influence, [force: _, positive_inputs: _, negative_inputs:Sources, sign;_, output:_]),
member((_*S), Sources).
is_target(S, Influence):-
item([kind: influence, item: Influence]),
influence(Influence, [force: _, positive_inputs: _, negative_inputs:_, sign;_, output: S]).
......@@ -2,7 +2,9 @@
rate_independence,
[
% Commands
test_rate_independence/0
test_rate_independence/0,
test_rate_independence/2,
rate_independence_reduction/2
]
).
......@@ -15,7 +17,7 @@
% Test the 3 graphical conditions (and not yet well-formedness) for rate independence.
test_rate_independence :-
biocham_command,
doc('Test the graphical sufficient conditions for rate independence of the current model.'),
doc('Test graphical sufficient conditions for rate independence of the current model for all output species.'),
reaction_graph,
get_current_graph(Id),
(
......@@ -32,6 +34,29 @@ test_rate_independence :-
).
% 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('Checks whether the current reaction model is rate independent
for computing the given set of molecular species \\argument{Outputs} from \\argument{Inputs}.
Warning: this command destroys the current model by applying the \\command{rate_independence_reduction}.'),
rate_independence_check(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 indepence for computing \\argument{Outputs} species from \\argument{Inputs} species.'),
simplify(Inputs, Outputs).
% API
%! is_synthesis_free is det.
%
% check for the absence of a reaction with reactants included in products
......@@ -114,3 +139,129 @@ dfs_no_loop(GraphId, Vertex) :-
)
),
assertz(seen(Vertex)).
% Graph rewriting for checking rate independence with respect to Inputs and Outputs
rate_independence_check(Inputs, Outputs):-
simplify(Inputs, Outputs),
check_rate_independence.
% Compression of species with no outgoing fork and trivial loops
simplify(Inputs, Outputs):-
enumerate_molecules(L),
subtract(L, Inputs, L1),
subtract(L1, Outputs, Molecules),
\+ simplify(Molecules),
\+ trivial_reaction,
\+ same_reaction.
simplify(Inputs, Outputs):-
simplify(Inputs, Outputs).
% Elimination of trivial reactions with products same as reactants
trivial_reaction:-
item([kind:reaction, id: Id, item: Reaction]),
reaction(Reaction, _Kinetics, Reactants, _Inhibitors, Reactants),
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]),
reaction(Reaction, _Kinetics, Reactants, _Inhibitors, Products),
item([kind:reaction, id: Id2, item: Reaction2]),
Id\=Id2,
reaction(Reaction2, _, Reactants, _, Products),
delete_item(Id2).
% Merge of reactions on molecules with no outgoing fork or in trivial loops
simplify(Molecules):-
member(M, Molecules),
simpl(M),
!.
% Elimination of species M if no outgoing fork rewritten by merging reactions
% alpha_i => beta_i + M and M + gamma => delta rewritten alpha_i + gamma => beta_i + delta
simpl(M):-
reactions_with_reactant([M], [Reaction1]),
reactions_with_product([M], Reactions),
subtract(Reactions, [Reaction1], Reactions2),
Reactions2\=[],
!,
delete_reaction(Reaction1),
member(Reaction2, Reactions2),
delete_reaction(Reaction2),
add_merge_reactions(M, Reaction1, Reaction2).
% Todo: accept reversible reactions with syntax <=>
% Elimination of M in 2 loops alpha 1<=>3 M (4<)=>2 gamma rewritten alpha (<)=> gamma
simpl(M):-
reactions_with_reactant([M], [Reaction1, Reaction2]),
reactions_with_product([M], [Reaction3 | Reactions]),
(
Reactions=[Reaction4]
->
reaction(Reaction4, _, Reactants4, _, [S*M])
;
Reactions=[]
),
reaction(Reaction1, _, [S*M], _, Products1),
reaction(Reaction2, _, [S*M], _, Products2),
reaction(Reaction3, _, Reactants3, _, [S*M]),
(
Products1=Reactants3
->
add_merge_reactions(M, Reaction3, Reaction2),
(
Reactions=[]
->
true
;
Products2=Reactants4,
add_merge_reactions(M, Reaction4, Reaction1),
delete_reaction(Reaction4)
)
;
Products2=Reactants3,
add_merge_reactions(M, Reaction3, Reaction1),
(
Reactions=[]
->
true
;
Products1=Reactants4,
add_merge_reactions(M, Reaction4, Reaction2),
delete_reaction(Reaction4)
)
),
delete_reaction(Reaction1),
delete_reaction(Reaction2),
delete_reaction(Reaction3).
add_merge_reactions(M, Reaction1, Reaction2):-
reaction(Reaction1, Kinetics1, Reactants1, Inhibitors1, Products1),
subtract(Reactants1, [_*M], Reac1),
subtract(Products1, [_*M], Prod1),
subtract(Inhibitors1, [M], Inhib1),
reaction(Reaction2, Kinetics2, Reactants2, Inhibitors2, Products2),
subtract(Reactants2, [_*M], Reac2),
subtract(Products2, [_*M], Prod2),
subtract(Inhibitors2, [M], Inhib2),
kinetics(Reac1, Inhib1, Kinetics1, Kin1),
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).
\ No newline at end of file
......@@ -4,6 +4,11 @@
% Commands
add_reaction/1,
list_reactions/0,
list_reactions_with_species/1,
list_reactions_with_reactant/1,
list_reactions_with_product/1,
list_reactions_with_catalyst/1,
list_reactions_with_inhibitor/1,
list_rules/0,
delete_reaction/1,
delete_reactions/0,
......@@ -11,12 +16,24 @@
delete_reaction_prefixed/1,
merge_reactions/2,
simplify_all_reactions/0,
% Public API
% Public API
reactions_with_species/2,
reactions_with_reactant/2,
reactions_with_product/2,
reactions_with_catalyst/2,
reactions_with_inhibitor/2,
is_reactant/2,
is_product/2,
is_catalyst/2,
is_inhibitor/2,
find_reaction_prefixed/2,
is_reaction_model/0,
check_reaction_model/0,
list_model_reactions/0,
reaction/2,
reaction/5,
add_reaction/5,
add_reaction/6,
solution_to_list/2,
list_to_solution/2,
compute_ode_for_reactions/0
......@@ -84,7 +101,7 @@ delete_reaction_prefixed(Prefix):-
biocham_command,
type(Prefix, name),
doc('
removes the reaction rules that match a name prefix.
removes all the reaction rules that match a name prefix.
'),
Item=(Name--_),
\+ (
......@@ -139,9 +156,113 @@ list_rules :-
list_items([kind: [reaction, influence]]).
list_reactions_with_species(ObjectSet):-
biocham_command,
type(ObjectSet, {object}),
doc('Lists reactions having a reactant, product or inhibitor in the set of molecular species \\argument{Objectset}.'),
reactions_with_species(ObjectSet, Reactions),
forall(member(R,Reactions), writeln(R)).
reactions_with_species(ObjectSet, Reactions):-
realsetof(Reaction, reaction_editor:species(ObjectSet, Reaction), Reactions).
species(ObjectSet, Reaction):-
member(M,ObjectSet),
(is_reactant(M, Reaction); is_product(M, Reaction); is_catalyst(M, Reaction); is_inhibitor(M, Reaction)).
list_reactions_with_reactant(ObjectSet):-
biocham_command,
type(ObjectSet, {object}),
doc('Lists reactions having a reactant in the set of molecular species \\argument{Objectset}.'),
reactions_with_reactant(ObjectSet, Reactions),
forall(member(R,Reactions), writeln(R)).
% setof with composite goals is bugged !
reactions_with_reactant(ObjectSet, Reactions):-
realsetof(Reaction, reaction_editor:reactant(ObjectSet, Reaction), Reactions).
reactant(ObjectSet, Reaction):-
member(M, ObjectSet),
is_reactant(M, Reaction).
list_reactions_with_product(ObjectSet):-
biocham_command,
type(ObjectSet, {object}),
doc('Lists reactions having a product in the set of molecular species \\argument{Objectset}.'),
reactions_with_product(ObjectSet, Reactions),
forall(member(R,Reactions), writeln(R)).
reactions_with_product(ObjectSet, Reactions):-
realsetof(Reaction, reaction_editor:product(ObjectSet, Reaction), Reactions).
product(ObjectSet, Reaction):-
member(M,ObjectSet),
is_product(M, Reaction).
list_reactions_with_catalyst(ObjectSet):-
biocham_command,
type(ObjectSet, {object}),
doc('Lists reactions with a catalyst in the set of molecular species \\argument{Objectset}.'),
reactions_with_catalyst(ObjectSet, Reactions),
forall(member(R,Reactions), writeln(R)).
reactions_with_catalyst(ObjectSet, Reactions):-
realsetof(Reaction, reaction_editor:catalyst(ObjectSet, Reaction), Reactions).
catalyst(ObjectSet, Reaction):-
member(M,ObjectSet),
is_catalyst(M, Reaction).
list_reactions_with_inhibitor(ObjectSet):-
biocham_command,
type(ObjectSet, {object}),
doc('Lists reactions with an inhibitor in the set of molecular species \\argument{Objectset}.'),
reactions_with_inhibitor(ObjectSet, Reactions),
forall(member(R,Reactions), writeln(R)).
reactions_with_inhibitor(ObjectSet, Reactions):-
realsetof(Reaction, reaction_editor:inhibitor(ObjectSet, Reaction), Reactions).
inhibitor(ObjectSet, Reaction):-
member(M,ObjectSet),
is_inhibitor(M, Reaction).
:- devdoc('\\section{Public API}').
is_reactant(S, Reaction):-
item([kind: reaction, item: Reaction]),
reaction(Reaction, [kinetics: _,reactants: Reactants,inhibitors: _, products: _]),
member((_*S), Reactants).
is_product(S, Reaction):-
item([kind: reaction, item: Reaction]),
reaction(Reaction, [kinetics: _,reactants: _,inhibitors: _, products: Products]),
member((_*S), Products).
is_catalyst(S, Reaction):-
item([kind: reaction, item: Reaction]),
reaction(Reaction, [kinetics: _,reactants: Reactants,inhibitors: _, products: Products]),
member((_*S), Reactants),
member((_*S), Products).
is_inhibitor(S, Reaction):-
item([kind: reaction, item: Reaction]),
reaction(Reaction, [kinetics: _,reactants: _,inhibitors: Inhibitors, products: _]),
member(S, Inhibitors).
is_reaction_model :-
devdoc('
succeeds if the current model is a reaction model
......@@ -778,3 +899,4 @@ get_stoichiometry_and_kinetics(Item, SimplifiedSolution, KineticsExpression) :-
maplist(negate_coefficient, Reactants, NegatedReactants),
append(NegatedReactants, Products, StoichiometricSolution),
simplify_solution(StoichiometricSolution, SimplifiedSolution).
......@@ -61,6 +61,8 @@
#** Other files
* Part II: Qualitative Analysis and Synthesis
* Static Analyses
** Graph properties
- graph_properties.pl
** Dimensional analysis
- units.pl
** Conservation laws and invariants
......
:- module(
util,
[
add_list/3,
multiset_inclusion/2,
add_list/3,
multiply_list/3,
name_variables_and_anonymous/2,
camel_case/2,
......@@ -55,12 +56,21 @@
write_condition_to_c/2,
write_condition_to_c_with_parentheses/2,
normal_random/1,
atom_to_int/2
atom_to_int/2,
realsetof/3
]).
:- use_module(library(process)).
:- use_module(doc).
multiset_inclusion([],_).
multiset_inclusion([(S*M)|L], Multiset):-
member((T*M), Multiset),
S=<T,
multiset_inclusion(L, Multiset).
name_variables(L) :-
doc('
instantiates a list of bindings returned by \\texttt{read_term/3}\'s
......@@ -908,3 +918,13 @@ normal_random(X) :-
atom_to_int(String, Int) :-
read_term_from_atom(String, Int, []),
integer(Int).
% what setof should have always been
% still setof with composite goals is bugged !
realsetof(X, Goal, List):-
setof(X, Goal, List),
!.
realsetof(_, _, []).
Supports Markdown
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