Commit e72db7ba authored by MARTINEZ Thierry 's avatar MARTINEZ Thierry

search_parameters

parent fe0c42d5
[submodule "modules/c-cmaes"]
path = modules/c-cmaes
url = https://github.com/CMA-ES/c-cmaes.git
......@@ -14,7 +14,7 @@ $(foreach var, CC PLBASE PLCFLAGS PLLDFLAGS PLLIB, \
INCLUDEDIRS=modules/graphviz modules/sbml modules/partialfrac $(PLBASE)/include
CFLAGS=$(addprefix -I, $(INCLUDEDIRS)) $(PLCFLAGS)
CFLAGS=-g $(addprefix -I, $(INCLUDEDIRS)) $(PLCFLAGS)
LDFLAGS=$(PLLDFLAGS) $(addprefix -L, $(dir $(wildcard $(PLBASE)/lib/*/)))
......@@ -26,13 +26,16 @@ LDLIBS=$(PLLIB) `pkg-config --libs libgvc` $(LIBSBML) -lgsl -lgslcblas -lm
SWIPL=$(PWD)/swipl-biocham
CMAES_LIB=library/cmaes.c library/cmaes.h library/cmaes_interface.h
all: biocham biocham_debug test doc jupyter
quick: unit_tests
.PHONY: all slow test unit_tests doc clean zip
biocham: platform/current $(SWIPL) $(MODULES) toc.org Makefile
biocham: platform/current $(SWIPL) $(MODULES) toc.org $(CMAES_LIB) \
library/gsl_solver.o library/cmaes.o Makefile
$(SWIPL) -q -o biocham \
--goal="\
set_prolog_flag(verbose, normal), \
......@@ -129,6 +132,7 @@ clean:
- rm biocham_full_tests
- rm swipl-biocham
- rm swipl-biocham.o
- rm library/gsl_solver.o
biocham-src.tar.gz: $(MODULES) $(TEST_MODULES) toc.org Makefile
mkdir tmp/
......@@ -174,3 +178,20 @@ biocham-src-%.zip:
git archive --prefix=biocham/ -o $@ HEAD
zip: biocham-src-$(shell git rev-parse --short HEAD).zip
library/cmaes.c: modules/c-cmaes/src/cmaes.c
cp $< $@
library/cmaes.h: modules/c-cmaes/src/cmaes.h
cp $< $@
library/cmaes_interface.h: modules/c-cmaes/src/cmaes_interface.h
cp $< $@
modules/c-cmaes/src/cmaes.c \
modules/c-cmaes/src/cmaes.h \
modules/c-cmaes/src/cmaes_interface.h: modules/c-cmaes
modules/c-cmaes:
git submodule init
git submodule update
:- module(
biocham,
[
% Public API
start/0,
initialize/0,
initial/1,
biocham_command/0,
biocham_command/1,
library_path/1
......@@ -56,10 +58,27 @@ initialize :-
nb_setval(graph_pdf, 0),
nb_setval(current_models, []),
load_biocham('library:initial'),
\+ (
initial_command(Command),
\+ (
command(Command)
)
),
new_model,
assertz(initialized).
:- dynamic(initial_command/1).
initial(Command) :-
devdoc('
Adds a command to run during initialization
(in the context of the "initial" model)
'),
assertz(initial_command(Command)).
:- dynamic(library_path/1).
......
:- module(
c_compiler,
[
% Public API
compile_c_program/3,
compile_cpp_program/3,
gsl_compile_options/1,
ppl_compile_options/1
]
).
% Public API
compile_c_program(SourceFiles, Options, ExecutableFilename) :-
compile_program('gcc', SourceFiles, Options, ExecutableFilename).
compile_cpp_program(SourceFiles, Options, ExecutableFilename) :-
compile_program(
'g++', SourceFiles, ['--std=c++1z' | Options], ExecutableFilename
).
gsl_compile_options(['-lgsl', '-lgslcblas', '-lm']).
ppl_compile_options(['-lppl_c', '-lppl', '-lgmpxx', '-lgmp']).
% Private predicates
:- dynamic(debug/0).
debug.
compile_program(Compiler, SourceFiles, Options, ExecutableFilename) :-
library_path(LibraryPath),
findall(
AbsoluteSourceFile,
(
member(SourceFile, SourceFiles),
absolute_file_name(
SourceFile, AbsoluteSourceFile, [relative_to(LibraryPath)]
)
),
AbsoluteSourceFiles
),
append(AbsoluteSourceFiles, Options, AllOptions),
(
debug
->
AllOptions0 = ['-g' | AllOptions]
;
AllOptions0 = AllOptions
),
call_subprocess(
path(Compiler), ['-o', ExecutableFilename, '-I.' | AllOptions0], []
).
......@@ -436,6 +436,15 @@ format_doc_html(prompt(Prompt, Command, Output)) :-
[Prompt, EscapedCommand, EscapedOutput]
).
format_doc_html(command(Command)) :-
escape_html(Command, EscapedCommand),
format(
'</div>
<div><kbd>~w.</kbd></div>
<div>',
[EscapedCommand]
).
format_doc_html(symbol(backslash)) :-
format('\\', []).
......@@ -773,10 +782,17 @@ format_doc_tex(prompt(Prompt, Command, Output)) :-
escape_tex(Command, EscapedCommand),
escape_tex(Output, EscapedOutput),
format(
'\\textbf{\\texttt{~a.}}\\texttt{~a}',
'\\textbf{\\texttt{~a.}}\\texttt{~a}\\\\\\texttt{~a}',
[Prompt, EscapedCommand, EscapedOutput]
).
format_doc_tex(command(Command)) :-
escape_tex(Command, EscapedCommand),
format(
'\\texttt{~a}.',
[EscapedCommand]
).
format_doc_tex(symbol(backslash)) :-
format('\\(\\backslash\\)', []).
......@@ -1106,6 +1122,7 @@ generate_body_item(section(Level, Title), _Type) :-
generate_body_item(file(File), Type) :-
format(user_error, 'File: ~a\n', [File]),
atom_concat('../', File, ParentFile),
setup_call_cleanup(
open(ParentFile, read, Stream),
......@@ -1178,11 +1195,14 @@ generate_clause(Clause, Variables, VariableNames, Type) :-
Type = devdoc
;
Clause = (:- devcom(Comment)),
((atom_chars(Comment,[A1,A2|_]),'A'@=<A1,A1@=<'Z','A'@=<A2,A2@=<'Z')
->
atom_chars(Author,[A1,A2])
;
Author='gray'
(
atom_chars(Comment,[A1,A2|_]),
uppercase_char(A1),
uppercase_char(A2)
->
atom_chars(Author,[A1,A2])
;
Author='gray'
),
atom_concat('\\begin{quote}\\quoteAuthor{',Author,A),
atom_concat(A,'}',B),
......@@ -1223,6 +1243,11 @@ generate_clause(Clause, Variables, VariableNames, Type) :-
->
instantiate_grammar_body(Body),
format_doc(grammar_item(Item))
;
Clause = (:- initial(Command))
->
close_grammar,
format_doc(command(Command))
;
Clause = (:- _)
->
......@@ -1495,12 +1520,12 @@ write_doc_item(devdoc(DocBody), Type) :-
).
write_doc_item(biocham_silent(Command), _Type) :-
command(Command).
echo_command(Command).
write_doc_item(biocham(Command), _Type) :-
with_output_to(
atom(Output),
command(Command)
echo_command(Command)
),
prompt(Prompt),
format_doc(prompt(Prompt, Command, Output)).
......@@ -1599,8 +1624,13 @@ escape_tex_chars([H | T0], [H | T1]) :-
write_doc(Atom) :-
atom_chars(Atom, Chars),
write_doc_chars(Chars).
watch_error(
(
atom_chars(Atom, Chars),
doc:write_doc_chars(Chars)
),
format(user_error, 'Error while formatting: ~a\n', [Atom])
).
write_doc_chars([]).
......@@ -1928,6 +1958,13 @@ biocham_reexecute(Trace, Prompt) :-
).
echo_command(Command) :-
watch_error(
command(Command),
format(user_error, 'Error while executing ~p\n', [Command])
).
close_opened_li :-
(
nb_getval(opened_li, true)
......
This diff is collapsed.
:- module(
gsl,
[
solve/2
% Public API
solve/2,
write_gsl_file/2,
write_field/1
]).
......@@ -23,7 +26,7 @@ solve(Options, Table) :-
ExecutableFilename = 'ode',
TableFilename = 'ode.csv',
write_gsl_file(CFilename, Options),
compile_c_program(ExecutableFilename),
compile_gsl_c_program(ExecutableFilename),
call_subprocess(ExecutableFilename, [], []),
csv_read_file(TableFilename, Table),
delete_file(CFilename),
......@@ -31,9 +34,6 @@ solve(Options, Table) :-
delete_file(TableFilename).
:- devdoc('\\section{Private predicates}').
write_gsl_file(CFilename, Options) :-
setup_call_cleanup(
open(CFilename, write, Stream),
......@@ -41,6 +41,23 @@ write_gsl_file(CFilename, Options) :-
close(Stream)
).
write_field(t) :-
write(t).
write_field(x(VariableIndex)) :-
format('x[~d]', [VariableIndex]).
write_field(p(ParameterIndex)) :-
format('p[~d]', [ParameterIndex]).
write_field(expression(Expr)) :-
write_arithmetic_expression(Expr).
:- devdoc('\\section{Private predicates}').
write_gsl_stream(Stream, Options) :-
with_output_to(
Stream,
......@@ -48,7 +65,7 @@ write_gsl_stream(Stream, Options) :-
).
write_gsl(Options) :-
write_headers(Options),
write_config(Options),
write_equations(Options),
write_jacobian(Options),
write_initial_values(Options),
......@@ -57,7 +74,7 @@ write_gsl(Options) :-
write_print_headers(Options),
write_print(Options).
write_headers(Options) :-
write_config(Options) :-
memberchk(equations: Equations, Options),
length(Equations, VariableCount),
memberchk(initial_parameter_values: InitialParameterValues, Options),
......@@ -69,25 +86,31 @@ write_headers(Options) :-
;
EventCount = 0
),
write_header('VARIABLE_COUNT', '~d', VariableCount),
write_header('PARAMETER_COUNT', '~d', ParameterCount),
write_header('EVENT_COUNT', '~d', EventCount),
write_option('METHOD', 'gsl_odeiv2_step_~a', method, Options),
write_option('ERROR_EPSILON_ABSOLUTE', '~g', error_epsilon_absolute, Options),
write_option('ERROR_EPSILON_RELATIVE', '~g', error_epsilon_relative, Options),
write_option('INITIAL_STEP_SIZE', '~g', initial_step_size, Options),
write_option('MAXIMUM_STEP_SIZE', '~g', maximum_step_size, Options),
write_option('TIME_INITIAL', '~g', time_initial, Options),
write_option('TIME_FINAL', '~g', time_final, Options).
write_header(Variable, Format, Value) :-
format('#define ~a ', [Variable]),
write('
void
init_gsl_solver_config(struct gsl_solver_config *config)
{
'),
write_config_field('variable_count', '~d', VariableCount),
write_config_field('parameter_count', '~d', ParameterCount),
write_config_field('event_count', '~d', EventCount),
write_option('method', 'gsl_odeiv2_step_~a', method, Options),
write_option('error_epsilon_absolute', '~g', error_epsilon_absolute, Options),
write_option('error_epsilon_relative', '~g', error_epsilon_relative, Options),
write_option('initial_step_size', '~g', initial_step_size, Options),
write_option('maximum_step_size', '~g', maximum_step_size, Options),
write_option('time_initial', '~g', time_initial, Options),
write_option('time_final', '~g', time_final, Options),
write('}\n').
write_config_field(Variable, Format, Value) :-
format(' config->~a = ', [Variable]),
format(Format, [Value]),
nl.
write(';\n').
write_option(Variable, Format, Option, Options) :-
memberchk(Option: Value, Options),
write_header(Variable, Format, Value).
write_config_field(Variable, Format, Value).
write_equations(Options) :-
memberchk(equations: Equations, Options),
......@@ -95,7 +118,7 @@ write_equations(Options) :-
int
functions(double t, const double x[], double f[], void *data)
{
double *p = data;
double *p = (double *) data;
'),
\+ (
nth0(I, Equations, Equation),
......@@ -122,7 +145,7 @@ write_jacobian(Options) :-
int
jacobian(double t, const double x[], double *dfdx, double dfdt[], void *data)
{
double *p = data;
double *p = (double *) data;
gsl_matrix_view dfdx_mat = gsl_matrix_view_array(dfdx, ~d, ~d);
gsl_matrix *m = &dfdx_mat.matrix;
', [VariableCount, VariableCount]),
......@@ -284,7 +307,7 @@ initial_parameter_values(double p[])
\+ (
nth0(ParameterIndex, InitialParameterValues, ParameterValue),
\+ (
format(' p[~d] = ~g;\n', [ParameterIndex, ParameterValue])
format(' p[~d] = ~g;\n', [ParameterIndex, ParameterValue])
)
),
write('\c
......@@ -322,19 +345,25 @@ write_events(Options) :-
),
write('
void
events(int e[], double p[], double x[], gsl_odeiv2_driver* d, double t,
double* t_upper)
events(struct gsl_solver_state *state, double *t_upper)
{
double f[VARIABLE_COUNT];
const struct gsl_solver *solver = state->solver;
const struct gsl_solver_config *config = solver->config;
gsl_odeiv2_driver *d = solver->d;
double f[config->variable_count];
double t = state->t;
double *p = solver->p;
double *x = state->x;
bool *e = state->e;
'),
\+ (
nth0(EventIndex, Events, (Condition -> ParameterValues)),
\+ (
format(' if (!e[~d] && (', [EventIndex]),
format(' if (!e[~d] && (', [EventIndex]),
write_condition(Condition),
format(
')) {
e[~d] = 1;',
e[~d] = true;',
[EventIndex]
),
\+ (
......@@ -342,7 +371,7 @@ events(int e[], double p[], double x[], gsl_odeiv2_driver* d, double t,
\+ (
format(
'
p[~d] = ',
p[~d] = ',
[ParameterIndex]
),
write_arithmetic_expression(Value),
......@@ -351,16 +380,16 @@ events(int e[], double p[], double x[], gsl_odeiv2_driver* d, double t,
),
format(
'
gsl_odeiv2_driver_reset(d);
}
else if (e[~d] && !(',
gsl_odeiv2_driver_reset(d);
}
else if (e[~d] && !(',
[EventIndex]
),
write_condition(Condition),
format(
')) {
e[~d] = 0;
}
e[~d] = false;
}
',
[EventIndex]
),
......@@ -370,16 +399,16 @@ events(int e[], double p[], double x[], gsl_odeiv2_driver* d, double t,
% \cite{EKP01hscc}
is_linear(Condition, Derivatives, Options, Normalized)
->
format(' else if (!e[~d] && !(', [EventIndex]),
format(' else if (!e[~d] && !(', [EventIndex]),
write_condition(Condition),
write(')) {
functions(t, x, f, p);
*t_upper = fmax(t + INITIAL_STEP_SIZE, fmin(*t_upper, t + fabs('),
functions(t, x, f, p);
*t_upper = fmax(t + config->initial_step_size, fmin(*t_upper, t + fabs('),
write_arithmetic_expression(Normalized),
write(')/'),
write_arithmetic_expression(Derivatives),
write('));
}
}
')
;
true
......@@ -392,103 +421,48 @@ events(int e[], double p[], double x[], gsl_odeiv2_driver* d, double t,
write_print_headers(Options) :-
(
memberchk(nonconstant_parameters: NCParameters, Options)
->
true
;
NCParameters = []
),
(
memberchk(nonparametric_functions: Functions, Options)
->
true
;
Functions = []
),
(
memberchk(headers: Headers, Options)
->
write('
memberchk(fields: [FirstField: _ | Fields], Options),
format('
void
print_headers(FILE *csv)
{
fprintf(csv, "#Time'),
fprintf(csv, "#\\"~a\\"', [FirstField]),
\+ (
member(Field: _, Fields),
\+ (
(
member(Header, Headers)
;
member((Header, _), NCParameters)
;
member((Header, _), Functions)
),
\+ (
format(',\\"~a\\"', [Header])
)
),
write('\\n"'),
write(');
}
')
;
write('
void
print_headers(FILE *csv)
{
format(',\\"~a\\"', [Field])
)
),
write('\\n");
}
')
).
').
write_print(Options) :-
memberchk(equations: Equations, Options),
length(Equations, VariableCount),
VariableMax is VariableCount - 1,
memberchk(precision: Precision, Options),
(
memberchk(nonconstant_parameters: NCParameters, Options)
->
true
;
NCParameters = []
),
(
memberchk(nonparametric_functions: Functions, Options)
->
true
;
Functions = []
),
length(NCParameters, NCParamCount),
length(Functions, FunctionsCount),
PlotCount is VariableCount + NCParamCount + FunctionsCount,
memberchk(fields: Fields, Options),
format('
void
print(FILE *csv, double t, double x[], double p[])
print(FILE *csv, struct gsl_solver_state *state)
{
fprintf(csv, "%.~dg', [Precision]),
double t = state->t;
double *p = state->solver->p;
double *x = state->x;
bool *e = state->e;
fprintf(csv, "%.~dg', [Precision]),
\+ (
between(1, PlotCount, _),
Fields = [_First | Others],
member(_Field, Others),
\+ (
format(',%.~dg', [Precision])
)
),
write('\\n", t'),
write('\\n"'),
\+ (
between(0, VariableMax, VariableIndex),
member(_Header: Field, Fields),
\+ (
format(', x[~d]', [VariableIndex])
)
),
forall(
member((_, ParamIndex), NCParameters),
format(', p[~d]', [ParamIndex])
),
forall(
member((_, Expr), Functions),
(
write(', '),
write_arithmetic_expression(Expr)
format(', '),
write_field(Field)
)
),
write(');
......@@ -496,15 +470,10 @@ print(FILE *csv, double t, double x[], double p[])
').
compile_c_program(ExecutableFilename) :-
library_path(LibraryPath),
absolute_file_name('gsl_solver.c', GslSolverC, [relative_to(LibraryPath)]),
call_subprocess(
path(gcc),
['-o', ExecutableFilename, '-I.', GslSolverC,
'-lgsl', '-lgslcblas', '-lm'],
[]
).
compile_gsl_c_program(ExecutableFilename) :-
gsl_compile_options(GslCompileOptions),
compile_c_program(
['gsl_loop.c', 'gsl_solver.o'], GslCompileOptions, ExecutableFilename).
% FIXME and=max and or=min but after normalization
......
......@@ -4,6 +4,7 @@
test('van_der_pol', [condition(flag(slow_test, true, true))]) :-
Options = [
fields: ['Time': t, 'x0': x(0), 'x1': x(1)],
equations: [[1], -[0] + p(0) * [1] * (1 - [0] ^ 2)],
initial_values: [1.0, 0.0],
initial_parameter_values: [10],
......
......@@ -103,6 +103,9 @@ check_influence_model :-
).
:- meta_predicate with_influence_model(0).
with_influence_model(G) :-
(
is_influence_model
......
......@@ -53,3 +53,18 @@ make_product([RS * RO | Tail], S, O, Expression, Product) :-
substitute([S, O], [RS, RO], Expression, SubstitutedExpression),
Product = SubstitutedExpression * TailProduct,
make_product(Tail, S, O, Expression, TailProduct).
% FF strange to me and not robust for summing kinetics
:- initial('function MA(k) = k * product(S * O in [reactants], O ^ S)').
:- initial('
function MM(Vm, Km) = Vm * single_reactant / (Km + single_reactant)
').
:- initial('
function H(Vm, Km, n) =
Vm * single_reactant ^ n / (Km ^ n + single_reactant ^ n).
').
#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>
#include <math.h>
#include "gsl_solver.h"
#include "ode.inc"
int
main(void)
{
FILE *csv = fopen("ode.csv", "w");
struct gsl_solver_config config;
struct gsl_solver solver;
struct gsl_solver_state state;
init_gsl_solver_config(&config);
init_gsl_solver(&solver, &config);
init_gsl_solver_state(&state, &solver);
print_headers(csv);
do {
print(csv, &state);
} while (gsl_solver_step(&state));
free_gsl_solver_state(&state);
free_gsl_solver(&solver);
return EXIT_SUCCESS;
}
#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
#include "gsl_solver.h"
#include "ode.inc"
void *
xcalloc(size_t count, size_t size)
{
void *result = calloc(count, size);
if (!result) {
fprintf(stderr, "Memory failure\n");
exit(EXIT_FAILURE);
}
return result;
}
int
main(void)
void
init_gsl_solver(
struct gsl_solver *solver, const struct gsl_solver_config *config)
{
FILE *csv = fopen("ode.csv", "w");
double p[PARAMETER_COUNT];
const gsl_odeiv2_step_type *T = METHOD;
gsl_odeiv2_step *s = gsl_odeiv2_step_alloc(T, VARIABLE_COUNT);
gsl_odeiv2_system sys = {functions, jacobian, VARIABLE_COUNT, p};
const gsl_odeiv2_step_type *method = config->method;
int variable_count = config->variable_count;
double *p = xcalloc(config->parameter_count, sizeof(double));
gsl_odeiv2_system sys = {functions, jacobian, variable_count, p};
solver->sys = sys;
gsl_odeiv2_step *s = gsl_odeiv2_step_alloc(method, variable_count);
gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new(
&sys, METHOD, INITIAL_STEP_SIZE, ERROR_EPSILON_ABSOLUTE,
ERROR_EPSILON_RELATIVE);
double t = TIME_INITIAL;
double delta = TIME_FINAL - TIME_INITIAL;
double maximum_step_size = delta * MAXIMUM_STEP_SIZE;
double x[VARIABLE_COUNT];
int e[EVENT_COUNT];