Commit 7573596d authored by Thierry Martinez's avatar Thierry Martinez

GSL Solver

parents
MODULES= \
biocham.pl \
about.pl \
toplevel.pl \
models.pl \
rules.pl \
ode.pl \
util.pl \
simplify_arithmetic.pl \
formal_derivation.pl \
formal_derivation.plt \
gsl.pl \
gsl.plt \
counters.pl
all: biocham biocham_debug tests
biocham: $(MODULES) Makefile
swipl -o biocham --goal=about --toplevel=toplevel -c $(MODULES)
biocham_debug: $(MODULES) Makefile
swipl -o biocham_debug -c $(MODULES)
tests: biocham_tests
./biocham_tests
biocham_tests: $(MODULES) Makefile
swipl -o biocham_tests \
--goal="call_cleanup((run_tests, halt(0)), halt(1))" -c $(MODULES)
:- module(
about,
[
version/1,
copyright/1,
license/1,
url/1,
about/0
]).
version('4.0').
copyright(
'Copyright (C) 2003-2015 Inria, EPI Lifeware, Paris-Rocquencourt, France'
).
license('GNU GPL 2').
url('http://lifeware.inria.fr/biocham/').
about :-
biocham_command,
version(Version),
copyright(Copyright),
license(License),
url(URL),
format('Biocham ~a\n', [Version]),
format('~a,\n', [Copyright]),
format('license ~a, ~a\n', [License, URL]).
:- module(
biocham,
[
biocham_command/0,
biocham_command/1
]).
biocham_command.
biocham_command(Functor/Arity) :-
functor(Head, Functor, Arity),
once(clause(Head, (biocham_command, _))).
:- module(
counters,
[
set_counter/2,
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).
:- module(
formal_derivation,
[
derivate/3,
jacobian/2
]).
:- use_module(library(error)).
derivate(Expression, VariableIndex, Result) :-
derivate_raw(Expression, VariableIndex, UnsimplifiedResult),
simplify(UnsimplifiedResult, Result).
derivate_raw(x(VariableIndex0), VariableIndex1, Result) :-
!,
(
VariableIndex0 = VariableIndex1
->
Result = v(1)
;
Result = v(0)
).
derivate_raw(p(_Value), _VariableIndex, v(0)) :-
!.
derivate_raw(v(_Value), _VariableIndex, v(0)) :-
!.
derivate_raw(- A, VariableIndex, - Aprime) :-
!,
derivate_raw(A, VariableIndex, Aprime).
derivate_raw(A + B, VariableIndex, Aprime + Bprime) :-
!,
derivate_raw(A, VariableIndex, Aprime),
derivate_raw(B, VariableIndex, Bprime).
derivate_raw(A - B, VariableIndex, Aprime - Bprime) :-
!,
derivate_raw(A, VariableIndex, Aprime),
derivate_raw(B, VariableIndex, Bprime).
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, VariableIndex, Aprime),
derivate_raw(B, VariableIndex, Bprime).
derivate_raw(A ^ B, VariableIndex, Result) :-
!,
derivate_raw(exp(B * log(A)), VariableIndex, Result).
derivate_raw(sqrt(A), VariableIndex, Result) :-
!,
derivate_raw(A ^ v(0.5), VariableIndex, Result).
derivate_raw(log(A), VariableIndex, Aprime / A) :-
!,
derivate_raw(A, VariableIndex, Aprime).
derivate_raw(exp(A), VariableIndex, Aprime * exp(A)) :-
!,
derivate_raw(A, VariableIndex, Aprime).
derivate_raw(sin(A), VariableIndex, Aprime * cos(A)) :-
!,
derivate_raw(A, VariableIndex, Aprime).
derivate_raw(cos(A), VariableIndex, - Aprime * sin(A)) :-
!,
derivate_raw(A, VariableIndex, Aprime).
derivate_raw(E, _VariableIndex, _Aprime) :-
type_error(arithmetic_expression, E).
jacobian(Equations, Jacobian) :-
length(Equations, VariableCount),
findall(
Line,
(
member(Equation, Equations),
jacobian_line(Equation, VariableCount, Line)
),
Jacobian).
jacobian_line(Equation, VariableCount, Line) :-
VariableMax is VariableCount - 1,
findall(
Derivative,
(
between(0, VariableMax, VariableIndex),
derivate(Equation, VariableIndex, Derivative)
),
Line).
:- use_module(library(plunit)).
:- begin_tests(formal_derivation).
test('dx/dt(x) = 1', [true(E == v(1))]) :-
derivate(x(0), 0, E).
test('dx/dt(y) = 0', [true(E == v(0))]) :-
derivate(x(1), 0, E).
test('dx/dt(x^2) = 2x', [true(E == v(2) * x(0))]) :-
derivate(x(0) ^ v(2), 0, 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).
:- end_tests(formal_derivation).
:- module(
ode,
[
solve/2
]).
:- use_module(library(csv)).
solve(Options, Table) :-
\+ memberchk(jacobian: _, Options),
!,
memberchk(equations: Equations, Options),
jacobian(Equations, Jacobian),
solve([jacobian: Jacobian | Options], Table).
solve(Options, Table) :-
CFilename = 'ode.inc',
ExecutableFilename = 'ode',
write_gsl_file(CFilename, Options),
compile_c_program(ExecutableFilename),
call_subprocess(ExecutableFilename, []),
csv_read_file('ode.csv', Table).
write_gsl_file(CFilename, Options) :-
setup_call_cleanup(
open(CFilename, write, Stream),
write_gsl_stream(Stream, Options),
close(Stream)
).
write_gsl_stream(Stream, Options) :-
with_output_to(
Stream,
write_gsl(Options)
).
write_gsl(Options) :-
write_headers(Options),
write_equations(Options),
write_jacobian(Options),
write_initial_values(Options),
write_initial_parameter_values(Options),
write_print(Options).
write_headers(Options) :-
memberchk(equations: Equations, Options),
length(Equations, VariableCount),
memberchk(initial_parameter_values: InitialParameterValues, Options),
length(InitialParameterValues, ParameterCount),
write_header('VARIABLE_COUNT', '~d', VariableCount),
write_header('PARAMETER_COUNT', '~d', ParameterCount),
write_option('METHOD', '~a', method, Options),
write_option('ERROR_EPSILON', '~f', error_epsilon, Options),
write_option('INITIAL_STEP_SIZE', '~f', initial_step_size, Options),
write_option('TIME_INITIAL', '~f', time_initial, Options),
write_option('TIME_FINAL', '~f', time_final, Options).
write_header(Variable, Format, Value) :-
format('#define ~a ', [Variable]),
format(Format, [Value]),
nl.
write_option(Variable, Format, Option, Options) :-
memberchk(Option: Value, Options),
write_header(Variable, Format, Value).
write_equations(Options) :-
memberchk(equations: Equations, Options),
write('
int
functions(double t, const double x[], double f[], void *data)
{
double *p = data;
'),
\+ (
nth0(I, Equations, Equation),
\+ (
write_equation(I, Equation)
)
),
write(' \c
return GSL_SUCCESS;
}
').
write_equation(I, Equation) :-
format(' f[~d] = ', [I]),
write_arithmetic_expression(Equation),
write(';\n').
write_jacobian(Options) :-
memberchk(jacobian: Jacobian, Options),
length(Jacobian, VariableCount),
format('
int
jacobian(double t, const double x[], double *dfdx, double dfdt[], void *data)
{
double *p = data;
gsl_matrix_view dfdx_mat = gsl_matrix_view_array(dfdx, ~d, ~d);
gsl_matrix *m = &dfdx_mat.matrix;
', [VariableCount, VariableCount]),
\+ (
nth0(I, Jacobian, JacobianLine),
nth0(J, JacobianLine, Expression),
\+ (
write_jacobian_cell(I, J, Expression)
)
),
write(' \c
dfdt[0] = 0.0;
dfdt[1] = 0.0;
return GSL_SUCCESS;
}
').
write_jacobian_cell(I, J, Expression) :-
format(' gsl_matrix_set(m, ~d, ~d, ', [I, J]),
write_arithmetic_expression(Expression),
write(');\n').
write_arithmetic_expression(x(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).
write_arithmetic_expression(A + B) :-
write_arithmetic_expression_with_parentheses(A),
write(' + '),
write_arithmetic_expression_with_parentheses(B).
write_arithmetic_expression(A - B) :-
write_arithmetic_expression_with_parentheses(A),
write(' - '),
write_arithmetic_expression_with_parentheses(B).
write_arithmetic_expression(A * B) :-
write_arithmetic_expression_with_parentheses(A),
write(' * '),
write_arithmetic_expression_with_parentheses(B).
write_arithmetic_expression(A / B) :-
write_arithmetic_expression_with_parentheses(A),
write(' / '),
write_arithmetic_expression_with_parentheses(B).
write_arithmetic_expression(A ^ B) :-
write('pow('),
write_arithmetic_expression(A),
write(', '),
write_arithmetic_expression(B),
write(')').
write_arithmetic_expression(sqrt(A)) :-
write('sqrt('),
write_arithmetic_expression(A),
write(')').
write_arithmetic_expression(cos(A)) :-
write('cos('),
write_arithmetic_expression(A),
write(')').
write_arithmetic_expression(sin(A)) :-
write('sin('),
write_arithmetic_expression(A),
write(')').
write_arithmetic_expression(exp(A)) :-
write('exp('),
write_arithmetic_expression(A),
write(')').
write_arithmetic_expression(log(A)) :-
write('log('),
write_arithmetic_expression(A),
write(')').
write_arithmetic_expression_with_parentheses(Expression) :-
write('('),
write_arithmetic_expression(Expression),
write(')').
write_initial_values(Options) :-
memberchk(initial_values: InitialValues, Options),
write('
void
initial_values(double x[])
{
'),
\+ (
nth0(VariableIndex, InitialValues, InitialValue),
\+ (
format(' x[~d] = ~f;\n', [VariableIndex, InitialValue])
)
),
write('\c
}
').
write_initial_parameter_values(Options) :-
memberchk(initial_parameter_values: InitialParameterValues, Options),
write('
void
initial_parameter_values(double p[])
{
'),
\+ (
nth0(ParameterIndex, InitialParameterValues, ParameterValue),
\+ (
format(' p[~d] = ~f;\n', [ParameterIndex, ParameterValue])
)
),
write('\c
}
').
write_print(Options) :-
memberchk(equations: Equations, Options),
length(Equations, VariableCount),
VariableMax is VariableCount - 1,
memberchk(precision: Precision, Options),
format('
void
print(FILE *csv, double t, double x[])
{
fprintf(csv, "%.~de', [Precision]),
\+ (
between(1, VariableCount, _),
\+ (
format(',%.~de', [Precision])
)
),
write('\\n", t'),
\+ (
between(0, VariableMax, VariableIndex),
\+ (
format(', x[~d]', [VariableIndex])
)
),
write(');
}
').
compile_c_program(ExecutableFilename) :-
call_subprocess(
path(gcc),
['-o', ExecutableFilename, 'ode_integrator.c',
'-lgsl', '-lgslcblas', '-lm']).
call_subprocess(ExecutableFilename, Arguments) :-
process_create(ExecutableFilename, Arguments, [process(Pid)]),
process_wait(Pid, exit(0)).
:- use_module(library(plunit)).
:- begin_tests(ode).
test('van_der_pol', [
true((
Table = [FirstRow | OtherRows],
append(_, [LastRow], OtherRows),
FirstRow = row(FirstTimeStamp, _, _),
LastRow = row(LastTimeStamp, _, _),
FirstTimeStamp < 1e-5,
LastTimeStamp == 100.0
))]) :-
Options = [
equations: [x(1), -x(0) + p(0) * x(1) * (v(1) - x(0) ^ v(2))],
initial_values: [1.0, 0.0],
initial_parameter_values: [10],
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).
:- end_tests(ode).
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
#include "ode.inc"
int
main(void)
{
FILE *csv = fopen("ode.csv", "w");
double parameters[PARAMETER_COUNT];
const gsl_odeiv2_step_type *T = METHOD;
gsl_odeiv2_step *s = gsl_odeiv2_step_alloc(T, VARIABLE_COUNT);
gsl_odeiv2_control *c = gsl_odeiv2_control_y_new(ERROR_EPSILON, 0.0);
gsl_odeiv2_evolve * e = gsl_odeiv2_evolve_alloc(VARIABLE_COUNT);
gsl_odeiv2_system sys = {functions, jacobian, VARIABLE_COUNT, parameters};
double t = TIME_INITIAL;
double h = INITIAL_STEP_SIZE;
double x[VARIABLE_COUNT];
initial_values(x);
initial_parameter_values(parameters);
while (t < TIME_FINAL) {
int status = gsl_odeiv2_evolve_apply(
e, c, s, &sys, &t, TIME_FINAL, &h, x);
if (status != GSL_SUCCESS) {
fprintf(stderr, "error, return value=%d\n", status);
return EXIT_FAILURE;
}
print(csv, t, x);
}
gsl_odeiv2_evolve_free(e);
gsl_odeiv2_control_free(c);
gsl_odeiv2_step_free(s);
return EXIT_SUCCESS;
}
:- module(
models,
[
add_item/3,
list_items/2,
remove_item/1
]
).
add_item(Model, Kind, Item) :-
create_item_id(Id),
assertz(item(Id, Model, Kind, Item)).
item(Options, Item) :-
optional(model: Model, Options),
optional(kind: Kind, Options),
optional(id: Id, Options),
item(Id, Model, Kind, Item).
optional(Item, List) :-
(
memberchk(Item, List)
->
true
;
true
).
list_items(Options) :-
set_counter(list_item_counter, 0),
\+ (
item(Options, Item),
count(list_item_counter, Counter),
\+ (
format('[~d] ~w\n', [Counter, Item])
)
).
remove_item(Id) :-
retract(item(Id, _Model, _Kind, _Item)).
create_item_id(Id) :-
count(item_id, Id).
:- module(
ode,
[
list_ODE/0
]).
:- dynamic(ode/2).
list_ODE :-
biocham_command,
compute_ode,
print_ode.
compute_ode :-
retractall(ode(_, _)),
\+ (
item([model: current_model, kind: rule], Item),
rule(Item, Kinetics, Reactants, Products)
\+ (
print_ode :-
write('Init \tODE\n'),
\+ (
ode(Molecule, Equation),
initial_state(Molecule, InitialState)
\+ (
format('~f \td[~w]/dt=~w\n', Equation)
)
).
:- module(
rules,
[
add_rule/1,
list_rules/0,
rule/4,
op(800, xfx, for),
op(600, xfx, =>),
op(600, xfx, <=>)
]
).
add_rule(Rule) :-
biocham_command,
add_item(current_model, rule, Rule).
list_rules :-
biocham_command,
list_items([model: current_model, kind: rule]).
rule(Item, Kinetics, Reactants, Products) :-
(
Item = PairKinetics for Body
->
true
;
Kinetics = 'MA'(1)
)
(
Item = (Left => Right)
->
Reversible = no,
Catalyst = '_'
;
Item = (Left <=> Right)
->
Reversibe = yes,
Catalyst = '_'
;
Item = (Left =[ Catalyst ]=> Right)
->
Reversible = no
;
Item = (Left <=[ Catalyst ]=> Right)
->
Reversible = yes
),
solution(Left, LeftMolecules),
solution(Right, RightMolecules),
solution(Catalyst, CatalystMolecules),
append(Catalyst, LeftMolecules, LeftMoleculesWithCatalyst),
append(Catalyst, RightMolecules, RightMoleculesWithCatalyst),
(
Reversible = yes
->
pair_kinetics(PairKinetics, ForwardKinetics, BackwardKinetics)
(
Kinetics = ForwardKinetics,
Reactants = LeftMoleculesWithCatalyst,
Products = RightMoleculesWithCatalyst
;
Kinetics = BackwardKinetics,
Reactants = RightMoleculesWithCatalyst,
Products = LeftMoleculesWithCatalyst
)
;
Kinetics = PairKinetics,
Reactants = LeftMoleculesWithCatalyst,
Products = RightMoleculesWithCatalyst
).
pair_kinetics(PairKinetics, ForwardKinetics, BackwardKinetics) :-
(
PairKinetics = (ForwardKinetics, BackwardKinetics)
->
true
;
ForwardKinetics = PairKinetics,
BackwardKinetics = PairKinetics
).
solution('_', []) :-
!.
solution(A + B, Solution) :-
!,
solution(A, SolutionA),
solution(B, SolutionB),
append(SolutionA, SolutionB, Solution).
solution(Coefficient * Object, Solution) :-
!,
Solution = [Coefficient * Object].
solution(Object, [1 * Object]).
:- module(
arithmetic_simplify,
[
simplify/2
]).
simplify(x(VariableIndex), x(VariableIndex)).
simplify(v(Value), v(Value)).
simplify(p(ParameterIndex), p(ParameterIndex)).
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).