support strength=0 PYPs, final notation clean-up

2 files changed, 115 insertions, 96 deletions

diff --git a/utils/ccrp.h b/utils/ccrp.h
index c883c027..5f9db7a6 100644
--- a/utils/ccrp.h
+++ b/utils/ccrp.h
@@ -18,27 +18,27 @@
 template <typename Dish, typename DishHash = boost::hash<Dish> >
 class CCRP {
  public:
-  CCRP(double disc, double alpha) :
+  CCRP(double disc, double strength) :
       num_tables_(),
       num_customers_(),
       discount_(disc),
-      alpha_(alpha),
-      discount_prior_alpha_(std::numeric_limits<double>::quiet_NaN()),
+      strength_(strength),
+      discount_prior_strength_(std::numeric_limits<double>::quiet_NaN()),
       discount_prior_beta_(std::numeric_limits<double>::quiet_NaN()),
-      alpha_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
-      alpha_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {
+      strength_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
+      strength_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {
     check_hyperparameters();
   }
 
-  CCRP(double d_alpha, double d_beta, double c_shape, double c_rate, double d = 0.9, double c = 1.0) :
+  CCRP(double d_strength, double d_beta, double c_shape, double c_rate, double d = 0.9, double c = 1.0) :
       num_tables_(),
       num_customers_(),
       discount_(d),
-      alpha_(c),
-      discount_prior_alpha_(d_alpha),
+      strength_(c),
+      discount_prior_strength_(d_strength),
       discount_prior_beta_(d_beta),
-      alpha_prior_shape_(c_shape),
-      alpha_prior_rate_(c_rate) {
+      strength_prior_shape_(c_shape),
+      strength_prior_rate_(c_rate) {
     check_hyperparameters();
   }
 
@@ -47,23 +47,23 @@ class CCRP {
       std::cerr << "Bad discount: " << discount_ << std::endl;
       abort();
     }
-    if (alpha_ <= -discount_) {
-      std::cerr << "Bad strength: " << alpha_ << " (discount=" << discount_ << ")" << std::endl;
+    if (strength_ <= -discount_) {
+      std::cerr << "Bad strength: " << strength_ << " (discount=" << discount_ << ")" << std::endl;
       abort();
     }
   }
 
   double discount() const { return discount_; }
-  double alpha() const { return alpha_; }
+  double strength() const { return strength_; }
   void set_discount(double d) { discount_ = d; check_hyperparameters(); }
-  void set_alpha(double a) { alpha_ = a; check_hyperparameters(); }
+  void set_strength(double a) { strength_ = a; check_hyperparameters(); }
 
   bool has_discount_prior() const {
-    return !std::isnan(discount_prior_alpha_);
+    return !std::isnan(discount_prior_strength_);
   }
 
-  bool has_alpha_prior() const {
-    return !std::isnan(alpha_prior_shape_);
+  bool has_strength_prior() const {
+    return !std::isnan(strength_prior_shape_);
   }
 
   void clear() {
@@ -97,7 +97,7 @@ class CCRP {
     DishLocations& loc = dish_locs_[dish];
     bool share_table = false;
     if (loc.total_dish_count_) {
-      const double p_empty = (alpha_ + num_tables_ * discount_) * p0;
+      const double p_empty = (strength_ + num_tables_ * discount_) * p0;
       const double p_share = (loc.total_dish_count_ - loc.table_counts_.size() * discount_);
       share_table = rng->SelectSample(p_empty, p_share);
     }
@@ -131,7 +131,7 @@ class CCRP {
     DishLocations& loc = dish_locs_[dish];
     bool share_table = false;
     if (loc.total_dish_count_) {
-      const T p_empty = T(alpha_ + num_tables_ * discount_) * p0;
+      const T p_empty = T(strength_ + num_tables_ * discount_) * p0;
       const T p_share = T(loc.total_dish_count_ - loc.table_counts_.size() * discount_);
       share_table = rng->SelectSample(p_empty, p_share);
     }
@@ -198,47 +198,47 @@ class CCRP {
 
   double prob(const Dish& dish, const double& p0) const {
     const typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.find(dish);
-    const double r = num_tables_ * discount_ + alpha_;
+    const double r = num_tables_ * discount_ + strength_;
     if (it == dish_locs_.end()) {
-      return r * p0 / (num_customers_ + alpha_);
+      return r * p0 / (num_customers_ + strength_);
     } else {
       return (it->second.total_dish_count_ - discount_ * it->second.table_counts_.size() + r * p0) /
-               (num_customers_ + alpha_);
+               (num_customers_ + strength_);
     }
   }
 
   template <typename T>
   T probT(const Dish& dish, const T& p0) const {
     const typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.find(dish);
-    const T r = T(num_tables_ * discount_ + alpha_);
+    const T r = T(num_tables_ * discount_ + strength_);
     if (it == dish_locs_.end()) {
-      return r * p0 / T(num_customers_ + alpha_);
+      return r * p0 / T(num_customers_ + strength_);
     } else {
       return (T(it->second.total_dish_count_ - discount_ * it->second.table_counts_.size()) + r * p0) /
-               T(num_customers_ + alpha_);
+               T(num_customers_ + strength_);
     }
   }
 
   double log_crp_prob() const {
-    return log_crp_prob(discount_, alpha_);
+    return log_crp_prob(discount_, strength_);
   }
 
   // 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& alpha) const {
+  double log_crp_prob(const double& discount, const double& strength) const {
     double lp = 0.0;
     if (has_discount_prior())
-      lp = Md::log_beta_density(discount, discount_prior_alpha_, discount_prior_beta_);
-    if (has_alpha_prior())
-      lp += Md::log_gamma_density(alpha + discount, alpha_prior_shape_, alpha_prior_rate_);
+      lp = Md::log_beta_density(discount, discount_prior_strength_, discount_prior_beta_);
+    if (has_strength_prior())
+      lp += Md::log_gamma_density(strength + discount, strength_prior_shape_, strength_prior_rate_);
     assert(lp <= 0.0);
     if (num_customers_) {
       if (discount > 0.0) {
         const double r = lgamma(1.0 - discount);
-        if (alpha)
-          lp += lgamma(alpha) - lgamma(alpha / discount);
-        lp += - lgamma(alpha + num_customers_)
-             + num_tables_ * log(discount) + lgamma(alpha / discount + num_tables_);
+        if (strength)
+          lp += lgamma(strength) - lgamma(strength / discount);
+        lp += - lgamma(strength + num_customers_)
+             + num_tables_ * log(discount) + lgamma(strength / discount + num_tables_);
         assert(std::isfinite(lp));
         for (typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.begin();
              it != dish_locs_.end(); ++it) {
@@ -247,8 +247,16 @@ class CCRP {
             lp += lgamma(*ti - discount) - r;
           }
         }
+      } else if (!discount) { // discount == 0.0
+        lp += lgamma(strength) + num_tables_ * log(strength) - lgamma(strength + num_tables_);
+        assert(std::isfinite(lp));
+        for (typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.begin();
+             it != dish_locs_.end(); ++it) {
+          const DishLocations& cur = it->second;
+          lp += lgamma(cur.table_counts_.size());
+        }
       } else {
-        assert(!"not implemented yet");
+        assert(!"discount less than 0 detected!");
       }
     }
     assert(std::isfinite(lp));
@@ -256,22 +264,22 @@ class CCRP {
   }
 
   void resample_hyperparameters(MT19937* rng, const unsigned nloop = 5, const unsigned niterations = 10) {
-    assert(has_discount_prior() || has_alpha_prior());
+    assert(has_discount_prior() || has_strength_prior());
     DiscountResampler dr(*this);
     StrengthResampler sr(*this);
     for (int iter = 0; iter < nloop; ++iter) {
-      if (has_alpha_prior()) {
-        alpha_ = slice_sampler1d(sr, alpha_, *rng, -discount_,
+      if (has_strength_prior()) {
+        strength_ = slice_sampler1d(sr, strength_, *rng, -discount_,
                                std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
       }
       if (has_discount_prior()) {
         double min_discount = std::numeric_limits<double>::min();
-        if (alpha_ < 0.0) min_discount = -alpha_;
+        if (strength_ < 0.0) min_discount = -strength_;
         discount_ = slice_sampler1d(dr, discount_, *rng, min_discount,
                                1.0, 0.0, niterations, 100*niterations);
       }
     }
-    alpha_ = slice_sampler1d(sr, alpha_, *rng, -discount_,
+    strength_ = slice_sampler1d(sr, strength_, *rng, -discount_,
                              std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
   }
 
@@ -279,15 +287,15 @@ class CCRP {
     DiscountResampler(const CCRP& crp) : crp_(crp) {}
     const CCRP& crp_;
     double operator()(const double& proposed_discount) const {
-      return crp_.log_crp_prob(proposed_discount, crp_.alpha_);
+      return crp_.log_crp_prob(proposed_discount, crp_.strength_);
     }
   };
 
   struct StrengthResampler {
     StrengthResampler(const CCRP& crp) : crp_(crp) {}
     const CCRP& crp_;
-    double operator()(const double& proposed_alpha) const {
-      return crp_.log_crp_prob(crp_.discount_, proposed_alpha);
+    double operator()(const double& proposed_strength) const {
+      return crp_.log_crp_prob(crp_.discount_, proposed_strength);
     }
   };
 
@@ -299,7 +307,7 @@ class CCRP {
   };
 
   void Print(std::ostream* out) const {
-    std::cerr << "PYP(d=" << discount_ << ",c=" << alpha_ << ") customers=" << num_customers_ << std::endl;
+    std::cerr << "PYP(d=" << discount_ << ",c=" << strength_ << ") customers=" << num_customers_ << std::endl;
     for (typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.begin();
          it != dish_locs_.end(); ++it) {
       (*out) << it->first << " (" << it->second.total_dish_count_ << " on " << it->second.table_counts_.size() << " tables): ";
@@ -324,15 +332,15 @@ class CCRP {
   std::tr1::unordered_map<Dish, DishLocations, DishHash> dish_locs_;
 
   double discount_;
-  double alpha_;
+  double strength_;
 
   // optional beta prior on discount_ (NaN if no prior)
-  double discount_prior_alpha_;
+  double discount_prior_strength_;
   double discount_prior_beta_;
 
-  // optional gamma prior on alpha_ (NaN if no prior)
-  double alpha_prior_shape_;
-  double alpha_prior_rate_;
+  // optional gamma prior on strength_ (NaN if no prior)
+  double strength_prior_shape_;
+  double strength_prior_rate_;
 };
 
 template <typename T,typename H>
diff --git a/utils/mfcr.h b/utils/mfcr.h
index df988f51..aeaf599d 100644
--- a/utils/mfcr.h
+++ b/utils/mfcr.h
@@ -39,37 +39,37 @@ template <typename Dish, typename DishHash = boost::hash<Dish> >
 class MFCR {
  public:
 
-  MFCR(unsigned num_floors, double d, double alpha) :
+  MFCR(unsigned num_floors, double d, double strength) :
     num_floors_(num_floors),
     num_tables_(),
     num_customers_(),
     discount_(d),
-    alpha_(alpha),
-    discount_prior_alpha_(std::numeric_limits<double>::quiet_NaN()),
+    strength_(strength),
+    discount_prior_strength_(std::numeric_limits<double>::quiet_NaN()),
     discount_prior_beta_(std::numeric_limits<double>::quiet_NaN()),
-    alpha_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
-    alpha_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {}
+    strength_prior_shape_(std::numeric_limits<double>::quiet_NaN()),
+    strength_prior_rate_(std::numeric_limits<double>::quiet_NaN()) {}
 
-  MFCR(unsigned num_floors, double discount_alpha, double discount_beta, double alpha_shape, double alpha_rate, double d = 0.9, double alpha = 10.0) :
+  MFCR(unsigned num_floors, double discount_strength, double discount_beta, double strength_shape, double strength_rate, double d = 0.9, double strength = 10.0) :
     num_floors_(num_floors),
     num_tables_(),
     num_customers_(),
     discount_(d),
-    alpha_(alpha),
-    discount_prior_alpha_(discount_alpha),
+    strength_(strength),
+    discount_prior_strength_(discount_strength),
     discount_prior_beta_(discount_beta),
-    alpha_prior_shape_(alpha_shape),
-    alpha_prior_rate_(alpha_rate) {}
+    strength_prior_shape_(strength_shape),
+    strength_prior_rate_(strength_rate) {}
 
   double discount() const { return discount_; }
-  double alpha() const { return alpha_; }
+  double strength() const { return strength_; }
 
   bool has_discount_prior() const {
-    return !std::isnan(discount_prior_alpha_);
+    return !std::isnan(discount_prior_strength_);
   }
 
-  bool has_alpha_prior() const {
-    return !std::isnan(alpha_prior_shape_);
+  bool has_strength_prior() const {
+    return !std::isnan(strength_prior_shape_);
   }
 
   void clear() {
@@ -122,7 +122,7 @@ class MFCR {
     int floor = -1;
     bool share_table = false;
     if (loc.total_dish_count_) {
-      const double p_empty = (alpha_ + num_tables_ * discount_) * marg_p0;
+      const double p_empty = (strength_ + num_tables_ * discount_) * marg_p0;
       const double p_share = (loc.total_dish_count_ - loc.table_counts_.size() * discount_);
       share_table = rng->SelectSample(p_empty, p_share);
     }
@@ -206,44 +206,53 @@ class MFCR {
     const double marg_p0 = std::inner_product(p0s.begin(), p0s.end(), lambdas.begin(), 0.0);
     assert(marg_p0 <= 1.0);
     const typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.find(dish);
-    const double r = num_tables_ * discount_ + alpha_;
+    const double r = num_tables_ * discount_ + strength_;
     if (it == dish_locs_.end()) {
-      return r * marg_p0 / (num_customers_ + alpha_);
+      return r * marg_p0 / (num_customers_ + strength_);
     } else {
       return (it->second.total_dish_count_ - discount_ * it->second.table_counts_.size() + r * marg_p0) /
-               (num_customers_ + alpha_);
+               (num_customers_ + strength_);
     }
   }
 
   double log_crp_prob() const {
-    return log_crp_prob(discount_, alpha_);
+    return log_crp_prob(discount_, strength_);
   }
 
   // taken from http://en.wikipedia.org/wiki/Chinese_restaurant_process
   // does not include draws from G_w's
-  double log_crp_prob(const double& d, const double& alpha) const {
+  double log_crp_prob(const double& discount, const double& strength) const {
     double lp = 0.0;
     if (has_discount_prior())
-      lp = Md::log_beta_density(d, discount_prior_alpha_, discount_prior_beta_);
-    if (has_alpha_prior())
-      lp += Md::log_gamma_density(alpha, alpha_prior_shape_, alpha_prior_rate_);
+      lp = Md::log_beta_density(discount, discount_prior_strength_, discount_prior_beta_);
+    if (has_strength_prior())
+      lp += Md::log_gamma_density(strength + discount, strength_prior_shape_, strength_prior_rate_);
     assert(lp <= 0.0);
     if (num_customers_) {
-      if (d > 0.0) {
-        const double r = lgamma(1.0 - d);
-        lp += lgamma(alpha) - lgamma(alpha + num_customers_)
-             + num_tables_ * log(d) + lgamma(alpha / d + num_tables_)
-             - lgamma(alpha / d);
+      if (discount > 0.0) {
+        const double r = lgamma(1.0 - discount);
+        if (strength)
+          lp += lgamma(strength) - lgamma(strength / discount);
+        lp += - lgamma(strength + num_customers_)
+             + num_tables_ * log(discount) + lgamma(strength / discount + num_tables_);
         assert(std::isfinite(lp));
         for (typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.begin();
              it != dish_locs_.end(); ++it) {
           const DishLocations& cur = it->second;
           for (std::list<TableCount>::const_iterator ti = cur.table_counts_.begin(); ti != cur.table_counts_.end(); ++ti) {
-            lp += lgamma(ti->count - d) - r;
+            lp += lgamma(ti->count - discount) - r;
           }
         }
+      } else if (!discount) { // discount == 0.0
+        lp += lgamma(strength) + num_tables_ * log(strength) - lgamma(strength + num_tables_);
+        assert(std::isfinite(lp));
+        for (typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.begin();
+             it != dish_locs_.end(); ++it) {
+          const DishLocations& cur = it->second;
+          lp += lgamma(cur.table_counts_.size());
+        }
       } else {
-        assert(!"not implemented yet");
+        assert(!"discount less than 0 detected!");
       }
     }
     assert(std::isfinite(lp));
@@ -251,20 +260,22 @@ class MFCR {
   }
 
   void resample_hyperparameters(MT19937* rng, const unsigned nloop = 5, const unsigned niterations = 10) {
-    assert(has_discount_prior() || has_alpha_prior());
+    assert(has_discount_prior() || has_strength_prior());
     DiscountResampler dr(*this);
-    ConcentrationResampler cr(*this);
+    StrengthResampler sr(*this);
     for (int iter = 0; iter < nloop; ++iter) {
-      if (has_alpha_prior()) {
-        alpha_ = slice_sampler1d(cr, alpha_, *rng, 0.0,
+      if (has_strength_prior()) {
+        strength_ = slice_sampler1d(sr, strength_, *rng, -discount_,
                                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(),
+        double min_discount = std::numeric_limits<double>::min();
+        if (strength_ < 0.0) min_discount = -strength_;
+        discount_ = slice_sampler1d(dr, discount_, *rng, min_discount,
                                1.0, 0.0, niterations, 100*niterations);
       }
     }
-    alpha_ = slice_sampler1d(cr, alpha_, *rng, 0.0,
+    strength_ = slice_sampler1d(sr, strength_, *rng, -discount_,
                              std::numeric_limits<double>::infinity(), 0.0, niterations, 100*niterations);
   }
 
@@ -272,15 +283,15 @@ class MFCR {
     DiscountResampler(const MFCR& crp) : crp_(crp) {}
     const MFCR& crp_;
     double operator()(const double& proposed_d) const {
-      return crp_.log_crp_prob(proposed_d, crp_.alpha_);
+      return crp_.log_crp_prob(proposed_d, crp_.strength_);
     }
   };
 
-  struct ConcentrationResampler {
-    ConcentrationResampler(const MFCR& crp) : crp_(crp) {}
+  struct StrengthResampler {
+    StrengthResampler(const MFCR& crp) : crp_(crp) {}
     const MFCR& crp_;
-    double operator()(const double& proposediscount_alpha) const {
-      return crp_.log_crp_prob(crp_.discount_, proposediscount_alpha);
+    double operator()(const double& proposediscount_strength) const {
+      return crp_.log_crp_prob(crp_.discount_, proposediscount_strength);
     }
   };
 
@@ -292,7 +303,7 @@ class MFCR {
   };
 
   void Print(std::ostream* out) const {
-    (*out) << "MFCR(d=" << discount_ << ",alpha=" << alpha_ << ") customers=" << num_customers_ << std::endl;
+    (*out) << "MFCR(d=" << discount_ << ",strength=" << strength_ << ") customers=" << num_customers_ << std::endl;
     for (typename std::tr1::unordered_map<Dish, DishLocations, DishHash>::const_iterator it = dish_locs_.begin();
          it != dish_locs_.end(); ++it) {
       (*out) << it->first << " (" << it->second.total_dish_count_ << " on " << it->second.table_counts_.size() << " tables): ";
@@ -318,15 +329,15 @@ class MFCR {
   std::tr1::unordered_map<Dish, DishLocations, DishHash> dish_locs_;
 
   double discount_;
-  double alpha_;
+  double strength_;
 
   // optional beta prior on discount_ (NaN if no prior)
-  double discount_prior_alpha_;
+  double discount_prior_strength_;
   double discount_prior_beta_;
 
-  // optional gamma prior on alpha_ (NaN if no prior)
-  double alpha_prior_shape_;
-  double alpha_prior_rate_;
+  // optional gamma prior on strength_ (NaN if no prior)
+  double strength_prior_shape_;
+  double strength_prior_rate_;
 };
 
 template <typename T,typename H>