Une MAJ de sécurité est nécessaire sur notre version actuelle. Elle sera effectuée lundi 02/08 entre 12h30 et 13h. L'interruption de service devrait durer quelques minutes (probablement moins de 5 minutes).

rosenbrock.pl 31.7 KB
Newer Older
1 2 3 4 5
:- module(rosenbrock,
   [
      ode_solver/1
   ] 
).
Paul Remondeau's avatar
Test  
Paul Remondeau committed
6

7
% Only for separate compilation and linting
HEMERY Mathieu's avatar
HEMERY Mathieu committed
8 9
:- use_module(arithmetic_rules).
:- use_module(counters).
10
:- use_module(doc).
11
:- use_module(util).
12

Mathieu Hemery's avatar
Mathieu Hemery committed
13 14
%! rosenbrock_init
%
15 16
% Initialize coefficients for the rosenbrock method, the errors and the minimal step size
% for events handling
17
 
18 19
rosenbrock_init(Epsilon_abs) :-
  nb_setval(not_zero_error, Epsilon_abs),
Mathieu Hemery's avatar
Mathieu Hemery committed
20
  nb_setval(k_method,stiff),
21 22
  nb_setval(event_step_size, 1e-3),
  nb_setval(step_doubling_error, Epsilon_abs), % Control the method's precision
23
  SAFE = 0.9, GROW = 2.5, PGROW = -0.25,
Mathieu Hemery's avatar
Mathieu Hemery committed
24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
  SHRINK = 0.5, PSHRINK is -1.0/3.0, ERRCON = 0.1296,
  GAM is 1.0/2.0, A21=2.0, A31 is 48.0/25.0, A32 is 6.0/25.0, C21 = -8.0,
  C31 is 372.0/25.0, C32 is 12.0/5.0, C41 is -112.0/125.0,
  C42 is -54.0/125.0, C43 is -2.0/5.0,
  B1 is 19.0/9.0, B2 is 1.0/2.0, B3 is 25.0/108.0, B4 is 125.0/108.0,
  E1 is 17.0/54.0, E2 is 7.0/36.0, E3 is 0.0, E4 is 125.0/108.0,
  nb_setval(k_rbk_safe,SAFE),
  nb_setval(k_rbk_grow,GROW),
  nb_setval(k_rbk_pgrow,PGROW),
  nb_setval(k_rbk_shrink,SHRINK),
  nb_setval(k_rbk_pshrink,PSHRINK),
  nb_setval(k_rbk_errcon,ERRCON),
  nb_setval(k_rbk_gam,GAM),
  nb_setval(k_rbk_a21,A21),
  nb_setval(k_rbk_a31,A31),
  nb_setval(k_rbk_a32,A32),
  nb_setval(k_rbk_c21,C21),
  nb_setval(k_rbk_c31,C31),
  nb_setval(k_rbk_c32,C32),
  nb_setval(k_rbk_c41,C41),
  nb_setval(k_rbk_c42,C42),
  nb_setval(k_rbk_c43,C43),
  nb_setval(k_rbk_b1,B1),
  nb_setval(k_rbk_b2,B2),
  nb_setval(k_rbk_b3,B3),
  nb_setval(k_rbk_b4,B4),
  nb_setval(k_rbk_e1,E1),
  nb_setval(k_rbk_e2,E2),
  nb_setval(k_rbk_e3,E3),
  nb_setval(k_rbk_e4,E4).
Paul Remondeau's avatar
Test  
Paul Remondeau committed
54

55
:- dynamic(k_fail/0).
56

57 58
:- dynamic(rsbk_continue/0).

Mathieu Hemery's avatar
Mathieu Hemery committed
59 60 61 62

%! ode_solver(+Options)
%
% Start of the ode solver, cf numerical_simulation.pl for Options
63 64

ode_solver([
65
    fields: Fields,
66 67
    equations: Equations,
    initial_values: InitialState,
68
    initial_parameter_values: InitialParameters,
69
    events: Events,
70
    method: _,
71
    error_epsilon_absolute: Epsilon_abs,
72 73
    error_epsilon_relative: _,
    initial_step_size: _,
74 75
    maximum_step_size: MaxStSz,
    minimum_step_size: MinStSz,
76 77
    precision: _,
    filter: _,
78
    time_initial: InitialTime,
79
    time_final: Duration,
80 81
    jacobian: Jacobian
  ]) :- 
Mathieu Hemery's avatar
Mathieu Hemery committed
82
  (
83
    nb_current(rsbk_continue,yes)
Mathieu Hemery's avatar
Mathieu Hemery committed
84
  ->
85
    get_last_row(InitialTime2,TrueInitialState),
Mathieu Hemery's avatar
Mathieu Hemery committed
86 87 88 89 90 91 92 93 94 95 96 97
    prepare_rosenbrock_events(TrueInitialState,InitialTime2),
    prepare_conditional_events(TrueInitialState,InitialTime2),
    TrueInitialTime is InitialTime2
  ;
    nb_setval(parameters_list,InitialParameters),
    nb_setval(variable_list,[]),
    length(Equations,NumberOfVariable),
    NNumberOfVariable is NumberOfVariable - 1,
    forall(
             between(0,NNumberOfVariable,VariableIndex),
             (
               nb_getval(variable_list,VariableList),
98 99 100
               numerical_simulation:convert_identifier(Molecule, [VariableIndex]),
               append(VariableList,[Molecule],NewVariableList),
               nb_setval(variable_list,NewVariableList)
Mathieu Hemery's avatar
Mathieu Hemery committed
101 102 103
             )
          ),
    nb_getval(variable_list,VariableList),
104
    initialize_functions_list,
105 106 107
    % add_functions_names(VariableList, Var_Func_List),
    initialize_names(Fields, Names),
    assertz(saved_row(Names)), %initialize the label of the numerical table
Mathieu Hemery's avatar
Mathieu Hemery committed
108
    eval_state(InitialState,InitialParameters,TrueInitialState),
109 110 111
    split_events(Events, RegularEvents, TimeEvents),
    nb_setval(events_list,RegularEvents),
    nb_setval(time_events_list,TimeEvents),
112
    nb_setval(last_time_event, -1),
Mathieu Hemery's avatar
Mathieu Hemery committed
113 114
    prepare_rosenbrock_events(TrueInitialState,InitialTime),
    prepare_conditional_events(TrueInitialState,InitialTime),
115 116
    add_functions_values(TrueInitialState, InitialTime, FullState),
    SecondRow =.. [row,InitialTime|FullState],
117
    nb_setval(last_row, SecondRow),
118
    log_current_row(InitialTime, FullState, InitialParameters, Fields),
Mathieu Hemery's avatar
Mathieu Hemery committed
119 120
    TrueInitialTime is InitialTime
  ),
121
  rosenbrock_init(Epsilon_abs),
122
  rosenbrock(Equations,TrueInitialTime,Duration,MaxStSz,MinStSz,Jacobian,Fields),
Mathieu Hemery's avatar
Mathieu Hemery committed
123 124 125 126 127 128
  retractall(k_fail).


%! initialize_functions_list
%
% Create a global variable: a list which contains all the functions defined in the .bc file
129

130
initialize_functions_list :-
Mathieu Hemery's avatar
Mathieu Hemery committed
131 132 133 134 135 136 137 138 139
  all_items([kind: function],Functions),
  initialize_functions_list(Functions,FunctionsList),
  (
    FunctionsList = []
  ->
    nb_setval(functions_number,0)
  ;
    nb_setval(functions_list,FunctionsList),
    length(FunctionsList,Number),
140
    nb_setval(functions_number,Number)
Mathieu Hemery's avatar
Mathieu Hemery committed
141
  ).
142

Mathieu Hemery's avatar
Mathieu Hemery committed
143 144 145 146 147 148 149 150 151 152 153 154 155 156 157
initialize_functions_list([],[]) :- !.

initialize_functions_list([Function|Functions], FunctionList) :- 
  Function =.. [function,TrueFunction],
  TrueFunction =.. [=,Head,Expression],
  initialize_functions_list(Functions,NextList),
  ( % skip over kinetic functions
    is_kinetics_function(Head)
  ->
    FunctionList = NextList
  ;
    numerical_simulation:convert_expression(Expression,IndexedExpression),
    NewFunction = function(Head,IndexedExpression),
    append(NextList,[NewFunction],FunctionList)
  ).
158 159


Mathieu Hemery's avatar
Mathieu Hemery committed
160 161 162 163 164 165 166 167 168 169 170
%! is_kinetics_function(F)
%
% Manually check if F is a kinetic function
% Rq: It'll have to be updated if others kinetics are added

is_kinetics_function(Head) :-
  Head = 'MA'(k);
  Head = 'MAI'(k);
  Head = 'MM'('Vm', 'Km');
  Head = 'Hill'('Vm', 'Km', n);
  Head = 'HillI'('Vm', 'Km', n).
171

172

173
%! add_functions_names(+VarName, -VarFuncName)
Mathieu Hemery's avatar
Mathieu Hemery committed
174
%
175
% Add the name of all functions to construct the FullNames list
176

177 178 179 180 181 182 183 184 185 186
add_functions_names(VarName, VarFuncName) :-
  (
    nb_current(functions_list,_)
  ->
    nb_getval(functions_list, Functions),
    maplist(give_function_name, Functions, FuncName),
    append(VarName, FuncName, VarFuncName)
  ;
    VarFuncName = VarName
  ).
Mathieu Hemery's avatar
Mathieu Hemery committed
187

188
give_function_name(function(Name, _Expr), Name).
Mathieu Hemery's avatar
Mathieu Hemery committed
189

190

191 192 193 194 195 196 197 198 199 200
%! initialize_names(+Fields, -Names)
%
% Extract the name of the different fields to be plotted

initialize_names_int([], []).

initialize_names_int([Name:_|FTail], [Name|NTail]) :-
  initialize_names_int(FTail, NTail).

initialize_names([Name:_|FTail], [NName|NTail]) :-
201
  atom_concat('#', Name, NName),
202 203 204
  initialize_names_int(FTail, NTail).


205
%! log_current_row(+Time, +State, +Parameters, +Fields)
206

207 208
log_current_row(Time, State, Parameters, Fields) :-
  prepare_row(Time, State, Parameters, Fields, ToSave),
209 210
  assertz(saved_row(ToSave)).

211
prepare_row(_Time, _State, _Param, [], []) :- !.
212

213
prepare_row(Time, State, Parameters, [_N:t|FTail], [Time|VTail]) :-
214
  !,
215
  prepare_row(Time, State, Parameters, FTail, VTail).
216

217 218
prepare_row(Time, State, Parameters, [_N:x(N)|FTail], [Val|VTail]) :-
  !,
219
  nth0(N, State, Val),
220 221 222 223 224 225
  prepare_row(Time, State, Parameters, FTail, VTail).

prepare_row(Time, State, Parameters, [_N:expression(Expr)|FTail], [Val|VTail]) :-
  !,
  eval_coeff(Expr, State, Parameters, Time, Val),
  prepare_row(Time, State, Parameters, FTail, VTail).
226 227


228
%! add_functions_values(+State, +Time, -FullState)
Mathieu Hemery's avatar
Mathieu Hemery committed
229
%
230
% Add the value of all functions to construct the FullState list
231

232 233 234 235 236 237 238 239 240 241 242 243 244 245
add_functions_values(State, Time, FullState) :-
  (
    nb_current(functions_list,_)
  ->
    nb_getval(parameters_list, Parameters),
    nb_getval(functions_list, Functions),
    maplist(eval_function(State, Parameters, Time), Functions, Values),
    append(State, Values, FullState)
  ;
    FullState = State
  ).

eval_function(State, Parameters, Time, function(_Name, Expr), Value) :-
  eval_coeff(Expr, State, Parameters, Time, Value).
246

247
%! eval_matrix(+Matrix, +States, +Parameters, -Evaluated_Matrix)
248
%
249 250
% insert all the numerical value in a symbolic matrix
% States and Parameters are simple list of value
251

252 253
eval_matrix([],_,_,[]). 
eval_matrix([Row|T],CurrentState,Parameters,[RowEvaluated|ET]) :-
254 255
   eval_row(Row,CurrentState,Parameters,RowEvaluated),
   eval_matrix(T,CurrentState,Parameters,ET).
256

257
eval_row([],_,_,[]).
258 259
eval_row([Element|T],CurrentState,Parameters,[Evaluation|ET]) :-
  eval_coeff(Element, CurrentState, Parameters, _, Evaluation),
260 261 262
  eval_row(T,CurrentState,Parameters,ET).


263 264 265 266 267 268 269 270
%! eval_coeff(+Expr, +States, +Parameters, +Time, -Evaluation)
%
% Evaluate an expression by rewritting the occurence of the concentration and parameters,
% simplyfying the resulting expression and calling is.
% The rewritting is delegated to rw_coeff

eval_coeff(Expr, States, Parameters, Time, Evaluation) :-
  rw_coeff(Expr, States, Parameters, Time, Rw_expr),!,
Mathieu Hemery's avatar
Mathieu Hemery committed
271 272
  % simplify(Rw_expr, Simple_expr), % skipped here for huge gains, trouble may arrive
  Evaluation is Rw_expr.
273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292


rw_coeff(Element,_,_,_,Element) :-
  number(Element),!.

rw_coeff(Element1 ^ Element2,CurrentState,Parameters,Time,ElementEvaluated1 ^ ElementEvaluated2) :-
  rw_coeff(Element1,CurrentState,Parameters,Time,ElementEvaluated1),
  rw_coeff(Element2,CurrentState,Parameters,Time,ElementEvaluated2).

rw_coeff(Element1 + Element2,CurrentState,Parameters,Time,ElementEvaluated1 + ElementEvaluated2) :-
	rw_coeff(Element1,CurrentState,Parameters,Time,ElementEvaluated1),
	rw_coeff(Element2,CurrentState,Parameters,Time,ElementEvaluated2).

rw_coeff(Element1 - Element2,CurrentState,Parameters,Time,ElementEvaluated1 - ElementEvaluated2) :-
	rw_coeff(Element1,CurrentState,Parameters,Time,ElementEvaluated1),
	rw_coeff(Element2,CurrentState,Parameters,Time,ElementEvaluated2).	

rw_coeff(Element1 * Element2,CurrentState,Parameters,Time,ElementEvaluated1 * ElementEvaluated2) :-
	rw_coeff(Element1,CurrentState,Parameters,Time,ElementEvaluated1),
	rw_coeff(Element2,CurrentState,Parameters,Time,ElementEvaluated2).
293

294 295 296
rw_coeff(Element1 / Element2,CurrentState,Parameters,Time,ElementEvaluated1 / ElementEvaluated2) :-
   rw_coeff(Element1,CurrentState,Parameters,Time,ElementEvaluated1),
   rw_coeff(Element2,CurrentState,Parameters,Time,ElementEvaluated2).
297

298 299
rw_coeff(-Element,CurrentState,Parameters,Time,-ElementEvaluated) :-
	rw_coeff(Element,CurrentState,Parameters,Time,ElementEvaluated).
300

301 302
rw_coeff([IndexVariable],CurrentState,_,_,ElementEvaluated) :-
	nth0(IndexVariable,CurrentState,ElementEvaluated).
303

304 305 306
rw_coeff(p(ParameterIndex),CurrentState,Parameters,Time,ElementEvaluated) :-
   nth0(ParameterIndex,Parameters,Element2),
   rw_coeff(Element2,CurrentState,Parameters,Time,ElementEvaluated).
307

308 309 310 311 312 313 314 315
rw_coeff(floor(Element),CurrentState,Parameters,Time,Result) :-
   !,
   rw_coeff(Element,CurrentState,Parameters,Time,ElementEval),
   Result is floor(ElementEval).

rw_coeff(random,_,_,_,Result) :-
   !,Result is random_float.

SOLIMAN Sylvain's avatar
merge  
SOLIMAN Sylvain committed
316 317 318 319 320 321 322 323 324
rw_coeff(infinity,_,_,_,Result) :-
  !,
  catch(
    Result is inf,
    % old SWIPL (7.2.3) doesn't have inf/0
    error(type_error(evaluable, inf/0), _),
    % use something big
    current_flag(max_tagged_integer, Result)
  ).
325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365

rw_coeff(t,_,_,Time,Time) :-
   !.

rw_coeff(step_size,_,_,_,Result) :-
   !,
   nb_getval(hdid,Result).

rw_coeff(min(Element1, Element2),CurrentState,Parameters,Time,min(ElementEvaluated1,ElementEvaluated2)) :-
   !,
   rw_coeff(Element1,CurrentState,Parameters,Time,ElementEvaluated1),
   rw_coeff(Element2,CurrentState,Parameters,Time,ElementEvaluated2).

rw_coeff(max(Element1, Element2),CurrentState,Parameters,Time,max(ElementEvaluated1,ElementEvaluated2)) :-
   !,
   rw_coeff(Element1,CurrentState,Parameters,Time,ElementEvaluated1),
   rw_coeff(Element2,CurrentState,Parameters,Time,ElementEvaluated2).

rw_coeff(log(Element),CurrentState,Parameters,Time,log(ElementEvaluated)) :-
   !,
   rw_coeff(Element,CurrentState,Parameters,Time,ElementEvaluated).

rw_coeff(exp(Element),CurrentState,Parameters,Time,exp(ElementEvaluated)) :-
   !,
   rw_coeff(Element,CurrentState,Parameters,Time,ElementEvaluated).

rw_coeff(cos(Element),CurrentState,Parameters,Time,cos(ElementEvaluated)) :-
   !,
   rw_coeff(Element,CurrentState,Parameters,Time,ElementEvaluated).

rw_coeff(sin(Element),CurrentState,Parameters,Time,sin(ElementEvaluated)) :-
   !,
   rw_coeff(Element,CurrentState,Parameters,Time,ElementEvaluated).


%! eval_state(+State, +Parameters, -Evaluated_state)
%
% evaluate the parameter in a state

eval_state(State, Parameters, Evaluated) :-
  eval_row(State, [], Parameters, Evaluated).
366 367


HEMERY Mathieu's avatar
HEMERY Mathieu committed
368 369 370 371
%! update_state(+Old, +StepSize, +Grad, -New) 
%
% Compute New = Old + StepSize x Grad
% This a tool for the rosenbrock method
372

Mathieu Hemery's avatar
Mathieu Hemery committed
373 374 375
update_state(Old, Stepsize, Grad, New) :-
  multiply_list(Grad, Stepsize, Move),
  add_list(Old, Move, New).
376 377


HEMERY Mathieu's avatar
HEMERY Mathieu committed
378 379 380 381
%! get_last_row(+Data, -Time, -CurrentState)
%
% 

382 383
get_last_row(Time,CurrentState):-
  nb_getval(last_row, LastRow),
HEMERY Mathieu's avatar
HEMERY Mathieu committed
384
  LastRow=..[row,Time|CurrentState].
385

386

387
/*********************/
388
/* Rosenbrock method */
389
/*********************/
HEMERY Mathieu's avatar
HEMERY Mathieu committed
390

391
%! rosenbrock(+Equations, +InitialTime, +Duration, +Max, +Min, +Jacobian, +Fields)
HEMERY Mathieu's avatar
HEMERY Mathieu committed
392 393 394 395
%
% Parameters are stored within the global variable parameters_list
% The current state of the simulation if obtained in the final row of the global
% variable numerical_table
396

397
rosenbrock(Equations, InitialTime, Duration, MaxStSz, MinStSz, Jacobian, Fields) :-
HEMERY Mathieu's avatar
HEMERY Mathieu committed
398
  FinalTime is InitialTime + Duration,
399 400 401
  StepSizeMax is Duration*MaxStSz,
  StepSizeMin is Duration*MinStSz,
  nb_setval(event_step_size, StepSizeMin),
HEMERY Mathieu's avatar
HEMERY Mathieu committed
402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417
  (
    StepSizeMax =:= 0.0,
    !,
    debug(error, "Simulation length cannot be 0", []),
    fail
  ;
    true
  ),
  nb_setval(hmax,StepSizeMax),
  (
    nb_current(h, _) % same as nb_getval but fail if 'h' do not exist instead of crashing
  ->
    true
  ;
    StepSize is StepSizeMax/5,
    nb_setval(h,StepSize),
418
    nb_setval(hdid,StepSize),
419
    nb_setval(h_init,StepSize),
HEMERY Mathieu's avatar
HEMERY Mathieu committed
420 421
    nb_setval(k_time_units,2)
  ),
422
  nb_setval(event_handling, false),% flag to know if we are currently dealing with events 
HEMERY Mathieu's avatar
HEMERY Mathieu committed
423
  repeat,
424
    get_last_row(Time, State_full),
HEMERY Mathieu's avatar
HEMERY Mathieu committed
425 426
    debug(rosenbrock, "entering step for ~p", [Time]),
    length(Equations, NumberTrueVariable), % Cut the function values at the end of State_full 
427
    separate_list(State_full, NumberTrueVariable, CurrentState, _FunctionState),
428
    handle_time_event(CurrentState, Time, NextEventTime),
429 430 431 432 433 434 435
    % As time events may change parameters, we need to fetch them now !
    nb_getval(parameters_list, Parameters),
    % Tailor cutting for the integration time step
    nb_getval(event_step_size, SmallestStep),
    EventStep is max((NextEventTime - Time), SmallestStep),
    FinalStep is (FinalTime - Time),
    MaxStep is min(EventStep, FinalStep),
HEMERY Mathieu's avatar
HEMERY Mathieu committed
436
    % Here is the true rosenbrock step
437
    rosenbrock_step(Equations, Jacobian, Time, CurrentState, Parameters, MaxStep, Next_state),
HEMERY Mathieu's avatar
HEMERY Mathieu committed
438
    nb_getval(hdid,NewStepSize), 
439
    Time2 is Time+NewStepSize,
440 441 442
    add_functions_values(Next_state, Time2, Next_state_full),
    NewLine =.. [row,Time2|Next_state_full],
    nb_setval(last_row, NewLine),
443
    log_current_row(Time2, Next_state_full, Parameters, Fields),
444
    % assertz(saved_row([Time2|Next_state_full])),
445 446
    nb_setval(last_time,Time2),
  Time2 >= FinalTime,!,
447 448 449
  %Collect the data for numerical table
  findall(Row, (saved_row(Data), Row =.. [row|Data]), RowList),
  retractall(saved_row(Data)),
450 451 452 453 454 455 456 457 458
  ( % patch to append the simulation when continue
    nb_current(rsbk_continue,yes)
  ->
    nb_getval(numerical_table, PreviousList),
    append(PreviousList, RowList, Complete_list)
  ;
    Complete_list = RowList
  ),
  nb_setval(numerical_table, Complete_list).
HEMERY Mathieu's avatar
HEMERY Mathieu committed
459 460


461
%! rosenbrock_step(+Equations, +Jacobian, +Time, +CurrentState, +Parameters, +MaxStep, -NextState)
HEMERY Mathieu's avatar
HEMERY Mathieu committed
462 463 464
%
% perform one step of rosenbrock integration, take also care of adjusting the step_size
% (stored in global variable nb_getval(hdid, _))
465
% MaxStep avoid overtaking the final step
466
% The next stepsize being stored in the 'h' global variable
467

468
rosenbrock_step(Equations,Jacobian,Time,CurrentState,Parameters,MaxStep,NewCurrentState) :-
HEMERY Mathieu's avatar
HEMERY Mathieu committed
469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495
  nb_getval(k_rbk_safe,SAFE),
  nb_getval(k_rbk_grow,GROW),
  nb_getval(k_rbk_pgrow,PGROW),
  nb_getval(k_rbk_shrink,SHRINK),
  nb_getval(k_rbk_pshrink,PSHRINK),
  nb_getval(k_rbk_errcon,ERRCON),
  nb_getval(k_rbk_gam,GAM),
  nb_getval(k_rbk_a21,A21),
  nb_getval(k_rbk_a31,A31),
  nb_getval(k_rbk_a32,A32),
  nb_getval(k_rbk_c21,C21),
  nb_getval(k_rbk_c31,C31),
  nb_getval(k_rbk_c32,C32),
  nb_getval(k_rbk_c41,C41),
  nb_getval(k_rbk_c42,C42),
  nb_getval(k_rbk_c43,C43),
  nb_getval(k_rbk_b1,B1),
  nb_getval(k_rbk_b2,B2),
  nb_getval(k_rbk_b3,B3),
  nb_getval(k_rbk_b4,B4),
  nb_getval(k_rbk_e1,E1),
  nb_getval(k_rbk_e2,E2),
  nb_getval(k_rbk_e3,E3),
  nb_getval(k_rbk_e4,E4),
  repeat,
    % evaluate Jacobian at current time point
    eval_state(CurrentState,Parameters,TrueCurrentState),
496
    eval_matrix(Jacobian,TrueCurrentState,Parameters,JacobianCurrentState),
HEMERY Mathieu's avatar
HEMERY Mathieu committed
497 498 499 500
    debug(rosenbrock,"Jacobian evaluated",[]),
    debug(rosenbrock,"Jacobian = ~w     CurrentState = ~w",[Jacobian, TrueCurrentState]),
    debug(rosenbrock,"JacobianCurrentState = ~w",[JacobianCurrentState]),

501 502
    nb_getval(h,StepSizeRaw),
    StepSize is min(StepSizeRaw, MaxStep),
HEMERY Mathieu's avatar
HEMERY Mathieu committed
503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591
    nb_setval(hdid,StepSize),
    length(Equations,Length),

    % fill the left hand side of AX=B
    GGAM is 1.0/(GAM*StepSize),
    create_diag_matrix(Length,GGAM,_),
    sub_matrix(JacobianCurrentState),

    % LU decompose A
    ludcmp_matrix,
    debug(rosenbrock,"Left side evaluated",[]), 

    % fill the right hand side, G1
    eval_row(Equations,TrueCurrentState,Parameters,G1),
    debug(rosenbrock,"Right side evaluated",[]),
    debug(rosenbrock,"Equations = ~w     CurrentState = ~w",[Equations, TrueCurrentState]),
    debug(rosenbrock,"G1 = ~w",[G1]),

    % compute by backsubstitution
    lusubst(G1,G1X),
    debug(rosenbrock,"G1=~p",[G1X]),
    % implement inc_vector/4 and inc_val_list/4 (which will be the same predicate I think)

    % get the new values at time + A2X * StepSize
    update_state(TrueCurrentState,A21,G1X,CurrentState1),
    debug(rosenbrock,"CurrentState1 = ~w",[CurrentState1]),
    eval_row(Equations,CurrentState1,Parameters,Equations1),
    debug(rosenbrock,"Equations at +A2X = ~p",[Equations1]),
    % fill rhs again for G2
    C21H is C21/StepSize,
    update_state(Equations1,C21H,G1X,G2),
    debug(rosenbrock,"Right side evaluated",[]),

    % compute by backsubstitution
    lusubst(G2,G2X),
    debug(rosenbrock,"G2=~p",[G2X]),

    % get the new values at time + (A2X+A3X) * StepSize
    update_state(TrueCurrentState,A31,G1X,CurrentState21),
    update_state(CurrentState21,A32,G2X,CurrentState2),
    eval_row(Equations,CurrentState2,Parameters,Equations2),
    debug(rosenbrock,"Equations at +A2X+A3X = ~p",[Equations2]),
    % fill rhs again for G3
    C31H is C31/StepSize,
    C32H is C32/StepSize,
    update_state(Equations2,C31H,G1X,G31),
    update_state(G31,C32H,G2X,G3),
    debug(rosenbrock,"Right side evaluated",[]),

    % compute by backsubstitution
    lusubst(G3,G3X),
    debug(rosenbrock,"G3=~p",[G3X]),

    % for G4 no need to recompute CurrentState and Equations
    C41H is C41/StepSize,
    C42H is C42/StepSize,
    C43H is C43/StepSize,
    update_state(Equations2,C41H,G1X,G41),
    update_state(G41,C42H,G2X,G42),
    update_state(G42,C43H,G3X,G4),
    debug(rosenbrock,"Right side evaluated",[]),

    % compute by backsubstitution
    lusubst(G4,G4X),
    debug(rosenbrock,"G4=~p",[G4X]),

    % now get 4th order estimate for CurrentState
    update_state(TrueCurrentState,B1,G1X,CurrentState41),
    update_state(CurrentState41,B2,G2X,CurrentState42),
    update_state(CurrentState42,B3,G3X,CurrentState43),
    update_state(CurrentState43,B4,G4X,NewCurrentState),
    debug(rosenbrock,"CurrentState 4th order evaluated",[]),

    % Estimate error
    constant_list(Length,0,ErrorListInit), nb_setval(error,ErrorListInit),
    nb_getval(error,ErrorList0),
    update_state(ErrorList0,E1,G1X,ErrorList1),
    update_state(ErrorList1,E2,G2X,ErrorList2),
    update_state(ErrorList2,E3,G3X,ErrorList3),
    update_state(ErrorList3,E4,G4X,ErrorList4),
    scale_error(ErrorList4,TrueCurrentState,ErrorList),
    max_list(ErrorList,Erro),
    debug(rosenbrock,"Error computed",[]),
    debug(rosenbrock,"Erro = ~w",[Erro]),

    % check error
    nb_getval(step_doubling_error,ErrMax),
    Error is Erro/ErrMax,
    debug(rosenbrock,"Error = ~w",[Error]),
592

HEMERY Mathieu's avatar
HEMERY Mathieu committed
593
    (
594
      (
HEMERY Mathieu's avatar
HEMERY Mathieu committed
595 596 597 598 599 600 601 602 603
        Error > 1
      ;
        (Error = nan)   % catch the case where Error=nan
      )
    ->
      % failure: reduce the step
      % TODO if not yet at maxtry
      (
        StepSize < ErrMax/1000
604
      ->
HEMERY Mathieu's avatar
HEMERY Mathieu committed
605 606 607 608 609 610 611 612
        (
          k_fail
        ->
          true
        ;
          assertz(k_fail),
          debug(warning, "Failure with minimal step, error will not be guaranteed.", [])
        )
613
      ;
HEMERY Mathieu's avatar
HEMERY Mathieu committed
614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629
        StepSize1 is SAFE*StepSize*(Error**PSHRINK),
        StepSize2 is SHRINK*StepSize,
        max_list([StepSize1,StepSize2],StepSizeNext),
        nb_setval(h,StepSizeNext),
        debug(rosenbrock,"error ~p, candidates ~p ~p",[Erro,StepSize1,StepSize2]),
        debug(rosenbrock,"step size ~p failed, trying ~p...",[StepSize,StepSizeNext]),
        fail
      )
    ;
      % success, try increasing the step size
      (
        Error > ERRCON*ErrMax
      ->
        StepSizeNext is SAFE*StepSize*(Error**PGROW)
      ;
        StepSizeNext is StepSize*GROW
630
      ),
HEMERY Mathieu's avatar
HEMERY Mathieu committed
631 632 633 634 635 636 637 638 639 640
      debug(rosenbrock,"step size ~p ok, trying ~p...",[StepSize,StepSizeNext]),
      nb_getval(hmax,StepSizeMAX),
      (
        StepSizeNext < StepSizeMAX
      ->
        nb_setval(h,StepSizeNext)
      ;
        true
      )
    ),
641
    %%% EVENTS HANDLING PART %%% 
642
    (
643
      NewTime is Time+StepSize,
644 645 646 647
      check_events(NewCurrentState, NewTime),
      check_conditional_events(NewCurrentState, NewTime),
      nb_getval(events_number,NumberOfEvents),
      NumberOfEvents > 0
648 649 650
    ->
      nb_setval(event_handling, true),
      (
651
        nb_getval(event_step_size, ESS),
652
        StepSize =< ESS
653 654 655 656
      ->
        handle_rosenbrock_events(NewCurrentState, NewTime),
        handle_conditional_events(NewCurrentState, NewTime),
        nb_setval(event_handling, false),
657 658
        nb_getval(h_init, NaturalStepSize),
        nb_setval(h, NaturalStepSize)
659 660 661 662 663 664
      ;
        StepSizeRed is SHRINK*StepSize,
        nb_setval(h, StepSizeRed),
        fail
      )
    ;
665 666 667 668 669 670 671 672 673
      (
        nb_getval(event_handling, true)
      ->
        nb_getval(event_step_size, ESS),
        StepSizeRed is max(SHRINK*StepSize, ESS),
        nb_setval(h, StepSizeRed)
      ;
        true
      )
674
    ),
675
  !,
HEMERY Mathieu's avatar
HEMERY Mathieu committed
676
  true.
677

HEMERY Mathieu's avatar
HEMERY Mathieu committed
678 679 680
%! scale_error(+Errors, +Values, -Scaled_errors)
% 
% Scale error vector with absolute value of variable (if |value| is > 1)
681

682 683 684
scale_error([],[],[]).

scale_error([E|Err],[V|L],[SE|SErr]) :-
HEMERY Mathieu's avatar
HEMERY Mathieu committed
685 686 687 688 689 690 691 692 693
  (
    abs(V) < 1
  ->
    SE is abs(E)
  ;
    SE is abs(E/V)
  ),
  scale_error(Err,L,SErr).

694

695
/*******************/
696
/* Events handling */
697
/*******************/
698

Mathieu Hemery's avatar
Mathieu Hemery committed
699 700
/* Conditionnal event is seen as true when it goes form False to True only */

701
:- dynamic(conditional_events_status/2).
702

703
:- dynamic(rosenbrock_events_status/2).
704

705 706 707 708 709 710 711
%! split_events(+Events, -RegularEvents, -TimeEvents)
%
% Split the events between the time and regular ones, would be deprecated when time_events
% will be correctly handled.

split_events([], [], []).

712 713 714 715 716
split_events([t>=Thr->Do|Tail], RegTail, [t>=Thr->Do|TimeTail]) :-
  !,
  split_events(Tail, RegTail, TimeTail).

split_events([t>Thr->Do|Tail], RegTail, [t>=Thr->Do|TimeTail]) :-
717 718 719 720 721 722 723
  !,
  split_events(Tail, RegTail, TimeTail).

split_events([Head|Tail], [Head|RegTail], TimeTail) :-
  split_events(Tail, RegTail, TimeTail).


Mathieu Hemery's avatar
Mathieu Hemery committed
724 725 726
%! get_events_lists(-EventsList, -ConditionalEventsList)
%
% Gives acces to the events list and the kinetics conditionnal events list
727

728
get_events_lists(EventsList,ConditionalEventsList) :-
Mathieu Hemery's avatar
Mathieu Hemery committed
729 730 731 732 733 734
  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).
735

Mathieu Hemery's avatar
Mathieu Hemery committed
736 737 738 739

%! prepare_conditional_events(+State, +Time)
%
% Creates the initial facts concerning theses events
740

741
prepare_conditional_events(State,Time) :-
Mathieu Hemery's avatar
Mathieu Hemery committed
742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770
  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)
%
%
771

772
prepare_rosenbrock_events(State,Time) :-
Mathieu Hemery's avatar
Mathieu Hemery committed
773 774
  get_events_lists(Events,_),
  forall(
775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791
    (member(Event, Events),Event =.. [->,Condition,AssignmentList]),
    (
      test_condition(State,Time,Condition,Result),
      nth0(EventIndex,Events,Event),
      assertz(rosenbrock_events_status(EventIndex,Result)),
      %format("event new state is ~w ~n",[Result]),
      Result
    ->
      maplist(apply_changes(State,Time), AssignmentList),
      test_condition(State,Time,Condition,NextResult), % In case the condition have changed (hybrid events for example)
      retract(rosenbrock_events_status(EventIndex,_)),
      assertz(rosenbrock_events_status(EventIndex,NextResult))
    ;
      true
    )
  ).

Mathieu Hemery's avatar
Mathieu Hemery committed
792 793 794
%! check_conditional_events(+State, +Time)
%
% Check if there is some conditional events which become true during the time step of the
795
% rosenbrock stimulation.
796

797
check_conditional_events(State,Time) :-
Mathieu Hemery's avatar
Mathieu Hemery committed
798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821
  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)
%
%
822

823
check_events(State,Time) :-
Mathieu Hemery's avatar
Mathieu Hemery committed
824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845
  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
    )
  ).
846

847 848 849 850 851 852 853 854

%! handle_time_event(+State, +CurrentTime, -NextTime)
%
% Apply all time events that should be and indicate the next event time.

handle_time_event(State, CurrentTime, NextTime) :-
  nb_getval(time_events_list, TimeEvents),
  nb_getval(parameters_list, Parameters),
855
  compute_next_event(TimeEvents, State, Parameters, CurrentTime, Time, Effect),
856 857 858 859
  (
    Time =< CurrentTime
  ->
    maplist(apply_changes(State,CurrentTime), Effect),
860 861
    nb_setval(last_time_event, Time),
    handle_time_event(State, CurrentTime, NextTime)
862 863 864 865 866
  ;
    NextTime = Time
  ).


867 868 869 870 871 872 873
%! compute_next_event(-TimeEvents, -State, -Parameters, -Time, +NextTime, +NextEvent)
%
% Compute the next time and type of time_event, return a big number if no event

compute_next_event([], _State, _Parameters, _CurrentTime, 999999, -1).

compute_next_event([(t>=ExprTime->Effect2)|Tail], State, Parameters, CurrentTime, Time0, Effect0) :-
874
  rw_coeff(ExprTime, State, Parameters, CurrentTime, Time2),
875
  compute_next_event(Tail, State, Parameters, CurrentTime, Time1, Effect1),
876
  (
877 878
    nb_getval(last_time_event, LastTime),
    Time2 > LastTime,
879 880 881 882 883 884 885 886 887 888
    Time2 < Time1
  ->
    Time0 = Time2,
    Effect0 = Effect2
  ;
    Time0 = Time1,
    Effect0 = Effect1
  ).


889
handle_conditional_events(State,Time) :-
890
   get_events_lists(_,Events),
891 892 893
   forall(
    (member(Event, Events),Event =.. [->,Condition,AssignmentList]),
    (
894
      test_condition(State,Time,Condition,Result),
895 896 897
      nth0(EventIndex,Events,Event),
      retract(conditional_events_status(EventIndex,EventInfo)),
      assertz(conditional_events_status(EventIndex,Result)),
898
      %format("Eventtt last state is ~w event new state is ~w ~n",[EventInfo,Result]),
899
      \+(EventInfo),
900
      Result
901
    ->
902
      maplist(apply_changes(State,Time), AssignmentList),
903 904 905 906 907 908
      nb_getval(events_number,EventsNumber),
      NewEventsNumber is EventsNumber - 1,
      nb_setval(events_number,NewEventsNumber)
    ;
      true
    )
909
  ).
910

911

912
handle_rosenbrock_events(State,Time) :-
913
   get_events_lists(Events,_),
914
   forall(
915
    (member(Event, Events),Event =.. [->,Condition,AssignmentList]),
916
    (
917
      test_condition(State,Time,Condition,Result),
918 919 920 921 922
      nth0(EventIndex,Events,Event),
      retract(rosenbrock_events_status(EventIndex,EventInfo)),
      assertz(rosenbrock_events_status(EventIndex,Result)),
      %format("Eventtt last state is ~w event new state is ~w ~n",[EventInfo,Result]),
      \+(EventInfo),
923
      Result
924
    ->
925
      maplist(apply_changes(State,Time), AssignmentList),
926 927
      nb_getval(events_number,EventsNumber),
      NewEventsNumber is EventsNumber - 1,
928 929 930 931
      nb_setval(events_number,NewEventsNumber),
      test_condition(State,Time,Condition,NextResult), % In case the condition have changed (hybrid events for example)
      retract(rosenbrock_events_status(EventIndex,_)),
      assertz(rosenbrock_events_status(EventIndex,NextResult))
932 933 934 935 936
    ;
      true
    )
  ).

937 938 939 940 941

%! test_condition(+State, +Time, +Condition, +Result)
%
% Check if the condition to fire an event is true
% Can also be used if conditional changes in events
942

943
test_condition(State,Time,Condition,Result) :-
Mathieu Hemery's avatar
Mathieu Hemery committed
944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977
  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
  ).
978

979
test_condition(State,Time,not Condition,Result) :-
Mathieu Hemery's avatar
Mathieu Hemery committed
980 981 982 983 984 985 986 987
  test_condition(State,Time,Condition,NotResult),
  (
    NotResult
  ->
    Result = false
  ;
    Result = true
  ).
988

989
test_condition(State,Time,Condition,Result) :-
Mathieu Hemery's avatar
Mathieu Hemery committed
990 991
  nb_getval(parameters_list,Parameters), 
  rw_coeff(Condition,State,Parameters,Time,Result).
992

993 994 995 996 997
test_condition(_,_,true,true).

test_condition(_,_,false,false).

test_condition(State,Time,Condition1 and Condition2,Result) :-
Mathieu Hemery's avatar
Mathieu Hemery committed
998 999 1000 1001 1002 1003 1004 1005 1006 1007
  test_condition(State,Time,Condition1,Result1),
  test_condition(State,Time,Condition2,Result2),
  (
    Result1,
    Result2
  ->
    Result = true
  ;
    Result = false
  ).
1008 1009

test_condition(State,Time,Condition1 or Condition2,Result) :-
Mathieu Hemery's avatar
Mathieu Hemery committed
1010 1011 1012 1013 1014 1015 1016 1017 1018
  test_condition(State,Time,Condition1,Result1),
  test_condition(State,Time,Condition2,Result2),
  (
    (Result1;Result2)
  ->
    Result = true
  ;
    Result = false
  ).
1019

Mathieu Hemery's avatar
Mathieu Hemery committed
1020 1021 1022 1023 1024

%! apply_changes(State, Time, Condition)
%
% to be called when an events is triggered

1025
apply_changes(State,Time, Parameter = (if Condition then IfTrue else IfFalse)) :-
Mathieu Hemery's avatar
Mathieu Hemery committed
1026 1027 1028 1029 1030 1031 1032 1033
  test_condition(State,_,Condition,Result),
  (
    Result
  ->
    apply_changes(State,Time,Parameter = IfTrue)
  ;
    apply_changes(State,Time,Parameter = IfFalse)
  ).
1034
   
1035
apply_changes(State,Time, Parameter = Expression) :-
Mathieu Hemery's avatar
Mathieu Hemery committed
1036 1037 1038 1039 1040
  IndexParameter is Parameter + 1,
  nb_getval(parameters_list,Parameters),
  eval_coeff(Expression,State,Parameters,Time,Value),
  matrix:replace(Parameters,IndexParameter,Value,NewParameters),
  nb_setval(parameters_list,NewParameters).
1041

1042
/*:- doc('\\begin{example}Ball simulation using event:').
1043 1044 1045
:- biocham_silent(clear_model).
:- biocham(load_biocham('library:examples/RobustMonitoring/ball.bc')).
:- biocham(numerical_simulation(method: rsbk)).
Mathieu Hemery's avatar
Mathieu Hemery committed
1046 1047
:- doc('\\end{example}').*/