Commit af524052 authored by MARTINEZ Thierry 's avatar MARTINEZ Thierry

First class ODE systems

parent 5e4ba131
ADDITIONAL_MODULES=modules/sbml/sbml_utils.pl
MODULES=$(shell sed -n -E 's/^- (.*\.pl)$$/\1/p' toc.org) $(ADDITIONAL_MODULES)
MODULES=$(shell sed -n -E 's/^[+-] (.*\.pl)$$/\1/p' toc.org) \
$(ADDITIONAL_MODULES)
# load_test_files/1 should make this useless, but I cannot find how to use it
TEST_MODULES=$(wildcard $(MODULES:.pl=.plt))
......
......@@ -62,7 +62,7 @@ check_conservations :-
'checks all conservation laws against reactions, and if necessary kinetics
(see also \\command{add_conservation/1}).'
),
compute_ode,
ode_system,
forall(
(
% model: current_model? compute_ode already supposes a single model
......
......@@ -245,6 +245,10 @@ read_toc(Stream, Type) :-
;
atom_codes(Atom, Line),
(
atom_concat('%', _, Atom)
->
true
;
atom_concat('* ', Section, Atom)
->
add_section(1, Section)
......@@ -257,7 +261,7 @@ read_toc(Stream, Type) :-
->
add_section(3, Section)
;
Atom = '+ index'
Atom = '- index'
->
assertz(toc(index))
;
......
......@@ -10,6 +10,7 @@
delete_vertex/1,
edge/1,
add_edge/1,
edgeref/1,
delete_edge/1,
list_edges/0,
list_isolated_vertices/0,
......@@ -116,15 +117,22 @@ add_edge(EdgeList) :-
).
delete_edge(EdgeList) :-
:- grammar(edgeref).
edgeref(Edge) :-
edge(Edge).
delete_edge(EdgeRefList) :-
biocham_command(*),
type(EdgeList, '*'(edge)),
type(EdgeRefList, '*'(edgeref)),
doc('Deletes a set of edges from the current graph.'),
get_current_graph(GraphId),
\+ (
member(Edge, EdgeList),
member(EdgeRef, EdgeRefList),
\+ (
delete_item([parent: GraphId, kind: edge, key: Edge])
delete_item([parent: GraphId, kind: edge, key: EdgeRef])
)
).
......@@ -258,11 +266,11 @@ list_attributes(GraphObject) :-
get_current_graph(Id) :-
get_selection(current_graph, [Id]).
get_selection(current_model, current_graph, [Id]).
set_current_graph(Id) :-
set_selection(current_graph, [Id]).
set_selection(current_model, current_graph, [Id]).
get_graph_name(Id, Name) :-
......
......@@ -37,8 +37,8 @@
inherits/1,
inherits/2,
begin_command/0,
get_selection/2,
set_selection/2,
get_selection/3,
set_selection/3,
at_delete/2,
get_parent/2
]
......@@ -338,7 +338,7 @@ load_biocham_stream(Stream) :-
current_models(Models) :-
get_selection(current_models, Models).
get_selection(top, current_models, Models).
set_current_models(Models) :-
......@@ -347,7 +347,7 @@ set_current_models(Models) :-
->
throw(error(cannot_select_no_models))
;
set_selection(current_models, Models)
set_selection(top, current_models, Models)
).
......@@ -360,7 +360,7 @@ new_model(Id) :-
;
add_item([parent: top, kind: model, key: new_model, id: Id]),
(
item([parent:top, kind: model, key: initial, id: InitialId])
item([parent: top, kind: model, key: initial, id: InitialId])
->
inherits(Id, InitialId)
;
......@@ -567,7 +567,7 @@ all_ids(Options, Ids) :-
optional(Item, List) :-
(
memberchk(Item, List)
member(Item, List)
->
true
;
......@@ -646,7 +646,7 @@ list_ids_aux(Options, Ids) :-
item(Id, _, _, Item),
indent(Indent),
(
selection(_, Id)
selection(_, _, Id)
->
Selected = '*'
;
......@@ -710,12 +710,12 @@ delete_item(Id) :-
)
)
),
retractall(selection(_, Id)),
retractall(selection(_, _, Id)),
(
KindItem = model
->
(
selection(current_models, _)
selection(top, current_models, _)
->
true
;
......@@ -829,18 +829,21 @@ inherits_from(Id, AncestorId) :-
inherits_from(IntermediateAncestorId, AncestorId).
:- dynamic(selection/2).
:- dynamic(selection/3).
get_selection(SelectionName, Ids) :-
get_selection(ParentOrCurrentModel, SelectionName, Ids) :-
get_selection_parent(ParentOrCurrentModel, Parent),
findall(
Id,
selection(SelectionName, Id),
Ids).
selection(Parent, SelectionName, Id),
Ids
).
set_selection(SelectionName, IdsOrOptions) :-
retractall(selection(SelectionName, _)),
set_selection(ParentOrCurrentModel, SelectionName, IdsOrOptions) :-
get_selection_parent(ParentOrCurrentModel, Parent),
retractall(selection(Parent, SelectionName, _)),
\+ (
(
IdsOrOptions = [_: _| _]
......@@ -850,13 +853,24 @@ set_selection(SelectionName, IdsOrOptions) :-
member(Id, IdsOrOptions)
),
\+ (
assertz(selection(SelectionName, Id))
assertz(selection(Parent, SelectionName, Id))
)
).
get_selection_parent(ParentOrCurrentModel, Parent) :-
(
ParentOrCurrentModel = current_model
->
single_model(Parent)
;
Parent = ParentOrCurrentModel
).
:- dynamic(at_delete_goal/2).
at_delete(Id, Goal) :-
assertz(at_delete_goal(Id, Goal)).
......
......@@ -28,7 +28,7 @@ numerical_simulation(Time) :-
biocham_command,
type(Time, time),
doc('performs a numerical simulation up to a given time.'),
compute_ode,
ode_system,
make_ode_system,
solve(Time).
......
:- module(
ode,
[
list_ODE/0,
compute_ode/0,
new_ode_system/0,
delete_ode_system/1,
set_ode_system_name/1,
list_ode_systems/0,
select_ode_system/1,
ode/1,
add_ode/1,
oderef/1,
delete_ode/1,
list_ode/0,
ode_system/0,
import_ode/1,
import_reactions_from_ode_system/0,
get_current_ode_system/1,
all_odes/1,
ode/2,
load_ode/1
op(1000, fx, init)
]).
:- dynamic(ode/2).
list_ODE :-
:- devdoc('\\section{Commands}').
new_ode_system :-
biocham_command,
doc('creates an ODE system.'),
add_item([kind: ode_system, key: new_ode_system, id: Id]),
set_current_ode_system(Id).
delete_ode_system(Name) :-
biocham_command,
type(Name, name),
doc('deletes an ODE system.'),
delete_items([kind: ode_system, key: Name]).
set_ode_system_name(Name) :-
biocham_command,
type(Name, name),
doc('sets the name of the current ODE system.'),
get_current_ode_system(Id),
replace_item(Id, ode_system, Name, Name).
list_ode_systems :-
biocham_command,
doc('lists the ODE systems of the current model.'),
list_items([kind: ode_system]).
select_ode_system(Name) :-
biocham_command,
type(Name, name),
doc('selects an ODE system'),
find_item([kind: ode_system, key: Name, id: Id]),
set_current_ode_system(Id).
:- grammar(ode).
ode(d(X)/dt = E) :-
name(X),
arithmetic_expression(E).
ode(DX/dt = E) :-
name(DX),
arithmetic_expression(E).
add_ode(ODE) :-
ode(ODE),
!,
(
ODE = (d(X)/dt = _)
->
Item = ODE
;
ODE = (DX/dt = E),
atom_concat(d, X, DX)
->
Item = (d(X)/dt = E)
;
throw(error(illformed_ode(ODE)))
),
get_current_ode_system(Id),
add_item([parent: Id, kind: ode, key: X, item: Item]).
add_ode(ODEList) :-
biocham_command(*),
type(ODEList, '*'(ode)),
doc('adds the given set of ODEs to the current ODE system.'),
\+ (
member(ODE, ODEList),
\+ (
add_ode(ODE)
)
).
import_reactions_from_ode_system :-
biocham_command,
doc('adds the reactions that match the current ODE system.'),
\+ (
ode(X, E),
monomial(Monomial, E),
\+ (
(
Monomial = - OppositeMonomial,
coefficient_species(OppositeMonomial, Coefficient, Species)
->
add_reaction(Coefficient * Species for X => '_')
;
coefficient_species(Monomial, Coefficient, Species)
->
add_reaction(Coefficient * Species for '_' => X)
)
)
).
:- grammar(oderef).
oderef(X) :-
name(X).
delete_ode(ODERefList) :-
biocham_command(*),
type(ODERefList, '*'(oderef)),
doc('removes the given set of ODEs from the current ODE system.'),
\+ (
member(Name, ODERefList),
\+ (
delete_single_ode(Name)
)
).
list_ode :-
biocham_command,
doc('
returns the set of ordinary differential equations
and initial concentrations (one line per molecule).
lists the set of ODES of the current ODE system.
\\begin{example}
'),
biocham_silent(clear_model),
biocham(a => b),
biocham(present(a)),
biocham(list_ODE),
biocham(ode_system),
biocham(list_ode),
doc('
\\end{example}
'),
get_current_ode_system(Id),
list_items([parent: Id, kind: ode]).
ode_system :-
biocham_command,
doc('builds the set of ODES corresponding to the reactions of the current model.'),
delete_ode_system('ode_system'),
new_ode_system,
set_ode_system_name('ode_system'),
clean_draft_ode,
compute_ode,
print_ode.
simplify_ode,
put_ode_into_system,
clean_draft_ode.
import_ode(InputFile) :-
biocham_command,
type(InputFile, input_file),
doc('imports a set of ODEs.'),
new_ode_system,
set_ode_system_name(InputFile),
setup_call_cleanup(
open(InputFile, read, Stream),
read_ode(Stream),
close(Stream)
).
:- devdoc('\\section{Public API}').
get_current_ode_system(Id) :-
get_selection(current_model, current_ode_system, [Id]).
all_odes(ODEs) :-
get_current_ode_system(Id),
all_items([parent: Id, kind: ode], ODEs).
ode(X, E) :-
get_current_ode_system(Id),
(
var(X)
->
item([parent: Id, kind: ode, item: (d(X)/dt = E)])
;
find_item([parent: Id, kind: ode, key: X, item: (d(X)/dt = E)])
).
:- devdoc('\\section{Private predicates}').
set_current_ode_system(Id) :-
set_selection(current_model, current_ode_system, [Id]).
:- dynamic(draft_ode/2).
clean_draft_ode :-
retractall(draft_ode(_, _)).
compute_ode :-
retractall(ode(_, _)),
\+ (
item([kind: reaction, item: Item]),
reaction(Item, Kinetics, Reactants, Products),
......@@ -37,16 +232,15 @@ compute_ode :-
add_molecules(Reactants, -KineticsExpression),
add_molecules(Products, KineticsExpression)
)
),
simplify_ode.
).
simplify_ode :-
\+ (
retract(ode(Variable, Expression)),
retract(draft_ode(Variable, Expression)),
\+ (
simplify(Expression, SimplifiedExpression),
assertz(ode(Variable, SimplifiedExpression))
assertz(draft_ode(Variable, SimplifiedExpression))
)
).
......@@ -62,58 +256,153 @@ add_molecules(Molecules, Kinetics) :-
add_molecule(Coefficient * Molecule, Kinetics) :-
(
retract(ode(Molecule, Expression))
retract(draft_ode(Molecule, Expression))
->
true
;
Expression = 0
),
asserta(ode(Molecule, Expression + Coefficient * Kinetics )).
asserta(draft_ode(Molecule, Expression + Coefficient * Kinetics)).
print_ode :-
write('Init \tODE\n'),
put_ode_into_system :-
\+ (
ode(Molecule, Equation),
get_initial_concentration(Molecule, InitialConcentration),
draft_ode(X, Expression),
\+ (
format('~f \td[~w]/dt=~w\n', [InitialConcentration, Molecule, Equation])
add_ode(d(X)/dt = Expression)
)
).
load_ode(InputFile) :-
biocham_command,
type(InputFile, input_file),
doc('
creates a new Biocham model (reaction rules, initial state, parameters and
macros) inferred from an xppaut .ode file (preferably written with
parameters).
'),
external_convert_ode(InputFile, 'external_ode.bc'),
load_biocham('external_ode.bc').
delete_single_ode(Name) :-
get_current_ode_system(Id),
delete_item([parent: Id, kind: ode, key: Name]).
add_ode(InputFile) :-
biocham_command,
type(InputFile, input_file),
doc('
adds reaction rules, initial state, parameters and
macros inferred from an xppaut .ode file (preferably written with
parameters).
'),
external_convert_ode(InputFile, 'external_ode.bc'),
add_biocham('external_ode.bc').
read_ode(Stream) :-
get_current_ode_system(Id),
repeat,
read_line_to_codes(Stream, Codes),
(
Codes = end_of_file
->
!
;
read_term_from_codes(Codes, Expression, []),
(
Expression = end_of_file
->
fail
;
Expression = done
->
!
;
(
Expression = (init X = _E)
->
add_item([parent: Id, kind: initial_concentration, key: X, item: Expression])
;
Expression = (_DX/dt = _E)
->
add_ode(Expression)
)
->
fail
;
throw(error(syntax_error_in_ode_system(Codes, Expression)))
)
).
external_convert_ode(InputFile, OutputFile) :-
setup_call_cleanup(
open('external_load_ode.bc', write, Ode_bc),
monomial(Monomial, A + B) :-
!,
(
portray_clause(Ode_bc, load_ode(InputFile)),
portray_clause(Ode_bc, export_biocham(OutputFile)),
portray_clause(Ode_bc, 'quit')
),
close(Ode_bc)
),
call_subprocess(path(biocham), ['external_load_ode.bc'], [stdout(pipe(_))]).
monomial(Monomial, A)
;
monomial(Monomial, B)
).
monomial(Monomial, A - B) :-
!,
(
monomial(Monomial, A)
;
monomial_opposite(Monomial, B)
).
monomial(A, A).
monomial_opposite(Monomial, A + B) :-
!,
(
monomial_opposite(Monomial, A)
;
monomial_opposite(Monomial, B)
).
monomial_opposite(Monomial, A - B) :-
!,
(
monomial_opposite(Monomial, A)
;
monomial(Monomial, B)
).
monomial_opposite(A, - A).
coefficient_species(Coefficient * Species, Coefficient, Species) :-
!.
coefficient_species(Species, 1, Species).
%
%
%print_ode :-
% write('Init \tODE\n'),
% \+ (
% ode(Molecule, Equation),
% get_initial_concentration(Molecule, InitialConcentration),
% \+ (
% format('~f \td[~w]/dt=~w\n', [InitialConcentration, Molecule, Equation])
% )
% ).
%
%
%load_ode(InputFile) :-
% biocham_command,
% type(InputFile, input_file),
% doc('
% creates a new Biocham model (reaction rules, initial state, parameters and
% macros) inferred from an xppaut .ode file (preferably written with
% parameters).
% '),
% external_convert_ode(InputFile, 'external_ode.bc'),
% load_biocham('external_ode.bc').
%
%
%add_ode(InputFile) :-
% biocham_command,
% type(InputFile, input_file),
% doc('
% adds reaction rules, initial state, parameters and
% macros inferred from an xppaut .ode file (preferably written with
% parameters).
% '),
% external_convert_ode(InputFile, 'external_ode.bc'),
% add_biocham('external_ode.bc').
%
%
%external_convert_ode(InputFile, OutputFile) :-
% setup_call_cleanup(
% open('external_load_ode.bc', write, Ode_bc),
% (
% portray_clause(Ode_bc, load_ode(InputFile)),
% portray_clause(Ode_bc, export_biocham(OutputFile)),
% portray_clause(Ode_bc, 'quit')
% ),
% close(Ode_bc)
% ),
% call_subprocess(path(biocham), ['external_load_ode.bc'], [stdout(pipe(_))]).
......@@ -3,18 +3,16 @@
:- begin_tests(ode).
test('compute_ode', [true(ODEs == [ode(b, a), ode(a, -a)])]) :-
test('ode_system', [true(ODEs == [d(b)/dt = a, d(a)/dt = -a])]) :-
clear_model,
command(a => b),
compute_ode,
findall(
ode(Variable, Expression),
ode(Variable, Expression),
ODEs).
ode_system,
all_odes(ODEs).
test('load_ode', [true(Reactions = [[a] for a => b])]) :-
clear_model,
test(
'import_ode',
[true(ODEs == [init a = 0, init b = 0, d(a)/dt = -a, d(b)/dt = a])]
) :-
setup_call_cleanup(
open('test.ode', write, TestOde),
write(TestOde, '
......@@ -26,7 +24,19 @@ done
'),
close(TestOde)
),
load_ode('test.ode'),
import_ode('test.ode'),
get_current_ode_system(Id),
all_items([parent: Id], ODEs).
test(
'import_reactions_from_ode_system',
[true(Reactions == [1*a for a=>'_', 1*a for '_'=>b])]
) :-
clear_model,
new_ode_system,
add_ode(da/dt = -a),
add_ode(db/dt = a),
import_reactions_from_ode_system,
all_items([kind: reaction], Reactions).
......
......@@ -102,6 +102,9 @@ draw_reactions :-
draw_graph.
:- devdoc('\\section{Private predicates}').
:- dynamic(vertex_transition/1).
......
......@@ -6,7 +6,7 @@
solution/1,
patch_solution/2,
kinetics/1,
simple_kinetics/1,
arithmetic_expression/1,
op(800, xfx, for),
op(700, xfx, <=),
op(750, xfx, =>),
......@@ -106,71 +106,71 @@ patch_solution(Solution0, Solution1, Solution0 - Solution1).
:- grammar(kinetics).
kinetics(SimpleKinetics) :-
simple_kinetics(SimpleKinetics).
kinetics(ArithmeticExpression) :-
arithmetic_expression(ArithmeticExpression).
kinetics((SimpleKinetics0, SimpleKinetics1)) :-
simple_kinetics(SimpleKinetics0),
simple_kinetics(SimpleKinetics1).
kinetics((ArithmeticExpression0, ArithmeticExpression1)) :-
arithmetic_expression(ArithmeticExpression0),
arithmetic_expression(ArithmeticExpression1).
:- grammar(simple_kinetics).
:- grammar(arithmetic_expression).
simple_kinetics([A]) :-
simple_kinetics(A).
arithmetic_expression([A]) :-
arithmetic_expression(A).
simple_kinetics(SimpleKinetics0 + SimpleKinetics1) :-
simple_kinetics(SimpleKinetics0),
simple_kinetics(SimpleKinetics1).
arithmetic_expression(ArithmeticExpression0 + ArithmeticExpression1) :-