njtree.h

Go to the documentation of this file.

#ifndef NJTREE_H

#define NJTREE_H
//njtree.h  the class of njtree

#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <fstream>
#include <cstring>

// these are bad C habbits, should get rid of them!
//#define PREC 10              /* precision of branch-lengths  */
#define PRC  20 
#define LEN  50              /* length of taxon names        */

struct word
{
         word() : next(0) { }
    word(const char n[]) { strcpy(name, n); }
    char name[LEN];
    struct word *next;   // to make a chain word->word
};

struct wordList
{
        wordList() :head(0), tail(0) {}
   word *head;
   word *tail;
};

class njtree
{
        public:
      enum inputformat { MATRIX, NUMTAXON };
                njtree() : trees(0), delta(0), n(0) { };   
                ~njtree();
                bool load(std::istream &ins, inputformat ipf);  
      bool load_generic(std::istream &ins);

                int taxonCount() const { return n; }
                int remainingCount() const { return r; }
                bool symmetric();
                void compute_sums_Sx();
                void best_pair(int &a, int &b) const;
                void reduction(int a, int b);
                void finish(std::ostream &ous);
                void output(int i, std::ostream &ous);

        private:
                bool empty(int i) const { return trees[i].head == 0; }
                                                  //ith line is emptied or not
                double distance(int i, int j) const;
                double variance(int i, int j) const;

                double agglomerative_criterion(int i, int j) const;
                double branch_length(int a, int b) const;
                double lamda(int a, int b) const;
                double reduction4(int a, double la, int b, double lb, int i, double lam) const;
                double reduction10(int a, int b, int i, double lam, double vab) const;
                void concatenate(const char chn1[], int ind, int post);
                double finish_branch_length(int i, int j, int k) const;
                void clear();


                double **delta;
                int n;               // size of the matrix
                int r;               // taxons remaining
                //int a, b;            // best pair for reduction
                wordList *trees;
// static members to control the whole class behavior
      static int PREC; 
};


inline double njtree::agglomerative_criterion(int i, int j) const
{
    return  (r-2)*distance(i,j) - delta[i][i] - delta[j][j]; 
         // Formula (1)
}

inline double njtree::distance(int i, int j) const
{ 
        return i>j? delta[i][j] : delta[j][i];
}  //the lower half left trangle

inline double njtree::variance(int i, int j) const
{
        return (i > j)? (delta[j][i]) : (delta[i][j]);
}

inline double njtree::branch_length(int a, int b) const
{
          /* Formula (2) */
    return 0.5*(distance(a,b) + (delta[a][a]-delta[b][b])/(r-2));
}

inline double njtree::reduction4(int a, double la, int b, double lb, int i, double lam) const
{
    return lam*(distance(a,i)-la) +(1-lam)*(distance(b,i)-lb);
            /* Formula (4) */
}

inline double njtree::reduction10(int a, int b, int i, double lam, double vab) const
{
    return lam*variance(a,i)+(1-lam)*variance(b,i) -lam*(1-lam)*vab;
    /* Formula (10)  double Vci */
}

inline double njtree::finish_branch_length(int i, int j, int k) const
{
         return 0.5*(distance(i,j) + distance(i,k) -distance(j,k));
         // branch length of i to common ancestor
}
#endif