Commit cfdda700 authored by Thierry Martinez's avatar Thierry Martinez

Fixed events

parent 3f81cc05
......@@ -79,10 +79,13 @@ simplify_blocks(Block, Index, RebuildCoef, Blocks, CoefSubBlocks) :-
gather_loop(Block, Blocks, SubBlocks),
maplist(Index, SubBlocks, IndexedSubBlocks),
sort(4, @=<, IndexedSubBlocks, SortedSubBlocks),
check_cleaned(arithmetic_rules:canonical/3),
gather_indexed(SortedSubBlocks),
reduce_blocks(IndexedSubBlocks, ReducedSubBlocks),
clean(arithmetic_rules:canonical/3),
with_clean(
[arithmetic_rules:canonical/3],
(
arithmetic_rules:gather_indexed(SortedSubBlocks),
arithmetic_rules:reduce_blocks(IndexedSubBlocks, ReducedSubBlocks)
)
),
map_blocks(RebuildCoef, ReducedSubBlocks, CoefSubBlocks).
......
......@@ -15,18 +15,30 @@
add_event(Condition, ParameterValues) :-
biocham_command(*),
type(Condition, condition),
type(ParameterValues, '*'(parameter = simple_kinetics)),
type(ParameterValues, '*'(parameter = arithmetic_expression)),
doc('
sets up an event that will be fired each time the condition given as first
argument goes from false to true.
This command is effective in numerical simulations only.
Upon firing, the parameters receive new values
computed from the expression.
The initial values of the parameters are restored after the simulation.'),
The initial values of the parameters are restored after the simulation.
\\begin{example}
'),
biocham_silent(clear_model),
biocham('MA'(k) for a => b),
biocham(set_parameter(k = 1)),
biocham(add_event(b > 0.5, k = 0)),
biocham(present(a)),
biocham(numerical_simulation(2, method: msbdf, maximum_step_size: 1e-3)),
biocham(plot),
doc('
\\end{example}
'),
\+ (
member(ParameterValue, ParameterValues),
\+ (
add_item(event, event(Condition, ParameterValue))
add_item([kind: event, item: event(Condition, ParameterValue)])
)
).
......
......@@ -145,7 +145,7 @@ write_jacobian_cell(I, J, Expression) :-
write_arithmetic_expression(Value) :-
number(Value),
!,
format('~d', [Value]).
format('~f', [Value]).
write_arithmetic_expression([VariableIndex]) :-
format('x[~d]', [VariableIndex]).
......@@ -313,7 +313,7 @@ write_events(Options) :-
),
write('
void
events(int e[], double p[])
events(int e[], double p[], double x[])
{
'),
\+ (
......@@ -321,18 +321,27 @@ events(int e[], double p[])
\+ (
format(' if (!e[~d] && (', [EventIndex]),
write_condition(Condition),
format(')) {
format(
')) {
e[~d] = 1;
p[~d] = ', [EventIndex, ParameterIndex]),
p[~d] = ',
[EventIndex, ParameterIndex]
),
write_arithmetic_expression(Value),
format(';
format(
';
}
else if (e[~d] && !(', [EventIndex]),
else if (e[~d] && !(',
[EventIndex]
),
write_condition(Condition),
format(')) {
format(
')) {
e[~d] = 0;
}
')
',
[EventIndex]
)
)
),
write('\c
......
......@@ -42,7 +42,7 @@ main(void)
if (t_upper >= TIME_FINAL) {
t_upper = TIME_FINAL;
}
events(e, p);
events(e, p, x);
status = gsl_odeiv2_driver_apply(d, &t, t_upper, x);
if (status != GSL_SUCCESS && status != GSL_EMAXITER) {
fprintf(stderr, "error, return value=%d\n", status);
......
......@@ -82,35 +82,51 @@ numerical_simulation(Time) :-
),
doc('performs a numerical simulation up to a given time.'),
with_current_ode_system((
numerical_simulation:make_ode_system,
numerical_simulation:solve(
Time, Method, ErrorEpsilonAbsolute, ErrorEpsilonRelative, InitialStepSize,
MaximumStepSize, Precision
with_clean(
[
numerical_simulation:variable/2,
numerical_simulation:equation/2,
numerical_simulation:parameter/2
],
(
numerical_simulation:solve(
Time, Method, ErrorEpsilonAbsolute, ErrorEpsilonRelative,
InitialStepSize, MaximumStepSize, Precision
)
)
)
)).
make_ode_system :-
enumerate_variables,
convert_ode.
:- devdoc('\\section{Private predicate}').
:- devdoc('\\section{Private predicate}').
:- dynamic(variable/2).
:- dynamic(equation/2).
:- dynamic(parameter/2).
solve(
Time, Method, ErrorEpsilonAbsolute, ErrorEpsilonRelative, InitialStepSize,
MaximumStepSize, Precision
) :-
enumerate_variables,
convert_ode,
gather_headers(Headers),
gather_equations(Equations),
gather_initial_values(InitialValues),
gather_initial_parameter_values(InitialParameterValues),
gather_events(Events),
Options = [
headers: Headers,
equations: Equations,
initial_values: InitialValues,
initial_parameter_values: InitialParameterValues,
events: Events,
method: Method,
error_epsilon_absolute: ErrorEpsilonAbsolute,
error_epsilon_relative: ErrorEpsilonRelative,
......@@ -180,17 +196,20 @@ gather_initial_parameter_values(InitialParameterValues) :-
gather_events(Events) :-
% findall(
% Event,
% (
% item([kind: event; item: event(Condition, Parameter = Value)]),
Events = [].
findall(
(ConditionIndexed -> ParameterIndex = ValueIndexed),
(
item([kind: event, item: event(Condition, Parameter = Value)]),
convert_condition(Condition, ConditionIndexed),
parameter(Parameter, ParameterIndex),
convert_expression(Value, ValueIndexed)
),
Events
).
:- dynamic(variable/2).
enumerate_variables :-
retractall(variable(_, _)),
set_counter(variable_counter, 0),
\+ (
ode(Molecule, _),
......@@ -202,16 +221,7 @@ enumerate_variables :-
).
:- dynamic(equation/2).
:- dynamic(parameter/2).
convert_ode :-
print(convert_ode),
nl,
retractall(equation(_, _)),
set_counter(parameter_counter, 0),
\+ (
ode(Molecule, Expression),
......@@ -220,9 +230,7 @@ convert_ode :-
convert_expression(Expression, IndexedExpression),
assertz(equation(VariableIndex, IndexedExpression))
)
),
print(end_convert_ode),
nl.
).
convert_condition(Expression, IndexedExpression) :-
......
......@@ -8,5 +8,12 @@ test('mapk') :-
command(load(library:examples/mapk/mapk)),
numerical_simulation(100).
test('events') :-
clear_model,
command('MA'(k) for a => b),
command(set_parameter(k = 1)),
command(add_event(b > 0.5, k = 0)),
command(present(a)),
numerical_simulation(20).
:- end_tests(numerical_simulation).
......@@ -173,11 +173,14 @@ ode_system :-
delete_ode_system('ode_system'),
new_ode_system,
set_ode_system_name('ode_system'),
check_cleaned(ode:assoc/2),
compute_ode,
simplify_ode,
put_ode_into_system,
clean(ode:assoc/2).
with_clean(
[ode:assoc/2],
(
ode:compute_ode,
ode:simplify_ode,
ode:put_ode_into_system
)
).
import_ode(InputFile) :-
......
......@@ -29,7 +29,8 @@ set_parameter(ParameterList) :-
get_parameter_value(Parameter, Value) :-
find_item([kind: parameter, key: Parameter, item: parameter(_, Value)]).
find_item([kind: parameter, key: Parameter, item: parameter(_ = Value0)]),
Value = Value0.
show_parameter(Parameter) :-
......
......@@ -228,15 +228,15 @@ arithmetic_expression(if Condition then IfTrue else IfFalse) :-
arithmetic_expression(IfTrue),
arithmetic_expression(IfFalse).
arithmetic_expression(FunctionApplication) :-
function_application(arithmetic_expression, FunctionApplication).
arithmetic_expression(Number) :-
number(Number).
arithmetic_expression(Name) :-
name(Name).
arithmetic_expression(FunctionApplication) :-
function_application(arithmetic_expression, FunctionApplication).
:- devdoc('\\section{Public API}').
......
......@@ -97,6 +97,7 @@ command(Command) :-
append(Prefix, [Tail], NewArguments),
Command0 =.. [Functor | NewArguments]
;
Arity >= 1,
predicate_info(Functor/1, ArgumentTypes, CommandOptions, yes, _)
->
CommandWithoutOptions =.. [Functor | Arguments],
......
......@@ -51,7 +51,9 @@ predicate_info(Functor/Arity, ArgumentTypes, Options, BiochamCommand, Doc) :-
catch(
(
clause(Head, Body),
predicate_info((Head :- Body), ArgumentTypes, Options, BiochamCommand, Doc),
predicate_info(
(Head :- Body), ArgumentTypes, Options, BiochamCommand, Doc
),
(
BiochamCommand = yes
;
......@@ -91,7 +93,9 @@ collect_info(
) :-
!.
collect_info(DocItem, _Arguments, _ArgumentTypes, [], _BiochamCommand, [DocItem]) :-
collect_info(
DocItem, _Arguments, _ArgumentTypes, [], _BiochamCommand, [DocItem]
) :-
(
DocItem = doc(_)
;
......@@ -178,13 +182,13 @@ check_type(Grammar, Item, NewItem) :-
copy_term((NewHead, NewBody), (Head, Body)),
catch(
check_grammar_body(Body, NewBody),
error(expected(_)),
error(expected(_, _)),
fail
)
->
true
;
throw(error(expected(Grammar)))
throw(error(expected(Grammar, Item)))
).
......@@ -234,4 +238,5 @@ check_grammar_body(Goal, NewGoal) :-
function_application(Kind, Term) :-
Term =.. [_Function | Arguments],
Arguments = [_ | _],
maplist(Kind, Arguments).
......@@ -23,6 +23,7 @@
rewrite/3,
grammar_map/4,
call_subprocess/3,
with_clean/2,
clean/1,
check_cleaned/1
]).
......@@ -297,6 +298,14 @@ user:message_hook(_Msg, warning, _Lines) :-
fail.
with_clean(Prototypes, Goal) :-
setup_call_cleanup(
maplist(check_cleaned, Prototypes),
Goal,
maplist(clean, Prototypes)
).
clean(F / N) :-
functor(H, F, N),
retractall(H).
......
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