/**
A rewriting of the lazy negation module to be used with the gpac representation of a
PIVP, the purpose was to be able to negate the PIVP after the binomial reduction.
For the record, gpac PIVP are represented as pivp_list that is:
pivp_list - a list description of the PIVP used for internal computation
in the format [N_species,[reactions_s1,reactions_s2,...],Initial_concentration]
where reactions_s1 is a list of lists of the form
[k, [exponent s1,..., exponent sn]]
@author M. Hemery
@license GNU-GPL v.2
*/
:- module(
lng,
[
rewrite_PIVP/3
]).
:- use_module(util).
%! rewrite_PIVP(+PIVP_list, -PIVP_rewrited, -VarNeg)
%
% rewrite the ode set with variables negated is needed
% VarNeg is the list of variable that have been negated given as their position in
% the pivp_list (before rewriting).
rewrite_PIVP([N,ODE,Init], [NewN,NewODE,NewInit], VarNeg) :-
negative_initial_concentration(Init, VarNeg_Init),
find_troubling_variables_main(ODE, VarNeg_Init, VarNeg),
rewrite_derivative_main(ODE, VarNeg, ODE_Tempo),
clean_ODE(ODE_Tempo,VarNeg,NewODE),
rewrite_initial_concentration(Init, VarNeg, NewInit),
length(VarNeg, NVN), NewN is N+NVN.
%! rewrite_PIVP_all_negated(+PIVP_list, -PIVP_rewrited, -VarNeg)
%
% rewrite the ode set with all variables negated
rewrite_PIVP_all_negated([N,ODE,Init], [NewN,NewODE,NewInit], VarNeg) :-
numlist(1,N,VarNeg), %VarNeg contains all the variables
rewrite_derivative_main(ODE, VarNeg, ODE_Tempo),
clean_ODE(ODE_Tempo,VarNeg,NewODE),
rewrite_initial_concentration(Init, VarNeg, NewInit),
NewN is 2*N.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Detection part %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%! negative_initial_concentration(+Initial_conc, -VarNeg_list)
%
% detect the variables with a negative initial concentration
negative_initial_concentration(Init,NVL) :-
negative_initial_concentration(Init,NVL,1).
negative_initial_concentration([],[],_N).
negative_initial_concentration([Conc|ConcTail],NVL,N):-
(
substitute([input], [1.0], Conc, ConcA),
ConcA < 0
->
NVL = [N|NVL_tempo]
;
NVL = NVL_tempo
),
Np is N+1,
negative_initial_concentration(ConcTail,NVL_tempo,Np).
%! find_troubling_variables_main(+ODE, +VarNeg_Current, -VarNeg_Update)
%
% Detect the variables that need to be negated through recursive call to
% find_troubling_variables_sr
find_troubling_variables_main(ODE, VarNeg_Cur, VarNeg_Update) :-
find_troubling_variables_sr(ODE, VarNeg_Cur, New_VarNeg),
(
New_VarNeg = []
->
VarNeg_Cur = VarNeg_Update,!
;
append(VarNeg_Cur, New_VarNeg, Both),
find_troubling_variables_main(ODE, Both, VarNeg_Update)
).
%! find_troubling_variables_sr(+ODE, +VarNeg_Current, -VarNeg_New)
%
% Given a list of already negated variables, check whether other variables need
% also to be negated.
find_troubling_variables_sr(ODE,ListOfNegatedVariables,NewVarToNegated) :-
find_troubling_variables_sr(ODE,ListOfNegatedVariables,NewVarToNegated,1).
find_troubling_variables_sr([], _LNVar, [], _N).
% skip the already negated variables
find_troubling_variables_sr([_Poly|Tail], LNVar, New_LNvar, N) :-
member(N,LNVar),!,
Np is N+1,
find_troubling_variables_sr(Tail, LNVar, New_LNvar, Np).
find_troubling_variables_sr([Poly|Tail], LNVar, New_LNvar, N) :-
(
find_troubling_polynomial(N,Poly,LNVar)
->
New_LNvar = [N|New_LNvar_Tempo]
;
New_LNvar = New_LNvar_Tempo
),
Np is N+1,
find_troubling_variables_sr(Tail, LNVar, New_LNvar_Tempo, Np).
%! find_troubling_polynomial(+Variable, +Polynomial, +VarNeg_Current)
%
% true if the variable has to be negated at the polynomial level
find_troubling_polynomial(X, [Monom|Poly], LVar) :-
find_troubling_monomial(X, Monom, LVar),!;
find_troubling_polynomial(X, Poly, LVar).
%! find_troubling_monomial(+Variable,+Monomial,+VarNeg_Current)
%
% true if the variable has to be negated at the monomial level
find_troubling_monomial(X,[Rate,Exponent],_LVar) :-
Rate < 0,
nth1(X, Exponent, 0),
!.
find_troubling_monomial(X,[_Rate,Exponent],LVar) :-
nth1(X,Exponent,0),
member(Y,LVar),
nth1(Y,Exponent,N),
N > 0,!.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Rewriting part %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%! rewrite_initial_concentration(+Init, +VarNeg, -Init_New)
%
% Well, what did you expect?
% Delegate to rewrite_initial_concentration/4
rewrite_initial_concentration(Init, VarNeg, Init_New) :-
rewrite_initial_concentration(Init, VarNeg, Init_New, 1).
rewrite_initial_concentration([], _VarNeg, [], _N).
rewrite_initial_concentration([RawConc|Tail], VarNeg, Init_New, N) :-
(
member(N,VarNeg)
->
(
parameter_value(input,Input_Value)
->
substitute([input], [Input_Value], RawConc, Conc)
;
Conc = RawConc
),
(
Conc >= 0
->
append([RawConc,0], Init_New_Tail, Init_New)
;
(
catch(ConcN is -RawConc, error(_A,_B), fail)
->
append([0,ConcN], Init_New_Tail, Init_New)
;
append([0,-RawConc], Init_New_Tail, Init_New)
)
)
;
append([RawConc], Init_New_Tail, Init_New)
),
Np is N+1,
rewrite_initial_concentration(Tail, VarNeg, Init_New_Tail, Np).
%! rewrite_derivative_main(+Derivative,+VarNeg,-Derivatives_New)
%
% Loop the rewrite_derivative_sr on all the derivatives of the PIVP starting from
% the highest one to avoid offset errors.
rewrite_derivative_main(ODE, [], ODE) :- !.
rewrite_derivative_main(ODE, ListVarNeg, NewODE) :-
max_list(ListVarNeg, Max),
rewrite_derivative_sr(ODE, Max, ODE_Tempo),
delete(ListVarNeg, Max, ListVarNegRed),
rewrite_derivative_main(ODE_Tempo, ListVarNegRed, NewODE).
%! rewrite_derivative_sr(+ODE, +Negated_Var, -ODE_New)
%
% rewrite all the ODE of a PIVP by spliting Negated_Var in two
% At this level, the equation of the negative parts are left empty, they are filled with
% the negative rates through the clean_ODE/3 predicate.
rewrite_derivative_sr(ODE, Negated_Var, ODE_New) :-
rewrite_derivative_sr(ODE, Negated_Var, ODE_New, 1).
rewrite_derivative_sr([], _Negated_Var, [], _N).
rewrite_derivative_sr([Poly|Tail], Negated_Var, [PolyMod|TailMod2], N) :-
rewrite_polynomial(Poly, Negated_Var, PolyMod),
(
N is Negated_Var
->
append([[]],TailMod,TailMod2)
;
TailMod2 = TailMod
),
Np is N+1,
rewrite_derivative_sr(Tail, Negated_Var, TailMod, Np).
%! rewrite_polynomial(+Poly,+Negated_Var,-PolyMod)
%
% rewrite a polynomial by succesive delegation to rewrite_monomial
rewrite_polynomial([], _Negated_Var, []).
rewrite_polynomial([Monom|Poly], Negated_Var, PolyMod) :-
rewrite_monomial(Monom, Negated_Var, MonoPoly),
rewrite_polynomial(Poly, Negated_Var, PolyTail),
append(MonoPoly, PolyTail, PolyMod).
%! rewrite_monomial(+Monom,+Negated_Var,-Poly)
%
% return the polynomial associated to a monomial through variable negation
% Note that the parity of the exponent determine the sign of the second term
% Note also that as we allways have one of the two var that is absent, the cross terms
% may be ignored.
rewrite_monomial([Rate,Expo], Negated_Var, [[Rate,Expand]]) :-
nth1(Negated_Var, Expo, 0),
rewrite_exponent(Expo, Negated_Var, Expand, _Dummy),!.
rewrite_monomial([Rate,Expo], Negated_Var, [M1,M2]) :-
nth1(Negated_Var,Expo,EN),
rewrite_exponent(Expo, Negated_Var, Expo1, Expo2),
(
mod(EN,2) is 0
->
M1 = [Rate,Expo1],
M2 = [Rate,Expo2]
;
NegRate is -Rate,
M1 = [Rate,Expo1],
M2 = [NegRate,Expo2]
).
%! rewrite_exponent(+Expo, +Negated_Var, -Expo_p, -Expo_n)
%
% Return the two exponents corresponding to the separation of Negated_Var
% e.g.: rewrite_exponent([1,2,3], 2, [1, 2, 0, 3], [1, 0, 2, 3]).
rewrite_exponent([E|Tail],Negated_Var,[Ep|Tailp],[En|Tailn]) :-
(
Negated_Var > 1
->
Negated_Varm is Negated_Var-1,
Ep = E, En = E,
rewrite_exponent(Tail,Negated_Varm,Tailp,Tailn)
;
Ep = E, append([0],Tail,Tailp),
En = 0, append([E],Tail,Tailn)
).
%! clean_ODE(+ODE, +List_Var_Neg, -NewODE)
%
% Move the monomials with a negative rate to the negative variable if necessary to form
% an equivalent and positive ode.
% Work by scannig the different variables and when it founds a pair of npos/neg variable,
% it gathers all the monomials with gather_poly/3 and dispatch them in the convenient
% way with dispatch_poly/3
clean_ODE(ODE, List_Var_Neg, NewODE) :-
clean_ODE(ODE, List_Var_Neg, NewODE, 1).
clean_ODE([], _VarNeg, [], _N) :- !.
clean_ODE(ODE, List_Var_Neg, NewODE, N) :-
(
member(N, List_Var_Neg)
->
[P1,P2|Tail] = ODE,
[NP1,NP2|NewTail] = NewODE,
gather_poly(P1,P2,P12),
dispatch_poly(P12,NP1,NP2)
;
[Poly|Tail] = ODE,
[Poly|NewTail] = NewODE
),
Np is N+1,
clean_ODE(Tail, List_Var_Neg, NewTail, Np).
gather_poly(Poly1,[],Poly1) :- !.
gather_poly(Poly1,[[R,E]|Tail],Poly12) :-
NR is -R,
append([[NR,E]], Poly12_Tail, Poly12),
gather_poly(Poly1, Tail, Poly12_Tail).
dispatch_poly([],[],[]) :- !.
dispatch_poly([[R,E]|Tail],NP1,NP2) :-
(
R > 0
->
NP1 = [[R,E]|NP1_Tail],
NP2 = NP2_Tail
;
NR is -R,
NP1 = NP1_Tail,
NP2 = [[NR,E]|NP2_Tail]
),
dispatch_poly(Tail, NP1_Tail, NP2_Tail).