+++ /dev/null
-#include "Util/types.h"
-#include "Util/rng.h"
-#include <cgds/Vector.h>
-#include <stdlib.h>
-#include <math.h>
-
-// auxiliary function to obtain a random sample of 1..n with k elements
-static void sample(uint32_t* v, uint32_t n, uint32_t k)
-{
- // refVect = (0,1,...,n-1,n)
- uint32_t* refVect = (uint32_t*) malloc(n * sizeof(uint32_t));
- for (uint32_t i = 0; i < n; i++)
- refVect[i] = i;
-
- uint32_t curSize = n; // current size of the sampling set
- for (uint32_t j = 0; j < k; j++)
- {
- // pick an index in sampling set:
- uint32_t index = get_rand_int() % curSize;
- v[j] = refVect[index];
- // move this index outside of sampling set:
- refVect[index] = refVect[--curSize];
- }
-
- free(refVect);
-}
-
-// assign a vector (represented by its dissimilarities to others, as dissimilarities[index,])
-// to a cluster, represented by its center ==> output is integer in 0..K-1
-static uint32_t assignCluster(uint32_t index, float* dissimilarities,
- uint32_t* centers, uint32_t n, uint32_t K)
-{
- uint32_t minIndex = 0;
- float minDist = dissimilarities[index * n + centers[0]];
-
- for (uint32_t j = 1; j < K; j++)
- {
- if (dissimilarities[index * n + centers[j]] < minDist)
- {
- minDist = dissimilarities[index * n + centers[j]];
- minIndex = j;
- }
- }
-
- return minIndex;
-}
-
-// assign centers given a clustering, and also compute corresponding distortion
-static void assign_centers(uint32_t nbClusters, Vector** clusters, float* dissimilarities,
- uint32_t nbItems, uint32_t* ctrs, float* distor)
-{
- *distor = 0.0;
- // TODO [heuristic]: checking only a neighborhood of the former center ?
- for (uint32_t j = 0; j < nbClusters; j++)
- {
- // If the cluster is empty, choose a center at random (pathological case...)
- uint32_t minIndex = get_rand_int() % nbItems;
- float minSumDist = INFINITY;
- for (uint32_t i = 0; i < vector_size(clusters[j]); i++)
- {
- uint32_t index1;
- vector_get(clusters[j], i, index1);
- // attempt to use current index as center
- float sumDist = 0.0;
- for (uint32_t ii = 0; ii < vector_size(clusters[j]); ii++)
- {
- uint32_t index2;
- vector_get(clusters[j], ii, index2);
- sumDist += dissimilarities[index1 * nbItems + index2];
- }
- if (sumDist < minSumDist)
- {
- minSumDist = sumDist;
- minIndex = index1;
- }
- }
- ctrs[j] = minIndex;
- *distor += minSumDist;
- }
-}
-
-// Core PAM algorithm from a dissimilarity matrix; (e.g. nstart=10, maxiter=100)
-void pam(float* dissimilarities, uint32_t nbItems, uint32_t nbClusters, int clustOnMedoids,
- uint32_t nbStart, uint32_t maxNbIter, Result_t* result)
-{
- uint32_t* ctrs = result->medoids_ranks; //shorthand
- uint32_t* oldCtrs = (uint32_t*) malloc(nbClusters * sizeof(uint32_t));
- Vector** clusters = (Vector**) malloc(nbClusters * sizeof(Vector*));
- Vector** bestClusts = (Vector**) malloc(nbClusters * sizeof(Vector*));
- for (uint32_t j = 0; j < nbClusters; j++)
- {
- clusters[j] = vector_new(uint32_t);
- bestClusts[j] = vector_new(uint32_t);
- }
-
- float lastDistor, distor, bestDistor = INFINITY;
- for (uint32_t startKount = 0; startKount < nbStart; startKount++)
- {
- // centers (random) [re]initialization
- if (clustOnMedoids)
- {
- // In this special case (clustering groups of medoids), we can improve the sampling
- uint32_t randomMedoidRank = get_rand_int() % (nbItems / nbClusters);
- uint32_t startIndex = randomMedoidRank * nbClusters;
- for (uint32_t index = startIndex; index < startIndex + nbClusters; index++)
- ctrs[index - startIndex] = index;
- }
- else
- // No information: complete random sampling
- sample(ctrs, nbItems, nbClusters);
- for (uint32_t j = 0; j < nbClusters; j++)
- oldCtrs[j] = 0;
- uint32_t kounter = 0;
-
- distor = INFINITY;
- do // while 'distortion' decreases...
- {
- // (re)initialize clusters to empty sets
- for (uint32_t j = 0; j < nbClusters; j++)
- vector_clear(clusters[j]);
-
- // estimate clusters belongings
- for (uint32_t i = 0; i < nbItems; i++)
- {
- uint32_t affectation = assignCluster(i, dissimilarities, ctrs, nbItems, nbClusters);
- vector_push(clusters[affectation], i);
- }
-
- // copy current centers to old centers
- for (uint32_t j = 0; j < nbClusters; j++)
- oldCtrs[j] = ctrs[j];
-
- // recompute centers
- lastDistor = distor;
- assign_centers(nbClusters, clusters, dissimilarities, nbItems, ctrs, &distor);
- }
- while (distor < lastDistor && kounter++ < maxNbIter);
-
- if (distor < bestDistor)
- {
- // copy current clusters into bestClusts
- for (uint32_t j = 0; j < nbClusters; j++)
- {
- vector_clear(bestClusts[j]);
- for (uint32_t i = 0; i < vector_size(clusters[j]); i++)
- {
- uint32_t index;
- vector_get(clusters[j], i, index);
- vector_push(bestClusts[j], index);
- }
- }
- bestDistor = distor;
- }
- }
-
- // Assign centers from bestClusts
- assign_centers(nbClusters, bestClusts, dissimilarities, nbItems, ctrs, &distor);
-
- free(oldCtrs);
- for (uint32_t j = 0; j < nbClusters; j++)
- {
- vector_destroy(clusters[j]);
- vector_destroy(bestClusts[j]);
- }
- free(clusters);
- free(bestClusts);
-}
projects / epclust.git / blobdiff