Commit 01691014 authored by Thierry Martinez's avatar Thierry Martinez

Conditional and events

parent cc8c520c
:- module(
events,
[
% Public API
add_event/2,
list_events/0
list_events/0,
list_model_events/0
]
).
:- devdoc('\\section{Commands}').
add_event(Condition, ParameterValues) :-
biocham_command(*),
type(Condition, condition),
......@@ -30,3 +35,19 @@ list_events :-
biocham_command,
doc('lists all the declared events.'),
list_items([kind: event]).
:- devdoc('\\section{Private predicates}').
list_model_events :-
devdoc('
lists all the events in a loadable syntax
(auxiliary predicate of list_model).
'),
\+ (
item([no_inheritance, kind: event, item: event(Condition, ParameterValue)]),
\+ (
format('add_event(~w, ~w).\n', [Condition, ParameterValue])
)
).
......@@ -8,6 +8,9 @@
:- use_module(library(csv)).
:- devdoc('\\section{Public API}').
solve(Options, Table) :-
\+ memberchk(jacobian: _, Options),
!,
......@@ -24,6 +27,9 @@ solve(Options, Table) :-
csv_read_file('ode.csv', Table).
:- devdoc('\\section{Private predicates}').
write_gsl_file(CFilename, Options) :-
setup_call_cleanup(
open(CFilename, write, Stream),
......@@ -43,6 +49,7 @@ write_gsl(Options) :-
write_jacobian(Options),
write_initial_values(Options),
write_initial_parameter_values(Options),
write_events(Options),
write_print_headers(Options),
write_print(Options).
......@@ -51,8 +58,16 @@ write_headers(Options) :-
length(Equations, VariableCount),
memberchk(initial_parameter_values: InitialParameterValues, Options),
length(InitialParameterValues, ParameterCount),
(
memberchk(events: Events, Options)
->
length(Events, EventCount)
;
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', '~f', error_epsilon_absolute, Options),
write_option('ERROR_EPSILON_RELATIVE', '~f', error_epsilon_relative, Options),
......@@ -88,11 +103,13 @@ functions(double t, const double x[], double f[], void *data)
}
').
write_equation(I, Equation) :-
format(' f[~d] = ', [I]),
write_arithmetic_expression(Equation),
write(';\n').
write_jacobian(Options) :-
memberchk(jacobian: Jacobian, Options),
length(Jacobian, VariableCount),
......@@ -118,6 +135,7 @@ jacobian(double t, const double x[], double *dfdx, double dfdt[], void *data)
}
').
write_jacobian_cell(I, J, Expression) :-
format(' gsl_matrix_set(m, ~d, ~d, ', [I, J]),
write_arithmetic_expression(Expression),
......@@ -196,6 +214,58 @@ write_arithmetic_expression_with_parentheses(Expression) :-
write(')').
write_condition(A = B) :-
write_arithmetic_expression_with_parentheses(A),
write(' == '),
write_arithmetic_expression_with_parentheses(B).
write_condition(A < B) :-
write_arithmetic_expression_with_parentheses(A),
write(' < '),
write_arithmetic_expression_with_parentheses(B).
write_condition(A > B) :-
write_arithmetic_expression_with_parentheses(A),
write(' > '),
write_arithmetic_expression_with_parentheses(B).
write_condition(A <= B) :-
write_arithmetic_expression_with_parentheses(A),
write(' <= '),
write_arithmetic_expression_with_parentheses(B).
write_condition(A >= B) :-
write_arithmetic_expression_with_parentheses(A),
write(' >= '),
write_arithmetic_expression_with_parentheses(B).
write_condition(true) :-
write('true').
write_condition(false) :-
write('false').
write_condition(A or B) :-
write_condition_with_parentheses(A),
write(' || '),
write_condition_with_parentheses(B).
write_condition(A and B) :-
write_condition_with_parentheses(A),
write(' && '),
write_condition_with_parentheses(B).
write_condition(not(A)) :-
write('!'),
write_condition_with_parentheses(A).
write_condition_with_parentheses(Condition) :-
write('('),
write_condition(Condition),
write(')').
write_initial_values(Options) :-
memberchk(initial_values: InitialValues, Options),
write('
......@@ -232,6 +302,43 @@ initial_parameter_values(double p[])
').
write_events(Options) :-
(
memberchk(events: Events, Options)
->
true
;
Events = []
),
write('
void
events(double e[], double p[])
{
'),
\+ (
nth0(EventIndex, Events, (Condition -> ParameterIndex = Value)),
\+ (
format(' if (!e[~d] && (', [EventIndex]),
write_condition(Condition),
format(')) {
e[~d] = 1;
p[~d] = ', [EventIndex, ParameterIndex]),
write_arithmetic_expression(Value),
format(';
}
else if (e[~d] && !(', [EventIndex]),
write_condition(Condition),
format(')) {
e[~d] = 0;
}
')
)
),
write('\c
}
').
write_print_headers(Options) :-
(
memberchk(headers: Headers, Options)
......
#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
......@@ -11,21 +12,27 @@ int
main(void)
{
FILE *csv = fopen("ode.csv", "w");
double parameters[PARAMETER_COUNT];
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, parameters};
gsl_odeiv2_system sys = {functions, jacobian, VARIABLE_COUNT, p};
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 x[VARIABLE_COUNT];
int e[EVENT_COUNT];
int i;
gsl_odeiv2_step_set_driver(s, d);
gsl_odeiv2_driver_set_nmax(d, 1);
initial_values(x);
initial_parameter_values(parameters);
initial_parameter_values(p);
print_headers(csv);
for (i = 0; i < EVENT_COUNT; i++) {
e[i] = 0;
}
while (t < TIME_FINAL) {
events(e, p);
int status = gsl_odeiv2_driver_apply(d, &t, TIME_FINAL, x);
if (status != GSL_SUCCESS && status != GSL_EMAXITER) {
fprintf(stderr, "error, return value=%d\n", status);
......
......@@ -198,6 +198,7 @@ list_model :-
),
list_model_initial_state,
list_model_parameters,
list_model_events,
list_model_functions,
list_model_options.
......
......@@ -8,6 +8,13 @@
patch_solution/2,
kinetics/1,
arithmetic_expression/1,
condition/1,
op(720, fy, not),
op(740, xfy, and),
op(760, xfy, or),
op(780, fy, if),
op(780, xfy, then),
op(780, xfy, else),
op(1100, xfx, for),
op(700, xfx, <=),
op(750, xfx, =>),
......@@ -120,6 +127,45 @@ kinetics((ArithmeticExpression0, ArithmeticExpression1)) :-
arithmetic_expression(ArithmeticExpression1).
:- grammar(condition).
condition(true).
condition(false).
condition(ArithmeticExpression0 = ArithmeticExpression1) :-
arithmetic_expression(ArithmeticExpression0),
arithmetic_expression(ArithmeticExpression1).
condition(ArithmeticExpression0 < ArithmeticExpression1) :-
arithmetic_expression(ArithmeticExpression0),
arithmetic_expression(ArithmeticExpression1).
condition(ArithmeticExpression0 <= ArithmeticExpression1) :-
arithmetic_expression(ArithmeticExpression0),
arithmetic_expression(ArithmeticExpression1).
condition(ArithmeticExpression0 > ArithmeticExpression1) :-
arithmetic_expression(ArithmeticExpression0),
arithmetic_expression(ArithmeticExpression1).
condition(ArithmeticExpression0 >= ArithmeticExpression1) :-
arithmetic_expression(ArithmeticExpression0),
arithmetic_expression(ArithmeticExpression1).
condition(Condition0 or Condition1) :-
condition(Condition0),
condition(Condition1).
condition(Condition0 and Condition1) :-
condition(Condition0),
condition(Condition1).
condition(not Condition) :-
condition(Condition).
:- grammar(arithmetic_expression).
......@@ -172,6 +218,11 @@ arithmetic_expression(product(Pattern in List, ArithmeticExpression)) :-
term(List),
arithmetic_expression(ArithmeticExpression).
arithmetic_expression(if Condition then IfTrue else IfFalse) :-
condition(Condition),
arithmetic_expression(IfTrue),
arithmetic_expression(IfFalse).
arithmetic_expression(Number) :-
number(Number).
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment