Biostatistics 615/815 Winter 2011

From Genome Analysis Wiki
Jump to navigationJump to search

(This page is the Biostatistics 615/815 course web page in Winter 2011 semester)

Objective

In this winter, Biostatistics 615/815 aims for providing students with a practical understanding of computational aspects in implementing statistical methods. Although C++ language will be used throughout the course, using Java programming language for homework and project will be acceptable.

Target Audience

Students in Biostatistics 615 should be comfortable with simple algebra and basic statistics including probability distribution, linear model, and hypothesis testing. Previous experience in programming is not required, but those who do not have previous programming experience should expect to spend additional time studying and learning to be familiar with a programming language during the coursework. Most students registering for the course are Masters or Doctoral students in Biostatistics, Statistics, Bioinformatics or Human Genetics.

Students in Biostatistics 815 should be familiar with programming languages so that they can complete the class project tackling an advanced statistical problem during the semester. Project will be carried out in teams of 2. The details of the possible projects will be announced soon.

Textbook

  • Required Textbook : Cormen, Leiserson, Rivest, and Stein, "Introduction to Algorithms", Third Edition, The MIT Press, 2009 [Official Book Web Site]
  • Optional Textbook : Press, Teukolsky, Vetterling, Flannery, "Numerical Recipes", 3rd Edition, Cambridge University Press, 2007 [Official Book Web Site]

Class Schedule

Classes are scheduled for Tuesday and Thursdays, 8:30 - 10:00 am at SPH II M4318

Topics

The following contents are planned to be covered.

Part I : Algorithms 101

  • Understanding of Computational Time Complexity
  • Sorting
  • Divide and Conquer Algorithms
  • Searching
  • Key Data Structure
  • Dynamic Programming

Part II : Matrix Operations and Numerical Optimizations

  • Matrix decomposition (LU, QR, SVD)
  • Implementation of Linear Models
  • Numerical Optimizations

Part III : Advanced Statistical Methods

  • Hidden Markov Models
  • Expectation Maximization
  • Markov-Chain Monte Carlo Methods

Class Notes

Problem Sets


Source code for Homework 6 - simplex615.h

#ifndef __SIMPLEX_615_H
#define __SIMPLEX_615_H

#include <vector>
#include <cmath>
#include <iostream>

#define ZEPS 1e-10

class optFunc {
 public:
  virtual double operator() (std::vector<double>& x) = 0;
};

// Simplex contains (dim+1)*dim points
class simplex615 {
 protected:
  std::vector<std::vector<double> > X;
  std::vector<double> Y;
  std::vector<double> midPoint;
  std::vector<double> thruLine;

  int dim, idxLo, idxHi, idxNextHi;

  void evaluateFunction(optFunc& foo);
  void evaluateExtremes();
  void prepareUpdate();
  bool updateSimplex(optFunc& foo, double scale);
  void contractSimplex(optFunc& foo);
  static int check_tol(double fmax, double fmin, double ftol);

 public:
  simplex615(double* p, int d);
  void amoeba(optFunc& foo, double tol);
  std::vector<double>& xmin();
  double ymin();
};

simplex615::simplex615(double* p, int d) : dim(d) {
  X.resize(dim+1);
  Y.resize(dim+1);
  midPoint.resize(dim);
  thruLine.resize(dim);
  for(int i=0; i < dim+1; ++i) {
    X[i].resize(dim);
  }

  // set every point the same
  for(int i=0; i < dim+1; ++i) {
    for(int j=0; j < dim; ++j) {
      X[i][j] = p[j];
    }
  }
  
  // then increase each dimension by one except for the last point
  for(int i=0; i < dim; ++i) {
    X[i][i] += 1.;
  }
}

void simplex615::evaluateFunction(optFunc& foo) {
  for(int i=0; i < dim+1; ++i) {
    Y[i] = foo(X[i]);
  }
}

void simplex615::evaluateExtremes() {
  if ( Y[0] > Y[1] ) {
    idxHi = 0; 
    idxLo = idxNextHi = 1;
  }
  else {
    idxHi = 1;
    idxLo = idxNextHi = 0;
  }
  
  for(int i=2; i < dim+1; ++i) {
    if ( Y[i] <= Y[idxLo] ) {
      idxLo = i;
    }
    else if ( Y[i] > Y[idxHi] ) {
      idxNextHi = idxHi;
      idxHi = i;
    }
    else if ( Y[i] > Y[idxNextHi] ) {
      idxNextHi = i;
    }
  }
}

void simplex615::prepareUpdate() {
  for(int j=0; j < dim; ++j) {
    midPoint[j] = 0;
  }
  for(int i=0; i < dim+1; ++i) {
    if ( i != idxHi ) {
      for(int j=0; j < dim; ++j) {
	midPoint[j] += X[i][j];
      }
    }
  }
  for(int j=0; j < dim; ++j) {
    midPoint[j] /= dim;
    thruLine[j] = X[idxHi][j] - midPoint[j];
  }
}

bool simplex615::updateSimplex(optFunc& foo, double scale) {
  std::vector<double> nextPoint;
  nextPoint.resize(dim);
  for(int i=0; i < dim; ++i) {
    nextPoint[i] = midPoint[i] + scale * thruLine[i];
  }
  double fNext = foo(nextPoint);
  if ( fNext < Y[idxHi] ) { // exchange with maximum
    for(int i=0; i < dim; ++i) {
      X[idxHi][i] = nextPoint[i];
    }
    Y[idxHi] = fNext;
    return true;
  }
  else {
    return false;
  }
}

void simplex615::contractSimplex(optFunc& foo) {
  for(int i=0; i < dim+1; ++i) {
    if ( i != idxLo ) {
      for(int j=0; j < dim; ++j) {
	X[i][j] = 0.5*( X[idxLo][j] + X[i][j] );
	Y[i] = foo(X[i]);
      }
    }
  }
}

void simplex615::amoeba(optFunc& foo, double tol) {
  evaluateFunction(foo);
  while(true) {
    evaluateExtremes();
    prepareUpdate();
    
    if ( check_tol(Y[idxHi],Y[idxLo],tol) ) break;

    updateSimplex(foo, -1.0); // reflection
    
    if ( Y[idxHi] < Y[idxLo] ) {
      updateSimplex(foo, -2.0); // expansion
    }
    else if ( Y[idxHi] >= Y[idxNextHi] ) {
      if ( !updateSimplex(foo, 0.5) ) {
	contractSimplex(foo);
      }
    }
  }
}

std::vector<double>& simplex615::xmin() {
  return X[idxLo];
}

double simplex615::ymin() {
  return Y[idxLo];
}

int simplex615::check_tol(double fmax, double fmin, double ftol) {
  double delta = fabs(fmax - fmin); double accuracy = (fabs(fmax) + fabs(fmin)) * ftol;
  return (delta < (accuracy + ZEPS));
}

#endif // __SIMPLEX_615_H


Source code for Homework 4 - Matrix615.h

#ifndef __BIOSTAT615_MATRIX_H  // avoid including the same header twice
#define __BIOSTAT615_MATRIX_H

#include <iostream>
#include <fstream>
#include <boost/tokenizer.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/foreach.hpp>

#include <Eigen/Core>

// a generic class for Matrix
template<class T>
class Matrix615 {
protected:   // internal data
  std::vector<T> data;    // using std::vector - object copy is now possible
  int nr, nc;             // # rows and cols
  bool hasMissing;
  T valueMissing;
  std::string strMissing;
public:
  // default constructor
 Matrix615() : nr(0), nc(0), hasMissing(false) {}
 Matrix615(const char* filename) : nr(0), nc(0), hasMissing(false) {
    readFromFile(filename);
  }
 Matrix615(const T value, const char* str) : nr(0), nc(0) {
    enableMissingValue(value, str);
  }

  // Allow missing value as a pair of actual value and string value
  void enableMissingValue(const T value, const char* str) {
    hasMissing = true;
    valueMissing = value;
    strMissing = str;
  }

  // resize the dimension of the matrix
  void resize(int nrows, int ncols) {
    nr = nrows;
    nc = ncols;
    data.resize(nr*nc);
  }

  // fill the content 
  void fill(T defaultValue) {
    std::fill( data.begin(), data.end(), defaultValue );
  }

  void copyTo(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m) {
    m.resize(nr,nc);
    for(int j=0; j < nc; ++j) {
      for(int i=0; i < nr; ++i) {
	m(i,j) = data[i*nc+j];
      }
    }
  }

  // access individual element
  T& at(int r, int c) { return data[r*nc+c]; }

  int numRows() { return nr; }
  int numCols() { return nc; }

  // print the content of the matrix 
  void print(std::ostream& o) {
    for(int i=0; i < nr; ++i) {
      for(int j=0; j < nc; ++j) {
	if ( j > 0 ) o << "\t";
	if ( hasMissing && ( valueMissing == at(i,j) ) )
	  o << strMissing;
	else
	  o << at(i,j);
      }
      o << std::endl;
    }
  }

  // opens a file to fill the matrix
  void readFromFile(const char* file) {
    std::ifstream ifs(file);
    if ( ! ifs.is_open() ) {
      std::cerr << "Cannot open file " << file << std::endl;
      abort();
    }
    
    std::string line;
    boost::char_separator<char> sep(" \t");
    typedef boost::tokenizer< boost::char_separator<char> > wsTokenizer;
  
    data.clear();  
    nr = nc = 0;
    while( std::getline(ifs, line) ) {
      wsTokenizer t(line,sep);
      for(wsTokenizer::iterator i=t.begin(); i != t.end(); ++i) {
	// if hasMissing is set, convert string "Missing" into special value for Missing
	if ( hasMissing && ( i->compare(strMissing) == 0 ) )
	    data.push_back(valueMissing);
	// Otherwise, convert the string to a particular type
	else 
	  data.push_back(boost::lexical_cast<T>(i->c_str()));  
	if ( nr == 0 ) ++nc;  // count # of columns at the first row
      }
      ++nr;
      // when reading each line, make sure that the # of columns match to expectation;
      if ( (int)data.size() != nr*nc ) {
	std::cerr << "The input file is not rectangle at line " << nr << std::endl;
	abort();
      }
    }
  }
};

#endif

Source code for Homework 4 - Problem 1

#include "Matrix615.h"

#define NIL 999999

int main(int argc, char** argv) {
  if ( argc != 2 ) {
    std::cerr << "Usage: floydWarshall [weight matrix]" << std::endl;
    abort();
  }

  Matrix615<int> W;
  W.enableMissingValue(NIL,"NIL");
  W.readFromFile(argv[1]);
  if ( W.numRows() != W.numCols() ) {
    std::cerr << "The input file is not exactly square" << std::endl;
    abort();
  }

  int n = W.numRows();
  Matrix615<int> P;
  P.enableMissingValue(NIL,"NIL");
  P.resize(n,n);

  // Fill out the code in the part marked as *** [FILL HERE] ***

  // *** FILL HERE ** initialize P and W matrix as described in page 695-696 of the textbook

  Matrix615<int> D = W;  // object copy is valid because no pointer is involved

  // print initial D and W matrix
  std::cout << "---------- Initial D  Matrix ---------------------------------" << std::endl;
  D.print(std::cout);
  std::cout << "---------- Initial PI Matrix ---------------------------------" << std::endl;
  P.print(std::cout);

  for(int k=0; k < n; ++k) {
    // *** FILL HERE *** Run floyd-Warshall algorithm for each k and 

    // print out D and PI matrix for each iteration
    std::cout << "---------- D  Matrix after covering node " << k << "----------" << std::endl;
    D.print(std::cout);
    std::cout << "---------- PI Matrix after covering node " << k << "----------" << std::endl;
    P.print(std::cout);
  }

  // *** EXTRA POINTS ***  Extra points (10\%) will be given if the optimal path 
  //                       for every pair of node is printed below

  return 0;
}

Source code for Homework 4 - Problem 2

#include <cstdlib>
#include <iomanip>
#include <iostream>

#include "Matrix615.h"

#define N_STATES 2
#define N_DATA 2

const char* stateLabels[N_STATES] = {"F","B"};
const char* dataLabels[N_DATA] = {"H","T"};

class BiasedCoinHMM {
protected:
  int T;
  // HMM initial parameters
  std::vector<double> pi; // N_STATES * 1 vector
  Matrix615<double> A;    // N_STATES * N_STATES matrix 
                          //  : A_{ij} is \Pr(q_t=i|q_{t-1}=j)
  Matrix615<double> B;    // N_OBS * N_STATES matrix
                          //  : B_{ij} = b_j(o_i) = \Pr(o_t=i|q_t=j)

  // Data (observations)
  std::vector<int> o; // vector observations 

  // forward-backward probability
  Matrix615<double> alphas; // T * N_STATES : alphas[i,j] = alpha_i(i)
  Matrix615<double> betas;  // T * N_STATES : betas[i,j]  = beta_j(i)
  Matrix615<double> gammas; // T * N_STATES : gammas[i,j] = gammas_j(i) = Pr(q_t=i|o_t=i)

  // viterbi probability and paths
  Matrix615<double> deltas; // T * N_STATES
  Matrix615<int> phis;   // T * N_STATES
  std::vector<int> mleStates; // vector of MLE states

public:
  BiasedCoinHMM(double priors[N_STATES], double trans[N_STATES][N_STATES], double emis[N_DATA][N_STATES]) {
    for(int i=0; i < N_STATES; ++i) {
      pi.push_back(priors[i]);
    }

    A.resize(N_STATES,N_STATES);

    for(int i=0; i < N_STATES; ++i) {
      for(int j=0; j < N_STATES; ++j) {
        A.at(i,j) = trans[i][j];
      }
    }

    B.resize(N_DATA,N_STATES);
    for(int i=0; i < N_DATA; ++i) {
      for(int j=0; j < N_STATES; ++j) {
        B.at(i,j) = emis[i][j];
      }
    }

    T = 0;
  }

  void loadObservations(std::vector<int>& obs) {
    o = obs;
    T = o.size();

    alphas.resize(T,N_STATES);
    betas.resize(T,N_STATES);
    gammas.resize(T,N_STATES);
    deltas.resize(T,N_STATES);
    phis.resize(T,N_STATES);
    mleStates.resize(T);
  }

  void computeForwardBackward() {
    // *** FILL OUT *** to compute 
    // forward and backward probabilities into alphas, betas, and gammas 
  }

  void computeViterbiPath() {
    // *** FILL OUT *** to compute
    // viterbi path and likelihood into deltas, phis, and mleStates,
  }

  void outputResults(std::ostream& os, std::vector<int>& trueStates) {
    if ( T != (int)trueStates.size() ) {
      std::cerr << "True states are in different length with HMM" << std::endl;
      abort();
    }

    os << "#SEQ\tTRUE_S\tOBS\tP(q|o)\tMLE_S" << std::endl << std::fixed << std::setprecision(4);
    for(int t = 0; t < T; ++t) {
      os << t+1 << "\t" 
         << stateLabels[trueStates[t]] << "\t"
         << dataLabels[o[t]] << "\t"
         << gammas.at(t,1) << "\t"  // prints Pr(q_t=1|o)
         << stateLabels[mleStates[t]] << "\n";
    }
  }
};

int main(int argc, char** argv) {
  if ( argc != 2 ) {
    std::cerr << "Usage: coinHMM [inputFile]" << std::endl;
    return -1;
  }

  std::vector<int> trueStates;
  std::vector<int> observations;
  std::ifstream ifs(argv[1]);

  if ( !ifs.is_open() ) {
    std::cerr << "Cannot open file " << argv[1] << std::endl;
    return -1;    
  }

  std::string tok;
  for(int i=0; ifs >> tok; ++i) {
    if ( i % 3 == 1 ) {
      if ( tok.compare("F") == 0 ) {
        trueStates.push_back(0);
      }
      else if ( tok.compare("B") == 0 ) {
        trueStates.push_back(1);
      }
      else {
        std::cerr << "Cannot recognize state " << tok << std::endl;     
      }
    }
    else if ( i % 3 == 2 ) {
      if ( tok.compare("H") == 0 ) {
        observations.push_back(0);
      }
      else if ( tok.compare("T") == 0 ) {
        observations.push_back(1);
      }
      else {
        std::cerr << "Cannot recognize observation " << tok << std::endl;       
      }
    }
  }

  std::cout << "Finished reading " << trueStates.size() << " states and observations.." << std::endl;

  double trans[N_STATES][N_STATES] = { {0.95,0.2}, {0.05,0.8} };
  double emis[N_DATA][N_STATES] = { {0.5,0.9}, {0.5,0.1} };
  double pi[N_STATES] = {0.9,0.1};

  BiasedCoinHMM bcHMM(pi, trans, emis);

  bcHMM.loadObservations(observations);
  bcHMM.computeForwardBackward();
  bcHMM.computeViterbiPath();
  bcHMM.outputResults(std::cout, trueStates);
  return 0;
}

Additional example of containing negative weights

$ cat ./sampleInput2.txt
0	3	8	NIL	-4
NIL	0	NIL	1	7
NIL	4	0	NIL	NIL
2	NIL	-5	0	NIL
NIL	NIL	NIL	6	0

$ ./floydWarshall ./sampleInput2.txt
---------- Initial D  Matrix ---------------------------------
0	3	8	NIL	-4
NIL	0	NIL	1	7
NIL	4	0	NIL	NIL
2	NIL	-5	0	NIL
NIL	NIL	NIL	6	0
---------- Initial PI Matrix ---------------------------------
NIL	0	0	NIL	0
NIL	NIL	NIL	1	1
NIL	2	NIL	NIL	NIL
3	NIL	3	NIL	NIL
NIL	NIL	NIL	4	NIL
---------- D  Matrix after covering node 0----------
0	3	8	NIL	-4
NIL	0	NIL	1	7
NIL	4	0	NIL	999995
2	5	-5	0	-2
NIL	NIL	NIL	6	0
---------- PI Matrix after covering node 0----------
NIL	0	0	NIL	0
NIL	NIL	NIL	1	1
NIL	2	NIL	NIL	0
3	0	3	NIL	0
NIL	NIL	NIL	4	NIL
---------- D  Matrix after covering node 1----------
0	3	8	4	-4
NIL	0	NIL	1	7
NIL	4	0	5	11
2	5	-5	0	-2
NIL	NIL	NIL	6	0
---------- PI Matrix after covering node 1----------
NIL	0	0	1	0
NIL	NIL	NIL	1	1
NIL	2	NIL	1	1
3	0	3	NIL	0
NIL	NIL	NIL	4	NIL
---------- D  Matrix after covering node 2----------
0	3	8	4	-4
NIL	0	NIL	1	7
NIL	4	0	5	11
2	-1	-5	0	-2
NIL	NIL	NIL	6	0
---------- PI Matrix after covering node 2----------
NIL	0	0	1	0
NIL	NIL	NIL	1	1
NIL	2	NIL	1	1
3	2	3	NIL	0
NIL	NIL	NIL	4	NIL
---------- D  Matrix after covering node 3----------
0	3	-1	4	-4
3	0	-4	1	-1
7	4	0	5	3
2	-1	-5	0	-2
8	5	1	6	0
---------- PI Matrix after covering node 3----------
NIL	0	3	1	0
3	NIL	3	1	0
3	2	NIL	1	0
3	2	3	NIL	0
3	2	3	4	NIL
---------- D  Matrix after covering node 4----------
0	1	-3	2	-4
3	0	-4	1	-1
7	4	0	5	3
2	-1	-5	0	-2
8	5	1	6	0
---------- PI Matrix after covering node 4----------
NIL	2	3	4	0
3	NIL	3	1	0
3	2	NIL	1	0
3	2	3	NIL	0
3	2	3	4	NIL
Optimal path for 1 <- 0 (d = 1) : 1 <-(4)-- 2 <-(-5)-- 3 <-(6)-- 4 <-(-4)-- 0
Optimal path for 2 <- 0 (d = -3) : 2 <-(-5)-- 3 <-(6)-- 4 <-(-4)-- 0
Optimal path for 3 <- 0 (d = 2) : 3 <-(6)-- 4 <-(-4)-- 0
Optimal path for 4 <- 0 (d = -4) : 4 <-(-4)-- 0
Optimal path for 0 <- 1 (d = 3) : 0 <-(2)-- 3 <-(1)-- 1
Optimal path for 2 <- 1 (d = -4) : 2 <-(-5)-- 3 <-(1)-- 1
Optimal path for 3 <- 1 (d = 1) : 3 <-(1)-- 1
Optimal path for 4 <- 1 (d = -1) : 4 <-(-4)-- 0 <-(2)-- 3 <-(1)-- 1
Optimal path for 0 <- 2 (d = 7) : 0 <-(2)-- 3 <-(1)-- 1 <-(4)-- 2
Optimal path for 1 <- 2 (d = 4) : 1 <-(4)-- 2
Optimal path for 3 <- 2 (d = 5) : 3 <-(1)-- 1 <-(4)-- 2
Optimal path for 4 <- 2 (d = 3) : 4 <-(-4)-- 0 <-(2)-- 3 <-(1)-- 1 <-(4)-- 2
Optimal path for 0 <- 3 (d = 2) : 0 <-(2)-- 3
Optimal path for 1 <- 3 (d = -1) : 1 <-(4)-- 2 <-(-5)-- 3
Optimal path for 2 <- 3 (d = -5) : 2 <-(-5)-- 3
Optimal path for 4 <- 3 (d = -2) : 4 <-(-4)-- 0 <-(2)-- 3
Optimal path for 0 <- 4 (d = 8) : 0 <-(2)-- 3 <-(6)-- 4
Optimal path for 1 <- 4 (d = 5) : 1 <-(4)-- 2 <-(-5)-- 3 <-(6)-- 4
Optimal path for 2 <- 4 (d = 1) : 2 <-(-5)-- 3 <-(6)-- 4
Optimal path for 3 <- 4 (d = 6) : 3 <-(6)-- 4

Source code from Homework 2 - Problem 3

#include <iostream>
#include <vector>
#include <ctime>
#include <fstream>
#include <set>
#include "mySortedArray.h"
#include "myTree.h"
#include "myList.h"                                            

int main(int argc, char** argv) {
  int tok;
  std::vector<int> v;
  if ( argc > 1 ) {
    std::ifstream fin(argv[1]);
    while( fin >> tok ) { v.push_back(tok); }
    fin.close();
  }
  else {
    while( std::cin >> tok ) { v.push_back(tok); }
  }                                             
  mySortedArray<int> c1;
  myList<int> c2;
  myTree<int> c3;
  std::set<int> s;

  clock_t start = clock();
  for(int i=0; i < (int)v.size(); ++i) {
    c1.insert(v[i]);
  }
  clock_t finish = clock();
  double duration = (double)(finish-start)/CLOCKS_PER_SEC;  
  std::cout << "Sorted Array (Insert) " << duration << std::endl;

  start = clock();
  for(int i=0; i < (int)v.size(); ++i) {
    c2.insert(v[i]);
  }
  finish = clock();
  duration = (double)(finish-start)/CLOCKS_PER_SEC;  
  std::cout << "List (Insert) " << duration << std::endl;
 
  start = clock();
  for(int i=0; i < (int)v.size(); ++i) {
    c3.insert(v[i]);
  }
  finish = clock();
  duration = (double)(finish-start)/CLOCKS_PER_SEC;  
  std::cout << "Tree (Insert) " << duration << std::endl;

  start = clock();
  for(int i=0; i < (int)v.size(); ++i) {
    s.insert(v[i]);
  }
  finish = clock();
  duration = (double)(finish-start)/CLOCKS_PER_SEC;  
  std::cout << "std::set (Insert) " << duration << std::endl;

  start = clock();
  for(int i=0; i < (int)v.size(); ++i) {
    c1.search(v[i]);
  }
  finish = clock();
  duration = (double)(finish-start)/CLOCKS_PER_SEC;  
  std::cout << "Sorted Array (Search) " << duration << std::endl;

  start = clock();
  for(int i=0; i < (int)v.size(); ++i) {
    c2.search(v[i]);
  }
  finish = clock();
  duration = (double)(finish-start)/CLOCKS_PER_SEC;  
  std::cout << "List (Search) " << duration << std::endl;

  start = clock();
  for(int i=0; i < (int)v.size(); ++i) {
    c3.search(v[i]);
  }
  finish = clock();
  duration = (double)(finish-start)/CLOCKS_PER_SEC;  
  std::cout << "Tree (Search) " << duration << std::endl;

  start = clock();
  for(int i=0; i < (int)v.size(); ++i) {
    s.find(v[i]);
  }
  finish = clock();
  duration = (double)(finish-start)/CLOCKS_PER_SEC;  
  std::cout << "std::set (Search) " << duration << std::endl;
}

The header file of mySortedArray.h (For Homework 2 Problem 3)

#include <iostream>
#define DEFAULT_ALLOC 1024
template <class T> // template supporting a generic type
class mySortedArray {
 protected:     // member variables hidden from outside
   T *data;    // array of the genetic type
   int size;   // number of elements in the container
   int nalloc; // # of objects allocated in the memory
  mySortedArray(mySortedArray& a) {};
  int search(const T& x, int begin, int end);
public:      // abstract interface visible to outside
   mySortedArray();         // default constructor
   ~mySortedArray();        // desctructor
   void insert(const T& x); // insert an element x
   int search(const T& x);  // search for an element x and return its location
   bool remove(const T& x); // delete a particular element
};

template <class T>
mySortedArray<T>::mySortedArray() {     // default constructor
  size = 0;              // array do not have element initially
  nalloc = DEFAULT_ALLOC;
  data = new T[nalloc];  // allocate default # of objects in memory
}

template <class T>
mySortedArray<T>::~mySortedArray() {    // destructor
  if ( data != NULL ) {
    delete [] data;      // delete the allocated memory before destorying
  }                      // the object. otherwise, memory leak happens
}

template <class T>
void mySortedArray<T>::insert(const T& x) {
  if ( size >= nalloc ) {  // if container has more elements than allocated
    T* newdata = new T[nalloc*2];   // make an array at doubled size
    for(int i=0; i < nalloc; ++i) {
      newdata[i] = data[i];         // copy the contents of array
    }
    delete [] data;                 // delete the original array
    data = newdata;                 // and reassign data ptr
  }

  int i;
  for(i=size-1; (i >= 0) && (data[i] > x); --i) {
    data[i+1] = data[i];            // insert the list into right position
  }
  data[i+1] = x;
  ++size;                           // increase the size
}

template <class T>
int mySortedArray<T>::search(const T& x) {
  return search(x, 0, size-1);
}

template <class T>
int mySortedArray<T>::search(const T& x, int begin, int end) {
  if ( begin > end ) {
    return -1;
  }
  else {
    int mid = (begin+end)/2;
    if ( data[mid] == x ) {
      return mid;
    }
    else if ( data[mid] < x ) {
      return search(x, mid+1, end);
    }
    else {
      return search(x, begin, mid);
    }
  }
}

template <class T>
bool mySortedArray<T>::remove(const T& x) {
  int i = search(x);  // try to find the element
  if ( i >= 0 ) {      // if found
    for(int j=i; j < size-1; ++j) {
      data[i] = data[i+1];  // shift all the elements by one
    }
    --size;           // and reduce the array size
    return true;      // successfully removed the value
  }
  else {
    return false;     // cannot find the value to remove
  }
}

Example code to generare all possible permutations

Note that http://www.cplusplus.com/reference/algorithm/next_permutation/ provides much simpler solution for generating permutation. (Thanks to Matthew Snyder) 

// Code to generate permutations from 1..n
// These code was adopted from
// http://geeksforgeeks.org/?p=767 by Hyun Min Kang
#include <iostream>
#include <vector>

// swaps two elements
void swap(std::vector<int>& v, int i, int j) {
  int tmp = v[i];
  v[i] = v[j];
  v[j] = tmp;
  return;
}

// actual engine - recursively calls permutation
void permute(std::vector<int>& v, int from, int to) {
  if ( from == to ) {           // print the permutation
    for(int i=0; i < v.size(); ++i) {
      std::cout << " " << v[i];
    }
    std::cout << std::endl;
  }
  else {
    for(int i = from; i <= to; ++i) {
      swap(v,from,i);           // swaps two elements
      permute(v, from+1, to);   // permute the rest
      swap(v,from,i);           // recover the vector to the original state
    }
  }
}

// Permutation Example
int main(int argc, char** argv) {
  if ( argc != 2 ) {
    std::cerr << "Usage: " << argv[0] << " " << "[# of samples]" << std::endl;
    return -1;
  }

  // takes input from 1..n
  int n = atoi(argv[1]);
  std::vector<int> input;
  for(int i=1; i <= n; ++i) {
    input.push_back(i);
  }
  permute(input, 0, n-1); // print all possible permutations
  return 0;
}

#include <iostream>                       // for input/output
#include <boost/graph/adjacency_list.hpp> // for using graph type 
#include <boost/graph/dijkstra_shortest_paths.hpp> // for dijkstra algorithm
using namespace std;   // allow to omit prefix 'std::'
using namespace boost; // allow to omit prefix 'boost::'

int getNodeID(int node) {
  return ((node/5)+1)*10+(node%5+1);
}

Code for Problem 3 - Problem 2

#include <iostream>                       // for input/output
#include <boost/graph/adjacency_list.hpp> // for using graph type 
#include <boost/graph/dijkstra_shortest_paths.hpp> // for dijkstra algorithm
using namespace std;   // allow to omit prefix 'std::'
using namespace boost; // allow to omit prefix 'boost::'

int getNodeID(int node) {
  return ((node/5)+1)*10+(node%5+1);
}

int main(int argc, char** argv) {
  // defining a graph type
  // 1. edges are stored as std::list internally
  // 2. verticies are stored as std::vector internally
  // 3. the graph is directed (undirectedS, bidirectionalS can be used)
  // 4. vertices do not carry particular properties
  // 5. edges contains weight property as integer value
  typedef adjacency_list< listS, vecS, directedS, no_property, property< edge_weight_t, int> > graph_t;

  // vertex_descriptor is a type for representing vertices
  typedef graph_traits< graph_t >::vertex_descriptor vertex_descriptor;
  // a nodes is represented as an integer, and an edge is a pair of integers
  typedef std::pair<int, int> E;

  // Connect between verticies as in the Manhattan Tourist Problem
  // Each node is numbered as a two-digit integer of [#row] and [#col]
  enum { N11, N12, N13, N14, N15,
	 N21, N22, N23, N24, N25,
	 N31, N32, N33, N34, N35,
	 N41, N42, N43, N44, N45,
	 N51, N52, N53, N54, N55 };

  E edges [] = { E(N11,N12), E(N12,N13), E(N13,N14), E(N14,N15),
		 E(N21,N22), E(N22,N23), E(N23,N24), E(N24,N25),
		 E(N31,N32), E(N32,N33), E(N33,N34), E(N34,N35),
		 E(N41,N42), E(N42,N43), E(N43,N44), E(N44,N45),
		 E(N51,N52), E(N52,N53), E(N53,N54), E(N54,N55),
		 E(N11,N21), E(N12,N22), E(N13,N23), E(N14,N24), E(N15,N25),
		 E(N21,N31), E(N22,N32), E(N23,N33), E(N24,N34), E(N25,N35),
		 E(N31,N41), E(N32,N42), E(N33,N43), E(N34,N44), E(N35,N45),
		 E(N41,N51), E(N42,N52), E(N43,N53), E(N44,N54), E(N45,N55) };
  // Assign weights for each edge
  int weight [] = { 4, 2, 0, 7,     // horizontal weights
		    7, 4, 5, 9,
		    6, 8, 1, 0,
		    1, 6, 4, 7,
		    1, 5, 8, 5,
		    0, 6, 6, 2, 4,  // vertical weights
		    9, 7, 1, 0, 6,
		    1, 8, 4, 8, 9,
		    3, 6, 6, 0, 7 };

  // define a graph as an array of edges and weights
#if defined(BOOST_MSVC) && BOOST_MSVC <= 1300 // special routine for MS VC++
  typedef graph_traits < graph_t >::edge_descriptor edge_descriptor;
  graph_t g(25);
  property_map<graph_t, edge_weight_t>::type weightmap = get(edge_weight, g);
  for (std::size_t j = 0; j < sizeof(edges)/sizeof(E); ++j) {
    edge_descriptor e; 
    bool inserted;
    boost::tie(e, inserted) = add_edge(edge_array[j].first, edge_array[j].second, g);
    weightmap[e] = weights[j];
  }
#else  // for regulat compilers
  graph_t g(edges, edges + sizeof(edges) / sizeof(E), weight, 25);
#endif

  std::vector<vertex_descriptor> p(num_vertices(g));
  std::vector<int> d(num_vertices(g));
  vertex_descriptor s = vertex(N11, g);

#if defined(BOOST_MSVC) && BOOST_MSVC <= 1300 // for VC++
  property_map<graph_t, vertex_index_t>::type indexmap = get(vertex_index, g);
  dijkstra_shortest_paths(g, s, &p[0], &d[0], weightmap, indexmap, 
                          std::less<int>(), closed_plus<int>(), 
                          (std::numeric_limits<int>::max)(), 0,
                          default_dijkstra_visitor());
#else  // for regular compilers
  dijkstra_shortest_paths(g, s, predecessor_map(&p[0]).distance_map(&d[0]));
#endif
  graph_traits < graph_t >::vertex_iterator vi, vend;

  std::cout << "Backtracking the optimal path from the destination to source" << std::endl;
  for(int node = N55; node != N11; node = p[node]) {
    std::cout << "Path: N" << getNodeID(p[node]) << " -> N" << getNodeID(node) << ", Distance from origin is " << d[node] << std::endl;
  }
  return 0;
}

Office Hours

  • Friday 9:30AM-12:30PM

Information on Biostatistics Cluster

Standards of Academic Conduct

The following is an extract from the School of Public Health's Student Code of Conduct [1]:

Student academic misconduct includes behavior involving plagiarism, cheating, fabrication, falsification of records or official documents, intentional misuse of equipment or materials, and aiding and abetting the perpetration of such acts. The preparation of reports, papers, and examinations, assigned on an individual basis, must represent each student’s own effort. Reference sources should be indicated clearly. The use of assistance from other students or aids of any kind during a written examination, except when the use of books or notes has been approved by an instructor, is a violation of the standard of academic conduct.

In the context of this course, any work you hand-in should be your own.

Course History

Goncalo Abecasis taught it in several academic years previously. For previous course notes, see [Goncalo's older class notes].

Retrieved from "http://genome.sph.umich.edu/w/index.php?title=Biostatistics_615/815_Winter_2011&oldid=3648"