Commit 370b31b0 authored by SOLIMAN Sylvain's avatar SOLIMAN Sylvain
Browse files

tropicotine is almost in

parent de8d9550
......@@ -93,6 +93,8 @@
* Quantitative Behavior Specification, Verification and Synthesis
** Numerical data time series
- tables.pl
** Tropical equilibrations
- tropical.pl
** Temporal logic FO-LTL(Rlin)
- foltl.pl
** Robustness of FO-LTL(Rlin) properties
......
:- module(
tropical,
[
tropicalize/0
% grammar
%Public API
]
).
% Only for separate compilation/linting
:- use_module(doc).
:- use_module(biocham).
:- use_module(library(clpfd)).
:- doc('Tropical equilibration search as per \\cite{SFR14almob}.').
% store monomials of -inf degree
:- dynamic(is_zero/1).
% equilibrate(?LPos, ?LNeg)
%
% To equilibrate, the positive term of min degree should be of same degree as
% the neg term of min degree
equilibrate([], []) :-
!.
equilibrate(LPos, LNeg) :-
get_min(LPos, Min),
get_min(LNeg, Min).
% equilibrate_rec(-L)
%
% call equilibrate for each pair (Pos, Neg) of list L
equilibrate_rec([]).
equilibrate_rec([(Pos, Neg) | E]) :-
equilibrate(Pos, Neg),
equilibrate_rec(E).
% get_min(?L, ?M)
%
% set M as smaller than all variables of L
% C is the number of variables of L actually equal to M
%
% since here we use C #>=1, we rely on swipl's min fd expr
% get_min([], +inf).
get_min([A], M) :-
!,
M #= A.
get_min([H | T], M) :-
M #= min(H, MM),
get_min(T, MM).
% get_min(L, M) :-
% get_min_and_count(L, M, C),
% % here we enforce at least 2-term equilibration (i.e. max of positive is
% % equal to max of neg and all others negligible)
% % for strict 2-term equilibration one would have C #= 1.
% C #>= 1.
% get_min_and_count(?L, ?M, ?C)
%
% sets M as smaller than all variables of the list given as first arg
% sets C to the count of variables of that list actually equal to M
get_min_and_count([H | T], M, C) :-
M #=< H,
B #<==> M #= H,
C #= B + CC,
get_min_and_count(T, M, CC).
get_min_and_count([], _, 0).
% conserve(?L, -D)
%
% a conservation law implies that when epsilon goes to infinity, the remaining
% terms (those of minimal degree) should be of degree DD (since the sum is
% equal to a constant D) and at least two (since otherwise it would force one
% variable to be always constant)
conserve(_, 0) :-
% FIXME no idea what to do :-/
!.
conserve(L, D) :-
get_degree(D, DD),
% we do not use get_min here, even with C #>= 1 since L is a long list and
% deeply nested min is not that efficient...
get_min_and_count(L, DD, C),
C #>= 1.
% C #>= 2. % oscillations FIXME add an option to toggle on/off
% set_conservations(?Ai, +Vars)
%
% given a list of variables for the species of the system, use the basis of
% P-invariants already computed to enforce constraints on the exponents Ai
set_conservations(Ai, Places) :-
nb_getval(base, B),
(
B \= 0
->
invariants:vars(Vars),
length(Ai, N),
length(Bi, N),
reorder_place_vars(Ai, Bi, Places, Vars),
get_initial_values(Vars, Di),
set_conservations_rec(B, Bi, Di)
;
true
),
set_constant_vars.
%% reorder_place_vars(-L1, -L2, +L3, +L4)
%
% reorders the list of vars L1, with order L3, as L2 with order L4
reorder_place_vars([], _, [], _).
reorder_place_vars([H | T], L, [HP | P], V) :-
nth0(N, V, HP),
nth0(N, L, H),
!,
reorder_place_vars(T, L, P, V).
%% set_constant_vars
%
% add fake parameters corresponding to initial values of constant variables
set_constant_vars :-
invariants:vars(Vars),
enumerate_molecules(M),
member(V, M),
get_initial_concentration(V, I),
\+ member(V, Vars),
get_degree(I, DI),
assertz(degree(V, DI)),
fail.
set_constant_vars.
% set_conservations_rec(-B, ?Ai, -Di)
%
% same as above but for a given basis B and already computed initial values Di
set_conservations_rec([], _, _).
set_conservations_rec([B | BB], Ai, Di) :-
set_conservation(B, Ai, [], Di, 0),
set_conservations_rec(BB, Ai, Di).
% set_conservation(-B, ?Ai, L, -Di, ?D)
%
% same as above but for a given vector B of the basis and using L as
% accumulator for storing the exponents that are in the support of B
% and D as accumulator for the initial value
set_conservation([0 | BB], [_ | AA], L, [_ | DD], D) :-
!,
set_conservation(BB, AA, L, DD, D).
set_conservation([N | BB], [Ai | AA], L, [Di | DD], D) :-
N > 0,
Dj is N * Di + D,
set_conservation(BB, AA, [Ai | L], DD, Dj).
set_conservation([], [], L, [], D) :-
conserve(L, D).
:- dynamic(degree/2).
% set_degrees(+Epsilon)
%
% given a small Epsilon, compute rescaling degrees for parameters
set_degrees(Epsilon) :-
Epsilon < 1,
nb_setval(epsilon, Epsilon),
E is log(Epsilon),
nb_setval(log_epsilon, E),
retractall(degree(_, _)),
parameter_value(K, V),
get_degree(V, Deg),
assertz(degree(K, Deg)),
fail.
set_degrees(_).
get_degree(N, -inf) :-
N =< 0,
!.
get_degree(D, DD) :-
nb_getval(log_epsilon, LEpsilon),
nb_getval(ovidiu_denominator, O),
% DD is round((1 + O) * log(D)/LEpsilon).
DD is (1 + O) * round(log(D)/LEpsilon).
% subst_deg_and_vars(+E, +Vars, ?Ai, ?D)
%
% substitute parameters by their degree, species by their exponent and
% products by sums in E to get D
subst_deg_and_vars(A * B, Vars, Ai, AA + BB) :-
!,
subst_deg_and_vars(A, Vars, Ai, AA),
subst_deg_and_vars(B, Vars, Ai, BB).
subst_deg_and_vars(A ^ N, Vars, Ai, N * AA) :-
integer(N),
!,
subst_deg_and_vars(A, Vars, Ai, AA).
subst_deg_and_vars(A ^ N, Vars, Ai, NN * AA) :-
float(N),
NN is round(N),
NN =:= N,
!,
subst_deg_and_vars(A, Vars, Ai, AA).
subst_deg_and_vars(K, _Vars, _Ai, D) :-
degree(K, D),
!.
subst_deg_and_vars(V, Vars, Ai, A) :-
nth0(N, Vars, V),
!,
nth0(N, Ai, A).
% TODO use get_degree for numbers
subst_deg_and_vars(E, _, _, _) :-
format_to_atom(A, "We do not understand the kinetic term ~p~n", [E]),
throw(error(A)).
% split_pos_neg_and_subst(+O, +Vars, ?Ai, ?Pos, ?Neg)
%
% split expression O in positive and negative terms of the ODE
% for each term apply substitutions to keep the epsilon exponent
split_pos_neg_and_subst(O, Vars, Ai, Pos, Neg) :-
split_pos_neg_and_subst_rec(O, Vars, Ai, Pos, Neg).
split_pos_neg_and_subst_rec(E + N * T, Vars, Ai, Pos, Neg) :-
!,
split_pos_neg_and_subst_rec(N * T, Vars, Ai, Pos1, Neg1),
append(Pos1, PPos, Pos),
append(Neg1, NNeg, Neg),
split_pos_neg_and_subst_rec(E, Vars, Ai, PPos, NNeg).
split_pos_neg_and_subst_rec(N * (A - B), Vars, Ai, Pos, Neg) :-
!,
M is -N,
split_pos_neg_and_subst_rec(N * A + M * B, Vars, Ai, Pos, Neg).
split_pos_neg_and_subst_rec(N * (A + B), Vars, Ai, Pos, Neg) :-
!,
split_pos_neg_and_subst_rec(N * A + N * B, Vars, Ai, Pos, Neg).
split_pos_neg_and_subst_rec(N * T, Vars, Ai, Pos, Neg) :-
number(N),
!,
ode:substitute_functions(T, TT),
subst_deg_and_vars(TT, Vars, Ai, S),
(
atom_is_subterm(inf, S)
->
assertz(is_zero(T)),
Pos = [],
Neg = []
;
(
N > 0
->
Neg = [],
Pos = [S]
;
Pos = [],
Neg = [S]
)
).
split_pos_neg_and_subst_rec(K, Vars, Ai, Pos, Neg) :-
split_pos_neg_and_subst_rec(1*K, Vars, Ai, Pos, Neg).
% get_equilibrations(+Vars, ?Ai, ?E, ?DE)
%
% E is the list of equilibrations coming from the precomputed odes and taking
% into account that Ai are the exponents of species of Vars
% DE also substracts degree of left side for correct rescaled system
get_equilibrations(Vars, Ai, E, DE) :-
ode_system,
get_current_ode_system(Id),
findall(
(S, O),
(
item([parent: Id, kind: ode, item: I]),
I = (d(S)/dt = O)
),
L
),
get_equilibrations_rec(L, Vars, Ai, E, DE).
get_equilibrations_rec([], _, _, [], []).
get_equilibrations_rec([(S, O) | L], Vars, Ai, [E | EE], [DE | DEE]) :-
nth0(N, Vars, S),
nth0(N, Ai, A),
split_pos_neg_and_subst(O, Vars, Ai, Pos, Neg),
% since we will equalize mins, we do not really need to substract left and
% right...
substract_rec(Pos, A, PPos),
substract_rec(Neg, A, NNeg),
E = (Pos, Neg),
DE = (PPos, NNeg),
debug(tropical, "Monomials for ~w: ~w~n", [S, E]),
get_equilibrations_rec(L, Vars, Ai, EE, DEE).
substract_rec(L, A, LA) :-
maplist(L, reverse_minus(A), LA).
reverse_minus(A, B, B - A).
show_rescaling([], []) :-
nl.
show_rescaling([A | AA], [V | VV]) :-
epsilon_degree(A, D),
format("~p_ = ~p~p~n", [V, D, V]),
show_rescaling(AA, VV).
show_rescaling_ovidiu([D]) :-
!,
nb_getval(ovidiu_denominator, O),
DD is D / (1 + O),
write(ovidiu_file, DD),
nl(ovidiu_file).
show_rescaling_ovidiu([H | T]) :-
nb_getval(ovidiu_denominator, O),
HH is H / (1 + O),
write(ovidiu_file, HH),
write(ovidiu_file, '\t'),
show_rescaling_ovidiu(T).
epsilon_degree(D, A) :-
(
D =\= 0
->
format(atom(A), "ε^~p * ", [D])
;
A = ''
).
rescale_system(E, S) :-
eval_rec(E, EE),
get_rescaled_odes(EE, O),
sort_by_scale_rec(O, OO),
predsort(compare_first_degree, OO, S).
print_system([]) :-
nl.
print_system([(Species, ODE) | L]) :-
format("d~p/dt = ", [Species]),
print_ode(ODE),
print_system(L).
print_ode([]) :-
write('0\n').
print_ode([(D, T) | O]) :-
epsilon_degree(D, A),
print(A),
(
A \= '',
T = _ + _
->
format("(~p)", [T])
;
print(T)
),
(
O == []
->
nl
;
print(' + '),
print_ode(O)
).
compare_first_degree(Result, (S1, [(D1, _) | _]), (S2, [(D2, _) | _])) :-
compare(Result, (D1, S1), (D2, S2)).
get_rescaled_odes([], []).
get_rescaled_odes([(LPos, LNeg) | E], [(Species, O) | OO]) :-
retract(ode(Species, T)),
rescale_to(T, LPos, LNeg, O),
get_rescaled_odes(E, OO).
rescale_to(E + N * T, Pos, Neg, EEE) :-
!,
rescale_to(N * T, Pos1, Neg1, RT),
append(RT, EE, EEE),
append(Pos1, LPos, Pos),
append(Neg1, LNeg, Neg),
rescale_to(E, LPos, LNeg, EE).
rescale_to(N * (A + B), Pos, Neg, EE) :-
!,
rescale_to(N * A + N * B, Pos, Neg, EE).
rescale_to(N * (A - B), Pos, Neg, EE) :-
!,
M is -N,
rescale_to(N * A + M * B, Pos, Neg, EE).
rescale_to(_ * T, [], [], []) :-
is_zero(T),
!.
rescale_to(N * T, LPos, LNeg, [RT]) :-
!,
(
N > 0
->
LPos = [Exp],
LNeg = []
;
LNeg = [Exp],
LPos = []
),
RT = (Exp, N * T).
sort_by_scale_rec([], []).
sort_by_scale_rec([(Species, Ode) | L], [(Species, SOde) | LL]) :-
sort_by_scale(Ode, SOde),
sort_by_scale_rec(L, LL).
sort_by_scale(L, SL) :-
sort_by_scale_aux(L, SL, []).
sort_by_scale_aux([], L, L).
sort_by_scale_aux([H | T], SL, L) :-
merge_by_scale(L, H, LL),
sort_by_scale_aux(T, SL, LL).
merge_by_scale([(N, T1) | L], (N, T2), [(N, T1 + T2) | L]) :-
!.
merge_by_scale([(N, T1) | L], (M, T2), [(N, T1) | LL]) :-
N < M,
!,
merge_by_scale(L, (M, T2), LL).
merge_by_scale(L, (M, T2), [(M, T2) | L]).
eval_rec([], []) :-
!.
eval_rec([H | T], [HH | TT]) :-
!,
eval_rec(H, HH),
eval_rec(T, TT).
eval_rec((A, B), (AA, BB)) :-
!,
eval_rec(A, AA),
eval_rec(B, BB).
eval_rec(A, B) :-
B is A.
% FIXME add rescaled conservation to S1
reduce_system([], [], []).
reduce_system([(Species, [(D, T) | _]) | ODEs], [T = 0 | S1], S2) :-
D < 0,
format("~p is a fast variable with invariant manifold ~p = 0~n", [Species, T]),
reduce_system(ODEs, S1, S2).
reduce_system([(Species, [(D, T) | L]) | ODEs], S1, [(Species, [(D, T) | L]) | S2]) :-
D >= 0,
format("~p is a slow variable~n", [Species]),
reduce_system(ODEs, S1, S2).
% look for an atom as a leaf of a term, not trying to unify FD vars
atom_is_subterm(_, Y) :-
var(Y),
!,
fail.
atom_is_subterm(X, X) :-
!.
atom_is_subterm(_, Y) :-
atomic(Y),
!,
fail.
atom_is_subterm(X, F) :-
F =.. [_ | Args],
member(A, Args),
atom_is_subterm(X, A).
%% molecules_vars(-L)
%
% returns variables of vars(L), i.e. found as constrained by invariants
% but ordered as in molecules(L).
molecules_vars(L) :-
enumerate_molecules(P),
invariants:vars(V),
filter(P, V, L).
%% filter(+L1, +L2, +L3)
%
% L3 is elements of L1 that are also in L2
filter([], _, []).
filter([H | T], V, [H | L]) :-
member(H, V),
!,
filter(T, V, L).
filter([_ | T], V, L) :-
filter(T, V, L).
ignore_unbalanced([], []).
ignore_unbalanced([([_ | _], []) | L], LL) :-
!,
ignore_unbalanced(L, LL).
ignore_unbalanced([([], [_ | _]) | L], LL) :-
!,
ignore_unbalanced(L, LL).
ignore_unbalanced([H | T], [H | TT]) :-
!,
ignore_unbalanced(T, TT).
:- initial(option(tropical_epsilon: 0.1)).
:- initial(option(tropical_max_degree: 3)).
:- initial(option(tropical_partial: 0)).
:- initial(option(tropical_denominator: 0)).
% tropicalize(+F, ?Max, ?Eps)
%
% main predicate for tropical equilibration
% if Max is a variable, use the number of variables as bound on the absolute
% value of the degree in epsilon, otherwise iterate powers of 2 up to 2^Max
tropicalize :-
biocham_command,
doc('Try to solve a tropical equilibration problem.'),
option(
tropical_epsilon,
number,
Eps,
'Used for computing degrees (must be < 1)'
),
option(
tropical_max_degree,
number,
Max,
'Max absolute value of the degree in epsilon as a power of 2'
),
option(
tropical_partial,
number,
Partial,
'Nubmer of a variable to ignore for equilibration (0 to ignore none, -1
to ignore all unbalanced)'
),
option(
tropical_denominator,
number,
Denominator,
'Look for solutions that are integers/(tropical_denominator+1)'
),
nb_setval(ovidiu_denominator, Denominator),
set_degrees(Eps),
retractall(is_zero(_)),
(
% might not be defined so we use nb_current instead of nb_getval
nb_current(ovidiu_file, Ovidiu)
->
open(Ovidiu, write, _, [alias(ovidiu_file)]),
find_pinvar(5)
;
find_pinvar(4)
),