EM.hpp

#pragma once
/* (C) 2020 Roman Werpachowski. */
#include <memory>
#include <random>
#include <vector>
#include <Eigen/Cholesky>
#include <Eigen/Core>
#include "Clustering.hpp"
#include "dll.hpp"
 
 
namespace ml 
{   
    class EM: public Clustering::Model
    {
    public:     
        DLL_DECLSPEC EM(unsigned int number_components);
 
        DLL_DECLSPEC void set_seed(unsigned int seed);
 
        DLL_DECLSPEC void set_absolute_tolerance(double absolute_tolerance);
 
        DLL_DECLSPEC void set_relative_tolerance(double relative_tolerance);
 
        DLL_DECLSPEC void set_maximum_steps(unsigned int maximum_steps);
 
        DLL_DECLSPEC void set_means_initialiser(std::shared_ptr<const Clustering::CentroidsInitialiser> means_initialiser);
 
        DLL_DECLSPEC void set_responsibilities_initialiser(std::shared_ptr<const Clustering::ResponsibilitiesInitialiser> responsibilities_initialiser);
 
        void set_verbose(bool verbose)
        {
            verbose_ = verbose;
        }
 
        void set_maximise_first(bool maximise_first)
        {
            maximise_first_ = maximise_first;
        }
 
        DLL_DECLSPEC bool fit(Eigen::Ref<const Eigen::MatrixXd> data) override;
 
        auto number_components() const
        {
            return number_components_;
        }
 
        unsigned int number_clusters() const override
        {           
            return number_components();
        }
 
        const auto& means() const 
        {
            return means_;
        }
 
        const Eigen::MatrixXd& centroids() const override
        {           
            return means();
        }
 
        const auto& covariances() const
        {
            return covariances_;
        }
 
        DLL_DECLSPEC const Eigen::MatrixXd& covariance(unsigned int k) const;
 
        const auto& mixing_probabilities() const 
        {
            return mixing_probabilities_;
        }
 
        const auto& responsibilities() const 
        {
            return responsibilities_;
        }
 
        double log_likelihood() const 
        {
            return log_likelihood_;
        }       
 
        std::shared_ptr<const Clustering::CentroidsInitialiser> means_initialiser() const
        {
            return means_initialiser_;
        }
 
        DLL_DECLSPEC void assign_responsibilities(Eigen::Ref<const Eigen::VectorXd> x, Eigen::Ref<Eigen::VectorXd> u) const;        
 
        const std::vector<unsigned int>& labels() const override
        {
            return labels_;
        }
 
        bool converged() const override
        {
            return converged_;
        }
    private:
        std::default_random_engine prng_;
        std::shared_ptr<const Clustering::CentroidsInitialiser> means_initialiser_;
        std::shared_ptr<const Clustering::ResponsibilitiesInitialiser> responsibilities_initialiser_;
        Eigen::VectorXd mixing_probabilities_;
        Eigen::MatrixXd means_; 
        Eigen::MatrixXd responsibilities_; 
        Eigen::VectorXd work_vector_;
        std::vector<Eigen::MatrixXd> covariances_; 
        std::vector<Eigen::MatrixXd> inverse_covariances_; 
        std::vector<Eigen::LLT<Eigen::MatrixXd>> covariance_decompositions_; 
        Eigen::VectorXd sqrt_covariance_determinants_; 
        std::vector<unsigned int> labels_; 
        double absolute_tolerance_;
        double relative_tolerance_;
        double log_likelihood_;
        unsigned int number_components_;
        unsigned int maximum_steps_;
        bool verbose_;
        bool maximise_first_;
        bool converged_;
 
        static Eigen::MatrixXd calculate_sample_covariance(Eigen::Ref<const Eigen::MatrixXd> data);
 
        void process_covariances(Eigen::Index number_dimensions);
 
        void expectation_step(Eigen::Ref<const Eigen::MatrixXd> data);
 
        void maximisation_step(Eigen::Ref<const Eigen::MatrixXd> data);
 
        void calculate_labels();
    };
}