%%% SICStus Prolog Code for Protein Folding Prediction %%% in CLP(FD) %%% Version 030430 :-use_module(library(clpfd)). :-use_module(library(lists)). :-use_module(library(terms)). :-consult('data.txt'). :-consult('table.txt'). %%% (MAIN) PREDICATE: fcc_pf %%% INPUT: ID, the code of the protein according to names in file data.pl %%% OUTPUT: the minimum Energy associated to the conformation %%% represented by Tertiary %%% EXAMPLE: fcc_pf('1ed0'). %%% NOTE: "protein" is defined in data.pl fcc_pf( ID):- clean, %%% remove garbage if any protein(ID, Primary, Secondary), initial_time, constrain(Primary, Secondary, Indexes, Tertiary, Energy,M), writetime('Constraint time: '),nl,!, local_labeling(Primary, Secondary, Indexes, Tertiary, Energy,M), potential_list(Primary, Tertiary, List), sum(List,#=,EnergyTot), writetime('Global time: '),nl,!, write(Indexes),nl, write('*'),nl, write(ID),nl, write(EnergyTot),nl, print_results(Primary,Tertiary), write('*'),nl. %%% shut down Sicstus when a solution is found. fcc_pf( ID, 1):- fcc_pf( ID ), halt. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% CONSTRAINT DEFINITION PART %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% constrain(Primary, Secondary, Indexes, Tertiary, Energy, M):- write('Constrain:'),nl, length(Primary,N), checkamino(Primary), generate_indexes(N, Indexes), generate_tertiary(N, Primary, Tertiary), domain_bounds(N, Tertiary), avoid_self_loops(Tertiary,N), next_constraints(Tertiary), distance_constraints(Tertiary), secondary_info(Secondary, Indexes, Tertiary), avoid_symmetries(Indexes, Tertiary, Secondary), indexes_to_coordinates(Indexes, Tertiary, _Dir), energy_constraints(Primary, Secondary, Tertiary, Indexes, Energy, _ListaPot, M, _MD), write('Secondary:'),write(Secondary),nl, write('Indexes:'),write(Indexes),nl, write('Constrain end.'),nl. %%% PREDICATE: generate_tertiary %%% INPUT: N, a number greater than or equal to 0, %%% Primary, a list of aminoacids %%% OUTPUT: Ter, a list of 3*N variables to be used as coordinates generate_tertiary(N,Amino,Ter):- length(Amino,N), N3 is N*3, length(Ter,N3). %%% check correctness of primary input checkamino([T|R]):- aminoacid(T),!, checkamino(R). checkamino([_T|_R]):- write('Bad aminoacid sequence!.'),nl, fail. checkamino([]). %%% PREDICATE: domain_bounds %%% INPUT: N, a number greater than or equal to 0, %%% Tertiary %%% EFFECT: adding FD constraints to those variables %%% in order they keep value in 2*n fcc domain_bounds(N, [X,Y,Z|Rest]) :- CUBESIZE is N * 2, domain([X,Y,Z],0,CUBESIZE), sum([X,Y,Z], #=, Sum), even(Sum), domain_bounds(N, Rest). domain_bounds(_N, []). %%% PREDICATE: avoid_self_loops (aux.: positions_to_integers) %%% INPUT: Tertiary, the list of variables, N size of protein %%% EFFECT: adding FD constraints to those variables %%% in order to avoid equal triples avoid_self_loops(Tertiary, N):- CubeSize is 2*N+1, positions_to_integers(Tertiary, ListaInteri, CubeSize), all_different(ListaInteri). positions_to_integers([X,Y,Z|R], [I|S], CubeSize):- I #= X*CubeSize*CubeSize+Y*CubeSize+Z, %%% weigthed sum !!! positions_to_integers( R , S, CubeSize). positions_to_integers([],[],_). %%% PREDICATE: next_constraints (aux.: next) %%% INPUT: Tertiary, the list of variables %%% EFFECT: adding FD constraints to those variables %%% in order to set that the triples are consecutive next_constraints([X1,Y1,Z1,X2,Y2,Z2|C]) :- next(X1,Y1,Z1,X2,Y2,Z2), next_constraints([X2,Y2,Z2|C]). next_constraints([_,_,_]). next_constraints([]). next(X1,Y1,Z1,X2,Y2,Z2):- domain([Dx,Dy,Dz],0,1), Dx #= abs(X1-X2), Dy #= abs(Y1-Y2), Dz #= abs(Z1-Z2), Dx+Dy+Dz #= 2. %%% PREDICATE: distance_constraints (aux.: distance, non_next) %%% INPUT: Tertiary, the list of variables' triples %%% EFFECT: adding FD constraints to those variables %%% in order to set that non consecutive triples %%% has lattice distance greater than 1. %%% NOTE: It introduces an explicit disjunction distance_constraints([X1,Y1,Z1,X2,Y2,Z2|R]):- distance(X1,Y1,Z1,R), distance_constraints([X2,Y2,Z2|R]). distance_constraints([_,_,_]). distance_constraints([]). distance(X1,Y1,Z1,[X2,Y2,Z2|R]):- non_next(X1,Y1,Z1,X2,Y2,Z2), distance(X1,Y1,Z1,R). distance(_,_,_,[]). non_next(X1,Y1,Z1,X2,Y2,Z2):- %%% j > i + 1 !!! Dx #= (X1-X2)*(X1-X2), Dy #= (Y1-Y2)*(Y1-Y2), Dz #= (Z1-Z2)*(Z1-Z2), Dx+Dy+Dz #>2. %%% PREDICATE: generate_indexes %%% INPUT: N, a number greater than or equal to 0, %%% OUTPUT: a list of N-2 variables as \theta-indexes %%% 0-5 six torsion angles (0,60,120,180,240,300) generate_indexes(N, Indexes):- M is N-3, M > 0,!, length(Indexes,M), domain(Indexes,0,5). generate_indexes(0, []). generate_indexes(1, []). generate_indexes(2, []). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% AUXILIARY PREDICATES FOR SECONDARY STRUCTURE CONSTRAINTS: %%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% PREDICATE: secondary_info/3 %%% INPUT: The list Secondary, %%% The list Indexes, %%% The list Tertiary %%% EFFECT: Adding constraints to the variables in Indexes and Tertiary %%% according to the information in Secondary %%% NOTE: helix(1,..) or strand(1,..) are not allowed!!! %%% NOTE: helix(_,N) or strand(_,N) are not allowed!!! %%% We do not add +1 because in this way we leave free %%% the beginning and the ending of a chain/strand %%%%%%%%%%%%%%% helix and strand already processed in first_constrain %%% Empty list secondary_info([], _ , _ ). %%% Case 1: ssbond secondary_info([ssbond(A1,A2)|SecRest], Indexes, Tertiary):- !, ssbond_propagate(A1,A2,Tertiary), secondary_info(SecRest, Indexes, Tertiary). %%% Case 2: helix secondary_info([helix( Start, End)|SecRest], Indexes, Tertiary):- !, S is Start - 1, Helix_length #= End - Start, % We do not add + 1 sublist(S, Helix_length, Indexes, Subindexes), helix(Subindexes), secondary_info(SecRest,Indexes,Tertiary). %%% Case 3: strand secondary_info([strand( Start, End)|SecRest], Indexes, Tertiary):- !, S is Start - 1, Helix_length #= End - Start, % We do not add + 1 sublist(S, Helix_length, Indexes, Subindexes), strand(Subindexes), secondary_info(SecRest,Indexes, Tertiary). %%% PREDICATE: helix %%% INPUT: A sublist of Indexes, %%% EFFECT: Adding constraints to the variables involved in the sublist helix(Indexes):- domain(Indexes,5,5). %%% PREDICATE: strand %%% INPUT: A sublist of Indexes, %%% EFFECT: Adding constraints to the variables involved in the sublist strand(Indexes):- domain(Indexes,3,3). isodd(L,1):- odd(L). isodd(_,0). ssbond_propagate(A1,A2,Tertiary):- extract_ter(A1,A2,Tertiary,P1,P2), ssbond(P1,P2,6). %%% PREDICATE: ssbond %%% INPUT: Two triples of variables, C is the constant to increase the distance %%% EFFECT: Adding a "contact" constraint to the variables involved %%% in the sublist ssbond([X1,Y1,Z1],[X2,Y2,Z2],C):- (X1-X2)*(X1-X2) + (Y1-Y2)*(Y1-Y2) + (Z1-Z2)*(Z1-Z2) #=< C*C. avoid_symmetries(Indexes,Tertiary,Secondary):- length(Tertiary,N), N3 #= N/3, isodd(N3,Odd), P is N3+Odd, P1 is P+1, P2 is P+2, choose_subindexes(Indexes, Su1, L, Secondary), sublist(Su1,L,Indexes,Sub), getmaxrun(Sub,Su2,_E), S is Su1 + Su2-1, S1 is S+1, S2 is S+2, extract_ter(S,S1,Tertiary,[P,P,P],[P1,_P1,P1]), extract_ter(S,S2,Tertiary,[P,P,P],[P2,_P1,_P1]). %% selects in Tertiary the triples corresponding %% to amino A1 and A2 extract_ter(A1,A2,Tertiary,[PA1xa,PA1ya,PA1za],[PA2xa,PA2ya,PA2za]):- Ind11a is (A1-1)*3+1, Ind12a is (A1-1)*3+2, Ind13a is (A1-1)*3+3, Ind21a is (A2-1)*3+1, Ind22a is (A2-1)*3+2, Ind23a is (A2-1)*3+3, nth(Ind11a,Tertiary,PA1xa), nth(Ind12a,Tertiary,PA1ya), nth(Ind13a,Tertiary,PA1za), nth(Ind21a,Tertiary,PA2xa), nth(Ind22a,Tertiary,PA2ya), nth(Ind23a,Tertiary,PA2za). getmaxrun(Indexes,S,E):- indexes2segments(Indexes,Runs,1), getmax(Runs,0,0,S,E). getmax([],S,E,S,E). getmax([S,E|R],S1,E1,Q,W):- E-S > E1-S1,!, getmax(R,S,E,Q,W). getmax([_S,_E|R],Me,Ms,Q,W):- getmax(R,Me,Ms,Q,W). %%% PREDICATE: indexes_to_coordinates %%% (aux. indexes_to_directions, directions_to_coordinates, %%% vector_sum) %%% INPUT: The list of Indexes %%% The list of variables's triples Tertiary %%% EFFECT: Relates the \theta-constraints on Indexes to the %%% coordinates variables in Tertiary indexes_to_coordinates([], _, _). indexes_to_coordinates(Indexes, Tertiary, Directions):- directions_to_coordinates(Directions, Tertiary), indexes_to_directions(Indexes, Directions). directions_to_coordinates([D1 | Directions],[X11,X12,X13,X21,X22,X23|Coordinate]):- vector_sum(D1, [X11,X12,X13], [X21,X22,X23]), directions_to_coordinates(Directions,[X21,X22,X23|Coordinate]). directions_to_coordinates([], [_,_,_] ). directions_to_coordinates([], [] ). indexes_to_directions([I1|Indexes],[D0,D1,D2|Directions]):- new_direction(D0,D1,D2,I1), indexes_to_directions(Indexes,[D1,D2|Directions]). indexes_to_directions([], _). %%% PREDICATE: new_direction %%% INPUT: An index,The list of Indexes %%% The list of variables's triples Tertiary %%% EFFECT: Moving the \theta-constraints on Indexes to the %%% coordinates variables in Tertiary new_direction([Ax,Ay,Az], [Bx,By,Bz], [Cx,Cy,Cz], I):- domain([Rx,Ry,Rz,Tipo,Ax,Ay,Az,Bx,By,Bz,Cx,Cy,Cz],-1,1), domain([ANx,ANy,ANz],0,1), domain([Nx,Ny,Nz,T1,T2],-2,2), domain([I],0,5), %%% avoid forbidden bend angles (0 and 180) 10 * Ax + 5 * Ay + Az #\= 10 * Bx + 5 * By + Bz, 10 * Ax + 5 * Ay + Az #\= (-10 * Bx - 5 * By -Bz), %%% avoid forbidden bend angles (60) scalar_product([Ax,Ay,Az],[Bx,By,Bz],Ang), Ang #>=0, scalar_product([Cx,Cy,Cz],[Bx,By,Bz],Ang1), Ang1 #>=0, Ax*Ay*Az #= 0, Ax*Ay + Ax*Az + Ay*Az #\= 0, Bx*By*Bz #= 0, Bx*By + Bx*Bz + By*Bz #\= 0, Cx*Cy*Cz #= 0, Cx*Cy + Cx*Cz + Cy*Cz #\= 0, vector_product([Bx,By,Bz],[Ax,Ay,Az],[Nx,Ny,Nz]), vector_product([Ax,Ay,Az],[ANx,ANy,ANz],[Rx,Ry,Rz]), scalar_product([Rx,Ry,Rz],[Bx,By,Bz],T1), scalar_product([ANx,ANy,ANz],[Bx,By,Bz],T2), ANx #= 1-Ax*Ax, ANy #= 1-Ay*Ay, ANz #= 1-Az*Az, ANx * ANy #=0, ANx * ANz #=0, ANy * ANz #=0, ANx + ANy + ANz #= 1, Tipo #= -1*T1*T2, type_constraint(I,Tipo,Nx,Ny,Nz,Ax,Ay,Az,Bx,By,Bz,Cx,Cy,Cz). type_constraint(I,Tipo,Nx,Ny,Nz,Ax,Ay,Az,Bx,By,Bz,Cx,Cy,Cz) :- %tipo=0 caso 90 gradi %tipo=1 caso 60 orario rispetto al 90 %tipo=-1 caso 60 antiorario rispetto al 90 % Tipo -1 0 1 % K1 0 0 1 % K2 0 1 1 % K3 0 1 0 domain([K1,K2,K3],0,1), K1 #= (Tipo+1)/2, % K2 #= (Tipo+2)/2, % K3 #= K2-K1, % (I#=0 #/\ Cx#= -1*Ax+Bx - K3*Bx #/\ Cy#= -1*Ay+By - K3*By #/\ Cz#= -1*Az+Bz - K3*Bz ) #\/ (I#=1 #/\ (3-K3)*Cx #= (2-K3)*Nx+ (1+K1)*Bx+ (-2+K2)*Ax #/\ (3-K3)*Cy #= (2-K3)*Ny+ (1+K1)*By+ (-2+K2)*Ay #/\ (3-K3)*Cz #= (2-K3)*Nz+ (1+K1)*Bz+ (-2+K2)*Az ) #\/ (I#=2 #/\ (3-K3)*Cx #= (2-K3)*Nx+ (1-2*K1)*Bx+ (1+K1)*Ax #/\ (3-K3)*Cy #= (2-K3)*Ny+ (1-2*K1)*By+ (1+K1)*Ay #/\ (3-K3)*Cz #= (2-K3)*Nz+ (1-2*K1)*Bz+ (1+K1)*Az ) #\/ (I#=3 #/\ Cx#=Ax #/\ Cy#=Ay #/\ Cz#=Az) #\/ (I#=4 #/\ (-3+5*K2+K1)*Cx #= (2-3*K2-K1)*Nx+ Bx+ (-2+3*K2)*Ax #/\ (-3+5*K2+K1)*Cy #= (2-3*K2-K1)*Ny+ By+ (-2+3*K2)*Ay #/\ (-3+5*K2+K1)*Cz #= (2-3*K2-K1)*Nz+ Bz+ (-2+3*K2)*Az ) #\/ (I#=5 #/\ (3-1*K2-5*K1)*Cx #= (-2+1*K2+3*K1)*Nx+ (2-1*K2-2*K1)*Bx+ (-1+3*K1)*Ax #/\ (3-1*K2-5*K1)*Cy #= (-2+1*K2+3*K1)*Ny+ (2-1*K2-2*K1)*By+ (-1+3*K1)*Ay #/\ (3-1*K2-5*K1)*Cz #= (-2+1*K2+3*K1)*Nz+ (2-1*K2-2*K1)*Bz+ (-1+3*K1)*Az ). % K1 0 0 1 % K2 0 1 1 % K3 0 1 0 %%% PREDICATE: energy_constraints %%% INPUT: The Primary list of aminoacids, %%% The Tertiary list of triples, %%% OUTPUT: Computes in Energy (as a constraint!) %%% the value of energy associated to a spatial %%% conformation and the list of potentials energy_constraints(_Primary, _,[], _, 0, _, _). energy_constraints(_Primary, _,[_,_,_], _, 0, _, _). energy_constraints(_Primary, _,[_,_,_,_,_,_], _, 0, _, _). energy_constraints(Primary, _Secondary, Tertiary, Indexes, Energy, Potential_list,Matrix,Matrixdist):- indexes2segments(Indexes,Segm,1), %%creates the full matrix of pots length(Primary,N), Matrixlen is N*N, length(Matrix,Matrixlen), length(Matrixdist,Matrixlen), % constrain_secondary_contact(Matrix,N,7,18,23,30,Tertiary), constrain_matrix(1,Primary, Segm, Tertiary, Matrix, Matrixdist, N), print_matrix(Matrix), create_potlist(1,Matrix,N,Potential_list), sum(Potential_list,#=,Energy). constrain_secondary_contact(_Matrix,_N,A,B,C,D,T):- auxfor(C,A,B,C,D,T). auxfor(N,A,B,C,D,T):- N= I+2, J =< I+3,!, matrixel(Matrix,Size,I,J,0) . constrain_mindist(_Matrix,_Size,_I,_J). constrain_contact(Matrix, N, I, J, _A1, _A2, _, _, Segm):- included(I,J,Segm),!, matrixel(Matrix,N,I,J,0). constrain_contact(Matrix,N,I,J,A1,A2,[X1,Y1,Z1],[X2,Y2,Z2],_Segm):- %%% skip contacts for aa with distance of 2 in the sequence J >= I+2,!, matrixel(Matrix,N,I,J,El), table(A1, A2, Pot), El in {0,Pot}, DX #= abs(X1 - X2), DY #= abs(Y1 - Y2), DZ #= abs(Z1 - Z2), 2 #= DX + DY + DZ #<=> El #= Pot. %%% CONTACT CONTRIBUTION constrain_contact(Matrix, N, I, J, _A1, _A2, _, _, _):- matrixel(Matrix,N,I,J,0). %%%converts runs of ground indexes in intervals for matrix indexes2segments(Indexes,R,P):- vars(Indexes, Nv, I1 ), Nv>0,!, P1 is P+Nv, indexes2segments(I1,R,P1). indexes2segments(Indexes,[A,B|R1],P):- grounds(Indexes, Ng, I1), Ng >0,!, A is P+1, B is P+Ng+1, P1 is P+Ng, indexes2segments(I1,R1,P1). indexes2segments([],[],_). %%returns the nth amino terelem(Tertiary,A1,[X,Y,Z]):- Ind11a is (A1-1)*3+1, Ind12a is (A1-1)*3+2, Ind13a is (A1-1)*3+3, nth(Ind11a,Tertiary,X), nth(Ind12a,Tertiary,Y), nth(Ind13a,Tertiary,Z). matrixel(Matrix,N,I,J,El):- J > I, I >0, I < N+1, J > 0, J < N+1,!, Ind is J+(I-1)*N, nth(Ind, Matrix, El). %%%if I,J not in bounds, return a free variable matrixel(_Matrix,_N,_I,_J,_El). sortindex(I,J,I,J):- I =< J,!. sortindex(I,J,J,I). %%% selects from the matrix, only the elems not grounds (== 0) create_potlist(N, Matrix,Size,PotList):- N =< Size,!, M is N+1, create_potlist(N, M, Matrix, Size, List), create_potlist(M, Matrix, Size, List1), append(List,List1,PotList). create_potlist(_N, _Matrix, _Size, []). create_potlist(N, M, Matrix, Size, List):- M =< Size, matrixel(Matrix,Size,N,M,El), ground(El),!, %% El is ground, so is not considered M1 is M+1, create_potlist(N, M1, Matrix, Size, List). create_potlist(N, M, Matrix, Size, [El|List]):- M =< Size,!, matrixel(Matrix,Size,N,M,El), var(El),!, M1 is M+1, create_potlist(N, M1, Matrix, Size, List). create_potlist(_N, _M, _Matrix, _Size, []). %%% selects from the matrix, only the elems not grounds (== 0) create_subpotlist(N, Matrix,Size,PotList,Y):- N =< Y,!, M is N+1, create_subpotlist(N, M, Matrix, Size, List,Y), create_subpotlist(M, Matrix, Size, List1,Y), append(List,List1,PotList). create_subpotlist(_N, _Matrix, _Size, [],_). create_subpotlist(N, M, Matrix, Size, List,Y):- M =< Y, matrixel(Matrix,Size,N,M,El), ground(El),!, M1 is M+1, create_subpotlist(N, M1, Matrix, Size, List,Y). create_subpotlist(N, M, Matrix, Size, [El|List],Y):- M =< Y,!, matrixel(Matrix,Size,N,M,El), var(El),!, M1 is M+1, create_subpotlist(N, M1, Matrix, Size, List,Y). create_subpotlist(_N, _M, _Matrix, _Size, [],_). included(A,B,[H1,H2|_R]):- H1 =< A, B =< H2,!. included(A,B,[_H1,_H2|R]):- included(A,B,R). distances([T1,T2,T3|R],D):- distances([T1,T2,T3],R,D). distances([P1,P2,P3],[X1,X2,X3|R],[D1,D2,D3|R1]):- vector_sum([P1,P2,P3], [D1,D2,D3], [X1,X2,X3]), distances([P1,P2,P3],R,R1). distances(_,[],[]). %%% considers all contacts for tertiary potential_list([],[],[]). potential_list([_],[_,_,_],[]). potential_list( [A1,A2|ListaAmin], [P1x,P1y,P1z,P2x,P2y,P2z|ListaPos], Lnew):- potential_list([A2|ListaAmin], [P2x,P2y,P2z|ListaPos], L), list_pot_amino(A1,[P1x,P1y,P1z],ListaAmin,ListaPos,Lnew,L). list_pot_amino(Amino,[X1,Y1,Z1],[Amino2|ListaAmin],[X2,Y2,Z2|ListaPos],[C|L],L1):- table(Amino, Amino2, Pot), C in {0,Pot}, DX #= abs(X1 - X2), DY #= abs(Y1 - Y2), DZ #= abs(Z1 - Z2), 2 #= DX + DY + DZ #<=> C #= Pot, %%% CONTACT CONTRIBUTION list_pot_amino(Amino,[X1,Y1,Z1],ListaAmin,ListaPos,L,L1). list_pot_amino(_Amino,_Posizione,[],[],L,L). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% END OF CONSTRAINT DEFINITION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% LABELING PHASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% PREDICATE: local_labeling %%% INPUT: The list Primary of aminoacids, The list of Secondary %%% structure Information, The list of Indexes, %%% The flat list Tertiary_flat of coordinates %%% the Matrix of contacts %%% OUTPUT: The ground values for Tertiary variables! %%% NOTE: If Indexes are all known, perform an immediate standard labeling. %%% Otherwise, look at Indexes, which is a sequence of alternations: %%% [(V1),G1,V2,G2,...] of variables and ground values. %%% Find the "best" sublist [G1,V2,G2] (or [V1,G1], or [Gn,Vn]) %%% The best is that having |G1|+|V2|+|G2| - 2 * |V2| bigger %%% (i.e. less variables than constant part). %%% Perform the labeling on that sublist!!!! local_labeling(_Primary, _Secondary, Indexes, Tertiary, Energy,_M):- write('Welcome to local labeling end? '),write(Indexes),nl, ground( Indexes ), ground( Tertiary), ! , nl,write('Got all indexes labelled!'), write(' Index : '),write(Indexes),nl, write(' Ter : '),write(Tertiary),nl, write(' Energy : '),write(Energy),nl . local_labeling(_Primary, _Secondary, Indexes, Tertiary,Energy,_M):- write('Welcome to local only labeling? '),write(Indexes),nl, ground( Indexes ), ! , append(Indexes,Tertiary,AA), labeling([ff], AA), %%% Base case: Built-in labeling write(' Index : '),write(Indexes),nl, write(' Ter : '),write(Tertiary),nl, write(' Energy : '),write(Energy),nl . local_labeling(Primary, Secondary, Indexes, Tertiary, Energy,M):- nl,write('Welcome to local labeling for'), write(' Index : '),write(Indexes),nl, %% if subsequencing choose_subindexes(Indexes, S, L, Secondary), %% if no subsequencing % S=1, length(Indexes,L), %%% Debugging: write('local labeling: starting point and substring length:'), write(S), write(' , '), write(L),nl, %%% L_prim is 3+L, sublist(S,L,Indexes,SubIndexes), sublist(S,L_prim,Primary,SubPri), St is (S-1)*3+1, Lt is (L_prim)*3, sublist(St,Lt,Tertiary,SubTertiary), length(M,LM), Size is integer(sqrt(LM)), Y is S + L_prim-1, create_subpotlist(S, M,Size,ListaPot,Y), sum(ListaPot, #=, SubEnergy), write(' Index: '),write(SubIndexes),nl, write(' ter : '),write(Tertiary),nl, write(' Subter : '),write(SubTertiary),nl, (choose_labeling( SubIndexes, SubEnergy, ListaPot, SubTertiary,SubPri,M); !,fail), write(' Listap: '),write(ListaPot),nl, write(' Energy: '),write(SubEnergy),nl, write(' ter : '),write(SubTertiary),nl, print_matrix(M), local_labeling(Primary, Secondary, Indexes, Tertiary, Energy,M). global_tertiary(Primary,Tertiary):- length(Primary,N), N3 is N*3, length(Tertiary,N3). %%% PREDICATE: choose_subindexes (aux. admissible_sub, best_sub) %%% INPUT: The list of Indexes, %%% OUTPUT: Two indexes Start and Length (starting point and %%% length of the sum of variables and non variables (2 alternations) %%% at the beginning choose_subindexes(Indexes, Start, Length, Secondary):- admissible_sub(Indexes, 1, ListaPoss), write('Possibilities: '),write(ListaPoss),nl, best_sub(ListaPoss, Start, Length, Secondary,Indexes). %%% PREDICATE: choose_labeling %%% INPUT: The Sublist of the Indexes and the Energy, %%% OUTPUT: The List of partial potentials, and the associated Energy %%% NOTE: If the number of vars is little, use built-in labeling. %%% Otherwise, use labeling_new (with Heuristics). choose_labeling( SubIndexes, Energy, _ , SubTertiary,_SubPrimary,_M):- %%% (1): built-in num_var(SubIndexes, Nvars), length(SubIndexes, Sublength), %%% Test: (values can be changed) Nvars < 4, Sublength < 20,!, writetime('Choosing built-in labeling...'),nl, append(SubIndexes,SubTertiary,AA), ( labeling([ff,minimize(Energy)], AA);!,fail), %%% Debugging writetime('### built-in labeling time and value: '), write(Energy),nl, write(SubIndexes),nl. choose_labeling( SubIndexes, Energy, ListaPot, SubTertiary,SubPrimary, M):- %%% (2): ad hoc (with heuristics) write('Choosing ad hoc labeling...'),nl, % generate 0 matrix length(ListaPot,P), length(M0,P), domain(M0,0,0), assert( best_solution(10000) ), clean_bestlistapot, assert( best_listapot(M0) ), append(SubIndexes,SubTertiary,AA), length(AA,To_do), print_matrix(M), minimum_domain(ff, Min, AA, OtherVars), %% ff or leftmost or min minimize( (labeling_new(To_do, 0, 1, [Min|OtherVars], Energy, ListaPot, M,SubIndexes,SubTertiary) %%% Debugging ,writetime('### ad-hoc labeling time and value: '), nl,print_results(SubPrimary,SubTertiary), write(ListaPot),nl, write(SubIndexes),nl, write(Energy),nl ), Energy ), retract(best_solution(_)), retract(best_listapot(_)). %%% PREDICATE: labeling_new %%% INPUT: In the call to minimum_domain %%% "ff" first-fail, "leftmost" leftmost strategy %%% To be better explained %%% A first-fail or leftmost, self-programmed strategy: %%% (5 is the number of variables to be instantiated before %%% the next heuristic control) labeling_new( To_do, Done, MOD, [ A | R], Energy, ListaPot, M,I,Ter) :- Done < To_do, ground(A), Now_Done is Done + 1, labeling_new(To_do, Now_Done, MOD, R, Energy, ListaPot, M,I,Ter). labeling_new( To_do, Done, MOD, [ A | R], Energy, ListaPot, M,I,Ter) :- Done < To_do, MOD < 5, var(A), Now_Done is Done + 1, MOD1 is MOD + 1, minimum_domain(ff, Min, [ A | R], OtherVars), %% ff or leftmost indomain(Min), labeling_new(To_do, Now_Done, MOD1, OtherVars, Energy, ListaPot, M,I,Ter). labeling_new( To_do, Done, 5, [ A | R], Energy, ListaPot, M,I,Ter) :- Done < To_do, var(A), best_listapot(BestPL), minimum_domain(ff, Min, [ A | R], OtherVars), %% ff or leftmost indomain(Min), term_variables(ListaPot,LpV), length(LpV,LVarPot), length(ListaPot,LTotPot), sumgrounds( ListaPot, BestPL, Valint, ValintBest), Now_Done is Done + 1, %%% beginning end CoefS1 is (0.5 * (LVarPot) + 1.1 * (LTotPot-LVarPot)) / (LTotPot), S1 is integer( CoefS1 * ValintBest)+1, Valint < S1, write(Valint), write(' '), write(S1), write(' '), write(ValintBest), write(' '),write( LVarPot),writetime(' '), write(I),nl, labeling_new(To_do, Now_Done, 1, OtherVars, Energy, ListaPot, M,I,Ter). labeling_new( Done, Done, _, [] , Energy, ListaPot, _M, I,Ter) :- best_solution( Val ), write(Ter),nl, write(I),nl, ( Energy < Val, !, retract( best_solution(Val)), assert( best_solution(Energy)), retract( best_listapot(_)), assert( best_listapot(ListaPot)) ; true). %%%% PREDICATE: minimum_domain(Options,A, VarList, OrderedList). %%%% (aux. all_domain_size,remove_key; built.in: keysort,fd_size) %%%% INPUT: The list VarList of finite domain variables %%%% OUTPUT: One of the variable (A) with smaller domain, and %%%% the list of the other variables, ordered %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% minimum_domain(ff,A,VarList,SortedVarList) :- all_domain_size(VarList,DomainSizeList), keysort(DomainSizeList,DomainSizeListSorted), remove_key([A|SortedVarList],DomainSizeListSorted). minimum_domain(leftmost,A,[A|R], R). %%% Note: keysort sorts lists of pairs Key-Value all_domain_size([],[]). all_domain_size([A|R],[ASize-A|S]) :- fd_size(A,ASize), all_domain_size(R,S). remove_key([],[]). remove_key([A|R],[_-A|S]) :- remove_key(R,S). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% AUXILIARY PREDICATES FOR THE LABELING PHASE %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% PREDICATE: local_ssbond %%% INPUT: The list of Secondary Structure, The starting point Sn, the length Ln %%% OUTPUT: A new list of SSBONDS, parametric w.r.t. the sublist involved %%% NOTE: Currently, it works only if the ssbonds are put first in data.pl!!! local_ssbond([ssbond(A,B)|Rest], Sn, Ln, [ssbond(C,D)|Sec]):- %%% assumed that A < B A >= Sn + 1, B =< Sn + Ln +1, !, C is A - Sn + 1, D is B - Sn + 1, local_ssbond(Rest, Sn, Ln, Sec). local_ssbond([ssbond(_,_) | Rest], Sn, Ln, Sec):- !, local_ssbond(Rest, Sn, Ln, Sec). local_ssbond(_, _, _, []). local_secondary([ssbond(A,B)|Rest], Sn, Ln, [ssbond(C,D)|Sec]):- A >= Sn + 1, B =< Sn + Ln +1, !, C is A - Sn + 1, D is B - Sn + 1, local_secondary(Rest, Sn, Ln, Sec). local_secondary([ssbond(_,_) | Rest], Sn, Ln, Sec):- !, local_secondary(Rest, Sn, Ln, Sec). local_secondary([strand(A,B)|Rest], Sn, Ln, [strand(C,D)|Sec]):- C is A - Sn + 1, D is B - Sn + 1, local_secondary(Rest, Sn, Ln, Sec). local_secondary([helix(A,B)|Rest], Sn, Ln, [helix(C,D)|Sec]):- C is A - Sn + 1, D is B - Sn + 1, local_secondary(Rest, Sn, Ln, Sec). local_secondary(_, _, _, []). %%% PREDICATE: admissible_sub %%% INPUT: The list of Indexes, a number i %%% OUTPUT: A list of triples: [Starting index, %%% #vars in the sublist + #non vars in the sublist, %%% #vars in the sublist ] %%% Note: very compact code to avoid multiple calls of %%% auxiliary procedure %%% Nuova quintupla: con una struttura [c,X,d,y,e] fornisco le %%% cinque lunghezze. admissible_sub(Indexes, 1, [[1, [0,Nvars1,L2,Nvars2,L3,Nvars3], Length, Nvars]|ListaPoss]):- vars(Indexes, Nvars1, R1 ), Nvars1 > 0, !, grounds(R1, L2, R2), vars(R2, Nvars2, R3 ), grounds(R3, L3, R4 ), vars(R4, Nvars3, _R5 ), (%%Case 1: [ Vars ] R1 = [], !, ListaPoss = [], Length = Nvars1, Nvars = Nvars1, L2 = 0; %%Case 2: [ Vars, Grounds | ...] R2=[],!, Nvars=Nvars1, Length is L2 + Nvars1, NextStart is Nvars1 + 1, ListaPoss = []; %%Case 2/bis: [ Vars, Grounds, Vars| ...] R3=[],!, Nvars is Nvars1+Nvars2, Length is L2 + Nvars1 + Nvars2, NextStart is Nvars1 + 1, admissible_sub(R1, NextStart, ListaPoss); %%Case 2/ter: [ Vars, Grounds, Vars, Grounds| ...] R4=[],!, Nvars is Nvars1+Nvars2, Length is L2 + Nvars1 + Nvars2+L3, NextStart is Nvars1 + 1, admissible_sub(R1, NextStart, ListaPoss); %%Case 2/q: [ Vars, Grounds, Vars, Grounds| ...] Nvars is Nvars1+Nvars2+Nvars3, Length is L2 + Nvars1 + Nvars2+L3+Nvars3, NextStart is Nvars1 + 1, admissible_sub(R1, NextStart, ListaPoss) ). admissible_sub(Indexes, Start, [[Start,[L1,Nvars1,L2,Nvars2,L3,Nvars3],Length, Nvars]| ListaPoss]):- grounds(Indexes, L1, R1), vars(R1,Nvars1,R2), Nvars1 > 0, !, grounds(R2,L2,R3), vars(R3,Nvars2,R4), grounds(R4,L3,R5), vars(R5,Nvars3,R6), (%%% Case 4: [ Grounds , Vars ] R2 = [], !, ListaPoss = [], Nvars = Nvars1, Length is Nvars + L1; %%% Case 5: [ Ground , Vars , Ground] R3 = [], !, ListaPoss = [], Nvars = Nvars1, Length is Nvars1 + L1 + L2; %%% Case6: [ Grounds , Vars , Grounds , Vars ] R4 = [], !, % ListaPoss = [], Nvars is Nvars1 + Nvars2, Length is L1 + Nvars1 + L2 + Nvars2; %%% Case7: [ Grounds , Vars , Ground , Vars , Ground ] R5 = [], !, ListaPoss = [], Length is L1 + Nvars1 + L2 + Nvars2 + L3, Nvars is Nvars1 + Nvars2; %%%% Case8: [ Grounds , Vars , Ground , Vars , Ground | ...] R6 = [], !, Length is L1 + Nvars1 + L2 + Nvars2 + L3, Nvars is Nvars1 + Nvars2, NewStart is Start + L1 + Nvars1, admissible_sub(R2,NewStart,ListaPoss); %%%% Case9: [ Grounds , Vars , Ground , Vars , Ground, V | ...] Length is L1 + Nvars1 + L2 + Nvars2 + L3 + Nvars3, Nvars is Nvars1 + Nvars2 + Nvars3, NewStart is Start + L1 + Nvars1, admissible_sub(R2,NewStart,ListaPoss)). %%% Case 3: Indexes = [ Ground ] admissible_sub(_, _, []). %%%%%%%%%%% NOTE:further cases can be analyzed to optimize!!!! %%% PREDICATE: best_sub/3 (aux. best_subaux/3) %%% INPUT: A list of triples: [Starting index, L , Nvar ] %%% identifying sublists [G1,V1,G2,V2,G3] where %%% L = |G1|+|V1|+|G2|+ |V2| + |G3|, Nvar = |V1| + |V2|. %%% generated by admissible_sub %%% OUTPUT: The best tuple (the last value is changed, %%% for optimization) %%% NOTE: We maximize L - 2 * Nvar %%% We use best_sub/3 for optimizing %%% If more than 10 variables are in the selected sublist, %%% we cut it, using possible_cut best_sub(List, Start, Length, Secondary ,I) :- best_subaux(List, [_,_,_,-1000,_],Tuple, Secondary,I), write(Tuple),nl, possible_cut(Tuple, Start, Length, Secondary). %%%if included a ssbond the score is increased by 100 best_subaux([ [Start,[C1,V1,C2,V2,C3,V3],L,Nvar] | R], [_,_,_,Pbest,_],Triple, Secondary,I) :- check_include_ssbond(Start,L,SS, Secondary,I,A,B), %%% if less than 10 variables, set NV=10 NV #=10 #<=> Nvar #=< 10, NV #=Nvar #<=> Nvar#> 10, P1 is L - 2*NV+100*SS, write(Start),write(' '),write(L),write(' '),write(SS),write(' '),write(P1),write(' '),write(Pbest),nl, P1 > Pbest,!, best_subaux( R , [Start, [C1,V1,C2,V2,C3,V3], L, P1,[A,B]] ,Triple, Secondary,I). best_subaux([ [_,_,_,_] | R],[StartBest, Qbest,Lbest,NvarBest,Bond],Triple, Secondary,I) :- best_subaux( R , [StartBest, Qbest,Lbest,NvarBest,Bond] ,Triple, Secondary,I). best_subaux([] , BestTriple,BestTriple, _Secondary,_). possible_cut([Start,[C1,V1,C2,V2,C3,V3],Len,_,[A,B]], S, L, _Secondary) :- V1+V2+V3 =<12, A>0, B>0,!, Ai is A-1, Bi is B-2, %begin ( Ai= 1,!, SS is 1. check_include_ssbond(Start,L,SS, [ssbond(_A,_B)|R],Ind,A,B):- check_include_ssbond(Start,L,SS, R,Ind,A,B). check_include_ssbond(_Start,_L,0, _,_,-1,-1). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% AUXILIARY (TECHNICAL AND EASY) PREDICATES %%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% num_var(t) returns the number of non ground %%% elements of the list t %%% uses the built-in term_variables (returns the list of variables) num_var(L,N):- term_variables(L,LV), length(LV,N). %%% vars([A,B,C,a,b], X , Y ) %%% returns X = 3 (number of vars at the beginning) and %%% Y = [a,b] (the list but the first initial vars) vars([],0,[]). vars([X|Y], 0, [X|Y]):- ground(X). vars([X|Y], L, R):- var(X), vars(Y, M, R), L is M+1. %%% grounds([a,b,c,A,B,C], X , Y ) %%% returns X = 3 (number of non vars at the beginning) and %%% Y = [A,B,C] (the list bust the first initial non vars) grounds([],0,[]). grounds([X|Y], 0, [X|Y]):- var(X). grounds([X|Y], L, R):- ground(X), grounds(Y, M, R), L is M+1. %%% sumgrounds([2,X,3,C,5,1],C ) %%% returns C = 11 (the sum of the non ground values) sumgrounds( [ A | R ] , [_B|Best], S, S1 ) :- var(A), !, sumgrounds( R , Best, S, S1 ). sumgrounds( [ A | R ] , [B|Best], S, S1 ) :- sumgrounds( R , Best, Si, Si1 ), S is A + Si, S1 is B + Si1. sumgrounds( [], [], 0, 0). sumgroundster( [ A,B,C | R ] , [Rx,Ry,Rz], T ) :- sumgroundster( [ A,B,C | R ] , [Rx,Ry,Rz], T1 ,N ), T is T1 /N. sumgroundster( [ A,B,C | R ] , [Rx,Ry,Rz], T , N) :- ground(A), ground(B), ground(C), !, sumgroundster(R,[Rx,Ry,Rz], T1 ,N1), N is N1+1, T is T1 + abs(Rx-A)+abs(Ry-B)+abs(Rz-C). sumgroundster( [ _ , _, _ | R ] , Ref, T ,N ) :- sumgroundster( R , Ref, T ,N). sumgroundster( [], _ , 0,0). sumgrounds( [ A | R ] , C ) :- var(A), !, sumgrounds( R , C ). sumgrounds( [ A | R ] , C ) :- sumgrounds( R , C1 ), C is A + C1. sumgrounds( [], 0). %%% sublist(2,3,[a,b,c,d,e,f],X) %%% returns X = [b,c,d] (the list of 3 elements starting from 2) sublist(_Inizio, 0,[],[]). sublist(1, N, List, Sublist):- length( Sublist, N), prefix( Sublist, List). sublist(Inizio, N, [_Testa|Coda], Sublist):- Inizio > 1, NuovoInizio is Inizio-1, sublist( NuovoInizio, N, Coda, Sublist). %%% even/odd predicate: even(X):- 0 #= X mod 2. odd(X):- 1 #= X mod 2. %%% VECTORS (3D) operations: vector_sum([A,B,C],[D,E,F],[G,H,I]):- G #= A + D, H #= B + E, I #= C + F. vector_product([Ax,Ay,Az],[Bx,By,Bz],[Cx,Cy,Cz]):- Cx #= Ay*Bz - Az*By, Cy #= Az*Bx - Ax*Bz, Cz #= Ax*By - Ay*Bx. scalar_product([Ax,Ay,Az],[Bx,By,Bz],S):- S #= Ax*Bx + Ay*By+ Az*Bz. vector_minus([X,Y,Z],[X1,Y1,Z1]):- X1 #= 0 - X, Y1 #= 0 - Y, Z1 #= 0 - Z. extract([],[]). extract([A,B,C|R],[A,B,C|R1]):- extract(R,R1). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% The list of the 20 aminoacids %%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% aminoacid(a). aminoacid(c). aminoacid(d). aminoacid(e). aminoacid(f). aminoacid(g). aminoacid(h). aminoacid(i). aminoacid(k). aminoacid(l). aminoacid(m). aminoacid(n). aminoacid(p). aminoacid(q). aminoacid(r). aminoacid(s). aminoacid(t). aminoacid(v). aminoacid(w). aminoacid(y). %%%%%% Print predicate print_results( Primary, Tertiary ) :- print_results( 1 , Primary, Tertiary ). print_results( I , [Amino|R], [X,Y,Z|N]) :- write(I), write(' '), write(Amino), write(' '), write(X), write(' '), write(Y), write(' '), write(Z), write(' '), nl, I1 is I + 1, print_results( I1, R, N). print_results(_ , [], []). %%%%%%%%% Timing Predicates: initial_time :- statistics(runtime,[T|_]), assert(start(T)). writetime(Str) :- start(InitialTime), statistics(runtime,[CurrentTime|_]), Delta is CurrentTime - InitialTime, write(Str), print_time(Delta), write('\t'). print_time( K ) :- K < 60000,!, S is integer(K/1000), M is K mod 1000, write(S), write(','), write(M), write('s'). print_time( K ) :- K < 3600000,!, M is integer(K / 60000), S is K mod 60000, write(M), write('min.'), print_time(S). print_time( K ) :- H is integer(K / 3600000), M is K mod 3600000, write(H), write('hours.'), print_time(M). print_matrix(M,I):- length(M,LM), N is integer(sqrt(LM)), print_matrix(M,N,N,[*,*,*|I]). print_matrixt(M,I,T):- length(M,LM), N is integer(sqrt(LM)), print_matrix(M,N,N,[*,*,*|I],T). print_matrix(M):- length(M,LM), N is integer(sqrt(LM)), print_matrix(M,N,N). print_matrix(_Matrix,0,_N). print_matrix(Matrix,M,N):- A is (N-M)*N+1, Chi is (N-M)+1, sublist(A,N,Matrix,List), write(Chi),write(' '), print_line_matrix(List),nl, M1 is M-1, print_matrix(Matrix,M1,N). print_matrix(_Matrix,0,_N,_). print_matrix(Matrix,M,N,[Ii|I]):- A is (N-M)*N+1, Chi is (N-M)+1, sublist(A,N,Matrix,List), write(Chi),write(' '), print_line_matrix(List), write('|'),write(Ii),nl, M1 is M-1, print_matrix(Matrix,M1,N,I). print_matrix(_Matrix,0,_N,_,_). print_matrix(Matrix,M,N,[Ii|I],[T1,T2,T3|Tr]):- A is (N-M)*N+1, Chi is (N-M)+1, sublist(A,N,Matrix,List), write(Chi),write(' '), print_line_matrix(List), write('|'),write(Ii), write(' |'),write(T1), write(' '),write(T2), write(' '),write(T3),nl, M1 is M-1, print_matrix(Matrix,M1,N,I,Tr). print_line_matrix([]). print_line_matrix([A|R]):- var(A),!, write('.'),write(','), print_line_matrix(R). print_line_matrix([A|R]):- A=0,!, write(A),write(','), print_line_matrix(R). print_line_matrix([_|R]):- write('*,'), print_line_matrix(R). print_matrixd(M):- length(M,LM), N is integer(sqrt(LM)), print_matrixd(M,N,N). print_matrixd(_Matrix,0,_N). print_matrixd(Matrix,M,N):- A is (N-M)*N+1, Chi is (N-M)+1, sublist(A,N,Matrix,List), write(Chi),write(' '), print_line_matrixd(List),nl, M1 is M-1, print_matrixd(Matrix,M1,N). print_line_matrixd([]). print_line_matrixd([A|R]):- fd_max(A,AA), fd_min(A,AA1), fd_size(A,AA2), write(AA),write('-'), write(AA1),write('!'), write(AA2),write(','), print_line_matrixd(R). %%%%%%%% %%% For removing asserted facts, %%% remained in memory due to execution abort clean :- retract(start(_)), clean. clean :- retract(best_solution(_)), clean. clean :- retract(best_listapot(_)), clean. clean. clean_bestlistapot :- retract(best_listapot(_)), clean_bestlistapot. clean_bestlistapot. cleanmatrix :- retract(matrix(_)). cleanmatrix. %%% Execute: :-clean. :-cleanmatrix. printtable:- table(A,B,C), write('table('),write(A),write(' ,'),write(B),write(', '),write(C),write(').'),nl, fail.