Commit 2e43869d authored by Mathieu Hemery's avatar Mathieu Hemery

Merge branch 'feature/well_formed' into develop

parents fb10a043 b821e8b4
......@@ -11,6 +11,7 @@
]
).
:- use_module(doc).
:- devdoc('\\section{Commands}').
......
......@@ -6,6 +6,7 @@
additive_normal_form/2,
always_negative/1,
always_positive/1,
determine_sign/2,
normalize_number/2,
is_null/1,
modulo/3
......@@ -665,105 +666,132 @@ arithmetic_operation(- _).
arithmetic_operation(_ ^ _).
always_negative(A) :-
rewrite(rewrite_simplify, A, B),
B \= A,
!,
always_negative(B).
always_negative(A) :-
number(A),
!,
A =< 0.
%! always_negative(+Expr)
%! always_positive(+Expr)
%
% Check if Expr is assure to be always negative/positive, fail if it can change sign
% or if it is impossible to decide
always_negative(- A) :-
always_positive(A).
always_negative(Expr) :-
determine_sign(Expr, zero);
determine_sign(Expr, neg).
always_negative(A + B) :-
always_negative(A),
always_negative(B).
always_positive(Expr) :-
determine_sign(Expr, zero);
determine_sign(Expr, pos).
always_negative(A - B) :-
always_negative(A),
always_positive(B).
always_negative(A * B) :-
always_negative(A),
always_positive(B).
%! determine_sign(+Expr, -Sign)
%
% Try to infer the sign of Expr, Sign will be unified to: neg/pos/zero/dnk (do not know)
always_negative(0 * _) :-
!.
% TODO : Verify if it needs rewriting
%always_negative(A) :-
% rewrite(rewrite_simplify, A, B),
% B \= A,
% !,
% always_negative(B).
always_negative(A * B) :-
always_positive(A),
always_negative(B).
determine_sign(A, zero) :-
is_null(A),
!.
always_negative(A / B) :-
always_negative(A),
always_positive(B).
determine_sign(A, Sign) :-
number(A),
!,
(
A = 0
->
Sign = zero
;
A < 0
->
Sign = neg
;
Sign = pos
).
always_negative(A / B) :-
always_positive(A),
always_negative(B).
determine_sign(- A, NSign) :-
determine_sign(A, Sign),
!,
(
Sign = pos
->
NSign = neg
;
Sign = neg
->
NSign = pos
;
NSign = Sign
).
determine_sign(A + B, Sign) :-
!,
(
determine_sign(A, SignAB),
determine_sign(B, SignAB)
->
Sign = SignAB
;
Sign = dnk
).
always_positive(A) :-
rewrite(rewrite_simplify, A, B),
B \= A,
determine_sign(A - B, Sign) :-
!,
always_positive(B).
determine_sign(A + (- B), Sign).
always_positive(A) :-
number(A),
determine_sign(A * B, Sign) :-
!,
A >= 0.
determine_sign(A, SignA),
determine_sign(B, SignB),
(
(SignA = zero ; SignB = zero)
->
Sign = zero
;
(SignA = dnk ; SignB = dnk)
->
Sign = dnk
;
SignA = SignB
->
Sign = pos
;
Sign = neg
).
always_positive(A):- % assumed for concentrations, parameters
atom(A),
!.
determine_sign(A / B, Sign) :-
!,
(
is_null(B)
->
throw(error(division_by_zero, arithmetic_rules_determine_sign))
;
true
),
determine_sign(A * B, Sign).
always_positive(- A) :-
always_negative(A).
determine_sign(A, pos) :- % concentrations and parameters are assumed positive
atom(A),
!.
always_positive(_A ^ B) :-
determine_sign(_A ^ B, pos) :-
number(B),
(
modulo(B, 2, 0)
->
!,
true
).
always_positive(A ^ _B) :-
always_positive(A).
always_positive(A + B) :-
always_positive(A),
always_positive(B).
always_positive(A - B) :-
always_positive(A),
always_negative(B).
always_positive(0 * _) :-
!.
always_positive(A * B) :-
always_negative(A),
always_negative(B).
always_positive(A * B) :-
always_positive(A),
always_positive(B).
determine_sign(A ^ _B, Sign) :-
determine_sign(A, Sign).
always_positive(A / B) :-
always_positive(A),
always_positive(B).
determine_sign(exp(_A), pos).
always_positive(A / B) :-
always_negative(A),
always_negative(B).
always_positive(exp(_A)).
normalize_number(N, Norm) :-
(
......
......@@ -9,6 +9,14 @@ test('distribute2', []) :-
distribute((a+b)^2, a*a+a*b+(b*a+b*b)),
distribute((a+b)^ -2, (a+b)^ -2).
test('determine_sign', []) :-
determine_sign((a+1)^2, pos),
determine_sign((a-1)^3, dnk),
determine_sign(-exp(a), neg),
determine_sign(a/(a+1), pos),
determine_sign((2+a)*0, zero).
test('is_null', []) :-
is_null(a*0*1/2),
\+(is_null(a*b/0)).
......
......@@ -101,6 +101,7 @@ commands = [
"hybrid_static_simulation",
"import_ode",
"import_reactions_from_graph",
"infer_hidden_molecules",
"influence_graph",
"influence_hypergraph",
"influence_model",
......
#!/usr/bin/bash
# quadratization
swipl --goal=main --stand_alone=true -DSWIPL_LINKER_FLAGS="-static" -DSWIPL_SHARED_LIB=OFF\
--foreign=save -o quadratization -c biochamlib_quadratization.pl
# ODE2SBML
swipl --goal=main --stand_alone=true -DSWIPL_LINKER_FLAGS="-static" -DSWIPL_SHARED_LIB=OFF\
--foreign=save -o ODE2SBML -c biochamlib_ode2sbml.pl
#!/usr/bin/prolog
/*
* The prolog script to create the stand-alone saved state for conversion of an ODE to an
* SBML file.
* Coder: M. Hemery
*/
% Loading libraries
:- use_module(aliases).
:- use_module(arithmetic_rules).
:- use_module(biocham).
:- use_module(counters).
:- use_module(doc).
:- use_module(filename).
:- use_module(formal_derivation).
:- use_module(functions).
:- use_module(initial_state).
:- use_module(kinetics).
:- use_module(models).
:- use_module(molecules).
:- use_module(namespace).
:- use_module(objects).
:- use_module(ode).
:- use_module(parameters).
:- use_module(reaction_editor).
:- use_module(reaction_rules).
:- use_module(sbml_files).
:- use_module(toplevel).
:- use_module(types).
:- use_module(util).
:- use_module(xpp_parser).
set_my_option(Option, Value) :-
change_item([], option, Option, option(Option: Value)).
% parse_arguments(+Arguments, -Input)
%
% Handle the various option of the standalone:
% -i: input file
parse_arguments([], none).
parse_arguments(['-i', Input| Argv], Input) :-
parse_arguments(Argv, _I).
parse_arguments(['--help'| _Argv], _I, _M, _T) :-
print_help.
% print_help
print_help :-
format("help~n", []),
halt.
% main(+argv)
main(Argv) :-
(Argv = [] -> print_help ; true),
% Handling naming and options
parse_arguments(Argv, Input),
file_name_extension(InputBase, 'ode', Input),
file_name_extension(InputBase, 'sbml', Output),
format("Input: ~w~nOutput: ~w~n", [Input, Output]),
% Initialization (see biocham:initialize)
set_prolog_flag(allow_variable_name_as_functor, true),
set_counter(list_item_counter, 0),
set_counter(item_id, 0),
set_counter(model, 0),
set_counter(molecule_id, 0),
set_counter(parameter_id, 1),
nb_setval(ode_viewer, inline),
nb_setval(current_models, []),
new_model,
set_my_option(import_reactions_with_inhibitors, yes),
% Main procedure
read_xpp(Input),
list_ode,
add_reactions_from_ode_system,
export_sbml(Output).
#!/usr/bin/prolog
/*
* The prolog script to create the stand-alone saved state for quadratization, it reads
* and returns an .ode file.
* Coder: M. Hemery
*/
% Loading libraries
:- use_module(arithmetic_rules).
:- use_module(biocham).
:- use_module(counters).
:- use_module(doc).
:- use_module(kinetics).
:- use_module(models).
:- use_module(namespace).
:- use_module(ode).
:- use_module(quadratic_reduction).
:- use_module(sat).
:- use_module(toplevel).
:- use_module(types).
:- use_module(util).
:- use_module(xpp_parser).
set_my_option(Option, Value) :-
change_item([], option, Option, option(Option: Value)).
% parse_arguments(+Arguments, -Input, -Method, -Timeout)
%
% Handle the various option of the standalone:
% --method: method used to perform the reduction
% --timeout: timeout for sat call
% --help: print help and quit
% -i: input file
parse_arguments([], none, sat_species, 120).
parse_arguments(['--method', Method| Argv], Input, Method, Timeout) :-
parse_arguments(Argv, Input, _M, Timeout).
parse_arguments(['--timeout', Timeout| Argv], Input, Method, Timeout) :-
parse_arguments(Argv, Input, Method, _T).
parse_arguments(['-i', Input| Argv], Input, Method, Timeout) :-
parse_arguments(Argv, _I, Method, Timeout).
parse_arguments(['--help'| _Argv], _I, _M, _T) :-
print_help.
% print_help
print_help :-
format("help~n", []),
halt.
% main(+argv)
main(Argv) :-
(Argv = [] -> print_help ; true),
% Handling naming and options
parse_arguments(Argv, Input, Method, Timeout),
file_name_extension(InputBase, 'ode', Input),
string_concat(InputBase, '_quad', OutputBase),
file_name_extension(OutputBase, 'ode', Output),
format("Input: ~w~nOutput: ~w~n", [Input, Output]),
% Initialization (see biocham:initialize)
set_prolog_flag(allow_variable_name_as_functor, true),
set_counter(list_item_counter, 0),
set_counter(item_id, 0),
set_counter(model, 0),
set_counter(molecule_id, 0),
set_counter(parameter_id, 1),
nb_setval(ode_viewer, inline),
nb_setval(current_models, []),
new_model,
set_my_option(quadratic_reduction, Method),
set_my_option(timeout, Timeout),
set_my_option(stats, no),
% Main procedure
read_xpp(Input),
list_ode,
quadratic_reduction_ODE,
select_ode_system(quadratic_ode),
write_xpp(Output).
......@@ -7,6 +7,9 @@
]
).
:- use_module(objects).
:- use_module(reaction_editor).
kinetics(Reactants, Inhibitors, Kinetics, Value) :-
eval_kinetics(Reactants, Inhibitors, Kinetics, ValueToSimplify),
simplify(ValueToSimplify, Value).
......
......@@ -35,7 +35,7 @@ correct_model :-
(
remove_null_kinetic(Reaction),!
;
split_bidirectional_reaction(List_Molecule, Reaction),!
split_reaction(List_Molecule, Reaction),!
;
reverse_reaction(Reaction),!
;
......@@ -61,32 +61,51 @@ remove_null_kinetic(Expression for Reaction) :-
delete_reaction(Expression for Reaction)
).
%! split_bidirectional_reaction(+List_Mol, +Reaction)
%! split_reaction(+List_Mol, +Reaction)
%
% Split a bidirection reaction in two
split_bidirectional_reaction(List_Mol, Expression for Reactant => Product) :-
split_reaction(List_Mol, Expression for Reactant => Product) :-
(
ode:substitute_functions(Expression, Expr),
distribute(Expr, DistExpr),
DistExpr = Up - Down,
% test if it is effectively a bidir. reaction
member(P, List_Mol),
models:formal_product(P, Reactant => Product),
present_in_kinetics(P, Down)
->
delete_reaction(Expression for Reactant => Product),
(
Reactant = Reac/Inhib
DistExpr = _A + _B
->
add_reaction(Up for Reactant => Product),
add_reaction(Down for Product/Inhib => Reac)
delete_reaction(Expression for Reactant => Product),
split_reaction_sr(List_Mol, DistExpr for Reactant => Product)
;
add_reaction(Up for Reactant => Product),
add_reaction(Down for Product => Reactant)
DistExpr = _Up - Down,
% test if it is effectively a bidir. reaction
member(P, List_Mol),
models:formal_product(P, Reactant => Product),
present_in_kinetics(P, Down)
->
delete_reaction(Expression for Reactant => Product),
split_reaction_sr(List_Mol, DistExpr for Reactant => Product)
)
).
split_reaction_sr(List_Mol, Expression + Monomial for Reactant => Product) :-
!,
add_reaction(Monomial for Reactant => Product),
split_reaction_sr(List_Mol, Expression for Reactant => Product).
split_reaction_sr(List_Mol, Expression - Monomial for Reactant => Product) :-
!,
(
Reactant = Reac/Inhib
->
add_reaction(Monomial for Product/Inhib => Reac)
;
add_reaction(Monomial for Product => Reactant)
),
split_reaction_sr(List_Mol, Expression for Reactant => Product).
split_reaction_sr(_List_Mol, Expression for Reactant => Product) :-
add_reaction(Expression for Reactant => Product).
%! reverse_reaction(+Reaction)
%
% Reverse a reaction if its kinetic is negative
......
......@@ -327,12 +327,17 @@ check_model :-
(
forall(
member(Reaction, List_Reaction),
test_wellformed(List_Molecule, Reaction)
(
test_wellformed(List_Molecule, Reaction)
;
format("Due to: ~w ; ", [Reaction]),
fail
)
)
->
format("The model is wellformed.~n", [])
;
format("The model is NOT wellformed.~n", [])
format("the model is NOT wellformed.~n", [])
),
(
test_strictness(List_Molecule, List_Reaction)
......@@ -1357,11 +1362,11 @@ formal_inhibitor(Mol, _Reactants/Inhibitors => _Products) :- !,
% is_present(Mol, Expr)
is_present(A, A) :- !.
is_present(A, _K*A) :- !.
is_present(A, _Tail+A) :- !.
is_present(A, _Tail+_K*A) :- !.
is_present(A, Tail+_H) :- is_present(A, Tail).
% is_present(A, A) :- !.
% is_present(A, _K*A) :- !.
% is_present(A, _Tail+A) :- !.
% is_present(A, _Tail+_K*A) :- !.
% is_present(A, Tail+_H) :- is_present(A, Tail).
%! test_wellformed(+List_Molecule, +Reaction)
......
......@@ -23,6 +23,7 @@
add_influences_from_ode_system/0,
(init)/1,
remove_fraction/0,
infer_hidden_molecules/0,
% load_ode_system/1,
load_reactions_from_ode/1,
% add_ode_system/1,
......@@ -49,6 +50,7 @@
import_influences_from_ode_system/1,
with_current_ode_system/1,
edit_current_ode_system/1,
infer_hidden_molecules/1,
all_odes/1,
ode/2,
ode/3,
......@@ -222,6 +224,8 @@ add_reactions_from_ode(InputFile) :-
'),
option(import_reactions_with_inhibitors, yesno, _,
'Add inhibitors when inferring reactions.'),
option(infer_hidden_molecules, yesno, _,
'Run infer_hidden_molecules (1-X => X_m) before loading reactions'),
retractall(flag_influences),
models:add_all('ode', InputFile).
......@@ -381,6 +385,8 @@ init(InitList) :-
)
).
:- initial(option(infer_hidden_molecules: no)).
add_reactions_from_ode_system :-
biocham_command,
doc('adds a set of reactions equivalent to the current ODE system.'),
......@@ -389,6 +395,14 @@ add_reactions_from_ode_system :-
get_current_ode_system(Id),
!,
retractall(flag_influences),
(
get_option(infer_hidden_molecules, IHM),
IHM = yes
->
infer_hidden_molecules(Id)
;
true
),
import_reactions_from_ode_system(Id).
add_reactions_from_ode_system :-
print_message(warning,'there is no current ode system').
......@@ -413,6 +427,8 @@ load_reactions_from_ode_system :-
\\command{add_reactions_from_ode_system/0}.'),
option(import_reactions_with_inhibitors, yesno, _,
'Add inhibitors when inferring reactions.'),
option(infer_hidden_molecules, yesno, _,
'Run infer_hidden_molecules (1-X => X_m) before loading reactions'),
delete_reactions,
add_reactions_from_ode_system.
......@@ -453,7 +469,15 @@ remove_fraction :-
)),
select_ode_system(without_fraction_ode).
% For now a biocham command mainly to debug, remove also ode_system when switching to API
infer_hidden_molecules :-
biocham_command,
doc('Blabla'),
ode_system,
with_current_ode_system((
get_current_ode_system(Id),
infer_hidden_molecules(Id)
)).
:- devdoc('\\section{Public API}').
......@@ -667,18 +691,16 @@ parse_ode(ODE, X, Item) :-
enumerate_terms_in_odes(Id) :-
\+ (
forall(
ode(Id, X, E),