kmeans.c

Go to the documentation of this file.

/*
 * This work is licensed under a Creative Commons 
 * Attribution-Noncommercial-Share Alike 3.0 United States License.
 * 
 *    http://creativecommons.org/licenses/by-nc-sa/3.0/us/
 * 
 * You are free:
 * 
 *    to Share - to copy, distribute, display, and perform the work
 *    to Remix - to make derivative works
 * 
 * Under the following conditions:
 * 
 *    Attribution. You must attribute the work in the manner specified by the
 *    author or licensor (but not in any way that suggests that they endorse you
 *    or your use of the work).
 * 
 *    Noncommercial. You may not use this work for commercial purposes.
 * 
 *    Share Alike. If you alter, transform, or build upon this work, you may
 *    distribute the resulting work only under the same or similar license to
 *    this one.
 * 
 * For any reuse or distribution, you must make clear to others the license
 * terms of this work. The best way to do this is by including this header.
 * 
 * Any of the above conditions can be waived if you get permission from the
 * copyright holder.
 * 
 * Apart from the remix rights granted under this license, nothing in this
 * license impairs or restricts the author's moral rights.
 */


#include <jwsc/config.h>

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <limits.h>
#include <inttypes.h>

#include "jwsc/vector/vector.h"
#include "jwsc/matrix/matrix.h"
#include "jwsc/prob/pdf.h"
#include "jwsc/stat/kmeans.h"


void kmeans_d
(
    Matrix_d**      means_out, 
    Vector_i32**    members_out, 
    const Matrix_d* features,
    uint32_t        K,
    uint32_t        N
)
{
    uint32_t mean, mean2, feature;
    uint32_t num_features = features->num_cols;
    uint32_t num_records = features->num_rows;
    uint32_t n;
    int64_t record;

    double objective = INT_MAX;
    double delta_objective = 1.0;
    double d, sum, min;

    Matrix_d*   means;
    Vector_i32* members;
    Vector_i32* member_counts;

    assert(N <= floor(features->num_rows / (double)K));

    if (*means_out == NULL)
    {
        init_kmeans_randomly_d(means_out, features, K);
    }
    means = *means_out;

    create_zero_vector_i32(members_out, num_records);
    members = *members_out;

    member_counts = NULL;

    /* Until there is no change in the means, iterate. */
    while (delta_objective > 1.0e-10)
    {
        /* Assign membership of each record to its closest mean. */
        for (record = 0; record < num_records; record++)
        {
            min = INT_MAX;

            for (mean = 0; mean < K; mean++)
            {
                sum = 0.0;

                for (feature = 0; feature < num_features; feature++)
                {
                    d = means->elts[ mean ][ feature ];
                    d -= features->elts[ record ][ feature ];
                    sum += d*d;
                }

                if  (sum < min)
                {
                    min = sum;
                    members->elts[ record ] = mean;
                }
            }
        }

        create_zero_vector_i32(&member_counts, K);

        /* Ensure each mean has at least N member records. */
        for (record = 0; record < num_records; record++) 
        {
            mean = members->elts[ record ];
            member_counts->elts[ mean ]++;
        }
        for (mean = 0; mean < K; mean++)
        {
            if (member_counts->elts[ mean ] < N)
            {
                n = N - member_counts->elts[ mean ];

                while (n > 0)
                {
                    record = floor(sample_uniform_pdf_d(0, num_records*0.9999));
                    mean2 = members->elts[ record ];

                    if (member_counts->elts[ mean2 ] > N)
                    {
                        members->elts[ record ] = mean;
                        member_counts->elts[ mean2 ]--;
                        member_counts->elts[ mean ]++;
                        n--;
                    }
                }
            }

            assert(member_counts->elts[ mean ] >= N);
        }

        create_zero_matrix_d(&means, means->num_rows, means->num_cols);

        /* Re-calculate means from their members. */
        for (record = 0; record < num_records; record++)
        {
            mean = members->elts[ record ];

            for (feature = 0; feature < num_features; feature++)
            {
                means->elts[ mean ][ feature ] += 
                    features->elts[ record ][ feature ];
            }
        }
        for (mean = 0; mean < K; mean++)
        {
            for (feature = 0; feature < num_features; feature++)
            {
                means->elts[ mean ][ feature ] /= member_counts->elts[ mean ];
            }
        }

        sum = 0.0;

        /* Re-caclculate value of the objective function. */
        for (record = 0; record < num_records; record++)
        {
            mean = members->elts[ record ];

            for (feature = 0; feature < num_features; feature++)
            {
                d = means->elts[ mean ][ feature];
                d -= features->elts[ record ][ feature ];
                sum += d*d;
            }
        }

        /*assert((delta_objective = objective - sum) >= 0);*/
        delta_objective = fabs(objective - sum);
        objective = sum;

        /*fprintf(stderr, "Objective: %e\n", objective);*/
    }

    free_vector_i32(member_counts);
}


void init_kmeans_randomly_d
(
    Matrix_d**      means_out, 
    const Matrix_d* features,
    uint32_t        K
)
{
    uint32_t mean, feature;
    uint32_t num_features = features->num_cols;
    uint32_t num_records = features->num_rows;
    int64_t  record;
    double r;

    Vector_i32* choosen_records = NULL;
    Matrix_d*   means;

    create_zero_vector_i32(&choosen_records, num_records);

    create_zero_matrix_d(means_out, K, num_features);
    means = *means_out;

    /* Randomly choose k records to initialize the means with. */
    for (mean = 0; mean < K; mean++)
    {
        record = -1;

        while (record < 0)
        {
            r = rand() / (double)RAND_MAX;
            record = floor((num_records-1)*r + 0.5);

            if (!choosen_records->elts[ record ])
            {
                choosen_records->elts[ record ] = 1;
            }
            else
            {
                record = -1;
            }
        }

        for (feature = 0; feature < num_features; feature++)
        {
            means->elts[ mean ][ feature ] = 
                features->elts[ record ][ feature ];
        }
    }

    free_vector_i32(choosen_records);
}