%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Agostino Dovier, August 30, 2012
%% TESTED with BProlog
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Input = [ a, g, c, u, a, ...]

rna_align(Input,Pairing,Energy) :-
     length(Input,N),
     length(Pairing,N),
     rnadomain(Pairing,1,Input), 
     constraint(Input,Pairing),
     energy(Pairing,Input,Energy, N),
     statistics(runtime,_),
     labeling([ff,minimize(Energy)],Pairing),
     statistics(runtime,[_,T]),
     format("Running time: ~d ms~n",T),
     myprint(Input,Pairing).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% SET THE DOMAIN, ACCODING TO THE NUCLEOTIDES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rnadomain([],_,_).
rnadomain([A|R],I,Input) :- %% 0 stands for no pair
     nth1(I,Input,Base),
     (Base==a -> 
         findall(J,(nth1(J,Input,u),abs(I-J)>1),DOMLIST );
      Base==c -> 
         findall(J,(nth1(J,Input,g),abs(I-J)>1),DOMLIST );
      Base==g -> 
         findall(J,(or(J,Input,c,u),abs(I-J)>1),DOMLIST );
      Base==u -> 
         findall(J,(or(J,Input,a,g),abs(I-J)>1),DOMLIST );
      true -> (write('wrong input!'),nl,fail)),
     A in [0|DOMLIST],
     I1 is I+1,
     rnadomain(R,I1,Input).
         
or(V,LIST,A,_) :- nth1(V,LIST,A).
or(V,LIST,_,B) :- nth1(V,LIST,B).
         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% CONSTRAINTS (main procedure)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         
constraint(_,Pairing) :-
     different(Pairing),
     symmetric(Pairing),     
     %contiguity(Pairing),
     no_knots(Pairing),
     write('Constraints added. Starting search.'),nl.
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Injectivity of contacts
%%%% (it is not all_different due to 0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

different([]).
different([_]).
different([A,B|R]) :-
     different_aux(A,[B|R]),
     different([B|R]). 
different_aux(_A,[]).
different_aux(A,[B|R]) :-
     A #> 0 #=> A #\= B,
     different_aux(A,R).
              
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Symmetry of contacts
%%%% If Pairing[i]=j > 0 Then Pairing[j]=i :
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Expressed using "element"
%%% Assume Pairing[i] = 0  (P below)
%%% Then element(1,[i,...],i) is true
%%% Assume Pairing[i] > 0 (e.g. 2 = P)
%%% Then element(3,[i,_,Pairing[2],_, ...],i) is true iff Pairing[2]=i
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

symmetric(Pairing) :-
    symmetric(Pairing,Pairing,1).
symmetric([],_,_).
symmetric([P|Airing],Pairing,I) :-
     P1  #= P+1, element(P1,[I|Pairing],I), %%% <-Symmetry     
     I1 is I + 1,
     symmetric(Airing,Pairing,I1).
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Contiguity (partial) of contacts
%%%% If Pairing[i]=j > 0 
%%%%     Then Pairing[i+1]=0 or it is j+1 or j-1 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Observe that it is symmetrical...
%%% Difficult to be stated with the monotonic encoding
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

contiguity([]).
contiguity([_]).
contiguity([A,B|Pairing]) :-
    A #>0 #=>  (B #= A-1 #\/ B #= A+1 #\/ B #= 0),
    contiguity([B|Pairing]).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% no_pseudo_knots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

no_knots(Pairings) :-
   no_knots(Pairings,1,Pairings).

no_knots(L,_,_) :- length(L,NL),NL<3,!.
no_knots([End|Pairing],Init,Pairs) :-
    %% I only consider increasing pairs (not the symmetrical) 
    truthlist(Init,End,Pairs,Truth),
    (End #> Init) #=> Truth #= 1,
    I1 is Init + 1,
    no_knots(Pairing,I1, Pairs).

truthlist(Init,End,Pairs,Flag) :-
    myselect(Init,Pairs,InterestingPart),
    J is Init+2,
    loopcontrol(InterestingPart,J,Init,End,Truthlist),
    Flag #= min(Truthlist).


loopcontrol([],_,_,_,[]).
loopcontrol([P|Rest],J,Init,End,[T|Truth]) :-
     %%% (we know that (Init,End) in P), with Init < End
     %%% and J > Init.
     %%% We do not consider here backwards jumps
     P #< J                   #=> T #= 1,
    (P #>= J #/\ P #<  End)   #=> T #= 1, 
    (P #>= J #/\ P #>= End)   #=> T #= 0,         
    J1 is J+1,
    loopcontrol(Rest, J1, Init, End, Truth).
    
    
myselect(1,[_,_|Pairs],Pairs).
myselect(N,[_|Pairs],InterestingPart) :-
     N>1,
     N1 is N-1,
     myselect(N1,Pairs,InterestingPart).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% energy function (used a statistical one, as
%%% done by Dahl/Bavarian )
%%% AU~35%, CG~53%, UG~12%
%%% NIL = number of zeros
%%% Energy ~ k1 | AU / (N-NIL) - 0.35 | +
%%%          k2 | CG / (N-NIL) - 0.53 | +
%%%          k3 NIL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                       
energy(Pairing,Input,Energy, N) :-
     findall(J,nth1(J,Input,a),Alist),
     findall(J,nth1(J,Input,c),Clist),
     findall(J,nth1(J,Input,g),Glist),
     findall(J,nth1(J,Input,u),Ulist),
     %%%%%%%     
     cons_energy(Pairing,Input,Pairing, 1,  0,AU, 0,CG,   Alist,Clist,Glist,Ulist),
     count(0,Pairing,#=,NonContacts),
     %%%% PARAMETERS
     C1=400,
     C2=5,
     C3=3,
     %%%%
     Energy #= C1*NonContacts +
               C2*abs(100*AU -  70*N  +  70*NonContacts) + 
               C3*abs(100*CG - 106*N  + 106*NonContacts).

cons_energy([],_,_,_, AU,AU,CG,CG,_,_,_,_).
cons_energy([P|Airing],Input,Ass, I, AUo,AUn,CGo,CGn, 
            Aset,Cset,Gset,Uset) :-
     % write('#######'),write(I), 
     nth1(I,Input,Base),
     % write('#######'),write(Base),nl,
     (Base==a -> 
         (P #= 0 #=> AUi #= AUo,
          P #> 0 #=> AUi #= AUo+1,
          CGo #= CGi);
      Base==c -> 
         (P #= 0 #=> CGi #= CGo,
          P #> 0 #=> CGi #= CGo+1,
          AUo #= AUi);
      Base==g -> 
         (P #= 0 #=> CGi #= CGo,
          P in  Cset #=> CGi #= CGo+1,
          P in  Uset #=> CGi #= CGo,
          AUo #= AUi);
      Base==u -> 
         (P #= 0 #=> AUi #= AUo,
          P in Aset #=> AUi #= AUo+1,
          P in Gset #=> AUi #= AUo,
          CGo #= CGi)),
     I1 is I+1,     
     cons_energy(Airing,Input,Ass, I1, AUi,AUn,CGi,CGn, 
                 Aset,Cset,Gset,Uset).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Printing predicate
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

myprint(Input,Pairing) :-
      myprint(1,Input,Pairing,Input).
myprint(_,[],[],_) :- nl.
myprint(N,[I|Nput],[0|Pairing], Bases) :-
     !,
     write(N),write(':'),write(I),nl, 
     N1 is N+1,
     myprint(N1,Nput,Pairing,Bases).
myprint(N,[I|Nput],[J|Pairing], Bases) :-
     nth1(J,Bases,B),
     write(N),write(':'),write(I),write('--'),
     write(B),write(':'),write(J),nl,
     N1 is N+1,
     myprint(N1,Nput,Pairing,Bases).
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% call as
%%% :- rna(0,L),rna_align(L,Pairing,Energy). 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rna(0,[g,c,g,a,a,a,c,g,c]).
rna(1,[u,c,g,a,a,a,c,g,c,c,c,a]).
rna(2,[a,c,g,a,a,a,u,c,g,a,a,a,c,g,c,c,c,a,u,u,u,u,g,u]).
%%% prion
rna(3,[c,c,a,a,g,a,u,g,u,g,g,a,g,g,c,u,g,g,g,g,u,c,a,g,c,c,c,u,u,g,g,u,g,g,a,g,g,u,u,g,c]).
%%% top riboswitch:
rna(4,[g,a,g,c,a,a,g,g,g,g,u,u,c,u,a,g,u,g,c,a,g,c,u,c,g,a,c,u,a,g,a,a,c,c,c,u,g,c,a,a,u,g,g,c,a,c,g,g,a,c,a,u,c,g,a,u,a,a,g,u,g,u,g,u,u,c,g,u,g,c,u,u,u,u,g,a,c,a,u]).
