check in some experimental particle filtering code, some gitignore fixes

26 files changed, 6141 insertions, 3 deletions

diff --git a/.gitignore b/.gitignore
index 2a287bbc..5efe37b0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,5 +1,27 @@
+pro-train/.deps
+pro-train/mr_pro_map
+pro-train/mr_pro_reduce
+utils/reconstruct_weights
+decoder/.libs
+training/augment_grammar
+training/mpi_batch_optimize
+training/mpi_compute_cllh
+training/mpi_em_optimize
+training/mpi_extract_features
+training/mpi_extract_reachable
 klm/lm/build_binary
 extools/extractor_monolingual
+gi/pf/.deps
+gi/pf/brat
+gi/pf/cbgi
+gi/pf/dpnaive
+gi/pf/itg
+gi/pf/libpf.a
+gi/pf/pfbrat
+gi/pf/pfdist
+gi/pf/pfnaive
+gi/markov_al/.deps
+gi/markov_al/ml
 gi/posterior-regularisation/prjava/lib/*.jar
 klm/lm/libklm.a
 klm/util/.deps
@@ -120,4 +142,4 @@ gi/posterior-regularisation/prjava/lib/prjava-20100715.jar
 *.dvi
 *.ps
 *.toc
-*~
\ No newline at end of file
+*~
diff --git a/Makefile.am b/Makefile.am
index 98b4bac7..59c2fc0a 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -1,7 +1,7 @@
 # warning - the subdirectories in the following list should
 # be kept in topologically sorted order. Also, DO NOT introduce
 # cyclic dependencies between these directories!
-SUBDIRS = utils mteval klm/util klm/lm decoder phrasinator training mira vest pro-train extools
+SUBDIRS = utils mteval klm/util klm/lm decoder phrasinator training mira vest pro-train extools gi/pf gi/markov_al
 
 #gi/pyp-topics/src gi/clda/src gi/posterior-regularisation/prjava
 
diff --git a/configure.ac b/configure.ac
index 2e9cc36d..131a1705 100644
--- a/configure.ac
+++ b/configure.ac
@@ -113,4 +113,4 @@ then
   AM_CONDITIONAL([GLC], true)
 fi
 
-AC_OUTPUT(Makefile utils/Makefile mteval/Makefile extools/Makefile decoder/Makefile phrasinator/Makefile training/Makefile vest/Makefile pro-train/Makefile klm/util/Makefile klm/lm/Makefile mira/Makefile gi/pyp-topics/src/Makefile gi/clda/src/Makefile gi/cbgi/Makefile gi/ml/Makefile)
+AC_OUTPUT(Makefile utils/Makefile mteval/Makefile extools/Makefile decoder/Makefile phrasinator/Makefile training/Makefile vest/Makefile pro-train/Makefile klm/util/Makefile klm/lm/Makefile mira/Makefile gi/pyp-topics/src/Makefile gi/clda/src/Makefile gi/pf/Makefile gi/markov_al/Makefile)
diff --git a/gi/markov_al/Makefile.am b/gi/markov_al/Makefile.am
new file mode 100644
index 00000000..fe3e3349
--- /dev/null
+++ b/gi/markov_al/Makefile.am
@@ -0,0 +1,6 @@
+bin_PROGRAMS = ml
+
+ml_SOURCES = ml.cc
+
+AM_CPPFLAGS = -W -Wall -Wno-sign-compare -funroll-loops -I$(top_srcdir)/utils $(GTEST_CPPFLAGS) -I$(top_srcdir)/decoder
+AM_LDFLAGS = $(top_srcdir)/decoder/libcdec.a $(top_srcdir)/utils/libutils.a -lz
diff --git a/gi/markov_al/README b/gi/markov_al/README
new file mode 100644
index 00000000..9c10f7cd
--- /dev/null
+++ b/gi/markov_al/README
@@ -0,0 +1,2 @@
+Experimental translation models with Markovian dependencies.
+
diff --git a/gi/markov_al/ml.cc b/gi/markov_al/ml.cc
new file mode 100644
index 00000000..1e71edd6
--- /dev/null
+++ b/gi/markov_al/ml.cc
@@ -0,0 +1,470 @@
+#include <iostream>
+#include <tr1/unordered_map>
+
+#include <boost/shared_ptr.hpp>
+#include <boost/functional.hpp>
+#include <boost/program_options.hpp>
+#include <boost/program_options/variables_map.hpp>
+
+#include "tdict.h"
+#include "filelib.h"
+#include "sampler.h"
+#include "ccrp_onetable.h"
+#include "array2d.h"
+
+using namespace std;
+using namespace std::tr1;
+namespace po = boost::program_options;
+
+void PrintTopCustomers(const CCRP_OneTable<WordID>& crp) {
+  for (CCRP_OneTable<WordID>::const_iterator it = crp.begin(); it != crp.end(); ++it) {
+    cerr << "  " << TD::Convert(it->first) << " = " << it->second << endl;
+  }
+}
+
+void PrintAlignment(const vector<WordID>& src, const vector<WordID>& trg, const vector<unsigned char>& a) {
+  cerr << TD::GetString(src) << endl << TD::GetString(trg) << endl;
+  Array2D<bool> al(src.size(), trg.size());
+  for (int i = 0; i < a.size(); ++i)
+    if (a[i] != 255) al(a[i], i) = true;
+  cerr << al << endl;
+}
+
+void InitCommandLine(int argc, char** argv, po::variables_map* conf) {
+  po::options_description opts("Configuration options");
+  opts.add_options()
+        ("samples,s",po::value<unsigned>()->default_value(1000),"Number of samples")
+        ("input,i",po::value<string>(),"Read parallel data from")
+        ("random_seed,S",po::value<uint32_t>(), "Random seed");
+  po::options_description clo("Command line options");
+  clo.add_options()
+        ("config", po::value<string>(), "Configuration file")
+        ("help,h", "Print this help message and exit");
+  po::options_description dconfig_options, dcmdline_options;
+  dconfig_options.add(opts);
+  dcmdline_options.add(opts).add(clo);
+  
+  po::store(parse_command_line(argc, argv, dcmdline_options), *conf);
+  if (conf->count("config")) {
+    ifstream config((*conf)["config"].as<string>().c_str());
+    po::store(po::parse_config_file(config, dconfig_options), *conf);
+  }
+  po::notify(*conf);
+
+  if (conf->count("help") || (conf->count("input") == 0)) {
+    cerr << dcmdline_options << endl;
+    exit(1);
+  }
+}
+
+struct Unigram;
+struct Bigram {
+  Bigram() : trg(), cond() {}
+  Bigram(WordID prev, WordID cur, WordID t) : trg(t) { cond.first = prev; cond.second = cur; }
+  const pair<WordID,WordID>& ConditioningPair() const {
+    return cond;
+  }
+  WordID& prev_src() { return cond.first; }
+  WordID& cur_src() { return cond.second; }
+  const WordID& prev_src() const { return cond.first; }
+  const WordID& cur_src() const { return cond.second; }
+  WordID trg;
+ private:
+  pair<WordID, WordID> cond;
+};
+
+struct Unigram {
+  Unigram() : cur_src(), trg() {}
+  Unigram(WordID s, WordID t) : cur_src(s), trg(t) {}
+  WordID cur_src;
+  WordID trg;
+};
+
+ostream& operator<<(ostream& os, const Bigram& b) {
+  os << "( " << TD::Convert(b.trg) << " | " << TD::Convert(b.prev_src()) << " , " << TD::Convert(b.cur_src()) << " )";
+  return os;
+}
+
+ostream& operator<<(ostream& os, const Unigram& u) {
+  os << "( " << TD::Convert(u.trg) << " | " << TD::Convert(u.cur_src) << " )";
+  return os;
+}
+
+bool operator==(const Bigram& a, const Bigram& b) {
+  return a.trg == b.trg && a.cur_src() == b.cur_src() && a.prev_src() == b.prev_src();
+}
+
+bool operator==(const Unigram& a, const Unigram& b) {
+  return a.trg == b.trg && a.cur_src == b.cur_src;
+}
+
+size_t hash_value(const Bigram& b) {
+  size_t h = boost::hash_value(b.prev_src());
+  boost::hash_combine(h, boost::hash_value(b.cur_src()));
+  boost::hash_combine(h, boost::hash_value(b.trg));
+  return h;
+}
+
+size_t hash_value(const Unigram& u) {
+  size_t h = boost::hash_value(u.cur_src);
+  boost::hash_combine(h, boost::hash_value(u.trg));
+  return h;
+}
+
+void ReadParallelCorpus(const string& filename,
+                vector<vector<WordID> >* f,
+                vector<vector<WordID> >* e,
+                set<WordID>* vocab_f,
+                set<WordID>* vocab_e) {
+  f->clear();
+  e->clear();
+  vocab_f->clear();
+  vocab_e->clear();
+  istream* in;
+  if (filename == "-")
+    in = &cin;
+  else
+    in = new ifstream(filename.c_str());
+  assert(*in);
+  string line;
+  const WordID kDIV = TD::Convert("|||");
+  vector<WordID> tmp;
+  while(*in) {
+    getline(*in, line);
+    if (line.empty() && !*in) break;
+    e->push_back(vector<int>());
+    f->push_back(vector<int>());
+    vector<int>& le = e->back();
+    vector<int>& lf = f->back();
+    tmp.clear();
+    TD::ConvertSentence(line, &tmp);
+    bool isf = true;
+    for (unsigned i = 0; i < tmp.size(); ++i) {
+      const int cur = tmp[i];
+      if (isf) {
+        if (kDIV == cur) { isf = false; } else {
+          lf.push_back(cur);
+          vocab_f->insert(cur);
+        }
+      } else {
+        assert(cur != kDIV);
+        le.push_back(cur);
+        vocab_e->insert(cur);
+      }
+    }
+    assert(isf == false);
+  }
+  if (in != &cin) delete in;
+}
+
+struct UnigramModel {
+  UnigramModel(size_t src_voc_size, size_t trg_voc_size) :
+    unigrams(TD::NumWords() + 1, CCRP_OneTable<WordID>(1,1,1,1)),
+    p0(1.0 / trg_voc_size) {}
+
+  void increment(const Bigram& b) {
+    unigrams[b.cur_src()].increment(b.trg);
+  }
+
+  void decrement(const Bigram& b) {
+    unigrams[b.cur_src()].decrement(b.trg);
+  }
+
+  double prob(const Bigram& b) const {
+    const double q0 = unigrams[b.cur_src()].prob(b.trg, p0);
+    return q0;
+  }
+
+  double LogLikelihood() const {
+    double llh = 0;
+    for (unsigned i = 0; i < unigrams.size(); ++i) {
+      const CCRP_OneTable<WordID>& crp = unigrams[i];
+      if (crp.num_customers() > 0) {
+        llh += crp.log_crp_prob();
+        llh += crp.num_tables() * log(p0);
+      }
+    }
+    return llh;
+  }
+
+  void ResampleHyperparameters(MT19937* rng) {
+    for (unsigned i = 0; i < unigrams.size(); ++i)
+      unigrams[i].resample_hyperparameters(rng);
+  }
+
+  vector<CCRP_OneTable<WordID> > unigrams;  // unigrams[src].prob(trg, p0) = p(trg|src)
+
+  const double p0;
+};
+
+struct BigramModel {
+  BigramModel(size_t src_voc_size, size_t trg_voc_size) :
+    unigrams(TD::NumWords() + 1, CCRP_OneTable<WordID>(1,1,1,1)),
+    p0(1.0 / trg_voc_size) {}
+
+  void increment(const Bigram& b) {
+    BigramMap::iterator it = bigrams.find(b.ConditioningPair());
+    if (it == bigrams.end()) {
+      it = bigrams.insert(make_pair(b.ConditioningPair(), CCRP_OneTable<WordID>(1,1,1,1))).first;
+    }
+    if (it->second.increment(b.trg))
+      unigrams[b.cur_src()].increment(b.trg);
+  }
+
+  void decrement(const Bigram& b) {
+    BigramMap::iterator it = bigrams.find(b.ConditioningPair());
+    assert(it != bigrams.end());
+    if (it->second.decrement(b.trg)) {
+      unigrams[b.cur_src()].decrement(b.trg);
+      if (it->second.num_customers() == 0)
+        bigrams.erase(it);
+    }
+  }
+
+  double prob(const Bigram& b) const {
+    const double q0 = unigrams[b.cur_src()].prob(b.trg, p0);
+    const BigramMap::const_iterator it = bigrams.find(b.ConditioningPair());
+    if (it == bigrams.end()) return q0;
+    return it->second.prob(b.trg, q0);
+  }
+
+  double LogLikelihood() const {
+    double llh = 0;
+    for (unsigned i = 0; i < unigrams.size(); ++i) {
+      const CCRP_OneTable<WordID>& crp = unigrams[i];
+      if (crp.num_customers() > 0) {
+        llh += crp.log_crp_prob();
+        llh += crp.num_tables() * log(p0);
+      }
+    }
+    for (BigramMap::const_iterator it = bigrams.begin(); it != bigrams.end(); ++it) {
+      const CCRP_OneTable<WordID>& crp = it->second;
+      const WordID cur_src = it->first.second;
+      llh += crp.log_crp_prob();
+      for (CCRP_OneTable<WordID>::const_iterator bit = crp.begin(); bit != crp.end(); ++bit) {
+        llh += log(unigrams[cur_src].prob(bit->second, p0));
+      }
+    }
+    return llh;
+  }
+
+  void ResampleHyperparameters(MT19937* rng) {
+    for (unsigned i = 0; i < unigrams.size(); ++i)
+      unigrams[i].resample_hyperparameters(rng);
+    for (BigramMap::iterator it = bigrams.begin(); it != bigrams.end(); ++it)
+      it->second.resample_hyperparameters(rng);
+  }
+
+  typedef unordered_map<pair<WordID,WordID>, CCRP_OneTable<WordID>, boost::hash<pair<WordID,WordID> > > BigramMap;
+  BigramMap bigrams;  // bigrams[(src-1,src)].prob(trg, q0) = p(trg|src,src-1)
+  vector<CCRP_OneTable<WordID> > unigrams;  // unigrams[src].prob(trg, p0) = p(trg|src)
+
+  const double p0;
+};
+
+struct BigramAlignmentModel {
+  BigramAlignmentModel(size_t src_voc_size, size_t trg_voc_size) : bigrams(TD::NumWords() + 1, CCRP_OneTable<WordID>(1,1,1,1)), p0(1.0 / src_voc_size) {}
+  void increment(WordID prev, WordID next) {
+    bigrams[prev].increment(next);  // hierarchy?
+  }
+  void decrement(WordID prev, WordID next) {
+    bigrams[prev].decrement(next);  // hierarchy?
+  }
+  double prob(WordID prev, WordID next) {
+    return bigrams[prev].prob(next, p0);
+  }
+  double LogLikelihood() const {
+    double llh = 0;
+    for (unsigned i = 0; i < bigrams.size(); ++i) {
+      const CCRP_OneTable<WordID>& crp = bigrams[i];
+      if (crp.num_customers() > 0) {
+        llh += crp.log_crp_prob();
+        llh += crp.num_tables() * log(p0);
+      }
+    }
+    return llh;
+  }
+
+  vector<CCRP_OneTable<WordID> > bigrams;  // bigrams[prev].prob(next, p0) = p(next|prev)
+  const double p0;
+};
+
+struct Alignment {
+  vector<unsigned char> a;
+};
+
+int main(int argc, char** argv) {
+  po::variables_map conf;
+  InitCommandLine(argc, argv, &conf);
+  const unsigned samples = conf["samples"].as<unsigned>();
+
+  boost::shared_ptr<MT19937> prng;
+  if (conf.count("random_seed"))
+    prng.reset(new MT19937(conf["random_seed"].as<uint32_t>()));
+  else
+    prng.reset(new MT19937);
+  MT19937& rng = *prng;
+
+  vector<vector<WordID> > corpuse, corpusf;
+  set<WordID> vocabe, vocabf;
+  cerr << "Reading corpus...\n";
+  ReadParallelCorpus(conf["input"].as<string>(), &corpusf, &corpuse, &vocabf, &vocabe);
+  cerr << "F-corpus size: " << corpusf.size() << " sentences\t (" << vocabf.size() << " word types)\n";
+  cerr << "E-corpus size: " << corpuse.size() << " sentences\t (" << vocabe.size() << " word types)\n";
+  assert(corpusf.size() == corpuse.size());
+  const size_t corpus_len = corpusf.size();
+  const WordID kNULL = TD::Convert("<eps>");
+  const WordID kBOS = TD::Convert("<s>");
+  const WordID kEOS = TD::Convert("</s>");
+  Bigram TT(kBOS, TD::Convert("我"), TD::Convert("i"));
+  Bigram TT2(kBOS, TD::Convert("要"), TD::Convert("i"));
+
+  UnigramModel model(vocabf.size(), vocabe.size());
+  vector<Alignment> alignments(corpus_len);
+  for (unsigned ci = 0; ci < corpus_len; ++ci) {
+    const vector<WordID>& src = corpusf[ci];
+    const vector<WordID>& trg = corpuse[ci];
+    vector<unsigned char>& alg = alignments[ci].a;
+    alg.resize(trg.size());
+    int lenp1 = src.size() + 1;
+    WordID prev_src = kBOS;
+    for (int j = 0; j < trg.size(); ++j) {
+      int samp = lenp1 * rng.next();
+      --samp;
+      if (samp < 0) samp = 255;
+      alg[j] = samp;
+      WordID cur_src = (samp == 255 ? kNULL : src[alg[j]]);
+      Bigram b(prev_src, cur_src, trg[j]);
+      model.increment(b);
+      prev_src = cur_src;
+    }
+    Bigram b(prev_src, kEOS, kEOS);
+    model.increment(b);
+  }
+  cerr << "Initial LLH: " << model.LogLikelihood() << endl;
+
+  SampleSet<double> ss;
+  for (unsigned si = 0; si < 50; ++si) {
+    for (unsigned ci = 0; ci < corpus_len; ++ci) {
+      const vector<WordID>& src = corpusf[ci];
+      const vector<WordID>& trg = corpuse[ci];
+      vector<unsigned char>& alg = alignments[ci].a;
+      WordID prev_src = kBOS;
+      for (unsigned j = 0; j < trg.size(); ++j) {
+        unsigned char& a_j = alg[j];
+        WordID cur_e_a_j = (a_j == 255 ? kNULL : src[a_j]);
+        Bigram b(prev_src, cur_e_a_j, trg[j]);
+        //cerr << "DEC: " << b << "\t" << nextb << endl;
+        model.decrement(b);
+        ss.clear();
+        for (unsigned i = 0; i <= src.size(); ++i) {
+          const WordID cur_src = (i ? src[i-1] : kNULL);
+          b.cur_src() = cur_src;
+          ss.add(model.prob(b));
+        }
+        int sampled_a_j = rng.SelectSample(ss);
+        a_j = (sampled_a_j ? sampled_a_j - 1 : 255);
+        cur_e_a_j = (a_j == 255 ? kNULL : src[a_j]);
+        b.cur_src() = cur_e_a_j;
+        //cerr << "INC: " << b << "\t" << nextb << endl;
+        model.increment(b);
+        prev_src = cur_e_a_j;
+      }
+    }
+    cerr << '.' << flush;
+    if (si % 10 == 9) {
+      cerr << "[LLH prev=" << model.LogLikelihood();
+      //model.ResampleHyperparameters(&rng);
+      cerr << " new=" << model.LogLikelihood() << "]\n";
+      //pair<WordID,WordID> xx = make_pair(kBOS, TD::Convert("我"));
+      //PrintTopCustomers(model.bigrams.find(xx)->second);
+      cerr << "p(" << TT << ") = " << model.prob(TT) << endl;
+      cerr << "p(" << TT2 << ") = " << model.prob(TT2) << endl;
+      PrintAlignment(corpusf[0], corpuse[0], alignments[0].a);
+    }
+  }
+  {
+  // MODEL 2
+  BigramModel model(vocabf.size(), vocabe.size());
+  BigramAlignmentModel amodel(vocabf.size(), vocabe.size());
+  for (unsigned ci = 0; ci < corpus_len; ++ci) {
+    const vector<WordID>& src = corpusf[ci];
+    const vector<WordID>& trg = corpuse[ci];
+    vector<unsigned char>& alg = alignments[ci].a;
+    WordID prev_src = kBOS;
+    for (int j = 0; j < trg.size(); ++j) {
+      WordID cur_src = (alg[j] == 255 ? kNULL : src[alg[j]]);
+      Bigram b(prev_src, cur_src, trg[j]);
+      model.increment(b);
+      amodel.increment(prev_src, cur_src);
+      prev_src = cur_src;
+    }
+    amodel.increment(prev_src, kEOS);
+    Bigram b(prev_src, kEOS, kEOS);
+    model.increment(b);
+  }
+  cerr << "Initial LLH: " << model.LogLikelihood() << " " << amodel.LogLikelihood() << endl;
+
+  SampleSet<double> ss;
+  for (unsigned si = 0; si < samples; ++si) {
+    for (unsigned ci = 0; ci < corpus_len; ++ci) {
+      const vector<WordID>& src = corpusf[ci];
+      const vector<WordID>& trg = corpuse[ci];
+      vector<unsigned char>& alg = alignments[ci].a;
+      WordID prev_src = kBOS;
+      for (unsigned j = 0; j < trg.size(); ++j) {
+        unsigned char& a_j = alg[j];
+        WordID cur_e_a_j = (a_j == 255 ? kNULL : src[a_j]);
+        Bigram b(prev_src, cur_e_a_j, trg[j]);
+        WordID next_src = kEOS;
+        WordID next_trg = kEOS;
+        if (j < (trg.size() - 1)) {
+          next_src = (alg[j+1] == 255 ? kNULL : src[alg[j + 1]]);
+          next_trg = trg[j + 1];
+        }
+        Bigram nextb(cur_e_a_j, next_src, next_trg);
+        //cerr << "DEC: " << b << "\t" << nextb << endl;
+        model.decrement(b);
+        model.decrement(nextb);
+        amodel.decrement(prev_src, cur_e_a_j);
+        amodel.decrement(cur_e_a_j, next_src);
+        ss.clear();
+        for (unsigned i = 0; i <= src.size(); ++i) {
+          const WordID cur_src = (i ? src[i-1] : kNULL);
+          b.cur_src() = cur_src;
+          ss.add(model.prob(b) * model.prob(nextb) * amodel.prob(prev_src, cur_src) * amodel.prob(cur_src, next_src));
+          //cerr << log(ss[ss.size() - 1]) << "\t" << b << endl;
+        }
+        int sampled_a_j = rng.SelectSample(ss);
+        a_j = (sampled_a_j ? sampled_a_j - 1 : 255);
+        cur_e_a_j = (a_j == 255 ? kNULL : src[a_j]);
+        b.cur_src() = cur_e_a_j;
+        nextb.prev_src() = cur_e_a_j;
+        //cerr << "INC: " << b << "\t" << nextb << endl;
+        //exit(1);
+        model.increment(b);
+        model.increment(nextb);
+        amodel.increment(prev_src, cur_e_a_j);
+        amodel.increment(cur_e_a_j, next_src);
+        prev_src = cur_e_a_j;
+      }
+    }
+    cerr << '.' << flush;
+    if (si % 10 == 9) {
+      cerr << "[LLH prev=" << (model.LogLikelihood() + amodel.LogLikelihood());
+      //model.ResampleHyperparameters(&rng);
+      cerr << " new=" << model.LogLikelihood() << "]\n";
+      pair<WordID,WordID> xx = make_pair(kBOS, TD::Convert("我"));
+      cerr << "p(" << TT << ") = " << model.prob(TT) << endl;
+      cerr << "p(" << TT2 << ") = " << model.prob(TT2) << endl;
+      pair<WordID,WordID> xx2 = make_pair(kBOS, TD::Convert("要"));
+      PrintTopCustomers(model.bigrams.find(xx)->second);
+      //PrintTopCustomers(amodel.bigrams[TD::Convert("<s>")]);
+      //PrintTopCustomers(model.unigrams[TD::Convert("<eps>")]);
+      PrintAlignment(corpusf[0], corpuse[0], alignments[0].a);
+    }
+  }
+  }
+  return 0;
+}
+
diff --git a/gi/pf/Makefile.am b/gi/pf/Makefile.am
new file mode 100644
index 00000000..c9764ad5
--- /dev/null
+++ b/gi/pf/Makefile.am
@@ -0,0 +1,21 @@
+bin_PROGRAMS = cbgi brat dpnaive pfbrat pfdist itg pfnaive
+
+noinst_LIBRARIES = libpf.a
+libpf_a_SOURCES = base_measures.cc reachability.cc cfg_wfst_composer.cc
+
+itg_SOURCES = itg.cc
+
+dpnaive_SOURCES = dpnaive.cc
+
+pfdist_SOURCES = pfdist.cc
+
+pfnaive_SOURCES = pfnaive.cc
+
+cbgi_SOURCES = cbgi.cc
+
+brat_SOURCES = brat.cc
+
+pfbrat_SOURCES = pfbrat.cc
+
+AM_CPPFLAGS = -W -Wall -Wno-sign-compare -funroll-loops -I$(top_srcdir)/utils $(GTEST_CPPFLAGS) -I$(top_srcdir)/decoder
+AM_LDFLAGS = libpf.a $(top_srcdir)/decoder/libcdec.a $(top_srcdir)/utils/libutils.a -lz
diff --git a/gi/pf/README b/gi/pf/README
new file mode 100644
index 00000000..62e47541
--- /dev/null
+++ b/gi/pf/README
@@ -0,0 +1,2 @@
+Experimental Bayesian alignment tools. Nothing to see here.
+
diff --git a/gi/pf/base_measures.cc b/gi/pf/base_measures.cc
new file mode 100644
index 00000000..f8ddfd32
--- /dev/null
+++ b/gi/pf/base_measures.cc
@@ -0,0 +1,112 @@
+#include "base_measures.h"
+
+#include <iostream>
+
+#include "filelib.h"
+
+using namespace std;
+
+void Model1::LoadModel1(const string& fname) {
+  cerr << "Loading Model 1 parameters from " << fname << " ..." << endl;
+  ReadFile rf(fname);
+  istream& in = *rf.stream();
+  string line;
+  unsigned lc = 0;
+  while(getline(in, line)) {
+    ++lc;
+    int cur = 0;
+    int start = 0;
+    while(cur < line.size() && line[cur] != ' ') { ++cur; }
+    assert(cur != line.size());
+    line[cur] = 0;
+    const WordID src = TD::Convert(&line[0]);
+    ++cur;
+    start = cur;
+    while(cur < line.size() && line[cur] != ' ') { ++cur; }
+    assert(cur != line.size());
+    line[cur] = 0;
+    WordID trg = TD::Convert(&line[start]);
+    const double logprob = strtod(&line[cur + 1], NULL);
+    if (src >= ttable.size()) ttable.resize(src + 1);
+    ttable[src][trg].logeq(logprob);
+  }
+  cerr << "  read " << lc << " parameters.\n";
+}
+
+prob_t PhraseConditionalBase::p0(const vector<WordID>& vsrc,
+                                 const vector<WordID>& vtrg,
+                                 int start_src, int start_trg) const {
+  const int flen = vsrc.size() - start_src;
+  const int elen = vtrg.size() - start_trg;
+  prob_t uniform_src_alignment; uniform_src_alignment.logeq(-log(flen + 1));
+  prob_t p;
+  p.logeq(log_poisson(elen, flen + 0.01));       // elen | flen          ~Pois(flen + 0.01)
+  for (int i = 0; i < elen; ++i) {               // for each position i in e-RHS
+    const WordID trg = vtrg[i + start_trg];
+    prob_t tp = prob_t::Zero();
+    for (int j = -1; j < flen; ++j) {
+      const WordID src = j < 0 ? 0 : vsrc[j + start_src];
+      tp += kM1MIXTURE * model1(src, trg);
+      tp += kUNIFORM_MIXTURE * kUNIFORM_TARGET;
+    }
+    tp *= uniform_src_alignment;                 //     draw a_i         ~uniform
+    p *= tp;                                     //     draw e_i         ~Model1(f_a_i) / uniform
+  }
+  if (p.is_0()) {
+    cerr << "Zero! " << vsrc << "\nTRG=" << vtrg << endl;
+    abort();
+  }
+  return p;
+}
+
+prob_t PhraseJointBase::p0(const vector<WordID>& vsrc,
+                           const vector<WordID>& vtrg,
+                           int start_src, int start_trg) const {
+  const int flen = vsrc.size() - start_src;
+  const int elen = vtrg.size() - start_trg;
+  prob_t uniform_src_alignment; uniform_src_alignment.logeq(-log(flen + 1));
+  prob_t p;
+  p.logeq(log_poisson(flen, 1.0));               // flen                 ~Pois(1)
+                                                 // elen | flen          ~Pois(flen + 0.01)
+  prob_t ptrglen; ptrglen.logeq(log_poisson(elen, flen + 0.01));
+  p *= ptrglen;
+  p *= kUNIFORM_SOURCE.pow(flen);                // each f in F ~Uniform
+  for (int i = 0; i < elen; ++i) {               // for each position i in E
+    const WordID trg = vtrg[i + start_trg];
+    prob_t tp = prob_t::Zero();
+    for (int j = -1; j < flen; ++j) {
+      const WordID src = j < 0 ? 0 : vsrc[j + start_src];
+      tp += kM1MIXTURE * model1(src, trg);
+      tp += kUNIFORM_MIXTURE * kUNIFORM_TARGET;
+    }
+    tp *= uniform_src_alignment;                 //     draw a_i         ~uniform
+    p *= tp;                                     //     draw e_i         ~Model1(f_a_i) / uniform
+  }
+  if (p.is_0()) {
+    cerr << "Zero! " << vsrc << "\nTRG=" << vtrg << endl;
+    abort();
+  }
+  return p;
+}
+
+JumpBase::JumpBase() : p(200) {
+  for (unsigned src_len = 1; src_len < 200; ++src_len) {
+    map<int, prob_t>& cpd = p[src_len];
+    int min_jump = 1 - src_len;
+    int max_jump = src_len;
+    prob_t z;
+    for (int j = min_jump; j <= max_jump; ++j) {
+      prob_t& cp = cpd[j];
+      if (j < 0)
+        cp.logeq(log_poisson(1.5-j, 1));
+      else if (j > 0)
+        cp.logeq(log_poisson(j, 1));
+      cp.poweq(0.2);
+      z += cp;
+    }
+    for (int j = min_jump; j <= max_jump; ++j) {
+      cpd[j] /= z;
+    }
+  }
+}
+
diff --git a/gi/pf/base_measures.h b/gi/pf/base_measures.h
new file mode 100644
index 00000000..df17aa62
--- /dev/null
+++ b/gi/pf/base_measures.h
@@ -0,0 +1,116 @@
+#ifndef _BASE_MEASURES_H_
+#define _BASE_MEASURES_H_
+
+#include <vector>
+#include <map>
+#include <string>
+#include <cmath>
+#include <iostream>
+
+#include "trule.h"
+#include "prob.h"
+#include "tdict.h"
+
+inline double log_poisson(unsigned x, const double& lambda) {
+  assert(lambda > 0.0);
+  return log(lambda) * x - lgamma(x + 1) - lambda;
+}
+
+inline std::ostream& operator<<(std::ostream& os, const std::vector<WordID>& p) {
+  os << '[';
+  for (int i = 0; i < p.size(); ++i)
+    os << (i==0 ? "" : " ") << TD::Convert(p[i]);
+  return os << ']';
+}
+
+struct Model1 {
+  explicit Model1(const std::string& fname) :
+      kNULL(TD::Convert("<eps>")),
+      kZERO() {
+    LoadModel1(fname);
+  }
+
+  void LoadModel1(const std::string& fname);
+
+  // returns prob 0 if src or trg is not found
+  const prob_t& operator()(WordID src, WordID trg) const {
+    if (src == 0) src = kNULL;
+    if (src < ttable.size()) {
+      const std::map<WordID, prob_t>& cpd = ttable[src];
+      const std::map<WordID, prob_t>::const_iterator it = cpd.find(trg);
+      if (it != cpd.end())
+        return it->second;
+    }
+    return kZERO;
+  }
+
+  const WordID kNULL;
+  const prob_t kZERO;
+  std::vector<std::map<WordID, prob_t> > ttable;
+};
+
+struct PhraseConditionalBase {
+  explicit PhraseConditionalBase(const Model1& m1, const double m1mixture, const unsigned vocab_e_size) :
+      model1(m1),
+      kM1MIXTURE(m1mixture),
+      kUNIFORM_MIXTURE(1.0 - m1mixture),
+      kUNIFORM_TARGET(1.0 / vocab_e_size) {
+    assert(m1mixture >= 0.0 && m1mixture <= 1.0);
+    assert(vocab_e_size > 0);
+  }
+
+  // return p0 of rule.e_ | rule.f_
+  prob_t operator()(const TRule& rule) const {
+    return p0(rule.f_, rule.e_, 0, 0);
+  }
+
+  prob_t p0(const std::vector<WordID>& vsrc, const std::vector<WordID>& vtrg, int start_src, int start_trg) const;
+
+  const Model1& model1;
+  const prob_t kM1MIXTURE;  // Model 1 mixture component
+  const prob_t kUNIFORM_MIXTURE; // uniform mixture component
+  const prob_t kUNIFORM_TARGET;
+};
+
+struct PhraseJointBase {
+  explicit PhraseJointBase(const Model1& m1, const double m1mixture, const unsigned vocab_e_size, const unsigned vocab_f_size) :
+      model1(m1),
+      kM1MIXTURE(m1mixture),
+      kUNIFORM_MIXTURE(1.0 - m1mixture),
+      kUNIFORM_SOURCE(1.0 / vocab_f_size),
+      kUNIFORM_TARGET(1.0 / vocab_e_size) {
+    assert(m1mixture >= 0.0 && m1mixture <= 1.0);
+    assert(vocab_e_size > 0);
+  }
+
+  // return p0 of rule.e_ | rule.f_
+  prob_t operator()(const TRule& rule) const {
+    return p0(rule.f_, rule.e_, 0, 0);
+  }
+
+  prob_t p0(const std::vector<WordID>& vsrc, const std::vector<WordID>& vtrg, int start_src, int start_trg) const;
+
+  const Model1& model1;
+  const prob_t kM1MIXTURE;  // Model 1 mixture component
+  const prob_t kUNIFORM_MIXTURE; // uniform mixture component
+  const prob_t kUNIFORM_SOURCE;
+  const prob_t kUNIFORM_TARGET;
+};
+
+// base distribution for jump size multinomials
+// basically p(0) = 0 and then, p(1) is max, and then
+// you drop as you move to the max jump distance
+struct JumpBase {
+  JumpBase();
+
+  const prob_t& operator()(int jump, unsigned src_len) const {
+    assert(jump != 0);
+    const std::map<int, prob_t>::const_iterator it = p[src_len].find(jump);
+    assert(it != p[src_len].end());
+    return it->second;
+  }
+  std::vector<std::map<int, prob_t> > p;
+};
+
+
+#endif
diff --git a/gi/pf/brat.cc b/gi/pf/brat.cc
new file mode 100644
index 00000000..4c6ba3ef
--- /dev/null
+++ b/gi/pf/brat.cc
@@ -0,0 +1,554 @@
+#include <iostream>
+#include <tr1/memory>
+#include <queue>
+
+#include <boost/functional.hpp>
+#include <boost/multi_array.hpp>
+#include <boost/program_options.hpp>
+#include <boost/program_options/variables_map.hpp>
+
+#include "viterbi.h"
+#include "hg.h"
+#include "trule.h"
+#include "tdict.h"
+#include "filelib.h"
+#include "dict.h"
+#include "sampler.h"
+#include "ccrp_nt.h"
+#include "cfg_wfst_composer.h"
+
+using namespace std;
+using namespace tr1;
+namespace po = boost::program_options;
+
+static unsigned kMAX_SRC_PHRASE;
+static unsigned kMAX_TRG_PHRASE;
+struct FSTState;
+
+size_t hash_value(const TRule& r) {
+  size_t h = 2 - r.lhs_;
+  boost::hash_combine(h, boost::hash_value(r.e_));
+  boost::hash_combine(h, boost::hash_value(r.f_));
+  return h;
+}
+
+bool operator==(const TRule& a, const TRule& b) {
+  return (a.lhs_ == b.lhs_ && a.e_ == b.e_ && a.f_ == b.f_);
+}
+
+double log_poisson(unsigned x, const double& lambda) {
+  assert(lambda > 0.0);
+  return log(lambda) * x - lgamma(x + 1) - lambda;
+}
+
+struct ConditionalBase {
+  explicit ConditionalBase(const double m1mixture, const unsigned vocab_e_size, const string& model1fname) :
+      kM1MIXTURE(m1mixture),
+      kUNIFORM_MIXTURE(1.0 - m1mixture),
+      kUNIFORM_TARGET(1.0 / vocab_e_size),
+      kNULL(TD::Convert("<eps>")) {
+    assert(m1mixture >= 0.0 && m1mixture <= 1.0);
+    assert(vocab_e_size > 0);
+    LoadModel1(model1fname);
+  }
+
+  void LoadModel1(const string& fname) {
+    cerr << "Loading Model 1 parameters from " << fname << " ..." << endl;
+    ReadFile rf(fname);
+    istream& in = *rf.stream();
+    string line;
+    unsigned lc = 0;
+    while(getline(in, line)) {
+      ++lc;
+      int cur = 0;
+      int start = 0;
+      while(cur < line.size() && line[cur] != ' ') { ++cur; }
+      assert(cur != line.size());
+      line[cur] = 0;
+      const WordID src = TD::Convert(&line[0]);
+      ++cur;
+      start = cur;
+      while(cur < line.size() && line[cur] != ' ') { ++cur; }
+      assert(cur != line.size());
+      line[cur] = 0;
+      WordID trg = TD::Convert(&line[start]);
+      const double logprob = strtod(&line[cur + 1], NULL);
+      if (src >= ttable.size()) ttable.resize(src + 1);
+      ttable[src][trg].logeq(logprob);
+    }
+    cerr << "  read " << lc << " parameters.\n";
+  }
+
+  // return logp0 of rule.e_ | rule.f_
+  prob_t operator()(const TRule& rule) const {
+    const int flen = rule.f_.size();
+    const int elen = rule.e_.size();
+    prob_t uniform_src_alignment; uniform_src_alignment.logeq(-log(flen + 1));
+    prob_t p;
+    p.logeq(log_poisson(elen, flen + 0.01));       // elen | flen          ~Pois(flen + 0.01)
+    for (int i = 0; i < elen; ++i) {               // for each position i in e-RHS
+      const WordID trg = rule.e_[i];
+      prob_t tp = prob_t::Zero();
+      for (int j = -1; j < flen; ++j) {
+        const WordID src = j < 0 ? kNULL : rule.f_[j];
+        const map<WordID, prob_t>::const_iterator it = ttable[src].find(trg);
+        if (it != ttable[src].end()) {
+          tp += kM1MIXTURE * it->second;
+        }
+        tp += kUNIFORM_MIXTURE * kUNIFORM_TARGET;
+      }
+      tp *= uniform_src_alignment;                 //     draw a_i         ~uniform
+      p *= tp;                                     //     draw e_i         ~Model1(f_a_i) / uniform
+    }
+    return p;
+  }
+
+  const prob_t kM1MIXTURE;  // Model 1 mixture component
+  const prob_t kUNIFORM_MIXTURE; // uniform mixture component
+  const prob_t kUNIFORM_TARGET;
+  const WordID kNULL;
+  vector<map<WordID, prob_t> > ttable;
+};
+
+void InitCommandLine(int argc, char** argv, po::variables_map* conf) {
+  po::options_description opts("Configuration options");
+  opts.add_options()
+        ("samples,s",po::value<unsigned>()->default_value(1000),"Number of samples")
+        ("input,i",po::value<string>(),"Read parallel data from")
+        ("max_src_phrase",po::value<unsigned>()->default_value(3),"Maximum length of source language phrases")
+        ("max_trg_phrase",po::value<unsigned>()->default_value(3),"Maximum length of target language phrases")
+        ("model1,m",po::value<string>(),"Model 1 parameters (used in base distribution)")
+        ("model1_interpolation_weight",po::value<double>()->default_value(0.95),"Mixing proportion of model 1 with uniform target distribution")
+        ("random_seed,S",po::value<uint32_t>(), "Random seed");
+  po::options_description clo("Command line options");
+  clo.add_options()
+        ("config", po::value<string>(), "Configuration file")
+        ("help,h", "Print this help message and exit");
+  po::options_description dconfig_options, dcmdline_options;
+  dconfig_options.add(opts);
+  dcmdline_options.add(opts).add(clo);
+  
+  po::store(parse_command_line(argc, argv, dcmdline_options), *conf);
+  if (conf->count("config")) {
+    ifstream config((*conf)["config"].as<string>().c_str());
+    po::store(po::parse_config_file(config, dconfig_options), *conf);
+  }
+  po::notify(*conf);
+
+  if (conf->count("help") || (conf->count("input") == 0)) {
+    cerr << dcmdline_options << endl;
+    exit(1);
+  }
+}
+
+void ReadParallelCorpus(const string& filename,
+                vector<vector<WordID> >* f,
+                vector<vector<int> >* e,
+                set<int>* vocab_f,
+                set<int>* vocab_e) {
+  f->clear();
+  e->clear();
+  vocab_f->clear();
+  vocab_e->clear();
+  istream* in;
+  if (filename == "-")
+    in = &cin;
+  else
+    in = new ifstream(filename.c_str());
+  assert(*in);
+  string line;
+  const WordID kDIV = TD::Convert("|||");
+  vector<WordID> tmp;
+  while(*in) {
+    getline(*in, line);
+    if (line.empty() && !*in) break;
+    e->push_back(vector<int>());
+    f->push_back(vector<int>());
+    vector<int>& le = e->back();
+    vector<int>& lf = f->back();
+    tmp.clear();
+    TD::ConvertSentence(line, &tmp);
+    bool isf = true;
+    for (unsigned i = 0; i < tmp.size(); ++i) {
+      const int cur = tmp[i];
+      if (isf) {
+        if (kDIV == cur) { isf = false; } else {
+          lf.push_back(cur);
+          vocab_f->insert(cur);
+        }
+      } else {
+        assert(cur != kDIV);
+        le.push_back(cur);
+        vocab_e->insert(cur);
+      }
+    }
+    assert(isf == false);
+  }
+  if (in != &cin) delete in;
+}
+
+struct UniphraseLM {
+  UniphraseLM(const vector<vector<int> >& corpus,
+              const set<int>& vocab,
+              const po::variables_map& conf) :
+    phrases_(1,1),
+    gen_(1,1),
+    corpus_(corpus),
+    uniform_word_(1.0 / vocab.size()),
+    gen_p0_(0.5),
+    p_end_(0.5),
+    use_poisson_(conf.count("poisson_length") > 0) {}
+
+  void ResampleHyperparameters(MT19937* rng) {
+    phrases_.resample_hyperparameters(rng);
+    gen_.resample_hyperparameters(rng);
+    cerr << " " << phrases_.concentration();
+  }
+
+  CCRP_NoTable<vector<int> > phrases_;
+  CCRP_NoTable<bool> gen_;
+  vector<vector<bool> > z_;   // z_[i] is there a phrase boundary after the ith word
+  const vector<vector<int> >& corpus_;
+  const double uniform_word_;
+  const double gen_p0_;
+  const double p_end_; // in base length distribution, p of the end of a phrase
+  const bool use_poisson_;
+};
+
+struct Reachability {
+  boost::multi_array<bool, 4> edges;  // edges[src_covered][trg_covered][x][trg_delta] is this edge worth exploring?
+  boost::multi_array<short, 2> max_src_delta; // msd[src_covered][trg_covered] -- the largest src delta that's valid
+
+  Reachability(int srclen, int trglen, int src_max_phrase_len, int trg_max_phrase_len) :
+      edges(boost::extents[srclen][trglen][src_max_phrase_len+1][trg_max_phrase_len+1]),
+      max_src_delta(boost::extents[srclen][trglen]) {
+    ComputeReachability(srclen, trglen, src_max_phrase_len, trg_max_phrase_len);
+  }
+
+ private:
+  struct SState {
+    SState() : prev_src_covered(), prev_trg_covered() {}
+    SState(int i, int j) : prev_src_covered(i), prev_trg_covered(j) {}
+    int prev_src_covered;
+    int prev_trg_covered;
+  };
+
+  struct NState {
+    NState() : next_src_covered(), next_trg_covered() {}
+    NState(int i, int j) : next_src_covered(i), next_trg_covered(j) {}
+    int next_src_covered;
+    int next_trg_covered;
+  };
+
+  void ComputeReachability(int srclen, int trglen, int src_max_phrase_len, int trg_max_phrase_len) {
+    typedef boost::multi_array<vector<SState>, 2> array_type;
+    array_type a(boost::extents[srclen + 1][trglen + 1]);
+    a[0][0].push_back(SState());
+    for (int i = 0; i < srclen; ++i) {
+      for (int j = 0; j < trglen; ++j) {
+        if (a[i][j].size() == 0) continue;
+        const SState prev(i,j);
+        for (int k = 1; k <= src_max_phrase_len; ++k) {
+          if ((i + k) > srclen) continue;
+          for (int l = 1; l <= trg_max_phrase_len; ++l) {
+            if ((j + l) > trglen) continue;
+            a[i + k][j + l].push_back(prev);
+          }
+        }
+      }
+    }
+    a[0][0].clear();
+    cerr << "Final cell contains " << a[srclen][trglen].size() << " back pointers\n";
+    assert(a[srclen][trglen].size() > 0);
+
+    typedef boost::multi_array<bool, 2> rarray_type;
+    rarray_type r(boost::extents[srclen + 1][trglen + 1]);
+//    typedef boost::multi_array<vector<NState>, 2> narray_type;
+//    narray_type b(boost::extents[srclen + 1][trglen + 1]);
+    r[srclen][trglen] = true;
+    for (int i = srclen; i >= 0; --i) {
+      for (int j = trglen; j >= 0; --j) {
+        vector<SState>& prevs = a[i][j];
+        if (!r[i][j]) { prevs.clear(); }
+//        const NState nstate(i,j);
+        for (int k = 0; k < prevs.size(); ++k) {
+          r[prevs[k].prev_src_covered][prevs[k].prev_trg_covered] = true;
+          int src_delta = i - prevs[k].prev_src_covered;
+          edges[prevs[k].prev_src_covered][prevs[k].prev_trg_covered][src_delta][j - prevs[k].prev_trg_covered] = true;
+          short &msd = max_src_delta[prevs[k].prev_src_covered][prevs[k].prev_trg_covered];
+          if (src_delta > msd) msd = src_delta;
+//          b[prevs[k].prev_src_covered][prevs[k].prev_trg_covered].push_back(nstate);
+        }
+      }
+    }
+    assert(!edges[0][0][1][0]);
+    assert(!edges[0][0][0][1]);
+    assert(!edges[0][0][0][0]);
+    cerr << "  MAX SRC DELTA[0][0] = " << max_src_delta[0][0] << endl;
+    assert(max_src_delta[0][0] > 0);
+    //cerr << "First cell contains " << b[0][0].size() << " forward pointers\n";
+    //for (int i = 0; i < b[0][0].size(); ++i) {
+    //  cerr << "  -> (" << b[0][0][i].next_src_covered << "," << b[0][0][i].next_trg_covered << ")\n";
+    //}
+  }
+};
+
+ostream& operator<<(ostream& os, const FSTState& q);
+struct FSTState {
+  explicit FSTState(int src_size) :
+      trg_covered_(),
+      src_covered_(),
+      src_coverage_(src_size) {}
+
+  FSTState(short trg_covered, short src_covered, const vector<bool>& src_coverage, const vector<short>& src_prefix) :
+      trg_covered_(trg_covered),
+      src_covered_(src_covered),
+      src_coverage_(src_coverage),
+      src_prefix_(src_prefix) {
+    if (src_coverage_.size() == src_covered) {
+      assert(src_prefix.size() == 0);
+    }
+  }
+
+  // if we extend by the word at src_position, what are
+  // the next states that are reachable and lie on a valid
+  // path to the final state?
+  vector<FSTState> Extensions(int src_position, int src_len, int trg_len, const Reachability& r) const {
+    assert(src_position < src_coverage_.size());
+    if (src_coverage_[src_position]) {
+      cerr << "Trying to extend " << *this << " with position " << src_position << endl;
+      abort();
+    }
+    vector<bool> ncvg = src_coverage_;
+    ncvg[src_position] = true;
+
+    vector<FSTState> res;
+    const int trg_remaining = trg_len - trg_covered_;
+    if (trg_remaining <= 0) {
+      cerr << "Target appears to have been covered: " << *this << " (trg_len=" << trg_len << ",trg_covered=" << trg_covered_ << ")" << endl;
+      abort();
+    }
+    const int src_remaining = src_len - src_covered_;
+    if (src_remaining <= 0) {
+      cerr << "Source appears to have been covered: " << *this << endl;
+      abort();
+    }
+
+    for (int tc = 1; tc <= kMAX_TRG_PHRASE; ++tc) {
+      if (r.edges[src_covered_][trg_covered_][src_prefix_.size() + 1][tc]) {
+        int nc = src_prefix_.size() + 1 + src_covered_;
+        res.push_back(FSTState(trg_covered_ + tc, nc, ncvg, vector<short>()));
+      }
+    }
+
+    if ((src_prefix_.size() + 1) < r.max_src_delta[src_covered_][trg_covered_]) {
+      vector<short> nsp = src_prefix_;
+      nsp.push_back(src_position);
+      res.push_back(FSTState(trg_covered_, src_covered_, ncvg, nsp));
+    }
+
+    if (res.size() == 0) {
+      cerr << *this << " can't be extended!\n";
+      abort();
+    }
+    return res;
+  }
+
+  short trg_covered_, src_covered_;
+  vector<bool> src_coverage_;
+  vector<short> src_prefix_;
+};
+bool operator<(const FSTState& q, const FSTState& r) {
+  if (q.trg_covered_ != r.trg_covered_) return q.trg_covered_ < r.trg_covered_;
+  if (q.src_covered_!= r.src_covered_) return q.src_covered_ < r.src_covered_;
+  if (q.src_coverage_ != r.src_coverage_) return q.src_coverage_ < r.src_coverage_;
+  return q.src_prefix_ < r.src_prefix_;
+}
+
+ostream& operator<<(ostream& os, const FSTState& q) {
+  os << "[" << q.trg_covered_ << " : ";
+  for (int i = 0; i < q.src_coverage_.size(); ++i)
+    os << q.src_coverage_[i];
+  os << " : <";
+  for (int i = 0; i < q.src_prefix_.size(); ++i) {
+    if (i != 0) os << ' ';
+    os << q.src_prefix_[i];
+  }
+  return os << ">]";
+}
+
+struct MyModel {
+  MyModel(ConditionalBase& rcp0) : rp0(rcp0) {}
+  typedef unordered_map<vector<WordID>, CCRP_NoTable<TRule>, boost::hash<vector<WordID> > > SrcToRuleCRPMap;
+
+  void DecrementRule(const TRule& rule) {
+    SrcToRuleCRPMap::iterator it = rules.find(rule.f_);
+    assert(it != rules.end());
+    it->second.decrement(rule);
+    if (it->second.num_customers() == 0) rules.erase(it);
+  }
+
+  void IncrementRule(const TRule& rule) {
+    SrcToRuleCRPMap::iterator it = rules.find(rule.f_);
+    if (it == rules.end()) {
+      CCRP_NoTable<TRule> crp(1,1);
+      it = rules.insert(make_pair(rule.f_, crp)).first;
+    }
+    it->second.increment(rule);
+  }
+
+  // conditioned on rule.f_
+  prob_t RuleConditionalProbability(const TRule& rule) const {
+    const prob_t base = rp0(rule);
+    SrcToRuleCRPMap::const_iterator it = rules.find(rule.f_);
+    if (it == rules.end()) {
+      return base;
+    } else {
+      const double lp = it->second.logprob(rule, log(base));
+      prob_t q; q.logeq(lp);
+      return q;
+    }
+  }
+
+  const ConditionalBase& rp0;
+  SrcToRuleCRPMap rules;
+};
+
+struct MyFST : public WFST {
+  MyFST(const vector<WordID>& ssrc, const vector<WordID>& strg, MyModel* m) :
+      src(ssrc), trg(strg),
+      r(src.size(),trg.size(),kMAX_SRC_PHRASE, kMAX_TRG_PHRASE),
+      model(m) {
+    FSTState in(src.size());
+    cerr << " INIT: " << in << endl;
+    init = GetNode(in);
+    for (int i = 0; i < in.src_coverage_.size(); ++i) in.src_coverage_[i] = true;
+    in.src_covered_ = src.size();
+    in.trg_covered_ = trg.size();
+    cerr << "FINAL: " << in << endl;
+    final = GetNode(in);
+  }
+  virtual const WFSTNode* Final() const;
+  virtual const WFSTNode* Initial() const;
+
+  const WFSTNode* GetNode(const FSTState& q);
+  map<FSTState, boost::shared_ptr<WFSTNode> > m;
+  const vector<WordID>& src;
+  const vector<WordID>& trg;
+  Reachability r;
+  const WFSTNode* init;
+  const WFSTNode* final;
+  MyModel* model;
+};
+
+struct MyNode : public WFSTNode {
+  MyNode(const FSTState& q, MyFST* fst) : state(q), container(fst) {}
+  virtual vector<pair<const WFSTNode*, TRulePtr> > ExtendInput(unsigned srcindex) const;
+  const FSTState state;
+  mutable MyFST* container;
+};
+
+vector<pair<const WFSTNode*, TRulePtr> > MyNode::ExtendInput(unsigned srcindex) const {
+  cerr << "EXTEND " << state << " with " << srcindex << endl;
+  vector<FSTState> ext = state.Extensions(srcindex, container->src.size(), container->trg.size(), container->r);
+  vector<pair<const WFSTNode*,TRulePtr> > res(ext.size());
+  for (unsigned i = 0; i < ext.size(); ++i) {
+    res[i].first = container->GetNode(ext[i]);
+    if (ext[i].src_prefix_.size() == 0) {
+      const unsigned trg_from = state.trg_covered_;
+      const unsigned trg_to = ext[i].trg_covered_;
+      const unsigned prev_prfx_size = state.src_prefix_.size();
+      res[i].second.reset(new TRule);
+      res[i].second->lhs_ = -TD::Convert("X");
+      vector<WordID>& src = res[i].second->f_;
+      vector<WordID>& trg = res[i].second->e_;
+      src.resize(prev_prfx_size + 1);
+      for (unsigned j = 0; j < prev_prfx_size; ++j)
+        src[j] = container->src[state.src_prefix_[j]];
+      src[prev_prfx_size] = container->src[srcindex];
+      for (unsigned j = trg_from; j < trg_to; ++j)
+        trg.push_back(container->trg[j]);
+      res[i].second->scores_.set_value(FD::Convert("Proposal"), log(container->model->RuleConditionalProbability(*res[i].second)));
+    }
+  }
+  return res;
+}
+
+const WFSTNode* MyFST::GetNode(const FSTState& q) {
+  boost::shared_ptr<WFSTNode>& res = m[q];
+  if (!res) {
+    res.reset(new MyNode(q, this));
+  }
+  return &*res;
+}
+
+const WFSTNode* MyFST::Final() const {
+  return final;
+}
+
+const WFSTNode* MyFST::Initial() const {
+  return init;
+}
+
+int main(int argc, char** argv) {
+  po::variables_map conf;
+  InitCommandLine(argc, argv, &conf);
+  kMAX_TRG_PHRASE = conf["max_trg_phrase"].as<unsigned>();
+  kMAX_SRC_PHRASE = conf["max_src_phrase"].as<unsigned>();
+
+  if (!conf.count("model1")) {
+    cerr << argv[0] << "Please use --model1 to specify model 1 parameters\n";
+    return 1;
+  }
+  shared_ptr<MT19937> prng;
+  if (conf.count("random_seed"))
+    prng.reset(new MT19937(conf["random_seed"].as<uint32_t>()));
+  else
+    prng.reset(new MT19937);
+  MT19937& rng = *prng;
+
+  vector<vector<int> > corpuse, corpusf;
+  set<int> vocabe, vocabf;
+  ReadParallelCorpus(conf["input"].as<string>(), &corpusf, &corpuse, &vocabf, &vocabe);
+  cerr << "f-Corpus size: " << corpusf.size() << " sentences\n";
+  cerr << "f-Vocabulary size: " << vocabf.size() << " types\n";
+  cerr << "f-Corpus size: " << corpuse.size() << " sentences\n";
+  cerr << "f-Vocabulary size: " << vocabe.size() << " types\n";
+  assert(corpusf.size() == corpuse.size());
+
+  ConditionalBase lp0(conf["model1_interpolation_weight"].as<double>(),
+                      vocabe.size(),
+                      conf["model1"].as<string>());
+  MyModel m(lp0);
+
+  TRule x("[X] ||| kAnwntR myN ||| at the convent ||| 0");
+  m.IncrementRule(x);
+  TRule y("[X] ||| nY dyN ||| gave ||| 0");
+  m.IncrementRule(y);
+
+
+  MyFST fst(corpusf[0], corpuse[0], &m);
+  ifstream in("./kimura.g");
+  assert(in);
+  CFG_WFSTComposer comp(fst);
+  Hypergraph hg;
+  bool succeed = comp.Compose(&in, &hg);
+  hg.PrintGraphviz();
+  if (succeed) { cerr << "SUCCESS.\n"; } else { cerr << "FAILURE REPORTED.\n"; }
+
+#if 0
+  ifstream in2("./amnabooks.g");
+  assert(in2);
+  MyFST fst2(corpusf[1], corpuse[1], &m);
+  CFG_WFSTComposer comp2(fst2);
+  Hypergraph hg2;
+  bool succeed2 = comp2.Compose(&in2, &hg2);
+  if (succeed2) { cerr << "SUCCESS.\n"; } else { cerr << "FAILURE REPORTED.\n"; }
+#endif
+
+  SparseVector<double> w; w.set_value(FD::Convert("Proposal"), 1.0);
+  hg.Reweight(w);
+  cerr << ViterbiFTree(hg) << endl;
+  return 0;
+}
+
diff --git a/gi/pf/cbgi.cc b/gi/pf/cbgi.cc
new file mode 100644
index 00000000..20204e8a
--- /dev/null
+++ b/gi/pf/cbgi.cc
@@ -0,0 +1,340 @@
+#include <queue>
+#include <sstream>
+#include <iostream>
+
+#include <boost/unordered_map.hpp>
+#include <boost/functional/hash.hpp>
+
+#include "sampler.h"
+#include "filelib.h"
+#include "hg_io.h"
+#include "hg.h"
+#include "ccrp_nt.h"
+#include "trule.h"
+#include "inside_outside.h"
+
+using namespace std;
+using namespace std::tr1;
+
+double log_poisson(unsigned x, const double& lambda) {
+  assert(lambda > 0.0);
+  return log(lambda) * x - lgamma(x + 1) - lambda;
+}
+
+double log_decay(unsigned x, const double& b) {
+  assert(b > 1.0);
+  assert(x > 0);
+  return log(b - 1) - x * log(b);
+}
+
+size_t hash_value(const TRule& r) {
+  // TODO fix hash function
+  size_t h = boost::hash_value(r.e_) * boost::hash_value(r.f_) * r.lhs_;
+  return h;
+}
+
+bool operator==(const TRule& a, const TRule& b) {
+  return (a.lhs_ == b.lhs_ && a.e_ == b.e_ && a.f_ == b.f_);
+}
+
+struct SimpleBase {
+  SimpleBase(unsigned esize, unsigned fsize, unsigned ntsize = 144) :
+    uniform_e(-log(esize)),
+    uniform_f(-log(fsize)),
+    uniform_nt(-log(ntsize)) {
+  }
+
+  // binomial coefficient
+  static double choose(unsigned n, unsigned k) {
+    return exp(lgamma(n + 1) - lgamma(k + 1) - lgamma(n - k + 1));
+  }
+
+  // count the number of patterns of terminals and NTs in the rule, given elen and flen
+  static double log_number_of_patterns(const unsigned flen, const unsigned elen) {
+    static vector<vector<double> > counts;
+    if (elen >= counts.size()) counts.resize(elen + 1);
+    if (flen >= counts[elen].size()) counts[elen].resize(flen + 1);
+    double& count = counts[elen][flen];
+    if (count) return log(count);
+    const unsigned max_arity = min(elen, flen);
+    for (unsigned a = 0; a <= max_arity; ++a)
+      count += choose(elen, a) * choose(flen, a);
+    return log(count);
+  }
+
+  // return logp0 of rule | LHS
+  double operator()(const TRule& rule) const {
+    const unsigned flen = rule.f_.size();
+    const unsigned elen = rule.e_.size();
+#if 0
+    double p = 0;
+    p += log_poisson(flen, 0.5);                   // flen                 ~Pois(0.5)
+    p += log_poisson(elen, flen);                  // elen | flen          ~Pois(flen)
+    p -= log_number_of_patterns(flen, elen);       // pattern | flen,elen  ~Uniform
+    for (unsigned i = 0; i < flen; ++i) {          // for each position in f-RHS
+      if (rule.f_[i] <= 0)                         //   according to pattern
+        p += uniform_nt;                           //     draw NT          ~Uniform
+      else
+        p += uniform_f;                            //     draw f terminal  ~Uniform
+    }
+    p -= lgamma(rule.Arity() + 1);                 // draw permutation     ~Uniform 
+    for (unsigned i = 0; i < elen; ++i) {          // for each position in e-RHS
+      if (rule.e_[i] > 0)                          //   according to pattern
+        p += uniform_e;                            //     draw e|f term    ~Uniform
+        // TODO this should prob be model 1
+    }
+#else
+    double p = 0;
+    bool is_abstract = rule.f_[0] <= 0;
+    p += log(0.5);
+    if (is_abstract) {
+      if (flen == 2) p += log(0.99); else p += log(0.01);
+    } else {
+      p += log_decay(flen, 3);
+    }
+
+    for (unsigned i = 0; i < flen; ++i) {          // for each position in f-RHS
+      if (rule.f_[i] <= 0)                         //   according to pattern
+        p += uniform_nt;                           //     draw NT          ~Uniform
+      else
+        p += uniform_f;                            //     draw f terminal  ~Uniform
+    }
+#endif
+    return p;
+  }
+  const double uniform_e;
+  const double uniform_f;
+  const double uniform_nt;
+  vector<double> arities;
+};
+
+MT19937* rng = NULL;
+
+template <typename Base>
+struct MHSamplerEdgeProb {
+  MHSamplerEdgeProb(const Hypergraph& hg,
+                  const map<int, CCRP_NoTable<TRule> >& rdp,
+                  const Base& logp0,
+                  const bool exclude_multiword_terminals) : edge_probs(hg.edges_.size()) {
+    for (int i = 0; i < edge_probs.size(); ++i) {
+      const TRule& rule = *hg.edges_[i].rule_;
+      const map<int, CCRP_NoTable<TRule> >::const_iterator it = rdp.find(rule.lhs_);
+      assert(it != rdp.end());
+      const CCRP_NoTable<TRule>& crp = it->second;
+      edge_probs[i].logeq(crp.logprob(rule, logp0(rule)));
+      if (exclude_multiword_terminals && rule.f_[0] > 0 && rule.f_.size() > 1)
+        edge_probs[i] = prob_t::Zero();
+    }
+  }
+  inline prob_t operator()(const Hypergraph::Edge& e) const {
+    return edge_probs[e.id_];
+  }
+  prob_t DerivationProb(const vector<int>& d) const {
+    prob_t p = prob_t::One();
+    for (unsigned i = 0; i < d.size(); ++i)
+      p *= edge_probs[d[i]];
+    return p;
+  }
+  vector<prob_t> edge_probs;
+};
+
+template <typename Base>
+struct ModelAndData {
+  ModelAndData() :
+     base_lh(prob_t::One()),
+     logp0(10000, 10000),
+     mh_samples(),
+     mh_rejects() {}
+
+  void SampleCorpus(const string& hgpath, int i);
+  void ResampleHyperparameters() {
+    for (map<int, CCRP_NoTable<TRule> >::iterator it = rules.begin(); it != rules.end(); ++it)
+      it->second.resample_hyperparameters(rng);
+  }
+
+  CCRP_NoTable<TRule>& RuleCRP(int lhs) {
+    map<int, CCRP_NoTable<TRule> >::iterator it = rules.find(lhs);
+    if (it == rules.end()) {
+      rules.insert(make_pair(lhs, CCRP_NoTable<TRule>(1,1)));
+      it = rules.find(lhs);
+    }
+    return it->second;
+  }
+
+  void IncrementRule(const TRule& rule) {
+    CCRP_NoTable<TRule>& crp = RuleCRP(rule.lhs_);
+    if (crp.increment(rule)) {
+      prob_t p; p.logeq(logp0(rule));
+      base_lh *= p;
+    }
+  }
+
+  void DecrementRule(const TRule& rule) {
+    CCRP_NoTable<TRule>& crp = RuleCRP(rule.lhs_);
+    if (crp.decrement(rule)) {
+      prob_t p; p.logeq(logp0(rule));
+      base_lh /= p;
+    }
+  }
+
+  void DecrementDerivation(const Hypergraph& hg, const vector<int>& d) {
+    for (unsigned i = 0; i < d.size(); ++i) {
+      const TRule& rule = *hg.edges_[d[i]].rule_;
+      DecrementRule(rule);
+    }
+  }
+
+  void IncrementDerivation(const Hypergraph& hg, const vector<int>& d) {
+    for (unsigned i = 0; i < d.size(); ++i) {
+      const TRule& rule = *hg.edges_[d[i]].rule_;
+      IncrementRule(rule);
+    }
+  }
+
+  prob_t Likelihood() const {
+    prob_t p = prob_t::One();
+    for (map<int, CCRP_NoTable<TRule> >::const_iterator it = rules.begin(); it != rules.end(); ++it) {
+      prob_t q; q.logeq(it->second.log_crp_prob());
+      p *= q;
+    }
+    p *= base_lh;
+    return p;
+  }
+
+  void ResampleDerivation(const Hypergraph& hg, vector<int>* sampled_derivation);
+
+  map<int, CCRP_NoTable<TRule> > rules;  // [lhs] -> distribution over RHSs
+  prob_t base_lh;
+  SimpleBase logp0;
+  vector<vector<int> > samples;   // sampled derivations
+  unsigned int mh_samples;
+  unsigned int mh_rejects;
+};
+
+template <typename Base>
+void ModelAndData<Base>::SampleCorpus(const string& hgpath, int n) {
+  vector<Hypergraph> hgs(n); hgs.clear();
+  boost::unordered_map<TRule, unsigned> acc;
+  map<int, unsigned> tot;
+  for (int i = 0; i < n; ++i) {
+    ostringstream os;
+    os << hgpath << '/' << i << ".json.gz";
+    if (!FileExists(os.str())) continue;
+    hgs.push_back(Hypergraph());
+    ReadFile rf(os.str());
+    HypergraphIO::ReadFromJSON(rf.stream(), &hgs.back());
+  }
+  cerr << "Read " << hgs.size() << " alignment hypergraphs.\n";
+  samples.resize(hgs.size());
+  const unsigned SAMPLES = 2000;
+  const unsigned burnin = 3 * SAMPLES / 4;
+  const unsigned every = 20;
+  for (unsigned s = 0; s < SAMPLES; ++s) {
+    if (s % 10 == 0) {
+      if (s > 0) { cerr << endl; ResampleHyperparameters(); }
+      cerr << "[" << s << " LLH=" << log(Likelihood()) << " REJECTS=" << ((double)mh_rejects / mh_samples) << " LHS's=" << rules.size() << " base=" << log(base_lh) << "] ";
+    }
+    cerr << '.';
+    for (unsigned i = 0; i < hgs.size(); ++i) {
+      ResampleDerivation(hgs[i], &samples[i]);
+      if (s > burnin && s % every == 0) {
+        for (unsigned j = 0; j < samples[i].size(); ++j) {
+          const TRule& rule = *hgs[i].edges_[samples[i][j]].rule_;
+          ++acc[rule];
+          ++tot[rule.lhs_];
+        }
+      }
+    }
+  }
+  cerr << endl;
+  for (boost::unordered_map<TRule,unsigned>::iterator it = acc.begin(); it != acc.end(); ++it) {
+    cout << it->first << " MyProb=" << log(it->second)-log(tot[it->first.lhs_]) << endl;
+  }
+}
+
+template <typename Base>
+void ModelAndData<Base>::ResampleDerivation(const Hypergraph& hg, vector<int>* sampled_deriv) {
+  vector<int> cur;
+  cur.swap(*sampled_deriv);
+
+  const prob_t p_cur = Likelihood();
+  DecrementDerivation(hg, cur);
+  if (cur.empty()) {
+    // first iteration, create restaurants
+    for (int i = 0; i < hg.edges_.size(); ++i)
+      RuleCRP(hg.edges_[i].rule_->lhs_);
+  }
+  MHSamplerEdgeProb<SimpleBase> wf(hg, rules, logp0, cur.empty());
+//  MHSamplerEdgeProb<SimpleBase> wf(hg, rules, logp0, false);
+  const prob_t q_cur = wf.DerivationProb(cur);
+  vector<prob_t> node_probs;
+  Inside<prob_t, MHSamplerEdgeProb<SimpleBase> >(hg, &node_probs, wf);
+  queue<unsigned> q;
+  q.push(hg.nodes_.size() - 3);
+  while(!q.empty()) {
+    unsigned cur_node_id = q.front();
+//    cerr << "NODE=" << cur_node_id << endl;
+    q.pop();
+    const Hypergraph::Node& node = hg.nodes_[cur_node_id];
+    const unsigned num_in_edges = node.in_edges_.size();
+    unsigned sampled_edge = 0;
+    if (num_in_edges == 1) {
+      sampled_edge = node.in_edges_[0];
+    } else {
+      prob_t z;
+      assert(num_in_edges > 1);
+      SampleSet<prob_t> ss;
+      for (unsigned j = 0; j < num_in_edges; ++j) {
+        const Hypergraph::Edge& edge = hg.edges_[node.in_edges_[j]];
+        prob_t p = wf.edge_probs[edge.id_];             // edge proposal prob
+        for (unsigned k = 0; k < edge.tail_nodes_.size(); ++k)
+          p *= node_probs[edge.tail_nodes_[k]];
+        ss.add(p);
+//        cerr << log(ss[j]) << " ||| " << edge.rule_->AsString() << endl;
+        z += p;
+      }
+//      for (unsigned j = 0; j < num_in_edges; ++j) {
+//        const Hypergraph::Edge& edge = hg.edges_[node.in_edges_[j]];
+//        cerr << exp(log(ss[j] / z)) << " ||| " << edge.rule_->AsString() << endl;
+//      }
+//      cerr << " --- \n";
+      sampled_edge = node.in_edges_[rng->SelectSample(ss)];
+    }
+    sampled_deriv->push_back(sampled_edge);
+    const Hypergraph::Edge& edge = hg.edges_[sampled_edge];
+    for (unsigned j = 0; j < edge.tail_nodes_.size(); ++j) {
+      q.push(edge.tail_nodes_[j]);
+    }
+  }
+  IncrementDerivation(hg, *sampled_deriv);
+
+//  cerr << "sampled derivation contains " << sampled_deriv->size() << " edges\n";
+//  cerr << "DERIV:\n";
+//  for (int i = 0; i < sampled_deriv->size(); ++i) {
+//    cerr << "  " << hg.edges_[(*sampled_deriv)[i]].rule_->AsString() << endl;
+//  }
+
+  if (cur.empty()) return;  // accept first sample
+
+  ++mh_samples;
+  // only need to do MH if proposal is different to current state
+  if (cur != *sampled_deriv) {
+    const prob_t q_prop = wf.DerivationProb(*sampled_deriv);
+    const prob_t p_prop = Likelihood();
+    if (!rng->AcceptMetropolisHastings(p_prop, p_cur, q_prop, q_cur)) {
+      ++mh_rejects;
+      DecrementDerivation(hg, *sampled_deriv);
+      IncrementDerivation(hg, cur);
+      swap(cur, *sampled_deriv);
+    }
+  }
+}
+
+int main(int argc, char** argv) {
+  rng = new MT19937;
+  ModelAndData<SimpleBase> m;
+  m.SampleCorpus("./hgs", 50);
+  // m.SampleCorpus("./btec/hgs", 5000);
+  return 0;
+}
+
diff --git a/gi/pf/cfg_wfst_composer.cc b/gi/pf/cfg_wfst_composer.cc
new file mode 100644
index 00000000..a31b5be8
--- /dev/null
+++ b/gi/pf/cfg_wfst_composer.cc
@@ -0,0 +1,730 @@
+#include "cfg_wfst_composer.h"
+
+#include <iostream>
+#include <fstream>
+#include <map>
+#include <queue>
+#include <tr1/unordered_set>
+
+#include <boost/shared_ptr.hpp>
+#include <boost/program_options.hpp>
+#include <boost/program_options/variables_map.hpp>
+#include "fast_lexical_cast.hpp"
+
+#include "phrasetable_fst.h"
+#include "sparse_vector.h"
+#include "tdict.h"
+#include "hg.h"
+
+using boost::shared_ptr;
+namespace po = boost::program_options;
+using namespace std;
+using namespace std::tr1;
+
+WFSTNode::~WFSTNode() {}
+WFST::~WFST() {}
+
+// Define the following macro if you want to see lots of debugging output
+// when you run the chart parser
+#undef DEBUG_CHART_PARSER
+
+// A few constants used by the chart parser ///////////////
+static const int kMAX_NODES = 2000000;
+static const string kPHRASE_STRING = "X";
+static bool constants_need_init = true;
+static WordID kUNIQUE_START;
+static WordID kPHRASE;
+static TRulePtr kX1X2;
+static TRulePtr kX1;
+static WordID kEPS;
+static TRulePtr kEPSRule;
+
+static void InitializeConstants() {
+  if (constants_need_init) {
+    kPHRASE = TD::Convert(kPHRASE_STRING) * -1;
+    kUNIQUE_START = TD::Convert("S") * -1;
+    kX1X2.reset(new TRule("[X] ||| [X,1] [X,2] ||| [X,1] [X,2]"));
+    kX1.reset(new TRule("[X] ||| [X,1] ||| [X,1]"));
+    kEPSRule.reset(new TRule("[X] ||| <eps> ||| <eps>"));
+    kEPS = TD::Convert("<eps>");
+    constants_need_init = false;
+  }
+}
+////////////////////////////////////////////////////////////
+
+class EGrammarNode {
+  friend bool CFG_WFSTComposer::Compose(const Hypergraph& src_forest, Hypergraph* trg_forest);
+  friend void AddGrammarRule(const string& r, map<WordID, EGrammarNode>* g);
+ public:
+#ifdef DEBUG_CHART_PARSER
+  string hint;
+#endif
+  EGrammarNode() : is_some_rule_complete(false), is_root(false) {}
+  const map<WordID, EGrammarNode>& GetTerminals() const { return tptr; }
+  const map<WordID, EGrammarNode>& GetNonTerminals() const { return ntptr; }
+  bool HasNonTerminals() const { return (!ntptr.empty()); }
+  bool HasTerminals() const { return (!tptr.empty()); }
+  bool RuleCompletes() const {
+    return (is_some_rule_complete || (ntptr.empty() && tptr.empty()));
+  }
+  bool GrammarContinues() const {
+    return !(ntptr.empty() && tptr.empty());
+  }
+  bool IsRoot() const {
+    return is_root;
+  }
+  // these are the features associated with the rule from the start
+  // node up to this point.  If you use these features, you must
+  // not Extend() this rule.
+  const SparseVector<double>& GetCFGProductionFeatures() const {
+    return input_features;
+  }
+
+  const EGrammarNode* Extend(const WordID& t) const {
+    if (t < 0) {
+      map<WordID, EGrammarNode>::const_iterator it = ntptr.find(t);
+      if (it == ntptr.end()) return NULL;
+      return &it->second;
+    } else {
+      map<WordID, EGrammarNode>::const_iterator it = tptr.find(t);
+      if (it == tptr.end()) return NULL;
+      return &it->second;
+    }
+  }
+
+ private:
+  map<WordID, EGrammarNode> tptr;
+  map<WordID, EGrammarNode> ntptr;
+  SparseVector<double> input_features;
+  bool is_some_rule_complete;
+  bool is_root;
+};
+typedef map<WordID, EGrammarNode> EGrammar;    // indexed by the rule LHS
+
+// edges are immutable once created
+struct Edge {
+#ifdef DEBUG_CHART_PARSER
+  static int id_count;
+  const int id;
+#endif
+  const WordID cat;                   // lhs side of rule proved/being proved
+  const EGrammarNode* const dot;      // dot position
+  const WFSTNode* const q;             // start of span
+  const WFSTNode* const r;             // end of span
+  const Edge* const active_parent;    // back pointer, NULL for PREDICT items
+  const Edge* const passive_parent;   // back pointer, NULL for SCAN and PREDICT items
+  TRulePtr tps;   // translations
+  shared_ptr<SparseVector<double> > features; // features from CFG rule
+
+  bool IsPassive() const {
+    // when a rule is completed, this value will be set
+    return static_cast<bool>(features);
+  }
+  bool IsActive() const { return !IsPassive(); }
+  bool IsInitial() const {
+    return !(active_parent || passive_parent);
+  }
+  bool IsCreatedByScan() const {
+    return active_parent && !passive_parent && !dot->IsRoot();
+  }
+  bool IsCreatedByPredict() const {
+    return dot->IsRoot();
+  }
+  bool IsCreatedByComplete() const {
+    return active_parent && passive_parent;
+  }
+
+  // constructor for PREDICT
+  Edge(WordID c, const EGrammarNode* d, const WFSTNode* q_and_r) :
+#ifdef DEBUG_CHART_PARSER
+    id(++id_count),
+#endif
+    cat(c), dot(d), q(q_and_r), r(q_and_r), active_parent(NULL), passive_parent(NULL), tps() {}
+  Edge(WordID c, const EGrammarNode* d, const WFSTNode* q_and_r, const Edge* act_parent) :
+#ifdef DEBUG_CHART_PARSER
+    id(++id_count),
+#endif
+    cat(c), dot(d), q(q_and_r), r(q_and_r), active_parent(act_parent), passive_parent(NULL), tps() {}
+
+  // constructors for SCAN
+  Edge(WordID c, const EGrammarNode* d, const WFSTNode* i, const WFSTNode* j,
+       const Edge* act_par, const TRulePtr& translations) :
+#ifdef DEBUG_CHART_PARSER
+    id(++id_count),
+#endif
+    cat(c), dot(d), q(i), r(j), active_parent(act_par), passive_parent(NULL), tps(translations) {}
+
+  Edge(WordID c, const EGrammarNode* d, const WFSTNode* i, const WFSTNode* j,
+       const Edge* act_par, const TRulePtr& translations,
+       const SparseVector<double>& feats) :
+#ifdef DEBUG_CHART_PARSER
+    id(++id_count),
+#endif
+    cat(c), dot(d), q(i), r(j), active_parent(act_par), passive_parent(NULL), tps(translations),
+    features(new SparseVector<double>(feats)) {}
+
+  // constructors for COMPLETE
+  Edge(WordID c, const EGrammarNode* d, const WFSTNode* i, const WFSTNode* j,
+       const Edge* act_par, const Edge *pas_par) :
+#ifdef DEBUG_CHART_PARSER
+    id(++id_count),
+#endif
+    cat(c), dot(d), q(i), r(j), active_parent(act_par), passive_parent(pas_par), tps() {
+      assert(pas_par->IsPassive());
+      assert(act_par->IsActive());
+    }
+
+  Edge(WordID c, const EGrammarNode* d, const WFSTNode* i, const WFSTNode* j,
+       const Edge* act_par, const Edge *pas_par, const SparseVector<double>& feats) :
+#ifdef DEBUG_CHART_PARSER
+    id(++id_count),
+#endif
+    cat(c), dot(d), q(i), r(j), active_parent(act_par), passive_parent(pas_par), tps(),
+    features(new SparseVector<double>(feats)) {
+      assert(pas_par->IsPassive());
+      assert(act_par->IsActive());
+    }
+
+  // constructor for COMPLETE query
+  Edge(const WFSTNode* _r) :
+#ifdef DEBUG_CHART_PARSER
+    id(0),
+#endif
+    cat(0), dot(NULL), q(NULL),
+    r(_r), active_parent(NULL), passive_parent(NULL), tps() {}
+  // constructor for MERGE quere
+  Edge(const WFSTNode* _q, int) :
+#ifdef DEBUG_CHART_PARSER
+    id(0),
+#endif
+    cat(0), dot(NULL), q(_q),
+    r(NULL), active_parent(NULL), passive_parent(NULL), tps() {}
+};
+#ifdef DEBUG_CHART_PARSER
+int Edge::id_count = 0;
+#endif
+
+ostream& operator<<(ostream& os, const Edge& e) {
+  string type = "PREDICT";
+  if (e.IsCreatedByScan())
+    type = "SCAN";
+  else if (e.IsCreatedByComplete())
+    type = "COMPLETE";
+  os << "["
+#ifdef DEBUG_CHART_PARSER
+     << '(' << e.id << ") "
+#else
+     << '(' << &e << ") "
+#endif
+     << "q=" << e.q << ", r=" << e.r
+     << ", cat="<< TD::Convert(e.cat*-1) << ", dot="
+     << e.dot
+#ifdef DEBUG_CHART_PARSER
+     << e.dot->hint
+#endif
+     << (e.IsActive() ? ", Active" : ", Passive")
+     << ", " << type;
+#ifdef DEBUG_CHART_PARSER
+  if (e.active_parent) { os << ", act.parent=(" << e.active_parent->id << ')'; }
+  if (e.passive_parent) { os << ", psv.parent=(" << e.passive_parent->id << ')'; }
+#endif
+  if (e.tps) { os << ", tps=" << e.tps->AsString(); }
+  return os << ']';
+}
+
+struct Traversal {
+  const Edge* const edge;      // result from the active / passive combination
+  const Edge* const active;
+  const Edge* const passive;
+  Traversal(const Edge* me, const Edge* a, const Edge* p) : edge(me), active(a), passive(p) {}
+};
+
+struct UniqueTraversalHash {
+  size_t operator()(const Traversal* t) const {
+    size_t x = 5381;
+    x = ((x << 5) + x) ^ reinterpret_cast<size_t>(t->active);
+    x = ((x << 5) + x) ^ reinterpret_cast<size_t>(t->passive);
+    x = ((x << 5) + x) ^ t->edge->IsActive();
+    return x;
+  }
+};
+
+struct UniqueTraversalEquals {
+  size_t operator()(const Traversal* a, const Traversal* b) const {
+    return (a->passive == b->passive && a->active == b->active && a->edge->IsActive() == b->edge->IsActive());
+  }
+};
+
+struct UniqueEdgeHash {
+  size_t operator()(const Edge* e) const {
+    size_t x = 5381;
+    if (e->IsActive()) {
+      x = ((x << 5) + x) ^ reinterpret_cast<size_t>(e->dot);
+      x = ((x << 5) + x) ^ reinterpret_cast<size_t>(e->q);
+      x = ((x << 5) + x) ^ reinterpret_cast<size_t>(e->r);
+      x = ((x << 5) + x) ^ static_cast<size_t>(e->cat);
+      x += 13;
+    } else {  // with passive edges, we don't care about the dot
+      x = ((x << 5) + x) ^ reinterpret_cast<size_t>(e->q);
+      x = ((x << 5) + x) ^ reinterpret_cast<size_t>(e->r);
+      x = ((x << 5) + x) ^ static_cast<size_t>(e->cat);
+    }
+    return x;
+  }
+};
+
+struct UniqueEdgeEquals {
+  bool operator()(const Edge* a, const Edge* b) const {
+    if (a->IsActive() != b->IsActive()) return false;
+    if (a->IsActive()) {
+      return (a->cat == b->cat) && (a->dot == b->dot) && (a->q == b->q) && (a->r == b->r);
+    } else {
+      return (a->cat == b->cat) && (a->q == b->q) && (a->r == b->r);
+    }
+  }
+};
+
+struct REdgeHash {
+  size_t operator()(const Edge* e) const {
+    size_t x = 5381;
+    x = ((x << 5) + x) ^ reinterpret_cast<size_t>(e->r);
+    return x;
+  }
+};
+
+struct REdgeEquals {
+  bool operator()(const Edge* a, const Edge* b) const {
+    return (a->r == b->r);
+  }
+};
+
+struct QEdgeHash {
+  size_t operator()(const Edge* e) const {
+    size_t x = 5381;
+    x = ((x << 5) + x) ^ reinterpret_cast<size_t>(e->q);
+    return x;
+  }
+};
+
+struct QEdgeEquals {
+  bool operator()(const Edge* a, const Edge* b) const {
+    return (a->q == b->q);
+  }
+};
+
+struct EdgeQueue {
+  queue<const Edge*> q;
+  EdgeQueue() {}
+  void clear() { while(!q.empty()) q.pop(); }
+  bool HasWork() const { return !q.empty(); }
+  const Edge* Next() { const Edge* res = q.front(); q.pop(); return res; }
+  void AddEdge(const Edge* s) { q.push(s); }
+};
+
+class CFG_WFSTComposerImpl {
+ public:
+  CFG_WFSTComposerImpl(WordID start_cat,
+                       const WFSTNode* q_0,
+                       const WFSTNode* q_final) : start_cat_(start_cat), q_0_(q_0), q_final_(q_final) {}
+
+  // returns false if the intersection is empty
+  bool Compose(const EGrammar& g, Hypergraph* forest) {
+    goal_node = NULL;
+    EGrammar::const_iterator sit = g.find(start_cat_);
+    forest->ReserveNodes(kMAX_NODES);
+    assert(sit != g.end());
+    Edge* init = new Edge(start_cat_, &sit->second, q_0_);
+    assert(IncorporateNewEdge(init));
+    while (exp_agenda.HasWork() || agenda.HasWork()) {
+      while(exp_agenda.HasWork()) {
+        const Edge* edge = exp_agenda.Next();
+        FinishEdge(edge, forest);
+      }
+      if (agenda.HasWork()) {
+        const Edge* edge = agenda.Next();
+#ifdef DEBUG_CHART_PARSER
+        cerr << "processing (" << edge->id << ')' << endl;
+#endif
+        if (edge->IsActive()) {
+          if (edge->dot->HasTerminals())
+            DoScan(edge);
+          if (edge->dot->HasNonTerminals()) {
+            DoMergeWithPassives(edge);
+            DoPredict(edge, g);
+          }
+        } else {
+          DoComplete(edge);
+        }
+      }
+    }
+    if (goal_node) {
+      forest->PruneUnreachable(goal_node->id_);
+      forest->EpsilonRemove(kEPS);
+    }
+    FreeAll();
+    return goal_node;
+  }
+
+  void FreeAll() {
+    for (int i = 0; i < free_list_.size(); ++i)
+      delete free_list_[i];
+    free_list_.clear();
+    for (int i = 0; i < traversal_free_list_.size(); ++i)
+      delete traversal_free_list_[i];
+    traversal_free_list_.clear();
+    all_traversals.clear();
+    exp_agenda.clear();
+    agenda.clear();
+    tps2node.clear();
+    edge2node.clear();
+    all_edges.clear();
+    passive_edges.clear();
+    active_edges.clear();
+  }
+
+  ~CFG_WFSTComposerImpl() {
+    FreeAll();
+  }
+
+  // returns the total number of edges created during composition
+  int EdgesCreated() const {
+    return free_list_.size();
+  }
+
+ private:
+  void DoScan(const Edge* edge) {
+    // here, we assume that the FST will potentially have many more outgoing
+    // edges than the grammar, which will be just a couple.  If you want to
+    // efficiently handle the case where both are relatively large, this code
+    // will need to change how the intersection is done.  The best general
+    // solution would probably be the Baeza-Yates double binary search.
+
+    const EGrammarNode* dot = edge->dot;
+    const WFSTNode* r = edge->r;
+    const map<WordID, EGrammarNode>& terms = dot->GetTerminals();
+    for (map<WordID, EGrammarNode>::const_iterator git = terms.begin();
+         git != terms.end(); ++git) {
+
+      if (!(TD::Convert(git->first)[0] >= '0' && TD::Convert(git->first)[0] <= '9')) {
+        std::cerr << "TERMINAL SYMBOL: " << TD::Convert(git->first) << endl;
+        abort();
+      }
+      std::vector<std::pair<const WFSTNode*, TRulePtr> > extensions = r->ExtendInput(atoi(TD::Convert(git->first)));
+      for (unsigned nsi = 0; nsi < extensions.size(); ++nsi) {
+        const WFSTNode* next_r = extensions[nsi].first;
+        const EGrammarNode* next_dot = &git->second;
+        const bool grammar_continues = next_dot->GrammarContinues();
+        const bool rule_completes    = next_dot->RuleCompletes();
+        if (extensions[nsi].second)
+          cerr << "!!! " << extensions[nsi].second->AsString() << endl;
+        // cerr << "  rule completes: " << rule_completes << " after consuming " << TD::Convert(git->first) << endl;
+        assert(grammar_continues || rule_completes);
+        const SparseVector<double>& input_features = next_dot->GetCFGProductionFeatures();
+        if (rule_completes)
+          IncorporateNewEdge(new Edge(edge->cat, next_dot, edge->q, next_r, edge, extensions[nsi].second, input_features));
+        if (grammar_continues)
+          IncorporateNewEdge(new Edge(edge->cat, next_dot, edge->q, next_r, edge, extensions[nsi].second));
+      }
+    }
+  }
+
+  void DoPredict(const Edge* edge, const EGrammar& g) {
+    const EGrammarNode* dot = edge->dot;
+    const map<WordID, EGrammarNode>& non_terms = dot->GetNonTerminals();
+    for (map<WordID, EGrammarNode>::const_iterator git = non_terms.begin();
+         git != non_terms.end(); ++git) {
+      const WordID nt_to_predict = git->first;
+      //cerr << edge->id << " -- " << TD::Convert(nt_to_predict*-1) << endl;
+      EGrammar::const_iterator egi = g.find(nt_to_predict);
+      if (egi == g.end()) {
+        cerr << "[ERROR] Can't find any grammar rules with a LHS of type "
+             << TD::Convert(-1*nt_to_predict) << '!' << endl;
+        continue;
+      }
+      assert(edge->IsActive());
+      const EGrammarNode* new_dot = &egi->second;
+      Edge* new_edge = new Edge(nt_to_predict, new_dot, edge->r, edge);
+      IncorporateNewEdge(new_edge);
+    }
+  }
+
+  void DoComplete(const Edge* passive) {
+#ifdef DEBUG_CHART_PARSER
+    cerr << "  complete: " << *passive << endl;
+#endif
+    const WordID completed_nt = passive->cat;
+    const WFSTNode* q = passive->q;
+    const WFSTNode* next_r = passive->r;
+    const Edge query(q);
+    const pair<unordered_multiset<const Edge*, REdgeHash, REdgeEquals>::iterator,
+         unordered_multiset<const Edge*, REdgeHash, REdgeEquals>::iterator > p =
+      active_edges.equal_range(&query);
+    for (unordered_multiset<const Edge*, REdgeHash, REdgeEquals>::iterator it = p.first;
+         it != p.second; ++it) {
+      const Edge* active = *it;
+#ifdef DEBUG_CHART_PARSER
+      cerr << "    pos: " << *active << endl;
+#endif
+      const EGrammarNode* next_dot = active->dot->Extend(completed_nt);
+      if (!next_dot) continue;
+      const SparseVector<double>& input_features = next_dot->GetCFGProductionFeatures();
+      // add up to 2 rules
+      if (next_dot->RuleCompletes())
+        IncorporateNewEdge(new Edge(active->cat, next_dot, active->q, next_r, active, passive, input_features));
+      if (next_dot->GrammarContinues())
+        IncorporateNewEdge(new Edge(active->cat, next_dot, active->q, next_r, active, passive));
+    }
+  }
+
+  void DoMergeWithPassives(const Edge* active) {
+    // edge is active, has non-terminals, we need to find the passives that can extend it
+    assert(active->IsActive());
+    assert(active->dot->HasNonTerminals());
+#ifdef DEBUG_CHART_PARSER
+    cerr << "  merge active with passives: ACT=" << *active << endl;
+#endif
+    const Edge query(active->r, 1);
+    const pair<unordered_multiset<const Edge*, QEdgeHash, QEdgeEquals>::iterator,
+         unordered_multiset<const Edge*, QEdgeHash, QEdgeEquals>::iterator > p =
+      passive_edges.equal_range(&query);
+    for (unordered_multiset<const Edge*, QEdgeHash, QEdgeEquals>::iterator it = p.first;
+         it != p.second; ++it) {
+      const Edge* passive = *it;
+      const EGrammarNode* next_dot = active->dot->Extend(passive->cat);
+      if (!next_dot) continue;
+      const WFSTNode* next_r = passive->r;
+      const SparseVector<double>& input_features = next_dot->GetCFGProductionFeatures();
+      if (next_dot->RuleCompletes())
+        IncorporateNewEdge(new Edge(active->cat, next_dot, active->q, next_r, active, passive, input_features));
+      if (next_dot->GrammarContinues())
+        IncorporateNewEdge(new Edge(active->cat, next_dot, active->q, next_r, active, passive));
+    }
+  }
+
+  // take ownership of edge memory, add to various indexes, etc
+  // returns true if this edge is new
+  bool IncorporateNewEdge(Edge* edge) {
+    free_list_.push_back(edge);
+    if (edge->passive_parent && edge->active_parent) {
+      Traversal* t = new Traversal(edge, edge->active_parent, edge->passive_parent);
+      traversal_free_list_.push_back(t);
+      if (all_traversals.find(t) != all_traversals.end()) {
+        return false;
+      } else {
+        all_traversals.insert(t);
+      }
+    }
+    exp_agenda.AddEdge(edge);
+    return true;
+  }
+
+  bool FinishEdge(const Edge* edge, Hypergraph* hg) {
+    bool is_new = false;
+    if (all_edges.find(edge) == all_edges.end()) {
+#ifdef DEBUG_CHART_PARSER
+      cerr << *edge << " is NEW\n";
+#endif
+      all_edges.insert(edge);
+      is_new = true;
+      if (edge->IsPassive()) passive_edges.insert(edge);
+      if (edge->IsActive()) active_edges.insert(edge);
+      agenda.AddEdge(edge);
+    } else {
+#ifdef DEBUG_CHART_PARSER
+      cerr << *edge << " is NOT NEW.\n";
+#endif
+    }
+    AddEdgeToTranslationForest(edge, hg);
+    return is_new;
+  }
+
+  // build the translation forest
+  void AddEdgeToTranslationForest(const Edge* edge, Hypergraph* hg) {
+    assert(hg->nodes_.size() < kMAX_NODES);
+    Hypergraph::Node* tps = NULL;
+    // first add any target language rules
+    if (edge->tps) {
+      Hypergraph::Node*& node = tps2node[(size_t)edge->tps.get()];
+      if (!node) {
+        // cerr << "Creating phrases for " << edge->tps << endl;
+        const TRulePtr& rule = edge->tps;
+        node = hg->AddNode(kPHRASE);
+        Hypergraph::Edge* hg_edge = hg->AddEdge(rule, Hypergraph::TailNodeVector());
+        hg_edge->feature_values_ += rule->GetFeatureValues();
+        hg->ConnectEdgeToHeadNode(hg_edge, node);
+      }
+      tps = node;
+    }
+    Hypergraph::Node*& head_node = edge2node[edge];
+    if (!head_node)
+      head_node = hg->AddNode(kPHRASE);
+    if (edge->cat == start_cat_ && edge->q == q_0_ && edge->r == q_final_ && edge->IsPassive()) {
+      assert(goal_node == NULL || goal_node == head_node);
+      goal_node = head_node;
+    }
+    Hypergraph::TailNodeVector tail;
+    SparseVector<double> extra;
+    if (edge->IsCreatedByPredict()) {
+      // extra.set_value(FD::Convert("predict"), 1);
+    } else if (edge->IsCreatedByScan()) {
+      tail.push_back(edge2node[edge->active_parent]->id_);
+      if (tps) {
+        tail.push_back(tps->id_);
+      }
+      //extra.set_value(FD::Convert("scan"), 1);
+    } else if (edge->IsCreatedByComplete()) {
+      tail.push_back(edge2node[edge->active_parent]->id_);
+      tail.push_back(edge2node[edge->passive_parent]->id_);
+      //extra.set_value(FD::Convert("complete"), 1);
+    } else {
+      assert(!"unexpected edge type!");
+    }
+    //cerr << head_node->id_ << "<--" << *edge << endl;
+
+#ifdef DEBUG_CHART_PARSER
+      for (int i = 0; i < tail.size(); ++i)
+        if (tail[i] == head_node->id_) {
+          cerr << "ERROR: " << *edge << "\n   i=" << i << endl;
+          if (i == 1) { cerr << "\tP: " << *edge->passive_parent << endl; }
+          if (i == 0) { cerr << "\tA: " << *edge->active_parent << endl; }
+          assert(!"self-loop found!");
+        }
+#endif
+    Hypergraph::Edge* hg_edge = NULL;
+    if (tail.size() == 0) {
+      hg_edge = hg->AddEdge(kEPSRule, tail);
+    } else if (tail.size() == 1) {
+      hg_edge = hg->AddEdge(kX1, tail);
+    } else if (tail.size() == 2) {
+      hg_edge = hg->AddEdge(kX1X2, tail);
+    }
+    if (edge->features)
+      hg_edge->feature_values_ += *edge->features;
+    hg_edge->feature_values_ += extra;
+    hg->ConnectEdgeToHeadNode(hg_edge, head_node);
+  }
+
+  Hypergraph::Node* goal_node;
+  EdgeQueue exp_agenda;
+  EdgeQueue agenda;
+  unordered_map<size_t, Hypergraph::Node*> tps2node;
+  unordered_map<const Edge*, Hypergraph::Node*, UniqueEdgeHash, UniqueEdgeEquals> edge2node;
+  unordered_set<const Traversal*, UniqueTraversalHash, UniqueTraversalEquals> all_traversals;
+  unordered_set<const Edge*, UniqueEdgeHash, UniqueEdgeEquals> all_edges;
+  unordered_multiset<const Edge*, QEdgeHash, QEdgeEquals> passive_edges;
+  unordered_multiset<const Edge*, REdgeHash, REdgeEquals> active_edges;
+  vector<Edge*> free_list_;
+  vector<Traversal*> traversal_free_list_;
+  const WordID start_cat_;
+  const WFSTNode* const q_0_;
+  const WFSTNode* const q_final_;
+};
+
+#ifdef DEBUG_CHART_PARSER
+static string TrimRule(const string& r) {
+  size_t start = r.find(" |||") + 5;
+  size_t end = r.rfind(" |||");
+  return r.substr(start, end - start);
+}
+#endif
+
+void AddGrammarRule(const string& r, EGrammar* g) {
+  const size_t pos = r.find(" ||| ");
+  if (pos == string::npos || r[0] != '[') {
+    cerr << "Bad rule: " << r << endl;
+    return;
+  }
+  const size_t rpos = r.rfind(" ||| ");
+  string feats;
+  string rs = r;
+  if (rpos != pos) {
+    feats = r.substr(rpos + 5);
+    rs = r.substr(0, rpos);
+  }
+  string rhs = rs.substr(pos + 5);
+  string trule = rs + " ||| " + rhs + " ||| " + feats;
+  TRule tr(trule);
+  cerr << "X: " << tr.e_[0] << endl;
+#ifdef DEBUG_CHART_PARSER
+  string hint_last_rule;
+#endif
+  EGrammarNode* cur = &(*g)[tr.GetLHS()];
+  cur->is_root = true;
+  for (int i = 0; i < tr.FLength(); ++i) {
+    WordID sym = tr.f()[i];
+#ifdef DEBUG_CHART_PARSER
+    hint_last_rule = TD::Convert(sym < 0 ? -sym : sym);
+    cur->hint += " <@@> (*" + hint_last_rule + ") " + TrimRule(tr.AsString());
+#endif
+    if (sym < 0)
+      cur = &cur->ntptr[sym];
+    else
+      cur = &cur->tptr[sym];
+  }
+#ifdef DEBUG_CHART_PARSER
+  cur->hint += " <@@> (" + hint_last_rule + "*) " + TrimRule(tr.AsString());
+#endif
+  cur->is_some_rule_complete = true;
+  cur->input_features = tr.GetFeatureValues();
+}
+
+CFG_WFSTComposer::~CFG_WFSTComposer() {
+  delete pimpl_;
+}
+
+CFG_WFSTComposer::CFG_WFSTComposer(const WFST& wfst) {
+  InitializeConstants();
+  pimpl_ = new CFG_WFSTComposerImpl(kUNIQUE_START, wfst.Initial(), wfst.Final());
+}
+
+bool CFG_WFSTComposer::Compose(const Hypergraph& src_forest, Hypergraph* trg_forest) {
+  // first, convert the src forest into an EGrammar
+  EGrammar g;
+  const int nedges = src_forest.edges_.size();
+  const int nnodes = src_forest.nodes_.size();
+  vector<int> cats(nnodes);
+  bool assign_cats = false;
+  for (int i = 0; i < nnodes; ++i)
+    if (assign_cats) {
+      cats[i] = TD::Convert("CAT_" + boost::lexical_cast<string>(i)) * -1;
+    } else {
+      cats[i] = src_forest.nodes_[i].cat_;
+    }
+  // construct the grammar
+  for (int i = 0; i < nedges; ++i) {
+    const Hypergraph::Edge& edge = src_forest.edges_[i];
+    const vector<WordID>& src = edge.rule_->f();
+    EGrammarNode* cur = &g[cats[edge.head_node_]];
+    cur->is_root = true;
+    int ntc = 0;
+    for (int j = 0; j < src.size(); ++j) {
+      WordID sym = src[j];
+      if (sym <= 0) {
+        sym = cats[edge.tail_nodes_[ntc]];
+        ++ntc;
+        cur = &cur->ntptr[sym];
+      } else {
+        cur = &cur->tptr[sym];
+      }
+    }
+    cur->is_some_rule_complete = true;
+    cur->input_features = edge.feature_values_;
+  }
+  EGrammarNode& goal_rule = g[kUNIQUE_START];
+  assert((goal_rule.ntptr.size() == 1 && goal_rule.tptr.size() == 0) ||
+         (goal_rule.ntptr.size() == 0 && goal_rule.tptr.size() == 1));
+
+  return pimpl_->Compose(g, trg_forest);
+}
+
+bool CFG_WFSTComposer::Compose(istream* in, Hypergraph* trg_forest) {
+  EGrammar g;
+  while(*in) {
+    string line;
+    getline(*in, line);
+    if (line.empty()) continue;
+    AddGrammarRule(line, &g);
+  }
+
+  return pimpl_->Compose(g, trg_forest);
+}
diff --git a/gi/pf/cfg_wfst_composer.h b/gi/pf/cfg_wfst_composer.h
new file mode 100644
index 00000000..cf47f459
--- /dev/null
+++ b/gi/pf/cfg_wfst_composer.h
@@ -0,0 +1,46 @@
+#ifndef _CFG_WFST_COMPOSER_H_
+#define _CFG_WFST_COMPOSER_H_
+
+#include <iostream>
+#include <vector>
+#include <utility>
+
+#include "trule.h"
+#include "wordid.h"
+
+class CFG_WFSTComposerImpl;
+class Hypergraph;
+
+struct WFSTNode {
+  virtual ~WFSTNode();
+  // returns the next states reachable by consuming srcindex (which identifies a word)
+  // paired with the output string generated by taking that transition.
+  virtual std::vector<std::pair<const WFSTNode*,TRulePtr> > ExtendInput(unsigned srcindex) const = 0;
+};
+
+struct WFST {
+  virtual ~WFST();
+  virtual const WFSTNode* Final() const = 0;
+  virtual const WFSTNode* Initial() const = 0;
+};
+
+class CFG_WFSTComposer {
+ public:
+  ~CFG_WFSTComposer();
+  explicit CFG_WFSTComposer(const WFST& wfst);
+  bool Compose(const Hypergraph& in_forest, Hypergraph* trg_forest);
+
+  // reads the grammar from a file. There must be a single top-level
+  // S -> X rule.  Anything else is possible. Format is:
+  // [S] ||| [SS,1]
+  // [SS] ||| [NP,1] [VP,2] ||| Feature1=0.2 Feature2=-2.3
+  // [SS] ||| [VP,1] [NP,2] ||| Feature1=0.8
+  // [NP] ||| [DET,1] [N,2] ||| Feature3=2
+  // ...
+  bool Compose(std::istream* grammar_file, Hypergraph* trg_forest);
+
+ private:
+  CFG_WFSTComposerImpl* pimpl_;
+};
+
+#endif
diff --git a/gi/pf/dpnaive.cc b/gi/pf/dpnaive.cc
new file mode 100644
index 00000000..582d1be7
--- /dev/null
+++ b/gi/pf/dpnaive.cc
@@ -0,0 +1,349 @@
+#include <iostream>
+#include <tr1/memory>
+#include <queue>
+
+#include <boost/multi_array.hpp>
+#include <boost/program_options.hpp>
+#include <boost/program_options/variables_map.hpp>
+
+#include "base_measures.h"
+#include "trule.h"
+#include "tdict.h"
+#include "filelib.h"
+#include "dict.h"
+#include "sampler.h"
+#include "ccrp_nt.h"
+
+using namespace std;
+using namespace std::tr1;
+namespace po = boost::program_options;
+
+static unsigned kMAX_SRC_PHRASE;
+static unsigned kMAX_TRG_PHRASE;
+struct FSTState;
+
+size_t hash_value(const TRule& r) {
+  size_t h = 2 - r.lhs_;
+  boost::hash_combine(h, boost::hash_value(r.e_));
+  boost::hash_combine(h, boost::hash_value(r.f_));
+  return h;
+}
+
+bool operator==(const TRule& a, const TRule& b) {
+  return (a.lhs_ == b.lhs_ && a.e_ == b.e_ && a.f_ == b.f_);
+}
+
+void InitCommandLine(int argc, char** argv, po::variables_map* conf) {
+  po::options_description opts("Configuration options");
+  opts.add_options()
+        ("samples,s",po::value<unsigned>()->default_value(1000),"Number of samples")
+        ("input,i",po::value<string>(),"Read parallel data from")
+        ("max_src_phrase",po::value<unsigned>()->default_value(4),"Maximum length of source language phrases")
+        ("max_trg_phrase",po::value<unsigned>()->default_value(4),"Maximum length of target language phrases")
+        ("model1,m",po::value<string>(),"Model 1 parameters (used in base distribution)")
+        ("model1_interpolation_weight",po::value<double>()->default_value(0.95),"Mixing proportion of model 1 with uniform target distribution")
+        ("random_seed,S",po::value<uint32_t>(), "Random seed");
+  po::options_description clo("Command line options");
+  clo.add_options()
+        ("config", po::value<string>(), "Configuration file")
+        ("help,h", "Print this help message and exit");
+  po::options_description dconfig_options, dcmdline_options;
+  dconfig_options.add(opts);
+  dcmdline_options.add(opts).add(clo);
+  
+  po::store(parse_command_line(argc, argv, dcmdline_options), *conf);
+  if (conf->count("config")) {
+    ifstream config((*conf)["config"].as<string>().c_str());
+    po::store(po::parse_config_file(config, dconfig_options), *conf);
+  }
+  po::notify(*conf);
+
+  if (conf->count("help") || (conf->count("input") == 0)) {
+    cerr << dcmdline_options << endl;
+    exit(1);
+  }
+}
+
+void ReadParallelCorpus(const string& filename,
+                vector<vector<WordID> >* f,
+                vector<vector<int> >* e,
+                set<int>* vocab_e,
+                set<int>* vocab_f) {
+  f->clear();
+  e->clear();
+  vocab_f->clear();
+  vocab_e->clear();
+  istream* in;
+  if (filename == "-")
+    in = &cin;
+  else
+    in = new ifstream(filename.c_str());
+  assert(*in);
+  string line;
+  const WordID kDIV = TD::Convert("|||");
+  vector<WordID> tmp;
+  while(*in) {
+    getline(*in, line);
+    if (line.empty() && !*in) break;
+    e->push_back(vector<int>());
+    f->push_back(vector<int>());
+    vector<int>& le = e->back();
+    vector<int>& lf = f->back();
+    tmp.clear();
+    TD::ConvertSentence(line, &tmp);
+    bool isf = true;
+    for (unsigned i = 0; i < tmp.size(); ++i) {
+      const int cur = tmp[i];
+      if (isf) {
+        if (kDIV == cur) { isf = false; } else {
+          lf.push_back(cur);
+          vocab_f->insert(cur);
+        }
+      } else {
+        assert(cur != kDIV);
+        le.push_back(cur);
+        vocab_e->insert(cur);
+      }
+    }
+    assert(isf == false);
+  }
+  if (in != &cin) delete in;
+}
+
+shared_ptr<MT19937> prng;
+
+template <typename Base>
+struct ModelAndData {
+  explicit ModelAndData(const Base& b, const vector<vector<int> >& ce, const vector<vector<int> >& cf, const set<int>& ve, const set<int>& vf) :
+     rng(&*prng),
+     p0(b),
+     baseprob(prob_t::One()),
+     corpuse(ce),
+     corpusf(cf),
+     vocabe(ve),
+     vocabf(vf),
+     rules(1,1),
+     mh_samples(),
+     mh_rejects(),
+     kX(-TD::Convert("X")),
+     derivations(corpuse.size()) {}
+
+  void ResampleHyperparameters() {
+    rules.resample_hyperparameters(&*prng);
+  }
+
+  void InstantiateRule(const pair<short,short>& from,
+                       const pair<short,short>& to,
+                       const vector<int>& sentf,
+                       const vector<int>& sente,
+                       TRule* rule) const {
+    rule->f_.clear();
+    rule->e_.clear();
+    rule->lhs_ = kX;
+    for (short i = from.first; i < to.first; ++i)
+      rule->f_.push_back(sentf[i]);
+    for (short i = from.second; i < to.second; ++i)
+      rule->e_.push_back(sente[i]);
+  }
+
+  void DecrementDerivation(const vector<pair<short,short> >& d, const vector<int>& sentf, const vector<int>& sente) {
+    if (d.size() < 2) return;
+    TRule x;
+    for (int i = 1; i < d.size(); ++i) {
+      InstantiateRule(d[i], d[i-1], sentf, sente, &x);
+      //cerr << "REMOVE: " << x.AsString() << endl;
+      if (rules.decrement(x)) {
+        baseprob /= p0(x);
+        //cerr << "  (REMOVED ONLY INSTANCE)\n";
+      }
+    }
+  }
+
+  void PrintDerivation(const vector<pair<short,short> >& d, const vector<int>& sentf, const vector<int>& sente) {
+    if (d.size() < 2) return;
+    TRule x;
+    for (int i = 1; i < d.size(); ++i) {
+      InstantiateRule(d[i], d[i-1], sentf, sente, &x);
+      cerr << i << '/' << (d.size() - 1) << ": " << x << endl;
+    }
+  }
+
+  void IncrementDerivation(const vector<pair<short,short> >& d, const vector<int>& sentf, const vector<int>& sente) {
+    if (d.size() < 2) return;
+    TRule x;
+    for (int i = 1; i < d.size(); ++i) {
+      InstantiateRule(d[i], d[i-1], sentf, sente, &x);
+      if (rules.increment(x)) {
+        baseprob *= p0(x);
+      }
+    }
+  }
+
+  prob_t Likelihood() const {
+    prob_t p;
+    p.logeq(rules.log_crp_prob());
+    return p * baseprob;
+  }
+
+  prob_t DerivationProposalProbability(const vector<pair<short,short> >& d, const vector<int>& sentf, const vector<int>& sente) const {
+    prob_t p = prob_t::One();
+    if (d.size() < 2) return p;
+    TRule x;
+    for (int i = 1; i < d.size(); ++i) {
+      InstantiateRule(d[i], d[i-1], sentf, sente, &x);
+      prob_t rp; rp.logeq(rules.logprob(x, log(p0(x))));
+      p *= rp;
+    }
+    return p;
+  }
+
+  void Sample();
+
+  MT19937* rng;
+  const Base& p0;
+  prob_t baseprob; // cached value of generating the table table labels from p0
+                   // this can't be used if we go to a hierarchical prior!
+  const vector<vector<int> >& corpuse, corpusf;
+  const set<int>& vocabe, vocabf;
+  CCRP_NoTable<TRule> rules;
+  unsigned mh_samples, mh_rejects;
+  const int kX;
+  vector<vector<pair<short, short> > > derivations;
+};
+
+template <typename Base>
+void ModelAndData<Base>::Sample() {
+  unsigned MAXK = 4;
+  unsigned MAXL = 4;
+  TRule x;
+  x.lhs_ = -TD::Convert("X");
+  for (int samples = 0; samples < 1000; ++samples) {
+    if (samples % 1 == 0 && samples > 0) {
+      //ResampleHyperparameters();
+      cerr << " [" << samples << " LLH=" << log(Likelihood()) << " MH=" << ((double)mh_rejects / mh_samples) << "]\n";
+      for (int i = 0; i < 10; ++i) {
+        cerr << "SENTENCE: " << TD::GetString(corpusf[i]) << " ||| " << TD::GetString(corpuse[i]) << endl;
+        PrintDerivation(derivations[i], corpusf[i], corpuse[i]);
+      }
+    }
+    cerr << '.' << flush;
+    for (int s = 0; s < corpuse.size(); ++s) {
+      const vector<int>& sentf = corpusf[s];
+      const vector<int>& sente = corpuse[s];
+//      cerr << "  CUSTOMERS: " << rules.num_customers() << endl;
+//      cerr << "SENTENCE: " << TD::GetString(sentf) << " ||| " << TD::GetString(sente) << endl;
+
+      vector<pair<short, short> >& deriv = derivations[s];
+      const prob_t p_cur = Likelihood();
+      DecrementDerivation(deriv, sentf, sente);
+
+      boost::multi_array<prob_t, 2> a(boost::extents[sentf.size() + 1][sente.size() + 1]);
+      boost::multi_array<prob_t, 4> trans(boost::extents[sentf.size() + 1][sente.size() + 1][MAXK][MAXL]);
+      a[0][0] = prob_t::One();
+      for (int i = 0; i < sentf.size(); ++i) {
+        for (int j = 0; j < sente.size(); ++j) {
+          const prob_t src_a = a[i][j];
+          x.f_.clear();
+          for (int k = 1; k <= MAXK; ++k) {
+            if (i + k > sentf.size()) break;
+            x.f_.push_back(sentf[i + k - 1]);
+            x.e_.clear();
+            for (int l = 1; l <= MAXL; ++l) {
+              if (j + l > sente.size()) break;
+              x.e_.push_back(sente[j + l - 1]);
+              trans[i][j][k - 1][l - 1].logeq(rules.logprob(x, log(p0(x))));
+              a[i + k][j + l] += src_a * trans[i][j][k - 1][l - 1];
+            }
+          }
+        }
+      }
+//      cerr << "Inside: " << log(a[sentf.size()][sente.size()]) << endl;
+      const prob_t q_cur = DerivationProposalProbability(deriv, sentf, sente);
+
+      vector<pair<short,short> > newderiv;
+      int cur_i = sentf.size();
+      int cur_j = sente.size();
+      while(cur_i > 0 && cur_j > 0) {
+        newderiv.push_back(pair<short,short>(cur_i, cur_j));
+//        cerr << "NODE: (" << cur_i << "," << cur_j << ")\n";
+        SampleSet<prob_t> ss;
+        vector<pair<short,short> > nexts;
+        for (int k = 1; k <= MAXK; ++k) {
+          const int hyp_i = cur_i - k;
+          if (hyp_i < 0) break;
+          for (int l = 1; l <= MAXL; ++l) {
+            const int hyp_j = cur_j - l;
+            if (hyp_j < 0) break;
+            const prob_t& inside = a[hyp_i][hyp_j];
+            if (inside == prob_t::Zero()) continue;
+            const prob_t& transp = trans[hyp_i][hyp_j][k - 1][l - 1];
+            if (transp == prob_t::Zero()) continue;
+            const prob_t p = inside * transp;
+            ss.add(p);
+            nexts.push_back(pair<short,short>(hyp_i, hyp_j));
+//            cerr << "    (" << hyp_i << "," << hyp_j << ")  <--- " << log(p) << endl;
+          }
+        }
+//        cerr << "  sample set has " << nexts.size() << " elements.\n";
+        const int selected = rng->SelectSample(ss);
+        cur_i = nexts[selected].first;
+        cur_j = nexts[selected].second;
+      }
+      newderiv.push_back(pair<short,short>(0,0));
+      const prob_t q_new = DerivationProposalProbability(newderiv, sentf, sente);
+      IncrementDerivation(newderiv, sentf, sente);
+//      cerr << "SANITY: " << q_new << "  " <<log(DerivationProposalProbability(newderiv, sentf, sente)) << endl;
+      if (deriv.empty()) { deriv = newderiv; continue; }
+      ++mh_samples;
+
+      if (deriv != newderiv) {
+        const prob_t p_new = Likelihood();
+//        cerr << "p_cur=" << log(p_cur) << "\t p_new=" << log(p_new) << endl;
+//        cerr << "q_cur=" << log(q_cur) << "\t q_new=" << log(q_new) << endl;
+        if (!rng->AcceptMetropolisHastings(p_new, p_cur, q_new, q_cur)) {
+          ++mh_rejects;
+          DecrementDerivation(newderiv, sentf, sente);
+          IncrementDerivation(deriv, sentf, sente);
+        } else {
+//          cerr << "  ACCEPT\n";
+          deriv = newderiv;
+        }
+      }
+    }
+  }
+}
+
+int main(int argc, char** argv) {
+  po::variables_map conf;
+  InitCommandLine(argc, argv, &conf);
+  kMAX_TRG_PHRASE = conf["max_trg_phrase"].as<unsigned>();
+  kMAX_SRC_PHRASE = conf["max_src_phrase"].as<unsigned>();
+
+  if (!conf.count("model1")) {
+    cerr << argv[0] << "Please use --model1 to specify model 1 parameters\n";
+    return 1;
+  }
+  if (conf.count("random_seed"))
+    prng.reset(new MT19937(conf["random_seed"].as<uint32_t>()));
+  else
+    prng.reset(new MT19937);
+//  MT19937& rng = *prng;
+
+  vector<vector<int> > corpuse, corpusf;
+  set<int> vocabe, vocabf;
+  ReadParallelCorpus(conf["input"].as<string>(), &corpusf, &corpuse, &vocabf, &vocabe);
+  cerr << "f-Corpus size: " << corpusf.size() << " sentences\n";
+  cerr << "f-Vocabulary size: " << vocabf.size() << " types\n";
+  cerr << "f-Corpus size: " << corpuse.size() << " sentences\n";
+  cerr << "f-Vocabulary size: " << vocabe.size() << " types\n";
+  assert(corpusf.size() == corpuse.size());
+
+  Model1 m1(conf["model1"].as<string>());
+  PhraseJointBase lp0(m1, conf["model1_interpolation_weight"].as<double>(), vocabe.size(), vocabf.size());
+
+  ModelAndData<PhraseJointBase> posterior(lp0, corpuse, corpusf, vocabe, vocabf);
+  posterior.Sample();
+
+  return 0;
+}
+
diff --git a/gi/pf/itg.cc b/gi/pf/itg.cc
new file mode 100644
index 00000000..2c2a86f9
--- /dev/null
+++ b/gi/pf/itg.cc
@@ -0,0 +1,224 @@
+#include <iostream>
+#include <tr1/memory>
+#include <queue>
+
+#include <boost/functional.hpp>
+#include <boost/program_options.hpp>
+#include <boost/program_options/variables_map.hpp>
+
+#include "viterbi.h"
+#include "hg.h"
+#include "trule.h"
+#include "tdict.h"
+#include "filelib.h"
+#include "dict.h"
+#include "sampler.h"
+#include "ccrp_nt.h"
+#include "ccrp_onetable.h"
+
+using namespace std;
+using namespace tr1;
+namespace po = boost::program_options;
+
+ostream& operator<<(ostream& os, const vector<WordID>& p) {
+  os << '[';
+  for (int i = 0; i < p.size(); ++i)
+    os << (i==0 ? "" : " ") << TD::Convert(p[i]);
+  return os << ']';
+}
+
+size_t hash_value(const TRule& r) {
+  size_t h = boost::hash_value(r.e_);
+  boost::hash_combine(h, -r.lhs_);
+  boost::hash_combine(h, boost::hash_value(r.f_));
+  return h;
+}
+
+bool operator==(const TRule& a, const TRule& b) {
+  return (a.lhs_ == b.lhs_ && a.e_ == b.e_ && a.f_ == b.f_);
+}
+
+double log_poisson(unsigned x, const double& lambda) {
+  assert(lambda > 0.0);
+  return log(lambda) * x - lgamma(x + 1) - lambda;
+}
+
+struct Model1 {
+  explicit Model1(const string& fname) :
+      kNULL(TD::Convert("<eps>")),
+      kZERO() {
+    LoadModel1(fname);
+  }
+
+  void LoadModel1(const string& fname) {
+    cerr << "Loading Model 1 parameters from " << fname << " ..." << endl;
+    ReadFile rf(fname);
+    istream& in = *rf.stream();
+    string line;
+    unsigned lc = 0;
+    while(getline(in, line)) {
+      ++lc;
+      int cur = 0;
+      int start = 0;
+      while(cur < line.size() && line[cur] != ' ') { ++cur; }
+      assert(cur != line.size());
+      line[cur] = 0;
+      const WordID src = TD::Convert(&line[0]);
+      ++cur;
+      start = cur;
+      while(cur < line.size() && line[cur] != ' ') { ++cur; }
+      assert(cur != line.size());
+      line[cur] = 0;
+      WordID trg = TD::Convert(&line[start]);
+      const double logprob = strtod(&line[cur + 1], NULL);
+      if (src >= ttable.size()) ttable.resize(src + 1);
+      ttable[src][trg].logeq(logprob);
+    }
+    cerr << "  read " << lc << " parameters.\n";
+  }
+
+  // returns prob 0 if src or trg is not found!
+  const prob_t& operator()(WordID src, WordID trg) const {
+    if (src == 0) src = kNULL;
+    if (src < ttable.size()) {
+      const map<WordID, prob_t>& cpd = ttable[src];
+      const map<WordID, prob_t>::const_iterator it = cpd.find(trg);
+      if (it != cpd.end())
+        return it->second;
+    }
+    return kZERO;
+  }
+
+  const WordID kNULL;
+  const prob_t kZERO;
+  vector<map<WordID, prob_t> > ttable;
+};
+
+void InitCommandLine(int argc, char** argv, po::variables_map* conf) {
+  po::options_description opts("Configuration options");
+  opts.add_options()
+        ("samples,s",po::value<unsigned>()->default_value(1000),"Number of samples")
+        ("particles,p",po::value<unsigned>()->default_value(25),"Number of particles")
+        ("input,i",po::value<string>(),"Read parallel data from")
+        ("max_src_phrase",po::value<unsigned>()->default_value(7),"Maximum length of source language phrases")
+        ("max_trg_phrase",po::value<unsigned>()->default_value(7),"Maximum length of target language phrases")
+        ("model1,m",po::value<string>(),"Model 1 parameters (used in base distribution)")
+        ("inverse_model1,M",po::value<string>(),"Inverse Model 1 parameters (used in backward estimate)")
+        ("model1_interpolation_weight",po::value<double>()->default_value(0.95),"Mixing proportion of model 1 with uniform target distribution")
+        ("random_seed,S",po::value<uint32_t>(), "Random seed");
+  po::options_description clo("Command line options");
+  clo.add_options()
+        ("config", po::value<string>(), "Configuration file")
+        ("help,h", "Print this help message and exit");
+  po::options_description dconfig_options, dcmdline_options;
+  dconfig_options.add(opts);
+  dcmdline_options.add(opts).add(clo);
+  
+  po::store(parse_command_line(argc, argv, dcmdline_options), *conf);
+  if (conf->count("config")) {
+    ifstream config((*conf)["config"].as<string>().c_str());
+    po::store(po::parse_config_file(config, dconfig_options), *conf);
+  }
+  po::notify(*conf);
+
+  if (conf->count("help") || (conf->count("input") == 0)) {
+    cerr << dcmdline_options << endl;
+    exit(1);
+  }
+}
+
+void ReadParallelCorpus(const string& filename,
+                vector<vector<WordID> >* f,
+                vector<vector<WordID> >* e,
+                set<WordID>* vocab_f,
+                set<WordID>* vocab_e) {
+  f->clear();
+  e->clear();
+  vocab_f->clear();
+  vocab_e->clear();
+  istream* in;
+  if (filename == "-")
+    in = &cin;
+  else
+    in = new ifstream(filename.c_str());
+  assert(*in);
+  string line;
+  const WordID kDIV = TD::Convert("|||");
+  vector<WordID> tmp;
+  while(*in) {
+    getline(*in, line);
+    if (line.empty() && !*in) break;
+    e->push_back(vector<int>());
+    f->push_back(vector<int>());
+    vector<int>& le = e->back();
+    vector<int>& lf = f->back();
+    tmp.clear();
+    TD::ConvertSentence(line, &tmp);
+    bool isf = true;
+    for (unsigned i = 0; i < tmp.size(); ++i) {
+      const int cur = tmp[i];
+      if (isf) {
+        if (kDIV == cur) { isf = false; } else {
+          lf.push_back(cur);
+          vocab_f->insert(cur);
+        }
+      } else {
+        assert(cur != kDIV);
+        le.push_back(cur);
+        vocab_e->insert(cur);
+      }
+    }
+    assert(isf == false);
+  }
+  if (in != &cin) delete in;
+}
+
+int main(int argc, char** argv) {
+  po::variables_map conf;
+  InitCommandLine(argc, argv, &conf);
+  const size_t kMAX_TRG_PHRASE = conf["max_trg_phrase"].as<unsigned>();
+  const size_t kMAX_SRC_PHRASE = conf["max_src_phrase"].as<unsigned>();
+  const unsigned particles = conf["particles"].as<unsigned>();
+  const unsigned samples = conf["samples"].as<unsigned>();
+
+  if (!conf.count("model1")) {
+    cerr << argv[0] << "Please use --model1 to specify model 1 parameters\n";
+    return 1;
+  }
+  shared_ptr<MT19937> prng;
+  if (conf.count("random_seed"))
+    prng.reset(new MT19937(conf["random_seed"].as<uint32_t>()));
+  else
+    prng.reset(new MT19937);
+  MT19937& rng = *prng;
+
+  vector<vector<WordID> > corpuse, corpusf;
+  set<WordID> vocabe, vocabf;
+  cerr << "Reading corpus...\n";
+  ReadParallelCorpus(conf["input"].as<string>(), &corpusf, &corpuse, &vocabf, &vocabe);
+  cerr << "F-corpus size: " << corpusf.size() << " sentences\t (" << vocabf.size() << " word types)\n";
+  cerr << "E-corpus size: " << corpuse.size() << " sentences\t (" << vocabe.size() << " word types)\n";
+  assert(corpusf.size() == corpuse.size());
+
+  const int kLHS = -TD::Convert("X");
+  Model1 m1(conf["model1"].as<string>());
+  Model1 invm1(conf["inverse_model1"].as<string>());
+  for (int si = 0; si < conf["samples"].as<unsigned>(); ++si) {
+    cerr << '.' << flush;
+    for (int ci = 0; ci < corpusf.size(); ++ci) {
+      const vector<WordID>& src = corpusf[ci];
+      const vector<WordID>& trg = corpuse[ci];
+      for (int i = 0; i < src.size(); ++i) {
+        for (int j = 0; j < trg.size(); ++j) {
+          const int eff_max_src = min(src.size() - i, kMAX_SRC_PHRASE);
+          for (int k = 0; k < eff_max_src; ++k) {
+            const int eff_max_trg = (k == 0 ? 1 : min(trg.size() - j, kMAX_TRG_PHRASE));
+            for (int l = 0; l < eff_max_trg; ++l) {
+            }
+          }
+        }
+      }
+    }
+  }
+}
+
diff --git a/gi/pf/pfbrat.cc b/gi/pf/pfbrat.cc
new file mode 100644
index 00000000..4c6ba3ef
--- /dev/null
+++ b/gi/pf/pfbrat.cc
@@ -0,0 +1,554 @@
+#include <iostream>
+#include <tr1/memory>
+#include <queue>
+
+#include <boost/functional.hpp>
+#include <boost/multi_array.hpp>
+#include <boost/program_options.hpp>
+#include <boost/program_options/variables_map.hpp>
+
+#include "viterbi.h"
+#include "hg.h"
+#include "trule.h"
+#include "tdict.h"
+#include "filelib.h"
+#include "dict.h"
+#include "sampler.h"
+#include "ccrp_nt.h"
+#include "cfg_wfst_composer.h"
+
+using namespace std;
+using namespace tr1;
+namespace po = boost::program_options;
+
+static unsigned kMAX_SRC_PHRASE;
+static unsigned kMAX_TRG_PHRASE;
+struct FSTState;
+
+size_t hash_value(const TRule& r) {
+  size_t h = 2 - r.lhs_;
+  boost::hash_combine(h, boost::hash_value(r.e_));
+  boost::hash_combine(h, boost::hash_value(r.f_));
+  return h;
+}
+
+bool operator==(const TRule& a, const TRule& b) {
+  return (a.lhs_ == b.lhs_ && a.e_ == b.e_ && a.f_ == b.f_);
+}
+
+double log_poisson(unsigned x, const double& lambda) {
+  assert(lambda > 0.0);
+  return log(lambda) * x - lgamma(x + 1) - lambda;
+}
+
+struct ConditionalBase {
+  explicit ConditionalBase(const double m1mixture, const unsigned vocab_e_size, const string& model1fname) :
+      kM1MIXTURE(m1mixture),
+      kUNIFORM_MIXTURE(1.0 - m1mixture),
+      kUNIFORM_TARGET(1.0 / vocab_e_size),
+      kNULL(TD::Convert("<eps>")) {
+    assert(m1mixture >= 0.0 && m1mixture <= 1.0);
+    assert(vocab_e_size > 0);
+    LoadModel1(model1fname);
+  }
+
+  void LoadModel1(const string& fname) {
+    cerr << "Loading Model 1 parameters from " << fname << " ..." << endl;
+    ReadFile rf(fname);
+    istream& in = *rf.stream();
+    string line;
+    unsigned lc = 0;
+    while(getline(in, line)) {
+      ++lc;
+      int cur = 0;
+      int start = 0;
+      while(cur < line.size() && line[cur] != ' ') { ++cur; }
+      assert(cur != line.size());
+      line[cur] = 0;
+      const WordID src = TD::Convert(&line[0]);
+      ++cur;
+      start = cur;
+      while(cur < line.size() && line[cur] != ' ') { ++cur; }
+      assert(cur != line.size());
+      line[cur] = 0;
+      WordID trg = TD::Convert(&line[start]);
+      const double logprob = strtod(&line[cur + 1], NULL);
+      if (src >= ttable.size()) ttable.resize(src + 1);
+      ttable[src][trg].logeq(logprob);
+    }
+    cerr << "  read " << lc << " parameters.\n";
+  }
+
+  // return logp0 of rule.e_ | rule.f_
+  prob_t operator()(const TRule& rule) const {
+    const int flen = rule.f_.size();
+    const int elen = rule.e_.size();
+    prob_t uniform_src_alignment; uniform_src_alignment.logeq(-log(flen + 1));
+    prob_t p;
+    p.logeq(log_poisson(elen, flen + 0.01));       // elen | flen          ~Pois(flen + 0.01)
+    for (int i = 0; i < elen; ++i) {               // for each position i in e-RHS
+      const WordID trg = rule.e_[i];
+      prob_t tp = prob_t::Zero();
+      for (int j = -1; j < flen; ++j) {
+        const WordID src = j < 0 ? kNULL : rule.f_[j];
+        const map<WordID, prob_t>::const_iterator it = ttable[src].find(trg);
+        if (it != ttable[src].end()) {
+          tp += kM1MIXTURE * it->second;
+        }
+        tp += kUNIFORM_MIXTURE * kUNIFORM_TARGET;
+      }
+      tp *= uniform_src_alignment;                 //     draw a_i         ~uniform
+      p *= tp;                                     //     draw e_i         ~Model1(f_a_i) / uniform
+    }
+    return p;
+  }
+
+  const prob_t kM1MIXTURE;  // Model 1 mixture component
+  const prob_t kUNIFORM_MIXTURE; // uniform mixture component
+  const prob_t kUNIFORM_TARGET;
+  const WordID kNULL;
+  vector<map<WordID, prob_t> > ttable;
+};
+
+void InitCommandLine(int argc, char** argv, po::variables_map* conf) {
+  po::options_description opts("Configuration options");
+  opts.add_options()
+        ("samples,s",po::value<unsigned>()->default_value(1000),"Number of samples")
+        ("input,i",po::value<string>(),"Read parallel data from")
+        ("max_src_phrase",po::value<unsigned>()->default_value(3),"Maximum length of source language phrases")
+        ("max_trg_phrase",po::value<unsigned>()->default_value(3),"Maximum length of target language phrases")
+        ("model1,m",po::value<string>(),"Model 1 parameters (used in base distribution)")
+        ("model1_interpolation_weight",po::value<double>()->default_value(0.95),"Mixing proportion of model 1 with uniform target distribution")
+        ("random_seed,S",po::value<uint32_t>(), "Random seed");
+  po::options_description clo("Command line options");
+  clo.add_options()
+        ("config", po::value<string>(), "Configuration file")
+        ("help,h", "Print this help message and exit");
+  po::options_description dconfig_options, dcmdline_options;
+  dconfig_options.add(opts);
+  dcmdline_options.add(opts).add(clo);
+  
+  po::store(parse_command_line(argc, argv, dcmdline_options), *conf);
+  if (conf->count("config")) {
+    ifstream config((*conf)["config"].as<string>().c_str());
+    po::store(po::parse_config_file(config, dconfig_options), *conf);
+  }
+  po::notify(*conf);
+
+  if (conf->count("help") || (conf->count("input") == 0)) {
+    cerr << dcmdline_options << endl;
+    exit(1);
+  }
+}
+
+void ReadParallelCorpus(const string& filename,
+                vector<vector<WordID> >* f,
+                vector<vector<int> >* e,
+                set<int>* vocab_f,
+                set<int>* vocab_e) {
+  f->clear();
+  e->clear();
+  vocab_f->clear();
+  vocab_e->clear();
+  istream* in;
+  if (filename == "-")
+    in = &cin;
+  else
+    in = new ifstream(filename.c_str());
+  assert(*in);
+  string line;
+  const WordID kDIV = TD::Convert("|||");
+  vector<WordID> tmp;
+  while(*in) {
+    getline(*in, line);
+    if (line.empty() && !*in) break;
+    e->push_back(vector<int>());
+    f->push_back(vector<int>());
+    vector<int>& le = e->back();
+    vector<int>& lf = f->back();
+    tmp.clear();
+    TD::ConvertSentence(line, &tmp);
+    bool isf = true;
+    for (unsigned i = 0; i < tmp.size(); ++i) {
+      const int cur = tmp[i];
+      if (isf) {
+        if (kDIV == cur) { isf = false; } else {
+          lf.push_back(cur);
+          vocab_f->insert(cur);
+        }
+      } else {
+        assert(cur != kDIV);
+        le.push_back(cur);
+        vocab_e->insert(cur);
+      }
+    }
+    assert(isf == false);
+  }
+  if (in != &cin) delete in;
+}
+
+struct UniphraseLM {
+  UniphraseLM(const vector<vector<int> >& corpus,
+              const set<int>& vocab,
+              const po::variables_map& conf) :
+    phrases_(1,1),
+    gen_(1,1),
+    corpus_(corpus),
+    uniform_word_(1.0 / vocab.size()),
+    gen_p0_(0.5),
+    p_end_(0.5),
+    use_poisson_(conf.count("poisson_length") > 0) {}
+
+  void ResampleHyperparameters(MT19937* rng) {
+    phrases_.resample_hyperparameters(rng);
+    gen_.resample_hyperparameters(rng);
+    cerr << " " << phrases_.concentration();
+  }
+
+  CCRP_NoTable<vector<int> > phrases_;
+  CCRP_NoTable<bool> gen_;
+  vector<vector<bool> > z_;   // z_[i] is there a phrase boundary after the ith word
+  const vector<vector<int> >& corpus_;
+  const double uniform_word_;
+  const double gen_p0_;
+  const double p_end_; // in base length distribution, p of the end of a phrase
+  const bool use_poisson_;
+};
+
+struct Reachability {
+  boost::multi_array<bool, 4> edges;  // edges[src_covered][trg_covered][x][trg_delta] is this edge worth exploring?
+  boost::multi_array<short, 2> max_src_delta; // msd[src_covered][trg_covered] -- the largest src delta that's valid
+
+  Reachability(int srclen, int trglen, int src_max_phrase_len, int trg_max_phrase_len) :
+      edges(boost::extents[srclen][trglen][src_max_phrase_len+1][trg_max_phrase_len+1]),
+      max_src_delta(boost::extents[srclen][trglen]) {
+    ComputeReachability(srclen, trglen, src_max_phrase_len, trg_max_phrase_len);
+  }
+
+ private:
+  struct SState {
+    SState() : prev_src_covered(), prev_trg_covered() {}
+    SState(int i, int j) : prev_src_covered(i), prev_trg_covered(j) {}
+    int prev_src_covered;
+    int prev_trg_covered;
+  };
+
+  struct NState {
+    NState() : next_src_covered(), next_trg_covered() {}
+    NState(int i, int j) : next_src_covered(i), next_trg_covered(j) {}
+    int next_src_covered;
+    int next_trg_covered;
+  };
+
+  void ComputeReachability(int srclen, int trglen, int src_max_phrase_len, int trg_max_phrase_len) {
+    typedef boost::multi_array<vector<SState>, 2> array_type;
+    array_type a(boost::extents[srclen + 1][trglen + 1]);
+    a[0][0].push_back(SState());
+    for (int i = 0; i < srclen; ++i) {
+      for (int j = 0; j < trglen; ++j) {
+        if (a[i][j].size() == 0) continue;
+        const SState prev(i,j);
+        for (int k = 1; k <= src_max_phrase_len; ++k) {
+          if ((i + k) > srclen) continue;
+          for (int l = 1; l <= trg_max_phrase_len; ++l) {
+            if ((j + l) > trglen) continue;
+            a[i + k][j + l].push_back(prev);
+          }
+        }
+      }
+    }
+    a[0][0].clear();
+    cerr << "Final cell contains " << a[srclen][trglen].size() << " back pointers\n";
+    assert(a[srclen][trglen].size() > 0);
+
+    typedef boost::multi_array<bool, 2> rarray_type;
+    rarray_type r(boost::extents[srclen + 1][trglen + 1]);
+//    typedef boost::multi_array<vector<NState>, 2> narray_type;
+//    narray_type b(boost::extents[srclen + 1][trglen + 1]);
+    r[srclen][trglen] = true;
+    for (int i = srclen; i >= 0; --i) {
+      for (int j = trglen; j >= 0; --j) {
+        vector<SState>& prevs = a[i][j];
+        if (!r[i][j]) { prevs.clear(); }
+//        const NState nstate(i,j);
+        for (int k = 0; k < prevs.size(); ++k) {
+          r[prevs[k].prev_src_covered][prevs[k].prev_trg_covered] = true;
+          int src_delta = i - prevs[k].prev_src_covered;
+          edges[prevs[k].prev_src_covered][prevs[k].prev_trg_covered][src_delta][j - prevs[k].prev_trg_covered] = true;
+          short &msd = max_src_delta[prevs[k].prev_src_covered][prevs[k].prev_trg_covered];
+          if (src_delta > msd) msd = src_delta;
+//          b[prevs[k].prev_src_covered][prevs[k].prev_trg_covered].push_back(nstate);
+        }
+      }
+    }
+    assert(!edges[0][0][1][0]);
+    assert(!edges[0][0][0][1]);
+    assert(!edges[0][0][0][0]);
+    cerr << "  MAX SRC DELTA[0][0] = " << max_src_delta[0][0] << endl;
+    assert(max_src_delta[0][0] > 0);
+    //cerr << "First cell contains " << b[0][0].size() << " forward pointers\n";
+    //for (int i = 0; i < b[0][0].size(); ++i) {
+    //  cerr << "  -> (" << b[0][0][i].next_src_covered << "," << b[0][0][i].next_trg_covered << ")\n";
+    //}
+  }
+};
+
+ostream& operator<<(ostream& os, const FSTState& q);
+struct FSTState {
+  explicit FSTState(int src_size) :
+      trg_covered_(),
+      src_covered_(),
+      src_coverage_(src_size) {}
+
+  FSTState(short trg_covered, short src_covered, const vector<bool>& src_coverage, const vector<short>& src_prefix) :
+      trg_covered_(trg_covered),
+      src_covered_(src_covered),
+      src_coverage_(src_coverage),
+      src_prefix_(src_prefix) {
+    if (src_coverage_.size() == src_covered) {
+      assert(src_prefix.size() == 0);
+    }
+  }
+
+  // if we extend by the word at src_position, what are
+  // the next states that are reachable and lie on a valid
+  // path to the final state?
+  vector<FSTState> Extensions(int src_position, int src_len, int trg_len, const Reachability& r) const {
+    assert(src_position < src_coverage_.size());
+    if (src_coverage_[src_position]) {
+      cerr << "Trying to extend " << *this << " with position " << src_position << endl;
+      abort();
+    }
+    vector<bool> ncvg = src_coverage_;
+    ncvg[src_position] = true;
+
+    vector<FSTState> res;
+    const int trg_remaining = trg_len - trg_covered_;
+    if (trg_remaining <= 0) {
+      cerr << "Target appears to have been covered: " << *this << " (trg_len=" << trg_len << ",trg_covered=" << trg_covered_ << ")" << endl;
+      abort();
+    }
+    const int src_remaining = src_len - src_covered_;
+    if (src_remaining <= 0) {
+      cerr << "Source appears to have been covered: " << *this << endl;
+      abort();
+    }
+
+    for (int tc = 1; tc <= kMAX_TRG_PHRASE; ++tc) {
+      if (r.edges[src_covered_][trg_covered_][src_prefix_.size() + 1][tc]) {
+        int nc = src_prefix_.size() + 1 + src_covered_;
+        res.push_back(FSTState(trg_covered_ + tc, nc, ncvg, vector<short>()));
+      }
+    }
+
+    if ((src_prefix_.size() + 1) < r.max_src_delta[src_covered_][trg_covered_]) {
+      vector<short> nsp = src_prefix_;
+      nsp.push_back(src_position);
+      res.push_back(FSTState(trg_covered_, src_covered_, ncvg, nsp));
+    }
+
+    if (res.size() == 0) {
+      cerr << *this << " can't be extended!\n";
+      abort();
+    }
+    return res;
+  }
+
+  short trg_covered_, src_covered_;
+  vector<bool> src_coverage_;
+  vector<short> src_prefix_;
+};
+bool operator<(const FSTState& q, const FSTState& r) {
+  if (q.trg_covered_ != r.trg_covered_) return q.trg_covered_ < r.trg_covered_;
+  if (q.src_covered_!= r.src_covered_) return q.src_covered_ < r.src_covered_;
+  if (q.src_coverage_ != r.src_coverage_) return q.src_coverage_ < r.src_coverage_;
+  return q.src_prefix_ < r.src_prefix_;
+}
+
+ostream& operator<<(ostream& os, const FSTState& q) {
+  os << "[" << q.trg_covered_ << " : ";
+  for (int i = 0; i < q.src_coverage_.size(); ++i)
+    os << q.src_coverage_[i];
+  os << " : <";
+  for (int i = 0; i < q.src_prefix_.size(); ++i) {
+    if (i != 0) os << ' ';
+    os << q.src_prefix_[i];
+  }
+  return os << ">]";
+}
+
+struct MyModel {
+  MyModel(ConditionalBase& rcp0) : rp0(rcp0) {}
+  typedef unordered_map<vector<WordID>, CCRP_NoTable<TRule>, boost::hash<vector<WordID> > > SrcToRuleCRPMap;
+
+  void DecrementRule(const TRule& rule) {
+    SrcToRuleCRPMap::iterator it = rules.find(rule.f_);
+    assert(it != rules.end());
+    it->second.decrement(rule);
+    if (it->second.num_customers() == 0) rules.erase(it);
+  }
+
+  void IncrementRule(const TRule& rule) {
+    SrcToRuleCRPMap::iterator it = rules.find(rule.f_);
+    if (it == rules.end()) {
+      CCRP_NoTable<TRule> crp(1,1);
+      it = rules.insert(make_pair(rule.f_, crp)).first;
+    }
+    it->second.increment(rule);
+  }
+
+  // conditioned on rule.f_
+  prob_t RuleConditionalProbability(const TRule& rule) const {
+    const prob_t base = rp0(rule);
+    SrcToRuleCRPMap::const_iterator it = rules.find(rule.f_);
+    if (it == rules.end()) {
+      return base;
+    } else {
+      const double lp = it->second.logprob(rule, log(base));
+      prob_t q; q.logeq(lp);
+      return q;
+    }
+  }
+
+  const ConditionalBase& rp0;
+  SrcToRuleCRPMap rules;
+};
+
+struct MyFST : public WFST {
+  MyFST(const vector<WordID>& ssrc, const vector<WordID>& strg, MyModel* m) :
+      src(ssrc), trg(strg),
+      r(src.size(),trg.size(),kMAX_SRC_PHRASE, kMAX_TRG_PHRASE),
+      model(m) {
+    FSTState in(src.size());
+    cerr << " INIT: " << in << endl;
+    init = GetNode(in);
+    for (int i = 0; i < in.src_coverage_.size(); ++i) in.src_coverage_[i] = true;
+    in.src_covered_ = src.size();
+    in.trg_covered_ = trg.size();
+    cerr << "FINAL: " << in << endl;
+    final = GetNode(in);
+  }
+  virtual const WFSTNode* Final() const;
+  virtual const WFSTNode* Initial() const;
+
+  const WFSTNode* GetNode(const FSTState& q);
+  map<FSTState, boost::shared_ptr<WFSTNode> > m;
+  const vector<WordID>& src;
+  const vector<WordID>& trg;
+  Reachability r;
+  const WFSTNode* init;
+  const WFSTNode* final;
+  MyModel* model;
+};
+
+struct MyNode : public WFSTNode {
+  MyNode(const FSTState& q, MyFST* fst) : state(q), container(fst) {}
+  virtual vector<pair<const WFSTNode*, TRulePtr> > ExtendInput(unsigned srcindex) const;
+  const FSTState state;
+  mutable MyFST* container;
+};
+
+vector<pair<const WFSTNode*, TRulePtr> > MyNode::ExtendInput(unsigned srcindex) const {
+  cerr << "EXTEND " << state << " with " << srcindex << endl;
+  vector<FSTState> ext = state.Extensions(srcindex, container->src.size(), container->trg.size(), container->r);
+  vector<pair<const WFSTNode*,TRulePtr> > res(ext.size());
+  for (unsigned i = 0; i < ext.size(); ++i) {
+    res[i].first = container->GetNode(ext[i]);
+    if (ext[i].src_prefix_.size() == 0) {
+      const unsigned trg_from = state.trg_covered_;
+      const unsigned trg_to = ext[i].trg_covered_;
+      const unsigned prev_prfx_size = state.src_prefix_.size();
+      res[i].second.reset(new TRule);
+      res[i].second->lhs_ = -TD::Convert("X");
+      vector<WordID>& src = res[i].second->f_;
+      vector<WordID>& trg = res[i].second->e_;
+      src.resize(prev_prfx_size + 1);
+      for (unsigned j = 0; j < prev_prfx_size; ++j)
+        src[j] = container->src[state.src_prefix_[j]];
+      src[prev_prfx_size] = container->src[srcindex];
+      for (unsigned j = trg_from; j < trg_to; ++j)
+        trg.push_back(container->trg[j]);
+      res[i].second->scores_.set_value(FD::Convert("Proposal"), log(container->model->RuleConditionalProbability(*res[i].second)));
+    }
+  }
+  return res;
+}
+
+const WFSTNode* MyFST::GetNode(const FSTState& q) {
+  boost::shared_ptr<WFSTNode>& res = m[q];
+  if (!res) {
+    res.reset(new MyNode(q, this));
+  }
+  return &*res;
+}
+
+const WFSTNode* MyFST::Final() const {
+  return final;
+}
+
+const WFSTNode* MyFST::Initial() const {
+  return init;
+}
+
+int main(int argc, char** argv) {
+  po::variables_map conf;
+  InitCommandLine(argc, argv, &conf);
+  kMAX_TRG_PHRASE = conf["max_trg_phrase"].as<unsigned>();
+  kMAX_SRC_PHRASE = conf["max_src_phrase"].as<unsigned>();
+
+  if (!conf.count("model1")) {
+    cerr << argv[0] << "Please use --model1 to specify model 1 parameters\n";
+    return 1;
+  }
+  shared_ptr<MT19937> prng;
+  if (conf.count("random_seed"))
+    prng.reset(new MT19937(conf["random_seed"].as<uint32_t>()));
+  else
+    prng.reset(new MT19937);
+  MT19937& rng = *prng;
+
+  vector<vector<int> > corpuse, corpusf;
+  set<int> vocabe, vocabf;
+  ReadParallelCorpus(conf["input"].as<string>(), &corpusf, &corpuse, &vocabf, &vocabe);
+  cerr << "f-Corpus size: " << corpusf.size() << " sentences\n";
+  cerr << "f-Vocabulary size: " << vocabf.size() << " types\n";
+  cerr << "f-Corpus size: " << corpuse.size() << " sentences\n";
+  cerr << "f-Vocabulary size: " << vocabe.size() << " types\n";
+  assert(corpusf.size() == corpuse.size());
+
+  ConditionalBase lp0(conf["model1_interpolation_weight"].as<double>(),
+                      vocabe.size(),
+                      conf["model1"].as<string>());
+  MyModel m(lp0);
+
+  TRule x("[X] ||| kAnwntR myN ||| at the convent ||| 0");
+  m.IncrementRule(x);
+  TRule y("[X] ||| nY dyN ||| gave ||| 0");
+  m.IncrementRule(y);
+
+
+  MyFST fst(corpusf[0], corpuse[0], &m);
+  ifstream in("./kimura.g");
+  assert(in);
+  CFG_WFSTComposer comp(fst);
+  Hypergraph hg;
+  bool succeed = comp.Compose(&in, &hg);
+  hg.PrintGraphviz();
+  if (succeed) { cerr << "SUCCESS.\n"; } else { cerr << "FAILURE REPORTED.\n"; }
+
+#if 0
+  ifstream in2("./amnabooks.g");
+  assert(in2);
+  MyFST fst2(corpusf[1], corpuse[1], &m);
+  CFG_WFSTComposer comp2(fst2);
+  Hypergraph hg2;
+  bool succeed2 = comp2.Compose(&in2, &hg2);
+  if (succeed2) { cerr << "SUCCESS.\n"; } else { cerr << "FAILURE REPORTED.\n"; }
+#endif
+
+  SparseVector<double> w; w.set_value(FD::Convert("Proposal"), 1.0);
+  hg.Reweight(w);
+  cerr << ViterbiFTree(hg) << endl;
+  return 0;
+}
+
diff --git a/gi/pf/pfdist.cc b/gi/pf/pfdist.cc
new file mode 100644
index 00000000..18dfd03b
--- /dev/null
+++ b/gi/pf/pfdist.cc
@@ -0,0 +1,621 @@
+#include <iostream>
+#include <tr1/memory>
+#include <queue>
+
+#include <boost/functional.hpp>
+#include <boost/program_options.hpp>
+#include <boost/program_options/variables_map.hpp>
+
+#include "base_measures.h"
+#include "reachability.h"
+#include "viterbi.h"
+#include "hg.h"
+#include "trule.h"
+#include "tdict.h"
+#include "filelib.h"
+#include "dict.h"
+#include "sampler.h"
+#include "ccrp_nt.h"
+#include "ccrp_onetable.h"
+
+using namespace std;
+using namespace tr1;
+namespace po = boost::program_options;
+
+shared_ptr<MT19937> prng;
+
+size_t hash_value(const TRule& r) {
+  size_t h = boost::hash_value(r.e_);
+  boost::hash_combine(h, -r.lhs_);
+  boost::hash_combine(h, boost::hash_value(r.f_));
+  return h;
+}
+
+bool operator==(const TRule& a, const TRule& b) {
+  return (a.lhs_ == b.lhs_ && a.e_ == b.e_ && a.f_ == b.f_);
+}
+
+void InitCommandLine(int argc, char** argv, po::variables_map* conf) {
+  po::options_description opts("Configuration options");
+  opts.add_options()
+        ("samples,s",po::value<unsigned>()->default_value(1000),"Number of samples")
+        ("particles,p",po::value<unsigned>()->default_value(30),"Number of particles")
+        ("filter_frequency,f",po::value<unsigned>()->default_value(5),"Number of time steps between filterings")
+        ("input,i",po::value<string>(),"Read parallel data from")
+        ("max_src_phrase",po::value<unsigned>()->default_value(5),"Maximum length of source language phrases")
+        ("max_trg_phrase",po::value<unsigned>()->default_value(5),"Maximum length of target language phrases")
+        ("model1,m",po::value<string>(),"Model 1 parameters (used in base distribution)")
+        ("inverse_model1,M",po::value<string>(),"Inverse Model 1 parameters (used in backward estimate)")
+        ("model1_interpolation_weight",po::value<double>()->default_value(0.95),"Mixing proportion of model 1 with uniform target distribution")
+        ("random_seed,S",po::value<uint32_t>(), "Random seed");
+  po::options_description clo("Command line options");
+  clo.add_options()
+        ("config", po::value<string>(), "Configuration file")
+        ("help,h", "Print this help message and exit");
+  po::options_description dconfig_options, dcmdline_options;
+  dconfig_options.add(opts);
+  dcmdline_options.add(opts).add(clo);
+  
+  po::store(parse_command_line(argc, argv, dcmdline_options), *conf);
+  if (conf->count("config")) {
+    ifstream config((*conf)["config"].as<string>().c_str());
+    po::store(po::parse_config_file(config, dconfig_options), *conf);
+  }
+  po::notify(*conf);
+
+  if (conf->count("help") || (conf->count("input") == 0)) {
+    cerr << dcmdline_options << endl;
+    exit(1);
+  }
+}
+
+void ReadParallelCorpus(const string& filename,
+                vector<vector<WordID> >* f,
+                vector<vector<WordID> >* e,
+                set<WordID>* vocab_f,
+                set<WordID>* vocab_e) {
+  f->clear();
+  e->clear();
+  vocab_f->clear();
+  vocab_e->clear();
+  istream* in;
+  if (filename == "-")
+    in = &cin;
+  else
+    in = new ifstream(filename.c_str());
+  assert(*in);
+  string line;
+  const WordID kDIV = TD::Convert("|||");
+  vector<WordID> tmp;
+  while(*in) {
+    getline(*in, line);
+    if (line.empty() && !*in) break;
+    e->push_back(vector<int>());
+    f->push_back(vector<int>());
+    vector<int>& le = e->back();
+    vector<int>& lf = f->back();
+    tmp.clear();
+    TD::ConvertSentence(line, &tmp);
+    bool isf = true;
+    for (unsigned i = 0; i < tmp.size(); ++i) {
+      const int cur = tmp[i];
+      if (isf) {
+        if (kDIV == cur) { isf = false; } else {
+          lf.push_back(cur);
+          vocab_f->insert(cur);
+        }
+      } else {
+        assert(cur != kDIV);
+        le.push_back(cur);
+        vocab_e->insert(cur);
+      }
+    }
+    assert(isf == false);
+  }
+  if (in != &cin) delete in;
+}
+
+#if 0
+struct MyConditionalModel {
+  MyConditionalModel(PhraseConditionalBase& rcp0) : rp0(&rcp0), base(prob_t::One()), src_phrases(1,1), src_jumps(200, CCRP_NoTable<int>(1,1)) {}
+
+  prob_t srcp0(const vector<WordID>& src) const {
+    prob_t p(1.0 / 3000.0);
+    p.poweq(src.size());
+    prob_t lenp; lenp.logeq(log_poisson(src.size(), 1.0));
+    p *= lenp;
+    return p;
+  }
+
+  void DecrementRule(const TRule& rule) {
+    const RuleCRPMap::iterator it = rules.find(rule.f_);
+    assert(it != rules.end());
+    if (it->second.decrement(rule)) {
+      base /= (*rp0)(rule);
+      if (it->second.num_customers() == 0)
+        rules.erase(it);
+    }
+    if (src_phrases.decrement(rule.f_))
+      base /= srcp0(rule.f_);
+  }
+
+  void IncrementRule(const TRule& rule) {
+    RuleCRPMap::iterator it = rules.find(rule.f_);
+    if (it == rules.end())
+      it = rules.insert(make_pair(rule.f_, CCRP_NoTable<TRule>(1,1))).first;
+    if (it->second.increment(rule)) {
+      base *= (*rp0)(rule);
+    }
+    if (src_phrases.increment(rule.f_))
+      base *= srcp0(rule.f_);
+  }
+
+  void IncrementRules(const vector<TRulePtr>& rules) {
+    for (int i = 0; i < rules.size(); ++i)
+      IncrementRule(*rules[i]);
+  }
+
+  void DecrementRules(const vector<TRulePtr>& rules) {
+    for (int i = 0; i < rules.size(); ++i)
+      DecrementRule(*rules[i]);
+  }
+
+  void IncrementJump(int dist, unsigned src_len) {
+    assert(src_len > 0);
+    if (src_jumps[src_len].increment(dist))
+      base *= jp0(dist, src_len);
+  }
+
+  void DecrementJump(int dist, unsigned src_len) {
+    assert(src_len > 0);
+    if (src_jumps[src_len].decrement(dist))
+      base /= jp0(dist, src_len);
+  }
+
+  void IncrementJumps(const vector<int>& js, unsigned src_len) {
+    for (unsigned i = 0; i < js.size(); ++i)
+      IncrementJump(js[i], src_len);
+  }
+
+  void DecrementJumps(const vector<int>& js, unsigned src_len) {
+    for (unsigned i = 0; i < js.size(); ++i)
+      DecrementJump(js[i], src_len);
+  }
+
+  // p(jump = dist | src_len , z)
+  prob_t JumpProbability(int dist, unsigned src_len) {
+    const prob_t p0 = jp0(dist, src_len);
+    const double lp = src_jumps[src_len].logprob(dist, log(p0));
+    prob_t q; q.logeq(lp);
+    return q;
+  }
+
+  // p(rule.f_ | z) * p(rule.e_ | rule.f_ , z)
+  prob_t RuleProbability(const TRule& rule) const {
+    const prob_t p0 = (*rp0)(rule);
+    prob_t srcp; srcp.logeq(src_phrases.logprob(rule.f_, log(srcp0(rule.f_))));
+    const RuleCRPMap::const_iterator it = rules.find(rule.f_);
+    if (it == rules.end()) return srcp * p0;
+    const double lp = it->second.logprob(rule, log(p0));
+    prob_t q; q.logeq(lp);
+    return q * srcp;
+  }
+
+  prob_t Likelihood() const {
+    prob_t p = base;
+    for (RuleCRPMap::const_iterator it = rules.begin();
+         it != rules.end(); ++it) {
+      prob_t cl; cl.logeq(it->second.log_crp_prob());
+      p *= cl;
+    }
+    for (unsigned l = 1; l < src_jumps.size(); ++l) {
+      if (src_jumps[l].num_customers() > 0) {
+        prob_t q;
+        q.logeq(src_jumps[l].log_crp_prob());
+        p *= q;
+      }
+    }
+    return p;
+  }
+
+  JumpBase jp0;
+  const PhraseConditionalBase* rp0;
+  prob_t base;
+  typedef unordered_map<vector<WordID>, CCRP_NoTable<TRule>, boost::hash<vector<WordID> > > RuleCRPMap;
+  RuleCRPMap rules;
+  CCRP_NoTable<vector<WordID> > src_phrases;
+  vector<CCRP_NoTable<int> > src_jumps;
+};
+
+#endif
+
+struct MyJointModel {
+  MyJointModel(PhraseJointBase& rcp0) :
+    rp0(rcp0), base(prob_t::One()), rules(1,1), src_jumps(200, CCRP_NoTable<int>(1,1)) {}
+
+  void DecrementRule(const TRule& rule) {
+    if (rules.decrement(rule))
+      base /= rp0(rule);
+  }
+
+  void IncrementRule(const TRule& rule) {
+    if (rules.increment(rule))
+      base *= rp0(rule);
+  }
+
+  void IncrementRules(const vector<TRulePtr>& rules) {
+    for (int i = 0; i < rules.size(); ++i)
+      IncrementRule(*rules[i]);
+  }
+
+  void DecrementRules(const vector<TRulePtr>& rules) {
+    for (int i = 0; i < rules.size(); ++i)
+      DecrementRule(*rules[i]);
+  }
+
+  void IncrementJump(int dist, unsigned src_len) {
+    assert(src_len > 0);
+    if (src_jumps[src_len].increment(dist))
+      base *= jp0(dist, src_len);
+  }
+
+  void DecrementJump(int dist, unsigned src_len) {
+    assert(src_len > 0);
+    if (src_jumps[src_len].decrement(dist))
+      base /= jp0(dist, src_len);
+  }
+
+  void IncrementJumps(const vector<int>& js, unsigned src_len) {
+    for (unsigned i = 0; i < js.size(); ++i)
+      IncrementJump(js[i], src_len);
+  }
+
+  void DecrementJumps(const vector<int>& js, unsigned src_len) {
+    for (unsigned i = 0; i < js.size(); ++i)
+      DecrementJump(js[i], src_len);
+  }
+
+  // p(jump = dist | src_len , z)
+  prob_t JumpProbability(int dist, unsigned src_len) {
+    const prob_t p0 = jp0(dist, src_len);
+    const double lp = src_jumps[src_len].logprob(dist, log(p0));
+    prob_t q; q.logeq(lp);
+    return q;
+  }
+
+  // p(rule.f_ | z) * p(rule.e_ | rule.f_ , z)
+  prob_t RuleProbability(const TRule& rule) const {
+    prob_t p; p.logeq(rules.logprob(rule, log(rp0(rule))));
+    return p;
+  }
+
+  prob_t Likelihood() const {
+    prob_t p = base;
+    prob_t q; q.logeq(rules.log_crp_prob());
+    p *= q;
+    for (unsigned l = 1; l < src_jumps.size(); ++l) {
+      if (src_jumps[l].num_customers() > 0) {
+        prob_t q;
+        q.logeq(src_jumps[l].log_crp_prob());
+        p *= q;
+      }
+    }
+    return p;
+  }
+
+  JumpBase jp0;
+  const PhraseJointBase& rp0;
+  prob_t base;
+  CCRP_NoTable<TRule> rules;
+  vector<CCRP_NoTable<int> > src_jumps;
+};
+
+struct BackwardEstimate {
+  BackwardEstimate(const Model1& m1, const vector<WordID>& src, const vector<WordID>& trg) :
+      model1_(m1), src_(src), trg_(trg) {
+  }
+  const prob_t& operator()(const vector<bool>& src_cov, unsigned trg_cov) const {
+    assert(src_.size() == src_cov.size());
+    assert(trg_cov <= trg_.size());
+    prob_t& e = cache_[src_cov][trg_cov];
+    if (e.is_0()) {
+      if (trg_cov == trg_.size()) { e = prob_t::One(); return e; }
+      vector<WordID> r(src_.size() + 1); r.clear();
+      r.push_back(0);  // NULL word
+      for (int i = 0; i < src_cov.size(); ++i)
+        if (!src_cov[i]) r.push_back(src_[i]);
+      const prob_t uniform_alignment(1.0 / r.size());
+      e.logeq(log_poisson(trg_.size() - trg_cov, r.size() - 1)); // p(trg len remaining | src len remaining)
+      for (unsigned j = trg_cov; j < trg_.size(); ++j) {
+        prob_t p;
+        for (unsigned i = 0; i < r.size(); ++i)
+          p += model1_(r[i], trg_[j]);
+        if (p.is_0()) {
+          cerr << "ERROR: p(" << TD::Convert(trg_[j]) << " | " << TD::GetString(r) << ") = 0!\n";
+          abort();
+        }
+        p *= uniform_alignment;
+        e *= p;
+      }
+    }
+    return e;
+  }
+  const Model1& model1_;
+  const vector<WordID>& src_;
+  const vector<WordID>& trg_;
+  mutable unordered_map<vector<bool>, map<unsigned, prob_t>, boost::hash<vector<bool> > > cache_;
+};
+
+struct BackwardEstimateSym {
+  BackwardEstimateSym(const Model1& m1,
+                      const Model1& invm1, const vector<WordID>& src, const vector<WordID>& trg) :
+      model1_(m1), invmodel1_(invm1), src_(src), trg_(trg) {
+  }
+  const prob_t& operator()(const vector<bool>& src_cov, unsigned trg_cov) const {
+    assert(src_.size() == src_cov.size());
+    assert(trg_cov <= trg_.size());
+    prob_t& e = cache_[src_cov][trg_cov];
+    if (e.is_0()) {
+      if (trg_cov == trg_.size()) { e = prob_t::One(); return e; }
+      vector<WordID> r(src_.size() + 1); r.clear();
+      for (int i = 0; i < src_cov.size(); ++i)
+        if (!src_cov[i]) r.push_back(src_[i]);
+      r.push_back(0);  // NULL word
+      const prob_t uniform_alignment(1.0 / r.size());
+      e.logeq(log_poisson(trg_.size() - trg_cov, r.size() - 1)); // p(trg len remaining | src len remaining)
+      for (unsigned j = trg_cov; j < trg_.size(); ++j) {
+        prob_t p;
+        for (unsigned i = 0; i < r.size(); ++i)
+          p += model1_(r[i], trg_[j]);
+        if (p.is_0()) {
+          cerr << "ERROR: p(" << TD::Convert(trg_[j]) << " | " << TD::GetString(r) << ") = 0!\n";
+          abort();
+        }
+        p *= uniform_alignment;
+        e *= p;
+      }
+      r.pop_back();
+      const prob_t inv_uniform(1.0 / (trg_.size() - trg_cov + 1.0));
+      prob_t inv;
+      inv.logeq(log_poisson(r.size(), trg_.size() - trg_cov));
+      for (unsigned i = 0; i < r.size(); ++i) {
+        prob_t p;
+        for (unsigned j = trg_cov - 1; j < trg_.size(); ++j)
+          p += invmodel1_(j < trg_cov ? 0 : trg_[j], r[i]);
+        if (p.is_0()) {
+          cerr << "ERROR: p_inv(" << TD::Convert(r[i]) << " | " << TD::GetString(trg_) << ") = 0!\n";
+          abort();
+        }
+        p *= inv_uniform;
+        inv *= p;
+      }
+      prob_t x = pow(e * inv, 0.5);
+      e = x;
+      //cerr << "Forward: " << log(e) << "\tBackward: " << log(inv) << "\t prop: " << log(x) << endl;
+    }
+    return e;
+  }
+  const Model1& model1_;
+  const Model1& invmodel1_;
+  const vector<WordID>& src_;
+  const vector<WordID>& trg_;
+  mutable unordered_map<vector<bool>, map<unsigned, prob_t>, boost::hash<vector<bool> > > cache_;
+};
+
+struct Particle {
+  Particle() : weight(prob_t::One()), src_cov(), trg_cov(), prev_pos(-1) {}
+  prob_t weight;
+  prob_t gamma_last;
+  vector<int> src_jumps;
+  vector<TRulePtr> rules;
+  vector<bool> src_cv;
+  int src_cov;
+  int trg_cov;
+  int prev_pos;
+};
+
+ostream& operator<<(ostream& o, const vector<bool>& v) {
+  for (int i = 0; i < v.size(); ++i)
+    o << (v[i] ? '1' : '0');
+  return o;
+}
+ostream& operator<<(ostream& o, const Particle& p) {
+  o << "[cv=" << p.src_cv << "  src_cov=" << p.src_cov << " trg_cov=" << p.trg_cov << " last_pos=" << p.prev_pos << " num_rules=" << p.rules.size() << "  w=" << log(p.weight) << ']';
+  return o;
+}
+
+void FilterCrapParticlesAndReweight(vector<Particle>* pps) {
+  vector<Particle>& ps = *pps;
+  SampleSet<prob_t> ss;
+  for (int i = 0; i < ps.size(); ++i)
+    ss.add(ps[i].weight);
+  vector<Particle> nps; nps.reserve(ps.size());
+  const prob_t uniform_weight(1.0 / ps.size());
+  for (int i = 0; i < ps.size(); ++i) {
+    nps.push_back(ps[prng->SelectSample(ss)]);
+    nps[i].weight = uniform_weight;
+  }
+  nps.swap(ps);
+}
+
+int main(int argc, char** argv) {
+  po::variables_map conf;
+  InitCommandLine(argc, argv, &conf);
+  const unsigned kMAX_TRG_PHRASE = conf["max_trg_phrase"].as<unsigned>();
+  const unsigned kMAX_SRC_PHRASE = conf["max_src_phrase"].as<unsigned>();
+  const unsigned particles = conf["particles"].as<unsigned>();
+  const unsigned samples = conf["samples"].as<unsigned>();
+  const unsigned rejuv_freq = conf["filter_frequency"].as<unsigned>();
+
+  if (!conf.count("model1")) {
+    cerr << argv[0] << "Please use --model1 to specify model 1 parameters\n";
+    return 1;
+  }
+  if (conf.count("random_seed"))
+    prng.reset(new MT19937(conf["random_seed"].as<uint32_t>()));
+  else
+    prng.reset(new MT19937);
+  MT19937& rng = *prng;
+
+  vector<vector<WordID> > corpuse, corpusf;
+  set<WordID> vocabe, vocabf;
+  cerr << "Reading corpus...\n";
+  ReadParallelCorpus(conf["input"].as<string>(), &corpusf, &corpuse, &vocabf, &vocabe);
+  cerr << "F-corpus size: " << corpusf.size() << " sentences\t (" << vocabf.size() << " word types)\n";
+  cerr << "E-corpus size: " << corpuse.size() << " sentences\t (" << vocabe.size() << " word types)\n";
+  assert(corpusf.size() == corpuse.size());
+
+  const int kLHS = -TD::Convert("X");
+  Model1 m1(conf["model1"].as<string>());
+  Model1 invm1(conf["inverse_model1"].as<string>());
+
+#if 0
+  PhraseConditionalBase lp0(m1, conf["model1_interpolation_weight"].as<double>(), vocabe.size());
+  MyConditionalModel m(lp0);
+#else
+  PhraseJointBase lp0(m1, conf["model1_interpolation_weight"].as<double>(), vocabe.size(), vocabf.size());
+  MyJointModel m(lp0);
+#endif
+
+  cerr << "Initializing reachability limits...\n";
+  vector<Particle> ps(corpusf.size());
+  vector<Reachability> reaches; reaches.reserve(corpusf.size());
+  for (int ci = 0; ci < corpusf.size(); ++ci)
+    reaches.push_back(Reachability(corpusf[ci].size(),
+                                   corpuse[ci].size(),
+                                   kMAX_SRC_PHRASE,
+                                   kMAX_TRG_PHRASE));
+  cerr << "Sampling...\n"; 
+  vector<Particle> tmp_p(10000);  // work space
+  SampleSet<prob_t> pfss;
+  for (int SS=0; SS < samples; ++SS) {
+    for (int ci = 0; ci < corpusf.size(); ++ci) {
+      vector<int>& src = corpusf[ci];
+      vector<int>& trg = corpuse[ci];
+      m.DecrementRules(ps[ci].rules);
+      m.DecrementJumps(ps[ci].src_jumps, src.size());
+
+      //BackwardEstimate be(m1, src, trg);
+      BackwardEstimateSym be(m1, invm1, src, trg);
+      const Reachability& r = reaches[ci];
+      vector<Particle> lps(particles);
+
+      for (int pi = 0; pi < particles; ++pi) {
+        Particle& p = lps[pi];
+        p.src_cv.resize(src.size(), false);
+      }
+
+      bool all_complete = false;
+      while(!all_complete) {
+        SampleSet<prob_t> ss;
+
+        // all particles have now been extended a bit, we will reweight them now
+        if (lps[0].trg_cov > 0)
+          FilterCrapParticlesAndReweight(&lps);
+
+        // loop over all particles and extend them
+        bool done_nothing = true;
+        for (int pi = 0; pi < particles; ++pi) {
+          Particle& p = lps[pi];
+          int tic = 0;
+          while(p.trg_cov < trg.size() && tic < rejuv_freq) {
+            ++tic;
+            done_nothing = false;
+            ss.clear();
+            TRule x; x.lhs_ = kLHS;
+            prob_t z;
+            int first_uncovered = src.size();
+            int last_uncovered = -1;
+            for (int i = 0; i < src.size(); ++i) {
+              const bool is_uncovered = !p.src_cv[i];
+              if (i < first_uncovered && is_uncovered) first_uncovered = i;
+              if (is_uncovered && i > last_uncovered) last_uncovered = i;
+            }
+            assert(last_uncovered > -1);
+            assert(first_uncovered < src.size());
+
+            for (int trg_len = 1; trg_len <= kMAX_TRG_PHRASE; ++trg_len) {
+              x.e_.push_back(trg[trg_len - 1 + p.trg_cov]);
+              for (int src_len = 1; src_len <= kMAX_SRC_PHRASE; ++src_len) {
+                if (!r.edges[p.src_cov][p.trg_cov][src_len][trg_len]) continue;
+
+                const int last_possible_start = last_uncovered - src_len + 1;
+                assert(last_possible_start >= 0);
+                //cerr << src_len << "," << trg_len << " is allowed. E=" << TD::GetString(x.e_) << endl;
+                //cerr << "  first_uncovered=" << first_uncovered << "  last_possible_start=" << last_possible_start << endl;
+                for (int i = first_uncovered; i <= last_possible_start; ++i) {
+                  if (p.src_cv[i]) continue;
+                  assert(ss.size() < tmp_p.size());  // if fails increase tmp_p size
+                  Particle& np = tmp_p[ss.size()];
+                  np = p;
+                  x.f_.clear();
+                  int gap_add = 0;
+                  bool bad = false;
+                  prob_t jp = prob_t::One();
+                  int prev_pos = p.prev_pos;
+                  for (int j = 0; j < src_len; ++j) {
+                    if ((j + i + gap_add) == src.size()) { bad = true; break; }
+                    while ((i+j+gap_add) < src.size() && p.src_cv[i + j + gap_add]) { ++gap_add; }
+                    if ((j + i + gap_add) == src.size()) { bad = true; break; }
+                    np.src_cv[i + j + gap_add] = true;
+                    x.f_.push_back(src[i + j + gap_add]);
+                    jp *= m.JumpProbability(i + j + gap_add - prev_pos, src.size());
+                    int jump = i + j + gap_add - prev_pos;
+                    assert(jump != 0);
+                    np.src_jumps.push_back(jump);
+                    prev_pos = i + j + gap_add;
+                  }
+                  if (bad) continue;
+                  np.prev_pos = prev_pos;
+                  np.src_cov += x.f_.size();
+                  np.trg_cov += x.e_.size();
+                  if (x.f_.size() != src_len) continue;
+                  prob_t rp = m.RuleProbability(x);
+                  np.gamma_last = rp * jp;
+                  const prob_t u = pow(np.gamma_last * be(np.src_cv, np.trg_cov), 0.2);
+                  //cerr << "**rule=" << x << endl;
+                  //cerr << "  u=" << log(u) << "  rule=" << rp << " jump=" << jp << endl;
+                  ss.add(u);
+                  np.rules.push_back(TRulePtr(new TRule(x)));
+                  z += u;
+
+                  const bool completed = (p.trg_cov == trg.size());
+                  if (completed) {
+                    int last_jump = src.size() - p.prev_pos;
+                    assert(last_jump > 0);
+                    p.src_jumps.push_back(last_jump);
+                    p.weight *= m.JumpProbability(last_jump, src.size());
+                  }
+                }
+              }
+            }
+            cerr << "number of edges to consider: " << ss.size() << endl;
+            const int sampled = rng.SelectSample(ss);
+            prob_t q_n = ss[sampled] / z;
+            p = tmp_p[sampled];
+            //m.IncrementRule(*p.rules.back());
+            p.weight *= p.gamma_last / q_n;
+            cerr << "[w=" << log(p.weight) << "]\tsampled rule: " << p.rules.back()->AsString() << endl;
+            cerr << p << endl;
+          }
+        } // loop over particles (pi = 0 .. particles)
+        if (done_nothing) all_complete = true;
+      }
+      pfss.clear();
+      for (int i = 0; i < lps.size(); ++i)
+        pfss.add(lps[i].weight);
+      const int sampled = rng.SelectSample(pfss);
+      ps[ci] = lps[sampled];
+      m.IncrementRules(lps[sampled].rules);
+      m.IncrementJumps(lps[sampled].src_jumps, src.size());
+      for (int i = 0; i < lps[sampled].rules.size(); ++i) { cerr << "S:\t" << lps[sampled].rules[i]->AsString() << "\n"; }
+      cerr << "tmp-LLH: " << log(m.Likelihood()) << endl;
+    }
+    cerr << "LLH: " << log(m.Likelihood()) << endl;
+    for (int sni = 0; sni < 5; ++sni) {
+      for (int i = 0; i < ps[sni].rules.size(); ++i) { cerr << "\t" << ps[sni].rules[i]->AsString() << endl; }
+    }
+  }
+  return 0;
+}
+
diff --git a/gi/pf/pfdist.new.cc b/gi/pf/pfdist.new.cc
new file mode 100644
index 00000000..3169eb75
--- /dev/null
+++ b/gi/pf/pfdist.new.cc
@@ -0,0 +1,620 @@
+#include <iostream>
+#include <tr1/memory>
+#include <queue>
+
+#include <boost/functional.hpp>
+#include <boost/program_options.hpp>
+#include <boost/program_options/variables_map.hpp>
+
+#include "base_measures.h"
+#include "reachability.h"
+#include "viterbi.h"
+#include "hg.h"
+#include "trule.h"
+#include "tdict.h"
+#include "filelib.h"
+#include "dict.h"
+#include "sampler.h"
+#include "ccrp_nt.h"
+#include "ccrp_onetable.h"
+
+using namespace std;
+using namespace tr1;
+namespace po = boost::program_options;
+
+shared_ptr<MT19937> prng;
+
+size_t hash_value(const TRule& r) {
+  size_t h = boost::hash_value(r.e_);
+  boost::hash_combine(h, -r.lhs_);
+  boost::hash_combine(h, boost::hash_value(r.f_));
+  return h;
+}
+
+bool operator==(const TRule& a, const TRule& b) {
+  return (a.lhs_ == b.lhs_ && a.e_ == b.e_ && a.f_ == b.f_);
+}
+
+void InitCommandLine(int argc, char** argv, po::variables_map* conf) {
+  po::options_description opts("Configuration options");
+  opts.add_options()
+        ("samples,s",po::value<unsigned>()->default_value(1000),"Number of samples")
+        ("particles,p",po::value<unsigned>()->default_value(25),"Number of particles")
+        ("input,i",po::value<string>(),"Read parallel data from")
+        ("max_src_phrase",po::value<unsigned>()->default_value(5),"Maximum length of source language phrases")
+        ("max_trg_phrase",po::value<unsigned>()->default_value(5),"Maximum length of target language phrases")
+        ("model1,m",po::value<string>(),"Model 1 parameters (used in base distribution)")
+        ("inverse_model1,M",po::value<string>(),"Inverse Model 1 parameters (used in backward estimate)")
+        ("model1_interpolation_weight",po::value<double>()->default_value(0.95),"Mixing proportion of model 1 with uniform target distribution")
+        ("random_seed,S",po::value<uint32_t>(), "Random seed");
+  po::options_description clo("Command line options");
+  clo.add_options()
+        ("config", po::value<string>(), "Configuration file")
+        ("help,h", "Print this help message and exit");
+  po::options_description dconfig_options, dcmdline_options;
+  dconfig_options.add(opts);
+  dcmdline_options.add(opts).add(clo);
+  
+  po::store(parse_command_line(argc, argv, dcmdline_options), *conf);
+  if (conf->count("config")) {
+    ifstream config((*conf)["config"].as<string>().c_str());
+    po::store(po::parse_config_file(config, dconfig_options), *conf);
+  }
+  po::notify(*conf);
+
+  if (conf->count("help") || (conf->count("input") == 0)) {
+    cerr << dcmdline_options << endl;
+    exit(1);
+  }
+}
+
+void ReadParallelCorpus(const string& filename,
+                vector<vector<WordID> >* f,
+                vector<vector<WordID> >* e,
+                set<WordID>* vocab_f,
+                set<WordID>* vocab_e) {
+  f->clear();
+  e->clear();
+  vocab_f->clear();
+  vocab_e->clear();
+  istream* in;
+  if (filename == "-")
+    in = &cin;
+  else
+    in = new ifstream(filename.c_str());
+  assert(*in);
+  string line;
+  const WordID kDIV = TD::Convert("|||");
+  vector<WordID> tmp;
+  while(*in) {
+    getline(*in, line);
+    if (line.empty() && !*in) break;
+    e->push_back(vector<int>());
+    f->push_back(vector<int>());
+    vector<int>& le = e->back();
+    vector<int>& lf = f->back();
+    tmp.clear();
+    TD::ConvertSentence(line, &tmp);
+    bool isf = true;
+    for (unsigned i = 0; i < tmp.size(); ++i) {
+      const int cur = tmp[i];
+      if (isf) {
+        if (kDIV == cur) { isf = false; } else {
+          lf.push_back(cur);
+          vocab_f->insert(cur);
+        }
+      } else {
+        assert(cur != kDIV);
+        le.push_back(cur);
+        vocab_e->insert(cur);
+      }
+    }
+    assert(isf == false);
+  }
+  if (in != &cin) delete in;
+}
+
+#if 0
+struct MyConditionalModel {
+  MyConditionalModel(PhraseConditionalBase& rcp0) : rp0(&rcp0), base(prob_t::One()), src_phrases(1,1), src_jumps(200, CCRP_NoTable<int>(1,1)) {}
+
+  prob_t srcp0(const vector<WordID>& src) const {
+    prob_t p(1.0 / 3000.0);
+    p.poweq(src.size());
+    prob_t lenp; lenp.logeq(log_poisson(src.size(), 1.0));
+    p *= lenp;
+    return p;
+  }
+
+  void DecrementRule(const TRule& rule) {
+    const RuleCRPMap::iterator it = rules.find(rule.f_);
+    assert(it != rules.end());
+    if (it->second.decrement(rule)) {
+      base /= (*rp0)(rule);
+      if (it->second.num_customers() == 0)
+        rules.erase(it);
+    }
+    if (src_phrases.decrement(rule.f_))
+      base /= srcp0(rule.f_);
+  }
+
+  void IncrementRule(const TRule& rule) {
+    RuleCRPMap::iterator it = rules.find(rule.f_);
+    if (it == rules.end())
+      it = rules.insert(make_pair(rule.f_, CCRP_NoTable<TRule>(1,1))).first;
+    if (it->second.increment(rule)) {
+      base *= (*rp0)(rule);
+    }
+    if (src_phrases.increment(rule.f_))
+      base *= srcp0(rule.f_);
+  }
+
+  void IncrementRules(const vector<TRulePtr>& rules) {
+    for (int i = 0; i < rules.size(); ++i)
+      IncrementRule(*rules[i]);
+  }
+
+  void DecrementRules(const vector<TRulePtr>& rules) {
+    for (int i = 0; i < rules.size(); ++i)
+      DecrementRule(*rules[i]);
+  }
+
+  void IncrementJump(int dist, unsigned src_len) {
+    assert(src_len > 0);
+    if (src_jumps[src_len].increment(dist))
+      base *= jp0(dist, src_len);
+  }
+
+  void DecrementJump(int dist, unsigned src_len) {
+    assert(src_len > 0);
+    if (src_jumps[src_len].decrement(dist))
+      base /= jp0(dist, src_len);
+  }
+
+  void IncrementJumps(const vector<int>& js, unsigned src_len) {
+    for (unsigned i = 0; i < js.size(); ++i)
+      IncrementJump(js[i], src_len);
+  }
+
+  void DecrementJumps(const vector<int>& js, unsigned src_len) {
+    for (unsigned i = 0; i < js.size(); ++i)
+      DecrementJump(js[i], src_len);
+  }
+
+  // p(jump = dist | src_len , z)
+  prob_t JumpProbability(int dist, unsigned src_len) {
+    const prob_t p0 = jp0(dist, src_len);
+    const double lp = src_jumps[src_len].logprob(dist, log(p0));
+    prob_t q; q.logeq(lp);
+    return q;
+  }
+
+  // p(rule.f_ | z) * p(rule.e_ | rule.f_ , z)
+  prob_t RuleProbability(const TRule& rule) const {
+    const prob_t p0 = (*rp0)(rule);
+    prob_t srcp; srcp.logeq(src_phrases.logprob(rule.f_, log(srcp0(rule.f_))));
+    const RuleCRPMap::const_iterator it = rules.find(rule.f_);
+    if (it == rules.end()) return srcp * p0;
+    const double lp = it->second.logprob(rule, log(p0));
+    prob_t q; q.logeq(lp);
+    return q * srcp;
+  }
+
+  prob_t Likelihood() const {
+    prob_t p = base;
+    for (RuleCRPMap::const_iterator it = rules.begin();
+         it != rules.end(); ++it) {
+      prob_t cl; cl.logeq(it->second.log_crp_prob());
+      p *= cl;
+    }
+    for (unsigned l = 1; l < src_jumps.size(); ++l) {
+      if (src_jumps[l].num_customers() > 0) {
+        prob_t q;
+        q.logeq(src_jumps[l].log_crp_prob());
+        p *= q;
+      }
+    }
+    return p;
+  }
+
+  JumpBase jp0;
+  const PhraseConditionalBase* rp0;
+  prob_t base;
+  typedef unordered_map<vector<WordID>, CCRP_NoTable<TRule>, boost::hash<vector<WordID> > > RuleCRPMap;
+  RuleCRPMap rules;
+  CCRP_NoTable<vector<WordID> > src_phrases;
+  vector<CCRP_NoTable<int> > src_jumps;
+};
+
+#endif
+
+struct MyJointModel {
+  MyJointModel(PhraseJointBase& rcp0) :
+    rp0(rcp0), base(prob_t::One()), rules(1,1), src_jumps(200, CCRP_NoTable<int>(1,1)) {}
+
+  void DecrementRule(const TRule& rule) {
+    if (rules.decrement(rule))
+      base /= rp0(rule);
+  }
+
+  void IncrementRule(const TRule& rule) {
+    if (rules.increment(rule))
+      base *= rp0(rule);
+  }
+
+  void IncrementRules(const vector<TRulePtr>& rules) {
+    for (int i = 0; i < rules.size(); ++i)
+      IncrementRule(*rules[i]);
+  }
+
+  void DecrementRules(const vector<TRulePtr>& rules) {
+    for (int i = 0; i < rules.size(); ++i)
+      DecrementRule(*rules[i]);
+  }
+
+  void IncrementJump(int dist, unsigned src_len) {
+    assert(src_len > 0);
+    if (src_jumps[src_len].increment(dist))
+      base *= jp0(dist, src_len);
+  }
+
+  void DecrementJump(int dist, unsigned src_len) {
+    assert(src_len > 0);
+    if (src_jumps[src_len].decrement(dist))
+      base /= jp0(dist, src_len);
+  }
+
+  void IncrementJumps(const vector<int>& js, unsigned src_len) {
+    for (unsigned i = 0; i < js.size(); ++i)
+      IncrementJump(js[i], src_len);
+  }
+
+  void DecrementJumps(const vector<int>& js, unsigned src_len) {
+    for (unsigned i = 0; i < js.size(); ++i)
+      DecrementJump(js[i], src_len);
+  }
+
+  // p(jump = dist | src_len , z)
+  prob_t JumpProbability(int dist, unsigned src_len) {
+    const prob_t p0 = jp0(dist, src_len);
+    const double lp = src_jumps[src_len].logprob(dist, log(p0));
+    prob_t q; q.logeq(lp);
+    return q;
+  }
+
+  // p(rule.f_ | z) * p(rule.e_ | rule.f_ , z)
+  prob_t RuleProbability(const TRule& rule) const {
+    prob_t p; p.logeq(rules.logprob(rule, log(rp0(rule))));
+    return p;
+  }
+
+  prob_t Likelihood() const {
+    prob_t p = base;
+    prob_t q; q.logeq(rules.log_crp_prob());
+    p *= q;
+    for (unsigned l = 1; l < src_jumps.size(); ++l) {
+      if (src_jumps[l].num_customers() > 0) {
+        prob_t q;
+        q.logeq(src_jumps[l].log_crp_prob());
+        p *= q;
+      }
+    }
+    return p;
+  }
+
+  JumpBase jp0;
+  const PhraseJointBase& rp0;
+  prob_t base;
+  CCRP_NoTable<TRule> rules;
+  vector<CCRP_NoTable<int> > src_jumps;
+};
+
+struct BackwardEstimate {
+  BackwardEstimate(const Model1& m1, const vector<WordID>& src, const vector<WordID>& trg) :
+      model1_(m1), src_(src), trg_(trg) {
+  }
+  const prob_t& operator()(const vector<bool>& src_cov, unsigned trg_cov) const {
+    assert(src_.size() == src_cov.size());
+    assert(trg_cov <= trg_.size());
+    prob_t& e = cache_[src_cov][trg_cov];
+    if (e.is_0()) {
+      if (trg_cov == trg_.size()) { e = prob_t::One(); return e; }
+      vector<WordID> r(src_.size() + 1); r.clear();
+      r.push_back(0);  // NULL word
+      for (int i = 0; i < src_cov.size(); ++i)
+        if (!src_cov[i]) r.push_back(src_[i]);
+      const prob_t uniform_alignment(1.0 / r.size());
+      e.logeq(log_poisson(trg_.size() - trg_cov, r.size() - 1)); // p(trg len remaining | src len remaining)
+      for (unsigned j = trg_cov; j < trg_.size(); ++j) {
+        prob_t p;
+        for (unsigned i = 0; i < r.size(); ++i)
+          p += model1_(r[i], trg_[j]);
+        if (p.is_0()) {
+          cerr << "ERROR: p(" << TD::Convert(trg_[j]) << " | " << TD::GetString(r) << ") = 0!\n";
+          abort();
+        }
+        p *= uniform_alignment;
+        e *= p;
+      }
+    }
+    return e;
+  }
+  const Model1& model1_;
+  const vector<WordID>& src_;
+  const vector<WordID>& trg_;
+  mutable unordered_map<vector<bool>, map<unsigned, prob_t>, boost::hash<vector<bool> > > cache_;
+};
+
+struct BackwardEstimateSym {
+  BackwardEstimateSym(const Model1& m1,
+                      const Model1& invm1, const vector<WordID>& src, const vector<WordID>& trg) :
+      model1_(m1), invmodel1_(invm1), src_(src), trg_(trg) {
+  }
+  const prob_t& operator()(const vector<bool>& src_cov, unsigned trg_cov) const {
+    assert(src_.size() == src_cov.size());
+    assert(trg_cov <= trg_.size());
+    prob_t& e = cache_[src_cov][trg_cov];
+    if (e.is_0()) {
+      if (trg_cov == trg_.size()) { e = prob_t::One(); return e; }
+      vector<WordID> r(src_.size() + 1); r.clear();
+      for (int i = 0; i < src_cov.size(); ++i)
+        if (!src_cov[i]) r.push_back(src_[i]);
+      r.push_back(0);  // NULL word
+      const prob_t uniform_alignment(1.0 / r.size());
+      e.logeq(log_poisson(trg_.size() - trg_cov, r.size() - 1)); // p(trg len remaining | src len remaining)
+      for (unsigned j = trg_cov; j < trg_.size(); ++j) {
+        prob_t p;
+        for (unsigned i = 0; i < r.size(); ++i)
+          p += model1_(r[i], trg_[j]);
+        if (p.is_0()) {
+          cerr << "ERROR: p(" << TD::Convert(trg_[j]) << " | " << TD::GetString(r) << ") = 0!\n";
+          abort();
+        }
+        p *= uniform_alignment;
+        e *= p;
+      }
+      r.pop_back();
+      const prob_t inv_uniform(1.0 / (trg_.size() - trg_cov + 1.0));
+      prob_t inv;
+      inv.logeq(log_poisson(r.size(), trg_.size() - trg_cov));
+      for (unsigned i = 0; i < r.size(); ++i) {
+        prob_t p;
+        for (unsigned j = trg_cov - 1; j < trg_.size(); ++j)
+          p += invmodel1_(j < trg_cov ? 0 : trg_[j], r[i]);
+        if (p.is_0()) {
+          cerr << "ERROR: p_inv(" << TD::Convert(r[i]) << " | " << TD::GetString(trg_) << ") = 0!\n";
+          abort();
+        }
+        p *= inv_uniform;
+        inv *= p;
+      }
+      prob_t x = pow(e * inv, 0.5);
+      e = x;
+      //cerr << "Forward: " << log(e) << "\tBackward: " << log(inv) << "\t prop: " << log(x) << endl;
+    }
+    return e;
+  }
+  const Model1& model1_;
+  const Model1& invmodel1_;
+  const vector<WordID>& src_;
+  const vector<WordID>& trg_;
+  mutable unordered_map<vector<bool>, map<unsigned, prob_t>, boost::hash<vector<bool> > > cache_;
+};
+
+struct Particle {
+  Particle() : weight(prob_t::One()), src_cov(), trg_cov(), prev_pos(-1) {}
+  prob_t weight;
+  prob_t gamma_last;
+  vector<int> src_jumps;
+  vector<TRulePtr> rules;
+  vector<bool> src_cv;
+  int src_cov;
+  int trg_cov;
+  int prev_pos;
+};
+
+ostream& operator<<(ostream& o, const vector<bool>& v) {
+  for (int i = 0; i < v.size(); ++i)
+    o << (v[i] ? '1' : '0');
+  return o;
+}
+ostream& operator<<(ostream& o, const Particle& p) {
+  o << "[cv=" << p.src_cv << "  src_cov=" << p.src_cov << " trg_cov=" << p.trg_cov << " last_pos=" << p.prev_pos << " num_rules=" << p.rules.size() << "  w=" << log(p.weight) << ']';
+  return o;
+}
+
+void FilterCrapParticlesAndReweight(vector<Particle>* pps) {
+  vector<Particle>& ps = *pps;
+  SampleSet<prob_t> ss;
+  for (int i = 0; i < ps.size(); ++i)
+    ss.add(ps[i].weight);
+  vector<Particle> nps; nps.reserve(ps.size());
+  const prob_t uniform_weight(1.0 / ps.size());
+  for (int i = 0; i < ps.size(); ++i) {
+    nps.push_back(ps[prng->SelectSample(ss)]);
+    nps[i].weight = uniform_weight;
+  }
+  nps.swap(ps);
+}
+
+int main(int argc, char** argv) {
+  po::variables_map conf;
+  InitCommandLine(argc, argv, &conf);
+  const unsigned kMAX_TRG_PHRASE = conf["max_trg_phrase"].as<unsigned>();
+  const unsigned kMAX_SRC_PHRASE = conf["max_src_phrase"].as<unsigned>();
+  const unsigned particles = conf["particles"].as<unsigned>();
+  const unsigned samples = conf["samples"].as<unsigned>();
+
+  if (!conf.count("model1")) {
+    cerr << argv[0] << "Please use --model1 to specify model 1 parameters\n";
+    return 1;
+  }
+  if (conf.count("random_seed"))
+    prng.reset(new MT19937(conf["random_seed"].as<uint32_t>()));
+  else
+    prng.reset(new MT19937);
+  MT19937& rng = *prng;
+
+  vector<vector<WordID> > corpuse, corpusf;
+  set<WordID> vocabe, vocabf;
+  cerr << "Reading corpus...\n";
+  ReadParallelCorpus(conf["input"].as<string>(), &corpusf, &corpuse, &vocabf, &vocabe);
+  cerr << "F-corpus size: " << corpusf.size() << " sentences\t (" << vocabf.size() << " word types)\n";
+  cerr << "E-corpus size: " << corpuse.size() << " sentences\t (" << vocabe.size() << " word types)\n";
+  assert(corpusf.size() == corpuse.size());
+
+  const int kLHS = -TD::Convert("X");
+  Model1 m1(conf["model1"].as<string>());
+  Model1 invm1(conf["inverse_model1"].as<string>());
+
+#if 0
+  PhraseConditionalBase lp0(m1, conf["model1_interpolation_weight"].as<double>(), vocabe.size());
+  MyConditionalModel m(lp0);
+#else
+  PhraseJointBase lp0(m1, conf["model1_interpolation_weight"].as<double>(), vocabe.size(), vocabf.size());
+  MyJointModel m(lp0);
+#endif
+
+  cerr << "Initializing reachability limits...\n";
+  vector<Particle> ps(corpusf.size());
+  vector<Reachability> reaches; reaches.reserve(corpusf.size());
+  for (int ci = 0; ci < corpusf.size(); ++ci)
+    reaches.push_back(Reachability(corpusf[ci].size(),
+                                   corpuse[ci].size(),
+                                   kMAX_SRC_PHRASE,
+                                   kMAX_TRG_PHRASE));
+  cerr << "Sampling...\n"; 
+  vector<Particle> tmp_p(10000);  // work space
+  SampleSet<prob_t> pfss;
+  for (int SS=0; SS < samples; ++SS) {
+    for (int ci = 0; ci < corpusf.size(); ++ci) {
+      vector<int>& src = corpusf[ci];
+      vector<int>& trg = corpuse[ci];
+      m.DecrementRules(ps[ci].rules);
+      m.DecrementJumps(ps[ci].src_jumps, src.size());
+
+      //BackwardEstimate be(m1, src, trg);
+      BackwardEstimateSym be(m1, invm1, src, trg);
+      const Reachability& r = reaches[ci];
+      vector<Particle> lps(particles);
+
+      for (int pi = 0; pi < particles; ++pi) {
+        Particle& p = lps[pi];
+        p.src_cv.resize(src.size(), false);
+      }
+
+      bool all_complete = false;
+      while(!all_complete) {
+        SampleSet<prob_t> ss;
+
+        // all particles have now been extended a bit, we will reweight them now
+        if (lps[0].trg_cov > 0)
+          FilterCrapParticlesAndReweight(&lps);
+
+        // loop over all particles and extend them
+        bool done_nothing = true;
+        for (int pi = 0; pi < particles; ++pi) {
+          Particle& p = lps[pi];
+          int tic = 0;
+          const int rejuv_freq = 1;
+          while(p.trg_cov < trg.size() && tic < rejuv_freq) {
+            ++tic;
+            done_nothing = false;
+            ss.clear();
+            TRule x; x.lhs_ = kLHS;
+            prob_t z;
+            int first_uncovered = src.size();
+            int last_uncovered = -1;
+            for (int i = 0; i < src.size(); ++i) {
+              const bool is_uncovered = !p.src_cv[i];
+              if (i < first_uncovered && is_uncovered) first_uncovered = i;
+              if (is_uncovered && i > last_uncovered) last_uncovered = i;
+            }
+            assert(last_uncovered > -1);
+            assert(first_uncovered < src.size());
+
+            for (int trg_len = 1; trg_len <= kMAX_TRG_PHRASE; ++trg_len) {
+              x.e_.push_back(trg[trg_len - 1 + p.trg_cov]);
+              for (int src_len = 1; src_len <= kMAX_SRC_PHRASE; ++src_len) {
+                if (!r.edges[p.src_cov][p.trg_cov][src_len][trg_len]) continue;
+
+                const int last_possible_start = last_uncovered - src_len + 1;
+                assert(last_possible_start >= 0);
+                //cerr << src_len << "," << trg_len << " is allowed. E=" << TD::GetString(x.e_) << endl;
+                //cerr << "  first_uncovered=" << first_uncovered << "  last_possible_start=" << last_possible_start << endl;
+                for (int i = first_uncovered; i <= last_possible_start; ++i) {
+                  if (p.src_cv[i]) continue;
+                  assert(ss.size() < tmp_p.size());  // if fails increase tmp_p size
+                  Particle& np = tmp_p[ss.size()];
+                  np = p;
+                  x.f_.clear();
+                  int gap_add = 0;
+                  bool bad = false;
+                  prob_t jp = prob_t::One();
+                  int prev_pos = p.prev_pos;
+                  for (int j = 0; j < src_len; ++j) {
+                    if ((j + i + gap_add) == src.size()) { bad = true; break; }
+                    while ((i+j+gap_add) < src.size() && p.src_cv[i + j + gap_add]) { ++gap_add; }
+                    if ((j + i + gap_add) == src.size()) { bad = true; break; }
+                    np.src_cv[i + j + gap_add] = true;
+                    x.f_.push_back(src[i + j + gap_add]);
+                    jp *= m.JumpProbability(i + j + gap_add - prev_pos, src.size());
+                    int jump = i + j + gap_add - prev_pos;
+                    assert(jump != 0);
+                    np.src_jumps.push_back(jump);
+                    prev_pos = i + j + gap_add;
+                  }
+                  if (bad) continue;
+                  np.prev_pos = prev_pos;
+                  np.src_cov += x.f_.size();
+                  np.trg_cov += x.e_.size();
+                  if (x.f_.size() != src_len) continue;
+                  prob_t rp = m.RuleProbability(x);
+                  np.gamma_last = rp * jp;
+                  const prob_t u = pow(np.gamma_last * be(np.src_cv, np.trg_cov), 0.2);
+                  //cerr << "**rule=" << x << endl;
+                  //cerr << "  u=" << log(u) << "  rule=" << rp << " jump=" << jp << endl;
+                  ss.add(u);
+                  np.rules.push_back(TRulePtr(new TRule(x)));
+                  z += u;
+
+                  const bool completed = (p.trg_cov == trg.size());
+                  if (completed) {
+                    int last_jump = src.size() - p.prev_pos;
+                    assert(last_jump > 0);
+                    p.src_jumps.push_back(last_jump);
+                    p.weight *= m.JumpProbability(last_jump, src.size());
+                  }
+                }
+              }
+            }
+            cerr << "number of edges to consider: " << ss.size() << endl;
+            const int sampled = rng.SelectSample(ss);
+            prob_t q_n = ss[sampled] / z;
+            p = tmp_p[sampled];
+            //m.IncrementRule(*p.rules.back());
+            p.weight *= p.gamma_last / q_n;
+            cerr << "[w=" << log(p.weight) << "]\tsampled rule: " << p.rules.back()->AsString() << endl;
+            cerr << p << endl;
+          }
+        } // loop over particles (pi = 0 .. particles)
+        if (done_nothing) all_complete = true;
+      }
+      pfss.clear();
+      for (int i = 0; i < lps.size(); ++i)
+        pfss.add(lps[i].weight);
+      const int sampled = rng.SelectSample(pfss);
+      ps[ci] = lps[sampled];
+      m.IncrementRules(lps[sampled].rules);
+      m.IncrementJumps(lps[sampled].src_jumps, src.size());
+      for (int i = 0; i < lps[sampled].rules.size(); ++i) { cerr << "S:\t" << lps[sampled].rules[i]->AsString() << "\n"; }
+      cerr << "tmp-LLH: " << log(m.Likelihood()) << endl;
+    }
+    cerr << "LLH: " << log(m.Likelihood()) << endl;
+    for (int sni = 0; sni < 5; ++sni) {
+      for (int i = 0; i < ps[sni].rules.size(); ++i) { cerr << "\t" << ps[sni].rules[i]->AsString() << endl; }
+    }
+  }
+  return 0;
+}
+
diff --git a/gi/pf/pfnaive.cc b/gi/pf/pfnaive.cc
new file mode 100644
index 00000000..43c604c3
--- /dev/null
+++ b/gi/pf/pfnaive.cc
@@ -0,0 +1,385 @@
+#include <iostream>
+#include <tr1/memory>
+#include <queue>
+
+#include <boost/functional.hpp>
+#include <boost/program_options.hpp>
+#include <boost/program_options/variables_map.hpp>
+
+#include "base_measures.h"
+#include "reachability.h"
+#include "viterbi.h"
+#include "hg.h"
+#include "trule.h"
+#include "tdict.h"
+#include "filelib.h"
+#include "dict.h"
+#include "sampler.h"
+#include "ccrp_nt.h"
+#include "ccrp_onetable.h"
+
+using namespace std;
+using namespace tr1;
+namespace po = boost::program_options;
+
+shared_ptr<MT19937> prng;
+
+size_t hash_value(const TRule& r) {
+  size_t h = boost::hash_value(r.e_);
+  boost::hash_combine(h, -r.lhs_);
+  boost::hash_combine(h, boost::hash_value(r.f_));
+  return h;
+}
+
+bool operator==(const TRule& a, const TRule& b) {
+  return (a.lhs_ == b.lhs_ && a.e_ == b.e_ && a.f_ == b.f_);
+}
+
+void InitCommandLine(int argc, char** argv, po::variables_map* conf) {
+  po::options_description opts("Configuration options");
+  opts.add_options()
+        ("samples,s",po::value<unsigned>()->default_value(1000),"Number of samples")
+        ("particles,p",po::value<unsigned>()->default_value(30),"Number of particles")
+        ("filter_frequency,f",po::value<unsigned>()->default_value(5),"Number of time steps between filterings")
+        ("input,i",po::value<string>(),"Read parallel data from")
+        ("max_src_phrase",po::value<unsigned>()->default_value(5),"Maximum length of source language phrases")
+        ("max_trg_phrase",po::value<unsigned>()->default_value(5),"Maximum length of target language phrases")
+        ("model1,m",po::value<string>(),"Model 1 parameters (used in base distribution)")
+        ("inverse_model1,M",po::value<string>(),"Inverse Model 1 parameters (used in backward estimate)")
+        ("model1_interpolation_weight",po::value<double>()->default_value(0.95),"Mixing proportion of model 1 with uniform target distribution")
+        ("random_seed,S",po::value<uint32_t>(), "Random seed");
+  po::options_description clo("Command line options");
+  clo.add_options()
+        ("config", po::value<string>(), "Configuration file")
+        ("help,h", "Print this help message and exit");
+  po::options_description dconfig_options, dcmdline_options;
+  dconfig_options.add(opts);
+  dcmdline_options.add(opts).add(clo);
+  
+  po::store(parse_command_line(argc, argv, dcmdline_options), *conf);
+  if (conf->count("config")) {
+    ifstream config((*conf)["config"].as<string>().c_str());
+    po::store(po::parse_config_file(config, dconfig_options), *conf);
+  }
+  po::notify(*conf);
+
+  if (conf->count("help") || (conf->count("input") == 0)) {
+    cerr << dcmdline_options << endl;
+    exit(1);
+  }
+}
+
+void ReadParallelCorpus(const string& filename,
+                vector<vector<WordID> >* f,
+                vector<vector<WordID> >* e,
+                set<WordID>* vocab_f,
+                set<WordID>* vocab_e) {
+  f->clear();
+  e->clear();
+  vocab_f->clear();
+  vocab_e->clear();
+  istream* in;
+  if (filename == "-")
+    in = &cin;
+  else
+    in = new ifstream(filename.c_str());
+  assert(*in);
+  string line;
+  const WordID kDIV = TD::Convert("|||");
+  vector<WordID> tmp;
+  while(*in) {
+    getline(*in, line);
+    if (line.empty() && !*in) break;
+    e->push_back(vector<int>());
+    f->push_back(vector<int>());
+    vector<int>& le = e->back();
+    vector<int>& lf = f->back();
+    tmp.clear();
+    TD::ConvertSentence(line, &tmp);
+    bool isf = true;
+    for (unsigned i = 0; i < tmp.size(); ++i) {
+      const int cur = tmp[i];
+      if (isf) {
+        if (kDIV == cur) { isf = false; } else {
+          lf.push_back(cur);
+          vocab_f->insert(cur);
+        }
+      } else {
+        assert(cur != kDIV);
+        le.push_back(cur);
+        vocab_e->insert(cur);
+      }
+    }
+    assert(isf == false);
+  }
+  if (in != &cin) delete in;
+}
+
+struct MyJointModel {
+  MyJointModel(PhraseJointBase& rcp0) :
+    rp0(rcp0), base(prob_t::One()), rules(1,1) {}
+
+  void DecrementRule(const TRule& rule) {
+    if (rules.decrement(rule))
+      base /= rp0(rule);
+  }
+
+  void IncrementRule(const TRule& rule) {
+    if (rules.increment(rule))
+      base *= rp0(rule);
+  }
+
+  void IncrementRules(const vector<TRulePtr>& rules) {
+    for (int i = 0; i < rules.size(); ++i)
+      IncrementRule(*rules[i]);
+  }
+
+  void DecrementRules(const vector<TRulePtr>& rules) {
+    for (int i = 0; i < rules.size(); ++i)
+      DecrementRule(*rules[i]);
+  }
+
+  prob_t RuleProbability(const TRule& rule) const {
+    prob_t p; p.logeq(rules.logprob(rule, log(rp0(rule))));
+    return p;
+  }
+
+  prob_t Likelihood() const {
+    prob_t p = base;
+    prob_t q; q.logeq(rules.log_crp_prob());
+    p *= q;
+    for (unsigned l = 1; l < src_jumps.size(); ++l) {
+      if (src_jumps[l].num_customers() > 0) {
+        prob_t q;
+        q.logeq(src_jumps[l].log_crp_prob());
+        p *= q;
+      }
+    }
+    return p;
+  }
+
+  const PhraseJointBase& rp0;
+  prob_t base;
+  CCRP_NoTable<TRule> rules;
+  vector<CCRP_NoTable<int> > src_jumps;
+};
+
+struct BackwardEstimateSym {
+  BackwardEstimateSym(const Model1& m1,
+                      const Model1& invm1, const vector<WordID>& src, const vector<WordID>& trg) :
+      model1_(m1), invmodel1_(invm1), src_(src), trg_(trg) {
+  }
+  const prob_t& operator()(unsigned src_cov, unsigned trg_cov) const {
+    assert(src_cov <= src_.size());
+    assert(trg_cov <= trg_.size());
+    prob_t& e = cache_[src_cov][trg_cov];
+    if (e.is_0()) {
+      if (trg_cov == trg_.size()) { e = prob_t::One(); return e; }
+      vector<WordID> r(src_.size() + 1); r.clear();
+      for (int i = src_cov; i < src_.size(); ++i)
+        r.push_back(src_[i]);
+      r.push_back(0);  // NULL word
+      const prob_t uniform_alignment(1.0 / r.size());
+      e.logeq(log_poisson(trg_.size() - trg_cov, r.size() - 1)); // p(trg len remaining | src len remaining)
+      for (unsigned j = trg_cov; j < trg_.size(); ++j) {
+        prob_t p;
+        for (unsigned i = 0; i < r.size(); ++i)
+          p += model1_(r[i], trg_[j]);
+        if (p.is_0()) {
+          cerr << "ERROR: p(" << TD::Convert(trg_[j]) << " | " << TD::GetString(r) << ") = 0!\n";
+          abort();
+        }
+        p *= uniform_alignment;
+        e *= p;
+      }
+      r.pop_back();
+      const prob_t inv_uniform(1.0 / (trg_.size() - trg_cov + 1.0));
+      prob_t inv;
+      inv.logeq(log_poisson(r.size(), trg_.size() - trg_cov));
+      for (unsigned i = 0; i < r.size(); ++i) {
+        prob_t p;
+        for (unsigned j = trg_cov - 1; j < trg_.size(); ++j)
+          p += invmodel1_(j < trg_cov ? 0 : trg_[j], r[i]);
+        if (p.is_0()) {
+          cerr << "ERROR: p_inv(" << TD::Convert(r[i]) << " | " << TD::GetString(trg_) << ") = 0!\n";
+          abort();
+        }
+        p *= inv_uniform;
+        inv *= p;
+      }
+      prob_t x = pow(e * inv, 0.5);
+      e = x;
+      //cerr << "Forward: " << log(e) << "\tBackward: " << log(inv) << "\t prop: " << log(x) << endl;
+    }
+    return e;
+  }
+  const Model1& model1_;
+  const Model1& invmodel1_;
+  const vector<WordID>& src_;
+  const vector<WordID>& trg_;
+  mutable unordered_map<unsigned, map<unsigned, prob_t> > cache_;
+};
+
+struct Particle {
+  Particle() : weight(prob_t::One()), src_cov(), trg_cov() {}
+  prob_t weight;
+  prob_t gamma_last;
+  vector<TRulePtr> rules;
+  int src_cov;
+  int trg_cov;
+};
+
+ostream& operator<<(ostream& o, const vector<bool>& v) {
+  for (int i = 0; i < v.size(); ++i)
+    o << (v[i] ? '1' : '0');
+  return o;
+}
+ostream& operator<<(ostream& o, const Particle& p) {
+  o << "[src_cov=" << p.src_cov << " trg_cov=" << p.trg_cov << " num_rules=" << p.rules.size() << "  w=" << log(p.weight) << ']';
+  return o;
+}
+
+void FilterCrapParticlesAndReweight(vector<Particle>* pps) {
+  vector<Particle>& ps = *pps;
+  SampleSet<prob_t> ss;
+  for (int i = 0; i < ps.size(); ++i)
+    ss.add(ps[i].weight);
+  vector<Particle> nps; nps.reserve(ps.size());
+  const prob_t uniform_weight(1.0 / ps.size());
+  for (int i = 0; i < ps.size(); ++i) {
+    nps.push_back(ps[prng->SelectSample(ss)]);
+    nps[i].weight = uniform_weight;
+  }
+  nps.swap(ps);
+}
+
+int main(int argc, char** argv) {
+  po::variables_map conf;
+  InitCommandLine(argc, argv, &conf);
+  const unsigned kMAX_TRG_PHRASE = conf["max_trg_phrase"].as<unsigned>();
+  const unsigned kMAX_SRC_PHRASE = conf["max_src_phrase"].as<unsigned>();
+  const unsigned particles = conf["particles"].as<unsigned>();
+  const unsigned samples = conf["samples"].as<unsigned>();
+  const unsigned rejuv_freq = conf["filter_frequency"].as<unsigned>();
+
+  if (!conf.count("model1")) {
+    cerr << argv[0] << "Please use --model1 to specify model 1 parameters\n";
+    return 1;
+  }
+  if (conf.count("random_seed"))
+    prng.reset(new MT19937(conf["random_seed"].as<uint32_t>()));
+  else
+    prng.reset(new MT19937);
+  MT19937& rng = *prng;
+
+  vector<vector<WordID> > corpuse, corpusf;
+  set<WordID> vocabe, vocabf;
+  cerr << "Reading corpus...\n";
+  ReadParallelCorpus(conf["input"].as<string>(), &corpusf, &corpuse, &vocabf, &vocabe);
+  cerr << "F-corpus size: " << corpusf.size() << " sentences\t (" << vocabf.size() << " word types)\n";
+  cerr << "E-corpus size: " << corpuse.size() << " sentences\t (" << vocabe.size() << " word types)\n";
+  assert(corpusf.size() == corpuse.size());
+
+  const int kLHS = -TD::Convert("X");
+  Model1 m1(conf["model1"].as<string>());
+  Model1 invm1(conf["inverse_model1"].as<string>());
+
+#if 0
+  PhraseConditionalBase lp0(m1, conf["model1_interpolation_weight"].as<double>(), vocabe.size());
+  MyConditionalModel m(lp0);
+#else
+  PhraseJointBase lp0(m1, conf["model1_interpolation_weight"].as<double>(), vocabe.size(), vocabf.size());
+  MyJointModel m(lp0);
+#endif
+
+  cerr << "Initializing reachability limits...\n";
+  vector<Particle> ps(corpusf.size());
+  vector<Reachability> reaches; reaches.reserve(corpusf.size());
+  for (int ci = 0; ci < corpusf.size(); ++ci)
+    reaches.push_back(Reachability(corpusf[ci].size(),
+                                   corpuse[ci].size(),
+                                   kMAX_SRC_PHRASE,
+                                   kMAX_TRG_PHRASE));
+  cerr << "Sampling...\n"; 
+  vector<Particle> tmp_p(10000);  // work space
+  SampleSet<prob_t> pfss;
+  for (int SS=0; SS < samples; ++SS) {
+    for (int ci = 0; ci < corpusf.size(); ++ci) {
+      vector<int>& src = corpusf[ci];
+      vector<int>& trg = corpuse[ci];
+      m.DecrementRules(ps[ci].rules);
+
+      BackwardEstimateSym be(m1, invm1, src, trg);
+      const Reachability& r = reaches[ci];
+      vector<Particle> lps(particles);
+
+      bool all_complete = false;
+      while(!all_complete) {
+        SampleSet<prob_t> ss;
+
+        // all particles have now been extended a bit, we will reweight them now
+        if (lps[0].trg_cov > 0)
+          FilterCrapParticlesAndReweight(&lps);
+
+        // loop over all particles and extend them
+        bool done_nothing = true;
+        for (int pi = 0; pi < particles; ++pi) {
+          Particle& p = lps[pi];
+          int tic = 0;
+          while(p.trg_cov < trg.size() && tic < rejuv_freq) {
+            ++tic;
+            done_nothing = false;
+            ss.clear();
+            TRule x; x.lhs_ = kLHS;
+            prob_t z;
+
+            for (int trg_len = 1; trg_len <= kMAX_TRG_PHRASE; ++trg_len) {
+              x.e_.push_back(trg[trg_len - 1 + p.trg_cov]);
+              for (int src_len = 1; src_len <= kMAX_SRC_PHRASE; ++src_len) {
+                if (!r.edges[p.src_cov][p.trg_cov][src_len][trg_len]) continue;
+
+                int i = p.src_cov;
+                assert(ss.size() < tmp_p.size());  // if fails increase tmp_p size
+                Particle& np = tmp_p[ss.size()];
+                np = p;
+                x.f_.clear();
+                for (int j = 0; j < src_len; ++j)
+                  x.f_.push_back(src[i + j]);
+                np.src_cov += x.f_.size();
+                np.trg_cov += x.e_.size();
+                prob_t rp = m.RuleProbability(x);
+                np.gamma_last = rp;
+                const prob_t u = pow(np.gamma_last * pow(be(np.src_cov, np.trg_cov), 1.2), 0.1);
+                //cerr << "**rule=" << x << endl;
+                //cerr << "  u=" << log(u) << "  rule=" << rp << endl;
+                ss.add(u);
+                np.rules.push_back(TRulePtr(new TRule(x)));
+                z += u;
+              }
+            }
+            //cerr << "number of edges to consider: " << ss.size() << endl;
+            const int sampled = rng.SelectSample(ss);
+            prob_t q_n = ss[sampled] / z;
+            p = tmp_p[sampled];
+            //m.IncrementRule(*p.rules.back());
+            p.weight *= p.gamma_last / q_n;
+            //cerr << "[w=" << log(p.weight) << "]\tsampled rule: " << p.rules.back()->AsString() << endl;
+            //cerr << p << endl;
+          }
+        } // loop over particles (pi = 0 .. particles)
+        if (done_nothing) all_complete = true;
+      }
+      pfss.clear();
+      for (int i = 0; i < lps.size(); ++i)
+        pfss.add(lps[i].weight);
+      const int sampled = rng.SelectSample(pfss);
+      ps[ci] = lps[sampled];
+      m.IncrementRules(lps[sampled].rules);
+      for (int i = 0; i < lps[sampled].rules.size(); ++i) { cerr << "S:\t" << lps[sampled].rules[i]->AsString() << "\n"; }
+      cerr << "tmp-LLH: " << log(m.Likelihood()) << endl;
+    }
+    cerr << "LLH: " << log(m.Likelihood()) << endl;
+  }
+  return 0;
+}
+
diff --git a/gi/pf/reachability.cc b/gi/pf/reachability.cc
new file mode 100644
index 00000000..73dd8d39
--- /dev/null
+++ b/gi/pf/reachability.cc
@@ -0,0 +1,64 @@
+#include "reachability.h"
+
+#include <vector>
+#include <iostream>
+
+using namespace std;
+
+struct SState {
+  SState() : prev_src_covered(), prev_trg_covered() {}
+  SState(int i, int j) : prev_src_covered(i), prev_trg_covered(j) {}
+  int prev_src_covered;
+  int prev_trg_covered;
+};
+
+void Reachability::ComputeReachability(int srclen, int trglen, int src_max_phrase_len, int trg_max_phrase_len) {
+    typedef boost::multi_array<vector<SState>, 2> array_type;
+    array_type a(boost::extents[srclen + 1][trglen + 1]);
+    a[0][0].push_back(SState());
+    for (int i = 0; i < srclen; ++i) {
+      for (int j = 0; j < trglen; ++j) {
+        if (a[i][j].size() == 0) continue;
+        const SState prev(i,j);
+        for (int k = 1; k <= src_max_phrase_len; ++k) {
+          if ((i + k) > srclen) continue;
+          for (int l = 1; l <= trg_max_phrase_len; ++l) {
+            if ((j + l) > trglen) continue;
+            a[i + k][j + l].push_back(prev);
+          }
+        }
+      }
+    }
+    a[0][0].clear();
+    //cerr << "Final cell contains " << a[srclen][trglen].size() << " back pointers\n";
+    if (a[srclen][trglen].size() == 0) {
+      cerr << "Sentence with length (" << srclen << ',' << trglen << ") violates reachability constraints\n";
+      return;
+    }
+
+    typedef boost::multi_array<bool, 2> rarray_type;
+    rarray_type r(boost::extents[srclen + 1][trglen + 1]);
+    r[srclen][trglen] = true;
+    for (int i = srclen; i >= 0; --i) {
+      for (int j = trglen; j >= 0; --j) {
+        vector<SState>& prevs = a[i][j];
+        if (!r[i][j]) { prevs.clear(); }
+        for (int k = 0; k < prevs.size(); ++k) {
+          r[prevs[k].prev_src_covered][prevs[k].prev_trg_covered] = true;
+          int src_delta = i - prevs[k].prev_src_covered;
+          edges[prevs[k].prev_src_covered][prevs[k].prev_trg_covered][src_delta][j - prevs[k].prev_trg_covered] = true;
+          short &msd = max_src_delta[prevs[k].prev_src_covered][prevs[k].prev_trg_covered];
+          if (src_delta > msd) msd = src_delta;
+        }
+      }
+    }
+    assert(!edges[0][0][1][0]);
+    assert(!edges[0][0][0][1]);
+    assert(!edges[0][0][0][0]);
+    assert(max_src_delta[0][0] > 0);
+    //cerr << "First cell contains " << b[0][0].size() << " forward pointers\n";
+    //for (int i = 0; i < b[0][0].size(); ++i) {
+    //  cerr << "  -> (" << b[0][0][i].next_src_covered << "," << b[0][0][i].next_trg_covered << ")\n";
+    //}
+  }
+
diff --git a/gi/pf/reachability.h b/gi/pf/reachability.h
new file mode 100644
index 00000000..98450ec1
--- /dev/null
+++ b/gi/pf/reachability.h
@@ -0,0 +1,28 @@
+#ifndef _REACHABILITY_H_
+#define _REACHABILITY_H_
+
+#include "boost/multi_array.hpp"
+
+// determines minimum and maximum lengths of outgoing edges from all
+// coverage positions such that the alignment path respects src and
+// trg maximum phrase sizes
+//
+// runs in O(n^2 * src_max * trg_max) time but should be relatively fast
+//
+// currently forbids 0 -> n and n -> 0 alignments
+
+struct Reachability {
+  boost::multi_array<bool, 4> edges;  // edges[src_covered][trg_covered][x][trg_delta] is this edge worth exploring?
+  boost::multi_array<short, 2> max_src_delta; // msd[src_covered][trg_covered] -- the largest src delta that's valid
+
+  Reachability(int srclen, int trglen, int src_max_phrase_len, int trg_max_phrase_len) :
+      edges(boost::extents[srclen][trglen][src_max_phrase_len+1][trg_max_phrase_len+1]),
+      max_src_delta(boost::extents[srclen][trglen]) {
+    ComputeReachability(srclen, trglen, src_max_phrase_len, trg_max_phrase_len);
+  }
+
+ private:
+  void ComputeReachability(int srclen, int trglen, int src_max_phrase_len, int trg_max_phrase_len);
+};
+
+#endif
diff --git a/gi/pf/tpf.cc b/gi/pf/tpf.cc
new file mode 100644
index 00000000..7348d21c
--- /dev/null
+++ b/gi/pf/tpf.cc
@@ -0,0 +1,99 @@
+#include <iostream>
+#include <tr1/memory>
+#include <queue>
+
+#include "sampler.h"
+
+using namespace std;
+using namespace tr1;
+
+shared_ptr<MT19937> prng;
+
+struct Particle {
+  Particle() : weight(prob_t::One()) {}
+  vector<int> states;
+  prob_t weight;
+  prob_t gamma_last;
+};
+
+ostream& operator<<(ostream& os, const Particle& p) {
+  os << "[";
+  for (int i = 0; i < p.states.size(); ++i) os << p.states[i] << ' ';
+  os << "| w=" << log(p.weight) << ']';
+  return os;
+}
+
+void Rejuvenate(vector<Particle>& pps) {
+  SampleSet<prob_t> ss;
+  vector<Particle> nps(pps.size());
+  for (int i = 0; i < pps.size(); ++i) {
+//    cerr << pps[i] << endl;
+    ss.add(pps[i].weight);
+  }
+//  cerr << "REJUVINATING...\n";
+  for (int i = 0; i < pps.size(); ++i) {
+    nps[i] = pps[prng->SelectSample(ss)];
+    nps[i].weight = prob_t(1.0 / pps.size());
+//    cerr << nps[i] << endl;
+  }
+  nps.swap(pps);
+//  exit(1);
+}
+
+int main(int argc, char** argv) {
+  const unsigned particles = 100;
+  prng.reset(new MT19937);
+  MT19937& rng = *prng;
+
+  // q(a) = 0.8
+  // q(b) = 0.8
+  // q(c) = 0.4
+  SampleSet<double> ssq;
+  ssq.add(0.4);
+  ssq.add(0.6);
+  ssq.add(0);
+  double qz = 1;
+
+  // p(a) = 0.2
+  // p(b) = 0.8
+  vector<double> p(3);
+  p[0] = 0.2;
+  p[1] = 0.8;
+  p[2] = 0;
+
+  vector<int> counts(3);
+  int tot = 0;
+
+  vector<Particle> pps(particles);
+  SampleSet<prob_t> ppss;
+  int LEN = 12;
+  int PP = 1;
+  while (pps[0].states.size() < LEN) {
+    for (int pi = 0; pi < particles; ++pi) {
+      Particle& prt = pps[pi];
+
+      bool redo = true;
+      const Particle savedp = prt;
+      while (redo) {
+        redo = false;
+        for (int i = 0; i < PP; ++i) {
+          int s = rng.SelectSample(ssq);
+          double gamma_last = p[s];
+          if (!gamma_last) { redo = true; break; }
+          double q = ssq[s] / qz;
+          prt.states.push_back(s);
+          prt.weight *= prob_t(gamma_last / q);
+        }
+        if (redo) { prt = savedp; continue; }
+      }
+    }
+    Rejuvenate(pps);
+  }
+  ppss.clear();
+  for (int i = 0; i < particles; ++i) { ppss.add(pps[i].weight); }
+  int sp = rng.SelectSample(ppss);
+  cerr << pps[sp] << endl;
+
+  return 0;
+}
+
diff --git a/m4/acx_pthread.m4 b/m4/acx_pthread.m4
new file mode 100644
index 00000000..2cf20de1
--- /dev/null
+++ b/m4/acx_pthread.m4
@@ -0,0 +1,363 @@
+# This was retrieved from
+#    http://svn.0pointer.de/viewvc/trunk/common/acx_pthread.m4?revision=1277&root=avahi
+# See also (perhaps for new versions?)
+#    http://svn.0pointer.de/viewvc/trunk/common/acx_pthread.m4?root=avahi
+#
+# We've rewritten the inconsistency check code (from avahi), to work
+# more broadly.  In particular, it no longer assumes ld accepts -zdefs.
+# This caused a restructing of the code, but the functionality has only
+# changed a little.
+
+dnl @synopsis ACX_PTHREAD([ACTION-IF-FOUND[, ACTION-IF-NOT-FOUND]])
+dnl
+dnl @summary figure out how to build C programs using POSIX threads
+dnl
+dnl This macro figures out how to build C programs using POSIX threads.
+dnl It sets the PTHREAD_LIBS output variable to the threads library and
+dnl linker flags, and the PTHREAD_CFLAGS output variable to any special
+dnl C compiler flags that are needed. (The user can also force certain
+dnl compiler flags/libs to be tested by setting these environment
+dnl variables.)
+dnl
+dnl Also sets PTHREAD_CC to any special C compiler that is needed for
+dnl multi-threaded programs (defaults to the value of CC otherwise).
+dnl (This is necessary on AIX to use the special cc_r compiler alias.)
+dnl
+dnl NOTE: You are assumed to not only compile your program with these
+dnl flags, but also link it with them as well. e.g. you should link
+dnl with $PTHREAD_CC $CFLAGS $PTHREAD_CFLAGS $LDFLAGS ... $PTHREAD_LIBS
+dnl $LIBS
+dnl
+dnl If you are only building threads programs, you may wish to use
+dnl these variables in your default LIBS, CFLAGS, and CC:
+dnl
+dnl        LIBS="$PTHREAD_LIBS $LIBS"
+dnl        CFLAGS="$CFLAGS $PTHREAD_CFLAGS"
+dnl        CC="$PTHREAD_CC"
+dnl
+dnl In addition, if the PTHREAD_CREATE_JOINABLE thread-attribute
+dnl constant has a nonstandard name, defines PTHREAD_CREATE_JOINABLE to
+dnl that name (e.g. PTHREAD_CREATE_UNDETACHED on AIX).
+dnl
+dnl ACTION-IF-FOUND is a list of shell commands to run if a threads
+dnl library is found, and ACTION-IF-NOT-FOUND is a list of commands to
+dnl run it if it is not found. If ACTION-IF-FOUND is not specified, the
+dnl default action will define HAVE_PTHREAD.
+dnl
+dnl Please let the authors know if this macro fails on any platform, or
+dnl if you have any other suggestions or comments. This macro was based
+dnl on work by SGJ on autoconf scripts for FFTW (www.fftw.org) (with
+dnl help from M. Frigo), as well as ac_pthread and hb_pthread macros
+dnl posted by Alejandro Forero Cuervo to the autoconf macro repository.
+dnl We are also grateful for the helpful feedback of numerous users.
+dnl
+dnl @category InstalledPackages
+dnl @author Steven G. Johnson <stevenj@alum.mit.edu>
+dnl @version 2006-05-29
+dnl @license GPLWithACException
+dnl 
+dnl Checks for GCC shared/pthread inconsistency based on work by
+dnl Marcin Owsiany <marcin@owsiany.pl>
+
+
+AC_DEFUN([ACX_PTHREAD], [
+AC_REQUIRE([AC_CANONICAL_HOST])
+AC_LANG_SAVE
+AC_LANG_C
+acx_pthread_ok=no
+
+# We used to check for pthread.h first, but this fails if pthread.h
+# requires special compiler flags (e.g. on True64 or Sequent).
+# It gets checked for in the link test anyway.
+
+# First of all, check if the user has set any of the PTHREAD_LIBS,
+# etcetera environment variables, and if threads linking works using
+# them:
+if test x"$PTHREAD_LIBS$PTHREAD_CFLAGS" != x; then
+        save_CFLAGS="$CFLAGS"
+        CFLAGS="$CFLAGS $PTHREAD_CFLAGS"
+        save_LIBS="$LIBS"
+        LIBS="$PTHREAD_LIBS $LIBS"
+        AC_MSG_CHECKING([for pthread_join in LIBS=$PTHREAD_LIBS with CFLAGS=$PTHREAD_CFLAGS])
+        AC_TRY_LINK_FUNC(pthread_join, acx_pthread_ok=yes)
+        AC_MSG_RESULT($acx_pthread_ok)
+        if test x"$acx_pthread_ok" = xno; then
+                PTHREAD_LIBS=""
+                PTHREAD_CFLAGS=""
+        fi
+        LIBS="$save_LIBS"
+        CFLAGS="$save_CFLAGS"
+fi
+
+# We must check for the threads library under a number of different
+# names; the ordering is very important because some systems
+# (e.g. DEC) have both -lpthread and -lpthreads, where one of the
+# libraries is broken (non-POSIX).
+
+# Create a list of thread flags to try.  Items starting with a "-" are
+# C compiler flags, and other items are library names, except for "none"
+# which indicates that we try without any flags at all, and "pthread-config"
+# which is a program returning the flags for the Pth emulation library.
+
+acx_pthread_flags="pthreads none -Kthread -kthread lthread -pthread -pthreads -mthreads pthread --thread-safe -mt pthread-config"
+
+# The ordering *is* (sometimes) important.  Some notes on the
+# individual items follow:
+
+# pthreads: AIX (must check this before -lpthread)
+# none: in case threads are in libc; should be tried before -Kthread and
+#       other compiler flags to prevent continual compiler warnings
+# -Kthread: Sequent (threads in libc, but -Kthread needed for pthread.h)
+# -kthread: FreeBSD kernel threads (preferred to -pthread since SMP-able)
+# lthread: LinuxThreads port on FreeBSD (also preferred to -pthread)
+# -pthread: Linux/gcc (kernel threads), BSD/gcc (userland threads)
+# -pthreads: Solaris/gcc
+# -mthreads: Mingw32/gcc, Lynx/gcc
+# -mt: Sun Workshop C (may only link SunOS threads [-lthread], but it
+#      doesn't hurt to check since this sometimes defines pthreads too;
+#      also defines -D_REENTRANT)
+#      ... -mt is also the pthreads flag for HP/aCC
+# pthread: Linux, etcetera
+# --thread-safe: KAI C++
+# pthread-config: use pthread-config program (for GNU Pth library)
+
+case "${host_cpu}-${host_os}" in
+        *solaris*)
+
+        # On Solaris (at least, for some versions), libc contains stubbed
+        # (non-functional) versions of the pthreads routines, so link-based
+        # tests will erroneously succeed.  (We need to link with -pthreads/-mt/
+        # -lpthread.)  (The stubs are missing pthread_cleanup_push, or rather
+        # a function called by this macro, so we could check for that, but
+        # who knows whether they'll stub that too in a future libc.)  So,
+        # we'll just look for -pthreads and -lpthread first:
+
+        acx_pthread_flags="-pthreads pthread -mt -pthread $acx_pthread_flags"
+        ;;
+esac
+
+if test x"$acx_pthread_ok" = xno; then
+for flag in $acx_pthread_flags; do
+
+        case $flag in
+                none)
+                AC_MSG_CHECKING([whether pthreads work without any flags])
+                ;;
+
+                -*)
+                AC_MSG_CHECKING([whether pthreads work with $flag])
+                PTHREAD_CFLAGS="$flag"
+                ;;
+
+		pthread-config)
+		AC_CHECK_PROG(acx_pthread_config, pthread-config, yes, no)
+		if test x"$acx_pthread_config" = xno; then continue; fi
+		PTHREAD_CFLAGS="`pthread-config --cflags`"
+		PTHREAD_LIBS="`pthread-config --ldflags` `pthread-config --libs`"
+		;;
+
+                *)
+                AC_MSG_CHECKING([for the pthreads library -l$flag])
+                PTHREAD_LIBS="-l$flag"
+                ;;
+        esac
+
+        save_LIBS="$LIBS"
+        save_CFLAGS="$CFLAGS"
+        LIBS="$PTHREAD_LIBS $LIBS"
+        CFLAGS="$CFLAGS $PTHREAD_CFLAGS"
+
+        # Check for various functions.  We must include pthread.h,
+        # since some functions may be macros.  (On the Sequent, we
+        # need a special flag -Kthread to make this header compile.)
+        # We check for pthread_join because it is in -lpthread on IRIX
+        # while pthread_create is in libc.  We check for pthread_attr_init
+        # due to DEC craziness with -lpthreads.  We check for
+        # pthread_cleanup_push because it is one of the few pthread
+        # functions on Solaris that doesn't have a non-functional libc stub.
+        # We try pthread_create on general principles.
+        AC_TRY_LINK([#include <pthread.h>],
+                    [pthread_t th; pthread_join(th, 0);
+                     pthread_attr_init(0); pthread_cleanup_push(0, 0);
+                     pthread_create(0,0,0,0); pthread_cleanup_pop(0); ],
+                    [acx_pthread_ok=yes])
+
+        LIBS="$save_LIBS"
+        CFLAGS="$save_CFLAGS"
+
+        AC_MSG_RESULT($acx_pthread_ok)
+        if test "x$acx_pthread_ok" = xyes; then
+                break;
+        fi
+
+        PTHREAD_LIBS=""
+        PTHREAD_CFLAGS=""
+done
+fi
+
+# Various other checks:
+if test "x$acx_pthread_ok" = xyes; then
+        save_LIBS="$LIBS"
+        LIBS="$PTHREAD_LIBS $LIBS"
+        save_CFLAGS="$CFLAGS"
+        CFLAGS="$CFLAGS $PTHREAD_CFLAGS"
+
+        # Detect AIX lossage: JOINABLE attribute is called UNDETACHED.
+	AC_MSG_CHECKING([for joinable pthread attribute])
+	attr_name=unknown
+	for attr in PTHREAD_CREATE_JOINABLE PTHREAD_CREATE_UNDETACHED; do
+	    AC_TRY_LINK([#include <pthread.h>], [int attr=$attr; return attr;],
+                        [attr_name=$attr; break])
+	done
+        AC_MSG_RESULT($attr_name)
+        if test "$attr_name" != PTHREAD_CREATE_JOINABLE; then
+            AC_DEFINE_UNQUOTED(PTHREAD_CREATE_JOINABLE, $attr_name,
+                               [Define to necessary symbol if this constant
+                                uses a non-standard name on your system.])
+        fi
+
+        AC_MSG_CHECKING([if more special flags are required for pthreads])
+        flag=no
+        case "${host_cpu}-${host_os}" in
+            *-aix* | *-freebsd* | *-darwin*) flag="-D_THREAD_SAFE";;
+            *solaris* | *-osf* | *-hpux*) flag="-D_REENTRANT";;
+        esac
+        AC_MSG_RESULT(${flag})
+        if test "x$flag" != xno; then
+            PTHREAD_CFLAGS="$flag $PTHREAD_CFLAGS"
+        fi
+
+        LIBS="$save_LIBS"
+        CFLAGS="$save_CFLAGS"
+        # More AIX lossage: must compile with xlc_r or cc_r
+	if test x"$GCC" != xyes; then
+          AC_CHECK_PROGS(PTHREAD_CC, xlc_r cc_r, ${CC})
+        else
+          PTHREAD_CC=$CC
+	fi
+
+	# The next part tries to detect GCC inconsistency with -shared on some
+	# architectures and systems. The problem is that in certain
+	# configurations, when -shared is specified, GCC "forgets" to
+	# internally use various flags which are still necessary.
+	
+	#
+	# Prepare the flags
+	#
+	save_CFLAGS="$CFLAGS"
+	save_LIBS="$LIBS"
+	save_CC="$CC"
+	
+	# Try with the flags determined by the earlier checks.
+	#
+	# -Wl,-z,defs forces link-time symbol resolution, so that the
+	# linking checks with -shared actually have any value
+	#
+	# FIXME: -fPIC is required for -shared on many architectures,
+	# so we specify it here, but the right way would probably be to
+	# properly detect whether it is actually required.
+	CFLAGS="-shared -fPIC -Wl,-z,defs $CFLAGS $PTHREAD_CFLAGS"
+	LIBS="$PTHREAD_LIBS $LIBS"
+	CC="$PTHREAD_CC"
+	
+	# In order not to create several levels of indentation, we test
+	# the value of "$done" until we find the cure or run out of ideas.
+	done="no"
+	
+	# First, make sure the CFLAGS we added are actually accepted by our
+	# compiler.  If not (and OS X's ld, for instance, does not accept -z),
+	# then we can't do this test.
+	if test x"$done" = xno; then
+	   AC_MSG_CHECKING([whether to check for GCC pthread/shared inconsistencies])
+	   AC_TRY_LINK(,, , [done=yes])
+	
+	   if test "x$done" = xyes ; then
+	      AC_MSG_RESULT([no])
+	   else
+	      AC_MSG_RESULT([yes])
+	   fi
+	fi
+	
+	if test x"$done" = xno; then
+	   AC_MSG_CHECKING([whether -pthread is sufficient with -shared])
+	   AC_TRY_LINK([#include <pthread.h>],
+	      [pthread_t th; pthread_join(th, 0);
+	      pthread_attr_init(0); pthread_cleanup_push(0, 0);
+	      pthread_create(0,0,0,0); pthread_cleanup_pop(0); ],
+	      [done=yes])
+	   
+	   if test "x$done" = xyes; then
+	      AC_MSG_RESULT([yes])
+	   else
+	      AC_MSG_RESULT([no])
+	   fi
+	fi
+	
+	#
+	# Linux gcc on some architectures such as mips/mipsel forgets
+	# about -lpthread
+	#
+	if test x"$done" = xno; then
+	   AC_MSG_CHECKING([whether -lpthread fixes that])
+	   LIBS="-lpthread $PTHREAD_LIBS $save_LIBS"
+	   AC_TRY_LINK([#include <pthread.h>],
+	      [pthread_t th; pthread_join(th, 0);
+	      pthread_attr_init(0); pthread_cleanup_push(0, 0);
+	      pthread_create(0,0,0,0); pthread_cleanup_pop(0); ],
+	      [done=yes])
+	
+	   if test "x$done" = xyes; then
+	      AC_MSG_RESULT([yes])
+	      PTHREAD_LIBS="-lpthread $PTHREAD_LIBS"
+	   else
+	      AC_MSG_RESULT([no])
+	   fi
+	fi
+	#
+	# FreeBSD 4.10 gcc forgets to use -lc_r instead of -lc
+	#
+	if test x"$done" = xno; then
+	   AC_MSG_CHECKING([whether -lc_r fixes that])
+	   LIBS="-lc_r $PTHREAD_LIBS $save_LIBS"
+	   AC_TRY_LINK([#include <pthread.h>],
+	       [pthread_t th; pthread_join(th, 0);
+	        pthread_attr_init(0); pthread_cleanup_push(0, 0);
+	        pthread_create(0,0,0,0); pthread_cleanup_pop(0); ],
+	       [done=yes])
+	
+	   if test "x$done" = xyes; then
+	      AC_MSG_RESULT([yes])
+	      PTHREAD_LIBS="-lc_r $PTHREAD_LIBS"
+	   else
+	      AC_MSG_RESULT([no])
+	   fi
+	fi
+	if test x"$done" = xno; then
+	   # OK, we have run out of ideas
+	   AC_MSG_WARN([Impossible to determine how to use pthreads with shared libraries])
+	
+	   # so it's not safe to assume that we may use pthreads
+	   acx_pthread_ok=no
+	fi
+	
+	CFLAGS="$save_CFLAGS"
+	LIBS="$save_LIBS"
+	CC="$save_CC"
+else
+        PTHREAD_CC="$CC"
+fi
+
+AC_SUBST(PTHREAD_LIBS)
+AC_SUBST(PTHREAD_CFLAGS)
+AC_SUBST(PTHREAD_CC)
+
+# Finally, execute ACTION-IF-FOUND/ACTION-IF-NOT-FOUND:
+if test x"$acx_pthread_ok" = xyes; then
+        ifelse([$1],,AC_DEFINE(HAVE_PTHREAD,1,[Define if you have POSIX threads libraries and header files.]),[$1])
+        :
+else
+        acx_pthread_ok=no
+        $2
+fi
+AC_LANG_RESTORE
+])dnl ACX_PTHREAD
diff --git a/utils/ccrp_nt.h b/utils/ccrp_nt.h
new file mode 100644
index 00000000..63b6f4c2
--- /dev/null
+++ b/utils/ccrp_nt.h
@@ -0,0 +1,169 @@
+#ifndef _CCRP_NT_H_
+#define _CCRP_NT_H_
+
+#include <numeric>
+#include <cassert>
+#include <cmath>
+#include <list>
+#include <iostream>
+#include <vector>
+#include <tr1/unordered_map>
+#include <boost/functional/hash.hpp>
+#include "sampler.h"
+#include "slice_sampler.h"
+
+// Chinese restaurant process (1 parameter)
+template <typename Dish, typename DishHash = boost::hash<Dish> >
+class CCRP_NoTable {
+ public:
+  explicit CCRP_NoTable(double conc) :
+    num_customers_(),
+    concentration_(conc),
+    concentration_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
+    concentration_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {}
+
+  CCRP_NoTable(double c_shape, double c_rate, double c = 10.0) :
+    num_customers_(),
+    concentration_(c),
+    concentration_prior_shape_(c_shape),
+    concentration_prior_rate_(c_rate) {}
+
+  double concentration() const { return concentration_; }
+
+  bool has_concentration_prior() const {
+    return !std::isnan(concentration_prior_shape_);
+  }
+
+  void clear() {
+    num_customers_ = 0;
+    custs_.clear();
+  }
+
+  unsigned num_customers() const {
+    return num_customers_;
+  }
+
+  unsigned num_customers(const Dish& dish) const {
+    const typename std::tr1::unordered_map<Dish, unsigned, DishHash>::const_iterator it = custs_.find(dish);
+    if (it == custs_.end()) return 0;
+    return it->second;
+  }
+
+  int increment(const Dish& dish) {
+    int table_diff = 0;
+    if (++custs_[dish] == 1)
+      table_diff = 1;
+    ++num_customers_;
+    return table_diff;
+  }
+
+  int decrement(const Dish& dish) {
+    int table_diff = 0;
+    int nc = --custs_[dish];
+    if (nc == 0) {
+      custs_.erase(dish);
+      table_diff = -1;
+    } else if (nc < 0) {
+      std::cerr << "Dish counts dropped below zero for: " << dish << std::endl;
+      abort();
+    }
+    --num_customers_;
+    return table_diff;
+  }
+
+  double prob(const Dish& dish, const double& p0) const {
+    const unsigned at_table = num_customers(dish);
+    return (at_table + p0 * concentration_) / (num_customers_ + concentration_);
+  }
+
+  double logprob(const Dish& dish, const double& logp0) const {
+    const unsigned at_table = num_customers(dish);
+    return log(at_table + exp(logp0 + log(concentration_))) - log(num_customers_ + concentration_);
+  }
+
+  double log_crp_prob() const {
+    return log_crp_prob(concentration_);
+  }
+
+  static double log_gamma_density(const double& x, const double& shape, const double& rate) {
+    assert(x >= 0.0);
+    assert(shape > 0.0);
+    assert(rate > 0.0);
+    const double lp = (shape-1)*log(x) - shape*log(rate) - x/rate - lgamma(shape);
+    return lp;
+  }
+
+  // taken from http://en.wikipedia.org/wiki/Chinese_restaurant_process
+  // does not include P_0's
+  double log_crp_prob(const double& concentration) const {
+    double lp = 0.0;
+    if (has_concentration_prior())
+      lp += log_gamma_density(concentration, concentration_prior_shape_, concentration_prior_rate_);
+    assert(lp <= 0.0);
+    if (num_customers_) {
+      lp += lgamma(concentration) - lgamma(concentration + num_customers_) +
+        custs_.size() * log(concentration);
+      assert(std::isfinite(lp));
+      for (typename std::tr1::unordered_map<Dish, unsigned, DishHash>::const_iterator it = custs_.begin();
+             it != custs_.end(); ++it) {
+          lp += lgamma(it->second);
+      }
+    }
+    assert(std::isfinite(lp));
+    return lp;
+  }
+
+  void resample_hyperparameters(MT19937* rng, const unsigned nloop = 5, const unsigned niterations = 10) {
+    assert(has_concentration_prior());
+    ConcentrationResampler cr(*this);
+    for (int iter = 0; iter < nloop; ++iter) {
+        concentration_ = slice_sampler1d(cr, concentration_, *rng, 0.0,
+                               std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
+    }
+  }
+
+  struct ConcentrationResampler {
+    ConcentrationResampler(const CCRP_NoTable& crp) : crp_(crp) {}
+    const CCRP_NoTable& crp_;
+    double operator()(const double& proposed_concentration) const {
+      return crp_.log_crp_prob(proposed_concentration);
+    }
+  };
+
+  void Print(std::ostream* out) const {
+    (*out) << "DP(alpha=" << concentration_ << ") customers=" << num_customers_ << std::endl;
+    int cc = 0;
+    for (typename std::tr1::unordered_map<Dish, unsigned, DishHash>::const_iterator it = custs_.begin();
+         it != custs_.end(); ++it) {
+      (*out) << " " << it->first << "(" << it->second << " eating)";
+      ++cc;
+      if (cc > 10) { (*out) << " ..."; break; }
+    }
+    (*out) << std::endl;
+  }
+
+  unsigned num_customers_;
+  std::tr1::unordered_map<Dish, unsigned, DishHash> custs_;
+
+  typedef typename std::tr1::unordered_map<Dish, unsigned, DishHash>::const_iterator const_iterator;
+  const_iterator begin() const {
+    return custs_.begin();
+  }
+  const_iterator end() const {
+    return custs_.end();
+  }
+
+  double concentration_;
+
+  // optional gamma prior on concentration_ (NaN if no prior)
+  double concentration_prior_shape_;
+  double concentration_prior_rate_;
+};
+
+template <typename T,typename H>
+std::ostream& operator<<(std::ostream& o, const CCRP_NoTable<T,H>& c) {
+  c.Print(&o);
+  return o;
+}
+
+#endif
diff --git a/utils/ccrp_onetable.h b/utils/ccrp_onetable.h
new file mode 100644
index 00000000..a868af9a
--- /dev/null
+++ b/utils/ccrp_onetable.h
@@ -0,0 +1,241 @@
+#ifndef _CCRP_ONETABLE_H_
+#define _CCRP_ONETABLE_H_
+
+#include <numeric>
+#include <cassert>
+#include <cmath>
+#include <list>
+#include <iostream>
+#include <tr1/unordered_map>
+#include <boost/functional/hash.hpp>
+#include "sampler.h"
+#include "slice_sampler.h"
+
+// Chinese restaurant process (Pitman-Yor parameters) with one table approximation
+
+template <typename Dish, typename DishHash = boost::hash<Dish> >
+class CCRP_OneTable {
+  typedef std::tr1::unordered_map<Dish, unsigned, DishHash> DishMapType;
+ public:
+  CCRP_OneTable(double disc, double conc) :
+    num_tables_(),
+    num_customers_(),
+    discount_(disc),
+    concentration_(conc),
+    discount_prior_alpha_(std::numeric_limits<double>::quiet_NaN()),
+    discount_prior_beta_(std::numeric_limits<double>::quiet_NaN()),
+    concentration_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
+    concentration_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {}
+
+  CCRP_OneTable(double d_alpha, double d_beta, double c_shape, double c_rate, double d = 0.9, double c = 1.0) :
+    num_tables_(),
+    num_customers_(),
+    discount_(d),
+    concentration_(c),
+    discount_prior_alpha_(d_alpha),
+    discount_prior_beta_(d_beta),
+    concentration_prior_shape_(c_shape),
+    concentration_prior_rate_(c_rate) {}
+
+  double discount() const { return discount_; }
+  double concentration() const { return concentration_; }
+  void set_concentration(double c) { concentration_ = c; }
+  void set_discount(double d) { discount_ = d; }
+
+  bool has_discount_prior() const {
+    return !std::isnan(discount_prior_alpha_);
+  }
+
+  bool has_concentration_prior() const {
+    return !std::isnan(concentration_prior_shape_);
+  }
+
+  void clear() {
+    num_tables_ = 0;
+    num_customers_ = 0;
+    dish_counts_.clear();
+  }
+
+  unsigned num_tables() const {
+    return num_tables_;
+  }
+
+  unsigned num_tables(const Dish& dish) const {
+    const typename DishMapType::const_iterator it = dish_counts_.find(dish);
+    if (it == dish_counts_.end()) return 0;
+    return 1;
+  }
+
+  unsigned num_customers() const {
+    return num_customers_;
+  }
+
+  unsigned num_customers(const Dish& dish) const {
+    const typename DishMapType::const_iterator it = dish_counts_.find(dish);
+    if (it == dish_counts_.end()) return 0;
+    return it->second;
+  }
+
+  // returns +1 or 0 indicating whether a new table was opened
+  int increment(const Dish& dish) {
+    unsigned& dc = dish_counts_[dish];
+    ++dc;
+    ++num_customers_;
+    if (dc == 1) {
+      ++num_tables_;
+      return 1;
+    } else {
+      return 0;
+    }
+  }
+
+  // returns -1 or 0, indicating whether a table was closed
+  int decrement(const Dish& dish) {
+    unsigned& dc = dish_counts_[dish];
+    assert(dc > 0);
+    if (dc == 1) {
+      dish_counts_.erase(dish);
+      --num_tables_;
+      --num_customers_;
+      return -1;
+    } else {
+      assert(dc > 1);
+      --dc;
+      --num_customers_;
+      return 0;
+    }
+  }
+
+  double prob(const Dish& dish, const double& p0) const {
+    const typename DishMapType::const_iterator it = dish_counts_.find(dish);
+    const double r = num_tables_ * discount_ + concentration_;
+    if (it == dish_counts_.end()) {
+      return r * p0 / (num_customers_ + concentration_);
+    } else {
+      return (it->second - discount_ + r * p0) /
+               (num_customers_ + concentration_);
+    }
+  }
+
+  double log_crp_prob() const {
+    return log_crp_prob(discount_, concentration_);
+  }
+
+  static double log_beta_density(const double& x, const double& alpha, const double& beta) {
+    assert(x > 0.0);
+    assert(x < 1.0);
+    assert(alpha > 0.0);
+    assert(beta > 0.0);
+    const double lp = (alpha-1)*log(x)+(beta-1)*log(1-x)+lgamma(alpha+beta)-lgamma(alpha)-lgamma(beta);
+    return lp;
+  }
+
+  static double log_gamma_density(const double& x, const double& shape, const double& rate) {
+    assert(x >= 0.0);
+    assert(shape > 0.0);
+    assert(rate > 0.0);
+    const double lp = (shape-1)*log(x) - shape*log(rate) - x/rate - lgamma(shape);
+    return lp;
+  }
+
+  // taken from http://en.wikipedia.org/wiki/Chinese_restaurant_process
+  // does not include P_0's
+  double log_crp_prob(const double& discount, const double& concentration) const {
+    double lp = 0.0;
+    if (has_discount_prior())
+      lp = log_beta_density(discount, discount_prior_alpha_, discount_prior_beta_);
+    if (has_concentration_prior())
+      lp += log_gamma_density(concentration, concentration_prior_shape_, concentration_prior_rate_);
+    assert(lp <= 0.0);
+    if (num_customers_) {
+      if (discount > 0.0) {
+        const double r = lgamma(1.0 - discount);
+        lp += lgamma(concentration) - lgamma(concentration + num_customers_)
+             + num_tables_ * log(discount) + lgamma(concentration / discount + num_tables_)
+             - lgamma(concentration / discount);
+        assert(std::isfinite(lp));
+        for (typename DishMapType::const_iterator it = dish_counts_.begin();
+             it != dish_counts_.end(); ++it) {
+          const unsigned& cur = it->second;
+          lp += lgamma(cur - discount) - r;
+        }
+      } else {
+        assert(!"not implemented yet");
+      }
+    }
+    assert(std::isfinite(lp));
+    return lp;
+  }
+
+  void resample_hyperparameters(MT19937* rng, const unsigned nloop = 5, const unsigned niterations = 10) {
+    assert(has_discount_prior() || has_concentration_prior());
+    DiscountResampler dr(*this);
+    ConcentrationResampler cr(*this);
+    for (int iter = 0; iter < nloop; ++iter) {
+      if (has_concentration_prior()) {
+        concentration_ = slice_sampler1d(cr, concentration_, *rng, 0.0,
+                               std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
+      }
+      if (has_discount_prior()) {
+        discount_ = slice_sampler1d(dr, discount_, *rng, std::numeric_limits<double>::min(),
+                               1.0, 0.0, niterations, 100*niterations);
+      }
+    }
+    concentration_ = slice_sampler1d(cr, concentration_, *rng, 0.0,
+                             std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
+  }
+
+  struct DiscountResampler {
+    DiscountResampler(const CCRP_OneTable& crp) : crp_(crp) {}
+    const CCRP_OneTable& crp_;
+    double operator()(const double& proposed_discount) const {
+      return crp_.log_crp_prob(proposed_discount, crp_.concentration_);
+    }
+  };
+
+  struct ConcentrationResampler {
+    ConcentrationResampler(const CCRP_OneTable& crp) : crp_(crp) {}
+    const CCRP_OneTable& crp_;
+    double operator()(const double& proposed_concentration) const {
+      return crp_.log_crp_prob(crp_.discount_, proposed_concentration);
+    }
+  };
+
+  void Print(std::ostream* out) const {
+    (*out) << "PYP(d=" << discount_ << ",c=" << concentration_ << ") customers=" << num_customers_ << std::endl;
+    for (typename DishMapType::const_iterator it = dish_counts_.begin(); it != dish_counts_.end(); ++it) {
+      (*out) << "  " << it->first << " = " << it->second << std::endl;
+    }
+  }
+
+  typedef typename DishMapType::const_iterator const_iterator;
+  const_iterator begin() const {
+    return dish_counts_.begin();
+  }
+  const_iterator end() const {
+    return dish_counts_.end();
+  }
+
+  unsigned num_tables_;
+  unsigned num_customers_;
+  DishMapType dish_counts_;
+
+  double discount_;
+  double concentration_;
+
+  // optional beta prior on discount_ (NaN if no prior)
+  double discount_prior_alpha_;
+  double discount_prior_beta_;
+
+  // optional gamma prior on concentration_ (NaN if no prior)
+  double concentration_prior_shape_;
+  double concentration_prior_rate_;
+};
+
+template <typename T,typename H>
+std::ostream& operator<<(std::ostream& o, const CCRP_OneTable<T,H>& c) {
+  c.Print(&o);
+  return o;
+}
+
+#endif