Commit 4530aac1 authored by MARTINEZ Thierry 's avatar MARTINEZ Thierry

ode

parent a85dc7fa
......@@ -22,6 +22,7 @@ alias(ObjectList) :-
EquivalenceClassList),
add_equivalence_class_list(EquivalenceClassList).
add_equivalence_class_list(EquivalenceClassList) :-
(
EquivalenceClassList = [_]
......@@ -32,7 +33,8 @@ add_equivalence_class_list(EquivalenceClassList) :-
add_item([
kind: alias, key: EquivalenceClassList, item: alias(EquivalenceClass)
])
).
),
simplify_all_reactions.
instance_member_delete(Object, ObjectList) :-
......
......@@ -6,7 +6,6 @@ test('alias', [true(Reactions == [2 * a => c])]) :-
clear_model,
command(a + b => c),
command(alias(a = b)),
command(simplify_all_reactions),
all_items([kind: reaction], Reactions).
:- end_tests(aliases).
:- module(
arithmetic_simplify,
arithmetic_rules,
[
simplify/2
simplify/2,
distribute/2,
additive_normal_form/2
]).
simplify(Value, Value) :-
......@@ -300,3 +302,66 @@ simplify_cos(A, cos(A)).
simplify_sin(A, sin(A)).
distribute(In, Out) :-
rewrite(arithmetic_rules:rewrite_distribute, In, Out).
rewrite_distribute((A + B) * C, A * C + B * C).
rewrite_distribute(A * (B + C), A * B + A * C).
rewrite_distribute((A - B) * C, A * C - B * C).
rewrite_distribute(A * (B - C), A * B - A * C).
rewrite_distribute((A + B) / C, A / C + B / C).
rewrite_distribute((A - B) / C, A / C - B / C).
rewrite_distribute(A - (B + C), A - B - C).
rewrite_distribute(- (A + B), - A - B).
additive_normal_form(In, Out) :-
rewrite(arithmetic_rules:rewrite_additive_normal_form, In, Out).
rewrite_additive_normal_form(A, B) :-
rewrite_eval(A, B).
rewrite_additive_normal_form(A, B) :-
rewrite_distribute(A, B).
rewrite_additive_normal_form(A * B, B * A) :-
number(B).
rewrite_additive_normal_form((A * B) * C, A * (B * C)) :-
number(A).
rewrite_additive_normal_form(- A, (-1) * A).
rewrite_additive_normal_form(A - B, A + (-1) * B).
rewrite_eval(In, Out) :-
arithmetic_operation(In),
In =.. [_ | Args],
maplist(number, Args),
!,
Out is In.
arithmetic_operation(_ + _).
arithmetic_operation(_ - _).
arithmetic_operation(_ * _).
arithmetic_operation(_ / _).
arithmetic_operation(- _).
arithmetic_operation(_ ^ _).
:- use_module(library(plunit)).
:- begin_tests(arithmetic_rules).
test('distribute', [true(Out == a * a + a * b + c)]) :-
distribute(a * (a + b) + c, Out).
:- end_tests(arithmetic_rules).
......@@ -7,81 +7,74 @@
:- use_module(library(error)).
derivate(Expression, VariableIndex, Result) :-
derivate_raw(Expression, VariableIndex, UnsimplifiedResult),
derivate(Expression, Variable, Result) :-
derivate_raw(Expression, Variable, UnsimplifiedResult),
simplify(UnsimplifiedResult, Result).
derivate_raw([VariableIndex0], VariableIndex1, Result) :-
!,
(
VariableIndex0 = VariableIndex1
->
Result = 1
;
Result = 0
).
derivate_raw(Parameter, _VariableIndex, 0) :-
atom(Parameter),
derivate_raw(p(_ParameterIndex), _Variable, 0) :-
!.
derivate_raw(p(_ParameterIndex), _VariableIndex, 0) :-
!.
derivate_raw(Value, _VariableIndex, 0) :-
derivate_raw(Value, _Variable, 0) :-
number(Value),
!.
derivate_raw(- A, VariableIndex, - Aprime) :-
derivate_raw(- A, Variable, - Aprime) :-
!,
derivate_raw(A, VariableIndex, Aprime).
derivate_raw(A, Variable, Aprime).
derivate_raw(A + B, VariableIndex, Aprime + Bprime) :-
derivate_raw(A + B, Variable, Aprime + Bprime) :-
!,
derivate_raw(A, VariableIndex, Aprime),
derivate_raw(B, VariableIndex, Bprime).
derivate_raw(A, Variable, Aprime),
derivate_raw(B, Variable, Bprime).
derivate_raw(A - B, VariableIndex, Aprime - Bprime) :-
derivate_raw(A - B, Variable, Aprime - Bprime) :-
!,
derivate_raw(A, VariableIndex, Aprime),
derivate_raw(B, VariableIndex, Bprime).
derivate_raw(A, Variable, Aprime),
derivate_raw(B, Variable, Bprime).
derivate_raw(A * B, VariableIndex, A * Bprime + Aprime * B) :-
derivate_raw(A * B, Variable, A * Bprime + Aprime * B) :-
!,
derivate_raw(A, VariableIndex, Aprime),
derivate_raw(B, VariableIndex, Bprime).
derivate_raw(A, Variable, Aprime),
derivate_raw(B, Variable, Bprime).
derivate_raw(A / B, VariableIndex, (A * Bprime - Aprime * B) / (A ^ 2)) :-
derivate_raw(A / B, Variable, (A * Bprime - Aprime * B) / (A ^ 2)) :-
!,
derivate_raw(A, VariableIndex, Aprime),
derivate_raw(B, VariableIndex, Bprime).
derivate_raw(A, Variable, Aprime),
derivate_raw(B, Variable, Bprime).
derivate_raw(A ^ B, VariableIndex, Result) :-
derivate_raw(A ^ B, Variable, Result) :-
!,
derivate_raw(exp(B * log(A)), VariableIndex, Result).
derivate_raw(exp(B * log(A)), Variable, Result).
derivate_raw(sqrt(A), VariableIndex, Result) :-
derivate_raw(sqrt(A), Variable, Result) :-
!,
derivate_raw(A ^ 0.5, VariableIndex, Result).
derivate_raw(A ^ 0.5, Variable, Result).
derivate_raw(log(A), VariableIndex, Aprime / A) :-
derivate_raw(log(A), Variable, Aprime / A) :-
!,
derivate_raw(A, VariableIndex, Aprime).
derivate_raw(A, Variable, Aprime).
derivate_raw(exp(A), VariableIndex, Aprime * exp(A)) :-
derivate_raw(exp(A), Variable, Aprime * exp(A)) :-
!,
derivate_raw(A, VariableIndex, Aprime).
derivate_raw(A, Variable, Aprime).
derivate_raw(sin(A), VariableIndex, Aprime * cos(A)) :-
derivate_raw(sin(A), Variable, Aprime * cos(A)) :-
!,
derivate_raw(A, VariableIndex, Aprime).
derivate_raw(A, Variable, Aprime).
derivate_raw(cos(A), VariableIndex, - Aprime * sin(A)) :-
derivate_raw(cos(A), Variable, - Aprime * sin(A)) :-
!,
derivate_raw(A, VariableIndex, Aprime).
derivate_raw(A, Variable, Aprime).
derivate_raw(E, _VariableIndex, _Aprime) :-
type_error(arithmetic_expression, E).
derivate_raw(Variable0, Variable1, Result) :-
!,
(
Variable0 = Variable1
->
Result = 1
;
Result = 0
).
jacobian(Equations, Jacobian) :-
......@@ -101,6 +94,6 @@ jacobian_line(Equation, VariableCount, Line) :-
Derivative,
(
between(0, VariableMax, VariableIndex),
derivate(Equation, VariableIndex, Derivative)
derivate(Equation, [VariableIndex], Derivative)
),
Line).
......@@ -3,16 +3,16 @@
:- begin_tests(formal_derivation).
test('dx/dt(x) = 1', [true(E == 1)]) :-
derivate([x], x, E).
derivate(x, x, E).
test('dx/dt(y) = 0', [true(E == 0)]) :-
derivate([y], x, E).
derivate(y, x, E).
test('dx/dt(x^2) = 2x', [true(E == 2 * [x])]) :-
derivate([x] ^ 2, x, E).
test('dx/dt(x^2) = 2x', [true(E == 2 * x)]) :-
derivate(x ^ 2, x, E).
test('dx/dt(cos(sqrt(x))) = - 0.5 / sqrt(x) * sin(sqrt(x))',
[true(E == - (0.5 / sqrt([x])) * sin(sqrt([x])))]) :-
derivate(cos(sqrt([x])), x, E).
[true(E == - (0.5 / sqrt(x)) * sin(sqrt(x)))]) :-
derivate(cos(sqrt(x)), x, E).
:- end_tests(formal_derivation).
......@@ -5,7 +5,12 @@
]
).
kinetics(Reactants, single_reactant, Value) :-
kinetics(Reactants, Kinetics, Value) :-
eval_kinetics(Reactants, Kinetics, ValueToSimplify),
simplify(ValueToSimplify, Value).
eval_kinetics(Reactants, single_reactant, Value) :-
!,
(
Reactants = [Value]
......@@ -15,7 +20,7 @@ kinetics(Reactants, single_reactant, Value) :-
throw(error(not_single_reactant_reaction))
).
kinetics(Reactants, product(Pattern in List, Expression), Value) :-
eval_kinetics(Reactants, product(Pattern in List, Expression), Value) :-
!,
(
Pattern = S * O
......@@ -33,13 +38,13 @@ kinetics(Reactants, product(Pattern in List, Expression), Value) :-
),
make_product(Reactants, S, O, Expression, Value).
kinetics(Reactants, FunctionApplication, Expression) :-
eval_kinetics(Reactants, FunctionApplication, Expression) :-
function_apply(FunctionApplication, NewBody),
!,
kinetics(Reactants, NewBody, Expression).
eval_kinetics(Reactants, NewBody, Expression).
kinetics(Reactants, Callable, Expression) :-
term_morphism(kinetics(Reactants), Callable, Expression).
eval_kinetics(Reactants, Callable, Expression) :-
term_morphism(kinetics:eval_kinetics(Reactants), Callable, Expression).
make_product([], _S, _O, _Expression, 1).
......
......@@ -28,6 +28,7 @@
add_ode/2,
import_reactions_from_ode_system/1,
with_current_ode_system/1,
edit_current_ode_system/1,
all_odes/1,
ode/2,
ode/3,
......@@ -88,21 +89,25 @@ ode(DX/dt = E) :-
add_ode(ODE) :-
ode(ODE),
!,
get_current_ode_system(Id),
!,
add_ode(Id, ODE).
add_ode(ODEList) :-
biocham_command(*),
type(ODEList, '*'(ode)),
doc('adds the given set of ODEs to the current ODE system.'),
\+ (
member(ODE, ODEList),
doc('
If there is a current ODE system, adds the given set of ODEs to it.
If there is no current ODE system, imports the reactions inferred from the given set of ODEs.'),
edit_current_ode_system((
\+ (
add_ode(ODE)
member(ODE, ODEList),
\+ (
add_ode(ODE)
)
)
).
)).
import_reactions_from_ode_system :-
......@@ -158,11 +163,11 @@ ode_system :-
delete_ode_system('ode_system'),
new_ode_system,
set_ode_system_name('ode_system'),
check_draft_ode_cleaned,
check_cleaned(ode:assoc/2),
compute_ode,
simplify_ode,
put_ode_into_system,
clean_draft_ode.
clean(ode:assoc/2).
import_ode(InputFile) :-
......@@ -267,6 +272,15 @@ ode(Id, X, E) :-
).
ode_variables(Id, X) :-
ode(Id, X, _).
ode_variables(Id, X) :-
item([parent: Id, kind: initial_concentration, item: (init X = _)]),
\+ item([parent: Id, kind: ode, key: X]).
set_ode_initial_value(Variable, Value) :-
get_current_ode_system(Id),
set_ode_initial_value(Id, Variable, Value).
......@@ -293,41 +307,172 @@ add_ode(Id, ODE) :-
import_reactions_from_ode_system(Id) :-
check_cleaned(ode:assoc/2),
enumerate_terms_in_odes(Id),
add_terms_as_reactions(Id),
clean(ode:assoc/2).
:- devdoc('\\section{Private predicates}').
enumerate_terms_in_odes(Id) :-
\+ (
ode(Id, X, E),
monomial(Monomial, E),
\+ (
enumerate_terms_in_ode(X, E)
)
).
enumerate_terms_in_ode(X, E) :-
\+ (
non_decomposable_term_in_expression(Coefficient, Term, E),
\+ (
(
retract(assoc(Term, Occurrences))
->
true
;
Occurrences = []
),
(
Monomial = - OppositeMonomial,
coefficient_species(OppositeMonomial, Coefficient, Species)
select((X, OtherCoefficient), Occurrences, OtherOccurrences)
->
add_reaction(Coefficient * Species for X => '_')
NewCoefficient is Coefficient + OtherCoefficient
;
coefficient_species(Monomial, Coefficient, Species)
NewCoefficient = Coefficient,
OtherOccurrences = Occurrences
),
(
NewCoefficient = 0
->
add_reaction(Coefficient * Species for '_' => X)
)
NewOccurrences = OtherOccurrences
;
NewOccurrences = [(X, NewCoefficient) | OtherOccurrences]
),
assertz(assoc(Term, NewOccurrences))
)
).
:- devdoc('\\section{Private predicates}').
non_decomposable_term_in_expression(Coefficient, Term, E) :-
additive_normal_form(E, ENormal),
non_decomposable_term_in_additive_normal_form(Coefficient, Term, ENormal).
:- dynamic(draft_ode/2).
non_decomposable_term_in_additive_normal_form(Coefficient, Term, A + B) :-
!,
(
non_decomposable_term_in_additive_normal_form(Coefficient, Term, A)
;
non_decomposable_term_in_additive_normal_form(Coefficient, Term, B)
).
check_draft_ode_cleaned :-
non_decomposable_term_in_additive_normal_form(Coefficient, Term, Coefficient * Term) :-
number(Coefficient),
!.
non_decomposable_term_in_additive_normal_form(1, Term, Term).
:- dynamic(reactant/2).
:- dynamic(product/2).
:- dynamic(inhibitor/2).
add_terms_as_reactions(Id) :-
\+ (
assoc(Term, Occurrences),
\+ (
add_term_as_reactions(Id, Term, Occurrences)
)
).
add_term_as_reactions(Id, Term, Occurrences) :-
check_cleaned(ode:reactant/2),
check_cleaned(ode:product/2),
check_cleaned(ode:inhibitor/2),
\+ (
member((X, Coefficient), Occurrences),
\+ (
(
Coefficient < 0
->
OppCoefficient is - Coefficient,
assertz(reactant(X, OppCoefficient))
;
assertz(product(X, Coefficient))
)
)
),
\+ (
ode_variables(Id, X),
\+ (
(
\+ reactant(X, _),
partial_has_pos_val(Term, X)
->
assertz(reactant(X, 1)),
(
retract(product(X, OldCoefficient))
->
true
;
OldCoefficient = 0
),
NewCoefficient is OldCoefficient + 1,
assertz(product(X, NewCoefficient))
;
true
),
(
partial_has_pos_val(- Term, X)
->
assertz(inhibitor(X, 1))
;
true
)
)
),
findall(
Coefficient * Object,
reactant(Object, Coefficient),
Reactants
),
findall(
Coefficient * Object,
inhibitor(Object, Coefficient),
_Inhibitors
),
findall(
Coefficient * Object,
product(Object, Coefficient),
Products
),
add_reaction(Term, Reactants, Products, false),
clean(ode:reactant/2),
clean(ode:product/2),
clean(ode:inhibitor/2).
partial_has_pos_val(Term, X) :-
derivate(Term, X, DTerm),
(
draft_ode(_, _)
number(DTerm),
DTerm =< 0
->
throw(error(check_draft_ode_cleaned))
fail
;
true
).
clean_draft_ode :-
retractall(draft_ode(_, _)).
:- dynamic(assoc/2).
compute_ode :-
......@@ -344,10 +489,10 @@ compute_ode :-
simplify_ode :-
\+ (
retract(draft_ode(Variable, Expression)),
retract(assoc(Variable, Expression)),
\+ (
simplify(Expression, SimplifiedExpression),
assertz(draft_ode(Variable, SimplifiedExpression))
assertz(assoc(Variable, SimplifiedExpression))
)
).
......@@ -363,18 +508,18 @@ add_molecules(Molecules, Kinetics) :-
add_molecule(Coefficient * Molecule, Kinetics) :-
(
retract(draft_ode(Molecule, Expression))
retract(assoc(Molecule, Expression))
->
true
;
Expression = 0
),
asserta(draft_ode(Molecule, Expression + Coefficient * Kinetics)).
asserta(assoc(Molecule, Expression + Coefficient * Kinetics)).
put_ode_into_system :-
\+ (
draft_ode(X, Expression),
assoc(X, Expression),
\+ (
add_ode(d(X)/dt = Expression),
get_initial_concentration(X, Concentration),
......@@ -458,7 +603,7 @@ monomial_opposite(Monomial, A - B) :-
monomial(Monomial, B)
).
monomial_opposite(A, - A).
monomial_opposite(- A, A).
coefficient_species(Coefficient * Species, Coefficient, Species) :-
......@@ -509,3 +654,21 @@ with_current_ode_system(Goal) :-
)
)
).
edit_current_ode_system(Goal) :-
(
get_current_ode_system(_)
->
Goal
;
setup_call_cleanup(
new_ode_system,
Goal,
(
get_current_ode_system(Id),
import_reactions_from_ode_system(Id),
delete_item(Id)
)
)
).
......@@ -30,7 +30,7 @@ done
test(
'import_reactions_from_ode_system',
[true(Reactions == [1*a for a=>'_', 1*a for '_'=>b])]
[true(Reactions == [a => b])]
) :-
clear_model,
new_ode_system,
......
:- module(
odefunction,
[
'add_function'/1
add_function/1
]
).
......@@ -32,20 +32,20 @@ add_function(FunctionList) :-
export_ode(_, OdeSystem) :-
enumerate_molecules(L), % useless if item backtracks on Object (not clear to FF)
member(Object, L),
get_initial_concentration(Object, Concentration),
set_ode_initial_value(OdeSystem, Object, Concentration),
fail.
enumerate_molecules(L), % useless if item backtracks on Object (not clear to FF)
member(Object, L),
get_initial_concentration(Object, Concentration),
set_ode_initial_value(OdeSystem, Object, Concentration),
fail.
export_ode(FunctionList, OdeSystem) :-
member(Variable = _, FunctionList),
set_ode_initial_value(OdeSystem, Variable, 0),
fail.
member(Variable = _, FunctionList),
set_ode_initial_value(OdeSystem, Variable, 0),
fail.
export_ode(FunctionList, OdeSystem) :-
member(Variable = Value, FunctionList),
add_ode(OdeSystem, d(Variable)/dt = Value - Variable),
fail.
member(Variable = Value, FunctionList),
add_ode(OdeSystem, d(Variable)/dt = Value - Variable),
fail.
export_ode(_,_).
......@@ -2,16 +2,21 @@
reaction_editor,
[
add_reaction/1,
add_reaction/4,
list_reactions/0,
simplify_all_reactions/0,
add_reaction/4,
reaction/5,
reaction/4,
simplify_all_reactions/0,
enumerate_molecules/1,
solution_to_list/2
solution_to_list/2,
list_to_solution/2
]
).
:- devdoc('\\section{Commands}').
add_reaction(Reaction) :-
biocham_command,
type(Reaction, reaction),
......@@ -19,12 +24,28 @@ add_reaction(Reaction) :-
adds reaction rules to the current set of reactions.
This command is implicit: reaction rules can be added directly in models.
'),
add_item([kind: reaction, item: Reaction]).
\+ (
reaction(Reaction, Kinetics, Reactants, Products),
\+ (
add_reaction(Kinetics, Reactants, Products, false)
)
).
add_reaction(Kinetics, Left, Right, Reversible) :-
make_reaction(Kinetics, Left, Right, Reversible, Reaction),
add_reaction(Reaction).
%simplify_all_reactions :-
% biocham_command,
% doc('
% replaces each reaction by a simplified form, by grouping common molecules,
% identifying catalysts, and by using the canonical molecule for aliases.
% '),
% \+ (
% item([kind: reaction, id: Id, item: Reaction]),
% \+ (
% delete_item(Id),
% simplify_reaction(Reaction, SimplifiedReaction),
% add_item([kind: reaction, item: SimplifiedReaction])
% )
% ).
list_reactions :-
......@@ -33,12 +54,21 @@ list_reactions :-
list_items([kind: reaction]).
:- devdoc('\\section{Public API}').
add_reaction(Kinetics, Left, Right, Reversible) :-
make_reaction(Kinetics, Left, Right, Reversible, Reaction),
add_item([kind: reaction, item: Reaction]).
simplify_all_reactions :-
biocham_command,
doc('
replaces each reaction by a simplified form, by grouping common molecules,
identifying catalysts, and by using the canonical molecule for aliases.
'</