Commit d52bfcfc authored by dcoudrin's avatar dcoudrin
Browse files

Merge branch 'develop' of...

Merge branch 'develop' of https://scm.gforge.inria.fr/authscm/dcoudrin/git/biocham/biocham into develop
parents b3b60f65 5807ee7a
......@@ -42,6 +42,7 @@ Rq: The distinction between pivp_string and pivp_list is not always obvious.
]).
:- use_module(doc).
:- use_module(biocham).
:- use_module(objects).
:- use_module(reaction_rules).
:- use_module(util).
......@@ -55,6 +56,8 @@ Rq: The distinction between pivp_string and pivp_list is not always obvious.
%%% Main Tools of the module %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
:- initial(option(binomial_reduction: no)).
%! compile_from_expression(+Expr, +Input, +Output)
%
% biocham command
......@@ -75,8 +78,14 @@ compile_from_expression(Expr, Input, Output) :-
creates a biochemical reaction network such that Output = Expr(Input). Use \'time\' as
input to get Output = Expr(t).
'),
option(
binomial_reduction,
yesno,
Reduction,
'Determine if the binomial reduction has to be performed'
),
expression_to_PIVP(Expr, PIVP),
pivp_compile(PIVP, Input, Output).
pivp_compile(PIVP, Input, Output, Reduction).
:- doc('
\\begin{example} Compilation of the expression 4+time^2:\n
......@@ -226,8 +235,14 @@ compile_from_pivp(PIVP, Input, Output) :-
compiles a PIVP into biochemical reactions. Use \'time\' as input
to get output f(t), or use any other input x to get f(x) as output.
'),
option(
binomial_reduction,
yesno,
Reduction,
'Determine if the binomial reduction has to be performed'
),
format_pivp(PIVP, P),
pivp_compile(P, Input, Output).
pivp_compile(P, Input, Output, Reduction).
:- devdoc('\\begin{todo} Compilation of the cosine function as a function of time defined as the solution of a PIVP)\n
\\texttt{
......@@ -265,13 +280,21 @@ test_2908 :-
').
%! pivp_compile(+PIVP_list, +Input, +Output)
%! pivp_compile(+PIVP_list, +Input, +Output, +Reduction)
%
% Compile a PIVP_list, this function make a call to compile_PIVP to
% Compile a PIVP_list, this function first reduce the PIVP to a binomial form
% then make a call to compile_PIVP to
% make the real work but take charge of the modification of the model
% when Input is not time through a call to g_to_c_PIVP/3.
pivp_compile(PIVP, Input, Output):-
pivp_compile(PIVP_non_bin, Input, Output, Reduction):-
(
Reduction = yes
->
reduce_to_binomial(PIVP_non_bin,PIVP)
;
PIVP = PIVP_non_bin
),
clear_model,
get_option(fast_rate, Fast),
set_parameter(fast, Fast),
......@@ -720,17 +743,36 @@ reduce_to_binomial([N,PODE,IV],[NewN,NewPODE,NewIV]) :-
Order > 2
->
highest_exponant_vector(PODE,HighestExp),
index_max(HighestExp,Max),
NewN is Max - 1,
generate_new_variable(HighestExp,NewN,ListVariable),
generate_pode(ListVariable,PODE,HighestExp,NewPODE),
generate_variable(HighestExp,NewN,ListVariable),
once(generate_pode(ListVariable,PODE,ListVariable,HighestExp,NewPODE)),
generate_iv(ListVariable,IV,NewIV)
; % PIVP is already binomial
NewN = N,
NewPODE = PODE,
NewIV = IV
).
%! generate_variable(+HighestExp,-NewN,-ListVariable)
%
% Generate all the variables in the new frame
generate_variable(HighestExp,NewN,ListVariable) :-
compute_number_species(HighestExp,NewSpeciesP1,UntouchedSpecies),
NewSpecies is NewSpeciesP1 - 1,
NewN is NewSpecies + UntouchedSpecies,
length(HighestExp,N),
generate_new_variable(HighestExp,NewSpecies,ListNewVariable),
generate_old_variable(ListOldVariable,N,0),
( % trick to keep the first variable (output) in head
[[1|_Tail1]|TruncListNewVariable] = ListNewVariable
->
union(ListOldVariable,TruncListNewVariable,ListVariable)
;
union(ListOldVariable,ListNewVariable,ListVariable)
).
%! generate_new_variable(+HighestExp,+NewN,-ListVariable)
%
......@@ -746,15 +788,24 @@ generate_new_variable(HighestExp,N,List) :-
append(ListTempo,[Var],List).
%! exponant_index(+OldExp,+HighestExp,-Index)
%! generate_old_variable(List,Length,Current)
%
% Compute the index of OldExp in the new variable set
% Generate the list of old indexes
exponant_index([],[],0).
generate_old_variable([],Length,Length) :- !.
exponant_index([HeadOld|TailOld],[HeadMax|TailMax],Index) :-
exponant_index(TailOld,TailMax,IndexTempo),
Index is (IndexTempo * (HeadMax + 1)) + HeadOld.
generate_old_variable([Head|Tail],Length,Current) :-
NC is Current + 1,
Right is Length - Current - 1,
displace_exponant([1],Current,Right,Head),
generate_old_variable(Tail,Length,NC).
%! exponant_index(?OldExp,+ListVariable,?Index)
%
% Compute the index of OldExp in the new list of variable set
exponant_index(OldExp,ListVariable,Index) :-
nth1(Index,ListVariable,OldExp).
%! index_exponant(+Index,+HighestExp,-OldExp)
......@@ -774,49 +825,54 @@ index_exponant_acc(Index,[Head|Tail],Acc,ExpList) :-
append([Rest],AccTempo,Acc).
%! index_max(+HighestExp,-Max)
%! compute_number_species(+HighestExp,Acc,AccZero).
%
% return the highest index in this set of variable +1
% subroutine to return the number of variables in the new set of variable
index_max([],1).
compute_number_species([],1,0).
index_max([Head|Tail],Max) :-
index_max(Tail,MaxTempo),
Max is MaxTempo * (Head+1).
compute_number_species([Head|Tail],Acc,AccZero) :-
(
Head =:= 0
->
compute_number_species(Tail,Acc,AccZeroTempo),
AccZero is AccZeroTempo+1
;
compute_number_species(Tail,AccTempo,AccZero),
Acc is AccTempo * (Head+1)
).
%! change_exponant(+OldExp,+HighestExp,-NewExp)
%! change_exponant(+OldExp,+ListVariable,-NewExp)
%
% Switch from the old variable set to the new one,
% HighestExp is needed for indexing
% Switch from the old variable set to the new one
change_exponant(_OldExp,[],[]).
change_exponant(OldExp,HighestExp,NewExp) :-
index_max(HighestExp,Max),
exponant_index(OldExp,HighestExp,Index),
Left is Index - 1,
Right is Max - Index - 1,
displace_exponant([1],Left,Right,NewExp).
change_exponant(OldExp,[OldExp|Tail],[1|T]) :- change_exponant(OldExp,Tail,T),!.
change_exponant(OldExp,[_Head|Tail],[0|T]) :- change_exponant(OldExp,Tail,T),!.
%! generate_pode(+ListVariable,+PODE,+HighestExp,-NewPODE)
%! generate_pode(+ListVariable,+PODE,+ListVariable,+HighestExp,-NewPODE)
%
% Generate the MultivarVector corresponding to PODE in the new set of variables
generate_pode([],_PODE,_HighestExp,[]) :- !.
generate_pode([],_PODE,_ListVariable,_HighestExp,[]) :- !.
generate_pode([Var|TailVar],PODE,HighestExp,[Multivar|Tail]) :-
generate_pode([Var|TailVar],PODE,ListVariable,HighestExp,[Multivar|Tail]) :-
sum_list(Var,VarTot),
(
VarTot =:= 1
->
% The variable already exist in the old set
extract_from_vector(Var,PODE,OldMultivar),
set_to_new_variable(OldMultivar,HighestExp,Multivar)
set_to_new_variable(OldMultivar,ListVariable,Multivar)
;
% This is a new variable
derive_new_variable([],Var,PODE,HighestExp,Multivar)
derive_new_variable([],Var,PODE,ListVariable,HighestExp,Multivar)
),
generate_pode(TailVar,PODE,HighestExp,Tail).
generate_pode(TailVar,PODE,ListVariable,HighestExp,Tail).
%! extract_from_vector(+Mask,+Vector,-Result)
......@@ -832,26 +888,26 @@ extract_from_vector([0|Tail],[_Skip|Vector],Element) :-
extract_from_vector(Tail,Vector,Element).
%! set_to_new_variable(+Multivar,+HighestExp,-NewMultivar)
%! set_to_new_variable(+Multivar,+ListVariable,-NewMultivar)
%
% Rephrase an old pODE in the new variable set
set_to_new_variable([],_,[]).
set_to_new_variable([[K,Exp]|Tail],HighExp,[[K,NewExp]|NewTail]) :-
change_exponant(Exp,HighExp,NewExp),
set_to_new_variable(Tail,HighExp,NewTail).
set_to_new_variable([[K,Exp]|Tail],ListVariable,[[K,NewExp]|NewTail]) :-
change_exponant(Exp,ListVariable,NewExp),
set_to_new_variable(Tail,ListVariable,NewTail).
%! derive_new_variable(Done,+Exponent,+MultivarVector,+HighestExp,-Multivar)
%! derive_new_variable(Done,+Exponent,+MultivarVector,+ListVariable,+HighestExp,-Multivar)
%
% Gives the derivative of an old exponent set in the new variable setting
% The MultivarVector is the core of the old PIVP, that is the derivative of
% the old (basic) variables.
derive_new_variable(_,[],[],_HighestExp,[]).
derive_new_variable(_,[],[],_ListVariable,_HighestExp,[]).
derive_new_variable(Done,[Exp|Tail],[Deriv|TailDeriv],HighestExp,Multivar) :-
derive_new_variable(Done,[Exp|Tail],[Deriv|TailDeriv],ListVariable,HighestExp,Multivar) :-
(
Exp >= 1
->
......@@ -859,30 +915,30 @@ derive_new_variable(Done,[Exp|Tail],[Deriv|TailDeriv],HighestExp,Multivar) :-
append(Done,[Expm],Exp_Tempo),
append(Exp_Tempo,Tail,ExpDeriv),
append(Done,[Exp],DoneE),
derive_new_variable(DoneE,Tail,TailDeriv,HighestExp,MultivarTempo),
multiply_derivative(ExpDeriv,Deriv,HighestExp,ThisDerivative),
derive_new_variable(DoneE,Tail,TailDeriv,ListVariable,HighestExp,MultivarTempo),
multiply_derivative(ExpDeriv,Deriv,ListVariable,HighestExp,ThisDerivative),
scalar_multivar(Exp,ThisDerivative,K_time_ThisDerivative),
append(K_time_ThisDerivative,MultivarTempo,Multivar)
;
Exp =:= 0
->
append(Done,[Exp],DoneE),
derive_new_variable(DoneE,Tail,TailDeriv,HighestExp,Multivar)
derive_new_variable(DoneE,Tail,TailDeriv,ListVariable,HighestExp,Multivar)
).
%! multiply_derivative(ExpSet,Deriv,HighestExp,ThisDerivative)
%! multiply_derivative(ExpSet,Deriv,ListVariable,HighestExp,ThisDerivative)
%
% Given an old ExpSet and a derivative of a basic variable
% return the corresponding Multivar in the new variable set
multiply_derivative(_ExpSet,[],_HighestExp,[]).
multiply_derivative(_ExpSet,[],_ListVariable,_HighestExp,[]).
multiply_derivative(ExpSet,[[K,Exp]|Deriv],HighestExp,[[K,NewExp]|ThisDerivative]) :-
change_exponant(ExpSet,HighestExp,NewExpSet),
change_exponant(Exp,HighestExp,NewExpDeriv),
multiply_derivative(ExpSet,[[K,Exp]|Deriv],ListVariable,HighestExp,[[K,NewExp]|ThisDerivative]) :-
change_exponant(ExpSet,ListVariable,NewExpSet),
change_exponant(Exp,ListVariable,NewExpDeriv),
add_list(NewExpSet,NewExpDeriv,NewExp),
multiply_derivative(ExpSet,Deriv,HighestExp,ThisDerivative).
multiply_derivative(ExpSet,Deriv,ListVariable,HighestExp,ThisDerivative).
%! generate_iv(+ListVariable,+IV,-NewIV)
......@@ -901,7 +957,7 @@ generate_iv([Var|TailVar],IV,[NewValue|TailValue]) :-
% determine the concentration of one variable (corresponding to Exponent)
% in the new set
set_new_iv([],[],1) :- !.
set_new_iv([],_IV,1) :- !.
set_new_iv([Exp|TailE],[Value|TailV],NewValue) :-
set_new_iv(TailE,TailV,NewValueTempo),
......@@ -966,6 +1022,10 @@ scan_order_multivar([[_K,Exponent]|Tail],Order) :-
max_list([],[],[]).
max_list([],List,List).
max_list(List,[],List).
max_list([Head1|Tail1],[Head2|Tail2],[HeadMax|TailMax]) :-
HeadMax is max(Head1,Head2),
max_list(Tail1,Tail2,TailMax).
......
......@@ -19,8 +19,15 @@ test(format_pivp_higher_order) :-
%%% Test of binomial reduction %%%
test(generate_new_variable) :-
gpac:generate_new_variable([3,1],7,[[1,0],[2,0],[3,0],[0,1],[1,1],[2,1],[3,1]]).
test() :-
gpac:compute_number_species([3,1],8,0),
gpac:compute_number_species([0,3,0,1,0],8,3),
gpac:compute_number_species([2,2,0],9,1).
test(generate_variable) :-
gpac:generate_variable([3,1],7,[[1,0],[2,0],[3,0],[0,1],[1,1],[2,1],[3,1]]),
gpac:generate_variable([2,2],8,[[1,0],[2,0],[0,1],[1,1],[2,1],[0,2],[1,2],[2,2]]),
gpac:generate_variable([0,2,0,1],7,[[1,0,0,0],[0,0,1,0],[0,1,0,0],[0,2,0,0],[0,0,0,1],[0,1,0,1],[0,2,0,1]]).
test(reduce_to_binomial_nomodif) :-
once(gpac:reduce_to_binomial([2,[[[1,[1, 1]]],[[-1,[1, 0]]]],[1,2]],P)),
......
......@@ -88,6 +88,7 @@ influence_graph(GraphId) :-
inhibitors: Inhibitors,
products: Products
]),
Origin = Item,
substract_list(Reactants, Products, Difference),
(
member(_ * Input, Reactants),
......@@ -112,6 +113,7 @@ influence_graph(GraphId) :-
sign: RuleSign,
output: Output
]),
Origin = Item,
(
member(Input, PositiveInputs),
InputSign = +
......@@ -139,6 +141,13 @@ influence_graph(GraphId) :-
)
;
set_attribute(EdgeId, sign = Sign)
),
(
get_attribute(EdgeId, origin = OldOrigin)
->
set_attribute(EdgeId, origin = [(Sign, Origin) | OldOrigin])
;
set_attribute(EdgeId, origin = [(Sign, Origin)])
)
)
),
......
/** <module> Influence graph export in Lemon graph format
Lemon is a [C++ graph library](http://lemon.cs.elte.hu/trac/lemon), lgf is its
input format
@author Sylvain Soliman
@license GPL
@copyright Inria EPI Lifeware 2018
*/
:- module(
lemon,
[
% Public API
export_lemon_graph/1
]
).
% :- use_module(doc).
:- use_module(counters).
:- use_module(invariants).
%! export_lemon_graph(+OutputFile) is det.
%
% exports the current influence graph to OutputFile (adding '.lgf' extension
% if needed)
export_lemon_graph(OutputFile) :-
biocham_command,
doc('exports the current influence graph to \\argument{OutputFile} in Lemon
graph format (\\url{http://lemon.cs.elte.hu/trac/lemon}) (adding \\texttt{.lgf}
extension if needed). Computes the conservation laws of the model in order
to do so.'),
type(OutputFile, output_file),
automatic_suffix(OutputFile, '.lgf', write, FilenameLgf),
setup_call_cleanup(
open(FilenameLgf, write, Stream),
export_lemon_stream(Stream),
close(Stream)
).
%! export_lemon_stream(+Stream) is det.
%
% exports the current influence graph to the given Stream
export_lemon_stream(Stream) :-
with_output_to(atom(_), invariants:find_pinvar),
setup_call_cleanup(
new_graph(GraphId),
(
influence_graph(GraphId),
export_lemon_stream(GraphId, Stream)
),
delete_item(GraphId)
).
%! species(-Species, -Uid) is nondet.
%
% stores a unique identifier Uid for each Species
:- dynamic(species/2).
%! export_lemon_stream(+Id, +Stream) is det.
%
% exports the graph associated to Id to the stream Stream
export_lemon_stream(Id, Stream) :-
write(Stream, '@nodes\nlabel\tspecies\n'),
set_counter(species, 0),
retractall(species(_, _)),
forall(
item([parent: Id, kind: vertex, item: VertexName]),
(
count(species, Count),
assertz(species(VertexName, Count)),
format(Stream, '~w\t~w\n', [Count, VertexName])
)
),
write(Stream, '@arcs\n\t\tlabel\treaction\tsign\n'),
number_reactions,
set_counter(influences, 0),
forall(
item([parent: Id, kind: edge, item: Edge, id: EdgeId]),
(
Edge = (VertexA -> VertexB),
species(VertexA, IdA),
species(VertexB, IdB),
get_attribute(EdgeId, origin = Origin),
forall(
member((Sign, Reaction), Origin),
(
count(influences, Count),
normalize_reaction(Reaction, Normalized, _Negated),
reactions(Normalized, RId),
sign(Sign, SignNumber),
format(
Stream,
'~w\t~w\t~w\t~w\t~w\n',
[IdA, IdB, Count, RId, SignNumber]
)
)
)
)
),
% TODO conservation laws
write(Stream, '@attributes\nconserv '),
findall(
Conserv,
(
invariants:base_mol(List),
atomic_list_concat(List, ',', Conserv)
),
ConservList
),
atomic_list_concat(ConservList,';', ConservAtom),
write(Stream, ConservAtom),
write(Stream, ';\n').
%! reactions(-Reaction, -Uid) is nondet.
%
% stores a unique identifier Uid for each Reaction
:- dynamic(reactions/2).
%! number_reactions is det.
%
% associates a unique number to each reaction
% such that reaction that are stoichiometrically inverse get opposite numbers
% and stoichiometrically equivalent get same number
number_reactions :-
retractall(reactions(_, _)),
set_counter(reactions, 1),
forall(
item([kind: reaction, item: Item]),
(
normalize_reaction(Item, Normalized, Negated),
(
reactions(Normalized, _Count)
->
true
;
reactions(Negated, Index)
->
Count is -Index,
assertz(reactions(Normalized, Count))
;
count(reactions, Count),
assertz(reactions(Normalized, Count))
)
)
).
%! sign(+Atom, -Number) is det.
%
% associates number Number to a sign Atom
sign('+', 1).
sign('-', -1).
%! normalize_reaction(+Item, -NormalizedStoichiometry, -Negated) is det.
%
% compute a sorted stoichiometry list for the reaction Item as
% NormalizedStoichiometry and a sorted list for the inverse reaction as
% Negated
normalize_reaction(Item, NormalizedStoichiometry, Negated) :-
reaction_editor:get_stoichiometry_and_kinetics(Item, Stoichiometry, _),
sort(Stoichiometry, NormalizedStoichiometry),
maplist(reaction_editor:negate_coefficient, Stoichiometry, NegatedStoich),
sort(NegatedStoich, Negated).
......@@ -600,15 +600,9 @@ make_reaction(
compute_ode_for_reactions :-
forall(
item([kind: reaction, item: Item]),
(
item([kind: reaction, item: Item]),
reaction(Item, Kinetics, Reactants, Inhibitors, Products)
),
(
kinetics(Reactants, Inhibitors, Kinetics, KineticsExpression),
maplist(negate_coefficient, Reactants, NegatedReactants),
append(NegatedReactants, Products, StoichiometricSolution),
simplify_solution(StoichiometricSolution, SimplifiedSolution),
get_stoichiometry_and_kinetics(Item, SimplifiedSolution, KineticsExpression),
add_molecules(SimplifiedSolution, KineticsExpression)
)
).
......@@ -627,4 +621,19 @@ add_molecule(Coefficient * Molecule, Kinetics) :-
ode_add_expression_to_molecule(Coefficient * Kinetics, Molecule).
negate_coefficient(C*X, (-C)*X).
negate_coefficient(C*X, D*X) :-
(
number(C)
->
D is -C
;
D = -C
).
get_stoichiometry_and_kinetics(Item, SimplifiedSolution, KineticsExpression) :-
reaction(Item, Kinetics, Reactants, Inhibitors, Products),
kinetics(Reactants, Inhibitors, Kinetics, KineticsExpression),
maplist(negate_coefficient, Reactants, NegatedReactants),
append(NegatedReactants, Products, StoichiometricSolution),
simplify_solution(StoichiometricSolution, SimplifiedSolution).
......@@ -45,6 +45,7 @@
- influence_editor.pl
** Influence graph
- influence_graphs.pl
- lemon.pl
* Events
- events.pl
* Importing and Exporting BIOCHAM Models
......
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