Merge remote-tracking branch 'origin/tasks/JAL-3035_remove_dasobert_dependency' into...
[jalview.git] / src / jalview / analysis / TreeBuilder.java
1 /*
2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
4  * 
5  * This file is part of Jalview.
6  * 
7  * Jalview is free software: you can redistribute it and/or
8  * modify it under the terms of the GNU General Public License 
9  * as published by the Free Software Foundation, either version 3
10  * of the License, or (at your option) any later version.
11  *  
12  * Jalview is distributed in the hope that it will be useful, but 
13  * WITHOUT ANY WARRANTY; without even the implied warranty 
14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
15  * PURPOSE.  See the GNU General Public License for more details.
16  * 
17  * You should have received a copy of the GNU General Public License
18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
19  * The Jalview Authors are detailed in the 'AUTHORS' file.
20  */
21 package jalview.analysis;
22
23 import jalview.api.analysis.ScoreModelI;
24 import jalview.api.analysis.SimilarityParamsI;
25 import jalview.datamodel.AlignmentView;
26 import jalview.datamodel.CigarArray;
27 import jalview.datamodel.SeqCigar;
28 import jalview.datamodel.SequenceI;
29 import jalview.datamodel.SequenceNode;
30 import jalview.math.MatrixI;
31 import jalview.viewmodel.AlignmentViewport;
32
33 import java.util.BitSet;
34 import java.util.Vector;
35
36 public abstract class TreeBuilder
37 {
38   public static final String AVERAGE_DISTANCE = "AV";
39
40   public static final String NEIGHBOUR_JOINING = "NJ";
41
42   protected Vector<BitSet> clusters;
43
44   protected SequenceI[] sequences;
45
46   public AlignmentView seqData;
47
48   protected BitSet done;
49
50   protected int noseqs;
51
52   int noClus;
53
54   protected MatrixI distances;
55
56   protected int mini;
57
58   protected int minj;
59
60   protected double ri;
61
62   protected double rj;
63
64   SequenceNode maxdist;
65
66   SequenceNode top;
67
68   double maxDistValue;
69
70   double maxheight;
71
72   int ycount;
73
74   Vector<SequenceNode> node;
75
76   private AlignmentView seqStrings;
77
78   /**
79    * Constructor
80    * 
81    * @param av
82    * @param sm
83    * @param scoreParameters
84    */
85   public TreeBuilder(AlignmentViewport av, ScoreModelI sm,
86           SimilarityParamsI scoreParameters)
87   {
88     int start, end;
89     boolean selview = av.getSelectionGroup() != null
90             && av.getSelectionGroup().getSize() > 1;
91     seqStrings = av.getAlignmentView(selview);
92     if (!selview)
93     {
94       start = 0;
95       end = av.getAlignment().getWidth();
96       this.sequences = av.getAlignment().getSequencesArray();
97     }
98     else
99     {
100       start = av.getSelectionGroup().getStartRes();
101       end = av.getSelectionGroup().getEndRes() + 1;
102       this.sequences = av.getSelectionGroup()
103               .getSequencesInOrder(av.getAlignment());
104     }
105
106     init(seqStrings, start, end);
107
108     computeTree(sm, scoreParameters);
109   }
110
111   public SequenceI[] getSequences()
112   {
113     return sequences;
114   }
115
116   /**
117    * DOCUMENT ME!
118    * 
119    * @param nd
120    *          DOCUMENT ME!
121    * 
122    * @return DOCUMENT ME!
123    */
124   double findHeight(SequenceNode nd)
125   {
126     if (nd == null)
127     {
128       return maxheight;
129     }
130
131     if ((nd.left() == null) && (nd.right() == null))
132     {
133       nd.height = ((SequenceNode) nd.parent()).height + nd.dist;
134
135       if (nd.height > maxheight)
136       {
137         return nd.height;
138       }
139       else
140       {
141         return maxheight;
142       }
143     }
144     else
145     {
146       if (nd.parent() != null)
147       {
148         nd.height = ((SequenceNode) nd.parent()).height + nd.dist;
149       }
150       else
151       {
152         maxheight = 0;
153         nd.height = (float) 0.0;
154       }
155
156       maxheight = findHeight((SequenceNode) (nd.left()));
157       maxheight = findHeight((SequenceNode) (nd.right()));
158     }
159
160     return maxheight;
161   }
162
163   /**
164    * DOCUMENT ME!
165    * 
166    * @param nd
167    *          DOCUMENT ME!
168    */
169   void reCount(SequenceNode nd)
170   {
171     ycount = 0;
172     // _lycount = 0;
173     // _lylimit = this.node.size();
174     _reCount(nd);
175   }
176
177   /**
178    * DOCUMENT ME!
179    * 
180    * @param nd
181    *          DOCUMENT ME!
182    */
183   void _reCount(SequenceNode nd)
184   {
185     // if (_lycount<_lylimit)
186     // {
187     // System.err.println("Warning: depth of _recount greater than number of
188     // nodes.");
189     // }
190     if (nd == null)
191     {
192       return;
193     }
194     // _lycount++;
195
196     if ((nd.left() != null) && (nd.right() != null))
197     {
198
199       _reCount((SequenceNode) nd.left());
200       _reCount((SequenceNode) nd.right());
201
202       SequenceNode l = (SequenceNode) nd.left();
203       SequenceNode r = (SequenceNode) nd.right();
204
205       nd.count = l.count + r.count;
206       nd.ycount = (l.ycount + r.ycount) / 2;
207     }
208     else
209     {
210       nd.count = 1;
211       nd.ycount = ycount++;
212     }
213     // _lycount--;
214   }
215
216   /**
217    * DOCUMENT ME!
218    * 
219    * @return DOCUMENT ME!
220    */
221   public SequenceNode getTopNode()
222   {
223     return top;
224   }
225
226   /**
227    * 
228    * @return true if tree has real distances
229    */
230   public boolean hasDistances()
231   {
232     return true;
233   }
234
235   /**
236    * 
237    * @return true if tree has real bootstrap values
238    */
239   public boolean hasBootstrap()
240   {
241     return false;
242   }
243
244   public boolean hasRootDistance()
245   {
246     return true;
247   }
248
249   /**
250    * Form clusters by grouping sub-clusters, starting from one sequence per
251    * cluster, and finishing when only two clusters remain
252    */
253   void cluster()
254   {
255     while (noClus > 2)
256     {
257       findMinDistance();
258
259       joinClusters(mini, minj);
260
261       noClus--;
262     }
263
264     int rightChild = done.nextClearBit(0);
265     int leftChild = done.nextClearBit(rightChild + 1);
266
267     joinClusters(leftChild, rightChild);
268     top = (node.elementAt(leftChild));
269
270     reCount(top);
271     findHeight(top);
272     findMaxDist(top);
273   }
274
275   /**
276    * Returns the minimum distance between two clusters, and also sets the
277    * indices of the clusters in fields mini and minj
278    * 
279    * @return
280    */
281   protected abstract double findMinDistance();
282
283   /**
284    * Calculates the tree using the given score model and parameters, and the
285    * configured tree type
286    * <p>
287    * If the score model computes pairwise distance scores, then these are used
288    * directly to derive the tree
289    * <p>
290    * If the score model computes similarity scores, then the range of the scores
291    * is reversed to give a distance measure, and this is used to derive the tree
292    * 
293    * @param sm
294    * @param scoreOptions
295    */
296   protected void computeTree(ScoreModelI sm, SimilarityParamsI scoreOptions)
297   {
298     distances = sm.findDistances(seqData, scoreOptions);
299
300     makeLeaves();
301
302     noClus = clusters.size();
303
304     cluster();
305   }
306
307   /**
308    * Finds the node, at or below the given node, with the maximum distance, and
309    * saves the node and the distance value
310    * 
311    * @param nd
312    */
313   void findMaxDist(SequenceNode nd)
314   {
315     if (nd == null)
316     {
317       return;
318     }
319
320     if ((nd.left() == null) && (nd.right() == null))
321     {
322       double dist = nd.dist;
323
324       if (dist > maxDistValue)
325       {
326         maxdist = nd;
327         maxDistValue = dist;
328       }
329     }
330     else
331     {
332       findMaxDist((SequenceNode) nd.left());
333       findMaxDist((SequenceNode) nd.right());
334     }
335   }
336
337   /**
338    * Calculates and returns r, whatever that is
339    * 
340    * @param i
341    * @param j
342    * 
343    * @return
344    */
345   protected double findr(int i, int j)
346   {
347     double tmp = 1;
348
349     for (int k = 0; k < noseqs; k++)
350     {
351       if ((k != i) && (k != j) && (!done.get(k)))
352       {
353         tmp = tmp + distances.getValue(i, k);
354       }
355     }
356
357     if (noClus > 2)
358     {
359       tmp = tmp / (noClus - 2);
360     }
361
362     return tmp;
363   }
364
365   protected void init(AlignmentView seqView, int start, int end)
366   {
367     this.node = new Vector<SequenceNode>();
368     if (seqView != null)
369     {
370       this.seqData = seqView;
371     }
372     else
373     {
374       SeqCigar[] seqs = new SeqCigar[sequences.length];
375       for (int i = 0; i < sequences.length; i++)
376       {
377         seqs[i] = new SeqCigar(sequences[i], start, end);
378       }
379       CigarArray sdata = new CigarArray(seqs);
380       sdata.addOperation(CigarArray.M, end - start + 1);
381       this.seqData = new AlignmentView(sdata, start);
382     }
383
384     /*
385      * count the non-null sequences
386      */
387     noseqs = 0;
388
389     done = new BitSet();
390
391     for (SequenceI seq : sequences)
392     {
393       if (seq != null)
394       {
395         noseqs++;
396       }
397     }
398   }
399
400   /**
401    * Merges cluster(j) to cluster(i) and recalculates cluster and node distances
402    * 
403    * @param i
404    * @param j
405    */
406   void joinClusters(final int i, final int j)
407   {
408     double dist = distances.getValue(i, j);
409
410     ri = findr(i, j);
411     rj = findr(j, i);
412
413     findClusterDistance(i, j);
414
415     SequenceNode sn = new SequenceNode();
416
417     sn.setLeft((node.elementAt(i)));
418     sn.setRight((node.elementAt(j)));
419
420     SequenceNode tmpi = (node.elementAt(i));
421     SequenceNode tmpj = (node.elementAt(j));
422
423     findNewDistances(tmpi, tmpj, dist);
424
425     tmpi.setParent(sn);
426     tmpj.setParent(sn);
427
428     node.setElementAt(sn, i);
429
430     /*
431      * move the members of cluster(j) to cluster(i)
432      * and mark cluster j as out of the game
433      */
434     clusters.get(i).or(clusters.get(j));
435     clusters.get(j).clear();
436     done.set(j);
437   }
438
439   /*
440    * Computes and stores new distances for nodei and nodej, given the previous
441    * distance between them
442    */
443   protected abstract void findNewDistances(SequenceNode nodei,
444           SequenceNode nodej, double previousDistance);
445
446   /**
447    * Calculates and saves the distance between the combination of cluster(i) and
448    * cluster(j) and all other clusters. The form of the calculation depends on
449    * the tree clustering method being used.
450    * 
451    * @param i
452    * @param j
453    */
454   protected abstract void findClusterDistance(int i, int j);
455
456   /**
457    * Start by making a cluster for each individual sequence
458    */
459   void makeLeaves()
460   {
461     clusters = new Vector<BitSet>();
462
463     for (int i = 0; i < noseqs; i++)
464     {
465       SequenceNode sn = new SequenceNode();
466
467       sn.setElement(sequences[i]);
468       sn.setName(sequences[i].getName());
469       node.addElement(sn);
470       BitSet bs = new BitSet();
471       bs.set(i);
472       clusters.addElement(bs);
473     }
474   }
475
476   public AlignmentView getOriginalData()
477   {
478     return seqStrings;
479   }
480
481 }