Properly substitute functions for ode_systems

especially if not coming from ode_system/0 (which already does some substitution)

library/examples/modalINF473L/satode.pl

@@ -6,13 +6,15 @@
  * numerical_simulation.
  * plot(show: {x_1, x_2}).
  *
+ * Or "read_dimacs('sat.cnf'), pl2xpp('out.ode')." and then you can directly
+ * run "numerical_simulation. plot."
+ *
  * @author Sylvain Soliman
  * @license GPL2
  */
 
-:- module(satode, [sat/2]).
+:- module(satode, [sat/2, read_dimacs/1]).
 :- reexport(minibiocham).
-:- use_module(symbolic).
 
 %! sat(+N:integer, +M:integer) is det.
 %
@@ -40,8 +42,7 @@ sat(N, M) :-
       big_prod(J, N, Pj),
       assertz(ode_func(Kj, Pj / 2^N))   % Could normalize with 2^k for k-SAT
     )
-  ),
-  print_xpp.
+  ).
 
 
 %! make_name(+Prefix:atom, +I:integer, -Name:atom) is det.
@@ -94,3 +95,55 @@ big_sum(I, M, Gi) :-
 
 
 sum(X, Y, Y+X).
+
+
+read_dimacs(File) :-
+  setup_call_cleanup(
+    open(File, read, Stream),
+    (
+      read_header(Stream, N, M),
+      sat(N, M),
+      forall(
+        between(1, M, J),
+        read_clause(Stream, J)
+      )
+    ),
+    close(Stream)
+  ).
+
+
+read_header(Stream, N, M) :-
+  read_line_to_string(Stream, Line),
+  split_string(Line, " ", " ", ["p", "cnf", Nvar, Nclauses]),
+  string_to_number(Nvar, N),
+  string_to_number(Nclauses, M).
+
+
+read_clause(Stream, Index) :-
+  read_line_to_string(Stream, Line),
+  split_string(Line, " ", " ", VarStrings),
+  maplist(string_to_number, VarStrings, Vars),
+  maplist(set_clause_parameter(Index), Vars).
+
+
+string_to_number(String, Number) :-
+  read_term_from_atom(String, Number, []),
+  number(Number).
+
+
+set_clause_parameter(_Clause, 0) :-
+  !.
+
+set_clause_parameter(Clause, N) :-
+  N > 0,
+  !,
+  make_name('c', Clause, Cj),
+  make_name(Cj, N, Cji),
+  assertz(ode_par(Cji, 1)).
+
+set_clause_parameter(Clause, N) :-
+  N < 0,
+  NN is -N,
+  make_name('c', Clause, Cj),
+  make_name(Cj, NN, Cji),
+  assertz(ode_par(Cji, -1)).

numerical_simulation.pl

@@ -417,7 +417,7 @@ enumerate_nonparametric_functions(Header: expression(Expr)) :-
   % (item([kind: function, item: function(Header = Body)]) ;
   get_ode_function(Header, Body),
   atomic(Header),
-  kinetics([], [], Body, NoFuncBody),
+  ode:substitute_functions(Body, NoFuncBody),
   convert_expression(NoFuncBody, Expr).
 
 
@@ -535,8 +535,9 @@ convert_ode :-
   \+ (
     ode(Molecule, Expression),
     \+ (
+      ode:substitute_functions(Expression, NoFuncExpr),
       variable(Molecule, VariableIndex),
-      convert_expression(Expression, IndexedExpression),
+      convert_expression(NoFuncExpr, IndexedExpression),
       assertz(equation(VariableIndex, IndexedExpression)),
       % also consider as parameters those used in initial values
       % so that they can be searched