still working on sum of delta

include/fl/distribution/sum_of_deltas.hpp

@@ -23,7 +23,7 @@
 #ifndef FL__DISTRIBUTION__SUM_OF_DELTAS_HPP
 #define FL__DISTRIBUTION__SUM_OF_DELTAS_HPP
 
-#include <Eigen/Dense>
+#include <Eigen/Core>
 
 // std
 #include <vector>
@@ -86,7 +86,7 @@ struct Traits<SumOfDeltas<Var>>
     /**
      * \brief Weight vector associated with the deltas
      */
-    typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Probabilities;
+    typedef Eigen::Array<Scalar, Eigen::Dynamic, 1> Probabilities;
 
     /**
      * \brief Moments interface of the SumOfDeltas distribution
@@ -113,7 +113,9 @@ class SumOfDeltas
     typedef typename Traits<This>::Locations    Locations;
     typedef typename Traits<This>::Probabilities      Probabilities;
 
+
 public:
+    // constructor and destructor ----------------------------------------------
     /**
      * Creates a dynamic or fixed-size SumOfDeltas.
      *
@@ -128,71 +130,97 @@ class SumOfDeltas
     {
         locations_ = Locations(1, Variate::Zero(dim));
         log_probabilities_ = Probabilities::Zero(1);
+        cumulative_ = Probabilities::Ones(1);
     }
 
-    /**
-     * \brief Overridable default constructor
-     */
     virtual ~SumOfDeltas() { }
 
 
 
-//    virtual void set_log_unnormalized_probability(const Weights& log_prob)
-//    {
-//        double max = *(std::max_element(log_prob.begin(), log_prob.end()));
-
-//        Weights prob(log_prob.size());
-//        for(int i = 0; i < int(log_prob.size()); i++)
-//            prob[i] = std::exp(log_prob[i] - max);
 
-//        set_unnormalized_probability(prob);
-//    }
+    /// non-const functions ****************************************************
 
-//    virtual void set_unnormalized_probability(const Weights& probs)
-//    {
-//        Scalar sum = 0;
-//        for(size_t i = 0; i < probs.size(); i++)
-//        {
-//            sum += prob[i];
-//        }
-
-//        weights_.resize(probs.size());
-//        for(size_t i = 0; i < probs.size(); i++)
-//        {
-//            weights_[i] = prob[i] / sum;
-//        }
-//    }
+    // set ---------------------------------------------------------------------
+    virtual void log_unnormalized_probabilities(const Probabilities& log_probs)
+    {
+        std::cout << "log_probs.size()" << log_probs.size()  << std::endl;
+        // rescale for numeric stability
+        log_probabilities_ = log_probs;
+        set_max(log_probabilities_);
+\
+        // copy to probabilities
+        probabilities_ = log_probabilities_.exp();
+        Scalar sum = probabilities_.sum();
+
+        // normalize
+        probabilities_ /= sum;
+        log_probabilities_ -= std::log(sum);
+
+        // compute cumulative
+        cumulative_.resize(log_probs.size());
+        cumulative_[0] = probabilities_[0];
+        for(int i = 1; i < cumulative_.size(); i++)
+            cumulative_[i] = cumulative_[i-1] + probabilities_[i];
+    }    
 
     virtual Variate& location(size_t i)
     {
         return locations_[i];
     }
 
-    virtual Scalar& log_probability(size_t i)
+    virtual void resize(size_t dim)
     {
-        return log_probabilities_(i);
+        locations_.resize(dim);
+        log_probabilities_.resize(dim);
     }
 
-    virtual Scalar probability(size_t i) const
+
+    /// const functions ********************************************************
+
+    // sampling ----------------------------------------------------------------
+    Variate map_standard_normal(Scalar gaussian_sample) const
     {
-        return std::exp(log_probabilities_(i));
+        Scalar uniform_sample =
+                0.5 * (1.0 + std::erf(gaussian_sample / std::sqrt(2.0)));
+
+        return map_standard_uniform(uniform_sample);
     }
 
-    virtual void resize(size_t dim)
+    Variate map_standard_uniform(Scalar uniform_sample) const
     {
-        locations_.resize(dim);
-        log_probabilities_.resize(dim);
+        typename std::vector<Scalar>::const_iterator
+                iterator = std::lower_bound(cumulative_.begin(),
+                                            cumulative_.end(),
+                                            uniform_sample);
+
+        int index = iterator - cumulative_.begin();
+        return locations_[index];
+    }
+
+
+    // get ---------------------------------------------------------------------
+    virtual Scalar log_probability(const size_t& i) const
+    {
+        return log_probabilities_(i);
+    }
+
+    virtual Scalar probability(const size_t& i) const
+    {
+        return std::exp(log_probabilities_(i));
     }
 
     virtual size_t size() const
     {
         return locations_.size();
     }
 
-    /**
-     * \return The weighted mean of the locations, or simply the first moment of
-     *         the distribution.
-     */
+    virtual int dimension() const
+    {
+        return locations_[0].rows();
+    }
+
+
+    // compute properties ------------------------------------------------------
     virtual Variate mean() const
     {
         Variate mu(Variate::Zero(dimension()));
@@ -202,33 +230,25 @@ class SumOfDeltas
         return mu;
     }
 
-    /**
-     * \return The covariance or the second central moment of the distribution
-     */
     virtual Covariance covariance() const
     {
         Variate mu = mean();
         Covariance cov(Covariance::Zero(dimension(), dimension()));
         for(size_t i = 0; i < locations_.size(); i++)
-            cov += probability(i) * (locations_[i]-mu) * (locations_[i]-mu).transpose();
+            cov += probability(i) *
+                    (locations_[i]-mu) * (locations_[i]-mu).transpose();
 
         return cov;
     }
 
-    /**
-     * \return Dimension of the distribution variate
-     */
-    virtual int dimension() const
-    {
-        return locations_[0].rows();
-    }
-
     virtual Scalar entropy() const
     {
         Scalar ent = 0;
-        for(size_t i = 0; i < log_probabilities_.size(); i++)
+        for(int i = 0; i < log_probabilities_.size(); i++)
         {
-            Scalar summand = - log_probabilities_(i) * std::exp(log_probabilities_(i));
+            Scalar summand =
+                    - log_probabilities_(i) * std::exp(log_probabilities_(i));
+
             if(!std::isfinite(summand))
                 summand = 0; // the limit for weight -> 0 is equal to 0
             ent += summand;
@@ -237,11 +257,24 @@ class SumOfDeltas
         return ent;
     }
 
+
+
+
+
+protected:
+    virtual void set_max(Probabilities& p, const Scalar& max = 0) const
+    {
+        const Scalar old_max = p.maxCoeff();
+        p += max - old_max;
+    }
+
+
 protected:
     Locations  locations_;
-    Probabilities probabilities_;
 
     Probabilities log_probabilities_;
+    Probabilities probabilities_;
+    Probabilities cumulative_;
 };
 
 }

test/distribution/discrete_distribution_test.cpp

@@ -81,7 +81,7 @@ TEST(discrete_distribution, sampling)
     for(size_t i = 0; i < sum_of_delta.size(); i++)
     {
         sum_of_delta.location(i) = gaussian.sample();
-        sum_of_delta.weight(i) = 1. / sum_of_delta.size();
+        sum_of_delta.probability(i) = 1. / sum_of_delta.size();
     }
 
     // compare mean and covariance
@@ -96,7 +96,7 @@ TEST(discrete_distribution, sampling)
     // check entropy of uniform distribution
     for(size_t i = 0; i < sum_of_delta.size(); i++)
     {
-        sum_of_delta.weight(i) = 1. / sum_of_delta.size();
+        sum_of_delta.probability(i) = 1. / sum_of_delta.size();
     }
 
     EXPECT_TRUE(fabs(std::log(double(sum_of_delta.size()))
@@ -105,9 +105,9 @@ TEST(discrete_distribution, sampling)
     // check entropy of certain distribution
     for(size_t i = 0; i < sum_of_delta.size(); i++)
     {
-        sum_of_delta.weight(i) = 0;
+        sum_of_delta.probability(i) = 0;
     }
-    sum_of_delta.weight(0) = 1;
+    sum_of_delta.probability(0) = 1;
 
     EXPECT_TRUE(fabs(sum_of_delta.entropy()) < 0.0000001);
 }

test/distribution/sum_of_delta_test.cpp

@@ -49,6 +49,7 @@
 
 #include <fl/distribution/sum_of_deltas.hpp>
 #include <fl/distribution/gaussian.hpp>
+
 #include <Eigen/Core>
 #include <Eigen/Eigenvalues>
 
@@ -59,6 +60,10 @@ TEST(sum_of_delta, mean_and_covariance)
 {
     typedef Eigen::Vector3d Variate;
     typedef Eigen::Matrix3d Covariance;
+    typedef Variate::Scalar Scalar;
+
+    typedef fl::SumOfDeltas<Variate> DiscreteDistribution;
+    typedef DiscreteDistribution::Probabilities Function;
 
     // pick some mean and covariance
     Covariance covariance;
@@ -76,12 +81,14 @@ TEST(sum_of_delta, mean_and_covariance)
     gaussian.covariance(covariance);
 
     // generate a sum of delta from gaussian
-    fl::SumOfDeltas<Variate> sum_of_delta;
-    sum_of_delta.resize(100000);
+    DiscreteDistribution sum_of_delta;
+        sum_of_delta.resize(100000);
+
+    sum_of_delta.log_unnormalized_probabilities(Function::Zero(100000));
+
     for(size_t i = 0; i < sum_of_delta.size(); i++)
     {
         sum_of_delta.location(i) = gaussian.sample();
-        sum_of_delta.log_probability(i) = std::log(1. / sum_of_delta.size());
     }
 
     // compare mean and covariance
@@ -94,20 +101,17 @@ TEST(sum_of_delta, mean_and_covariance)
                                     (sum_of_delta.mean()-mean)).norm() < 0.1);
 
     // check entropy of uniform distribution
-    for(size_t i = 0; i < sum_of_delta.size(); i++)
-    {
-        sum_of_delta.log_probability(i) = std::log(1. / sum_of_delta.size());
-    }
+    sum_of_delta.log_unnormalized_probabilities(Function::Zero(sum_of_delta.size()));
 
     EXPECT_TRUE(fabs(std::log(double(sum_of_delta.size()))
                      - sum_of_delta.entropy()) < 0.0000001);
 
     // check entropy of certain distribution
-    for(size_t i = 0; i < sum_of_delta.size(); i++)
-    {
-        sum_of_delta.log_probability(i) = -std::numeric_limits<double>::max();
-    }
-    sum_of_delta.log_probability(0) = 0;
+    Function log_pmf = Function::Constant(sum_of_delta.size(),
+                                          -std::numeric_limits<double>::max());
+    log_pmf(0) = 0;
+    sum_of_delta.log_unnormalized_probabilities(log_pmf);
+
 
     EXPECT_TRUE(fabs(sum_of_delta.entropy()) < 0.0000001);
 }