18 files changed, 4415 insertions, 0 deletions
diff --git a/gi/pyp-topics/src/Makefile b/gi/pyp-topics/src/Makefile
new file mode 100644
index 00000000..1d1391ae
--- /dev/null
+++ b/gi/pyp-topics/src/Makefile
@@ -0,0 +1,30 @@
+-include makefile.darwin
+
+local_objs = ../obj/corpus.o ../obj/gzstream.o ../obj/mt19937ar.o ../obj/pyp-topics.o ../obj/gammadist.o
+
+all: ../bin/pyp-topics-train
+
+-include makefile.depend
+
+#-----------------------#
+# Local stuff
+#-----------------------#
+
+../bin/pyp-topics-train: ../obj/train.o $(local_objs)
+	$(CXX) -o $@ $^ $(LDFLAGS)
+
+../obj/%.o: %.cc
+	${CXX} $(CXXFLAGS) -c $< -o $@
+
+../obj/%.o: %.c
+	${CC} $(CFLAGS) -c $< -o $@
+
+.PHONY: depend
+depend:
+	$(CXX) -MM $(CXXFLAGS) *.cc | sed 's/^\(.*\.o:\)/obj\/\1/' > makefile.depend
+
+clean:
+	rm -f ../obj/*.o
+
+#clobber: clean
+#	rm makefile.depend ../bin/${ARCH}/*
diff --git a/gi/pyp-topics/src/corpus.cc b/gi/pyp-topics/src/corpus.cc
new file mode 100644
index 00000000..93910ea3
--- /dev/null
+++ b/gi/pyp-topics/src/corpus.cc
@@ -0,0 +1,103 @@
+#include <sstream>
+#include <iostream>
+#include <set>
+
+#include "corpus.hh"
+#include "gzstream.hh"
+
+using namespace std;
+
+//////////////////////////////////////////////////
+// Corpus
+//////////////////////////////////////////////////
+
+Corpus::Corpus() {}
+
+unsigned Corpus::read(const std::string &filename) {
+  m_num_terms = 0;
+  m_num_types = 0;
+  std::set<int> seen_types;
+
+  igzstream in(filename.c_str());
+
+  string buf;
+  int token;
+  unsigned count=0;
+  while (getline(in, buf)) {
+    Document* doc(new Document());
+    istringstream ss(buf);
+
+    ss >> token; // the number of unique terms
+
+    char delimeter;
+    int count;
+    while(ss >> token >> delimeter >> count) {
+      for (int i=0; i<count; ++i)
+        doc->push_back(token);
+      m_num_terms += count;
+      seen_types.insert(token);
+    }
+
+    m_documents.push_back(doc);
+    count++;
+  }
+
+  m_num_types = seen_types.size();
+
+  return count;
+}
+
+//////////////////////////////////////////////////
+// TestCorpus
+//////////////////////////////////////////////////
+
+TestCorpus::TestCorpus() {}
+
+void TestCorpus::read(const std::string &filename) {
+  igzstream in(filename.c_str());
+
+  string buf;
+  Term term;
+  DocumentId doc;
+  char delimeter;
+  while (getline(in, buf)) {
+    DocumentTerms* line(new DocumentTerms());
+    istringstream ss(buf);
+
+    while(ss >> doc >> delimeter >> term)
+      line->push_back(DocumentTerm(doc, term));
+
+    m_lines.push_back(line);
+  }
+}
+
+//////////////////////////////////////////////////
+// TermBackoff
+//////////////////////////////////////////////////
+
+void TermBackoff::read(const std::string &filename) {
+  igzstream in(filename.c_str());
+
+  string buf;
+  int num_terms;
+  getline(in, buf);
+  istringstream ss(buf); 
+  ss >> num_terms >> m_backoff_order;
+
+  m_dict.resize(num_terms, -1);
+  for (int i=0; i<m_backoff_order; ++i) {
+    int count; ss >> count;
+    m_terms_at_order.push_back(count);
+  }
+
+  Term term, backoff;
+  while (getline(in, buf)) {
+    istringstream ss(buf);
+    ss >> term >> backoff;
+
+    assert(term < num_terms);
+    assert(term >= 0);
+
+    m_dict[term] = backoff;
+  }
+}
diff --git a/gi/pyp-topics/src/corpus.hh b/gi/pyp-topics/src/corpus.hh
new file mode 100644
index 00000000..3dd17cf9
--- /dev/null
+++ b/gi/pyp-topics/src/corpus.hh
@@ -0,0 +1,101 @@
+#ifndef _CORPUS_HH
+#define _CORPUS_HH
+
+#include <vector>
+#include <string>
+#include <map>
+
+#include <boost/ptr_container/ptr_vector.hpp>
+
+////////////////////////////////////////////////////////////////
+// Corpus
+////////////////////////////////////////////////////////////////
+typedef int Term;
+
+typedef std::vector<Term> Document;
+typedef std::vector<Term> Terms;
+
+class Corpus {
+public:
+    typedef boost::ptr_vector<Document>::const_iterator const_iterator;
+
+public:
+    Corpus();
+    ~Corpus() {}
+
+    unsigned read(const std::string &filename);
+
+    const_iterator begin() const { return m_documents.begin(); }
+    const_iterator end() const { return m_documents.end(); }
+
+    int num_documents() const { return m_documents.size(); }
+    int num_terms() const { return m_num_terms; }
+    int num_types() const { return m_num_types; }
+
+protected:
+    int m_num_terms, m_num_types;
+    boost::ptr_vector<Document> m_documents; 
+};
+
+typedef int DocumentId;
+struct DocumentTerm {
+  DocumentTerm(DocumentId d, Term t) : term(t), doc(d) {}
+  Term term;
+  DocumentId doc;
+};
+typedef std::vector<DocumentTerm> DocumentTerms;
+
+class TestCorpus {
+public:
+    typedef boost::ptr_vector<DocumentTerms>::const_iterator const_iterator;
+
+public:
+    TestCorpus();
+    ~TestCorpus() {}
+
+    void read(const std::string &filename);
+
+    const_iterator begin() const { return m_lines.begin(); }
+    const_iterator end() const { return m_lines.end(); }
+
+    int num_instances() const { return m_lines.size(); }
+
+protected:
+    boost::ptr_vector<DocumentTerms> m_lines; 
+};
+
+class TermBackoff {
+public:
+    typedef std::vector<Term> dictionary_type;
+    typedef dictionary_type::const_iterator const_iterator;
+
+public:
+    TermBackoff() : m_backoff_order(-1) {}
+    ~TermBackoff() {}
+
+    void read(const std::string &filename);
+
+    const_iterator begin() const { return m_dict.begin(); }
+    const_iterator end() const { return m_dict.end(); }
+
+    const Term& operator[](const Term& t) const {
+      assert(t < static_cast<int>(m_dict.size()));
+      return m_dict[t];
+    }
+
+    int order() const { return m_backoff_order; }
+//    int levels() const { return m_terms_at_order.size(); }
+    bool is_null(const Term& term) const { return term < 0; }
+    int terms_at_level(int level) const { 
+      assert (level < (int)m_terms_at_order.size());
+      return m_terms_at_order[level];
+    }
+
+    int size() const { return m_dict.size(); }
+
+protected:
+    dictionary_type m_dict; 
+    int m_backoff_order;
+    std::vector<int> m_terms_at_order;
+};
+#endif // _CORPUS_HH
diff --git a/gi/pyp-topics/src/gammadist.c b/gi/pyp-topics/src/gammadist.c
new file mode 100644
index 00000000..4e260db8
--- /dev/null
+++ b/gi/pyp-topics/src/gammadist.c
@@ -0,0 +1,247 @@
+/* gammadist.c -- computes probability of samples under / produces samples from a Gamma distribution
+ *
+ * Mark Johnson, 22nd March 2008
+ *
+ * WARNING: you need to set the flag -std=c99 to compile
+ *
+ * gammavariate() was translated from random.py in Python library
+ *
+ * The Gamma distribution is:
+ *
+ *   Gamma(x | alpha, beta) = pow(x/beta, alpha-1) * exp(-x/beta) / (gamma(alpha)*beta)
+ *
+ * shape parameter alpha > 0 (also called c), scale parameter beta > 0 (also called s); 
+ * mean is alpha*beta, variance is alpha*beta**2
+ *
+ * Note that many parameterizations of the Gamma function are in terms of an _inverse_
+ * scale parameter beta, which is the inverse of the beta given here.
+ *
+ * To define a main() that tests the routines, uncomment the following #define:
+ */
+/* #define GAMMATEST */
+
+#include <assert.h>
+#include <math.h> 
+
+#include "gammadist.h"
+#include "mt19937ar.h"
+
+/* gammadist() returns the probability density of x under a Gamma(alpha,beta) 
+ * distribution
+ */
+
+long double gammadist(long double x, long double alpha, long double beta) {
+  assert(alpha > 0);
+  assert(beta > 0);
+  return  pow(x/beta, alpha-1) * exp(-x/beta) / (tgamma(alpha)*beta);
+}
+
+/* lgammadist() returns the log probability density of x under a Gamma(alpha,beta)
+ * distribution
+ */
+
+long double lgammadist(long double x, long double alpha, long double beta) {
+  assert(alpha > 0);
+  assert(beta > 0);
+  return (alpha-1)*log(x) - alpha*log(beta) - x/beta - lgamma(alpha);
+}
+
+/* This definition of gammavariate is from Python code in
+ * the Python random module.
+ */
+
+long double gammavariate(long double alpha, long double beta) {
+
+  assert(alpha > 0);
+  assert(beta > 0);
+
+  if (alpha > 1.0) {
+    
+    /* Uses R.C.H. Cheng, "The generation of Gamma variables with
+       non-integral shape parameters", Applied Statistics, (1977), 26,
+       No. 1, p71-74 */
+
+    long double ainv = sqrt(2.0 * alpha - 1.0);
+    long double bbb = alpha - log(4.0);
+    long double ccc = alpha + ainv;
+    
+    while (1) {
+      long double u1 = mt_genrand_real3();
+      if (u1 > 1e-7  || u1 < 0.9999999) {
+	long double u2 = 1.0 - mt_genrand_real3();
+	long double v = log(u1/(1.0-u1))/ainv;
+	long double x = alpha*exp(v);
+	long double z = u1*u1*u2;
+	long double r = bbb+ccc*v-x;
+	if (r + (1.0+log(4.5)) - 4.5*z >= 0.0 || r >= log(z))
+	  return x * beta;
+      }
+    }
+  }
+  else if (alpha == 1.0) {
+    long double u = mt_genrand_real3();
+    while (u <= 1e-7)
+      u = mt_genrand_real3();
+    return -log(u) * beta;
+  }
+  else { 
+    /* alpha is between 0 and 1 (exclusive) 
+       Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle */
+    
+    while (1) {
+      long double u = mt_genrand_real3();
+      long double b = (exp(1) + alpha)/exp(1);
+      long double p = b*u;
+      long double x = (p <= 1.0) ? pow(p, 1.0/alpha) : -log((b-p)/alpha);
+      long double u1 = mt_genrand_real3();
+      if (! (((p <= 1.0) && (u1 > exp(-x))) ||
+	     ((p > 1.0)  &&  (u1 > pow(x, alpha - 1.0)))))
+	return x * beta;
+    }
+  }
+}
+
+/* betadist() returns the probability density of x under a Beta(alpha,beta)
+ * distribution.
+ */
+
+long double betadist(long double x, long double alpha, long double beta) {
+  assert(x >= 0);
+  assert(x <= 1);
+  assert(alpha > 0);
+  assert(beta > 0);
+  return pow(x,alpha-1)*pow(1-x,beta-1)*tgamma(alpha+beta)/(tgamma(alpha)*tgamma(beta));
+}
+
+/* lbetadist() returns the log probability density of x under a Beta(alpha,beta)
+ * distribution.
+ */
+
+long double lbetadist(long double x, long double alpha, long double beta) {
+  assert(x > 0);
+  assert(x < 1);
+  assert(alpha > 0);
+  assert(beta > 0);
+  return (alpha-1)*log(x)+(beta-1)*log(1-x)+lgamma(alpha+beta)-lgamma(alpha)-lgamma(beta);
+}
+
+/* betavariate() generates a sample from a Beta distribution with
+ * parameters alpha and beta.
+ *
+ * 0 < alpha < 1, 0 < beta < 1, mean is alpha/(alpha+beta)
+ */
+
+long double betavariate(long double alpha, long double beta) {
+  long double x = gammavariate(alpha, 1);
+  long double y = gammavariate(beta, 1);
+  return x/(x+y);
+}
+
+#ifdef GAMMATEST
+#include <stdio.h>
+
+int main(int argc, char **argv) {
+  int iteration, niterations = 1000;
+
+  for (iteration = 0; iteration < niterations; ++iteration) {
+    long double alpha = 100*mt_genrand_real3();
+    long double gv = gammavariate(alpha, 1);
+    long double pgv = gammadist(gv, alpha, 1);
+    long double pgvl = exp(lgammadist(gv, alpha, 1));
+    fprintf(stderr, "iteration = %d, gammavariate(%lg,1) = %lg, gammadist(%lg,%lg,1) = %lg, exp(lgammadist(%lg,%lg,1) = %lg\n",
+	    iteration, alpha, gv, gv, alpha, pgv, gv, alpha, pgvl);
+  }
+  return 0;
+}
+
+#endif /* GAMMATEST */
+
+
+/* Other routines I tried, but which weren't as good as the ones above */
+
+#if 0
+
+/*! gammavariate() returns samples from a Gamma distribution
+ *! where alpha is the shape parameter and beta is the scale 
+ *! parameter, using the algorithm described on p. 94 of 
+ *! Gentle (1998) Random Number Generation and Monte Carlo Methods, 
+ *! Springer.
+ */
+
+long double gammavariate(long double alpha) {
+
+  assert(alpha > 0); 
+  
+  if (alpha > 1.0) {
+    while (1) {
+      long double u1 = mt_genrand_real3();
+      long double u2 = mt_genrand_real3();
+      long double v = (alpha - 1/(6*alpha))*u1/(alpha-1)*u2;
+      if (2*(u2-1)/(alpha-1) + v + 1/v <= 2 
+         || 2*log(u2)/(alpha-1) - log(v) + v <= 1)
+	return (alpha-1)*v;
+    }
+  } else if (alpha < 1.0) {  
+    while (1) {
+      long double t = 0.07 + 0.75*sqrt(1-alpha);
+      long double b = alpha + exp(-t)*alpha/t;
+      long double u1 = mt_genrand_real3();
+      long double u2 = mt_genrand_real3();
+      long double v = b*u1;
+      if (v <= 1) {
+	long double x = t*pow(v, 1/alpha);
+	if (u2 <= (2 - x)/(2 + x))
+	  return x;
+	if (u2 <= exp(-x))
+	  return x;
+      }
+      else {
+	long double x = log(t*(b-v)/alpha);
+	long double y = x/t;
+	if (u2*(alpha + y*(1-alpha)) <= 1)
+	  return x;
+	if (u2 <= pow(y,alpha-1))
+	  return x;
+      }
+    }
+  }
+  else  
+    return -log(mt_genrand_real3());
+} 
+
+
+/*! gammavariate() returns a deviate distributed as a gamma
+ *! distribution of order alpha, beta, i.e., a waiting time to the alpha'th
+ *! event in a Poisson process of unit mean.
+ *!
+ *! Code from Numerical Recipes
+ */
+
+long double nr_gammavariate(long double ia) {
+  int j;
+  long double am,e,s,v1,v2,x,y;
+  assert(ia > 0);
+  if (ia < 10) { 
+    x=1.0; 
+    for (j=1;j<=ia;j++) 
+      x *= mt_genrand_real3();
+    x = -log(x);
+  } else { 
+    do {
+      do {
+	do { 
+	  v1=mt_genrand_real3();
+	  v2=2.0*mt_genrand_real3()-1.0;
+	} while (v1*v1+v2*v2 > 1.0); 
+	y=v2/v1;
+	am=ia-1;
+	s=sqrt(2.0*am+1.0);
+	x=s*y+am;
+      } while (x <= 0.0);
+      e=(1.0+y*y)*exp(am*log(x/am)-s*y);
+    } while (mt_genrand_real3() > e);
+  }
+  return x;
+} 
+
+#endif
diff --git a/gi/pyp-topics/src/gammadist.h b/gi/pyp-topics/src/gammadist.h
new file mode 100644
index 00000000..b6ad6c40
--- /dev/null
+++ b/gi/pyp-topics/src/gammadist.h
@@ -0,0 +1,72 @@
+/* gammadist.h -- computes probability of samples under / produces samples from a Gamma distribution
+ *
+ * Mark Johnson, 22nd March 2008
+ *
+ * gammavariate() was translated from random.py in Python library
+ *
+ * The Gamma distribution is:
+ *
+ *   Gamma(x | alpha, beta) = pow(x/beta, alpha-1) * exp(-x/beta) / (gamma(alpha)*beta)
+ *
+ * shape parameter alpha > 0 (also called c), scale parameter beta > 0 (also called s); 
+ * mean is alpha*beta, variance is alpha*beta**2
+ *
+ * Note that many parameterizations of the Gamma function are in terms of an _inverse_
+ * scale parameter beta, which is the inverse of the beta given here.
+ */
+
+#ifndef GAMMADIST_H
+#define GAMMADIST_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+  
+  /* gammadist() returns the probability density of x under a Gamma(alpha,beta) 
+   * distribution
+   */
+
+  long double gammadist(long double x, long double alpha, long double beta);
+
+  /* lgammadist() returns the log probability density of x under a Gamma(alpha,beta)
+   * distribution
+   */
+
+  long double lgammadist(long double x, long double alpha, long double beta);
+
+  /* gammavariate() generates samples from a Gamma distribution
+   * conditioned on the parameters alpha and beta.
+   * 
+   * alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
+   *
+   * Warning: a few older sources define the gamma distribution in terms
+   * of alpha > -1.0
+   */
+
+  long double gammavariate(long double alpha, long double beta);
+
+  /* betadist() returns the probability density of x under a Beta(alpha,beta)
+   * distribution.
+   */
+
+  long double betadist(long double x, long double alpha, long double beta);
+
+  /* lbetadist() returns the log probability density of x under a Beta(alpha,beta)
+   * distribution.
+   */
+  
+  long double lbetadist(long double x, long double alpha, long double beta);
+
+  /* betavariate() generates a sample from a Beta distribution with
+   * parameters alpha and beta.
+   *
+   * 0 < alpha < 1, 0 < beta < 1, mean is alpha/(alpha+beta)
+   */
+
+  long double betavariate(long double alpha, long double beta);
+
+#ifdef __cplusplus
+};
+#endif
+
+#endif /* GAMMADIST_H */
diff --git a/gi/pyp-topics/src/gzstream.cc b/gi/pyp-topics/src/gzstream.cc
new file mode 100644
index 00000000..7c4d3a12
--- /dev/null
+++ b/gi/pyp-topics/src/gzstream.cc
@@ -0,0 +1,165 @@
+// ============================================================================
+// gzstream, C++ iostream classes wrapping the zlib compression library.
+// Copyright (C) 2001  Deepak Bandyopadhyay, Lutz Kettner
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+// ============================================================================
+//
+// File          : gzstream.C
+// Revision      : $Revision: 1.1 $
+// Revision_date : $Date: 2006/03/30 04:05:52 $
+// Author(s)     : Deepak Bandyopadhyay, Lutz Kettner
+// 
+// Standard streambuf implementation following Nicolai Josuttis, "The 
+// Standard C++ Library".
+// ============================================================================
+
+#include "gzstream.hh"
+#include <iostream>
+#include <string.h>  // for memcpy
+
+#ifdef GZSTREAM_NAMESPACE
+namespace GZSTREAM_NAMESPACE {
+#endif
+
+// ----------------------------------------------------------------------------
+// Internal classes to implement gzstream. See header file for user classes.
+// ----------------------------------------------------------------------------
+
+// --------------------------------------
+// class gzstreambuf:
+// --------------------------------------
+
+gzstreambuf* gzstreambuf::open( const char* name, int open_mode) {
+    if ( is_open())
+        return (gzstreambuf*)0;
+    mode = open_mode;
+    // no append nor read/write mode
+    if ((mode & std::ios::ate) || (mode & std::ios::app)
+        || ((mode & std::ios::in) && (mode & std::ios::out)))
+        return (gzstreambuf*)0;
+    char  fmode[10];
+    char* fmodeptr = fmode;
+    if ( mode & std::ios::in)
+        *fmodeptr++ = 'r';
+    else if ( mode & std::ios::out)
+        *fmodeptr++ = 'w';
+    *fmodeptr++ = 'b';
+    *fmodeptr = '\0';
+    file = gzopen( name, fmode);
+    if (file == 0)
+        return (gzstreambuf*)0;
+    opened = 1;
+    return this;
+}
+
+gzstreambuf * gzstreambuf::close() {
+    if ( is_open()) {
+        sync();
+        opened = 0;
+        if ( gzclose( file) == Z_OK)
+            return this;
+    }
+    return (gzstreambuf*)0;
+}
+
+int gzstreambuf::underflow() { // used for input buffer only
+    if ( gptr() && ( gptr() < egptr()))
+        return * reinterpret_cast<unsigned char *>( gptr());
+
+    if ( ! (mode & std::ios::in) || ! opened)
+        return EOF;
+    // Josuttis' implementation of inbuf
+    int n_putback = gptr() - eback();
+    if ( n_putback > 4)
+        n_putback = 4;
+    memcpy( buffer + (4 - n_putback), gptr() - n_putback, n_putback);
+
+    int num = gzread( file, buffer+4, bufferSize-4);
+    if (num <= 0) // ERROR or EOF
+        return EOF;
+
+    // reset buffer pointers
+    setg( buffer + (4 - n_putback),   // beginning of putback area
+          buffer + 4,                 // read position
+          buffer + 4 + num);          // end of buffer
+
+    // return next character
+    return * reinterpret_cast<unsigned char *>( gptr());    
+}
+
+int gzstreambuf::flush_buffer() {
+    // Separate the writing of the buffer from overflow() and
+    // sync() operation.
+    int w = pptr() - pbase();
+    if ( gzwrite( file, pbase(), w) != w)
+        return EOF;
+    pbump( -w);
+    return w;
+}
+
+int gzstreambuf::overflow( int c) { // used for output buffer only
+    if ( ! ( mode & std::ios::out) || ! opened)
+        return EOF;
+    if (c != EOF) {
+        *pptr() = c;
+        pbump(1);
+    }
+    if ( flush_buffer() == EOF)
+        return EOF;
+    return c;
+}
+
+int gzstreambuf::sync() {
+    // Changed to use flush_buffer() instead of overflow( EOF)
+    // which caused improper behavior with std::endl and flush(),
+    // bug reported by Vincent Ricard.
+    if ( pptr() && pptr() > pbase()) {
+        if ( flush_buffer() == EOF)
+            return -1;
+    }
+    return 0;
+}
+
+// --------------------------------------
+// class gzstreambase:
+// --------------------------------------
+
+gzstreambase::gzstreambase( const char* name, int mode) {
+    init( &buf);
+    open( name, mode);
+}
+
+gzstreambase::~gzstreambase() {
+    buf.close();
+}
+
+void gzstreambase::open( const char* name, int open_mode) {
+    if ( ! buf.open( name, open_mode))
+        clear( rdstate() | std::ios::badbit);
+}
+
+void gzstreambase::close() {
+    if ( buf.is_open())
+        if ( ! buf.close())
+            clear( rdstate() | std::ios::badbit);
+}
+
+#ifdef GZSTREAM_NAMESPACE
+} // namespace GZSTREAM_NAMESPACE
+#endif
+
+// ============================================================================
+// EOF //
diff --git a/gi/pyp-topics/src/gzstream.hh b/gi/pyp-topics/src/gzstream.hh
new file mode 100644
index 00000000..ad9785fd
--- /dev/null
+++ b/gi/pyp-topics/src/gzstream.hh
@@ -0,0 +1,121 @@
+// ============================================================================
+// gzstream, C++ iostream classes wrapping the zlib compression library.
+// Copyright (C) 2001  Deepak Bandyopadhyay, Lutz Kettner
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+// ============================================================================
+//
+// File          : gzstream.h
+// Revision      : $Revision: 1.1 $
+// Revision_date : $Date: 2006/03/30 04:05:52 $
+// Author(s)     : Deepak Bandyopadhyay, Lutz Kettner
+// 
+// Standard streambuf implementation following Nicolai Josuttis, "The 
+// Standard C++ Library".
+// ============================================================================
+
+#ifndef GZSTREAM_H
+#define GZSTREAM_H 1
+
+// standard C++ with new header file names and std:: namespace
+#include <iostream>
+#include <fstream>
+#include <zlib.h>
+
+#ifdef GZSTREAM_NAMESPACE
+namespace GZSTREAM_NAMESPACE {
+#endif
+
+// ----------------------------------------------------------------------------
+// Internal classes to implement gzstream. See below for user classes.
+// ----------------------------------------------------------------------------
+
+class gzstreambuf : public std::streambuf {
+private:
+    static const int bufferSize = 47+256;    // size of data buff
+    // totals 512 bytes under g++ for igzstream at the end.
+
+    gzFile           file;               // file handle for compressed file
+    char             buffer[bufferSize]; // data buffer
+    char             opened;             // open/close state of stream
+    int              mode;               // I/O mode
+
+    int flush_buffer();
+public:
+    gzstreambuf() : opened(0) {
+        setp( buffer, buffer + (bufferSize-1));
+        setg( buffer + 4,     // beginning of putback area
+              buffer + 4,     // read position
+              buffer + 4);    // end position      
+        // ASSERT: both input & output capabilities will not be used together
+    }
+    int is_open() { return opened; }
+    gzstreambuf* open( const char* name, int open_mode);
+    gzstreambuf* close();
+    ~gzstreambuf() { close(); }
+    
+    virtual int     overflow( int c = EOF);
+    virtual int     underflow();
+    virtual int     sync();
+};
+
+class gzstreambase : virtual public std::ios {
+protected:
+    gzstreambuf buf;
+public:
+    gzstreambase() { init(&buf); }
+    gzstreambase( const char* name, int open_mode);
+    ~gzstreambase();
+    void open( const char* name, int open_mode);
+    void close();
+    gzstreambuf* rdbuf() { return &buf; }
+};
+
+// ----------------------------------------------------------------------------
+// User classes. Use igzstream and ogzstream analogously to ifstream and
+// ofstream respectively. They read and write files based on the gz* 
+// function interface of the zlib. Files are compatible with gzip compression.
+// ----------------------------------------------------------------------------
+
+class igzstream : public gzstreambase, public std::istream {
+public:
+    igzstream() : std::istream( &buf) {} 
+    igzstream( const char* name, int open_mode = std::ios::in)
+        : gzstreambase( name, open_mode), std::istream( &buf) {}  
+    gzstreambuf* rdbuf() { return gzstreambase::rdbuf(); }
+    void open( const char* name, int open_mode = std::ios::in) {
+        gzstreambase::open( name, open_mode);
+    }
+};
+
+class ogzstream : public gzstreambase, public std::ostream {
+public:
+    ogzstream() : std::ostream( &buf) {}
+    ogzstream( const char* name, int mode = std::ios::out)
+        : gzstreambase( name, mode), std::ostream( &buf) {}  
+    gzstreambuf* rdbuf() { return gzstreambase::rdbuf(); }
+    void open( const char* name, int open_mode = std::ios::out) {
+        gzstreambase::open( name, open_mode);
+    }
+};
+
+#ifdef GZSTREAM_NAMESPACE
+} // namespace GZSTREAM_NAMESPACE
+#endif
+
+#endif // GZSTREAM_H
+// ============================================================================
+// EOF //
+
diff --git a/gi/pyp-topics/src/log_add.h b/gi/pyp-topics/src/log_add.h
new file mode 100644
index 00000000..e0620c5a
--- /dev/null
+++ b/gi/pyp-topics/src/log_add.h
@@ -0,0 +1,30 @@
+#ifndef log_add_hh
+#define log_add_hh
+
+#include <limits>
+#include <iostream>
+#include <cassert>
+#include <cmath>
+
+template <typename T>
+struct Log
+{
+    static T zero() { return -std::numeric_limits<T>::infinity(); } 
+
+    static T add(T l1, T l2)
+    {
+        if (l1 == zero()) return l2;
+        if (l1 > l2) 
+            return l1 + std::log(1 + exp(l2 - l1));
+        else
+            return l2 + std::log(1 + exp(l1 - l2));
+    }
+
+    static T subtract(T l1, T l2)
+    {
+        //std::assert(l1 >= l2);
+        return l1 + log(1 - exp(l2 - l1));
+    }
+};
+
+#endif
diff --git a/gi/pyp-topics/src/makefile.darwin b/gi/pyp-topics/src/makefile.darwin
new file mode 100644
index 00000000..af608fd8
--- /dev/null
+++ b/gi/pyp-topics/src/makefile.darwin
@@ -0,0 +1,15 @@
+CC = /usr/bin/gcc
+CXX = /usr/bin/g++
+LD = /usr/bin/g++
+FC=/usr/bin/g77
+
+ARCH=i686-m64
+CXXFLAGS = -m64 -Wall -I/Users/pblunsom/packages/include 
+CFLAGS = -m64 -Wall -I/Users/pblunsom/packages/include 
+FFLAGS = -m64 -Wall
+LDFLAGS = -L/Users/pblunsom/packages/lib -lboost_program_options -lm -lz
+
+FFLAGS += -g -O3 -funroll-loops #-pg
+CFLAGS += -g -O3 -funroll-loops #-pg
+CXXFLAGS += -g -O3 -funroll-loops #-pg
+LDFLAGS += -g -O3 -funroll-loops #-pg
diff --git a/gi/pyp-topics/src/makefile.depend b/gi/pyp-topics/src/makefile.depend
new file mode 100644
index 00000000..89bafcf6
--- /dev/null
+++ b/gi/pyp-topics/src/makefile.depend
@@ -0,0 +1,1250 @@
+obj/corpus.o: corpus.cc corpus.hh \
+  /Users/pblunsom/packages/include/boost/ptr_container/ptr_vector.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/ptr_sequence_adapter.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/reversible_ptr_container.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/throw_exception.hpp \
+  /Users/pblunsom/packages/include/boost/assert.hpp \
+  /Users/pblunsom/packages/include/boost/config.hpp \
+  /Users/pblunsom/packages/include/boost/config/user.hpp \
+  /Users/pblunsom/packages/include/boost/config/select_compiler_config.hpp \
+  /Users/pblunsom/packages/include/boost/config/compiler/gcc.hpp \
+  /Users/pblunsom/packages/include/boost/config/select_stdlib_config.hpp \
+  /Users/pblunsom/packages/include/boost/config/no_tr1/utility.hpp \
+  /Users/pblunsom/packages/include/boost/config/stdlib/libstdcpp3.hpp \
+  /Users/pblunsom/packages/include/boost/config/select_platform_config.hpp \
+  /Users/pblunsom/packages/include/boost/config/platform/macos.hpp \
+  /Users/pblunsom/packages/include/boost/config/posix_features.hpp \
+  /Users/pblunsom/packages/include/boost/config/suffix.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/scoped_deleter.hpp \
+  /Users/pblunsom/packages/include/boost/scoped_array.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/scoped_array.hpp \
+  /Users/pblunsom/packages/include/boost/checked_delete.hpp \
+  /Users/pblunsom/packages/include/boost/detail/workaround.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/detail/operator_bool.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/static_move_ptr.hpp \
+  /Users/pblunsom/packages/include/boost/compressed_pair.hpp \
+  /Users/pblunsom/packages/include/boost/detail/compressed_pair.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/remove_cv.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/broken_compiler_spec.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/lambda_support.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/lambda.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/ttp.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/msvc.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/gcc.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/workaround.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/ctps.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/cv_traits_impl.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/type_trait_def.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/template_arity_spec.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/int.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/int_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/adl_barrier.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/adl.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/intel.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/nttp_decl.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/nttp.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/integral_wrapper.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/integral_c_tag.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/static_constant.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/static_cast.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/cat.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/config/config.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/template_arity_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessor/params.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/preprocessor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comma_if.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/punctuation/comma_if.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/if.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/iif.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/bool.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/facilities/empty.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/punctuation/comma.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repeat.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/repeat.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/debug/error.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/detail/auto_rec.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple/eat.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/inc.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/inc.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/overload_resolution.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/type_trait_undef.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_empty.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_convertible.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/intrinsics.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/config.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/yes_no_type.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_array.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/bool_trait_def.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/integral_constant.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/bool.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/bool_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/integral_c.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/integral_c_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/bool_trait_undef.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/add_reference.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_reference.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/ice.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/ice_or.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/ice_and.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/ice_not.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/ice_eq.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_arithmetic.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_integral.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_float.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_void.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_abstract.hpp \
+  /Users/pblunsom/packages/include/boost/static_assert.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_class.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_union.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_same.hpp \
+  /Users/pblunsom/packages/include/boost/call_traits.hpp \
+  /Users/pblunsom/packages/include/boost/detail/call_traits.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_pointer.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_member_pointer.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_member_function_pointer.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/is_mem_fun_pointer_impl.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/default_deleter.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/if.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/value_wknd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/integral.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/eti.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/na_spec.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/lambda_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/void_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/na.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/na_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/lambda_arity_param.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/arity.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/dtp.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessor/enum.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessor/def_params_tail.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/limits/arity.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/and.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/bitand.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/identity.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/facilities/identity.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/empty.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/add.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/dec.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/while.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/fold_left.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/detail/fold_left.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/expr_iif.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/adt.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/detail/is_binary.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/detail/check.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/compl.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/fold_right.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/detail/fold_right.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/reverse.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/detail/while.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple/elem.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/sub.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/remove_bounds.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/is_convertible.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/and.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/use_preprocessed.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/nested_type_wknd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/include_preprocessed.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/compiler.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/stringize.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/and.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/identity.hpp \
+  /Users/pblunsom/packages/include/boost/utility/enable_if.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/move.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/exception.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/clone_allocator.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/nullable.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/eval_if.hpp \
+  /Users/pblunsom/packages/include/boost/range/functions.hpp \
+  /Users/pblunsom/packages/include/boost/range/begin.hpp \
+  /Users/pblunsom/packages/include/boost/range/config.hpp \
+  /Users/pblunsom/packages/include/boost/range/iterator.hpp \
+  /Users/pblunsom/packages/include/boost/range/mutable_iterator.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/iterator_traits.hpp \
+  /Users/pblunsom/packages/include/boost/detail/iterator.hpp \
+  /Users/pblunsom/packages/include/boost/range/const_iterator.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/remove_const.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_volatile.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_const.hpp \
+  /Users/pblunsom/packages/include/boost/range/end.hpp \
+  /Users/pblunsom/packages/include/boost/range/detail/implementation_help.hpp \
+  /Users/pblunsom/packages/include/boost/range/detail/common.hpp \
+  /Users/pblunsom/packages/include/boost/range/detail/sfinae.hpp \
+  /Users/pblunsom/packages/include/boost/range/size.hpp \
+  /Users/pblunsom/packages/include/boost/range/difference_type.hpp \
+  /Users/pblunsom/packages/include/boost/range/distance.hpp \
+  /Users/pblunsom/packages/include/boost/range/empty.hpp \
+  /Users/pblunsom/packages/include/boost/range/rbegin.hpp \
+  /Users/pblunsom/packages/include/boost/range/reverse_iterator.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/reverse_iterator.hpp \
+  /Users/pblunsom/packages/include/boost/iterator.hpp \
+  /Users/pblunsom/packages/include/boost/utility.hpp \
+  /Users/pblunsom/packages/include/boost/utility/addressof.hpp \
+  /Users/pblunsom/packages/include/boost/utility/base_from_member.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_binary_params.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple/rem.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_params.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/repeat_from_to.hpp \
+  /Users/pblunsom/packages/include/boost/utility/binary.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/deduce_d.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/cat.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/fold_left.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/seq.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/elem.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/size.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/transform.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/mod.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/detail/div_base.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison/less_equal.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/not.hpp \
+  /Users/pblunsom/packages/include/boost/next_prior.hpp \
+  /Users/pblunsom/packages/include/boost/noncopyable.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/iterator_adaptor.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/iterator_categories.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/detail/config_def.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/placeholders.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/arg.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/arg_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/na_assert.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/assert.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/not.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/yes_no.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/arrays.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/pp_counter.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/arity_spec.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/arg_typedef.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/arg.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/placeholders.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/detail/config_undef.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/iterator_facade.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/interoperable.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/or.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/or.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/detail/facade_iterator_category.hpp \
+  /Users/pblunsom/packages/include/boost/detail/indirect_traits.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_function.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/false_result.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/is_function_ptr_helper.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/remove_reference.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/remove_pointer.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/detail/enable_if.hpp \
+  /Users/pblunsom/packages/include/boost/implicit_cast.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/add_const.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/add_pointer.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_pod.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_scalar.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_enum.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/always.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/apply.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/apply_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/apply_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/apply_wrap.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/has_apply.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/has_xxx.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/type_wrapper.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/has_xxx.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/msvc_typename.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/has_apply.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/msvc_never_true.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/apply_wrap.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/lambda.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/bind.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/bind_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/bind.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/bind_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/next.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/next_prior.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/common_name_wknd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/protect.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/bind.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/full_lambda.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/quote.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/void.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/has_type.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/bcc.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/quote.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/template_arity.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/template_arity.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/full_lambda.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/apply.hpp \
+  /Users/pblunsom/packages/include/boost/range/rend.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/indirect_fun.hpp \
+  /Users/pblunsom/packages/include/boost/utility/result_of.hpp \
+  /Users/pblunsom/packages/include/boost/type.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/library.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/div.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/mul.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/data.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/elem.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/size.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/insert.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/push_back.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison/not_equal.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/pop_back.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/deduce_z.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/pop_front.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/push_front.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/remove.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/replace.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/reverse.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple/reverse.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison/equal.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison/greater.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison/less.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison/greater_equal.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/config/limits.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/expr_if.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/debug.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/debug/assert.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/debug/line.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration/iterate.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/slot/slot.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/slot/detail/def.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/facilities.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/facilities/apply.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/detail/is_unary.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/facilities/expand.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/facilities/intercept.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration/local.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration/self.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/append.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/at.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/rest_n.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/cat.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/enum.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/for_each_i.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/for.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/detail/for.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/filter.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/first_n.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/for_each.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/for_each_product.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/to_tuple.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple/to_list.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/size.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/transform.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/bitnor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/bitor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/bitxor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/nor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/or.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/xor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/punctuation.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/punctuation/paren.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/punctuation/paren_if.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/deduce_r.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_params_with_a_default.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_params_with_defaults.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_shifted.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_shifted_binary_params.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_shifted_params.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_trailing.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_trailing_binary_params.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_trailing_params.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/selection.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/selection/max.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/selection/min.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/enum.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/filter.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/first_n.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/detail/split.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/fold_right.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/reverse.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/for_each.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/for_each_i.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/for_each_product.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/insert.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/rest_n.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/pop_back.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/pop_front.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/push_back.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/push_front.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/remove.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/replace.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/subseq.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/to_array.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/to_tuple.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/slot.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple/to_seq.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration/detail/iter/forward1.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration/detail/bounds/lower1.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/slot/detail/shared.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration/detail/bounds/upper1.hpp \
+  /Users/pblunsom/packages/include/boost/utility/detail/result_of_iterate.hpp \
+  /Users/pblunsom/packages/include/boost/pointee.hpp \
+  /Users/pblunsom/packages/include/boost/detail/is_incrementable.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/void_ptr_iterator.hpp \
+  gzstream.hh
+obj/gzstream.o: gzstream.cc gzstream.hh
+obj/pyp-topics.o: pyp-topics.cc pyp-topics.hh pyp.hh log_add.h gammadist.h \
+  slice-sampler.h mt19937ar.h corpus.hh \
+  /Users/pblunsom/packages/include/boost/ptr_container/ptr_vector.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/ptr_sequence_adapter.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/reversible_ptr_container.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/throw_exception.hpp \
+  /Users/pblunsom/packages/include/boost/assert.hpp \
+  /Users/pblunsom/packages/include/boost/config.hpp \
+  /Users/pblunsom/packages/include/boost/config/user.hpp \
+  /Users/pblunsom/packages/include/boost/config/select_compiler_config.hpp \
+  /Users/pblunsom/packages/include/boost/config/compiler/gcc.hpp \
+  /Users/pblunsom/packages/include/boost/config/select_stdlib_config.hpp \
+  /Users/pblunsom/packages/include/boost/config/no_tr1/utility.hpp \
+  /Users/pblunsom/packages/include/boost/config/stdlib/libstdcpp3.hpp \
+  /Users/pblunsom/packages/include/boost/config/select_platform_config.hpp \
+  /Users/pblunsom/packages/include/boost/config/platform/macos.hpp \
+  /Users/pblunsom/packages/include/boost/config/posix_features.hpp \
+  /Users/pblunsom/packages/include/boost/config/suffix.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/scoped_deleter.hpp \
+  /Users/pblunsom/packages/include/boost/scoped_array.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/scoped_array.hpp \
+  /Users/pblunsom/packages/include/boost/checked_delete.hpp \
+  /Users/pblunsom/packages/include/boost/detail/workaround.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/detail/operator_bool.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/static_move_ptr.hpp \
+  /Users/pblunsom/packages/include/boost/compressed_pair.hpp \
+  /Users/pblunsom/packages/include/boost/detail/compressed_pair.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/remove_cv.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/broken_compiler_spec.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/lambda_support.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/lambda.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/ttp.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/msvc.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/gcc.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/workaround.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/ctps.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/cv_traits_impl.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/type_trait_def.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/template_arity_spec.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/int.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/int_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/adl_barrier.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/adl.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/intel.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/nttp_decl.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/nttp.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/integral_wrapper.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/integral_c_tag.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/static_constant.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/static_cast.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/cat.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/config/config.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/template_arity_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessor/params.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/preprocessor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comma_if.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/punctuation/comma_if.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/if.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/iif.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/bool.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/facilities/empty.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/punctuation/comma.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repeat.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/repeat.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/debug/error.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/detail/auto_rec.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple/eat.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/inc.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/inc.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/overload_resolution.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/type_trait_undef.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_empty.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_convertible.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/intrinsics.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/config.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/yes_no_type.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_array.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/bool_trait_def.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/integral_constant.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/bool.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/bool_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/integral_c.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/integral_c_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/bool_trait_undef.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/add_reference.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_reference.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/ice.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/ice_or.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/ice_and.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/ice_not.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/ice_eq.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_arithmetic.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_integral.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_float.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_void.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_abstract.hpp \
+  /Users/pblunsom/packages/include/boost/static_assert.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_class.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_union.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_same.hpp \
+  /Users/pblunsom/packages/include/boost/call_traits.hpp \
+  /Users/pblunsom/packages/include/boost/detail/call_traits.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_pointer.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_member_pointer.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_member_function_pointer.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/is_mem_fun_pointer_impl.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/default_deleter.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/if.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/value_wknd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/integral.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/eti.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/na_spec.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/lambda_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/void_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/na.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/na_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/lambda_arity_param.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/arity.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/dtp.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessor/enum.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessor/def_params_tail.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/limits/arity.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/and.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/bitand.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/identity.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/facilities/identity.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/empty.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/add.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/dec.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/while.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/fold_left.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/detail/fold_left.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/expr_iif.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/adt.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/detail/is_binary.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/detail/check.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/compl.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/fold_right.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/detail/fold_right.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/reverse.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/detail/while.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple/elem.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/sub.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/remove_bounds.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/is_convertible.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/and.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/use_preprocessed.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/nested_type_wknd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/include_preprocessed.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/compiler.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/stringize.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/and.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/identity.hpp \
+  /Users/pblunsom/packages/include/boost/utility/enable_if.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/move.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/exception.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/clone_allocator.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/nullable.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/eval_if.hpp \
+  /Users/pblunsom/packages/include/boost/range/functions.hpp \
+  /Users/pblunsom/packages/include/boost/range/begin.hpp \
+  /Users/pblunsom/packages/include/boost/range/config.hpp \
+  /Users/pblunsom/packages/include/boost/range/iterator.hpp \
+  /Users/pblunsom/packages/include/boost/range/mutable_iterator.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/iterator_traits.hpp \
+  /Users/pblunsom/packages/include/boost/detail/iterator.hpp \
+  /Users/pblunsom/packages/include/boost/range/const_iterator.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/remove_const.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_volatile.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_const.hpp \
+  /Users/pblunsom/packages/include/boost/range/end.hpp \
+  /Users/pblunsom/packages/include/boost/range/detail/implementation_help.hpp \
+  /Users/pblunsom/packages/include/boost/range/detail/common.hpp \
+  /Users/pblunsom/packages/include/boost/range/detail/sfinae.hpp \
+  /Users/pblunsom/packages/include/boost/range/size.hpp \
+  /Users/pblunsom/packages/include/boost/range/difference_type.hpp \
+  /Users/pblunsom/packages/include/boost/range/distance.hpp \
+  /Users/pblunsom/packages/include/boost/range/empty.hpp \
+  /Users/pblunsom/packages/include/boost/range/rbegin.hpp \
+  /Users/pblunsom/packages/include/boost/range/reverse_iterator.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/reverse_iterator.hpp \
+  /Users/pblunsom/packages/include/boost/iterator.hpp \
+  /Users/pblunsom/packages/include/boost/utility.hpp \
+  /Users/pblunsom/packages/include/boost/utility/addressof.hpp \
+  /Users/pblunsom/packages/include/boost/utility/base_from_member.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_binary_params.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple/rem.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_params.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/repeat_from_to.hpp \
+  /Users/pblunsom/packages/include/boost/utility/binary.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/deduce_d.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/cat.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/fold_left.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/seq.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/elem.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/size.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/transform.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/mod.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/detail/div_base.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison/less_equal.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/not.hpp \
+  /Users/pblunsom/packages/include/boost/next_prior.hpp \
+  /Users/pblunsom/packages/include/boost/noncopyable.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/iterator_adaptor.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/iterator_categories.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/detail/config_def.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/placeholders.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/arg.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/arg_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/na_assert.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/assert.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/not.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/yes_no.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/arrays.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/pp_counter.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/arity_spec.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/arg_typedef.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/arg.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/placeholders.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/detail/config_undef.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/iterator_facade.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/interoperable.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/or.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/or.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/detail/facade_iterator_category.hpp \
+  /Users/pblunsom/packages/include/boost/detail/indirect_traits.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_function.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/false_result.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/is_function_ptr_helper.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/remove_reference.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/remove_pointer.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/detail/enable_if.hpp \
+  /Users/pblunsom/packages/include/boost/implicit_cast.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/add_const.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/add_pointer.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_pod.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_scalar.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_enum.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/always.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/apply.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/apply_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/apply_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/apply_wrap.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/has_apply.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/has_xxx.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/type_wrapper.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/has_xxx.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/msvc_typename.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/has_apply.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/msvc_never_true.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/apply_wrap.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/lambda.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/bind.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/bind_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/bind.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/bind_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/next.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/next_prior.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/common_name_wknd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/protect.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/bind.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/full_lambda.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/quote.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/void.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/has_type.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/bcc.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/quote.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/template_arity.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/template_arity.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/full_lambda.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/apply.hpp \
+  /Users/pblunsom/packages/include/boost/range/rend.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/indirect_fun.hpp \
+  /Users/pblunsom/packages/include/boost/utility/result_of.hpp \
+  /Users/pblunsom/packages/include/boost/type.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/library.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/div.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/mul.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/data.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/elem.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/size.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/insert.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/push_back.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison/not_equal.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/pop_back.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/deduce_z.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/pop_front.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/push_front.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/remove.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/replace.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/reverse.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple/reverse.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison/equal.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison/greater.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison/less.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison/greater_equal.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/config/limits.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/expr_if.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/debug.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/debug/assert.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/debug/line.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration/iterate.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/slot/slot.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/slot/detail/def.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/facilities.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/facilities/apply.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/detail/is_unary.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/facilities/expand.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/facilities/intercept.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration/local.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration/self.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/append.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/at.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/rest_n.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/cat.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/enum.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/for_each_i.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/for.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/detail/for.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/filter.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/first_n.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/for_each.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/for_each_product.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/to_tuple.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple/to_list.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/size.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/transform.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/bitnor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/bitor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/bitxor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/nor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/or.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/xor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/punctuation.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/punctuation/paren.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/punctuation/paren_if.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/deduce_r.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_params_with_a_default.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_params_with_defaults.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_shifted.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_shifted_binary_params.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_shifted_params.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_trailing.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_trailing_binary_params.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_trailing_params.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/selection.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/selection/max.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/selection/min.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/enum.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/filter.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/first_n.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/detail/split.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/fold_right.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/reverse.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/for_each.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/for_each_i.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/for_each_product.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/insert.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/rest_n.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/pop_back.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/pop_front.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/push_back.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/push_front.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/remove.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/replace.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/subseq.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/to_array.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/to_tuple.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/slot.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple/to_seq.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration/detail/iter/forward1.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration/detail/bounds/lower1.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/slot/detail/shared.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration/detail/bounds/upper1.hpp \
+  /Users/pblunsom/packages/include/boost/utility/detail/result_of_iterate.hpp \
+  /Users/pblunsom/packages/include/boost/pointee.hpp \
+  /Users/pblunsom/packages/include/boost/detail/is_incrementable.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/void_ptr_iterator.hpp
+obj/train.o: train.cc \
+  /Users/pblunsom/packages/include/boost/program_options/parsers.hpp \
+  /Users/pblunsom/packages/include/boost/program_options/config.hpp \
+  /Users/pblunsom/packages/include/boost/config.hpp \
+  /Users/pblunsom/packages/include/boost/config/user.hpp \
+  /Users/pblunsom/packages/include/boost/config/select_compiler_config.hpp \
+  /Users/pblunsom/packages/include/boost/config/compiler/gcc.hpp \
+  /Users/pblunsom/packages/include/boost/config/select_stdlib_config.hpp \
+  /Users/pblunsom/packages/include/boost/config/no_tr1/utility.hpp \
+  /Users/pblunsom/packages/include/boost/config/stdlib/libstdcpp3.hpp \
+  /Users/pblunsom/packages/include/boost/config/select_platform_config.hpp \
+  /Users/pblunsom/packages/include/boost/config/platform/macos.hpp \
+  /Users/pblunsom/packages/include/boost/config/posix_features.hpp \
+  /Users/pblunsom/packages/include/boost/config/suffix.hpp \
+  /Users/pblunsom/packages/include/boost/version.hpp \
+  /Users/pblunsom/packages/include/boost/config/auto_link.hpp \
+  /Users/pblunsom/packages/include/boost/program_options/option.hpp \
+  /Users/pblunsom/packages/include/boost/program_options/detail/cmdline.hpp \
+  /Users/pblunsom/packages/include/boost/program_options/errors.hpp \
+  /Users/pblunsom/packages/include/boost/program_options/cmdline.hpp \
+  /Users/pblunsom/packages/include/boost/program_options/options_description.hpp \
+  /Users/pblunsom/packages/include/boost/program_options/value_semantic.hpp \
+  /Users/pblunsom/packages/include/boost/any.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/remove_reference.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/broken_compiler_spec.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/lambda_support.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/lambda.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/ttp.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/msvc.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/gcc.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/workaround.hpp \
+  /Users/pblunsom/packages/include/boost/detail/workaround.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/ctps.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/type_trait_def.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/template_arity_spec.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/int.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/int_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/adl_barrier.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/adl.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/intel.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/nttp_decl.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/nttp.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/integral_wrapper.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/integral_c_tag.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/static_constant.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/static_cast.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/cat.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/config/config.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/template_arity_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessor/params.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/preprocessor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comma_if.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/punctuation/comma_if.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/if.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/iif.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/bool.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/facilities/empty.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/punctuation/comma.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repeat.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/repeat.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/debug/error.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/detail/auto_rec.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple/eat.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/inc.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/inc.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/overload_resolution.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/type_trait_undef.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_reference.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/config.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/bool_trait_def.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/integral_constant.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/bool.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/bool_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/integral_c.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/integral_c_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/bool_trait_undef.hpp \
+  /Users/pblunsom/packages/include/boost/throw_exception.hpp \
+  /Users/pblunsom/packages/include/boost/exception/detail/attribute_noreturn.hpp \
+  /Users/pblunsom/packages/include/boost/exception/exception.hpp \
+  /Users/pblunsom/packages/include/boost/current_function.hpp \
+  /Users/pblunsom/packages/include/boost/static_assert.hpp \
+  /Users/pblunsom/packages/include/boost/function/function1.hpp \
+  /Users/pblunsom/packages/include/boost/function/detail/maybe_include.hpp \
+  /Users/pblunsom/packages/include/boost/function/function_template.hpp \
+  /Users/pblunsom/packages/include/boost/function/detail/prologue.hpp \
+  /Users/pblunsom/packages/include/boost/config/no_tr1/functional.hpp \
+  /Users/pblunsom/packages/include/boost/function/function_base.hpp \
+  /Users/pblunsom/packages/include/boost/assert.hpp \
+  /Users/pblunsom/packages/include/boost/integer.hpp \
+  /Users/pblunsom/packages/include/boost/integer_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/limits.hpp \
+  /Users/pblunsom/packages/include/boost/integer_traits.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/has_trivial_copy.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/intrinsics.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_volatile.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/cv_traits_impl.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_pod.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_void.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_scalar.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_arithmetic.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_integral.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_float.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/ice_or.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_enum.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/add_reference.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_convertible.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/yes_no_type.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_array.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/ice.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/ice_and.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/ice_not.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/ice_eq.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_abstract.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_class.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_union.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/remove_cv.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_function.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/false_result.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/is_function_ptr_helper.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_pointer.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_member_pointer.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_member_function_pointer.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/is_mem_fun_pointer_impl.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/has_trivial_destructor.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_const.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/composite_traits.hpp \
+  /Users/pblunsom/packages/include/boost/ref.hpp \
+  /Users/pblunsom/packages/include/boost/utility/addressof.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/if.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/value_wknd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/integral.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/eti.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/na_spec.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/lambda_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/void_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/na.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/na_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/lambda_arity_param.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/arity.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/dtp.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessor/enum.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessor/def_params_tail.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/limits/arity.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/and.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/bitand.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/identity.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/facilities/identity.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/empty.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/add.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/dec.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/while.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/fold_left.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/detail/fold_left.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/expr_iif.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/adt.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/detail/is_binary.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/detail/check.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/compl.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/fold_right.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/detail/fold_right.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/reverse.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/detail/while.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple/elem.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/sub.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/alignment_of.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/size_t_trait_def.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/size_t.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/size_t_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/detail/size_t_trait_undef.hpp \
+  /Users/pblunsom/packages/include/boost/utility/enable_if.hpp \
+  /Users/pblunsom/packages/include/boost/function_equal.hpp \
+  /Users/pblunsom/packages/include/boost/function/function_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mem_fn.hpp \
+  /Users/pblunsom/packages/include/boost/bind/mem_fn.hpp \
+  /Users/pblunsom/packages/include/boost/get_pointer.hpp \
+  /Users/pblunsom/packages/include/boost/config/no_tr1/memory.hpp \
+  /Users/pblunsom/packages/include/boost/bind/mem_fn_template.hpp \
+  /Users/pblunsom/packages/include/boost/bind/mem_fn_cc.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/enum.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple/rem.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/enum_params.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_params.hpp \
+  /Users/pblunsom/packages/include/boost/detail/no_exceptions_support.hpp \
+  /Users/pblunsom/packages/include/boost/lexical_cast.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/make_unsigned.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_signed.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_unsigned.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_same.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/add_const.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/add_volatile.hpp \
+  /Users/pblunsom/packages/include/boost/call_traits.hpp \
+  /Users/pblunsom/packages/include/boost/detail/call_traits.hpp \
+  /Users/pblunsom/packages/include/boost/detail/lcast_precision.hpp \
+  /Users/pblunsom/packages/include/boost/program_options/detail/value_semantic.hpp \
+  /Users/pblunsom/packages/include/boost/function.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iterate.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration/iterate.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/elem.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/data.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/size.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/slot/slot.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/slot/detail/def.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration/detail/iter/forward1.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration/detail/bounds/lower1.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/slot/detail/shared.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration/detail/bounds/upper1.hpp \
+  /Users/pblunsom/packages/include/boost/function/detail/function_iterate.hpp \
+  /Users/pblunsom/packages/include/boost/shared_ptr.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/shared_ptr.hpp \
+  /Users/pblunsom/packages/include/boost/checked_delete.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/detail/shared_count.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/bad_weak_ptr.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/detail/sp_counted_base.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/detail/sp_has_sync.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/detail/sp_counted_base_gcc_x86.hpp \
+  /Users/pblunsom/packages/include/boost/detail/sp_typeinfo.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/detail/sp_counted_impl.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/detail/sp_convertible.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/detail/spinlock_pool.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/detail/spinlock.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/detail/spinlock_sync.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/detail/yield_k.hpp \
+  /Users/pblunsom/packages/include/boost/memory_order.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/detail/operator_bool.hpp \
+  /Users/pblunsom/packages/include/boost/program_options/positional_options.hpp \
+  /Users/pblunsom/packages/include/boost/program_options/detail/parsers.hpp \
+  /Users/pblunsom/packages/include/boost/program_options/detail/convert.hpp \
+  /Users/pblunsom/packages/include/boost/program_options/variables_map.hpp \
+  /Users/pblunsom/packages/include/boost/scoped_ptr.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/scoped_ptr.hpp \
+  pyp-topics.hh pyp.hh log_add.h gammadist.h slice-sampler.h mt19937ar.h \
+  corpus.hh \
+  /Users/pblunsom/packages/include/boost/ptr_container/ptr_vector.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/ptr_sequence_adapter.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/reversible_ptr_container.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/throw_exception.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/scoped_deleter.hpp \
+  /Users/pblunsom/packages/include/boost/scoped_array.hpp \
+  /Users/pblunsom/packages/include/boost/smart_ptr/scoped_array.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/static_move_ptr.hpp \
+  /Users/pblunsom/packages/include/boost/compressed_pair.hpp \
+  /Users/pblunsom/packages/include/boost/detail/compressed_pair.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/is_empty.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/default_deleter.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/remove_bounds.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/is_convertible.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/and.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/use_preprocessed.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/nested_type_wknd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/include_preprocessed.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/compiler.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/stringize.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/and.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/identity.hpp \
+  /Users/pblunsom/packages/include/boost/utility/enable_if.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/move.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/exception.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/clone_allocator.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/nullable.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/eval_if.hpp \
+  /Users/pblunsom/packages/include/boost/range/functions.hpp \
+  /Users/pblunsom/packages/include/boost/range/begin.hpp \
+  /Users/pblunsom/packages/include/boost/range/config.hpp \
+  /Users/pblunsom/packages/include/boost/range/iterator.hpp \
+  /Users/pblunsom/packages/include/boost/range/mutable_iterator.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/iterator_traits.hpp \
+  /Users/pblunsom/packages/include/boost/detail/iterator.hpp \
+  /Users/pblunsom/packages/include/boost/range/const_iterator.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/remove_const.hpp \
+  /Users/pblunsom/packages/include/boost/range/end.hpp \
+  /Users/pblunsom/packages/include/boost/range/detail/implementation_help.hpp \
+  /Users/pblunsom/packages/include/boost/range/detail/common.hpp \
+  /Users/pblunsom/packages/include/boost/range/detail/sfinae.hpp \
+  /Users/pblunsom/packages/include/boost/range/size.hpp \
+  /Users/pblunsom/packages/include/boost/range/difference_type.hpp \
+  /Users/pblunsom/packages/include/boost/range/distance.hpp \
+  /Users/pblunsom/packages/include/boost/range/empty.hpp \
+  /Users/pblunsom/packages/include/boost/range/rbegin.hpp \
+  /Users/pblunsom/packages/include/boost/range/reverse_iterator.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/reverse_iterator.hpp \
+  /Users/pblunsom/packages/include/boost/iterator.hpp \
+  /Users/pblunsom/packages/include/boost/utility.hpp \
+  /Users/pblunsom/packages/include/boost/utility/base_from_member.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_binary_params.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/repeat_from_to.hpp \
+  /Users/pblunsom/packages/include/boost/utility/binary.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/deduce_d.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/cat.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/fold_left.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/seq.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/elem.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/size.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/transform.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/mod.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/detail/div_base.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison/less_equal.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/not.hpp \
+  /Users/pblunsom/packages/include/boost/next_prior.hpp \
+  /Users/pblunsom/packages/include/boost/noncopyable.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/iterator_adaptor.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/iterator_categories.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/detail/config_def.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/placeholders.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/arg.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/arg_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/na_assert.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/assert.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/not.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/yes_no.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/arrays.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/pp_counter.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/arity_spec.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/arg_typedef.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/arg.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/placeholders.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/detail/config_undef.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/iterator_facade.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/interoperable.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/or.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/or.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/detail/facade_iterator_category.hpp \
+  /Users/pblunsom/packages/include/boost/detail/indirect_traits.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/remove_pointer.hpp \
+  /Users/pblunsom/packages/include/boost/iterator/detail/enable_if.hpp \
+  /Users/pblunsom/packages/include/boost/implicit_cast.hpp \
+  /Users/pblunsom/packages/include/boost/type_traits/add_pointer.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/always.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/apply.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/apply_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/apply_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/apply_wrap.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/has_apply.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/has_xxx.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/type_wrapper.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/has_xxx.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/msvc_typename.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/has_apply.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/msvc_never_true.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/apply_wrap.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/lambda.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/bind.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/bind_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/bind.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/bind_fwd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/next.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/next_prior.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/common_name_wknd.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/protect.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/bind.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/full_lambda.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/quote.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/void.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/has_type.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/config/bcc.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/quote.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/template_arity.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/template_arity.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/full_lambda.hpp \
+  /Users/pblunsom/packages/include/boost/mpl/aux_/preprocessed/gcc/apply.hpp \
+  /Users/pblunsom/packages/include/boost/range/rend.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/indirect_fun.hpp \
+  /Users/pblunsom/packages/include/boost/utility/result_of.hpp \
+  /Users/pblunsom/packages/include/boost/type.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/library.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/div.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/arithmetic/mul.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/insert.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/push_back.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison/not_equal.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/pop_back.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/deduce_z.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/pop_front.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/push_front.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/remove.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/replace.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/array/reverse.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple/reverse.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison/equal.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison/greater.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison/less.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/comparison/greater_equal.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/config/limits.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/control/expr_if.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/debug.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/debug/assert.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/debug/line.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/facilities.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/facilities/apply.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/detail/is_unary.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/facilities/expand.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/facilities/intercept.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration/local.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/iteration/self.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/append.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/at.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/rest_n.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/cat.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/enum.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/for_each_i.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/for.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/detail/for.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/filter.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/first_n.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/for_each.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/for_each_product.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/to_tuple.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple/to_list.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/size.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/list/transform.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/bitnor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/bitor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/bitxor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/nor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/or.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/logical/xor.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/punctuation.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/punctuation/paren.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/punctuation/paren_if.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/deduce_r.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_params_with_a_default.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_params_with_defaults.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_shifted.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_shifted_binary_params.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_shifted_params.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_trailing.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_trailing_binary_params.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/repetition/enum_trailing_params.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/selection.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/selection/max.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/selection/min.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/enum.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/filter.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/first_n.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/detail/split.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/fold_right.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/reverse.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/for_each.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/for_each_i.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/for_each_product.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/insert.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/rest_n.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/pop_back.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/pop_front.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/push_back.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/push_front.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/remove.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/replace.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/subseq.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/to_array.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/seq/to_tuple.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/slot.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple.hpp \
+  /Users/pblunsom/packages/include/boost/preprocessor/tuple/to_seq.hpp \
+  /Users/pblunsom/packages/include/boost/utility/detail/result_of_iterate.hpp \
+  /Users/pblunsom/packages/include/boost/pointee.hpp \
+  /Users/pblunsom/packages/include/boost/detail/is_incrementable.hpp \
+  /Users/pblunsom/packages/include/boost/ptr_container/detail/void_ptr_iterator.hpp \
+  gzstream.hh
diff --git a/gi/pyp-topics/src/mt19937ar.c b/gi/pyp-topics/src/mt19937ar.c
new file mode 100644
index 00000000..6551ea39
--- /dev/null
+++ b/gi/pyp-topics/src/mt19937ar.c
@@ -0,0 +1,194 @@
+/* 
+   A C-program for MT19937, with initialization improved 2002/1/26.
+   Coded by Takuji Nishimura and Makoto Matsumoto.
+
+   Before using, initialize the state by using mt_init_genrand(seed)  
+   or mt_init_by_array(init_key, key_length).
+
+   Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
+   All rights reserved.                          
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions
+   are met:
+
+     1. Redistributions of source code must retain the above copyright
+        notice, this list of conditions and the following disclaimer.
+
+     2. Redistributions in binary form must reproduce the above copyright
+        notice, this list of conditions and the following disclaimer in the
+        documentation and/or other materials provided with the distribution.
+
+     3. The names of its contributors may not be used to endorse or promote 
+        products derived from this software without specific prior written 
+        permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+   Any feedback is very welcome.
+   http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
+   email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
+*/
+
+#include "mt19937ar.h"   /* XXX MJ 17th March 2006 */
+
+/* Period parameters */  
+#define N 624
+#define M 397
+#define MATRIX_A 0x9908b0dfUL   /* constant vector a */
+#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
+#define LOWER_MASK 0x7fffffffUL /* least significant r bits */
+
+static unsigned long mt[N]; /* the array for the state vector  */
+static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */
+
+/* initializes mt[N] with a seed */
+void mt_init_genrand(unsigned long s)
+{
+    mt[0]= s & 0xffffffffUL;
+    for (mti=1; mti<N; mti++) {
+        mt[mti] = 
+	    (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti); 
+        /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
+        /* In the previous versions, MSBs of the seed affect   */
+        /* only MSBs of the array mt[].                        */
+        /* 2002/01/09 modified by Makoto Matsumoto             */
+        mt[mti] &= 0xffffffffUL;
+        /* for >32 bit machines */
+    }
+}
+
+/* initialize by an array with array-length */
+/* init_key is the array for initializing keys */
+/* key_length is its length */
+/* slight change for C++, 2004/2/26 */
+void mt_init_by_array(unsigned long init_key[], int key_length)
+{
+    int i, j, k;
+    mt_init_genrand(19650218UL);
+    i=1; j=0;
+    k = (N>key_length ? N : key_length);
+    for (; k; k--) {
+        mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
+          + init_key[j] + j; /* non linear */
+        mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
+        i++; j++;
+        if (i>=N) { mt[0] = mt[N-1]; i=1; }
+        if (j>=key_length) j=0;
+    }
+    for (k=N-1; k; k--) {
+        mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
+          - i; /* non linear */
+        mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
+        i++;
+        if (i>=N) { mt[0] = mt[N-1]; i=1; }
+    }
+
+    mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ 
+}
+
+/* generates a random number on [0,0xffffffff]-interval */
+unsigned long mt_genrand_int32(void)
+{
+    unsigned long y;
+    static unsigned long mag01[2]={0x0UL, MATRIX_A};
+    /* mag01[x] = x * MATRIX_A  for x=0,1 */
+
+    if (mti >= N) { /* generate N words at one time */
+        int kk;
+
+        if (mti == N+1)   /* if mt_init_genrand() has not been called, */
+            mt_init_genrand(5489UL); /* a default initial seed is used */
+
+        for (kk=0;kk<N-M;kk++) {
+            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
+            mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
+        }
+        for (;kk<N-1;kk++) {
+            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
+            mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
+        }
+        y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
+        mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
+
+        mti = 0;
+    }
+  
+    y = mt[mti++];
+
+    /* Tempering */
+    y ^= (y >> 11);
+    y ^= (y << 7) & 0x9d2c5680UL;
+    y ^= (y << 15) & 0xefc60000UL;
+    y ^= (y >> 18);
+
+    return y;
+}
+
+/* generates a random number on [0,0x7fffffff]-interval */
+long mt_genrand_int31(void)
+{
+    return (long)( mt_genrand_int32()>>1);
+}
+
+/* generates a random number on [0,1]-real-interval */
+double mt_genrand_real1(void)
+{
+    return mt_genrand_int32()*(1.0/4294967295.0); 
+    /* divided by 2^32-1 */ 
+}
+
+/* generates a random number on [0,1)-real-interval */
+double mt_genrand_real2(void)
+{
+    return mt_genrand_int32()*(1.0/4294967296.0); 
+    /* divided by 2^32 */
+}
+
+/* generates a random number on (0,1)-real-interval */
+double mt_genrand_real3(void)
+{
+    return (((double) mt_genrand_int32()) + 0.5)*(1.0/4294967296.0); 
+    /* divided by 2^32 */
+}
+
+/* generates a random number on [0,1) with 53-bit resolution*/
+double mt_genrand_res53(void) 
+{ 
+    unsigned long a=mt_genrand_int32()>>5, b=mt_genrand_int32()>>6; 
+    return(a*67108864.0+b)*(1.0/9007199254740992.0); 
+} 
+/* These real versions are due to Isaku Wada, 2002/01/09 added */
+
+/*
+#include <stdio.h>
+
+int main(void)
+{
+    int i;
+    unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4;
+    mt_init_by_array(init, length);
+    printf("1000 outputs of genrand_int32()\n");
+    for (i=0; i<1000; i++) {
+      printf("%10lu ", mt_genrand_int32());
+      if (i%5==4) printf("\n");
+    }
+    printf("\n1000 outputs of genrand_real2()\n");
+    for (i=0; i<1000; i++) {
+      printf("%10.8f ", mt_genrand_real2());
+      if (i%5==4) printf("\n");
+    }
+    return 0;
+}
+*/
diff --git a/gi/pyp-topics/src/mt19937ar.h b/gi/pyp-topics/src/mt19937ar.h
new file mode 100644
index 00000000..caab4045
--- /dev/null
+++ b/gi/pyp-topics/src/mt19937ar.h
@@ -0,0 +1,44 @@
+/* mt19937ar.h
+ *
+ * Mark Johnson, 17th March 2006
+ */
+
+#ifndef MT19937AR_H
+#define MT19937AR_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+  /* initializes mt[N] with a seed */
+  void mt_init_genrand(unsigned long s);
+
+  /* initialize by an array with array-length */
+  /* init_key is the array for initializing keys */
+  /* key_length is its length */
+  /* slight change for C++, 2004/2/26 */
+  void mt_init_by_array(unsigned long init_key[], int key_length);
+
+  /* generates a random number on [0,0xffffffff]-interval */
+  unsigned long mt_genrand_int32(void);
+
+  /* generates a random number on [0,0x7fffffff]-interval */
+  long mt_genrand_int31(void);
+
+  /* generates a random number on [0,1]-real-interval */
+  double mt_genrand_real1(void);
+
+  /* generates a random number on [0,1)-real-interval */
+  double mt_genrand_real2(void);
+
+  /* generates a random number on (0,1)-real-interval */
+  double mt_genrand_real3(void);
+
+  /* generates a random number on [0,1) with 53-bit resolution*/
+  double mt_genrand_res53(void);
+
+#ifdef __cplusplus
+};
+#endif
+
+#endif /* MT19937AR_H */
diff --git a/gi/pyp-topics/src/pyp-topics.cc b/gi/pyp-topics/src/pyp-topics.cc
new file mode 100644
index 00000000..9b43d6d1
--- /dev/null
+++ b/gi/pyp-topics/src/pyp-topics.cc
@@ -0,0 +1,224 @@
+#include "pyp-topics.hh"
+
+void PYPTopics::sample(const Corpus& corpus, int samples) {
+  if (!m_backoff.get()) {
+    m_word_pyps.clear();
+    m_word_pyps.push_back(PYPs());
+  }
+
+  std::cerr << " Training with " << m_word_pyps.size()-1 << " backoff level" 
+    << (m_word_pyps.size()==2 ? ":" : "s:") << std::endl;
+
+  for (int i=0; i<(int)m_word_pyps.size(); ++i)
+    m_word_pyps.at(i).resize(m_num_topics, PYP<int>(0.5, 1.0));
+  std::cerr << std::endl;
+
+  m_document_pyps.resize(corpus.num_documents(), PYP<int>(0.5, 1.0));
+
+  m_topic_p0 = 1.0/m_num_topics;
+  m_term_p0 = 1.0/corpus.num_types();
+  m_backoff_p0 = 1.0/corpus.num_documents();
+
+  std::cerr << " Documents: " << corpus.num_documents() << " Terms: " 
+    << corpus.num_types() << std::endl;
+
+  // Initialisation pass
+  int document_id=0, topic_counter=0;
+  for (Corpus::const_iterator corpusIt=corpus.begin(); 
+       corpusIt != corpus.end(); ++corpusIt, ++document_id) {
+    m_corpus_topics.push_back(DocumentTopics(corpusIt->size(), 0));
+
+    int term_index=0;
+    for (Document::const_iterator docIt=corpusIt->begin();
+         docIt != corpusIt->end(); ++docIt, ++term_index) {
+      topic_counter++;
+      Term term = *docIt;
+
+      // sample a new_topic
+      //int new_topic = (topic_counter % m_num_topics);
+      int new_topic = (document_id % m_num_topics);
+
+      // add the new topic to the PYPs
+      m_corpus_topics[document_id][term_index] = new_topic;
+      increment(term, new_topic);
+      m_document_pyps[document_id].increment(new_topic, m_topic_p0);
+    }
+  }
+
+  // Sampling phase
+  for (int curr_sample=0; curr_sample < samples; ++curr_sample) {
+    std::cerr << "\n  -- Sample " << curr_sample << " "; std::cerr.flush();
+
+    // for each document in the corpus
+    int document_id=0;
+    for (Corpus::const_iterator corpusIt=corpus.begin(); 
+         corpusIt != corpus.end(); ++corpusIt, ++document_id) {
+      // for each term in the document
+      int term_index=0;
+      for (Document::const_iterator docIt=corpusIt->begin();
+           docIt != corpusIt->end(); ++docIt, ++term_index) {
+        Term term = *docIt;
+
+        // remove the prevous topic from the PYPs
+        int current_topic = m_corpus_topics[document_id][term_index];
+        decrement(term, current_topic);
+        m_document_pyps[document_id].decrement(current_topic);
+
+        // sample a new_topic
+        int new_topic = sample(document_id, term);
+
+        // add the new topic to the PYPs
+        m_corpus_topics[document_id][term_index] = new_topic;
+        increment(term, new_topic);
+        m_document_pyps[document_id].increment(new_topic, m_topic_p0);
+      }
+      if (document_id && document_id % 10000 == 0) {
+        std::cerr << "."; std::cerr.flush();
+      }
+    }
+
+    if (curr_sample != 0 && curr_sample % 10 == 0) {
+      std::cerr << " ||| Resampling hyperparameters "; std::cerr.flush();
+      // resample the hyperparamters
+      F log_p=0.0; int resample_counter=0;
+      for (std::vector<PYPs>::iterator levelIt=m_word_pyps.begin();
+           levelIt != m_word_pyps.end(); ++levelIt) {
+        for (PYPs::iterator pypIt=levelIt->begin();
+             pypIt != levelIt->end(); ++pypIt) {
+          pypIt->resample_prior();
+          log_p += pypIt->log_restaurant_prob();
+          if (resample_counter++ % 100 == 0) {
+            std::cerr << "."; std::cerr.flush();
+          }
+        }
+      }
+
+      for (PYPs::iterator pypIt=m_document_pyps.begin();
+           pypIt != m_document_pyps.end(); ++pypIt) {
+        pypIt->resample_prior();
+        log_p += pypIt->log_restaurant_prob();
+      }
+      std::cerr << " ||| LLH=" << log_p << std::endl;
+    }
+  }
+}
+
+void PYPTopics::decrement(const Term& term, int topic, int level) {
+  //std::cerr << "PYPTopics::decrement(" << term << "," << topic << "," << level << ")" << std::endl;
+  m_word_pyps.at(level).at(topic).decrement(term);
+  if (m_backoff.get()) {
+    Term backoff_term = (*m_backoff)[term];
+    if (!m_backoff->is_null(backoff_term))
+      decrement(backoff_term, topic, level+1);
+  }
+}
+
+void PYPTopics::increment(const Term& term, int topic, int level) {
+  //std::cerr << "PYPTopics::increment(" << term << "," << topic << "," << level << ")" << std::endl;
+  m_word_pyps.at(level).at(topic).increment(term, word_pyps_p0(term, topic, level));
+
+  if (m_backoff.get()) {
+    Term backoff_term = (*m_backoff)[term];
+    if (!m_backoff->is_null(backoff_term))
+      increment(backoff_term, topic, level+1);
+  }
+}
+
+int PYPTopics::sample(const DocumentId& doc, const Term& term) {
+  // First pass: collect probs
+  F sum=0.0;
+  std::vector<F> sums;
+  for (int k=0; k<m_num_topics; ++k) {
+    F p_w_k = prob(term, k);
+    F p_k_d = m_document_pyps[doc].prob(k, m_topic_p0);
+    sum += (p_w_k*p_k_d);
+    sums.push_back(sum);
+  }
+  // Second pass: sample a topic
+  F cutoff = mt_genrand_res53() * sum;
+  for (int k=0; k<m_num_topics; ++k) {
+    if (cutoff <= sums[k]) 
+      return k;
+  }
+  assert(false);
+}
+
+PYPTopics::F PYPTopics::word_pyps_p0(const Term& term, int topic, int level) const {
+  //for (int i=0; i<level+1; ++i) std::cerr << "  ";
+  //std::cerr << "PYPTopics::word_pyps_p0(" << term << "," << topic << "," << level << ")" << std::endl;
+
+  F p0 = m_term_p0;
+  if (m_backoff.get()) {
+    //static F fudge=m_backoff_p0; // TODO
+
+    Term backoff_term = (*m_backoff)[term];
+    if (!m_backoff->is_null(backoff_term)) {
+      assert (level < m_backoff->order());
+      p0 = m_backoff->terms_at_level(level)*prob(backoff_term, topic, level+1);
+    }
+    else
+      p0 = m_term_p0;
+  }
+  //for (int i=0; i<level+1; ++i) std::cerr << "  ";
+  //std::cerr << "PYPTopics::word_pyps_p0(" << term << "," << topic << "," << level << ") = " << p0 << std::endl;
+  return p0;
+}
+
+PYPTopics::F PYPTopics::prob(const Term& term, int topic, int level) const {
+  //for (int i=0; i<level+1; ++i) std::cerr << "  ";
+  //std::cerr << "PYPTopics::prob(" << term << "," << topic << "," << level << " " << factor << ")" << std::endl;
+
+  F p0 = word_pyps_p0(term, topic, level);
+  F p_w_k = m_word_pyps.at(level).at(topic).prob(term, p0);
+
+  //for (int i=0; i<level+1; ++i) std::cerr << "  ";
+  //std::cerr << "PYPTopics::prob(" << term << "," << topic << "," << level << ") = " << p_w_k << std::endl;
+
+  return p_w_k;
+}
+
+int PYPTopics::max(const DocumentId& doc, const Term& term) {
+  //std::cerr << "PYPTopics::max(" << doc << "," << term << ")" << std::endl;
+  // collect probs
+  F current_max=0.0;
+  int current_topic=-1;
+  for (int k=0; k<m_num_topics; ++k) {
+    F p_w_k = prob(term, k);
+    F p_k_d = m_document_pyps[doc].prob(k, m_topic_p0);
+    F prob = (p_w_k*p_k_d);
+    if (prob > current_max) {
+      current_max = prob;
+      current_topic = k;
+    }
+  }
+  assert(current_topic >= 0);
+  return current_topic;
+}
+
+std::ostream& PYPTopics::print_document_topics(std::ostream& out) const {
+  for (CorpusTopics::const_iterator corpusIt=m_corpus_topics.begin(); 
+       corpusIt != m_corpus_topics.end(); ++corpusIt) {
+    int term_index=0;
+    for (DocumentTopics::const_iterator docIt=corpusIt->begin();
+         docIt != corpusIt->end(); ++docIt, ++term_index) {
+      if (term_index) out << " ";
+      out << *docIt;
+    }
+    out << std::endl;
+  }
+  return out;
+}
+
+std::ostream& PYPTopics::print_topic_terms(std::ostream& out) const {
+  for (PYPs::const_iterator pypsIt=m_word_pyps.front().begin(); 
+       pypsIt != m_word_pyps.front().end(); ++pypsIt) {
+    int term_index=0;
+    for (PYP<int>::const_iterator termIt=pypsIt->begin();
+         termIt != pypsIt->end(); ++termIt, ++term_index) {
+      if (term_index) out << " ";
+      out << termIt->first << ":" << termIt->second;
+    }
+    out << std::endl;
+  }
+  return out;
+}
diff --git a/gi/pyp-topics/src/pyp-topics.hh b/gi/pyp-topics/src/pyp-topics.hh
new file mode 100644
index 00000000..4036985a
--- /dev/null
+++ b/gi/pyp-topics/src/pyp-topics.hh
@@ -0,0 +1,52 @@
+#ifndef PYP_TOPICS_HH
+#define PYP_TOPICS_HH
+
+#include <vector>
+#include <iostream>
+
+#include "pyp.hh"
+#include "corpus.hh"
+
+
+class PYPTopics {
+public:
+  typedef std::vector<int> DocumentTopics;
+  typedef std::vector<DocumentTopics> CorpusTopics;
+  typedef long double F;
+
+public:
+  PYPTopics(int num_topics) : m_num_topics(num_topics), m_word_pyps(1), m_backoff(0) {}
+
+  void sample(const Corpus& corpus, int samples);
+  int sample(const DocumentId& doc, const Term& term);
+  int max(const DocumentId& doc, const Term& term);
+
+  void set_backoff(const std::string& filename) {
+    m_backoff.reset(new TermBackoff);
+    m_backoff->read(filename);
+    m_word_pyps.clear();
+    m_word_pyps.resize(m_backoff->order(), PYPs());
+  }
+
+  F prob(const Term& term, int topic, int level=0) const;
+  void decrement(const Term& term, int topic, int level=0);
+  void increment(const Term& term, int topic, int level=0);
+
+  std::ostream& print_document_topics(std::ostream& out) const;
+  std::ostream& print_topic_terms(std::ostream& out) const;
+
+private:
+  F word_pyps_p0(const Term& term, int topic, int level) const;
+
+  int m_num_topics;
+  F m_term_p0, m_topic_p0, m_backoff_p0;
+
+  CorpusTopics m_corpus_topics;
+  typedef std::vector< PYP<int> > PYPs;
+  PYPs m_document_pyps;
+  std::vector<PYPs> m_word_pyps;
+
+  std::auto_ptr<TermBackoff> m_backoff;
+};
+
+#endif // PYP_TOPICS_HH
diff --git a/gi/pyp-topics/src/pyp.hh b/gi/pyp-topics/src/pyp.hh
new file mode 100644
index 00000000..b06c6021
--- /dev/null
+++ b/gi/pyp-topics/src/pyp.hh
@@ -0,0 +1,477 @@
+#ifndef _pyp_hh
+#define _pyp_hh
+
+#include <math.h>
+#include <map>
+#include <tr1/unordered_map>
+
+#include "log_add.h"
+#include "gammadist.h"
+#include "slice-sampler.h"
+#include "mt19937ar.h"
+
+//
+// Pitman-Yor process with customer and table tracking
+//
+
+template <typename Dish, typename Hash=std::tr1::hash<Dish> >
+class PYP : protected std::tr1::unordered_map<Dish, int, Hash>
+{
+public:
+  using std::tr1::unordered_map<Dish,int>::const_iterator;
+  using std::tr1::unordered_map<Dish,int>::iterator;
+  using std::tr1::unordered_map<Dish,int>::begin;
+  using std::tr1::unordered_map<Dish,int>::end;
+
+  PYP(long double a, long double b, Hash hash=Hash());
+
+  int increment(Dish d, long double p0);
+  int decrement(Dish d);
+
+  // lookup functions
+  int count(Dish d) const;
+  long double prob(Dish dish, long double p0) const;
+  long double prob(Dish dish, long double dcd, long double dca, 
+                   long double dtd, long double dta, long double p0) const;
+
+  int num_customers() const { return _total_customers; }
+  int num_types() const { return std::tr1::unordered_map<Dish,int>::size(); }
+  bool empty() const { return _total_customers == 0; }
+
+  long double log_prob(Dish dish, long double log_p0) const;
+  // nb. d* are NOT logs
+  long double log_prob(Dish dish, long double dcd, long double dca, 
+                       long double dtd, long double dta, long double log_p0) const;
+
+  int num_tables(Dish dish) const;
+  int num_tables() const;
+
+  long double a() const { return _a; }
+  void set_a(long double a) { _a = a; }
+
+  long double b() const { return _b; }
+  void set_b(long double b) { _b = b; }
+
+  void clear();
+  std::ostream& debug_info(std::ostream& os) const;
+
+  long double log_restaurant_prob() const;
+  long double log_prior() const;
+  static long double log_prior_a(long double a, long double beta_a, long double beta_b);
+  static long double log_prior_b(long double b, long double gamma_c, long double gamma_s);
+
+  void resample_prior();
+  void resample_prior_a();
+  void resample_prior_b();
+
+private:
+  long double _a, _b; // parameters of the Pitman-Yor distribution
+  long double _a_beta_a, _a_beta_b; // parameters of Beta prior on a
+  long double _b_gamma_s, _b_gamma_c; // parameters of Gamma prior on b
+
+  struct TableCounter
+  {
+    TableCounter() : tables(0) {};
+    int tables;
+    std::map<int, int> table_histogram; // num customers at table -> number tables
+  };
+  typedef std::tr1::unordered_map<Dish, TableCounter, Hash> DishTableType;
+  DishTableType _dish_tables;
+  int _total_customers, _total_tables;
+
+  // Function objects for calculating the parts of the log_prob for 
+  // the parameters a and b
+  struct resample_a_type {
+    int n, m; long double b, a_beta_a, a_beta_b;
+    const DishTableType& dish_tables;
+    resample_a_type(int n, int m, long double b, long double a_beta_a, 
+                    long double a_beta_b, const DishTableType& dish_tables)
+      : n(n), m(m), b(b), a_beta_a(a_beta_a), a_beta_b(a_beta_b), dish_tables(dish_tables) {}
+
+    long double operator() (long double proposed_a) const {
+      long double log_prior = log_prior_a(proposed_a, a_beta_a, a_beta_b);
+      long double log_prob = 0.0;
+      long double lgamma1a = lgamma(1.0 - proposed_a);
+      for (typename DishTableType::const_iterator dish_it=dish_tables.begin(); dish_it != dish_tables.end(); ++dish_it) 
+        for (std::map<int, int>::const_iterator table_it=dish_it->second.table_histogram.begin(); 
+             table_it !=dish_it->second.table_histogram.end(); ++table_it) 
+          log_prob += (table_it->second * (lgamma(table_it->first - proposed_a) - lgamma1a));
+
+      log_prob += (proposed_a == 0.0 ? (m-1.0)*log(b) 
+                   : ((m-1.0)*log(proposed_a) + lgamma((m-1.0) + b/proposed_a) - lgamma(b/proposed_a)));
+      assert(std::isfinite(log_prob));
+      return log_prob + log_prior;
+    }
+  };
+
+  struct resample_b_type {
+    int n, m; long double a, b_gamma_c, b_gamma_s;
+    resample_b_type(int n, int m, long double a, long double b_gamma_c, long double b_gamma_s)
+      : n(n), m(m), a(a), b_gamma_c(b_gamma_c), b_gamma_s(b_gamma_s) {}
+
+    long double operator() (long double proposed_b) const {
+      long double log_prior = log_prior_b(proposed_b, b_gamma_c, b_gamma_s);
+      long double log_prob = 0.0;
+      log_prob += (a == 0.0  ? (m-1.0)*log(proposed_b) 
+                  : ((m-1.0)*log(a) + lgamma((m-1.0) + proposed_b/a) - lgamma(proposed_b/a)));
+      log_prob += (lgamma(1.0+proposed_b) - lgamma(n+proposed_b));
+      return log_prob + log_prior;
+    }
+  };
+};
+
+template <typename Dish, typename Hash>
+PYP<Dish,Hash>::PYP(long double a, long double b, Hash cmp)
+: std::tr1::unordered_map<Dish, int, Hash>(), _a(a), _b(b), 
+  _a_beta_a(1), _a_beta_b(1), _b_gamma_s(1), _b_gamma_c(1),
+  //_a_beta_a(1), _a_beta_b(1), _b_gamma_s(10), _b_gamma_c(0.1),
+  _total_customers(0), _total_tables(0)
+{
+//  std::cerr << "\t##PYP<Dish,Hash>::PYP(a=" << _a << ",b=" << _b << ")" << std::endl;
+}
+
+template <typename Dish, typename Hash>
+long double 
+PYP<Dish,Hash>::prob(Dish dish, long double p0) const
+{
+  int c = count(dish), t = num_tables(dish);
+  long double r = num_tables() * _a + _b;
+  //std::cerr << "\t\t\t\tPYP<Dish,Hash>::prob(" << dish << "," << p0 << ") c=" << c << " r=" << r << std::endl;
+  if (c > 0)
+    return (c - _a * t + r * p0) / (num_customers() + _b);
+  else
+    return r * p0 / (num_customers() + _b);
+}
+
+template <typename Dish, typename Hash>
+long double 
+PYP<Dish,Hash>::prob(Dish dish, long double dcd, long double dca, 
+                     long double dtd, long double dta, long double p0)
+const
+{
+  int c = count(dish) + dcd, t = num_tables(dish) + dtd;
+  long double r = (num_tables() + dta) * _a + _b;
+  if (c > 0)
+    return (c - _a * t + r * p0) / (num_customers() + dca + _b);
+  else
+    return r * p0 / (num_customers() + dca + _b);
+}
+
+template <typename Dish, typename Hash>
+long double 
+PYP<Dish,Hash>::log_prob(Dish dish, long double log_p0) const
+{
+  using std::log;
+  int c = count(dish), t = num_tables(dish);
+  long double r = log(num_tables() * _a + b);
+  if (c > 0)
+    return Log<double>::add(log(c - _a * t), r + log_p0)
+      - log(num_customers() + _b);
+  else
+    return r + log_p0 - log(num_customers() + b);
+}
+
+template <typename Dish, typename Hash>
+long double 
+PYP<Dish,Hash>::log_prob(Dish dish, long double dcd, long double dca, 
+                         long double dtd, long double dta, long double log_p0)
+const
+{
+  using std::log;
+  int c = count(dish) + dcd, t = num_tables(dish) + dtd;
+  long double r = log((num_tables() + dta) * _a + b);
+  if (c > 0)
+    return Log<double>::add(log(c - _a * t), r + log_p0)
+      - log(num_customers() + dca + _b);
+  else
+    return r + log_p0 - log(num_customers() + dca + b);
+}
+
+template <typename Dish, typename Hash>
+int 
+PYP<Dish,Hash>::increment(Dish dish, long double p0) {
+  int delta = 0;
+  TableCounter &tc = _dish_tables[dish];
+
+  // seated on a new or existing table?
+  int c = count(dish), t = num_tables(dish), T = num_tables();
+  long double pshare = (c > 0) ? (c - _a*t) : 0.0;
+  long double pnew = (_b + _a*T) * p0;
+  assert (pshare >= 0.0);
+  //assert (pnew > 0.0);
+
+  if (mt_genrand_res53() < pnew / (pshare + pnew)) {
+    // assign to a new table
+    tc.tables += 1;
+    tc.table_histogram[1] += 1;
+    _total_tables += 1;
+    delta = 1;
+  }
+  else {
+    // randomly assign to an existing table
+    // remove constant denominator from inner loop
+    long double r = mt_genrand_res53() * (c - _a*t);
+    for (std::map<int,int>::iterator
+         hit = tc.table_histogram.begin();
+         hit != tc.table_histogram.end(); ++hit) {
+      r -= ((hit->first - _a) * hit->second);
+      if (r <= 0) {
+        tc.table_histogram[hit->first+1] += 1;
+        hit->second -= 1;
+        if (hit->second == 0)
+          tc.table_histogram.erase(hit);
+        break;
+      }
+    }
+    if (r > 0) {
+      std::cerr << r << " " << c << " " << _a << " " << t << std::endl;
+      assert(false);
+    }
+    delta = 0;
+  }
+
+  std::tr1::unordered_map<Dish,int,Hash>::operator[](dish) += 1;
+  _total_customers += 1;
+
+  return delta;
+}
+
+template <typename Dish, typename Hash>
+int 
+PYP<Dish,Hash>::count(Dish dish) const
+{
+  typename std::tr1::unordered_map<Dish, int>::const_iterator 
+    dcit = find(dish);
+  if (dcit != end())
+    return dcit->second;
+  else
+    return 0;
+}
+
+template <typename Dish, typename Hash>
+int 
+PYP<Dish,Hash>::decrement(Dish dish)
+{
+  typename std::tr1::unordered_map<Dish, int>::iterator dcit = find(dish);
+  if (dcit == end()) {
+    std::cerr << dish << std::endl;
+    assert(false);
+  } 
+
+  int delta = 0;
+
+  typename std::tr1::unordered_map<Dish, TableCounter>::iterator dtit = _dish_tables.find(dish);
+  if (dtit == _dish_tables.end()) {
+    std::cerr << dish << std::endl;
+    assert(false);
+  } 
+  TableCounter &tc = dtit->second;
+
+  //std::cerr << "\tdecrement for " << dish << "\n";
+  //std::cerr << "\tBEFORE histogram: " << tc.table_histogram << " ";
+  //std::cerr << "count: " << count(dish) << " ";
+  //std::cerr << "tables: " << tc.tables << "\n";
+
+  long double r = mt_genrand_res53() * count(dish);
+  for (std::map<int,int>::iterator hit = tc.table_histogram.begin();
+       hit != tc.table_histogram.end(); ++hit)
+  {
+    //r -= (hit->first - _a) * hit->second;
+    r -= (hit->first) * hit->second;
+    if (r <= 0)
+    {
+      if (hit->first > 1)
+        tc.table_histogram[hit->first-1] += 1;
+      else
+      {
+        delta = -1;
+        tc.tables -= 1;
+        _total_tables -= 1;
+      }
+
+      hit->second -= 1;
+      if (hit->second == 0) tc.table_histogram.erase(hit);
+      break;
+    }
+  }
+  if (r > 0) {
+    std::cerr << r << " " << count(dish) << " " << _a << " " << num_tables(dish) << std::endl;
+    assert(false);
+  }
+
+  // remove the customer
+  dcit->second -= 1;
+  _total_customers -= 1;
+  assert(dcit->second >= 0);
+  if (dcit->second == 0) {
+    erase(dcit);
+    _dish_tables.erase(dtit);
+    //std::cerr << "\tAFTER histogram: Empty\n";
+  }
+  else {
+    //std::cerr << "\tAFTER histogram: " << _dish_tables[dish].table_histogram << " ";
+    //std::cerr << "count: " << count(dish) << " ";
+    //std::cerr << "tables: " << _dish_tables[dish].tables << "\n";
+  }
+
+  return delta;
+}
+
+template <typename Dish, typename Hash>
+int 
+PYP<Dish,Hash>::num_tables(Dish dish) const
+{
+  typename std::tr1::unordered_map<Dish, TableCounter, Hash>::const_iterator 
+    dtit = _dish_tables.find(dish);
+
+  //assert(dtit != _dish_tables.end());
+  if (dtit == _dish_tables.end())
+    return 0;
+
+  return dtit->second.tables;
+}
+
+template <typename Dish, typename Hash>
+int 
+PYP<Dish,Hash>::num_tables() const
+{
+  return _total_tables;
+}
+
+template <typename Dish, typename Hash>
+std::ostream&
+PYP<Dish,Hash>::debug_info(std::ostream& os) const
+{
+  int hists = 0, tables = 0;
+  for (typename std::tr1::unordered_map<Dish, TableCounter, Hash>::const_iterator 
+       dtit = _dish_tables.begin(); dtit != _dish_tables.end(); ++dtit)
+  {
+    hists += dtit->second.table_histogram.size();
+    tables += dtit->second.tables;
+
+    assert(dtit->second.tables > 0);
+    assert(!dtit->second.table_histogram.empty());
+
+    for (std::map<int,int>::const_iterator 
+         hit = dtit->second.table_histogram.begin();
+         hit != dtit->second.table_histogram.end(); ++hit)
+      assert(hit->second > 0);
+  }
+
+  os << "restaurant has " 
+    << _total_customers << " customers; "
+    << _total_tables << " tables; " 
+    << tables << " tables'; " 
+    << num_types() << " dishes; "
+    << _dish_tables.size() << " dishes'; and "
+    << hists << " histogram entries\n";
+
+  return os;
+}
+
+template <typename Dish, typename Hash>
+void 
+PYP<Dish,Hash>::clear()
+{
+  this->std::tr1::unordered_map<Dish,int,Hash>::clear();
+  _dish_tables.clear();
+  _total_tables = _total_customers = 0;
+}
+
+// log_restaurant_prob returns the log probability of the PYP table configuration.
+// Excludes Hierarchical P0 term which must be calculated separately.
+template <typename Dish, typename Hash>
+long double 
+PYP<Dish,Hash>::log_restaurant_prob() const {
+  if (_total_customers < 1)
+    return (long double)0.0;
+
+  long double log_prob = 0.0;
+  long double lgamma1a = lgamma(1.0-_a);
+
+  //std::cerr << "-------------------\n" << std::endl;
+  for (typename DishTableType::const_iterator dish_it=_dish_tables.begin(); 
+       dish_it != _dish_tables.end(); ++dish_it) {
+    for (std::map<int, int>::const_iterator table_it=dish_it->second.table_histogram.begin(); 
+         table_it !=dish_it->second.table_histogram.end(); ++table_it) {
+      log_prob += (table_it->second * (lgamma(table_it->first - _a) - lgamma1a));
+      //std::cerr << "|" << dish_it->first->parent << " --> " << dish_it->first->rhs << " " << table_it->first << " " << table_it->second << " " << log_prob;
+    }
+  }
+  //std::cerr << std::endl;
+
+  log_prob += (_a == (long double)0.0 ? (_total_tables-1.0)*log(_b) : (_total_tables-1.0)*log(_a) + lgamma((_total_tables-1.0) + _b/_a) - lgamma(_b/_a));
+  //std::cerr << "\t\t" << log_prob << std::endl;
+  log_prob += (lgamma(1.0 + _b) - lgamma(_total_customers + _b));
+
+  //std::cerr << _total_customers << " " << _total_tables << " " << log_prob << " " << log_prior() << std::endl;
+  //std::cerr << _a << " " << _b << std::endl;
+  if (!std::isfinite(log_prob)) {
+    assert(false);
+  }
+  //return log_prob;
+  return log_prob + log_prior();
+}
+
+template <typename Dish, typename Hash>
+long double 
+PYP<Dish,Hash>::log_prior() const {
+  long double prior = 0.0;
+  if (_a_beta_a > 0.0 && _a_beta_b > 0.0 && _a > 0.0)
+    prior += log_prior_a(_a, _a_beta_a, _a_beta_b);
+  if (_b_gamma_s > 0.0 && _b_gamma_c > 0.0)
+    prior += log_prior_b(_b, _b_gamma_c, _b_gamma_s);
+
+  return prior;
+}
+
+template <typename Dish, typename Hash>
+long double 
+PYP<Dish,Hash>::log_prior_a(long double a, long double beta_a, long double beta_b) {
+  return lbetadist(a, beta_a, beta_b); 
+}
+
+template <typename Dish, typename Hash>
+long double 
+PYP<Dish,Hash>::log_prior_b(long double b, long double gamma_c, long double gamma_s) {
+  return lgammadist(b, gamma_c, gamma_s); 
+}
+
+template <typename Dish, typename Hash>
+void 
+PYP<Dish,Hash>::resample_prior() {
+  for (int num_its=5; num_its >= 0; --num_its) {
+    resample_prior_b();
+    resample_prior_a();
+  }
+  resample_prior_b();
+}
+
+template <typename Dish, typename Hash>
+void 
+PYP<Dish,Hash>::resample_prior_b() {
+  if (_total_tables == 0) 
+    return;
+
+  int niterations = 10;   // number of resampling iterations
+  //std::cerr << "\n## resample_prior_b(), initial a = " << _a << ", b = " << _b << std::endl;
+  resample_b_type b_log_prob(_total_customers, _total_tables, _a, _b_gamma_c, _b_gamma_s);
+  _b = slice_sampler1d(b_log_prob, _b, mt_genrand_res53, (long double) 0.0, std::numeric_limits<long double>::infinity(), 
+                       (long double) 0.0, niterations, 100*niterations);
+  //std::cerr << "\n## resample_prior_b(), final a = " << _a << ", b = " << _b << std::endl;
+}
+
+template <typename Dish, typename Hash>
+void 
+PYP<Dish,Hash>::resample_prior_a() {
+  if (_total_tables == 0) 
+    return;
+
+  int niterations = 10;
+  //std::cerr << "\n## Initial a = " << _a << ", b = " << _b << std::endl;
+  resample_a_type a_log_prob(_total_customers, _total_tables, _b, _a_beta_a, _a_beta_b, _dish_tables);
+  _a = slice_sampler1d(a_log_prob, _a, mt_genrand_res53, std::numeric_limits<long double>::min(), 
+                       (long double) 1.0, (long double) 0.0, niterations, 100*niterations);
+}
+
+#endif
diff --git a/gi/pyp-topics/src/slice-sampler.h b/gi/pyp-topics/src/slice-sampler.h
new file mode 100644
index 00000000..3108a0f7
--- /dev/null
+++ b/gi/pyp-topics/src/slice-sampler.h
@@ -0,0 +1,192 @@
+//! slice-sampler.h is an MCMC slice sampler
+//!
+//! Mark Johnson, 1st August 2008
+
+#ifndef SLICE_SAMPLER_H
+#define SLICE_SAMPLER_H
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <iostream>
+#include <limits>
+
+//! slice_sampler_rfc_type{} returns the value of a user-specified
+//! function if the argument is within range, or - infinity otherwise
+//
+template <typename F, typename Fn, typename U>
+struct slice_sampler_rfc_type {
+  F min_x, max_x;
+  const Fn& f;
+  U max_nfeval, nfeval;
+  slice_sampler_rfc_type(F min_x, F max_x, const Fn& f, U max_nfeval) 
+    : min_x(min_x), max_x(max_x), f(f), max_nfeval(max_nfeval), nfeval(0) { }
+    
+  F operator() (F x) {
+    if (min_x < x && x < max_x) {
+      assert(++nfeval <= max_nfeval);
+      F fx = f(x);
+      assert(std::isfinite(fx));
+      return fx;
+    }
+    else
+      return -std::numeric_limits<F>::infinity();
+  }
+};  // slice_sampler_rfc_type{}
+
+//! slice_sampler1d() implements the univariate "range doubling" slice sampler
+//! described in Neal (2003) "Slice Sampling", The Annals of Statistics 31(3), 705-767.
+//
+template <typename F, typename LogF, typename Uniform01>
+F slice_sampler1d(const LogF& logF0,               //!< log of function to sample
+		  F x,                             //!< starting point
+		  Uniform01& u01,                  //!< uniform [0,1) random number generator
+		  F min_x = -std::numeric_limits<F>::infinity(),  //!< minimum value of support
+		  F max_x = std::numeric_limits<F>::infinity(),   //!< maximum value of support
+		  F w = 0.0,                       //!< guess at initial width
+		  unsigned nsamples=1,             //!< number of samples to draw
+		  unsigned max_nfeval=200)         //!< max number of function evaluations
+{
+  typedef unsigned U;
+  slice_sampler_rfc_type<F,LogF,U> logF(min_x, max_x, logF0, max_nfeval);
+
+  assert(std::isfinite(x));
+
+  if (w <= 0.0) {                           // set w to a default width 
+    if (min_x > -std::numeric_limits<F>::infinity() && max_x < std::numeric_limits<F>::infinity())
+      w = (max_x - min_x)/4;
+    else
+      w = std::max(((x < 0.0) ? -x : x)/4, (F) 0.1);
+  }
+  assert(std::isfinite(w));
+
+  F logFx = logF(x);
+  for (U sample = 0; sample < nsamples; ++sample) {
+    F logY = logFx + log(u01()+1e-100);     //! slice logFx at this value
+    assert(std::isfinite(logY));
+
+    F xl = x - w*u01();                     //! lower bound on slice interval
+    F logFxl = logF(xl);
+    F xr = xl + w;                          //! upper bound on slice interval
+    F logFxr = logF(xr);
+
+    while (logY < logFxl || logY < logFxr)  // doubling procedure
+      if (u01() < 0.5) 
+	logFxl = logF(xl -= xr - xl);
+      else
+	logFxr = logF(xr += xr - xl);
+	
+    F xl1 = xl;
+    F xr1 = xr;
+    while (true) {                          // shrinking procedure
+      F x1 = xl1 + u01()*(xr1 - xl1);
+      if (logY < logF(x1)) {
+	F xl2 = xl;                         // acceptance procedure
+	F xr2 = xr; 
+	bool d = false;
+	while (xr2 - xl2 > 1.1*w) {
+	  F xm = (xl2 + xr2)/2;
+	  if ((x < xm && x1 >= xm) || (x >= xm && x1 < xm))
+	    d = true;
+	  if (x1 < xm)
+	    xr2 = xm;
+	  else
+	    xl2 = xm;
+	  if (d && logY >= logF(xl2) && logY >= logF(xr2))
+	    goto unacceptable;
+	}
+	x = x1;
+	goto acceptable;
+      }
+      goto acceptable;
+    unacceptable:
+      if (x1 < x)                           // rest of shrinking procedure
+	xl1 = x1;
+      else 
+	xr1 = x1;
+    }
+  acceptable:
+    w = (4*w + (xr1 - xl1))/5;              // update width estimate
+  }
+  return x;
+}
+
+/*
+//! slice_sampler1d() implements a 1-d MCMC slice sampler.
+//! It should be correct for unimodal distributions, but
+//! not for multimodal ones.
+//
+template <typename F, typename LogP, typename Uniform01>
+F slice_sampler1d(const LogP& logP,     //!< log of distribution to sample
+		  F x,                  //!< initial sample
+		  Uniform01& u01,       //!< uniform random number generator
+		  F min_x = -std::numeric_limits<F>::infinity(),  //!< minimum value of support
+		  F max_x = std::numeric_limits<F>::infinity(),   //!< maximum value of support
+		  F w = 0.0,            //!< guess at initial width
+		  unsigned nsamples=1,  //!< number of samples to draw
+		  unsigned max_nfeval=200)  //!< max number of function evaluations
+{
+  typedef unsigned U;
+  assert(std::isfinite(x));
+  if (w <= 0.0) {
+    if (min_x > -std::numeric_limits<F>::infinity() && max_x < std::numeric_limits<F>::infinity())
+      w = (max_x - min_x)/4;
+    else
+      w = std::max(((x < 0.0) ? -x : x)/4, 0.1);
+  }
+  // TRACE4(x, min_x, max_x, w);
+  F logPx = logP(x);
+  assert(std::isfinite(logPx));
+  U nfeval = 1;
+  for (U sample = 0; sample < nsamples; ++sample) {
+    F x0 = x;
+    F logU = logPx + log(u01()+1e-100);
+    assert(std::isfinite(logU));
+    F r = u01();
+    F xl = std::max(min_x, x - r*w);
+    F xr = std::min(max_x, x + (1-r)*w);
+    // TRACE3(x, logPx, logU);
+    while (xl > min_x && logP(xl) > logU) {
+      xl -= w;
+      w *= 2;
+      ++nfeval;
+      if (nfeval >= max_nfeval)
+	std::cerr << "## Error: nfeval = " << nfeval << ", max_nfeval = " << max_nfeval << ", sample = " << sample << ", nsamples = " << nsamples << ", r = " << r << ", w = " << w << ", xl = " << xl << std::endl;
+      assert(nfeval < max_nfeval);
+    }
+    xl = std::max(xl, min_x);
+    while (xr < max_x && logP(xr) > logU) {
+      xr += w;
+      w *= 2;
+      ++nfeval;
+      if (nfeval >= max_nfeval)
+	std::cerr << "## Error: nfeval = " << nfeval << ", max_nfeval = " << max_nfeval << ", sample = " << sample << ", nsamples = " << nsamples << ", r = " << r << ", w = " << w << ", xr = " << xr << std::endl;
+      assert(nfeval < max_nfeval);
+    }
+    xr = std::min(xr, max_x);
+    while (true) {
+      r = u01();
+      x = r*xl + (1-r)*xr;
+      assert(std::isfinite(x));
+      logPx = logP(x);
+      // TRACE4(logPx, x, xl, xr);
+      assert(std::isfinite(logPx));
+      ++nfeval;
+      if (nfeval >= max_nfeval)
+	std::cerr << "## Error: nfeval = " << nfeval << ", max_nfeval = " << max_nfeval << ", sample = " << sample << ", nsamples = " << nsamples << ", r = " << r << ", w = " << w << ", xl = " << xl << ", xr = " << xr << ", x = " << x << std::endl;
+      assert(nfeval < max_nfeval);
+      if (logPx > logU)
+        break;
+      else if (x > x0)
+          xr = x;
+        else
+          xl = x;
+    }
+    // w = (4*w + (xr-xl))/5;   // gradually adjust w
+  }
+  // TRACE2(logPx, x);
+  return x;
+}  // slice_sampler1d()
+*/
+
+#endif  // SLICE_SAMPLER_H
diff --git a/gi/pyp-topics/src/train.cc b/gi/pyp-topics/src/train.cc
new file mode 100644
index 00000000..0d107f11
--- /dev/null
+++ b/gi/pyp-topics/src/train.cc
@@ -0,0 +1,136 @@
+// STL
+#include <iostream>
+#include <fstream>
+
+// Boost
+#include <boost/program_options/parsers.hpp>
+#include <boost/program_options/variables_map.hpp>
+#include <boost/scoped_ptr.hpp>
+
+// Local
+#include "pyp-topics.hh"
+#include "corpus.hh"
+#include "gzstream.hh"
+#include "mt19937ar.h"
+
+static const char *REVISION = "$Revision: 0.1 $";
+
+// Namespaces
+using namespace boost;
+using namespace boost::program_options;
+using namespace std;
+
+int main(int argc, char **argv)
+{
+  std::cout << "Pitman Yor topic models: Copyright 2010 Phil Blunsom\n";
+  std::cout << REVISION << '\n' << std::endl;
+
+  ////////////////////////////////////////////////////////////////////////////////////////////
+  // Command line processing
+  variables_map vm; 
+
+  // Command line processing
+  options_description cmdline_specific("Command line specific options");
+  cmdline_specific.add_options()
+    ("help,h", "print help message")
+    ("config,c", value<string>(), "config file specifying additional command line options")
+    ;
+  options_description generic("Allowed options");
+  generic.add_options()
+    ("documents,d", value<string>(), "file containing the documents")
+    ("topics,t", value<int>()->default_value(50), "number of topics")
+    ("document-topics-out,o", value<string>(), "file to write the document topics to")
+    ("topic-words-out,w", value<string>(), "file to write the topic word distribution to")
+    ("samples,s", value<int>()->default_value(10), "number of sampling passes through the data")
+    ("test-corpus", value<string>(), "file containing the test data")
+    ("backoff-paths", value<string>(), "file containing the term backoff paths")
+    ;
+  options_description config_options, cmdline_options;
+  config_options.add(generic);
+  cmdline_options.add(generic).add(cmdline_specific);
+
+  store(parse_command_line(argc, argv, cmdline_options), vm); 
+  if (vm.count("config") > 0) {
+    ifstream config(vm["config"].as<string>().c_str());
+    store(parse_config_file(config, cmdline_options), vm); 
+  }
+  notify(vm);
+  ////////////////////////////////////////////////////////////////////////////////////////////
+
+  if (vm.count("help")) { 
+    cout << cmdline_options << "\n"; 
+    return 1; 
+  }
+
+  if (vm.count("documents") == 0) {
+    cerr << "Please specify a file containing the documents." << endl;
+    cout << cmdline_options << "\n"; 
+    return 1;
+  }
+
+  // seed the random number generator
+  //mt_init_genrand(time(0));
+
+  // read the data
+  Corpus corpus;
+  corpus.read(vm["documents"].as<string>());
+
+  // run the sampler
+  PYPTopics model(vm["topics"].as<int>());
+
+  // read the backoff dictionary
+  if (vm.count("backoff-paths"))
+    model.set_backoff(vm["backoff-paths"].as<string>());
+
+  // train the sampler
+  model.sample(corpus, vm["samples"].as<int>());
+
+  if (vm.count("document-topics-out")) {
+    ogzstream documents_out(vm["document-topics-out"].as<string>().c_str());
+    //model.print_document_topics(documents_out);
+
+    int document_id=0;
+    for (Corpus::const_iterator corpusIt=corpus.begin(); 
+         corpusIt != corpus.end(); ++corpusIt, ++document_id) {
+      std::vector<int> unique_terms;
+      for (Document::const_iterator docIt=corpusIt->begin();
+           docIt != corpusIt->end(); ++docIt) {
+        if (unique_terms.empty() || *docIt != unique_terms.back())
+          unique_terms.push_back(*docIt);
+      }
+      documents_out << unique_terms.size();
+      for (std::vector<int>::const_iterator termIt=unique_terms.begin();
+           termIt != unique_terms.end(); ++termIt)
+        documents_out << " " << *termIt << ":" << model.max(document_id, *termIt);
+      documents_out << std::endl;
+    }
+    documents_out.close();
+  }
+
+  if (vm.count("topic-words-out")) {
+    ogzstream topics_out(vm["topic-words-out"].as<string>().c_str());
+    model.print_topic_terms(topics_out);
+    topics_out.close();
+  }
+
+  if (vm.count("test-corpus")) {
+    TestCorpus test_corpus;
+    test_corpus.read(vm["test-corpus"].as<string>());
+    ogzstream topics_out((vm["test-corpus"].as<string>() + ".topics.gz").c_str());
+
+    for (TestCorpus::const_iterator corpusIt=test_corpus.begin(); 
+         corpusIt != test_corpus.end(); ++corpusIt) {
+      int index=0;
+      for (DocumentTerms::const_iterator instanceIt=corpusIt->begin();
+           instanceIt != corpusIt->end(); ++instanceIt, ++index) {
+        int topic = model.max(instanceIt->doc, instanceIt->term);
+        if (index != 0) topics_out << " ";
+        topics_out << topic;
+      }
+      topics_out << std::endl;
+    }
+    topics_out.close();
+  }
+
+  return 0;
+}
diff --git a/gi/pyp-topics/src/utility.h b/gi/pyp-topics/src/utility.h
new file mode 100644
index 00000000..405a5b0a
--- /dev/null
+++ b/gi/pyp-topics/src/utility.h
@@ -0,0 +1,962 @@
+// utility.h
+//
+// (c) Mark Johnson, 24th January 2005
+//
+// modified 6th May 2002 to ensure write/read consistency, fixed 18th July 2002
+// modified 14th July 2002 to include insert() (generic inserter)
+// modified 26th September 2003 to use mapped_type instead of data_type
+// 25th August 2004 added istream >> const char*
+// 24th January 2005 added insert_newkey()
+//
+// Defines:
+//  loop macros foreach, cforeach
+//  dfind (default find function)
+//  afind (find function that asserts key exists)
+//  insert_newkey (inserts a new key into a map)
+//  insert (generic inserter into standard data structures)
+//  disjoint (set operation)
+//  first_lessthan and second_lessthan (compares elements of pairs)
+//
+// Simplified interfaces to STL routines:
+//
+//  includes (simplified interface)
+//  set_intersection (simplified interface)
+//  inserter (simplified interface)
+//  max_element (simplified interface)
+//  min_element (simplified interface)
+//  hash functions for pairs, vectors, lists, slists and maps
+//  input and output for pairs and vectors
+//  resource_usage (interface improved)
+
+
+#ifndef UTILITY_H
+#define UTILITY_H
+
+#include <algorithm>
+// #include <boost/smart_ptr.hpp>    // Comment out this line if boost is not used
+#include <cassert>
+#include <cmath>
+#include <cctype>
+#include <cstdio>
+#include <unordered_map>
+#include <unordered_set>
+#include <ext/slist>
+#include <iostream>
+#include <iterator>
+#include <list>
+#include <map>
+#include <set>
+#include <string>
+#include <utility>
+#include <vector>
+#include <memory>
+
+#if (__GNUC__ > 3) || (__GNUC__ >= 3 && __GNUC_MINOR__ >= 1)
+#define EXT_NAMESPACE __gnu_cxx
+#else
+#define EXT_NAMESPACE std
+#endif
+
+namespace ext = EXT_NAMESPACE;
+
+inline float power(float x, float y) { return powf(x, y); }
+inline double power(double x, double y) { return pow(x, y); }
+inline long double power(long double x, long double y) { return powl(x, y); }
+
+typedef unsigned U;
+typedef long double F;  // slower than double, but underflows less
+
+///////////////////////////////////////////////////////////////////////////
+//                                                                       //
+//                              Looping constructs                       //
+//                                                                       //
+///////////////////////////////////////////////////////////////////////////
+
+// foreach is a simple loop construct
+//
+//   STORE should be an STL container
+//   TYPE is the typename of STORE
+//   VAR will be defined as a local variable of type TYPE::iterator
+//
+#define foreach(TYPE, VAR, STORE) \
+  for (TYPE::iterator VAR = (STORE).begin(); VAR != (STORE).end(); ++VAR)
+
+// cforeach is just like foreach, except that VAR is a const_iterator
+//
+//   STORE should be an STL container
+//   TYPE is the typename of STORE
+//   VAR will be defined as a local variable of type TYPE::const_iterator
+//
+#define cforeach(TYPE, VAR, STORE) \
+  for (TYPE::const_iterator VAR = (STORE).begin(); VAR != (STORE).end(); ++VAR)
+
+
+///////////////////////////////////////////////////////////////////////////
+//                                                                       //
+//                             Map searching                             //
+//                                                                       //
+// dfind(map, key) returns the key's value in map, or map's default      //
+//   value if no such key exists (the default value is not inserted)     //
+//                                                                       //
+// afind(map, key) returns a reference to the key's value in map, and    //
+//    asserts that this value exists                                     //
+//                                                                       //
+///////////////////////////////////////////////////////////////////////////
+
+// dfind(Map, Key) returns the value Map associates with Key, or the
+//  Map's default value if no such Key exists
+//
+template <class Map, class Key>
+inline typename Map::mapped_type dfind(Map& m, const Key& k)
+{
+  typename Map::iterator i = m.find(k);
+  if (i == m.end())
+    return typename Map::mapped_type();
+  else
+    return i->second;
+}
+
+template <class Map, class Key>
+inline const typename Map::mapped_type dfind(const Map& m, const Key& k)
+{
+  typename Map::const_iterator i = m.find(k);
+  if (i == m.end())
+    return typename Map::mapped_type();
+  else
+    return i->second;
+}
+
+
+// afind(map, key) returns a reference to the value associated
+//  with key in map.  It uses assert to check that the key's value
+//  is defined.
+//
+template <class Map, class Key>
+inline typename Map::mapped_type& afind(Map& m, const Key& k)
+{
+  typename Map::iterator i = m.find(k);
+  assert(i != m.end());
+  return i->second;
+}
+
+template <class Map, class Key>
+inline const typename Map::mapped_type& afind(const Map& m, const Key& k)
+{
+  typename Map::const_iterator i = m.find(k);
+  assert(i != m.end());
+  return i->second;
+}
+
+//! insert_newkey(map, key, value) checks that map does not contain
+//! key, and binds key to value.
+//
+template <class Map, class Key, class Value>
+inline typename Map::value_type& 
+insert_newkey(Map& m, const Key& k,const Value& v) 
+{
+  std::pair<typename Map::iterator, bool> itb 
+    = m.insert(Map::value_type(k, v));
+  assert(itb.second);
+  return *(itb.first);
+}  // insert_newkey()
+
+
+///////////////////////////////////////////////////////////////////////////
+//                                                                       //
+//                        Insert operations                              //
+//                                                                       //
+///////////////////////////////////////////////////////////////////////////
+
+
+template <typename T>
+void insert(std::list<T>& xs, const T& x) {
+  xs.push_back(x);
+}
+
+template <typename T>
+void insert(std::set<T>& xs, const T& x) {
+  xs.insert(x);
+}
+
+template <typename T>
+void insert(std::vector<T>& xs, const T& x) {
+  xs.push_back(x);
+}
+
+
+///////////////////////////////////////////////////////////////////////////
+//                                                                       //
+//                Additional versions of standard algorithms             //
+//                                                                       //
+///////////////////////////////////////////////////////////////////////////
+
+template <typename Set1, typename Set2>
+inline bool includes(const Set1& set1, const Set2& set2)
+{
+  return std::includes(set1.begin(), set1.end(), set2.begin(), set2.end());
+}
+
+template <typename Set1, typename Set2, typename Compare>
+inline bool includes(const Set1& set1, const Set2& set2, Compare comp)
+{
+  return std::includes(set1.begin(), set1.end(), set2.begin(), set2.end(), comp);
+}
+
+
+template <typename InputIter1, typename InputIter2>
+bool disjoint(InputIter1 first1, InputIter1 last1,
+              InputIter2 first2, InputIter2 last2)
+{
+  while (first1 != last1 && first2 != last2)
+    if (*first1 < *first2)
+      ++first1;
+    else if (*first2 < *first1)
+      ++first2;
+    else // *first1 == *first2
+      return false;
+  return true;
+}
+
+template <typename InputIter1, typename InputIter2, typename Compare>
+bool disjoint(InputIter1 first1, InputIter1 last1,
+              InputIter2 first2, InputIter2 last2, Compare comp)
+{
+  while (first1 != last1 && first2 != last2)
+    if (comp(*first1, *first2))
+      ++first1;
+    else if (comp(*first2, *first1))
+      ++first2;
+    else // *first1 == *first2
+      return false;
+  return true;
+}
+
+template <typename Set1, typename Set2>
+inline bool disjoint(const Set1& set1, const Set2& set2)
+{
+  return disjoint(set1.begin(), set1.end(), set2.begin(), set2.end());
+}
+
+template <typename Set1, typename Set2, typename Compare>
+inline bool disjoint(const Set1& set1, const Set2& set2, Compare comp)
+{
+  return disjoint(set1.begin(), set1.end(), set2.begin(), set2.end(), comp);
+}
+
+
+template <typename Set1, typename Set2, typename OutputIterator>
+inline OutputIterator set_intersection(const Set1& set1, const Set2& set2, 
+                                       OutputIterator result)
+{
+  return set_intersection(set1.begin(), set1.end(), set2.begin(), set2.end(), result);
+}
+
+template <typename Set1, typename Set2, typename OutputIterator, typename Compare>
+inline OutputIterator set_intersection(const Set1& set1, const Set2& set2, 
+                                       OutputIterator result, Compare comp)
+{
+  return set_intersection(set1.begin(), set1.end(), set2.begin(), set2.end(), result, comp);
+}
+
+
+template <typename Container>
+inline std::insert_iterator<Container> inserter(Container& container)
+{
+  return std::inserter(container, container.begin());
+}
+
+// max_element
+//
+template <class Es> inline typename Es::iterator max_element(Es& es)
+{
+  return std::max_element(es.begin(), es.end());
+}
+
+template <class Es> inline typename Es::const_iterator max_element(const Es& es)
+{
+  return std::max_element(es.begin(), es.end());
+}
+
+template <class Es, class BinaryPredicate> 
+inline typename Es::iterator max_element(Es& es, BinaryPredicate comp)
+{
+  return std::max_element(es.begin(), es.end(), comp);
+}
+
+template <class Es, class BinaryPredicate> 
+inline typename Es::const_iterator max_element(const Es& es, BinaryPredicate comp)
+{
+  return std::max_element(es.begin(), es.end(), comp);
+}
+
+// min_element
+//
+template <class Es> inline typename Es::iterator min_element(Es& es)
+{
+  return std::min_element(es.begin(), es.end());
+}
+
+template <class Es> inline typename Es::const_iterator min_element(const Es& es)
+{
+  return std::min_element(es.begin(), es.end());
+}
+
+template <class Es, class BinaryPredicate> 
+inline typename Es::iterator min_element(Es& es, BinaryPredicate comp)
+{
+  return std::min_element(es.begin(), es.end(), comp);
+}
+
+template <class Es, class BinaryPredicate> 
+inline typename Es::const_iterator min_element(const Es& es, BinaryPredicate comp)
+{
+  return std::min_element(es.begin(), es.end(), comp);
+}
+
+// first_lessthan and second_lessthan
+//
+struct first_lessthan {
+  template <typename T1, typename T2>
+    bool operator() (const T1& e1, const T2& e2) {
+      return e1.first < e2.first;
+    }
+};
+
+struct second_lessthan {
+  template <typename T1, typename T2>
+    bool operator() (const T1& e1, const T2& e2) {
+      return e1.second < e2.second;
+    }
+};
+
+// first_greaterthan and second_greaterthan
+//
+struct first_greaterthan {
+  template <typename T1, typename T2>
+    bool operator() (const T1& e1, const T2& e2) {
+      return e1.first > e2.first;
+    }
+};
+
+struct second_greaterthan {
+  template <typename T1, typename T2>
+    bool operator() (const T1& e1, const T2& e2) {
+      return e1.second > e2.second;
+    }
+};
+
+
+///////////////////////////////////////////////////////////////////////////
+//                                                                       //
+//                          hash<> specializations                       //
+//                                                                       //
+// These must be in namespace std.   They permit the corresponding STL   //
+// container to be used as a key in an STL hash table.                   //
+//                                                                       //
+///////////////////////////////////////////////////////////////////////////
+
+//namespace EXT_NAMESPACE {
+namespace std {
+  /*
+  // hash function for bool
+  //
+  template <> struct hash<bool>
+  {
+  size_t operator() (bool b) const 
+  {
+  return b;
+  } // operator()
+  }; // hash<bool>{}
+
+  // hash function for double
+  //
+  template <> struct hash<double>
+  {
+  size_t operator() (double d) const 
+  {
+  int exponent;
+  double fraction = frexp(d, &exponent);
+  return size_t(exponent) ^ size_t(1000000.0*(fabs(fraction-0.5)));
+  } // operator()
+  }; // hash<double>{}
+
+  // hash function for strings
+  //
+  template <> struct hash<std::string> 
+  {
+  size_t operator()(const std::string& s) const 
+  {
+  typedef std::string::const_iterator CI;
+
+  unsigned long h = 0; 
+  unsigned long g;
+  CI p = s.begin();
+  CI end = s.end();
+
+  while (p!=end) {
+  h = (h << 4) + (*p++);
+  if ((g = h&0xf0000000)) {
+  h = h ^ (g >> 24);
+  h = h ^ g;
+  }}
+  return size_t(h);
+  }  // operator()
+  };  // hash<string>{}
+
+*/
+  // hash function for arbitrary pairs
+  //
+  template<class T1, class T2> struct hash<std::pair<T1,T2> > {
+    size_t operator()(const std::pair<T1,T2>& p) const
+    {
+      size_t h1 = hash<T1>()(p.first);
+      size_t h2 = hash<T2>()(p.second);
+      return h1 ^ (h1 >> 1) ^ h2 ^ (h2 << 1);
+    }
+  };
+
+
+  // hash function for vectors
+  //
+  template<class T> struct hash<std::vector<T> > 
+  { //  This is the fn hashpjw of Aho, Sethi and Ullman, p 436.
+    size_t operator()(const std::vector<T>& s) const 
+    {
+      typedef typename std::vector<T>::const_iterator CI;
+
+      unsigned long h = 0; 
+      unsigned long g;
+      CI p = s.begin();
+      CI end = s.end();
+
+      while (p!=end) {
+        h = (h << 5) + hash<T>()(*p++);
+        if ((g = h&0xff000000)) {
+          h = h ^ (g >> 23);
+          h = h ^ g;
+        }}
+        return size_t(h);
+    }
+  };
+
+  // hash function for slists
+  //
+  template<class T> struct hash<ext::slist<T> > 
+  { //  This is the fn hashpjw of Aho, Sethi and Ullman, p 436.
+    size_t operator()(const ext::slist<T>& s) const 
+    {
+      typedef typename ext::slist<T>::const_iterator CI;
+
+      unsigned long h = 0; 
+      unsigned long g;
+      CI p = s.begin();
+      CI end = s.end();
+
+      while (p!=end) {
+        h = (h << 7) + hash<T>()(*p++);
+        if ((g = h&0xff000000)) {
+          h = h ^ (g >> 23);
+          h = h ^ g;
+        }}
+        return size_t(h);
+    }
+  };
+
+  // hash function for maps
+  //
+  template<typename T1, typename T2> struct hash<std::map<T1,T2> >
+  {
+    size_t operator()(const std::map<T1,T2>& m) const
+    {
+      typedef typename std::map<T1,T2> M;
+      typedef typename M::const_iterator CI;
+
+      unsigned long h = 0;
+      unsigned long g;
+      CI p = m.begin();
+      CI end = m.end();
+
+      while (p != end) {
+        h = (h << 11) + hash<typename M::value_type>()(*p++);
+        if ((g = h&0xff000000)) {
+          h = h ^ (g >> 23);
+          h = h ^ g;
+        }}
+        return size_t(h);
+    }
+  };
+
+} // namespace EXT_NAMESPACE
+
+
+
+///////////////////////////////////////////////////////////////////////////
+//                                                                       //
+//                           Write/Read code                             //
+//                                                                       //
+// These routines should possess write/read invariance IF their elements //
+// also have write-read invariance.  Whitespace, '(' and ')' are used as //
+// delimiters.                                                           //
+//                                                                       //
+///////////////////////////////////////////////////////////////////////////
+
+
+// Define istream >> const char* so that it consumes the characters from the
+// istream.  Just as in scanf, a space consumes an arbitrary amount of whitespace.
+//
+inline std::istream& operator>> (std::istream& is, const char* cp)
+{
+  if (*cp == '\0')
+    return is;
+  else if (*cp == ' ') {
+    char c;
+    if (is.get(c)) {
+      if (isspace(c))
+        return is >> cp;
+      else {
+        is.unget();
+        return is >> (cp+1);
+      }
+    }
+    else {
+      is.clear(is.rdstate() & ~std::ios::failbit);  // clear failbit
+      return is >> (cp+1);
+    }
+  }
+  else {
+    char c;
+    if (is.get(c)) {
+      if (c == *cp)
+        return is >> (cp+1);
+      else {
+        is.unget();
+        is.setstate(std::ios::failbit);
+      }
+    }
+    return is;
+  }
+}
+
+
+// Write out an auto_ptr object just as you would write out the pointer object
+//
+template <typename T> 
+inline std::ostream& operator<<(std::ostream& os, const std::auto_ptr<T>& sp)
+{
+  return os << sp.get();
+}
+
+
+// Pairs
+//
+template <class T1, class T2> 
+std::ostream& operator<< (std::ostream& os, const std::pair<T1,T2>& p)
+{
+  return os << '(' << p.first << ' ' << p.second << ')';
+}
+
+template <class T1, class T2>
+std::istream& operator>> (std::istream& is, std::pair<T1,T2>& p)
+{
+  char c;
+  if (is >> c) {
+    if (c == '(') {
+      if (is >> p.first >> p.second >> c && c == ')')
+        return is;
+      else 
+        is.setstate(std::ios::badbit);
+    }
+    else
+      is.putback(c);
+  }
+  is.setstate(std::ios::failbit);
+  return is;
+}
+
+// Lists
+//
+template <class T>
+std::ostream& operator<< (std::ostream& os, const std::list<T>& xs)
+{
+  os << '(';
+  for (typename std::list<T>::const_iterator xi = xs.begin(); xi != xs.end(); ++xi) {
+    if (xi != xs.begin())
+      os << ' ';
+    os << *xi;
+  }
+  return os << ')';
+}
+
+template <class T>
+std::istream& operator>> (std::istream& is, std::list<T>& xs)
+{
+  char c;                          // This code avoids unnecessary copy
+  if (is >> c) {                   // read the initial '('
+    if (c == '(') {
+      xs.clear();                  // clear the list
+      do {
+        xs.push_back(T());         // create a new elt in list
+        is >> xs.back();           // read element
+      }
+      while (is.good());           // read as long as possible
+      xs.pop_back();               // last read failed; pop last elt
+      is.clear(is.rdstate() & ~std::ios::failbit);  // clear failbit
+      if (is >> c && c == ')')     // read terminating ')'
+        return is;                 // successful return
+      else 
+        is.setstate(std::ios::badbit); // something went wrong, set badbit
+    }
+    else                           // c is not '('
+      is.putback(c);               //  put c back into input
+  }
+  is.setstate(std::ios::failbit);  // read failed, set failbit
+  return is;
+}
+
+// Vectors
+//
+template <class T>
+std::ostream& operator<< (std::ostream& os, const std::vector<T>& xs)
+{
+  os << '(';
+  for (typename std::vector<T>::const_iterator xi = xs.begin(); xi != xs.end(); ++xi) {
+    if (xi != xs.begin())
+      os << ' ';
+    os << *xi;
+  }
+  return os << ')';
+}
+
+template <class T>
+std::istream& operator>> (std::istream& is, std::vector<T>& xs)
+{
+  char c;                          // This code avoids unnecessary copy
+  if (is >> c) {                   // read the initial '('
+    if (c == '(') {
+      xs.clear();                  // clear the list
+      do {
+        xs.push_back(T());         // create a new elt in list
+        is >> xs.back();           // read element
+      }
+      while (is.good());           // read as long as possible
+      xs.pop_back();               // last read failed; pop last elt
+      is.clear(is.rdstate() & ~std::ios::failbit);  // clear failbit
+      if (is >> c && c == ')')     // read terminating ')'
+        return is;                 // successful return
+      else 
+        is.setstate(std::ios::badbit); // something went wrong, set badbit
+    }
+    else                           // c is not '('
+      is.putback(c);               //  put c back into input
+  }
+  is.setstate(std::ios::failbit);  // read failed, set failbit
+  return is;
+}
+
+// Slists
+//
+template <class T>
+std::ostream& operator<< (std::ostream& os, const ext::slist<T>& xs)
+{
+  os << '(';
+  for (typename ext::slist<T>::const_iterator xi = xs.begin(); xi != xs.end(); ++xi) {
+    if (xi != xs.begin())
+      os << ' ';
+    os << *xi;
+  }
+  return os << ')';
+}
+
+template <class T>
+std::istream& operator>> (std::istream& is, ext::slist<T>& xs)
+{
+  char c; 
+  if (is >> c) { 
+    if (c == '(') {
+      xs.clear();
+      T e;
+      if (is >> e) {   
+        xs.push_front(e);
+        typename ext::slist<T>::iterator xi = xs.begin();
+        while (is >> e)
+          xi = xs.insert_after(xi, e);
+        is.clear(is.rdstate() & ~std::ios::failbit);
+        if (is >> c && c == ')')
+          return is;
+        else 
+          is.setstate(std::ios::badbit); 
+      }
+      else { // empty list
+        is.clear(is.rdstate() & ~std::ios::failbit);
+        if (is >> c && c == ')')
+          return is;
+        else  // didn't see closing ')'
+          is.setstate(std::ios::badbit);
+      }
+    }
+    else  // didn't read '('
+      is.putback(c);
+  }
+  is.setstate(std::ios::failbit);
+  return is;
+}
+
+// Sets
+//
+template <class T>
+std::ostream& operator<< (std::ostream& os, const std::set<T>& s)
+{
+  os << '(';
+  for (typename std::set<T>::const_iterator i = s.begin(); i != s.end(); ++i) {
+    if (i != s.begin())
+      os << ' ';
+    os << *i;
+  }
+  return os << ')';
+}
+
+template <class T>
+std::istream& operator>> (std::istream& is, std::set<T>& s)
+{
+  char c;
+  if (is >> c) {
+    if (c == '(') {
+      s.clear();
+      T e;
+      while (is >> e)
+        s.insert(e);
+      is.clear(is.rdstate() & ~std::ios::failbit);
+      if (is >> c && c == ')')
+        return is;
+      else
+        is.setstate(std::ios::badbit);
+    }
+    else
+      is.putback(c);
+  }
+  is.setstate(std::ios::failbit);
+  return is;
+}
+
+// Hash_sets
+//
+template <class T>
+std::ostream& operator<< (std::ostream& os, const std::unordered_set<T>& s)
+{
+  os << '(';
+  for (typename std::unordered_set<T>::const_iterator i = s.begin(); i != s.end(); ++i) {
+    if (i != s.begin())
+      os << ' ';
+    os << *i;
+  }
+  return os << ')';
+}
+
+template <class T>
+std::istream& operator>> (std::istream& is, std::unordered_set<T>& s)
+{
+  char c;
+  if (is >> c) {
+    if (c == '(') {
+      s.clear();
+      T e;
+      while (is >> e)
+        s.insert(e);
+      is.clear(is.rdstate() & ~std::ios::failbit);
+      if (is >> c && c == ')')
+        return is;
+      else
+        is.setstate(std::ios::badbit);
+    }
+    else
+      is.putback(c);
+  }
+  is.setstate(std::ios::failbit);
+  return is;
+}
+
+
+// Maps
+//
+template <class Key, class Value>
+std::ostream& operator<< (std::ostream& os, const std::map<Key,Value>& m)
+{
+  typedef std::map<Key,Value> M;
+  os << '(';
+  for (typename M::const_iterator it = m.begin(); it != m.end(); ++it) {
+    if (it != m.begin())
+      os << ' ';
+    os << *it;
+  }
+  return os << ")";
+}
+
+template <class Key, class Value>
+std::istream& operator>> (std::istream& is, std::map<Key,Value>& m)
+{
+  char c;
+  if (is >> c) {
+    if (c == '(') {
+      m.clear();
+      std::pair<Key,Value> e;
+      while (is >> e)
+        m.insert(e);
+      is.clear(is.rdstate() & ~std::ios::failbit);
+      if (is >> c && c == ')')
+        return is;
+      else
+        is.setstate(std::ios::badbit);
+    }
+    else
+      is.putback(c);
+  }
+  is.setstate(std::ios::failbit);
+  return is;
+}
+
+// Hash_maps
+//
+template <class Key, class Value>
+std::ostream& operator<< (std::ostream& os, const std::unordered_map<Key,Value>& m)
+{
+  typedef std::unordered_map<Key,Value> M;
+  os << '(';
+  for (typename M::const_iterator it = m.begin(); it != m.end(); ++it) {
+    if (it != m.begin())
+      os << ' ';
+    os << *it;
+  }
+  return os << ")";
+}
+
+template <class Key, class Value>
+std::istream& operator>> (std::istream& is, std::unordered_map<Key,Value>& m)
+{
+  char c;
+  if (is >> c) {
+    if (c == '(') {
+      m.clear();
+      std::pair<Key,Value> e;
+      while (is >> e)
+        m.insert(e);
+      is.clear(is.rdstate() & ~std::ios::failbit);
+      if (is >> c && c == ')')
+        return is;
+      else
+        is.setstate(std::ios::badbit);
+    }
+    else
+      is.putback(c);
+  }
+  is.setstate(std::ios::failbit);
+  return is;
+}
+
+
+///////////////////////////////////////////////////////////////////////////
+//                                                                       //
+//                       Boost library additions                         //
+//                                                                       //
+///////////////////////////////////////////////////////////////////////////
+
+#ifdef BOOST_SHARED_PTR_HPP_INCLUDED
+
+// enhancements to boost::shared_ptr so it can be used with hash
+//
+namespace std {
+  template <typename T> struct equal_to<boost::shared_ptr<T> > 
+    : public binary_function<boost::shared_ptr<T>, boost::shared_ptr<T>, bool> {
+      bool operator() (const boost::shared_ptr<T>& p1, const boost::shared_ptr<T>& p2) const {
+        return equal_to<T*>()(p1.get(), p2.get());
+      }
+    };
+}  // namespace std
+
+//namespace EXT_NAMESPACE {
+namespace std {
+  template <typename T> struct hash<boost::shared_ptr<T> > {
+    size_t operator() (const boost::shared_ptr<T>& a) const {
+      return hash<T*>()(a.get());
+    }
+  };
+}  // namespace ext
+
+template <typename T> 
+inline std::ostream& operator<< (std::ostream& os, const boost::shared_ptr<T>& sp)
+{
+  return os << sp.get();
+}
+
+#endif  // BOOST_SHARED_PTR_HPP_INCLUDED
+
+struct resource_usage { };
+
+#ifndef __i386
+inline std::ostream& operator<< (std::ostream& os, resource_usage r)
+{
+  return os;
+}
+#else // Assume we are on a 586 linux
+inline std::ostream& operator<< (std::ostream& os, resource_usage r)
+{
+  FILE* fp = fopen("/proc/self/stat", "r");
+  assert(fp);
+  int utime;
+  int stime;
+  unsigned int vsize;
+  unsigned int rss;
+  int result = 
+    fscanf(fp, "%*d %*s %*c %*d %*d %*d %*d %*d %*u %*u %*u %*u %*u %d %d %*d %*d %*d %*d"
+           "%*u %*u %*d %u %u", &utime, &stime, &vsize, &rss);
+  assert(result == 4);
+  fclose(fp);
+  // s << "utime = " << utime << ", stime = " << stime << ", vsize = " << vsize << ", rss = " << rss
+  ;
+  // return s << "utime = " << utime << ", vsize = " << vsize;
+  return os << "utime " << float(utime)/1.0e2 << "s, vsize " 
+    << float(vsize)/1048576.0 << " Mb.";
+}
+#endif
+
+//! A default_value_type{} object is used to read an object from a stream,
+//! assigning a default value if the read fails.  Users should not need to
+//! construct such objects, but should use default_value() instead.
+//
+template <typename object_type, typename default_type>
+struct default_value_type {
+  object_type& object;
+  const default_type defaultvalue;
+  default_value_type(object_type& object, const default_type defaultvalue)
+    : object(object), defaultvalue(defaultvalue) { }
+};
+
+//! default_value() is used to read an object from a stream, assigning a
+//! default value if the read fails.  It returns a default_value_type{}
+//! object, which does the actual reading.
+//
+template <typename object_type, typename default_type>
+default_value_type<object_type,default_type>
+default_value(object_type& object, const default_type defaultvalue=default_type()) {
+  return default_value_type<object_type,default_type>(object, defaultvalue);
+}
+
+//! This version of operator>>() reads default_value_type{} from an input stream.
+//
+template <typename object_type, typename default_type>
+std::istream& operator>> (std::istream& is, 
+                          default_value_type<object_type, default_type> dv) {
+  if (is) {
+    if (is >> dv.object)
+      ;
+    else {
+      is.clear(is.rdstate() & ~std::ios::failbit);  // clear failbit
+      dv.object = dv.defaultvalue;
+    }
+  }
+  return is;
+}
+
+// inline F random1() { return rand()/(RAND_MAX+1.0); }
+inline F random1() { return mt_genrand_res53(); }
+
+#endif  // UTILITY_H