Commit 8d0282bf authored by Thierry Martinez's avatar Thierry Martinez

Kinetics

parent 7573596d
MODULES= \
biocham.pl \
about.pl \
commandline.pl \
toplevel.pl \
models.pl \
rules.pl \
objects.pl \
ode.pl \
initial_state.pl \
kinetics.pl \
numerical_simulation.pl \
util.pl \
simplify_arithmetic.pl \
formal_derivation.pl \
......@@ -16,10 +21,10 @@ MODULES= \
all: biocham biocham_debug tests
biocham: $(MODULES) Makefile
swipl -o biocham --goal=about --toplevel=toplevel -c $(MODULES)
swipl -o biocham --goal=start --toplevel=toplevel -c $(MODULES)
biocham_debug: $(MODULES) Makefile
swipl -o biocham_debug -c $(MODULES)
swipl -o biocham_debug --goal=initialize -c $(MODULES)
tests: biocham_tests
./biocham_tests
......
:- module(
biocham,
[
start/0,
initialize/0,
biocham_command/0,
biocham_command/1
biocham_command/1,
load/1,
load_biocham/1
]).
start :-
about,
initialize,
do_arguments.
initialize :-
prolog_history(enable),
set_counter(item_id, 0).
biocham_command.
biocham_command(Functor/Arity) :-
functor(Head, Functor, Arity),
once(clause(Head, (biocham_command, _))).
load(Filename) :-
biocham_command,
load_biocham(Filename).
load_biocham(Filename) :-
biocham_command,
setup_call_cleanup(
open(Filename, read, Stream),
load_biocham_stream(Stream),
close(Stream)
).
load_biocham_stream(Stream) :-
\+ (
repeat,
read_term(Stream, Command, [variable_names(VariableNames)]),
(
Command = end_of_file
->
!,
fail
;
name_vars_and_anonymous(Command, VariableNames),
execute_command(Command),
fail
)
).
:- module(
commandline,
[
do_arguments/0
]
).
do_arguments :-
current_prolog_flag(argv, Argv),
do_arguments(Argv).
:- dynamic(file/1).
do_arguments(Argv) :-
retractall(file(_)),
parse_arguments(Argv),
load_files.
load_files :-
\+ (
file(File),
\+ (
load(File)
)
).
parse_arguments([]).
parse_arguments([Arg | Argv]) :-
(
Arg = '--'
->
add_files(Argv)
;
atom_concat('--', _, Arg)
->
(
option(_, Arg, Argv, ArgTail)
->
parse_arguments(ArgTail)
;
throw(error(unknown_option(Arg)))
)
;
add_file(Arg),
parse_arguments(Argv)
).
option('t', '--trace', Arg, Arg) :-
leash(-all),
trace.
add_files(Files) :-
\+ (
member(File, Files),
\+ (
add_file(File)
)
).
add_file(File) :-
assertz(file(File)).
......@@ -2,14 +2,20 @@
counters,
[
set_counter/2,
count/2
count/2,
peek_count/2
]
).
set_counter(CounterName, InitialValue) :-
nb_setval(CounterName, InitialValue).
count(CounterName, Value) :-
nb_getval(CounterName, Value),
NextValue is Value + 1,
nb_setval(CounterName, NextValue).
peek_count(CounterName, Value) :-
nb_getval(CounterName, Value).
......@@ -11,20 +11,25 @@ derivate(Expression, VariableIndex, Result) :-
derivate_raw(Expression, VariableIndex, UnsimplifiedResult),
simplify(UnsimplifiedResult, Result).
derivate_raw(x(VariableIndex0), VariableIndex1, Result) :-
derivate_raw([VariableIndex0], VariableIndex1, Result) :-
!,
(
VariableIndex0 = VariableIndex1
->
Result = v(1)
Result = 1
;
Result = v(0)
Result = 0
).
derivate_raw(p(_Value), _VariableIndex, v(0)) :-
derivate_raw(Parameter, _VariableIndex, 0) :-
atom(Parameter),
!.
derivate_raw(v(_Value), _VariableIndex, v(0)) :-
derivate_raw(p(_ParameterIndex), _VariableIndex, 0) :-
!.
derivate_raw(Value, _VariableIndex, 0) :-
number(Value),
!.
derivate_raw(- A, VariableIndex, - Aprime) :-
......@@ -46,7 +51,7 @@ derivate_raw(A * B, VariableIndex, A * Bprime + Aprime * B) :-
derivate_raw(A, VariableIndex, Aprime),
derivate_raw(B, VariableIndex, Bprime).
derivate_raw(A / B, VariableIndex, (A * Bprime - Aprime * B) / (A ^ v(2))) :-
derivate_raw(A / B, VariableIndex, (A * Bprime - Aprime * B) / (A ^ 2)) :-
!,
derivate_raw(A, VariableIndex, Aprime),
derivate_raw(B, VariableIndex, Bprime).
......@@ -57,7 +62,7 @@ derivate_raw(A ^ B, VariableIndex, Result) :-
derivate_raw(sqrt(A), VariableIndex, Result) :-
!,
derivate_raw(A ^ v(0.5), VariableIndex, Result).
derivate_raw(A ^ 0.5, VariableIndex, Result).
derivate_raw(log(A), VariableIndex, Aprime / A) :-
!,
......
......@@ -2,17 +2,17 @@
:- begin_tests(formal_derivation).
test('dx/dt(x) = 1', [true(E == v(1))]) :-
derivate(x(0), 0, E).
test('dx/dt(x) = 1', [true(E == 1)]) :-
derivate([x], x, E).
test('dx/dt(y) = 0', [true(E == v(0))]) :-
derivate(x(1), 0, E).
test('dx/dt(y) = 0', [true(E == 0)]) :-
derivate([y], x, E).
test('dx/dt(x^2) = 2x', [true(E == v(2) * x(0))]) :-
derivate(x(0) ^ v(2), 0, E).
test('dx/dt(x^2) = 2x', [true(E == 2 * [x])]) :-
derivate([x] ^ 2, x, E).
test('dx/dt(cos(sqrt(x))) = - v(0.5) / sqrt(x) * sin(sqrt(x))',
[true(E == - (v(0.5) / sqrt(x(0))) * sin(sqrt(x(0))))]) :-
derivate(cos(sqrt(x(0))), 0, 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).
:- end_tests(formal_derivation).
:- module(
ode,
gsl,
[
solve/2
]).
......@@ -119,15 +119,17 @@ write_jacobian_cell(I, J, Expression) :-
write_arithmetic_expression(Expression),
write(');\n').
write_arithmetic_expression(x(VariableIndex)) :-
write_arithmetic_expression(Value) :-
number(Value),
!,
format('~d', [Value]).
write_arithmetic_expression([VariableIndex]) :-
format('x[~d]', [VariableIndex]).
write_arithmetic_expression(p(ParameterIndex)) :-
format('p[~d]', [ParameterIndex]).
write_arithmetic_expression(v(Value)) :-
format('~d', [Value]).
write_arithmetic_expression(- A) :-
write('- '),
write_arithmetic_expression_with_parentheses(A).
......@@ -257,7 +259,7 @@ print(FILE *csv, double t, double x[])
compile_c_program(ExecutableFilename) :-
call_subprocess(
path(gcc),
['-o', ExecutableFilename, 'ode_integrator.c',
['-o', ExecutableFilename, 'gsl_solver.c',
'-lgsl', '-lgslcblas', '-lm']).
......
:- use_module(library(plunit)).
:- begin_tests(ode).
:- begin_tests(gsl).
test('van_der_pol', [
true((
......@@ -12,7 +12,7 @@ test('van_der_pol', [
LastTimeStamp == 100.0
))]) :-
Options = [
equations: [x(1), -x(0) + p(0) * x(1) * (v(1) - x(0) ^ v(2))],
equations: [[1], -[0] + p(0) * [1] * (1 - [0] ^ 2)],
initial_values: [1.0, 0.0],
initial_parameter_values: [10],
method: gsl_odeiv2_step_rk8pd,
......@@ -20,7 +20,8 @@ test('van_der_pol', [
initial_step_size: 1e-6,
precision: 5,
time_initial: 0,
time_final: 100],
time_final: 100
],
solve(Options, Table).
:- end_tests(ode).
:- end_tests(gsl).
:- module(
initial_state,
[
get_initial_state/2
]
).
:- dynamic(initial_state/2).
get_initial_state(Molecule, State) :-
(
initial_state(Molecule, State)
->
true
;
State = 0
).
:- module(
kinetics,
[
kinetics/3
]
).
kinetics(Kinetics, Reactants, Expression) :-
(
Kinetics = 'MA'(Coefficient)
->
mass_action(Coefficient, Reactants, Expression)
).
mass_action(Coefficient, Reactants, Expression) :-
concentration_product(Reactants, ConcentrationProduct),
Expression = Coefficient * ConcentrationProduct.
concentration_product([], 1).
concentration_product([Coefficient * Reactant | Reactants], Product) :-
Product = Coefficient * [Reactant] * ProductTail,
concentration_product(Reactants, ProductTail).
......@@ -2,8 +2,9 @@
models,
[
add_item/3,
list_items/2,
remove_item/1
list_items/1,
remove_item/1,
item/2
]
).
......
:- module(
numerical_simulation,
[
numerical_simulation/0
]
).
numerical_simulation :-
biocham_command,
compute_ode,
make_ode_system,
solve.
make_ode_system :-
enumerate_variables,
convert_ode.
solve :-
gather_equations(Equations),
gather_initial_values(InitialValues),
gather_initial_parameter_values(InitialParameterValues),
Options = [
equations: Equations,
initial_values: InitialValues,
initial_parameter_values: InitialParameterValues,
method: gsl_odeiv2_step_rk8pd,
error_epsilon: 1e-6,
initial_step_size: 1e-6,
precision: 5,
time_initial: 0,
time_final: 100
],
solve(Options, Table),
print(Table).
gather_equations(Equations) :-
peek_count(variable_counter, VariableCount),
VariableMax is VariableCount - 1,
findall(
Equation,
(
between(0, VariableMax, VariableIndex),
equation(VariableIndex, Equation)
),
Equations
).
gather_initial_values(InitialValues) :-
peek_count(variable_counter, VariableCount),
VariableMax is VariableCount - 1,
findall(
InitialValue,
(
between(0, VariableMax, VariableIndex),
get_initial_state(VariableIndex, InitialValue)
),
InitialValues
).
gather_initial_parameter_values(InitialParameterValues) :-
peek_count(parameter_counter, ParameterCount),
ParameterMax is ParameterCount - 1,
findall(
InitialParameterValue,
(
between(0, ParameterMax, ParameterIndex),
parameter(ParameterIndex, Parameter),
get_parameter_value(Parameter, InitialParameterValue)
),
InitialParameterValues
).
:- dynamic(variable/2).
enumerate_variables :-
retractall(variable(_, _)),
set_counter(variable_counter, 0),
\+ (
ode(Molecule, _),
\+ (
count(variable_counter, VariableIndex),
assertz(variable(Molecule, VariableIndex))
)
).
:- dynamic(equation/2).
:- dynamic(parameter/2).
convert_ode :-
retractall(equation(_, _)),
set_counter(parameter_counter, 0),
\+ (
ode(Molecule, Expression),
\+ (
variable(Molecule, VariableIndex),
convert_expression(Expression, IndexedExpression),
assertz(equation(VariableIndex, IndexedExpression))
)
).
convert_expression(Value, Value) :-
number(Value),
!.
convert_expression([Molecule], x(VariableIndex)) :-
!,
variable(Molecule, VariableIndex).
convert_expression(Parameter, p(ParameterIndex)) :-
!,
(
parameter(Parameter, ParameterIndex)
->
true
;
count(parameter_counter, ParameterIndex),
assertz(parameter(Parameter, ParameterIndex))
).
convert_expression(- A, - AIndexed) :-
!,
convert_expression(A, AIndexed).
convert_expression(A + B, AIndexed + BIndexed) :-
!,
convert_expression(A, AIndexed),
convert_expression(B, BIndexed).
convert_expression(A - B, AIndexed - BIndexed) :-
!,
convert_expression(A, AIndexed),
convert_expression(B, BIndexed).
convert_expression(A * B, AIndexed * BIndexed) :-
!,
convert_expression(A, AIndexed),
convert_expression(B, BIndexed).
convert_expression(A / B, AIndexed / BIndexed) :-
!,
convert_expression(A, AIndexed),
convert_expression(B, BIndexed).
convert_expression(A ^ B, AIndexed ^ BIndexed) :-
!,
convert_expression(A, AIndexed),
convert_expression(B, BIndexed).
convert_expression(sqrt(A), sqrt(AIndexed)) :-
!,
convert_expression(A, AIndexed).
convert_expression(exp(A), exp(AIndexed)) :-
!,
convert_expression(A, AIndexed).
convert_expression(log(A), log(AIndexed)) :-
!,
convert_expression(A, AIndexed).
convert_expression(cos(A), cos(AIndexed)) :-
!,
convert_expression(A, AIndexed).
convert_expression(sin(A), sin(AIndexed)) :-
!,
convert_expression(A, AIndexed).
convert_expression(E, _Result) :-
type_error(arithmetic_expression, E).
:- module(
objects,
[
rule/1
]
).
rule(_ for _).
rule(_ => _).
rule(_ <=> _).
:- module(
ode,
[
list_ODE/0
list_ODE/0,
compute_ode/0,
ode/2
]).
:- dynamic(ode/2).
......@@ -11,20 +13,57 @@ list_ODE :-
compute_ode,
print_ode.
compute_ode :-
retractall(ode(_, _)),
\+ (
item([model: current_model, kind: rule], Item),
rule(Item, Kinetics, Reactants, Products)
rule(Item, Kinetics, Reactants, Products),
\+ (
kinetics(Kinetics, Reactants, KineticsExpression),
add_molecules(Reactants, -KineticsExpression),
add_molecules(Products, KineticsExpression)
)
),
simplify_ode.
simplify_ode :-
\+ (
retract(ode(Variable, Expression)),
\+ (
simplify(Expression, SimplifiedExpression),
assertz(ode(Variable, SimplifiedExpression))
)
).
add_molecules(Molecules, Kinetics) :-
\+ (
member(Molecule, Molecules),
\+ (
add_molecule(Molecule, Kinetics)
)
).
add_molecule(Coefficient * Molecule, Kinetics) :-
(
retract(ode(Molecule, Expression))
->
true
;
Expression = 0
),
asserta(ode(Molecule, Expression + Kinetics / Coefficient)).
print_ode :-
write('Init \tODE\n'),
\+ (
ode(Molecule, Equation),
initial_state(Molecule, InitialState)
get_initial_state(Molecule, InitialState),
\+ (
format('~f \td[~w]/dt=~w\n', Equation)
format('~f \td[~w]/dt=~w\n', [InitialState, Molecule, Equation])
)
).
......@@ -5,8 +5,9 @@
list_rules/0,
rule/4,
op(800, xfx, for),
op(600, xfx, =>),
op(600, xfx, <=>)
op(700, xfx, <=),
op(750, xfx, =>),
op(750, xfx, <=>)
]
).
......@@ -20,40 +21,41 @@ list_rules :-
rule(Item, Kinetics, Reactants, Products) :-
(
Item = PairKinetics for Body
Item = (PairKinetics for Body)
->
true
;
Kinetics = 'MA'(1)
)
PairKinetics = 'MA'(1),
Body = Item
),
(
Item = (Left => Right)
Body = (Left => Right)
->
Reversible = no,
Catalyst = '_'
;
Item = (Left <=> Right)
Body = (Left <=> Right)
->
Reversibe = yes,
Reversible = yes,
Catalyst = '_'
;
Item = (Left =[ Catalyst ]=> Right)
Body = (Left =[ Catalyst ]=> Right)
->
Reversible = no
;
Item = (Left <=[ Catalyst ]=> Right)
Body = (Left <=[ Catalyst ]=> Right)
->
Reversible = yes
),
solution(Left, LeftMolecules),
solution(Right, RightMolecules),
solution(Catalyst, CatalystMolecules),
append(Catalyst, LeftMolecules, LeftMoleculesWithCatalyst),
append(Catalyst, RightMolecules, RightMoleculesWithCatalyst),
append(CatalystMolecules, LeftMolecules, LeftMoleculesWithCatalyst),
append(CatalystMolecules, RightMolecules, RightMoleculesWithCatalyst),
(
Reversible = yes
->
pair_kinetics(PairKinetics, ForwardKinetics, BackwardKinetics)
pair_kinetics(PairKinetics, ForwardKinetics, BackwardKinetics),
(
Kinetics = ForwardKinetics,
Reactants = LeftMoleculesWithCatalyst,
......
......@@ -4,202 +4,240 @@
simplify/2
]).
simplify(x(VariableIndex), x(VariableIndex)).
simplify(Value, Value) :-
number(Value),
!.
simplify(Parameter, Parameter) :-
atom(Parameter),
!.
simplify(v(Value), v(Value)).
simplify(p(ParameterIndex), p(ParameterIndex)) :-
!.
simplify(p(ParameterIndex), p(ParameterIndex)).
simplify([Variable], [Variable]) :-
!.
simplify(- A, Result) :-
!,
simplify(A, Asimplified),
simplify_opposite(Asimplified, Result).
simplify(A + B, Result) :-
!,
simplify(A, Asimplified),
simplify(B, Bsimplified),
simplify_addition(Asimplified, Bsimplified, Result).
simplify(A - B, Result) :-
!,