Package com.nr.ci

Source Code of com.nr.ci.Phylagglom

package com.nr.ci;

import static com.nr.NRUtil.*;

import java.io.*;

import org.netlib.util.doubleW;

// Abstract base class for constructing an agglomerative phylogenetic ree
public abstract class Phylagglom {
  public int n, root, fsroot;  // No. of data points, root node, forced root
  public double seqmax, depmax;  // max. values of seq, dep over the tree
  public Phylagglomnode[] t;   // the tree
  public abstract void premin(double[][] d, int[] nextp); //  function called before minimum search 
  public abstract double dminfn(double[][] d, int i, int j); //  function to be minimized 
  public abstract double dbranchfn(double[][] d, int i, int j); //  branch length, node i to mother (j is sister) 
  public abstract double dnewfn(double[][] d, int k, int i, int j, int ni, int nj)// distance function for newly constructed nodes
  public abstract void drootbranchfn(double[][] d, int i, int j, int ni, int nj, doubleW bi, doubleW bj);    // sets branch lengths to the final root node

  public Phylagglom(final double[][] dist){
    this(dist,-1);
  }
 
  public Phylagglom(final double[][] dist, final int fsr) {
    n = dist.length;
    fsroot = fsr;
    t = new Phylagglomnode[2*n-1];
    for(int i=0;i<t.length;i++)
      t[i] = new Phylagglomnode();
  }

  // routine that actually constructs the tree, called by the constructor of a derived class
  public void makethetree(final double[][] dist) {
    int i, j, k, imin=0, jmin=0, node, ntask;
    double dd, dmin;
    double[][] d = buildMatrix(dist)// Matrix d is initialized with dist
    int[] tp = new int[n], nextp = new int[n], prevp = new int[n], tasklist = new int[2*n+1];
    double[] tmp = new double[n];
    for (i=0;i<n;i++) {  // initializations on leaf elements
      // nextp and prevp are for looping on the distance matrix even as it becomes sparse
      nextp[i] = i+1;  
      prevp[i] = i-1;  
     
      tp[i] = i;  // tp points from a distance mattttttttix row to a tree element
      t[i].ldau = t[i].rdau = -1;
      t[i].nel = 1;
    }  
    prevp[0] = nextp[n-1] = -1;   // Signifying end of loop
    // ncurr = n;
    for (node = n; node < 2*n-2; node++) {   // Main loop!
      premin(d,nextp)// Any calculations needed before min finding
      dmin = 9.99e99;
      for (i=0; i>=0; i=nextp[i]) {   // Find i,j pair with min distance
        if (tp[i] == fsroot) continue;
        for (j=nextp[i]; j>=0; j=nextp[j]) {
          if (tp[j] == fsroot) continue;
          if ((dd = dminfn(d,i,j)) < dmin) {
            dmin = dd;
            imin = i; jmin = j;
          }
        }
      }
      i = imin; j = jmin;
        // Now set properties of the parent and children
      t[tp[i]].mo = t[tp[j]].mo = node;
      t[tp[i]].modist = dbranchfn(d,i,j);
      t[tp[j]].modist = dbranchfn(d,j,i);
      t[node].ldau = tp[i];
      t[node].rdau = tp[j];
      t[node].nel = t[tp[i]].nel + t[tp[j]].nel;
      for (k=0; k>=0; k=nextp[k]) {  // Get new-node distances
        tmp[k] = dnewfn(d,k,i,j,t[tp[i]].nel,t[tp[j]].nel);
      }
      for (k=0; k>=0; k=nextp[k]) d[i][k] = d[k][i] = tmp[k];
      // new node replaces child i in dist. matrix, while child gets patched around
      tp[i] = node; 
      if (prevp[j] >= 0) nextp[prevp[j]] = nextp[j];
      if (nextp[j] >= 0) prevp[nextp[j]] = prevp[j];
     
    // End of main loop
   
    i = 0; j = nextp[0]// set properties of the root node
    root = node;
    t[tp[i]].mo = t[tp[j]].mo = t[root].mo = root;
   
    doubleW bi = new doubleW(t[tp[i]].modist);
    doubleW bj = new doubleW(t[tp[j]].modist);
    drootbranchfn(d,i,j,t[tp[i]].nel,t[tp[j]].nel, bi,bj);
      //t[tp[i]].modist,t[tp[j]].modist);
    t[tp[i]].modist = bi.val;
    t[tp[j]].modist = bj.val;
   
    t[root].ldau = tp[i];
    t[root].rdau = tp[j];
    t[root].modist = t[root].dep = 0.;
    t[root].nel = t[tp[i]].nel + t[tp[j]].nel;
   
    // we now traverse the tree computing seq and dep, hints for where to plot nodes in a two-dimensional representation. See Numerical Recipes book text
    ntask = 0;
    seqmax = depmax = 0.;
    tasklist[ntask++] = root;
    while (ntask > 0) {
      i = tasklist[--ntask];
      if (i >= 0) {
        t[i].dep = t[t[i].mo].dep + t[i].modist;
        if (t[i].dep > depmax) depmax = t[i].dep;
        if (t[i].ldau < 0) {
          t[i].seq = seqmax++;
        } else {
          tasklist[ntask++] = -i-1;
          tasklist[ntask++] = t[i].ldau;
          tasklist[ntask++] = t[i].rdau;
        }
      } else {
        i = -i-1;
        t[i].seq = 0.5*(t[t[i].ldau].seq + t[t[i].rdau].seq);
      }
    }
  }
 
  public int comancestor(final int leafa, final int leafb) {
    int i, j;
    for (i = leafa; i != root; i = t[i].mo) {
      for (j = leafb; j != root; j = t[j].mo) if (i == j) break;
      if (i == j) break;
    }
    return i;
  }
 
  /*
   Output a pylogenetic tree p in the "Newick" or "New Hampshire" standard format. Text labels
   for the leaf nodes are input as the rows of str, each a null terminated string. The output file name is specified by filename
   */
  public static void newick(Phylagglom p, char[][] str, String filename) throws IOException {
    java.io.PrintWriter OUT =new java.io.PrintWriter(new FileWriter(filename));
    int i, s, ntask = 0, n = p.n, root = p.root;
    int[] tasklist = new int[2*n+1];
    tasklist[ntask++] = (1 << 16) + root;
    while (ntask-- > 0) {
      s = tasklist[ntask] >> 16;
      i = tasklist[ntask] & 0xffff;
      if (s == 1 || s == 2) {
        tasklist[ntask++] = ((s+2) << 16) + p.t[i].mo;
        if (p.t[i].ldau >= 0) {
          OUT.printf("(");
          tasklist[ntask++] = (2 << 16) + p.t[i].rdau;     
          tasklist[ntask++] = (1 << 16) + p.t[i].ldau;     
        }
        else OUT.printf("%s:%f",new String(str[i], 0, str[i].length-1),p.t[i].modist);
      }
      else if (s == 3) {if (ntask > 0) OUT.printf(",\n");}
      else if (s == 4) {
        if (i == root) OUT.printf(");\n");
        else OUT.printf("):%f",p.t[i].modist);
      }
    }
    OUT.close();
  }
 
  public static void phyl2ps(String filename, Phylagglom ph, char[][] str, int extend,
    double xl, double xr, double yt, double ybthrows IOException{
    int i,j;
    double id,jd,xi,yi,xj,yj,seqmax,depmax;
    java.io.PrintWriter OUT =new java.io.PrintWriter(new FileWriter(filename));
    OUT.printf("%%!PS\n/Courier findfont 8 scalefont setfont\n");
    seqmax = ph.seqmax;
    depmax = ph.depmax;
    for (i=0; i<2*(ph.n)-1; i++) {
      j = ph.t[i].mo;
      id = ph.t[i].dep;
      jd = ph.t[j].dep;
      xi = xl + (xr-xl)*id/depmax;
      yi = yt - (yt-yb)*(ph.t[i].seq+0.5)/seqmax;
      xj = xl + (xr-xl)*jd/depmax;
      yj = yt - (yt-yb)*(ph.t[j].seq+0.5)/seqmax;
      OUT.printf("%f %f moveto %f %f lineto %f %f lineto 0 setgray stroke\n",
        xj,yj,xj,yi,xi,yi);
      if (extend!=0) {
        if (i < ph.n) {
          OUT.printf("%f %f moveto %f %f lineto 0.7 setgray stroke\n",
            xi,yi,xr,yi);
          OUT.printf("%f %f moveto (%s (%02d)) 0 setgray show\n",
            xr+3.,yi-2.,new String(str[i], 0, str[i].length-1),i);
        }
      } else {
        if (i < ph.n) OUT.printf("%f %f moveto (%s (%02d)) 0 setgray show\n",
          xi+3.,yi-2.,new String(str[i], 0, str[i].length-1),i);
      }
    }
    OUT.printf("showpage\n\004");
    OUT.close();
  }
}
TOP

Related Classes of com.nr.ci.Phylagglom

TOP
Copyright © 2018 www.massapi.com. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.