mergining

Merge branch 'develop' of gitlab.inria.fr:lifeware/biocham into develop

library/examples/modalINF473L/satode.pl

@@ -15,6 +15,7 @@
 
 :- module(satode, [sat/2, read_dimacs/1]).
 :- reexport(minibiocham).
+:- use_module(symbolic).
 
 %! sat(+N:integer, +M:integer) is det.
 %
@@ -26,8 +27,9 @@ sat(N, M) :-
     (
       make_name('x', I, Xi),
       assertz(ode_init(Xi, 0.5)),
-      big_sum(I, M, Gi),
-      assertz(ode_dxdt(Xi, Gi)),
+      big_sum(I, M, N, Gi),
+      simplify(Gi, GGi),
+      assertz(ode_dxdt(Xi, GGi)),
       make_name('s', I, Si),
       assertz(ode_func(Si, 2 * Xi - 1))
     )
@@ -39,8 +41,9 @@ sat(N, M) :-
       make_name('k', J, Kj),
       assertz(ode_init(Aj, 1)),
       assertz(ode_dxdt(Aj, Aj * Kj)),
-      big_prod(J, N, Pj),
-      assertz(ode_func(Kj, Pj / 2^N))   % Could normalize with 2^k for k-SAT
+      big_prod(J, N, 0, Pj),
+      simplify(Pj, PPj),
+      assertz(ode_func(Kj, PPj))
     )
   ).
 
@@ -52,15 +55,18 @@ make_name(Prefix, I, Name) :-
   format(atom(Name), "~w_~d", [Prefix, I]).
 
 
-%! big_prod(+J:integer, +N:integer, -Pj:term) is det.
+%! big_prod(+J:integer, +N:integer, +Exclude, -Pj:term) is det.
 %
-% Pj is the term corresponding to the prod of all 1 - CjiSi
-% i.e. unnormalized Kj
-big_prod(J, N, Pj) :-
+% Pj is the term corresponding to the prod of all 1 - CjiSi except Exclude
+% and normalized with 1/2^N
+big_prod(J, N, Exclude, Pj) :-
+  %% Start is 1 / (2^N),   % Could normalize with 2^k for k-SAT
+  Start is 1,
   findall(
     C,
     (
       between(1, N, I),
+      I \= Exclude,
       make_name('c', J, Cj),
       make_name(Cj, I, Cji),
       make_name('s', I, Si),
@@ -68,17 +74,16 @@ big_prod(J, N, Pj) :-
     ),
     L
   ),
-  foldl(prod, L, 1, Pj).
+  foldl(prod, L, Start, Pj).
 
 
 prod(X, Y, Y*X).
 
 
-%! big_sum(+I:integer, +M:integer, -Gj:term) is det.
+%! big_sum(+I:integer, +M:integer, +N:integer, -Gj:term) is det.
 %
 % Gj is the term corresponding to the sum of all AjCjiKjiKj
-big_sum(I, M, Gi) :-
-  make_name('s', I, Si),
+big_sum(I, M, N, Gi) :-
   findall(
     C,
     (
@@ -87,7 +92,8 @@ big_sum(I, M, Gi) :-
       make_name(Cj, I, Cji),
       make_name('a', J, Aj),
       make_name('k', J, Kj),
-      C = Aj * Cji * Kj^2 / (1 - Cji * Si)
+      big_prod(J, N, I, Kji),
+      C = Aj * Cji * Kj * Kji
     ),
     L
   ),