Python interface for more parsing factors

ad3/Factor.h

@@ -167,6 +167,7 @@ class Factor {
 
   void SetAdditionalLogPotentials(
       const vector<double> &additional_log_potentials) {
+    CheckAdditionalLogPotentials(additional_log_potentials);
     additional_log_potentials_ = additional_log_potentials;
   }
 
@@ -218,6 +219,13 @@ class Factor {
             &additional_posteriors_last_);
   }
 
+  // check whether the given potentials are valid and raise an exception
+  // if they are not
+  virtual const void CheckAdditionalLogPotentials(
+    const vector<double> &additional_log_potentials){
+    return;
+  }
+
  private:
   int id_; // Factor id.
 

ad3/logval.h

@@ -0,0 +1,226 @@
+#ifndef LOGVAL_H_
+#define LOGVAL_H_
+
+#define LOGVAL_CHECK_NEG 0
+
+#include <iostream>
+#include <cstdlib>
+#include <cmath>
+#include <limits>
+#include <cassert>
+
+#define LOGVAL_LOG0 -std::numeric_limits<T>::infinity()
+
+//TODO: template for supporting negation or not - most uses are for nonnegative "probs" only; probably some 10-20% speedup available
+template <class T>
+class LogVal {
+public:
+  void print(std::ostream &o) const {
+    if (s_) o << "(-)";
+    o << v_;
+  }
+  //PRINT_SELF(LogVal<T>)
+
+  typedef LogVal<T> Self;
+
+  LogVal() : s_(), v_(LOGVAL_LOG0) {}
+  LogVal(double x) : s_(std::signbit(x)), v_(s_ ? std::log(-x) : std::log(x)) {}
+  const Self& operator=(double x) { s_ = std::signbit(x); v_ = s_ ? std::log(-x) : std::log(x); return *this; }
+  LogVal(double lnx, bool sign) : s_(sign), v_(lnx) {}
+  static Self exp(T lnx) { return Self(lnx, false); }
+
+  // maybe the below are faster than == 1 and == 0.  i don't know.
+  bool is_1() const { return v_ == 0 && s_ == 0; }
+  bool is_0() const { return v_ == LOGVAL_LOG0; }
+
+  static Self One() { return Self(1); }
+  static Self Zero() { return Self(); }
+  static Self e() { return Self(1, false); }
+  void logeq(const T& v) { s_ = false; v_ = v; }
+
+  // just like std::signbit, negative means true.  weird, i know
+  bool signbit() const {
+    return s_;
+  }
+  friend inline bool signbit(Self const& x) { return x.signbit(); }
+
+  T logabs() const {
+    return v_;
+  }
+
+  Self& besteq(const Self& a) {
+    assert(!a.s_ && !s_);
+    if (a.v_ < v_)
+      v_ = a.v_;
+    return *this;
+  }
+
+  Self& operator+=(const Self& a) {
+    if (a.is_0()) return *this;
+    if (a.s_ == s_) {
+      if (a.v_ < v_) {
+        v_ = v_ + log1p(std::exp(a.v_ - v_));
+      } else {
+        v_ = a.v_ + log1p(std::exp(v_ - a.v_));
+      }
+    } else {
+      if (a.v_ < v_) {
+        v_ = v_ + log1p(-std::exp(a.v_ - v_));
+      } else {
+        v_ = a.v_ + log1p(-std::exp(v_ - a.v_));
+        s_ = !s_;
+      }
+    }
+    return *this;
+  }
+
+  Self& operator*=(const Self& a) {
+    s_ = (s_ != a.s_);
+    v_ += a.v_;
+    return *this;
+  }
+
+  Self& operator/=(const Self& a) {
+    s_ = (s_ != a.s_);
+    v_ -= a.v_;
+    return *this;
+  }
+
+  Self& operator-=(const Self& a) {
+    Self b = a;
+    b.negate();
+    return *this += b;
+  }
+
+  // Self(fabs(log(x)),x.s_)
+  friend Self abslog(Self x) {
+    if (x.v_ < 0) x.v_ = -x.v_;
+    return x;
+  }
+
+  Self& poweq(const T& power) {
+#if LOGVAL_CHECK_NEG
+    if (s_) {
+      std::cerr << "poweq(T) not implemented when s_ is true\n";
+      std::abort();
+    } else
+#endif
+      v_ *= power;
+    return *this;
+  }
+
+  //remember, s_ means negative.
+  inline bool lt(Self const& o) const {
+    return s_ == o.s_ ? v_ < o.v_ : s_ > o.s_;
+  }
+  inline bool gt(Self const& o) const {
+    return s_ == o.s_ ? o.v_ < v_ : s_ < o.s_;
+  }
+
+  Self operator-() const {
+    return Self(v_, !s_);
+  }
+  void negate() { s_ = !s_; }
+
+  Self inverse() const { return Self(-v_, s_); }
+
+  Self pow(const T& power) const {
+    Self res = *this;
+    res.poweq(power);
+    return res;
+  }
+
+  Self root(const T& root) const {
+    return pow(1 / root);
+  }
+
+  T as_float() const {
+    if (s_) return -std::exp(v_); else return std::exp(v_);
+  }
+
+  bool s_;
+  T v_;
+};
+
+// copy elision - as opposed to explicit copy of LogVal<T> const& o1, we should be able to construct Logval r=a+(b+c) as a single result in place in r.  todo: return std::move(o1) - C++0x
+template<class T>
+LogVal<T> operator+(LogVal<T> o1, const LogVal<T>& o2) {
+  o1 += o2;
+  return o1;
+}
+
+template<class T>
+LogVal<T> operator*(LogVal<T> o1, const LogVal<T>& o2) {
+  o1 *= o2;
+  return o1;
+}
+
+template<class T>
+LogVal<T> operator/(LogVal<T> o1, const LogVal<T>& o2) {
+  o1 /= o2;
+  return o1;
+}
+
+template<class T>
+LogVal<T> operator-(LogVal<T> o1, const LogVal<T>& o2) {
+  o1 -= o2;
+  return o1;
+}
+
+template<class T>
+T log(const LogVal<T>& o) {
+#ifdef LOGVAL_CHECK_NEG
+  if (o.s_) return log(-1.0);
+#endif
+  return o.v_;
+}
+
+template<class T>
+LogVal<T> abs(const LogVal<T>& o) {
+  if (o.s_) {
+    LogVal<T> res = o;
+    res.s_ = false;
+    return res;
+  } else { return o; }
+}
+
+template <class T>
+LogVal<T> pow(const LogVal<T>& b, const T& e) {
+  return b.pow(e);
+}
+
+template <class T>
+bool operator==(const LogVal<T>& lhs, const LogVal<T>& rhs) {
+  return (lhs.v_ == rhs.v_) && (lhs.s_ == rhs.s_);
+}
+
+template <class T>
+bool operator!=(const LogVal<T>& lhs, const LogVal<T>& rhs) {
+  return !(lhs == rhs);
+}
+
+template <class T>
+bool operator<(const LogVal<T>& lhs, const LogVal<T>& rhs) {
+  if (lhs.s_ == rhs.s_) {
+    return (lhs.v_ < rhs.v_);
+  } else {
+    return lhs.s_ > rhs.s_;
+  }
+}
+
+template <class T>
+bool operator<=(const LogVal<T>& lhs, const LogVal<T>& rhs) {
+  return (lhs < rhs) || (lhs == rhs);
+}
+
+template <class T>
+bool operator>(const LogVal<T>& lhs, const LogVal<T>& rhs) {
+  return !(lhs <= rhs);
+}
+
+template <class T>
+bool operator>=(const LogVal<T>& lhs, const LogVal<T>& rhs) {
+  return !(lhs < rhs);
+}
+
+#endif

ad3_multi.cpp

@@ -570,8 +570,11 @@ int LoadGraph(ifstream &file_graph,
       vector<vector<int> > index_siblings(length, vector<int>(length+1, -1));
       int total = 0;
       vector<Sibling*> siblings;
+      vector<Arc*> arcs;
       vector<double> additional_scores;
       for (int m = 0; m < length; ++m) {
+        Arc *arc = new Arc(0, m);
+        arcs.push_back(arc);
         for (int s = m+1; s <= length; ++s) {
           // Create a fake sibling.
           Sibling *sibling = new Sibling(0, m, s);
@@ -584,10 +587,13 @@ int LoadGraph(ifstream &file_graph,
       }
       factor = new FactorHeadAutomaton;
       factor_graph->DeclareFactor(factor, binary_variables, true);
-      static_cast<FactorHeadAutomaton*>(factor)->Initialize(length, siblings);
+      static_cast<FactorHeadAutomaton*>(factor)->Initialize(arcs, siblings);
       for (int r = 0; r < siblings.size(); ++r) {
         delete siblings[r];
       }
+      for (int r = 0; r < arcs.size(); ++r) {
+        delete arcs[r];
+      }
       factor->SetAdditionalLogPotentials(additional_scores);
       num_factor_log_potentials += additional_scores.size();
       cout << "Read head automaton factor." << endl;

examples/cpp/parsing/Decode.cpp

@@ -0,0 +1,110 @@
+#include <Eigen/Dense>
+#include "FactorTree.h"
+#include "Utils.h"
+#include "logval.h"
+
+// Define a matrix of doubles using Eigen.
+typedef LogVal<double> LogValD;
+namespace Eigen {
+  typedef Eigen::Matrix<LogValD, Dynamic, Dynamic> MatrixXlogd;
+}
+
+namespace AD3
+{
+  // decoder for an arc factored parser; it invokes the matrix-tree theorem.
+  // index is a n x (n - 1) matrix mapping (head, modifier) to either -1 (if
+  // it doesn't exist) or its index in the scores vector.
+  void DecodeMatrixTree(const vector<vector<int> > &index, const vector<Arc*> &arcs,
+                        const vector<double> &scores,
+                        vector<double> *predicted_output,
+                        double *log_partition_function, double *entropy) {
+    int length = index.size();
+
+    // Matrix for storing the potentials.
+    Eigen::MatrixXlogd potentials(length, length);
+    // Kirchhoff matrix.
+    Eigen::MatrixXlogd kirchhoff(length - 1, length - 1);
+
+    // Compute an offset to improve numerical stability. This is a constant that
+    // is subtracted from all scores.
+    int num_arcs = arcs.size();
+    double constant = 0.0;
+    for (int r = 0; r < num_arcs; ++r) {
+      constant += scores[r];
+    }
+    constant /= static_cast<double>(num_arcs);
+
+    // Set the potentials.
+    for (int h = 0; h < length; ++h) {
+      for (int m = 0; m < length; ++m) {
+      potentials(m, h) = LogValD::Zero();
+      int r = index[h][m];
+      if (r >= 0) {
+        potentials(m, h) = LogValD(scores[r] - constant, false);
+      }
+      }
+    }
+
+    // Set the Kirchhoff matrix.
+    for (int h = 0; h < length - 1; ++h) {
+      for (int m = 0; m < length - 1; ++m) {
+      kirchhoff(h, m) = -potentials(m + 1, h + 1);
+      }
+    }
+    for (int m = 1; m < length; ++m) {
+      LogValD sum = LogValD::Zero();
+      for (int h = 0; h < length; ++h) {
+      sum += potentials(m, h);
+      }
+      kirchhoff(m - 1, m - 1) = sum;
+    }
+
+    // Inverse of the Kirchoff matrix.
+    Eigen::FullPivLU<Eigen::MatrixXlogd> lu(kirchhoff);
+    Eigen::MatrixXlogd inverted_kirchhoff = lu.inverse();
+    *log_partition_function = lu.determinant().logabs() +
+        constant * (length - 1);
+
+    Eigen::MatrixXlogd marginals(length, length);
+    for (int h = 0; h < length; ++h) {
+      marginals(0, h) = LogValD::Zero();
+    }
+    for (int m = 1; m < length; ++m) {
+      marginals(m, 0) = potentials(m, 0) * inverted_kirchhoff(m - 1, m - 1);
+      for (int h = 1; h < length; ++h) {
+      marginals(m, h) = potentials(m, h) *
+          (inverted_kirchhoff(m - 1, m - 1) - inverted_kirchhoff(m - 1, h - 1));
+      }
+    }
+
+    // Compute the entropy.
+    predicted_output->resize(num_arcs);
+    *entropy = *log_partition_function;
+    for (int r = 0; r < num_arcs; ++r) {
+      int h = arcs[r]->head();
+      int m = arcs[r]->modifier();
+      if (marginals(m, h).signbit()) {
+      if (!NEARLY_ZERO_TOL(marginals(m, h).as_float(), 1e-6)) {
+        // LOG(INFO) << "Marginals truncated to zero (" << marginals(m, h).as_float() << ")";
+        cerr << "Marginals truncated to zero (" << marginals(m, h).as_float() << ")";
+      }
+      // CHECK(!std::isnan(marginals(m, h).as_float()));
+      } else if (marginals(m, h).logabs() > 0) {
+      if (!NEARLY_ZERO_TOL(marginals(m, h).as_float() - 1.0, 1e-6)) {
+        // LOG(INFO) << "Marginals truncated to one (" << marginals(m, h).as_float() << ")";
+        cerr << "Marginals truncated to one (" << marginals(m, h).as_float() << ")";
+      }
+      }
+      (*predicted_output)[r] = marginals(m, h).as_float();
+      *entropy -= (*predicted_output)[r] * scores[r];
+    }
+    if (*entropy < 0.0) {
+      if (!NEARLY_ZERO_TOL(*entropy, 1e-6)) {
+        // LOG(INFO) << "Entropy truncated to zero (" << *entropy << ")";
+        cerr << "Entropy truncated to zero (" << *entropy << ")";
+      }
+      *entropy = 0.0;
+    }
+  }
+
+} // AD3