Commit 64e549da authored by Mathieu Hemery's avatar Mathieu Hemery
Browse files

Cleaning rosenbrock

parent 5d3f1e1d
......@@ -253,8 +253,8 @@ eval_row([Element|T],CurrentState,Parameters,[Evaluation|ET]) :-
eval_coeff(Expr, States, Parameters, Time, Evaluation) :-
rw_coeff(Expr, States, Parameters, Time, Rw_expr),!,
simplify(Rw_expr, Simple_expr), % may be skipped, to check with better tests on rail
Evaluation is Simple_expr.
% simplify(Rw_expr, Simple_expr), % skipped here for huge gains, trouble may arrive
Evaluation is Rw_expr.
rw_coeff(Element,_,_,_,Element) :-
......@@ -348,11 +348,9 @@ eval_state(State, Parameters, Evaluated) :-
% Compute New = Old + StepSize x Grad
% This a tool for the rosenbrock method
update_state([],_,[],[]).
update_state([Element|CurrentState],Constant,[Value|Values],[NewElement|NewCurrentState]) :-
NewElement is Element + Constant * Value,
update_state(CurrentState,Constant,Values,NewCurrentState).
update_state(Old, Stepsize, Grad, New) :-
multiply_list(Grad, Stepsize, Move),
add_list(Old, Move, New).
%! get_last_row(+Data, -Time, -CurrentState)
......@@ -657,51 +655,63 @@ scale_error([E|Err],[V|L],[SE|SErr]) :-
/* Events handling */
/*******************/
/* Conditionnal event is seen as true when it goes form False to True only */
:- dynamic(conditional_events_status/2).
:- dynamic(rosenbrock_events_status/2).
/* This gives acces to the traditionnal events list and the kinetics conditionnal events list */
%! get_events_lists(-EventsList, -ConditionalEventsList)
%
% Gives acces to the events list and the kinetics conditionnal events list
get_events_lists(EventsList,ConditionalEventsList) :-
nb_getval(events_list,Events),
length(Events,NumberOfEvents),
nb_getval(conditional_events,Number),
NumberOfConditionalEvents is Number*2,
Separation is NumberOfEvents - NumberOfConditionalEvents,
separate_list(Events,Separation,EventsList,ConditionalEventsList).
nb_getval(events_list,Events),
length(Events,NumberOfEvents),
nb_getval(conditional_events,Number),
NumberOfConditionalEvents is Number*2,
Separation is NumberOfEvents - NumberOfConditionalEvents,
separate_list(Events,Separation,EventsList,ConditionalEventsList).
/* Conditionnal event is seen as true when it goes form False to True only
This creates the initial facts concerning theses events */
%! prepare_conditional_events(+State, +Time)
%
% Creates the initial facts concerning theses events
prepare_conditional_events(State,Time) :-
nb_getval(conditional_events,Number),
NumberOfEvents is Number - 1,
forall(
between(0,NumberOfEvents,ParameterIndex),
(nb_getval(parameters_list,Parameters),
EventIndex is 2*ParameterIndex,
numerical_simulation:conditional_event(IndexedCondition, ParameterIndex),
test_condition(State,Time,IndexedCondition,Result),
assertz(conditional_events_status(EventIndex,Result)),
TrueParameterIndex is ParameterIndex + 1,
OppositeEventIndex is EventIndex + 1,
(
Result
->
assertz(conditional_events_status(OppositeEventIndex,false)),
replace(Parameters,TrueParameterIndex,1,NewParameters)
;
assertz(conditional_events_status(OppositeEventIndex,true)),
replace(Parameters,TrueParameterIndex,0,NewParameters)
),
nb_setval(parameters_list,NewParameters)
)
).
nb_getval(conditional_events,Number),
NumberOfEvents is Number - 1,
forall(
between(0,NumberOfEvents,ParameterIndex),
(
nb_getval(parameters_list,Parameters),
EventIndex is 2*ParameterIndex,
numerical_simulation:conditional_event(IndexedCondition, ParameterIndex),
test_condition(State,Time,IndexedCondition,Result),
assertz(conditional_events_status(EventIndex,Result)),
TrueParameterIndex is ParameterIndex + 1,
OppositeEventIndex is EventIndex + 1,
(
Result
->
assertz(conditional_events_status(OppositeEventIndex,false)),
replace(Parameters,TrueParameterIndex,1,NewParameters)
;
assertz(conditional_events_status(OppositeEventIndex,true)),
replace(Parameters,TrueParameterIndex,0,NewParameters)
),
nb_setval(parameters_list,NewParameters)
)
).
%! prepare_rosenbrock_events(+State, +Time)
%
%
prepare_rosenbrock_events(State,Time) :-
get_events_lists(Events,_),
forall(
get_events_lists(Events,_),
forall(
(member(Event, Events),Event =.. [->,Condition,AssignmentList]),
(
test_condition(State,Time,Condition,Result),
......@@ -719,51 +729,63 @@ prepare_rosenbrock_events(State,Time) :-
)
).
/* Check if there is some conditional events which become true during the time step of the rosenbrock stimulation
If there are indeed events, the simulation goes one step backward and switch to an explicit Euler method */
%! check_conditional_events(+State, +Time)
%
% Check if there is some conditional events which become true during the time step of the
% rosenbrock stimulation. If there are, the simulation goes one step backward and switch
% to an explicit Euler method
% /!\ I would prefer to go backward and make another call to rosenbrock with adjusted
% step size.
check_conditional_events(State,Time) :-
get_events_lists(_,Events),
forall(
(member(Event, Events),Event =.. [->,Condition|_]),
( test_condition(State,Time,Condition,Result),
nth0(EventIndex,Events,Event),
retract(conditional_events_status(EventIndex,EventInfo)),
assertz(conditional_events_status(EventIndex,Result)),
%format("Event last state is ~w event new state is ~w ~n",[EventInfo,Result]),
\+(EventInfo),
Result,
retract(conditional_events_status(EventIndex,Result)), % If the event goes from false to true, since we'll get back to the time step before, we get back to old values
assertz(conditional_events_status(EventIndex,EventInfo))
->
count(events_number,_)
;
true
)
).
/* Same for traditionnal events */
get_events_lists(_,Events),
forall(
(member(Event, Events),Event =.. [->,Condition|_]),
(
test_condition(State,Time,Condition,Result),
nth0(EventIndex,Events,Event),
retract(conditional_events_status(EventIndex,EventInfo)),
assertz(conditional_events_status(EventIndex,Result)),
%format("Event last state is ~w event new state is ~w ~n",[EventInfo,Result]),
\+(EventInfo),
Result,
retract(conditional_events_status(EventIndex,Result)), % If the event goes from false to true, since we'll get back to the time step before, we get back to old values
assertz(conditional_events_status(EventIndex,EventInfo))
->
count(events_number,_)
;
true
)
).
%! check_events(+State, +Time)
%
%
check_events(State,Time) :-
get_events_lists(Events,_),
set_counter(events_number,0),
forall(
(member(Event, Events),Event =.. [->,Condition|_]),
( test_condition(State,Time,Condition,Result),
nth0(EventIndex,Events,Event),
retract(rosenbrock_events_status(EventIndex,EventInfo)),
assertz(rosenbrock_events_status(EventIndex,Result)),
%format("Event last state is ~w event new state is ~w ~n",[EventInfo,Result]),
\+(EventInfo),
Result,
retract(rosenbrock_events_status(EventIndex,Result)), % If the event goes from false to true, since we'll get back to the time step before, we get back to old values
assertz(rosenbrock_events_status(EventIndex,EventInfo))
->
count(events_number,_)
;
true
)
).
get_events_lists(Events,_),
set_counter(events_number,0),
forall(
(member(Event, Events),Event =.. [->,Condition|_]),
(
test_condition(State,Time,Condition,Result),
nth0(EventIndex,Events,Event),
retract(rosenbrock_events_status(EventIndex,EventInfo)),
assertz(rosenbrock_events_status(EventIndex,Result)),
%format("Event last state is ~w event new state is ~w ~n",[EventInfo,Result]),
\+(EventInfo),
Result,
% If the event goes from false to true, since we'll get back to the time step
% before, we get back to old values
retract(rosenbrock_events_status(EventIndex,Result)),
assertz(rosenbrock_events_status(EventIndex,EventInfo))
->
count(events_number,_)
;
true
)
).
/* Used during the Euler method phase, and when an event is true, its changes are made */
......@@ -821,106 +843,100 @@ handle_rosenbrock_events(State,Time) :-
% Can also be used if conditional changes in events
test_condition(State,Time,Condition,Result) :-
Condition =.. [Op, A, B],
memberchk(Op, ['=', '<', '<=', '>', '>=']),
!,
test_condition(State,Time,A,ResultA),
test_condition(State,Time,B,ResultB),
(
Op = '<='
->
Op2 = '=<'
;
true
),
(
Op = '=>'
->
Op2 = '>='
;
true
),
(
var(Op2)
->
TrueOp = Op
;
TrueOp = Op2
),
VCondition =.. [TrueOp, ResultA, ResultB],
(
VCondition
->
Result = true
;
Result = false
).
Condition =.. [Op, A, B],
memberchk(Op, ['=', '<', '<=', '>', '>=']),
!,
test_condition(State,Time,A,ResultA),
test_condition(State,Time,B,ResultB),
(
Op = '<='
->
Op2 = '=<'
;
true
),
(
Op = '=>'
->
Op2 = '>='
;
true
),
(
var(Op2)
->
TrueOp = Op
;
TrueOp = Op2
),
VCondition =.. [TrueOp, ResultA, ResultB],
(
VCondition
->
Result = true
;
Result = false
).
test_condition(State,Time,not Condition,Result) :-
test_condition(State,Time,Condition,NotResult),
(
NotResult
->
Result = false
;
Result = true
).
test_condition(State,Time,Condition,NotResult),
(
NotResult
->
Result = false
;
Result = true
).
test_condition(State,Time,Condition,Result) :-
nb_getval(parameters_list,Parameters),
rw_coeff(Condition,State,Parameters,Time,Result).
nb_getval(parameters_list,Parameters),
rw_coeff(Condition,State,Parameters,Time,Result).
test_condition(_,_,true,true).
test_condition(_,_,false,false).
test_condition(State,Time,Condition1 and Condition2,Result) :-
test_condition(State,Time,Condition1,Result1),
test_condition(State,Time,Condition2,Result2),
(
Result1,
Result2
->
Result = true
;
Result = false
).
test_condition(State,Time,Condition1,Result1),
test_condition(State,Time,Condition2,Result2),
(
Result1,
Result2
->
Result = true
;
Result = false
).
test_condition(State,Time,Condition1 or Condition2,Result) :-
test_condition(State,Time,Condition1,Result1),
test_condition(State,Time,Condition2,Result2),
(
(Result1;Result2)
->
Result = true
;
Result = false
).
test_condition(State,Time,Condition1,Result1),
test_condition(State,Time,Condition2,Result2),
(
(Result1;Result2)
->
Result = true
;
Result = false
).
apply_changes(State,Time, Parameter = (if Condition then IfTrue else IfFalse)) :-
test_condition(State,_,Condition,Result),
(
Result
->
apply_changes(State,Time,Parameter = IfTrue)
;
apply_changes(State,Time,Parameter = IfFalse)
).
test_condition(State,_,Condition,Result),
(
Result
->
apply_changes(State,Time,Parameter = IfTrue)
;
apply_changes(State,Time,Parameter = IfFalse)
).
apply_changes(State,Time, Parameter = Expression) :-
IndexParameter is Parameter + 1,
nb_getval(parameters_list,Parameters),
eval_coeff(Expression,State,Parameters,Time,Value),
matrix:replace(Parameters,IndexParameter,Value,NewParameters),
%format("AncientParameters are ~w and NewParameters are ~w ~n",[Parameters,NewParameters]),
nb_setval(parameters_list,NewParameters).
scalar_multiply([],_,[]).
IndexParameter is Parameter + 1,
nb_getval(parameters_list,Parameters),
eval_coeff(Expression,State,Parameters,Time,Value),
matrix:replace(Parameters,IndexParameter,Value,NewParameters),
%format("AncientParameters are ~w and NewParameters are ~w ~n",[Parameters,NewParameters]),
nb_setval(parameters_list,NewParameters).
scalar_multiply([Element|Vector],Scalar,[NewElement|Result]) :-
NewElement is Element * Scalar,
scalar_multiply(Vector,Scalar,Result).
euler(Equations,Time,CurrentState,CurrentTable) :-
format("One event occured ~n",[]),
......@@ -952,7 +968,7 @@ euler(Equations,Time,CurrentState,CurrentTable) :-
euler_step(Equations,CurrentState,StepSize,Parameters,NewCurrentState) :-
eval_state(CurrentState,Parameters,EvaluatedCurrentState),
eval_row(Equations,CurrentState,Parameters,EvaluatedEquations),
scalar_multiply(EvaluatedEquations,StepSize,TrueEvaluatedEquations),
multiply_list(EvaluatedEquations,StepSize,TrueEvaluatedEquations),
row_sum(EvaluatedCurrentState,TrueEvaluatedEquations,NewCurrentState).
/*:- doc('\\begin{example}Ball simulation using event:').
......
Supports Markdown
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