/*
Implement the automatic correction of models loaded from a sbml file
Coder: M. Hemery
*/

:- module(
  correction,
  [
  % Commands
  correct_model/0
  % API
  ]
).

% Only for separate compilation/linting
:- use_module(arithmetic_rules).
:- use_module(doc).
:- use_module(formal_derivation).
:- use_module(reaction_rules).
:- use_module(util).

correct_model :-
  biocham_command,
  doc(''),
  single_model(ModelId),
  setof(
    Molecule,
    identifier_kind(ModelId, Molecule, object),
    List_Molecule
  ),
  % First loop that add/delete/split reactions
  reactions_with_species(List_Molecule, List_Reaction),
  forall(
    member(Reaction, List_Reaction),
    (
      remove_null_kinetic(Reaction),!
    ;
      split_bidirectional_reaction(Reaction),!
    ;
      reverse_reaction(Reaction),!
    ;
      true
    )
  ),
  % Second loop to add catalyst/inhibitor
  reactions_with_species(List_Molecule, Sec_List_Reaction),
  forall(
    member(Reaction, Sec_List_Reaction),
    correct_modifiers(List_Molecule, Reaction)
  ).

%! remove_null_kinetic(+Reaction)
%
% Detect and remove reactions with identically null kinetics

remove_null_kinetic(Expression for Reaction) :-
  (
    ode:substitute_functions(Expression, Expr),
    is_null(Expr)
  ->
    delete_reaction(Expression for Reaction)
  ).

%! split_bidirectional_reaction(+Reaction)
%
% Split a bidirection reaction in two

split_bidirectional_reaction(Expression for Reactant => Product) :-
  (
    ode:substitute_functions(Expression, Expr),
    distribute(Expr, DistExpr),
    DistExpr = Up - Down,
    models:formal_product(P, Reactant => Product),
    present_in_kinetics(P, Down)
  ->
    delete_reaction(Expression for Reactant => Product),
    add_reaction(Up for Reactant => Product),
    add_reaction(Down for Product => Reactant)
  ).

%! reverse_reaction(+Reaction)
%
% Reverse a reaction if its kinetic is negative

reverse_reaction(Expression for Reactant => Product) :-
  (
    ode:substitute_functions(Expression, - RevExpr)
  ->
    delete_reaction(Expression for Reactant => Product),
    add_reaction(RevExpr for Product => Reactant)
  ).


%! correct_modifiers(+List_Molecule, +Reaction)
%
% Ensure that the modifiers are correctly formulated

correct_modifiers(List_Molecule, Expression_raw for Reaction) :-
  ode:substitute_functions(Expression_raw, Expression),
  detect_modifiers(List_Molecule, Expression, Catalysts, Inhibitors),
  (
    append(Catalysts, Inhibitors, [])
  ->
    true
  ;
    add_catalysts(Catalysts, Reaction, ReactionTempo),
    add_inhibitors(Inhibitors, ReactionTempo, NewReaction),
    delete_reaction(Expression_raw for Reaction),
    add_reaction(Expression_raw for NewReaction)
  ).


detect_modifiers([], _Expression, [], []).

detect_modifiers([Mol|TailM], Expression, NewC, NewI) :-
  derivate(Expression, Mol, Deriv),
  (
    is_null(Deriv)
  ->
    NewC = TailC,
    NewI = TailI
  ;
    always_negative(Deriv)
  ->
    NewC = TailC,
    NewI = [Mol|TailI]
  ;
    NewC = [Mol|TailC],
    NewI = TailI
  ),
  detect_modifiers(TailM, Expression, TailC, TailI).

add_catalysts([], Reaction, Reaction).
add_catalysts([Head|Tail], Reaction, NewReaction) :-
  models:formal_reactant(Head, Reaction),
  !,
  add_catalysts(Tail, Reaction, NewReaction). 
add_catalysts([Head|Tail], Reactants => Products, NewReactants => NProducts+Head) :-
  add_catalysts(Tail, Reactants => Products, NReactants => NProducts),
  (
    NReactants = Reac / Inhib
  ->
    NewReactants = Reac+Head / Inhib
  ;
    NewReactants = Reactants+Head
  ).

% TODO -- remove Head from Reactants & Products if necessary
add_inhibitors([], Reaction, Reaction).
add_inhibitors([Head|Tail], Reaction, NewReaction) :-
  models:formal_inhibitor(Head, Reaction),
  !,
  add_inhibitors(Tail, Reaction, NewReaction). 
add_inhibitors([Head|Tail], Reactants => Products, NewReactants => NProducts) :-
  add_inhibitors(Tail, Reactants => Products, ReactantsTempo => ProductsTempo),
  (
    remove_molecule(ReactantsTempo, Head, NReactants),
    remove_molecule(ProductsTempo, Head, NProducts)
  ;
    NReactants = ReactantsTempo,
    NProducts = ProductsTempo
  ),
  (
    NReactants = Reac / Inhib
  ->
    NewReactants = Reac / (Head,Inhib)
  ;
    NewReactants = NReactants / Head
  ).

% remove_molecule(Species, Molecule, Without)

remove_molecule(Molecule, Molecule, '_') :- !.
remove_molecule(P*Molecule, Molecule, Result) :- 
  (
    P = 1
  ->
    Result = '_'
  ;
    P = 2
  ->
    Result = Molecule
  ;
    Pm is P-1,
    Result = Pm*Molecule
  ),!.
remove_molecule(Species+Molecule, Molecule, Species) :- !.
remove_molecule(Species+P*Molecule, Molecule, Result) :-
  (
    P = 1
  ->
    Result = Species
  ;
    P = 2
  ->
    Result = Species+Molecule
  ;
    Pm is P-1,
    Result = Species+Pm*Molecule
  ),!.
remove_molecule(Molecule+Other, Molecule, Result) :- !,
  remove_molecule(Other+Molecule, Molecule, Result).
remove_molecule(P*Molecule+Other, Molecule, Result) :- !,
  remove_molecule(Other+P*Molecule, Molecule, Result).
remove_molecule(Species+Other, Molecule, Without+Other) :-
  remove_molecule(Species, Molecule, Without).