Commit c117f67e authored by Thierry Martinez's avatar Thierry Martinez
parents 9f66b109 10fad23b
:- use_module(library(clpfd)).
:- use_module(library(lists)).
fd_max_integer(268435455).
cpu_time(T) :-
statistics(walltime, [T, _]).
% avoid useless symmetrical rules
:- dynamic(invar_done/2).
% store one in, one out places/transitions
:- dynamic(singular/1).
% store alternative paths
:- dynamic(path/2).
% store variables searched for invariants
:- dynamic(vars/1).
% store base
:- dynamic(base_mol/1).
%% find_all_pinvar
%
% find the complete set of minimal P-invariants
find_all_pinvar :-
fd_max_integer(MaxInt),
find_pinvar(MaxInt).
find_pinvar :-
find_pinvar(4).
%% find_pinvar(+ForcedMax)
%
% find a set of minimal P-invariants with a given max value
find_pinvar(ForcedMax) :-
find_invar_aux(ForcedMax, 'place', 'transition', '#=').
find_invar_aux(ForcedMax, Type, OtherType, Operator) :-
cpu_time(Time1),
findall(
(Ins, Outs),
(
Goal =.. [OtherType, T],
Goal,
findall(
(W, P),
arc(P, T, W),
Ins
),
findall(
(W, P),
arc(T, P, W),
Outs
)
),
L
),
mark_singular(Type),
(
find_invar(L, ForcedMax, Operator),
!
;
write('No complex invariant found\n')
),
% singletons
vars(Vars),
writeall(
(
Goal2 =.. [Type, P],
Goal2,
\+(member(P, Vars)),
assertz(base_mol([P]))
),
[P]
),
cpu_time(Time2),
Time is Time2 - Time1,
format_debug(3, "Time: ~w ms~n", [Time]).
%% mark_singular(+Type)
%
% mark singular places/transitions (depending on Type)
% cf. Soliman WCB'08 or MATHMOD'09
mark_singular(Type) :-
retractall(singular(_)),
(
Goal =.. [Type, E],
Goal,
% weight 1 necessary?
findall(Before, arc(Before, E, 1), [_]),
findall(After, arc(E, After, 1), [_]),
assertz(singular(E)),
fail
;
true
).
%% find_invar(+IOList, +ForcedMax, +Operator)
%
% find invariants for a list of weighted (In, Out) and a given max value
% Operator defines the relation (=, =<, >=) between left and right side for
% sub/sur-invariants
find_invar(IOList, ForcedMax, Operator) :-
retractall(invar_done(_, _)),
remove_common(IOList, UIOList),
format_debug(8, "~w no_common -> ~w~n", [IOList, UIOList]),
% add constraints X.M=0
get_constraints(UIOList, [], [], Vars, VarList, 1, MaxDomain1, Operator),
retractall(vars(_)),
assertz(vars(Vars)),
MaxDomain is min(MaxDomain1, ForcedMax),
% add a large upper bound for minimal invariants
VarList ins 0..MaxDomain,
format_debug(4, "~w~n", [MaxDomain]),
format_debug(4, "~w~n~w~n", [Vars, VarList]),
build_sum(VarList, Sum),
format_debug(8, "Sum: ~w~n", [Sum]),
!,
% avoid the trivial solution
Sum #> 0,
format_debug(9, "Sum #> 0~n", [Sum]),
nb_setval(base, []),
format_debug(9, "looking for paths~n", []),
find_paths,
format_debug(9, "found paths~n", []),
% zero the alternative of each path
zero_paths(Vars, VarList),
format_debug(9, "zeroed paths~n", []),
% reorder vars by neighborhood
% FIXME compare with reverse(VarList, ListVar),
reorder(Vars, VarList, ListVar),
format_debug(9, "reordered vars~n", []),
(
repeat,
(
get_base_constr(VarList),
cpu_time(Time1),
% does NOT guarantee a base e.g.
% add_rule(2*Y + Z => 3*X). for labeling not ff
% my_labeling(VarList)
label(ListVar)
->
cpu_time(Time2),
nb_getval(base, B),
(
MaxDomain > 1
->
make_base(B, VarList, BB)
;
BB = [VarList | B],
format_debug(6, "New minimal invariant: ~w~n", [VarList])
),
Time is Time2 - Time1,
nb_setval(base, BB),
length(BB, LBB),
format_debug(2, "new base length: ~w~n", [LBB]),
format_debug(2, "labeling time: ~w~n", [Time]),
fail
;
retractall(base_mol(_)),
expand_base(Vars),
nb_setval(nb, 0),
writeall(
(
base_mol(Mol),
nb_getval(nb, N),
NN is N + 1,
nb_setval(nb, NN),
sort(Mol, SMol)
),
SMol),
nb_getval(nb, LB),
format("~w complex invariant(s)~n", [LB])
)
->
true
).
%% expand_base(+Vars)
%
% unfolds symmetries in the current base of invariants, using the list of
% variables Vars
expand_base(Vars) :-
nb_getval(base, B),
member(VL, B),
to_mols(VL, Vars, Mol),
format_debug(4, "Mol ~w~n", [Mol]),
assertz(base_mol(Mol)),
fail.
expand_base(_) :-
normalized_path(E, L),
format_debug(5, "Path ~w ~w~n", [E, L]),
findall(
[E | MMol],
(
base_mol(Mol),
format_debug(5, "applying ~w to ~w~n", [L, Mol]),
subtract(Mol, L, MMol),
format_debug(5, "-> ~w~n", [MMol])
),
NewBases
),
format_debug(5, "Newbases ~w~n", [NewBases]),
assert_to_base(NewBases),
fail.
expand_base(_).
%% normalized_path(-E, +P)
%
% get a path starting from E and such that it follows other pathes to their
% end.
normalized_path(E, P) :-
path(E, L),
norm_path_rec(L, P).
norm_path(E, P) :-
path(E, L),
!,
norm_path_rec(L, P).
norm_path(E, [E]).
norm_path_rec([], []).
norm_path_rec([H | T], L) :-
norm_path(H, P),
append(P, Q, L),
norm_path_rec(T, Q).
%% assert_to_base(+L)
%
% add to the base the list of molecules L
assert_to_base([]).
assert_to_base([H | T]) :-
assertz(base_mol(H)),
assert_to_base(T).
%% find_paths, find_paths(CurrentLength, MaximumLength)
%
% find parallel singular paths (of length between CurrentLength and
% MaximumLength
find_paths :-
retractall(path(_, _)),
find_paths(1, 8).
find_paths(Length, _Maxlength) :-
format_debug(8, "looking for paths of length ~w~n", [Length]),
singular(E),
arc(Before, E, 1),
arc(E, After, 1),
(
get_path(Before, After, E, Length, Path)
->
retract(singular(E)),
assertz(path(E, Path)),
format_debug(3, "~w -> ~w~n", [E, Path]),
fail
).
find_paths(Length, MaxLength) :-
(
Length < MaxLength
->
NewLength is Length + 1,
find_paths(NewLength, MaxLength)
;
true
).
%% get_path(Before, After, Current, Length, Path)
%
% look for a path
get_path(A, B, E, 1, [C]) :-
!,
arc(A, C, 1),
C \= E,
singular(C),
arc(C, B, 1).
get_path(A, B, E, N, [C | L]) :-
arc(A, C, 1),
C \= E,
singular(C),
arc(C, D, 1),
NN is N - 1,
get_path(D, B, E, NN, L).
%% zero_paths(+Vars, ?VarList), zero_rec(+Path, +Vars, ?VarList)
%
% force all components of a pth to be equal to 0
zero_paths(Vars, VarList) :-
findall(E, path(E, _), L),
zero_rec(L, Vars, VarList).
zero_rec([], _, _).
zero_rec([M | L], Vars, VarList) :-
nth1(N, Vars, M),
!,
nth1(N, VarList, V),
V = 0,
zero_rec(L, Vars, VarList).
zero_rec([_ | L], Vars, VarList) :-
zero_rec(L, Vars, VarList).
%% get_constraints(+L, -Vars, -VarList, -NVars, -NVarList, -MaxDomain,
%% -NMaxDomain, +Operator)
%
% From a list of rules (KR, KL) get involved molecules, associate a variable
% list and the product of highest coefficients (for upper bound)
get_constraints([], Vars, VarList, Vars, VarList, MaxDomain, MaxDomain, _Op).
% if handling equality, remove symmetrical rules
get_constraints([(KL, KR) | L], Vars, VarList, NVars, NVarList, MaxDomain, NMaxDomain, '#=') :-
invar_done(KR, KL),
!,
get_constraints(L, Vars, VarList, NVars, NVarList, MaxDomain, NMaxDomain, '#=').
get_constraints([(KL, KR) | L], Vars, VarList, NVars, NVarList, MaxDomain, NMaxDomain, Operator) :-
assertz(invar_done(KL, KR)),
make_sum(KL, SL, Vars, VarList, Vars1, VarList1, Max1),
make_sum(KR, SR, Vars1, VarList1, Vars2, VarList2, Max2),
% set upper bound
MaxDomain3 is MaxDomain * max(max(1, Max1), max(1, Max2)),
(
% overflow :(
MaxDomain3 < 0
->
fd_max_integer(N),
MaxDomain2 = N
;
MaxDomain2 = MaxDomain3
),
Constr =.. [Operator, SL, SR],
Constr,
% add constraint X.M = O
format_debug(3, "~w~n", [Constr]),
format_debug(3, "Max:~w = ~w * max(~w, ~w)~n~w~n",
[MaxDomain2, MaxDomain, Max1, Max2, Constr]),
get_constraints(L, Vars2, VarList2, NVars, NVarList, MaxDomain2,
NMaxDomain, Operator).
%% make_sum(+L, -Sum, +Vars, +VarList, -Vars1, -VarList1, -Max)
%
% get molecule and variable list with maximum coeff
make_sum([], 0, V, VL, V, VL, 0).
make_sum([(S, M)], VV, Vars, VarList, Vars1, VarList1, S) :-
(
S = 1
->
VV = V
;
VV = S * V
),
(
nth1(N, Vars, M)
->
nth1(N, VarList, V),
Vars1 = Vars,
VarList1 = VarList
;
fd_max_integer(MaxInt),
V in 0..MaxInt,
Vars1 = [M | Vars],
VarList1 = [V | VarList]
),
!.
make_sum([(S, M) | L], VV + Sum, Vars, VarList, Vars1, VarList1, Max) :-
(
S = 1
->
VV = V
;
VV = S * V
),
(
nth1(N, Vars, M)
->
nth1(N, VarList, V),
Vars2 = Vars,
VarList2 = VarList
;
fd_max_integer(MaxInt),
V in 0..MaxInt,
Vars2 = [M | Vars],
VarList2 = [V | VarList]
),
make_sum(L, Sum, Vars2, VarList2, Vars1, VarList1, Max1),
Max is Max1 + S,
!.
%% to_mols(+ValueList, +VarList, -MolList)
%
% replace instantiated variables by the corresponding named and stoechiometric
% object (a conservation law for a P-inv)
to_mols([], [], []).
to_mols([0 | T], [_ | Vars], Mols) :-
!,
to_mols(T, Vars, Mols).
to_mols([1 | T], [V | Vars], [V | Mols]) :-
!,
to_mols(T, Vars, Mols).
to_mols([N | T], [V | Vars], [N * V | Mols]) :-
to_mols(T, Vars, Mols).
%% get_base_constr(?VarList)
%
% Add constraints on the variables that no base vector should be overlapped
get_base_constr(VarList) :-
nb_getval(base, B),
format_debug(3, "Base: ~w~n", [B]),
get_base_constr_rec(B, VarList),
format_debug(3, "got constr...~n", []).
get_base_constr_rec([], _).
get_base_constr_rec([B | BL], VarList) :-
build_prod(B, VarList, _Prod, Supp, _NotSupp),
element(_, Supp, 0),
get_base_constr_rec(BL, VarList).
%% build_prod(ValueList, ?VarList, -Prod, -S, -NS) :-
%
% builds the product of the variables, forgetting stoichiometry
% records list of the variables participating in the product in S
% and not participating in the product in NS
build_prod([], [], 1, [], []).
build_prod([0 | B], [V | VarList], Prod, S, [V | NS]) :-
!,
build_prod(B, VarList, Prod, S, NS).
build_prod([_ | B], [V | VarList], Prod * V, [V | S], NS) :-
build_prod(B, VarList, Prod, S, NS).
%% make_base(+B, +V, -BB)
%
% remove in B vectors subsumed by V, then add it and return BB
make_base(T, H, [H | T]) :-
% if new vector has max value 1, cannot subsume an earlier vect
% FIXME if MaxDomain=1 this check is useless...
max_list(H, 1),
format_debug(5, "skipping subsumption test for ~w~n", [H]),
!.
make_base([], V, [V]).
make_base([H | T], V, B) :-
(
bigger_support(H, V)
->
format_debug(3, "removing a vector from the base~n", []),
make_base(T, V, B)
;
B = [H | B1],
make_base(T, V, B1)
).
%% bigger_support(V1, V2)
%
% succeeds if V1 has a bigger support (non-null components) than V2
bigger_support([], []).
bigger_support([_ | V], [0 | B]) :-
!,
bigger_support(V, B).
bigger_support([N | V], [_ | B]) :-
!,
N > 0,
bigger_support(V, B).
%% remove_common(+IOList, -UIOList)
%
% Remove (stoechiometrically) catalyzers in a list of (I, O) pairs
% eauch pair is a list of stoichiometric pairs
remove_common([], []).
remove_common([(I, O) | IOList], [(UI, UO) | UIOList]) :-
remove_common(I, O, UI, UO),
remove_common(IOList, UIOList).
remove_common([], KR, [], KR).
remove_common([(SL, M) | KL], KR, UKL, UKR) :-
(
select((SR, M), KR, KR1)
->
(
SL > SR
->
SL2 is SL - SR,
UKL = [(SL2, M) | KL1],
remove_common(KL, KR1, KL1, UKR)
;
(
SL < SR
->
SR2 is SR - SL,
UKR = [(SR2, M) | KR2],
remove_common(KL, KR1, UKL, KR2)
;
remove_common(KL, KR1, UKL, UKR)
)
)
;
UKL = [(SL, M) | KL1],
remove_common(KL, KR, KL1, UKR)
).
%% reorder(+Values, ?VarList1, ?VarList2)
%
% reorders VarList1 into VarList2 s.t. most neighbors (computed using Values)
% are taken first
reorder(Vars, L1, L2) :-
format_debug(8, "Original order: ~w~n", [Vars]),
% initialize neighbors to 0 for everyone
make_keylist(Vars, KV),
reorder_aux(KV, L1, L2).
%% reorder_aux(+KeyValues, ?VarList1, ?VarList2)
%
% same as above but with a Key-Value list representing the current neighbors
reorder_aux([], [], []).
reorder_aux(KL, VL, [V | L]) :-
keysort(KL, SKL),
% sort is ascending, so biggest is last...
last(SKL, _-H),
format_debug(9, "~w ", [H]),
nth1(N, KL, _-H),
% get corresponding variable
nth1(N, VL, V),
remove_nth(N, KL, KT),
remove_nth(N, VL, VT),
update_neighbors(KT, H, KTT),
reorder_aux(KTT, VT, L).
%% make_keylist(+List, -KeyList)
%
% transforms List into a Key-Value pairs list, with all keys equal to 0
make_keylist([], []).
make_keylist([H | T], [0-H | TT]) :-
make_keylist(T, TT).
%% update_neighbors(+L, A, -LL)
%
% udates the Key-value list L into LL adding neighbors of A
update_neighbors([], _, []).
update_neighbors([N-A | L], B, [M-A | LL]) :-
(
(
arc(A, X, _),
arc(B, X, _)
;
arc(A, X, _),
arc(X, B, _)
;
arc(X, A, _),
arc(B, X, _)
;
arc(X, A, _),
arc(X, B, _)
)
->
M is N + 1
;
M = N
),
update_neighbors(L, B, LL).
% store the Petri Net
:- dynamic(place/1).
:- dynamic(transition/1).
:- dynamic(arc/3).
%% format_debug(+Level, +String, +Args)
%
% print debug messages if debug above Level
format_debug(Level, String, Args):-
catch(
nb_getval(debug, N),
error(existence_error(variable, debug), _),
N = 0
),
(