lazy_negation.pl

:- module(
  ln,
  [
    rewrite_ode/2,
    %Commands
    add_reactions_from_pivp/1
  ]).

:- use_module(ode).
:- use_module(gpac).

%! add_reactions_from_pivp(+PIVP)
%
% biocham command
%
% Infer from a given PIVP a biochemical reaction network implementing it
% The input is the PIVP formated as for the gpac module
% 
% Typical usage will be:
% ==
% :- add_reactions_from_pivp((1,(d(y)/dt=z)); (0,(d(z)/dt= (-1*y)))).
% ==

add_reactions_from_pivp(PIVP) :-
   biocham_command,
   type(PIVP, pivp),
   doc('
     creates a reaction model to implement a given PIVP (same syntax as compile_from_pivp),
     variables that may became negative are splitted x => x_p - x_m, the resulting ODE
     is translated in CRN through add_reactions_from_ode_system (see above).
   '),
   rewrite_ode(PIVP,PIVP_conc),
   split_init_ode(PIVP_conc,Init,Ode),
   ode:new_ode_system,
   ode:add_ode(Ode),
   ode:init(Init),
   ode:add_reactions_from_ode_system,
   get_option(fast_rate, Fast),
   parameters:set_parameter(fast, Fast).

:- doc('
  \\begin{example} In this example dual species are introduced for x because of its initial negative value but not for y: \n
').
:- biocham_silent(clear_model).
:- biocham(add_reactions_from_pivp((-1,(d(x)/dt= y-x)); (0,(d(y)/dt=2-y)))).
:- biocham(list_model).
:- biocham(numerical_simulation(time:10)). %, method:msbdf)).
:- biocham(plot).
:- doc('
  \\end{example}
  ').

%! rewrite_ode(+PIVP_string, -PIVP_rewrited)
%
% rewrite the ode set with lazy negated variables

rewrite_ode(PIVP_input, PIVP_output) :-
   detect_negative_variables(PIVP_input, VarNeg),
   rewrite_derivative_loop(PIVP_input, VarNeg, PIVP_output).


%! split_init_ode(+PIVP,-Init,-Ode)
%
% separate a PIVP int the Ode and the Initial Concentration

split_init_ode((Init,Ode),[Name = Init],[Ode]) :-
   Ode = (d(Name)/dt = _Expr).

split_init_ode(PIVP_A;PIVP_B,Init,Ode) :-
   split_init_ode(PIVP_A,Init_A,Ode_A),
   split_init_ode(PIVP_B,Init_B,Ode_B),
   append(Init_A,Init_B,Init),
   append(Ode_A,Ode_B,Ode).


%! detect_negative_variables(+PIVP_string, -NegVar_list)
%
% Detect the variables that need a negative conjugate
% We first search for variables that start in the negative domain
% then recursively check if some variables may go in the negative region.

detect_negative_variables(PIVP_string, NegVar_List) :-
   find_variables(PIVP_string,Var_list,Initial_conc),
   negative_initial_concentration(Var_list,Initial_conc,NegVar_Init),
   find_troubling_variables_loop(PIVP_string, NegVar_Init, NegVar_List).


%! negative_initial_concentration(+Var_list,+Initial_conc,-NegVar_list)
%
% detect the variables the concentration of which are initially negative

negative_initial_concentration([],[],[]).

negative_initial_concentration([Var|VarTail],[Conc|ConcTail],NegVar_list):-
   (
      Conc < 0
   ->
      NegVar_list = [Var|NegVar_tempo]
   ;
      NegVar_list = NegVar_tempo
   ),
   negative_initial_concentration(VarTail,ConcTail,NegVar_tempo).


%! find_troubling_variables_loop(+PIVP,+ListOfNegatedVariables,-NewVarToNegated)
%
% Recursively add the variables that need to be negated

find_troubling_variables_loop(PIVP_string, NegVar_Now, NegVar_Then) :-
   find_troubling_variables_sr(PIVP_string, NegVar_Now, MoreNegVar),
   (
      MoreNegVar = []
   ->
      NegVar_Now = NegVar_Then
   ;
      append(NegVar_Now,MoreNegVar,Both),
      find_troubling_variables_loop(PIVP_string, Both, NegVar_Then)
   ).


%! find_troubling_variables_sr(+PIVP,+ListOfNegatedVariables,-NewVarToNegated)
%
% Given a list of already negated variables, check whether other variables need
% also to be negated.

find_troubling_variables_sr(A;B, LVar, New_Lvar) :-
   find_troubling_variables_sr(A, LVar, VarA),
   find_troubling_variables_sr(B, LVar, VarB),
   append(VarA,VarB,New_Lvar).

find_troubling_variables_sr((_I,(d(X)/dt= E)), LVar, New_Lvar) :-
   (
      \+(member(X,LVar)),
      find_troubling_polynomial(X,E,LVar)
   ->
      New_Lvar = [X]
   ;
      New_Lvar = []
   ).


%! find_troubling_polynomial(+Variable,+Polynomial,+ListOfNegatedVariables)
%
% this predicate check if the variable has to be negated at the polynomial level

find_troubling_polynomial(X,Poly+Monom,LVar) :-
   find_troubling_monomial(X,Monom,LVar);
   find_troubling_polynomial(X,Poly,LVar).

find_troubling_polynomial(X,Poly-Monom,LVar) :-
   find_troubling_monomial(X,-Monom,LVar);
   find_troubling_polynomial(X,Poly,LVar).

find_troubling_polynomial(X,Monom,LVar) :-
   find_troubling_monomial(X,Monom,LVar).


%! find_troubling_monomial(+Variable,+Monomial,+ListOfNegatedVariables)
%
% this predicate check if the variable has to be negated at the monomial level

find_troubling_monomial(X,-Monom,_LVar) :-
   \+(is_in(X,Monom)),!.

find_troubling_monomial(X,Monom,LVar) :-
   \+(is_in(X,Monom)),
   is_in(Y,Monom),
   (
      number(Y) % allow to write stuff like -1*y
   ->
      Y < 0
   ;
      member(Y,LVar)
   ).


%! is_in(+X,+Expression)
%
% check if the variable X is present in the given Expression

is_in(X,Exp1 * Exp2) :-
   is_in(X,Exp1),!;
   is_in(X,Exp2),!.

is_in(X,Var^_Coeff) :-
   is_in(X,Var).

is_in(X,Var) :-
   X = Var,!;
   -X = Var,!.


%! is_negated(+Expression,-UnNegated)
%
% check if the Expression contains a negated sign and return a positive version if necessary

is_negated(-1*A,A) :- !.

is_negated(A*B,Ap*B) :-
   is_negated(A,Ap),!.

is_negated(-A,A) :- !.

is_negated(X,Xp) :-
   number(X),
   X < 0,
   Xp is -X.


%! replace_in(+Var,+Expr,+New_Var,-New_Expr)
%
% replaces all the occurences of Var in Expr by New_Var

replace_in(Var,Var,New_Var,New_Var) :- !.

replace_in(Var,X*Y,New_Var,NX*NY):-
   replace_in(Var,X,New_Var,NX),
   replace_in(Var,Y,New_Var,NY),!.

replace_in(Var,X^Y,New_Var,NX^NY):-
   replace_in(Var,X,New_Var,NX),
   replace_in(Var,Y,New_Var,NY),!.

replace_in(Var,-X,New_Var,-NX):-
   replace_in(Var,X,New_Var,NX),!.

replace_in(_Var,Other_Var,_New_Var,Other_Var).


%! rewrite_derivative_loop(+Derivative,+VarNeg,-Derivatives_New)
%
% Loop the rewrite_derivative subroutine on all the derivatives of the PIVP

rewrite_derivative_loop(A;B,VarNeg,Ab;Bb) :-
   rewrite_derivative_loop(A,VarNeg,Ab),
   rewrite_derivative_loop(B,VarNeg,Bb).

rewrite_derivative_loop(A,VarNeg,Ab) :-
   rewrite_derivative(A,VarNeg,Ab).


%! rewrite_derivative(+Derivative,+VarNeg,-Derivatives_New)
%
% rewrite a single derivative of a PIVP, the two predicates correspond to the cases where
% the variable has to be splitted or not.

rewrite_derivative((I,(d(Var)/dt= E)),VarNeg,(Ip,(d(Varp)/dt= Ep));(In,(d(Varn)/dt= En))) :-
   member(Var,VarNeg),
   name_p_m(Var,Varp,Varn),
   (
      I >= 0
   ->
      Ip is I,
      In is 0
   ;
      Ip is 0,
      In is -I
   ),
   rewrite_polynomial(E,Var,VarNeg,Ep_tempo,En_tempo),
   concatenate_expression(Ep_tempo,(-fast*Varp*Varn),Ep),
   concatenate_expression(En_tempo,(-fast*Varp*Varn),En),!.

rewrite_derivative((I,(d(Var)/dt= E)),VarNeg,(I,(d(Var)/dt= Ep))) :-
   rewrite_polynomial(E,Var,VarNeg,Ep_tempo,En_tempo),
   concatenate_expression(Ep_tempo,-En_tempo,Ep).


%! rewrite_polynomial(+Poly,+Var,+VarNeg,-Poly_p,-Poly_m)
%
% rewrite a polynomial by splitting it and then delegate to rewrite_monomial

rewrite_polynomial(Poly-Monom, Var, VarNeg, Polyp_M, Polym_M) :-
   rewrite_monomial(-Monom, VarNeg, Monom_p, Monom_m),
   rewrite_polynomial(Poly, Var, VarNeg, Poly_p, Poly_m),
   concatenate_expression(Poly_p, Monom_p, Polyp_M),
   concatenate_expression(Poly_m, Monom_m, Polym_M).

rewrite_polynomial(Poly+Monom, Var, VarNeg, Polyp_M, Polym_M) :-
   rewrite_monomial(Monom, VarNeg, Monom_p, Monom_m),
   rewrite_polynomial(Poly, Var, VarNeg, Poly_p, Poly_m),
   concatenate_expression(Poly_p, Monom_p, Polyp_M),
   concatenate_expression(Poly_m, Monom_m, Polym_M).

rewrite_polynomial(Monom, _Var, VarNeg, Monom_p,Monom_n) :-
   rewrite_monomial(Monom, VarNeg, Monom_p, Monom_n).


%! rewrite_monomial(+Monom,+VarNeg,-Monom_pos,-Monom_neg)
%
% return the associated monomial or expression or the term none if it should be null
% rewrite_monomial/4 is the front face for rewrite_monomial/5

rewrite_monomial(Monom,VarNeg,Mp,Mn) :-
   rewrite_monomial(Monom,none,VarNeg,Mp,Mn).

rewrite_monomial(M,none,_VarNeg,Mp2,Mn2) :-
   number(M),!,
   (
      M > 0
   ->
      Mp2 is M,
      Mn2 = none
   ;
      Mp2 = none,
      Mn2 is -M
   ).

rewrite_monomial(Mp,none,VarNeg,Mn2,Mp2) :-
   is_negated(Mp,Mpp),!,
   rewrite_monomial(Mpp,none,VarNeg,Mp2,Mn2).

rewrite_monomial(Mp,Mn,[Var|Tail],Mp2,Mn2) :-
   (is_in(Var,Mp);is_in(Var,Mn)),
   rewrite_monomial(Mp,Mn,Tail,Mp_tempo,Mn_tempo),
   name_p_m(Var,Varp,Varn),
   replace_in(Var,Mp_tempo,Varp,Mp_tempo1),
   replace_in(Var,Mp_tempo,Varn,Mn_tempo1),
   replace_in(Var,Mn_tempo,Varp,Mn_tempo2),
   replace_in(Var,Mn_tempo,Varn,Mp_tempo2),
   concatenate_expression(Mp_tempo1,Mp_tempo2,Mp2),
   concatenate_expression(Mn_tempo1,Mn_tempo2,Mn2),!.

rewrite_monomial(Mp,Mn,[_Var|Tail],Mp2,Mn2) :-
   rewrite_monomial(Mp,Mn,Tail,Mp2,Mn2),!.

rewrite_monomial(Mp,Mn,[],Mp,Mn).


%! concatenate_expression(+Exp1,+Exp2,-Exp12)
%
% Concatenate two expressions in a single one, expressions are monomial, polynomial or the
% none term.

concatenate_expression(Exp1,Exp2,Exp12) :-
   (
      Exp1 = none
   ->
      Exp12 = Exp2
   ;
      Exp2 = none
   ->
      Exp12 = Exp1
   ;
      is_negated(Exp2,Exp2p)
   ->
      Exp12 = Exp1 - Exp2p
   ;
      Exp2 = +(Exp2sp)
   ->
      Exp12 = Exp1 + Exp2sp
   ;
   Exp12 = Exp1+Exp2
   ).


%%% /!\ Do not insert this part, it is already in GPAC %%%

%! find_variables(+PIVP_string, -ListVariableNames, -InitialConcentrations)
%
% Scan a PIVP_string to extract the variables and initial concentrations.

find_variables((I, (d(X)/dt=_)), [X], [I]).

find_variables((X;Y), L, Init):-
   find_variables(X, L1, I1),
   find_variables(Y, L2, I2),
   append(L1, L2, L),
   append(I1, I2, Init).


%! name_p_m(+Name,-Positive_name,-Negative_name)
%
% Produce names for positive and negative terms

name_p_m(X,Xp,Xm):-
   atom_concat(X,'_p',Xp),
   atom_concat(X,'_m',Xm).