vonMisesDistribution.hpp

Go to the documentation of this file.

#ifndef VONMISESDISTRIBUTION_HPP
#define VONMISESDISTRIBUTION_HPP

#include <cmath>
#include <unsupported/Eigen/NonLinearOptimization>
#include <boost/math/special_functions/bessel.hpp>
#include <stdexcept>
#include <iostream>
#include <canopy/circularRegressor/vonMisesKappaFunctor.hpp>

namespace canopy
{

class vonMisesDistribution
{
    public:
        // Methods
        // -------

        vonMisesDistribution()
        : mu(0.0), kappa(0.0), S(0.0), C(0.0), pdf_normaliser(1.0)
        {
            // Nothing to do here
        }

        void initialise()
        {
            mu = 0.0;
            kappa = 0.0;
            S = 0.0;
            C = 0.0;
            pdf_normaliser = 1.0;
        }

        template <class TLabelIterator, class TIdIterator>
        void fit(TLabelIterator first_label, const TLabelIterator last_label, TIdIterator /*unused*/)
        {
            S = 0.0, C = 0.0;

            const int n_data = std::distance(first_label,last_label);

            // Loop through the data points, accumulating statistics
            for( ; first_label != last_label; ++first_label)
            {
                S += std::sin(*first_label);
                C += std::cos(*first_label);
            }

            // Length of resultant vector
            const double R2 = S*S + C*C;

            // Mean direction
            mu = std::atan2(S,C);

            // Unbiased R statistic
            const double Re = std::sqrt(R2) / n_data;

            // Find the kappa parameter
            if(Re > 0.98)
            {
                // There appears to be no solution for kappa in this case (look into this further!)
                // Saturate at roughly the value for when Re = 0.98
                kappa = 25.0;
            }
            else
            {
                // Set up and solve the non-linear equation for kappa
                vonMisesKappaFunctor vmftrinstance(Re);
                Eigen::VectorXd kappa_vec(1);
                kappa_vec << 25.0;
                Eigen::HybridNonLinearSolver<vonMisesKappaFunctor> solver(vmftrinstance);
                solver.hybrj1(kappa_vec);
                kappa = kappa_vec(0,0);
                pdf_normaliser = 1.0/(2.0*M_PI*boost::math::cyl_bessel_i(0,kappa));
            }
        }


        template<class TId>
        float pdf(const float x, const TId /*id*/) const
        {
            return pdf_normaliser*std::exp(kappa*std::cos(x - mu));
        }

        float pdf(const float x) const
        {
            return pdf(x,0);
        }

        void printOut(std::ofstream& stream) const
        {
            stream << mu << " " << kappa;
        }

        void readIn(std::ifstream& stream)
        {
            stream >> mu;
            stream >> kappa;

            S = std::sin(mu);
            C = std::cos(mu);
            pdf_normaliser = 1.0/(2.0*M_PI*boost::math::cyl_bessel_i(0,kappa));
        }

        float entropy() const
        {
            using boost::math::cyl_bessel_i;
            return std::log(2.0*M_PI*cyl_bessel_i(0,kappa)) - kappa*cyl_bessel_i(1,kappa)/cyl_bessel_i(0,kappa);
        }

        template <class TId>
        void combineWith(const vonMisesDistribution& dist, const TId /*id*/)
        {
            // Add the weighted mu value to the sine and cosine sums
            S += (dist.kappa)*(dist.S);
            C += (dist.kappa)*(dist.C);
        }

        // Normalises a distribution to find the resulting mu and kappa
        // Uses the sensor fusion approach of Stienne 2011
        void normalise()
        {
            mu = std::atan2(S,C);
            kappa = std::hypot(S,C);

            // Don't want this to raise an exception if kappa is very large (can happen
            // if there are many trees)
            try
            {
                const double kappa_bessel = boost::math::cyl_bessel_i(0,kappa);
                pdf_normaliser = 1.0/(2.0*M_PI*kappa_bessel);
            }
            catch (std::overflow_error)
            {
                kappa = 500.0;
                pdf_normaliser = 6.35397e-217;
            }
        }

        void reset()
        {
            initialise();
        }

        float getMu() const
        {
            return mu;
        }

        float getKappa() const
        {
            return kappa;
        }

        friend std::ofstream& operator<< (std::ofstream& stream, const vonMisesDistribution& dist) { dist.printOut(stream); return stream;}

        friend std::ifstream& operator>> (std::ifstream& stream, vonMisesDistribution& dist) { dist.readIn(stream); return stream;}

    protected:
        // Data
        float mu; 
        float kappa; 
        double S; 
        double C; 
        float pdf_normaliser; 

};

} // end of namespace

#endif
// VONMISESDISTRIBUTION_HPP