RefineLloyd.hpp

Go to the documentation of this file.

#ifndef KMEANS_LLOYD_HPP
#define KMEANS_LLOYD_HPP
 
#include <vector>
#include <algorithm>
 
#include "sanisizer/sanisizer.hpp"
 
#include "Refine.hpp"
#include "Details.hpp"
#include "QuickSearch.hpp"
#include "is_edge_case.hpp"
#include "compute_centroids.hpp"
#include "parallelize.hpp"
 
namespace kmeans {
 
struct RefineLloydOptions {
    int max_iterations = 10;
 
    int num_threads = 1;
};
 
template<typename Index_, typename Data_, typename Cluster_, typename Float_, typename Matrix_ = Matrix<Index_, Data_> >
class RefineLloyd final : public Refine<Index_, Data_, Cluster_, Float_, Matrix_> {
private:
    RefineLloydOptions my_options;
 
public:
    RefineLloyd(RefineLloydOptions options) : my_options(std::move(options)) {}
 
    RefineLloyd() = default;
 
public:
    RefineLloydOptions& get_options() {
        return my_options;
    }
 
public:
    Details<Index_> run(const Matrix_& data, const Cluster_ ncenters, Float_* const centers, Cluster_* const clusters) const {
        const auto nobs = data.num_observations();
        if (internal::is_edge_case(nobs, ncenters)) {
            return internal::process_edge_case(data, ncenters, centers, clusters);
        }
 
        auto sizes = sanisizer::create<std::vector<Index_> >(ncenters);
        auto copy = sanisizer::create<std::vector<Cluster_> >(nobs);
 
        const auto ndim = data.num_dimensions();
        internal::QuickSearch<Float_, Cluster_> index;
 
        decltype(I(my_options.max_iterations)) iter = 0;
        for (; iter < my_options.max_iterations; ++iter) {
            index.reset(ndim, ncenters, centers);
            parallelize(my_options.num_threads, nobs, [&](const int, const Index_ start, const Index_ length) -> void {
                auto work = data.new_extractor(start, length);
                for (Index_ obs = start, end = start + length; obs < end; ++obs) {
                    const auto dptr = work->get_observation();
                    copy[obs] = index.find(dptr); 
                }
            });
 
            // Checking if it already converged.
            bool updated = false;
            for (Index_ obs = 0; obs < nobs; ++obs) {
                if (copy[obs] != clusters[obs]) {
                    updated = true;
                    break;
                }
            }
            if (!updated) {
                break;
            }
            std::copy(copy.begin(), copy.end(), clusters);
 
            std::fill(sizes.begin(), sizes.end(), 0);
            for (Index_ obs = 0; obs < nobs; ++obs) {
                ++sizes[clusters[obs]];
            }
            internal::compute_centroids(data, ncenters, centers, clusters, sizes);
        }
 
        int status = 0;
        if (iter == my_options.max_iterations) {
            status = 2;
        } else {
            ++iter; // make it 1-based.
        }
        return Details<Index_>(std::move(sizes), iter, status);
    }
};
 
}
 
#endif