Reflow of rosenbrock.pl

rosenbrock.pl

@@ -4,55 +4,59 @@
    ] 
 ).
 
-/* Initialize  coefficients for the rosenborkc method 
-As well as the errors and the step size for euler method */
+%! rosenbrock_init
+%
+% Initialize coefficients for the rosenbrock method, the errors and the step size
+% for euler method
  
 rosenbrock_init :-
-   nb_setval(not_zero_error,0.0001),
-   nb_setval(k_method,stiff),
-   nb_setval(euler_step_size,0.01),
-   nb_setval(step_doubling_error,0.0001),
-   SAFE = 0.9, GROW = 1.5, PGROW = -0.25,
-   SHRINK = 0.5, PSHRINK is -1.0/3.0, ERRCON = 0.1296,
-   %MAXTRY = 40,
-   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,
-   % C1X, C2X, C3X, C4X,
-   % A2X is 1.0, A3X is 3.0/5.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).
+  nb_setval(not_zero_error,0.0001),
+  nb_setval(k_method,stiff),
+  nb_setval(euler_step_size,0.01),
+  nb_setval(step_doubling_error,0.0001),
+  SAFE = 0.9, GROW = 1.5, PGROW = -0.25,
+  SHRINK = 0.5, PSHRINK is -1.0/3.0, ERRCON = 0.1296,
+  %MAXTRY = 40,
+  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,
+  % C1X, C2X, C3X, C4X,
+  % A2X is 1.0, A3X is 3.0/5.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).
 
 :- dynamic(k_fail).
 
 :- dynamic(rsbk_continue/0).
 
-/* Start of the ode solver */
-/* ode_solver(+Options) ; cf numerical_simulation.pl for Options*/
+
+%! ode_solver(+Options)
+%
+% Start of the ode solver, cf numerical_simulation.pl for Options
 
 ode_solver([
     fields: _,
@@ -71,145 +75,163 @@ ode_solver([
     time_final: Duration,
     jacobian: Jacobian
   ]) :- 
-   (
-      nb_current(rsbk_continue,yes) 
-   ->
-      nb_getval(numerical_table,CurrentTable),
-      get_last_row(CurrentTable,InitialTime2,TrueInitialState),
-      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),
+  (
+    nb_current(rsbk_continue,yes) 
+  ->
+    nb_getval(numerical_table,CurrentTable),
+    get_last_row(CurrentTable,InitialTime2,TrueInitialState),
+    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),
                numerical_simulation:convert_identifier(Molecule, [VariableIndex]),
                append(VariableList,[Molecule],NewVariableList),
                nb_setval(variable_list,NewVariableList)
-               )
-            ),
-      nb_getval(variable_list,VariableList),
-      FirstRow =.. [row,'#''Time'|VariableList],
-      eval_state(InitialState,InitialParameters,TrueInitialState),
-      nb_setval(events_list,Events),
-      prepare_rosenbrock_events(TrueInitialState,InitialTime),
-      prepare_conditional_events(TrueInitialState,InitialTime),
-      SecondRow =.. [row,InitialTime|TrueInitialState],
-      InitialTable = [FirstRow,SecondRow],
-      nb_setval(numerical_table,InitialTable),
-      initialize_functions_list,
-      TrueInitialTime is InitialTime
-   ),
-   statistics(runtime,_),
-   rosenbrock_init,
-   rosenbrock(Equations,TrueInitialTime,Duration,Jacobian),
-   statistics(runtime,[_,T]),
-   TT is T/1000,
-   format("Simulation time: ~ps~n",[TT]),
-   retractall(k_fail).
-
-/* Tool to take the first elements of a list */
-/* take(+SubListSize,+List,-SubList) */
-
-take(N, _, Xs) :- N =< 0, !, N =:= 0, Xs = [].
-take(_, [], []).
-take(N, [X|Xs], [X|Ys]) :- M is N-1, take(M, Xs, Ys).
-
-/* Create a global variable which is a list which contains all the functions defined in the .bc file */
+             )
+          ),
+    nb_getval(variable_list,VariableList),
+    FirstRow =.. [row,'#''Time'|VariableList],
+    eval_state(InitialState,InitialParameters,TrueInitialState),
+    nb_setval(events_list,Events),
+    prepare_rosenbrock_events(TrueInitialState,InitialTime),
+    prepare_conditional_events(TrueInitialState,InitialTime),
+    SecondRow =.. [row,InitialTime|TrueInitialState],
+    InitialTable = [FirstRow,SecondRow],
+    nb_setval(numerical_table,InitialTable),
+    initialize_functions_list,
+    TrueInitialTime is InitialTime
+  ),
+  statistics(runtime,_),
+  rosenbrock_init,
+  rosenbrock(Equations,TrueInitialTime,Duration,Jacobian),
+  statistics(runtime,[_,T]),
+  TT is T/1000,
+  format("Simulation time: ~ps~n",[TT]),
+  retractall(k_fail).
+
+%! take(+N, +List, -Sublist)
+%
+% Sublist is the N first element of List
+
+take(Size, List, Sublist) :-
+  prefix(Sublist, List),
+  length(Sublist, Size).
+
+
+%! initialize_functions_list
+%
+% Create a global variable: a list which contains all the functions defined in the .bc file
 
 initialize_functions_list :-
-   all_items([kind: function],Functions),
-   initialize_functions_list(Functions,FunctionsList),
-   (
-      FunctionsList = []
-   ->
-      nb_setval(functions_number,0)
-   ;
-      %format("Function List is ~w ~n",[FunctionsList]),
-      nb_setval(functions_list,FunctionsList),
-      length(FunctionsList,Number),
-      nb_setval(functions_number,Number),
-      initialize_functions
-   ).
+  all_items([kind: function],Functions),
+  initialize_functions_list(Functions,FunctionsList),
+  (
+    FunctionsList = []
+  ->
+    nb_setval(functions_number,0)
+  ;
+    nb_setval(functions_list,FunctionsList),
+    length(FunctionsList,Number),
+    nb_setval(functions_number,Number),
+    initialize_functions
+  ).
 
-initialize_functions_list([],[]).
+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)
+  ).
 
-initialize_functions_list([Function|Functions],FunctionList) :- 
-   Function =.. [function,TrueFunction],
-   TrueFunction =.. [=,Head,Expression],
-   initialize_functions_list(Functions,NextList),
-   (
-      is_kinetics_function(Head)
-   ->
-      FunctionList = NextList
-   ;
-      numerical_simulation:convert_expression(Expression,IndexedExpression),
-      NewFunction = function(Head,IndexedExpression),
-      append(NextList,[NewFunction],FunctionList)
-   ).
 
- /* This preventsthe kinetics functions to be added to the functions list
- Must be updated if others kinetics are added */
+%! 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).
 
-is_kinetics_function(Head) :- % This prevents kinetics function to be in this list
-   Head = 'MA'(k) ; Head = 'MAI'(k) ; Head = 'MM'('Vm', 'Km') ; Head = 'Hill'('Vm', 'Km', n) ;
-   Head = 'HillI'('Vm', 'Km', n).
 
-/* Functions are evaluated at each step of the simulation and the result is stored in the global variable numerical_table
-This predicate adds the name of the functions at the end of the first row of the numerical_table, and the initial values of the functions in the second row */
+%! initialize_functions
+%
+% This predicate adds the name of the functions at the end of the first row of the 
+% numerical_table, and the initial values of the functions in the second row.
 
 initialize_functions :-
-   nb_getval(numerical_table,Table),
-   nb_getval(parameters_list,Parameters),
-   nth0(0,Table,FirstRow),
-   nth0(1,Table,SecondRow),
-   FirstRow =.. [row|VariablesNames],
-   SecondRow =.. [row,0|InitialState],
-   nb_getval(functions_list,Functions),
-   findall(
-           Variable,
-           member(function(Variable,_),Functions),
-           FunctionsVariables
-   ),
-   append(VariablesNames,FunctionsVariables,NewVariables),
-   findall(
-           InitialValue,
-           (member(function(_,Expression),Functions),
-            eval_coeff(Expression,InitialState,Parameters,0,InitialValue1),
-            simplify(InitialValue1,InitialValue2),
-            InitialValue is InitialValue2
-           ),
-           InitialValues
-   ),
-   append(InitialState,InitialValues,NewInitialState),
-   NewFirstRow =.. [row|NewVariables],
-   NewSecondRow =.. [row,0|NewInitialState],
-   NewTable = [NewFirstRow,NewSecondRow],
-   nb_setval(numerical_table,NewTable).
+  nb_getval(numerical_table,Table),
+  nb_getval(parameters_list,Parameters),
+  nth0(0,Table,FirstRow),
+  nth0(1,Table,SecondRow),
+  FirstRow =.. [row|VariablesNames],
+  SecondRow =.. [row,0|InitialState],
+  nb_getval(functions_list,Functions),
+  findall(
+          Variable,
+          member(function(Variable,_),Functions),
+          FunctionsVariables
+  ),
+  append(VariablesNames,FunctionsVariables,NewVariables),
+  findall(
+          InitialValue,
+          (member(function(_,Expression),Functions),
+           eval_coeff(Expression,InitialState,Parameters,0,InitialValue1),
+           simplify(InitialValue1,InitialValue2),
+           InitialValue is InitialValue2
+          ),
+          InitialValues
+  ),
+  append(InitialState,InitialValues,NewInitialState),
+  NewFirstRow =.. [row|NewVariables],
+  NewSecondRow =.. [row,0|NewInitialState],
+  NewTable = [NewFirstRow,NewSecondRow],
+  nb_setval(numerical_table,NewTable).
+
+
+%! update_functions_state(+State, +Time)
+%
+% Evaluate the functions at each step and add the result at the end of the global
+% numerical_table
 
-/* Evaluate the functions at each step and add the result at the end of list which contains the current state of the simulation */
- 
 update_functions_state(CurrentState,Time) :-
-   nb_getval(numerical_table,Table),
-   nb_getval(parameters_list,Parameters),
-   nb_getval(functions_list,Functions),
-   findall(
-           FunctionValue,
-           (member(function(_,Expression),Functions),
+  nb_getval(numerical_table,Table),
+  nb_getval(parameters_list,Parameters),
+  nb_getval(functions_list,Functions),
+  findall(
+          FunctionValue,
+          (
+            member(function(_,Expression),Functions),
             eval_coeff(Expression,CurrentState,Parameters,Time,FunctionValue1),
             simplify(FunctionValue1,FunctionValue2),
             FunctionValue is FunctionValue2
-           ),
-           FunctionsValues
-   ),
-   append(CurrentState,FunctionsValues,NewCurrentState),
-   NewLastRow =.. [row,Time|NewCurrentState],
-   append(Table,[NewLastRow],NewTable),
-   nb_setval(numerical_table,NewTable).
+          ),
+          FunctionsValues
+  ),
+  append(CurrentState,FunctionsValues,NewCurrentState),
+  NewLastRow =.. [row,Time|NewCurrentState],
+  append(Table,[NewLastRow],NewTable),
+  nb_setval(numerical_table,NewTable).
    
 /* Evaluates a matrix/row/coefficient with informations on current state of the simulation and the current parameters */
 
@@ -377,7 +399,7 @@ rosenbrock(Equations,InitialTime,Duration,Jacobian) :-
       nb_setval(k_time_units,2)
    ),
    repeat,
-   counter(TTTTT), write(TTTTT), nl,
+   %   counter(TTTTT), write(TTTTT), nl,
    nb_getval(parameters_list,Parameters),
    nb_getval(numerical_table,CurrentTable), 
    get_last_row(CurrentTable,Time,CurrentStateInitial),
@@ -946,4 +968,5 @@ euler_step(Equations,CurrentState,StepSize,Parameters,NewCurrentState) :-
 :- biocham_silent(clear_model).
 :- biocham(load_biocham('library:examples/RobustMonitoring/ball.bc')).
 :- biocham(numerical_simulation(method: rsbk)).
-:- doc('\\end{example}').*/
\ No newline at end of file
+:- doc('\\end{example}').*/
+

rosenbrock.plt

+:- use_module(rosenbrock).
+
+:- begin_tests(rosenbrock).
+
+:- end_tests(rosenbrock).