3
\$\begingroup\$

[NOTE] This question can be depreciated in favor of version 0.2.

This code works well.

This my first attempt at creating a robust, computationally lean K-Means Cluster analyzer. I first saw this algorithm in an intermediate programming course, and it was bloated to the max: it created clusters object into which points were physically copied from a universe. Centroids were computed and a scoring mechanism was implemented. This was computationally very expensive - it took an i7-2600 approximately 17 minutes to analyze a 5,000 point, 3-dimensional dataset.

This updated version analyzes a 50,000 point, 3-dimensional dataset in 30 seconds on the same processor.

Instead of using cluster objects, points are imported to the universe and are never moved. They are inserted within a wrapper object that points them toward whichever cluster they belong. Cluster wrappers contain their respective centroid, a count of the number of points belonging, and a move flag to prohibit the recomputation of centroids belonging to clusters that have seen no point movement in the current iteration.

Instead of a scoring mechanism, the algorithm exits when one of three conditions are met:

  • The maximum number of user-defined iterations is reached.

  • The maximum amount of user-defined bold error (epsilon) is achieved.

  • No point movement occurs between clusters.

Future Implementation

  • I'm leaving reporting in this version for the time being, however will be deleted in future versions.
  • Currently no exporting to file support - I'm still deciding how I want to do this (delimiters in file; file bursting for separate clusters).

Notes

kmeans.h

#ifndef KMEANS_KMEANS_H
#define KMEANS_KMEANS_H

#include <fstream>
#include <list>
#include <vector>



struct cluster_wrapper {
    std::vector<double> centroid;
    unsigned int count;
    bool move_flag;   // Signifies if anything moved in or out of a cluster - avoids unneeded centroid recomputation
    cluster_wrapper() : count(0), move_flag(1) {}
};



struct point_wrapper {
    std::vector<double> point;
    cluster_wrapper * centroid_pointer;
    point_wrapper() : centroid_pointer(nullptr) {}
};



class kmeans_set {
    std::list<point_wrapper> universe;      // The universe of points, a vector of vectors (vector of points)
    std::list<cluster_wrapper> clusters;   // Clusters with pointers to each point

    double epsilon;            // The maximum level of error (default 0)
    unsigned int hard_limit;   // Hard limit on the number of iterations (defaults to max of unsigned int)

    void import_points(std::list<std::vector<double>> &);   // Used by both import functions

    double calculate_distance(std::vector<double> &, std::vector<double> &);   // 2-norm distance calculation
    void move_to_nearest();                                                    // "Moves" points to the nearest cluster/centroid via pointer manipulation
    double recompute_centroids();
    bool check_move_state();                                                   // Checks if points have moved in/out of a cluster - avoids centroid recomputing when no movement
    void print_point(std::ostream &, std::vector<double>, char);               // A quick way to observe output

public:
    kmeans_set(std::ifstream &, char);                // Import points from file with specified delimiter
    kmeans_set(std::list<std::vector<double>> &);   // Import points from vector on construction

    void compute_centroids(unsigned int);
    void print_centroids(std::ostream &, char);       // output stream, centroid delimiter
    void print_clusters(std::ostream &, char, char);  // output stream, point delimiter, centroid designator

    void set_epsilon(double eps) { this->epsilon = eps; }
    void set_hard_limit(unsigned int lim) { this->hard_limit = lim; }
};

#endif //KMEANS_KMEANS_H

kmeans.cpp

#include "kmeans.h"
#include <cstdlib>
#include <limits>
#include <sstream>
#include <cmath>
#include <iostream>


// Import points from file with specified delimiter into a temporary list - calls import_points()
kmeans_set::kmeans_set(std::ifstream & input_file, char delimiter) :
        epsilon(0),
        hard_limit(std::numeric_limits<unsigned int>::max()) {
    std::string line;
    std::list<std::vector<double>> temp_point_list;
    while (std::getline(input_file, line)) {
        while ((line.length() == 0) && !(input_file.eof())) {
            std::getline(input_file, line);   // Skips blank lines in file
        }
        std::string parameter;
        std::stringstream ss(line);
        std::vector<double> temp_point;
        if ((line.length() != 0)) {
            while (std::getline(ss, parameter, delimiter)) {
                temp_point.push_back(atof(parameter.c_str()));
            }
            temp_point_list.push_back(temp_point);
        }
    }

// REPORTING - STARTS
/*
    std::cout << "\nFILE INTAKE\ntemp_list ct: " << temp_point_list.size();   // Reporting
    unsigned int ct = 0;
    for(auto temp_list_iter = temp_point_list.begin(); temp_list_iter != temp_point_list.end(); ++temp_list_iter) {
        std::cout << "\n" << ++ct << ": ";
        this->print_point(std::cout, *temp_list_iter, ',');
    }
    */
// REPORTING - ENDS

    this->import_points(temp_point_list);
}



// Import points from list on construction
kmeans_set::kmeans_set(std::list<std::vector<double>> & point_list) :
        epsilon(0),
        hard_limit(std::numeric_limits<unsigned int>::max()) {
    this->import_points(point_list);
}



// Primary point import function - Assures dimensional integrity by comparing all intake to the first point in the [points] list
void kmeans_set::import_points(std::list<std::vector<double>> & point_list) {
    for (auto point_list_iter = point_list.begin(); point_list_iter != point_list.end(); ++point_list_iter) {
        point_wrapper new_point;
        new_point.point = *point_list_iter;
        if (this->universe.empty()) {
            this->universe.push_back(new_point);
        } else if (this->universe.front().point.size() == point_list_iter->size()) {
            // Assures dimensional integrity
            // Sorting on intake assures well-dispursed centroid selection
            // Well-dispursed centroid selection helps mitigate centroids "fighting" over points
            auto universe_list_iter = this->universe.begin();
            while(universe_list_iter != this->universe.end() && (universe_list_iter->point < *point_list_iter)) {   // PRECEDENCE MATTERS!!
                // Sorts points on insert - this ensures evenly distributed centroid picking
                    ++universe_list_iter;
            }
            this->universe.insert(universe_list_iter,new_point);
        }
    }

    // REPORTING
    /*
    std::cout << "\n\nIMPORTING\nuniverse ct: " << this->universe.size();   // Reporting
    unsigned int ct = 0;
    for(auto universe_iter = this->universe.begin(); universe_iter != this->universe.end(); ++universe_iter) {
        std::cout << "\n" << ++ct << ": ";
        this->print_point(std::cout, universe_iter->point, ',');
        std::cout << "   centroid ptr: " << universe_iter->centroid_pointer;
    }
     */
    // REPORTING - ENDS

}



// Primary algorithm
void kmeans_set::compute_centroids(unsigned int k) {

    // Error handling
    if (this->universe.empty() || k == 0) {   // Exits program if no points have been imported or if no centroids are requested
        return;
    } else if (k > this->universe.size()) {   // Forces a maximum amount of centroids in a low-population universe.
        k = this->universe.size();
    }

    // Distributes centroids evenly among points
    {   // Encapsulation for variable destruction
        unsigned int subset = this->universe.size() / k;   // Used for even distribution
        unsigned int rem = this->universe.size() % k;      // Used for even distribution
        unsigned int i = 0;                                // Used for even distribution
        for (auto centroid_pick_iter = this->universe.begin(); centroid_pick_iter != this->universe.end(); ++centroid_pick_iter, --i) {
            if (i == 0) {
                cluster_wrapper new_cluster;
                new_cluster.centroid = centroid_pick_iter->point;
                this->clusters.push_back(new_cluster);
                i = subset + (rem ? 1 : 0);
                if (rem) { --rem; }
            }
            centroid_pick_iter->centroid_pointer = &(this->clusters.back());
            ++centroid_pick_iter->centroid_pointer->count;
        }
    }

    // REPORTING
    /*
    {
        std::cout << "\n\nASSIGNING CENTROIDS\nuniverse ct: " << this->universe.size();   // Reporting
        unsigned int ct = 0;
        for (auto universe_iter = this->universe.begin(); universe_iter != this->universe.end(); ++universe_iter) {
            std::cout << "\n" << ++ct << ": ";
            this->print_point(std::cout, universe_iter->point, ',');
            std::cout << "   centroid ptr: " << universe_iter->centroid_pointer;
        }
        std::cout << "\n\ncluster ct: " << this->clusters.size();   // Reporting
        ct = 0;
        for (auto cluster_iter = this->clusters.begin(); cluster_iter != this->clusters.end(); ++cluster_iter) {
            std::cout << "\n" << ++ct << ": ";
            this->print_point(std::cout, cluster_iter->centroid, ',');
            std::cout << "   cluster ct: " << cluster_iter->count << "   Address: " << &(*cluster_iter);
        }
    }
     */
    // REPORTING

    // Primary algo
    unsigned int n = this->hard_limit;
    double delta_max_1 = std::numeric_limits<double>::max();
    double delta_max_2 = 0;

    unsigned int iter_ct = 0; // REPORTING
    recompute_centroids();

    while ((fabs(delta_max_2 - delta_max_1) > fabs(this->epsilon)) && n && check_move_state()) {
        delta_max_2 = delta_max_1;

        // Reset move flags
        for(auto cluster_iter = this->clusters.begin(); cluster_iter != this->clusters.end(); ++cluster_iter) {
            cluster_iter->move_flag = false;
        }

        // Reassign all points to nearest centroid
        this->move_to_nearest();

        // Recompute all centroids, returns the highest movement (delta)
        delta_max_1 = recompute_centroids();
        --n;

        // REPORTING
        {
            std::cout << "\n\nITERATION " << ++iter_ct << "   cluster ct: " << this->clusters.size();   // Reporting
            unsigned int ct = 0;
            for (auto cluster_iter = this->clusters.begin(); cluster_iter != this->clusters.end(); ++cluster_iter) {
                std::cout << "\n"  << "Move State: " << ((cluster_iter->move_flag)?("INVALID    "):("  valid    ")) << ++ct << ": ";
                this->print_point(std::cout, cluster_iter->centroid, ',');
                std::cout << "   cluster ct: " << cluster_iter->count << "   Address: " << &(*cluster_iter);
            }
        }
        // REPORTING

    }
}





void kmeans_set::move_to_nearest() {
    for (auto universe_iter = this->universe.begin(); universe_iter != this->universe.end(); ++universe_iter) {
        for (auto cluster_iter = this->clusters.begin(); cluster_iter != this->clusters.end(); ++cluster_iter) {

            // REPORTING
            /*
            std::cout << "\nComparing: (";
            this->print_point(std::cout,universe_iter->point,',');
            std::cout << ")-(";
            this->print_point(std::cout,cluster_iter->centroid,',');
            std::cout << ") dist: ";
            std::cout << calculate_distance(universe_iter->point, cluster_iter->centroid);
            std::cout << "\n       To: (";
            this->print_point(std::cout,universe_iter->point,',');
            std::cout << ")-(";
            this->print_point(std::cout,universe_iter->centroid_pointer->centroid,',');
            std::cout << ") dist: ";
            std::cout << calculate_distance(universe_iter->point, universe_iter->centroid_pointer->centroid);
             */
            // REPORTING

            if (calculate_distance(universe_iter->point, cluster_iter->centroid) < calculate_distance(universe_iter->point, universe_iter->centroid_pointer->centroid)) {

                // REPORTING
                /*
                std::cout << "\n-----------------------------------------\nSending: (";
                this->print_point(std::cout, universe_iter->point, ',');
                std::cout << ") From: " << universe_iter->centroid_pointer;
                 */
                // REPORTING

                --universe_iter->centroid_pointer->count;
                universe_iter->centroid_pointer->move_flag = true;
                universe_iter->centroid_pointer = &(*cluster_iter);  // This is why I'm using lists vs. vectors - iterator would be invalid
                ++universe_iter->centroid_pointer->count;
                universe_iter->centroid_pointer->move_flag = true;

                // REPORTING
                /*
                std::cout << "    To: " << universe_iter->centroid_pointer << "\n-----------------------------------------";
                // REPORTING
                 */
            }
        }
    }
}



double kmeans_set::calculate_distance(std::vector<double> & a_vec, std::vector<double> & b_vec) {
    auto b_iter = b_vec.begin();
    double r_val = 0;
    for (auto a_iter = a_vec.begin(); a_iter != a_vec.end(); ++a_iter, ++b_iter) {
        r_val += (*a_iter - *b_iter) * (*a_iter - *b_iter);
    }
    // return std::sqrt(r_val);
    return r_val;   // No need to root - actual distance not necessary
}



double kmeans_set::recompute_centroids() {

    // The following is used for delta calculation
    std::list<std::vector<double>> old_centroids;
    for (auto cluster_iter = this->clusters.begin(); cluster_iter != this->clusters.end(); ++cluster_iter) {
        old_centroids.push_back(cluster_iter->centroid);   // Copy the current centroids state
    }
    for (auto cluster_iter = this->clusters.begin(); cluster_iter != this->clusters.end(); ++cluster_iter) {
        std::fill(cluster_iter->centroid.begin(), cluster_iter->centroid.end(), 0);   // Zero out the primary centroid list
    }
    auto old_centroid_iter = old_centroids.begin();
    for (auto cluster_iter = this->clusters.begin(); cluster_iter != this->clusters.end(); ++cluster_iter, ++old_centroid_iter) {
        if (!(cluster_iter->move_flag)) {
            cluster_iter->centroid = *old_centroid_iter;  // Recopy back to the current cluster if there was no movement - will avoid recomputation.
        }
    }

    // Recompute centroids
    for (auto universe_iter = this->universe.begin(); universe_iter != this->universe.end(); ++universe_iter) {
        if (universe_iter->centroid_pointer->move_flag) {
            auto centroid_pointer_centroid_iter = universe_iter->centroid_pointer->centroid.begin();
            for (auto point_iter = universe_iter->point.begin(); point_iter != universe_iter->point.end(); ++point_iter, ++centroid_pointer_centroid_iter) {
                // Division on the fly useful for accuracy super-large dimensions - doubles computation
                *centroid_pointer_centroid_iter += *point_iter / universe_iter->centroid_pointer->count;
            }
        }
    }

    // Compute max delta
    double max_delta = 0;
    auto old_centroids_iter = old_centroids.begin();
    for (auto centroids_iter = this->clusters.begin(); centroids_iter != this->clusters.end(); ++centroids_iter) {
        double test_delta = calculate_distance(*old_centroids_iter, centroids_iter->centroid);
        if (max_delta < test_delta) { max_delta = test_delta; }
    }
    return max_delta;
}



void kmeans_set::print_centroids(std::ostream & output, char delimiter) {
    for(auto cluster_iter = this->clusters.begin(); cluster_iter != this->clusters.end(); ++cluster_iter) {
        for (auto centroid_iter = cluster_iter->centroid.begin(); centroid_iter != cluster_iter->centroid.end(); ++centroid_iter) {
            if(cluster_iter != this->clusters.begin() && centroid_iter == cluster_iter->centroid.begin()) {
                output << "\n";
            } else if(centroid_iter != cluster_iter->centroid.begin()) {
                output << delimiter;
            }
            output << *centroid_iter;
        }
    }
}



void kmeans_set::print_point(std::ostream & output, std::vector<double> vec, char delim) {
    for(auto iter = vec.begin(); iter != vec.end(); ++iter) {
        if(iter != vec.begin()) { output << delim; }
        output << *iter;
    }
}



bool kmeans_set::check_move_state() {
    for (auto cluster_iter = this->clusters.begin(); cluster_iter != this->clusters.end(); ++cluster_iter) {
        if (cluster_iter->move_flag) { return true; }
    }
    return false;
}

main.cpp

#include <iostream>
#include "kmeans.h"

int main() {

    std::ifstream input;
    input.open("test_data_1.dat");
    std::ostream &output_console = std::cout;


    kmeans_set km(input, ',');
    km.compute_centroids(5);

    std::cout << "\n\nOutput in main():\n";
    km.print_centroids(output_console, ',');


}
\$\endgroup\$
  • 1
    \$\begingroup\$ I believe it is better to put input and output files on github or similar, and leave a link here. The amount of data makes the post not possible to open on mobile. Also it doesn't add much value. May be just layout could be described. \$\endgroup\$ – Incomputable Aug 22 '17 at 4:00
  • \$\begingroup\$ @Incomputable - That's perfect, thank you! I wasn't sure what the best option would be for linking. I will modify momentarily... \$\endgroup\$ – Miller Aug 22 '17 at 4:05
  • \$\begingroup\$ Please note, this question can be depreciated in favor of version 0.2. \$\endgroup\$ – Miller Sep 2 '17 at 21:57
  • \$\begingroup\$ errors on all "_iter"s was not declared in this scope any idea where to declare the code? thanks \$\endgroup\$ – Sujacka Retno Dec 17 '18 at 20:02
  • \$\begingroup\$ also note the OPs comment about this code being deprecated in favor of version 0.2. \$\endgroup\$ – Sᴀᴍ Onᴇᴌᴀ Dec 17 '18 at 20:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.