Merge branch 'release/4.4.14' into develop

Makefile

@@ -199,7 +199,7 @@ gadagne_restart: gadagne
 	ssh -t gadagne 'cd /opt/tmpnb/ && sudo ./cleanup.sh && sudo ./tmpnb_biocham.sh'
 
 release:
-	echo "current version: " && grep '^version' about.pl
+	@echo "current version: " && grep '^version' about.pl
 	@read -p "Version number (e.g. 4.1.0): " version && \
 	git flow release start "$$version" && \
 	sed -i"" -e "s/^version('.*')/version('$$version')/" about.pl && \

about.pl

@@ -8,7 +8,7 @@
     about/0
   ]).
 
-version('4.4.13').
+version('4.4.14').
 
 copyright(
   'Copyright (C) 2003-2020 Inria, EPI Lifeware, Saclay-Île de France, France'

arithmetic_rules.pl

@@ -389,7 +389,7 @@ insert_coef(A, Coef, Result) :-
     Coef > 0,
     Coef < 1,
     CoefInv is 1 / Coef,
-    0.0 is float_fractional_part(CoefInv)
+    0.0 =:= float_fractional_part(CoefInv)
   ->
     normalize_number(CoefInv, CoefInvNorm),
     Result = A / CoefInvNorm

binder/Dockerfile

-FROM registry.gitlab.inria.fr/lifeware/biocham:v4.4.13
+FROM registry.gitlab.inria.fr/lifeware/biocham:v4.4.14

binomial_reduction.pl

@@ -37,6 +37,10 @@ reduction_methods(no).
 
 reduction_methods(native).
 
+reduction_methods(fast).
+
+reduction_methods(fastnSAT).
+
 reduction_methods(sat_species).
 
 reduction_methods(sat_reactions).
@@ -123,7 +127,6 @@ list_reduction([N,PODE,IV], Name_list, [NewN,NewPODE,NewIV], Modified_Names) :-
       Order > 2
    ->
       generate_sufficient_variables(PODE, ListAllVariable),
-      length(ListAllVariable, Len), writeln(Len),
       breadth_search(PODE, ListAllVariable, ListVariable),
       length(ListVariable,NewN),
       % binomial_branch_bound(ListAllVariable, PODE, ListVariable, NewN),
@@ -230,6 +233,8 @@ highest([],[],[]).
 
 
 % strictly_smaller(+ExpM,+Exp)
+%
+% Exp is smaller than ExpM
 
 strictly_smaller(ExpM,Exp) :-
    smaller(ExpM,Exp),
@@ -551,6 +556,21 @@ breadth_search(PODE, ListVariable, Redvar) :-
      atom_concat('sat', _, Reduction)
    ->
      maxsat_binomial_reduction(ListVariable, ListDerivative, First_Var, Redvar)
+   ;
+     member(Reduction, [fast, fastnSAT])
+   ->
+     fast_search(PODE, First_Var, ListVariable, Redvar1),
+     (
+       Reduction = fast
+     ->
+       Redvar = Redvar1
+     ;
+       construct_all_derivatives(PODE, Redvar1, ListDerivative1),
+       with_option(
+         binomial_reduction:sat_species,
+         maxsat_binomial_reduction(Redvar1, ListDerivative1, First_Var, Redvar)
+       )
+     )
    ;
      Test_list = [[First_Var]],
      breadth_search_subroutine(Test_list, ListVariable, ListDerivative, _Next_level, Redvar_t0),
@@ -746,7 +766,7 @@ convert_to_ode(_Id, [], [], _VL) :- !.
 convert_poly_to_ode([Monomial], VarList, Expr) :-
   !,
   convert_mono_to_ode(Monomial, VarList, Expr).
-convert_poly_to_ode([Monomial|Tail], VarList, Expr+Remainder) :-
+convert_poly_to_ode([Monomial|Tail], VarList, Remainder+Expr) :-
   !,
   convert_mono_to_ode(Monomial, VarList, Expr),
   convert_poly_to_ode(Tail, VarList, Remainder).
@@ -1061,3 +1081,138 @@ is_molecule(A) :-
 is_numeric(N) :-
   item([kind:parameter, key:N]);
   number(N).
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%% FAST BINOMIAL REDUCTION %%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% fast_search(+PODE, +First_var, +ListAllVariable, -BinomialSubSet)
+%
+% perform a fast but inoptimal reduction
+
+fast_search(PODE, V0, AllVariables, BinomialSubSet) :-
+  delete(AllVariables, V0, TailVariables),
+  % construct_all_derivatives(PODE, TailVariables, ListDerivative),
+  add_variable_to_exponents([], PODE, V0, [], Exponents),
+  fast_search_sr([V0], BinomialSubSet, Exponents, [V0|TailVariables], PODE).
+
+fast_search_sr(Acc, RAcc, [], _LAV, _LD) :-
+  reverse(Acc,RAcc).
+
+fast_search_sr(OldVar, Acc, Exponents, LAV, PODE) :-
+  heuristic_few_expo(Exponents, OldVar, LAV, NewVar, PODE),
+  add_variable_to_exponents(Exponents, PODE, NewVar, OldVar, NewExponents),
+  fast_search_sr([NewVar|OldVar], Acc, NewExponents, LAV, PODE).
+
+
+%! add_variable_to_exponents(+ExpOld, +PODE, +NewVar, +OldVar, -ExpNew)
+%
+% compute the new set of unreachable exponents (ExpNew) once we add NewVar to the set of
+% variables
+
+add_variable_to_exponents(ExpOld, PODE, NewVar, OldVar, ExpNew) :-
+  compute_derivative(NewVar, PODE, Derivative),
+  concatenate_exponents(ExpOld, Derivative, ExpTempo),
+  sieve_exponents(ExpTempo, [NewVar|OldVar], ExpNew).
+
+
+%! heuristic_few_expo(+Exponents, +OldVar, +AllVar, -NewVar, +PODE)
+%
+% determine the NewVar to be added to the set of variable by looking first at all the
+% variable that are useful alone (either because they are in the set exponent, or because
+% added to an old variable they form a useful exponent) if this set is empty, it try any
+% remainding variable. Among this set, it choose the one that will introduce as few new
+% exponents as possible.
+
+heuristic_few_expo(Exponents, OldVar, AllVar, NewVar, PODE) :-
+  maybe_useful(AllVar, Exponents, OldVar, Useful, PODE),
+  (
+    Useful = []
+  ->
+    maybe_useless(AllVar, Exponents, OldVar, Useless, PODE),
+    keysort(Useless, [_K-NewVar|_UselessSortedTail])
+  ;
+    keysort(Useful, [_K-NewVar|_UsefulSortedTail])
+  ).
+
+%! maybe_useful(+AllVar, +Exponents, +OldVar, -ChoiceList, +PODE)
+%
+% Determine the list and cost of all variable that are useful alone
+% Choicelist is a list of the form [Cost-Variable|Tail]
+
+maybe_useful([], _Exponents, _OldVar, [], _PODE).
+
+maybe_useful([Head|TailAll], Exponents, OldVar, [Key-Head|TailCut], PODE) :-
+  (
+    member(Head, Exponents)
+  ;
+    member(Other, OldVar),
+    add_list(Other, Head, Expo),
+    member(Expo, Exponents)
+  ),!,
+  estimate_value(Exponents, PODE, Head, OldVar, Key),
+  maybe_useful(TailAll, Exponents, OldVar, TailCut, PODE).
+
+maybe_useful([_Head|TailAll], Exponents, OldVar, TailCut, PODE) :-
+  maybe_useful(TailAll, Exponents, OldVar, TailCut, PODE).
+
+%! maybe_useless(+AllVar, +Exponents, +OldVar, -ChoiceList, +PODE)
+%
+% Determine the list and cost of all variable that are not already used 
+% Choicelist is a list of the form [Cost-Variable|Tail]
+
+maybe_useless([], _Exponents, _OldVar, [], _PODE).
+
+maybe_useless([Head|TailAll], Exponents, OldVar, [Key-Head|TailCut], PODE) :-
+  \+(member(Head, OldVar)),
+  member(Expo, Exponents),
+  strictly_smaller(Expo, Head),
+  !,
+  estimate_value(Exponents, PODE, Head, OldVar, Key),
+  maybe_useless(TailAll, Exponents, OldVar, TailCut, PODE).
+
+maybe_useless([_Head|TailAll], Exponents, OldVar, TailCut, PODE) :-
+  maybe_useless(TailAll, Exponents, OldVar, TailCut, PODE).
+
+
+%! estimate_value(+Exponents, +PODE, +NewVar, +OldVar, -Cost)
+%
+% Determine the cost (number of new exponents to be computed) of NewVar
+
+estimate_value(ExpOld, PODE, NewVar, OldVar, Value) :-
+  add_variable_to_exponents(ExpOld, PODE, NewVar, OldVar, ExpNew),
+  length(ExpNew, Value).
+
+%! concatenate_exponents(+Exponents, +Derivative, -Exponent2)
+%
+% add all the exponent of the derivative to the list of Exponents
+
+concatenate_exponents(Exponents, Derivative, Exponent2):-
+  extract_exponent(Derivative, TempoDer),
+  append(Exponents, TempoDer, Exponent2).
+
+extract_exponent([], []) :- !.
+
+extract_exponent([[_R,E]|TailD], [E|TailE]) :- extract_exponent(TailD, TailE).
+
+%! sieve_exponents(+Exponents, -LV, +NewExponents)
+%
+% remove from the list of exponents all the ones that are binomially reachable with the
+% set of variables LV
+
+sieve_exponents([], _Var_list, []) :- !.
+
+sieve_exponents([Exp|Tail1], Var_list, Tail2) :-
+  (
+    constant(Exp)
+  ;
+    member(Exp,Var_list)
+  ;
+    member(X, Var_list), member(Y, Var_list),
+    add_list(X, Y, Exp)
+  ),
+  !,
+  sieve_exponents(Tail1, Var_list, Tail2).
+
+sieve_exponents([Exp|Tail1], Var_list, [Exp|Tail2]) :-
+  sieve_exponents(Tail1, Var_list, Tail2).

binomial_reduction.plt

@@ -6,11 +6,11 @@ list_reduction:- use_module(library(plunit)).
 test(polynomial_ODE_True) :-
   command(a => b),
   command(3*a^2 for a+a => c),
-  reduction:polynomial_ODE.
+  with_output_to(atom(_), reduction:polynomial_ODE).
 
 test(polynomial_ODE_False) :-
   command(a/b for a => b),
-  \+( reduction:polynomial_ODE ).
+  with_output_to(atom(_), \+ reduction:polynomial_ODE).
 
 test(generate_names) :-
   reduction:generate_names([a,'Out'], [[1,2],[1,0]], [aOut2, a]).
@@ -84,6 +84,6 @@ test(clean_writing) :-
 
 test(convert_poly_to_ode) :-
   reduction:convert_poly_to_ode([[1, [1,1]]], [a,b], a*b),
-  reduction:convert_poly_to_ode([[1, [1,1]], [-2, [0,2]]], [a,b], a*b + -(2*b^2)).
+  reduction:convert_poly_to_ode([[1, [1,1]], [-2, [0,2]]], [a,b], -(2*b^2) + a*b).
 
 :- end_tests(binomial_reduction).

biocham_jupyter/guinbextension/package.json

 {
   "name": "gui",
-  "version": "4.4.13",
+  "version": "4.4.14",
   "description": "biocham gui in jupyter notebook",
   "main": "src/index.js",
   "scripts": {

biocham_jupyter/guinbextension/src/config/workflows/commands.py

@@ -147,6 +147,7 @@ commands = [
     "list_stable_states",
     "list_tables",
     "list_tscc_candidates",
+    "list_units",
     "load",
     "load_biocham",
     "load_ginml",
@@ -188,7 +189,9 @@ commands = [
     "rate_independence_reduction_inputs_sinks",
     "reaction_graph",
     "reaction_model",
+    "read_xpp",
     "reduce_model",
+    "remove_fraction",
     "rename_column",
     "rename_table",
     "reset_options",
@@ -213,6 +216,7 @@ commands = [
     "set_model_name",
     "set_ode_system_name",
     "set_p_m_rate",
+    "set_units",
     "test_rate_independence",
     "test_rate_independence_inputs_sinks",
     "test_rate_independence_invariants",
@@ -223,4 +227,5 @@ commands = [
     "variation",
     "which_p_m_mode",
     "with_timer",
+    "write_xpp",
 ]

biocham_jupyter/kernel/biocham_kernel/__init__.py

 """Example magic"""
-__version__ = '4.4.13'
+__version__ = '4.4.14'

biocham_jupyter/kernel/biocham_kernel/images.py

@@ -69,6 +69,7 @@ def get_data_and_meta(filename):
     to_show = metadata[1].split(' ')[1:]
     against = metadata[2].split(' ')[1]
     axes = metadata[3].split(' ')[1:5]
+    axes_names = (' '.join(metadata[4].split(' ')[1:]), ' '.join(metadata[5].split(' ')[1:]))
 
     with open(filename, newline='') as csvfile:
         lines = list(csv.reader(csvfile))
@@ -78,7 +79,8 @@ def get_data_and_meta(filename):
             'logscale': logscale,
             'to_show': to_show,
             'against': against,
-            'axes': axes}
+            'axes': axes,
+            'axes_names': axes_names}
 
 
 def create_bokeh_plot(data, gui=False):
@@ -100,7 +102,8 @@ def create_bokeh_plot(data, gui=False):
         bokeh_loaded = True
 
     fig = figure(
-        x_axis_label=against,
+        x_axis_label=data['axes_names'][0],
+        y_axis_label=data['axes_names'][1],
         x_axis_type='log' if 'x' in logscale else 'linear',
         y_axis_type='log' if 'y' in logscale else 'linear',
         tools='pan,wheel_zoom,box_zoom,reset,save,undo',

commandline.pl

@@ -35,7 +35,7 @@ The following files can also be used to export some Biocham objects:
 :- doc('\\section{Biocham call options}').
 
 do_arguments :-
-  current_prolog_flag(argv, [_Exec | Args]),
+  current_prolog_flag(os_argv, [_Exec | Args]),
   do_arguments(Args).
 
 

gpac.pl

@@ -673,9 +673,9 @@ extract_names((_IC, d(Name)/dt = _Deriv; Remainder), [Name|Tail]) :-
 % A monomial is a unique object of the form: [K,[Exp1,Exp2,Exp3]]
 % corresponding to the term K v1^Exp1 V2^Exp2 etc.
 
-:- dynamic(input_name/1,
-  output_name/1,
-  name_spec/0).
+:- dynamic(input_name/1).
+:- dynamic(output_name/1).
+:- dynamic(name_spec/0).
 
 %! compare_exponants(+ExponentSet1,+ExponentSet2,-Comparator)
 %

models.pl

@@ -136,7 +136,7 @@ load(InputFile) :-
   type(InputFile, input_file),
   doc('
     acts as the corresponding \\command{load_biocham/1} / \\command{load_sbml/1}
-    / \\command{load_ode/1} / \\command{load_table/1},
+    / \\command{import_ode/1} / \\command{load_table/1},
     depending on the file extension
     (respectively \\texttt{.bc}, \\texttt{.xml}, \\texttt{.ode}, \\texttt{.csv}
     -- assuming no extension is \\texttt{.bc}).'),

modules/partialfrac/partialfraction.pl

@@ -86,7 +86,7 @@ normalized([], []) :- !.
 
 normalized(A, Anorm) :-
     leading(A, Lc),
-    mult(float(1 rdiv Lc), A, Anorm).
+    mult(float(1 / Lc), A, Anorm).
 
 % subtracts two polynomials
 subtract(A, B, C) :-
@@ -123,7 +123,7 @@ long_division(A, B, Q, R) :-
      leading(A, LcA),
      leading(B, LcB),
      D is DegA - DegB,
-     C is float(LcA rdiv LcB),
+     C is float(LcA / LcB),
      monomial_mult(D, C, B, B1),
      subtract(A, B1, A1),
      long_division(A1, B, Q1, R), !,
@@ -163,7 +163,7 @@ simplify_fraction(A, B, Aprime, Bprime) :-
 normalized_euclid(A, B, Unorm, Vnorm, Dnorm) :-
     euclid(A, B, U, V, D),
     leading(D, LcD),
-    C is float(1 rdiv LcD),
+    C is float(1 / LcD),
     mult(C, U, Unorm),
     mult(C, V, Vnorm),
     mult(C, D, Dnorm), !.
@@ -271,7 +271,7 @@ partial_fraction(P, Q, A, QL, ML) :-
     % the following three lines ensure Q is unitary
     % without altering the fraction P/Q
     leading(Q, Lead),
-    PCoef is float(1 rdiv Lead),
+    PCoef is float(1 / Lead),
     mult(PCoef, P, P_normal),
     
     roots(Q, QLraw, ML),

numerical_simulation.pl

@@ -133,6 +133,7 @@ filter(only_extrema).
 :- initial(option(filter: no_filter)).
 :- initial(option(stats: no)).
 
+:- dynamic(last_method/1).
 
 :- devdoc('\\section{Commands}').
 
@@ -179,6 +180,9 @@ numerical_simulation :-
     'filtering function for the trace'
   ),
   doc('performs a numerical simulation from time 0 up to a given time.'),
+  retractall(last_method(_)),
+  debug(numsim, "numerical simulation method: ~w", [Method]),
+  assertz(last_method(Method)),
   statistics(runtime,_),
   (
     Method = ssa
@@ -203,7 +207,7 @@ numerical_simulation :-
   ;
     Method = rsbk
   ->
-  assertz(rosenbrock_running),
+    assertz(rosenbrock_running),
     with_current_ode_system(
       with_clean(
         [
@@ -218,7 +222,6 @@ numerical_simulation :-
         solve2
       )
     ),
-    assertz(last_method_rsbk),
     retractall(rosenbrock_running),
     retractall(rosenbrock_running2)
   ;
@@ -233,8 +236,7 @@ numerical_simulation :-
         ],
         solve
       )
-    ),
-    retractall(last_method_rsbk)
+    )
   ),
   statistics(runtime,[_,T]),
   (
@@ -258,7 +260,8 @@ continue :-
   get_initial_values(Variables, InitialValues),
   set_initial_values(Variables, Values),
   % FIXME uses global options, not the ones used for previous simulation
-  ( last_method_rsbk
+  (
+    last_method(rsbk)
   ->
     nb_setval(rsbk_continue,yes),
     with_option([time: Time,method: rsbk], numerical_simulation),
@@ -305,6 +308,7 @@ prepare_numerical_simulation_options(Options) :-
   enumerate_variables,
   convert_ode,
   gather_equations(Equations),
+  debug(numsim, "ODEs: ~w", [Equations]),
   (
     rosenbrock_running
   ->
@@ -319,12 +323,16 @@ prepare_numerical_simulation_options(Options) :-
     true
   ),
   gather_initial_values(InitialValues),
+  debug(numsim, "Initial values: ~w", [InitialValues]),
   gather_events(RegularEvents, event),
   gather_events(TimeEvents, time_event),
   append(RegularEvents, TimeEvents, Events),
+  debug(numsim, "Events: ~w", [Events]),
   gather_fields(Events, Fields),
   gather_initial_parameter_values(InitialParameterValues),
+  debug(numsim, "Parameter values: ~w", [InitialParameterValues]),
   jacobian(Equations, Jacobian),
+  debug(numsim, "Jacobian: ~w", [Jacobian]),
   Options = [
     fields: Fields,
     equations: Equations,
@@ -371,7 +379,9 @@ solve2 :-
 
 solve :-
   prepare_numerical_simulation_options(Options),
+  debug(numsim, "simulation options: ~w", [Options]),
   solve(Options, Table),
+  debug(numsim, "Done", []),
   add_table('numerical_simulation', Table).
 
 
@@ -408,6 +418,7 @@ enumerate_nonconstant_parameters(Events, Header: Parameter) :-
     ),
     NCParameters
   ),
+  debug(numsim, "NC parameters ~w", [NCParameters]),
   member(ParameterIndex, NCParameters),
   Parameter = p(ParameterIndex),
   convert_identifier(Header, Parameter).
@@ -417,7 +428,9 @@ enumerate_nonparametric_functions(Header: expression(Expr)) :-
   % (item([kind: function, item: function(Header = Body)]) ;
   get_ode_function(Header, Body),
   atomic(Header),
+  debug(numsim, "NP function ~w with body ~w", [Header, Body]),
   ode:substitute_functions(Body, NoFuncBody),
+  debug(numsim, "substituted body ~w", [NoFuncBody]),
   convert_expression(NoFuncBody, Expr).
 
 
@@ -473,11 +486,13 @@ gather_initial_values(InitialValues) :-
 gather_initial_parameter_values(InitialParameterValues) :-
   peek_count(parameter_counter, ParameterCount),
   ParameterMax is ParameterCount - 1,
+  debug(numsim, "~d parameters", [ParameterCount]),
   findall(
     InitialParameterValue,
     (
       between(0, ParameterMax, ParameterIndex),
       parameter_index(Parameter, ParameterIndex),
+      debug(numsim, "Parameter #~d: ~w", [ParameterIndex, Parameter]),
       (%parameter_value(Parameter, InitialParameterValue)
         ode_parameter_value(Parameter, InitialParameterValue)
       ; conditional_parameter_value(ParameterIndex,InitialParameterValue) % virtual event parameters for conditional kinetics

ode.pl

@@ -22,6 +22,7 @@
     add_reactions_from_ode_system/0,
     add_influences_from_ode_system/0,
     (init)/1,
+    remove_fraction/0,
 %    load_ode_system/1,
     load_reactions_from_ode/1,
 %    add_ode_system/1,
@@ -52,7 +53,9 @@
     ode/2,
     ode/3,
     ode_add_expression_to_molecule/2,
-    ode_predicate/1
+    ode_predicate/1,
+    remove_fraction/2,
+    normalize_ode/1
   ]).
 
 % Only for separate compilation/linking
@@ -124,7 +127,10 @@ list_ode :-
     with_current_ode_system(
       (
         get_current_ode_system(Id),
-        list_items([parent: Id, kind: ode])
+        list_items([parent: Id, kind: ode]),
+        list_items([parent: Id, kind: initial_concentration]),
+        list_items([parent: Id, kind: parameter]),
+        list_items([parent: Id, kind: function])
       )
     )
   ;
@@ -414,7 +420,31 @@ load_influences_from_ode_system :-
   add_influences_from_ode_system.
 
 
-
+remove_fraction :-
+  biocham_command,
+  doc('Remove the rational fraction in the current ODE system by multiplying the derivative
+  by the GCD. WARNING: the resulting ODE system is NOT equivalent to the starting one.'),
+  new_ode_system(NewId),
+  set_ode_system_name(NewId, without_fraction_ode),
+  normalize_ode(ListOde),
+  maplist(remove_fraction, ListOde, NewListOde),
+  maplist(add_ode(NewId), NewListOde),
+  with_current_ode_system((
+    get_current_ode_system(Id),
+    forall(
+      item([parent: Id, kind:parameter, item: Item]),
+      add_item([parent: NewId, kind:parameter, item: Item])
+    ),
+    forall(
+      item([parent: Id, kind: function, item: Item]),
+      add_item([parent: NewId, kind: function, item: Item])
+    ),
+    forall(
+      item([parent: Id, kind:initial_concentration, item: Item]),
+      add_item([parent: NewId, kind:initial_concentration, item: Item])
+    )
+  )),
+  select_ode_system(without_fraction_ode).
 
 
 :- devdoc('\\section{Public API}').
@@ -559,6 +589,18 @@ add_ode(Id, ODE) :-
   add_item([parent: Id, kind: ode, key: X, item: Item]).
 
 
+change_ode(Id, ODE) :-
+  (
+    parse_ode(ODE, X, Item)
+  ->
+    true
+  ;
+    throw(error(illformed_ode(ODE)))
+  ),
+  delete_item([parent: Id, kind: ode, key: X]),
+  add_item([parent: Id, kind: ode, key: X, item: Item]).
+
+
 import_reactions_from_ode_system(Id) :-
   check_cleaned(ode:assoc/2),
   enumerate_terms_in_odes(Id),
@@ -662,25 +704,10 @@ non_decomposable_term_in_expression(Coefficient, Term, E) :-
 
 
 substitute_functions(E, ENoFunc) :-
-  (
-    get_current_ode_system(Id)
-  ->
-    all_items([parent: Id, kind: function], Functions),
-    maplist(split_func, Functions, Definitions, Bodies),
-    maplist(
-      rewrite(substitute(Definitions, Bodies)),
-      Bodies,
-      NoFuncBods
-    ),
-    substitute(Definitions, NoFuncBods, E, ENoFunc)
-  ;
-    % case of add_function, there is no current system, but no functions
-    % either
-    ENoFunc = E
-  ).
-
+  kinetics:eval_kinetics([], [], E, ENoFunc), % with kinetics/4 get some errors on simplification
+  !.
 
-split_func(func(F = D), F, D).
+substitute_functions(E, E). % case of add_function, there is no current system, but no functions
 
 
 non_decomposable_term_in_additive_normal_form(Coefficient, Term, A + B) :-
@@ -1259,3 +1286,185 @@ ode_add_expression_to_molecule(NewExpression, Molecule) :-
   ;
     asserta(assoc(Molecule, NewExpression))