JAL-1159 Merging 2.8.1 patches and Anne Menards work from 2012 into development branch
[jalview.git] / src / jalview / analysis / SecStrConsensus.java
1 package jalview.analysis;
2
3 import java.util.ArrayList;
4 import java.util.Hashtable;
5
6
7 public class SecStrConsensus {
8  
9  
10  /**
11   * Internal class to represent a simple base-pair.
12   * @author Yawn
13   * [JBPNote: ^^ is that Anne Menard or Ya(w)nn Ponty, I wonder ! ]
14   */
15  public static class SimpleBP{
16   int bp5;
17   int bp3;
18   
19   public SimpleBP()
20   {
21           
22   }
23   public SimpleBP(int i5, int i3)
24   {
25    bp5=i5;
26    bp3=i3;
27   }
28   public void setBP5(int i5)
29   {
30            bp5=i5;
31   }
32   
33   public void setBP3(int i3)
34   {
35            bp3=i3;
36   }
37   
38   public int getBP5()
39   {
40            return bp5;
41   }
42   
43   public int getBP3()
44   {
45            return bp3;
46   }
47   
48   public String toString()
49   {
50           return "("+bp5+","+bp3+")";
51   }
52   
53  }
54  
55  public static int[] extractConsensus(ArrayList<ArrayList<SimpleBP>> bps)
56  {
57   // We do not currently know the length of the alignment
58   // => Estimate it as the biggest index of a base-pair plus one. 
59   int maxlength = 0;
60   for (ArrayList<SimpleBP> strs : bps)
61   {
62    for (SimpleBP bp : strs)
63    {
64         
65     maxlength = Math.max(1+Math.max(bp.bp5, bp.bp3), maxlength);
66
67    }
68   }
69   // Now we have a good estimate for length, allocate and initialize data
70   // to be fed to the dynamic programming procedure.
71   ArrayList<Hashtable<Integer,Double>> seq = new ArrayList<Hashtable<Integer,Double>>();
72   for (int i=0;i<maxlength;i++)
73   { seq.add(new Hashtable<Integer,Double>()); }
74   for (ArrayList<SimpleBP> strs : bps)
75   {
76    for (SimpleBP bp : strs)
77    {
78     int i = bp.bp5;
79     int j = bp.bp3;
80     Hashtable<Integer,Double> h = seq.get(i);
81     if (!h.containsKey(j))
82     {
83      h.put(j, 0.0);
84     }
85     h.put(j, h.get(j)+1.);
86    }
87   }
88   // At this point, seq contains, at each position i, a hashtable which associates,
89   // to each possible end j, the number of time a base-pair (i,j) occurs in the alignment
90   
91   // We can now run the dynamic programming procedure on this data
92   double[][] mat = fillMatrix(seq); 
93   ArrayList<SimpleBP> res = backtrack(mat,seq);
94   
95   // Convert it to an array, ie finalres[i] = j >= 0 iff a base-pair (i,j) is present 
96   // in the consensus, or -1 otherwise
97   int[] finalres = new int[seq.size()];
98   for (int i=0;i<seq.size();i++)
99   { finalres[i] = -1; }
100   for (SimpleBP bp : res)
101   {
102    finalres[bp.bp5] = bp.bp3;
103    finalres[bp.bp3] = bp.bp5;
104   }
105   
106   return finalres;
107  }
108  
109  private static boolean canBasePair(ArrayList<Hashtable<Integer,Double>> seq, int i, int k)
110  {
111   return seq.get(i).containsKey(k);
112  }
113  
114  // Returns the score of a potential base-pair, ie the number of structures in which it is found.
115  private static double basePairScore(ArrayList<Hashtable<Integer,Double>> seq, int i, int k)
116  {
117   return seq.get(i).get(k); 
118  }
119  
120  
121   private static double[][] fillMatrix(ArrayList<Hashtable<Integer,Double>> seq)
122   {
123   int n = seq.size();
124   double[][] tab = new double[n][n];
125   for(int m=1;m<=n;m++)
126   {
127    for(int i=0;i<n-m+1;i++)
128    {
129     int j = i+m-1;
130     tab[i][j] = 0;
131     if (i<j)
132     { 
133      tab[i][j] = Math.max(tab[i][j], tab[i+1][j]); 
134      for (int k=i+1;k<=j;k++)
135      {
136       if (canBasePair(seq,i,k))
137       {
138        double fact1 = 0;
139        if (k>i+1)
140        {
141         fact1 = tab[i+1][k-1];
142        }
143        double fact2 = 0;
144        if (k<j)
145        {
146         fact2 = tab[k+1][j];
147        }
148        tab[i][j] = Math.max(tab[i][j],basePairScore(seq,i,k)+fact1+fact2);
149       } 
150      }
151     }
152    }   
153   }
154    return tab;
155   } 
156  
157   private static  ArrayList<SimpleBP> backtrack(double[][] tab,ArrayList<Hashtable<Integer,Double>> seq)
158   {
159    return backtrack(tab,seq,0,seq.size()-1);
160   }
161  
162   private static ArrayList<SimpleBP> backtrack(double[][] tab,ArrayList<Hashtable<Integer,Double>> seq, int i, int j)
163   {
164    ArrayList<SimpleBP> result = new ArrayList<SimpleBP>();
165    if (i<j)
166   { 
167     ArrayList<Integer> indices = new ArrayList<Integer>();
168     indices.add(-1);
169    for (int k=i+1;k<=j;k++)
170    {
171     indices.add(k);
172    }
173    for (int k : indices)
174     {
175      if (k==-1)
176      {
177       if (tab[i][j] == tab[i+1][j])
178       {
179        result = backtrack(tab, seq, i+1,j);
180       }
181      }
182      else
183      {
184       if (canBasePair(seq,i,k))
185       {
186        double fact1 = 0;
187        if (k>i+1)
188        {
189         fact1 = tab[i+1][k-1];
190        }
191        double fact2 = 0;
192        if (k<j)
193        {
194         fact2 = tab[k+1][j];
195        }
196        if (tab[i][j]==basePairScore(seq,i,k)+fact1+fact2)
197        { 
198         result = backtrack(tab, seq, i+1,k-1);
199         result.addAll(backtrack(tab, seq, k+1,j));
200         result.add(new SimpleBP(i,k));
201        }
202       }       
203      }     
204    }
205   }
206   else if  (i==j)
207   {
208    
209   }
210   else 
211   {
212    
213   }
214    return result;
215   }
216  
217
218
219 }