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);
        const auto ndim = data.num_dimensions();
        internal::QuickSearch<Float_, Cluster_> index;
 
        auto num_diff_threads = sanisizer::create<std::vector<Index_> >(my_options.num_threads);
        std::fill_n(clusters, nobs, 0); // just to avoid using uninitialized values in the first iteration of the loop.
 
        I<decltype(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 t, const Index_ start, const Index_ length) -> void {
                auto work = data.new_extractor(start, length);
                Index_ num_diff = 0;
                for (Index_ obs = start, end = start + length; obs < end; ++obs) {
                    const auto dptr = work->get_observation();
                    const auto closest = index.find(dptr); 
                    auto& previous = clusters[obs];
                    num_diff += (closest != previous);
                    previous = closest;
                }
                num_diff_threads[t] = num_diff;
            });
 
            if (iter) {
                // Checking if it already converged.
                bool updated = false;
                for (const auto num_diff : num_diff_threads) {
                    if (num_diff) {
                        updated = true;
                        break;
                    }
                }
                if (!updated) {
                    break;
                }
            }
 
            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