src/jalview/analysis/AverageDistanceEngine.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 java.util.ArrayList;
  24 import java.util.BitSet;
  25 import java.util.List;
  26 import java.util.Vector;
  27
  28 import jalview.datamodel.AlignmentAnnotation;
  29 import jalview.datamodel.BinaryNode;
  30 import jalview.datamodel.ContactListI;
  31 import jalview.datamodel.ContactMatrixI;
  32 import jalview.math.Matrix;
  33 import jalview.viewmodel.AlignmentViewport;
  34
  35 /**
  36  * This class implements distance calculations used in constructing a Average
  37  * Distance tree (also known as UPGMA)
  38  */
  39 public class AverageDistanceEngine extends TreeEngine
  40 {
  41   ContactMatrixI cm;
  42
  43   AlignmentViewport av;
  44
  45   AlignmentAnnotation aa;
  46
  47   /**
  48    * compute cosine distance matrix for a given contact matrix and create a
  49    * UPGMA tree
  50    *
  51    * @param cm
  52    */
  53   public AverageDistanceEngine(AlignmentViewport av, AlignmentAnnotation aa,
  54           ContactMatrixI cm)
  55   {
  56     this.av = av;
  57     this.aa = aa;
  58     this.cm = cm;
  59     calculate(cm);
  60
  61   }
  62
  63   // 0 - normalised dot product
  64   // 1 - L1 - ie (abs(v_1-v_2)/dim(v))
  65   // L1 is more rational - since can reason about value of difference,
  66   // normalised dot product might give cleaner clusters, but more difficult to
  67   // understand.
  68
  69   int mode = 1;
  70
  71   public void calculate(ContactMatrixI cm)
  72   {
  73     this.cm = cm;
  74     node = new Vector<BinaryNode>();
  75     clusters = new Vector<BitSet>();
  76     distances = new Matrix(new double[cm.getWidth()][cm.getWidth()]);
  77     noseqs = cm.getWidth();
  78     done = new BitSet();
  79     double moduli[] = new double[cm.getWidth()];
  80     double max;
  81     if (mode == 0)
  82     {
  83       max = 1;
  84     }
  85     else
  86     {
  87       max = cm.getMax() * cm.getMax();
  88     }
  89
  90     for (int i = 0; i < cm.getWidth(); i++)
  91     {
  92       // init the tree engine node for this column
  93       BinaryNode cnode = new BinaryNode();
  94       cnode.setElement(Integer.valueOf(i));
  95       cnode.setName("c" + i);
  96       node.addElement(cnode);
  97       BitSet bs = new BitSet();
  98       bs.set(i);
  99       clusters.addElement(bs);
 100
 101       // compute distance matrix element
 102       ContactListI ith = cm.getContactList(i);
 103
 104       for (int j = 0; j < i; j++)
 105       {
 106         distances.setValue(i, i, 0);
 107         ContactListI jth = cm.getContactList(j);
 108         if (jth==null)
 109         {
 110           break;
 111         }
 112         double prd = 0;
 113         for (int indx = 0; indx < cm.getHeight(); indx++)
 114         {
 115           if (mode == 0)
 116           {
 117             if (j == 0)
 118             {
 119               moduli[i] += ith.getContactAt(indx) * ith.getContactAt(indx);
 120             }
 121             prd += ith.getContactAt(indx) * jth.getContactAt(indx);
 122           }
 123           else
 124           {
 125             prd += Math
 126                     .abs(ith.getContactAt(indx) - jth.getContactAt(indx));
 127           }
 128         }
 129         if (mode == 0)
 130         {
 131           if (j == 0)
 132           {
 133             moduli[i] = Math.sqrt(moduli[i]);
 134           }
 135           prd = (moduli[i] != 0 && moduli[j] != 0)
 136                   ? prd / (moduli[i] * moduli[j])
 137                   : 0;
 138           prd = 1 - prd;
 139         }
 140         else
 141         {
 142           prd /= cm.getHeight();
 143         }
 144         distances.setValue(i, j, prd);
 145         distances.setValue(j, i, prd);
 146       }
 147     }
 148
 149     noClus = clusters.size();
 150     cluster();
 151   }
 152
 153   /**
 154    * Calculates and saves the distance between the combination of cluster(i) and
 155    * cluster(j) and all other clusters. An average of the distances from
 156    * cluster(i) and cluster(j) is calculated, weighted by the sizes of each
 157    * cluster.
 158    *
 159    * @param i
 160    * @param j
 161    */
 162   @Override
 163   protected void findClusterDistance(int i, int j)
 164   {
 165     int noi = clusters.elementAt(i).cardinality();
 166     int noj = clusters.elementAt(j).cardinality();
 167
 168     // New distances from cluster i to others
 169     double[] newdist = new double[noseqs];
 170
 171     for (int l = 0; l < noseqs; l++)
 172     {
 173       if ((l != i) && (l != j))
 174       {
 175         newdist[l] = ((distances.getValue(i, l) * noi)
 176                 + (distances.getValue(j, l) * noj)) / (noi + noj);
 177       }
 178       else
 179       {
 180         newdist[l] = 0;
 181       }
 182     }
 183
 184     for (int ii = 0; ii < noseqs; ii++)
 185     {
 186       distances.setValue(i, ii, newdist[ii]);
 187       distances.setValue(ii, i, newdist[ii]);
 188     }
 189   }
 190
 191   /**
 192    * {@inheritDoc}
 193    */
 194   @Override
 195   protected double findMinDistance()
 196   {
 197     double min = Double.MAX_VALUE;
 198
 199     for (int i = 0; i < (noseqs - 1); i++)
 200     {
 201       for (int j = i + 1; j < noseqs; j++)
 202       {
 203         if (!done.get(i) && !done.get(j))
 204         {
 205           if (distances.getValue(i, j) < min)
 206           {
 207             mini = i;
 208             minj = j;
 209
 210             min = distances.getValue(i, j);
 211           }
 212         }
 213       }
 214     }
 215     return min;
 216   }
 217
 218   /**
 219    * {@inheritDoc}
 220    */
 221   @Override
 222   protected void findNewDistances(BinaryNode nodei, BinaryNode nodej,
 223           double dist)
 224   {
 225     double ih = 0;
 226     double jh = 0;
 227
 228     BinaryNode sni = nodei;
 229     BinaryNode snj = nodej;
 230
 231     while (sni != null)
 232     {
 233       ih = ih + sni.dist;
 234       sni = (BinaryNode) sni.left();
 235     }
 236
 237     while (snj != null)
 238     {
 239       jh = jh + snj.dist;
 240       snj = (BinaryNode) snj.left();
 241     }
 242
 243     nodei.dist = ((dist / 2) - ih);
 244     nodej.dist = ((dist / 2) - jh);
 245   }
 246
 247   /***
 248    * not the right place - OH WELL!
 249    */
 250
 251   /**
 252    * Makes a list of groups, where each group is represented by a node whose
 253    * height (distance from the root node), as a fraction of the height of the
 254    * whole tree, is greater than the given threshold. This corresponds to
 255    * selecting the nodes immediately to the right of a vertical line
 256    * partitioning the tree (if the tree is drawn with root to the left). Each
 257    * such node represents a group that contains all of the sequences linked to
 258    * the child leaf nodes.
 259    *
 260    * @param threshold
 261    * @see #getGroups()
 262    */
 263   public List<BinaryNode> groupNodes(float threshold)
 264   {
 265     List<BinaryNode> groups = new ArrayList<BinaryNode>();
 266     _groupNodes(groups, getTopNode(), threshold);
 267     return groups;
 268   }
 269
 270   protected void _groupNodes(List<BinaryNode> groups, BinaryNode nd,
 271           float threshold)
 272   {
 273     if (nd == null)
 274     {
 275       return;
 276     }
 277
 278     if ((nd.height / maxheight) > threshold)
 279     {
 280       groups.add(nd);
 281     }
 282     else
 283     {
 284       _groupNodes(groups, nd.left(), threshold);
 285       _groupNodes(groups, nd.right(), threshold);
 286     }
 287   }
 288
 289   /**
 290    * DOCUMENT ME!
 291    *
 292    * @param nd
 293    *          DOCUMENT ME!
 294    *
 295    * @return DOCUMENT ME!
 296    */
 297   public double findHeight(BinaryNode nd)
 298   {
 299     if (nd == null)
 300     {
 301       return maxheight;
 302     }
 303
 304     if ((nd.left() == null) && (nd.right() == null))
 305     {
 306       nd.height = ((BinaryNode) nd.parent()).height + nd.dist;
 307
 308       if (nd.height > maxheight)
 309       {
 310         return nd.height;
 311       }
 312       else
 313       {
 314         return maxheight;
 315       }
 316     }
 317     else
 318     {
 319       if (nd.parent() != null)
 320       {
 321         nd.height = ((BinaryNode) nd.parent()).height + nd.dist;
 322       }
 323       else
 324       {
 325         maxheight = 0;
 326         nd.height = (float) 0.0;
 327       }
 328
 329       maxheight = findHeight((BinaryNode) (nd.left()));
 330       maxheight = findHeight((BinaryNode) (nd.right()));
 331     }
 332
 333     return maxheight;
 334   }
 335
 336   /**
 337    * Search for leaf nodes below (or at) the given node
 338    *
 339    * @param top2
 340    *          root node to search from
 341    *
 342    * @return
 343    */
 344   public Vector<BinaryNode> findLeaves(BinaryNode top2)
 345   {
 346     Vector<BinaryNode> leaves = new Vector<BinaryNode>();
 347     findLeaves(top2, leaves);
 348     return leaves;
 349   }
 350
 351   /**
 352    * Search for leaf nodes.
 353    *
 354    * @param nd
 355    *          root node to search from
 356    * @param leaves
 357    *          Vector of leaves to add leaf node objects too.
 358    *
 359    * @return Vector of leaf nodes on binary tree
 360    */
 361   Vector<BinaryNode> findLeaves(BinaryNode nd, Vector<BinaryNode> leaves)
 362   {
 363     if (nd == null)
 364     {
 365       return leaves;
 366     }
 367
 368     if ((nd.left() == null) && (nd.right() == null)) // Interior node
 369     // detection
 370     {
 371       leaves.addElement(nd);
 372
 373       return leaves;
 374     }
 375     else
 376     {
 377       /*
 378        * TODO: Identify internal nodes... if (node.isSequenceLabel()) {
 379        * leaves.addElement(node); }
 380        */
 381       findLeaves(nd.left(), leaves);
 382       findLeaves(nd.right(), leaves);
 383     }
 384
 385     return leaves;
 386   }
 387
 388 }