backup (very) old stuff

dataset.m

// datasets of 2D points:
// http://people.sc.fsu.edu/~jburkardt/datasets/spaeth/spaeth.html
// tr '\n' ';'

// small clouds are consist from 2d points
points = [41 45; 39 44; 42 43; 44 43; 10 42; 38 42; 8 41; 41 41; 
	  13 40; 45 40; 7 39; 38 39; 42 39; 9 38; 12 38; 19 38; 
	  25 38; 6 37; 13 35; 9 34; 12 34; 32 27; 26 25; 39 24;  
	  34 23; 37 23; 22 22; 38 21; 35 20; 31 18; 26 16; 38 13; 
	  29 11; 34 11; 37 10; 40 9; 42  9];

points2 = [16 50;37 50;40 50;20 49;37 48;34 46;38 45;
	   21 44;35 43;24 42;32 42;32 39;24 38;
	   30 38;27 36;30 35;25 34;28 34;27 32; 
	   1 31;29 31;25 30; 3 29;23 29;29 29;24 28;
	   2 26; 6 26 ;19 26;22 26;32 26;10 25;20 24;
	   5 22;12 22;16 22;34 22; 9 20;15 20;
	   13 19;18 19;15 18;13 17;22 17;11 15;
	   15 15;18 15;24 15;36 15;39 14;12 13;10 12;
	   15 12;19 12;24 12;29 12;37 11; 9 10; 7 9;
	   22  9;26  9;41  9;29  8;34  8;39  8; 
	   6  7; 4  5;26  5;37  5;42  5;44  5;32 4;
	   2  3;37  2;47  2; 1  1;30  1;41  1;44  1;46  1];

points3 = [38 50;40 50;44 50;46 50;49 50;42 49;
	      45 49;47 49;14 48;44 48;48 48;50 48;29 47;
	      43 47;46 47;49 47; 7 45;41 45;48 45;50 45;
	      45 44;38 43;42 43;47 43;50 43;44 41;48 41;
	      33 40;47 40;50 40;42 39;10 38;46 38;21 36; 
	      2 35;34 34;42 34;48 27;22 26;36 25; 7 14;
	      26 14;41  7;17  5];

points4 = [24 49;38 45;19 44;23 44;22 43;24 43;
	      26 43;42 43;23 42;40 42;10 41;41 41;44 41;
	      40 40;44 40;20 39; 3 38; 6 37; 8 37;26 37;
	      43 37; 7 36;13 36;38 36; 6 35; 8 35;10 34; 
	      5 33;26 33;38 33; 9 32;31 30;35 30;46 30;
	      11 28;28 28;34 27;36 27;18 26;37 26;35 25;
	      39 25;37 24; 9 23;41 23;31 21;38 21;23 19;
	      50 19;11 18;31 17;39 16;44 15; 7 14;21 13;
	      39 11;12 10;16  9; 5  8;44  8; 9  7;11  7;
	      42  7;10  6;13  6;41  6;43  6; 9  5;11  5;
	      43  5;10  3;27  3;14  1];

# two long thin ellipsoid
points5 = [45 50;22 49;20 48;18 46;21 47;46 47;19 44;
	      41 44;15 42;17 42;19 42;24 41;14 39;17 39;
	      23 39;42 39;21 38;38 38;14 37;12 36;21 36;
	      19 35;34 35;13 34;11 33;19 33;39 33;17 32;
	      35 32; 9 31;15 30;17 30; 9 29;32 29;15 28;
	      29 28;13 27;35 27;11 25;13 25;32 25;24 24;
	      11 23;27 22; 9 21; 7 18;20 17;28 17; 6 16;
	      23 16;18 12;22 11;17  7;12  6; 9  1];


points6 =  [25 42;23 40;25 40;21 39;22 37;18 36;
	       46 35;16 34;19 34;43 34;41 32;14 31;16 31;
	       38 31;34 29;37 29;16 28;36 28;12 27;14 27;
	       31 27;11 25;28 25;30 24;12 23; 9 22;22 22;
	       25 22; 8 20;19 20;22 20;19 18;16 17;14 16;
	       11 15; 6 13; 9 12;12  8;17  8; 9  7;21  7;
	       25  7;15  6;29  6;19  5;24  5;28  4;32  4;
	       37  3;41  3];

# three long thin ellipsoids
points7 = [11 50;16 48;23 47;13 47; 18 43;26 43;17 39;	
	27 39;34 39;23 38;31 35;22 34;37 34;42 34;
	27 33;34 32;37 31;40 31;24 29;39 29;32 28;
	35 28;37 27;26 26;29 25;43 25;33 24;23 23;
	27 23;38 23;19 22;29 22;13 21;25 20;34 20;
	19 19;45 19;15 18;22 18; 7 17;10 16;18 16;
	38 16;45 16;15 13;44 13; 3 12; 5 12;11 12;
	49 11; 1  9; 5  8];

kmeans.hs

function centroids = center(points, clustering, clusterCount)
  centroids = zeros(clusterCount, size(points, 2));
  inClusterCnt = zeros(clusterCount, 1);

  for i = 1:size(points, 1)
    ix = clustering(i);
    centroids(ix, :) = centroids(ix, :) + points(i, :);
    inClusterCnt(ix)++;
  end
  
  for i = 1:size(centroids, 1)
    if (inClusterCnt(i) != 0)
      centroids(i, :) = centroids(i, :) / inClusterCnt(i);
  end
end

function clustering = cluster(points, centroids)
  clustering = zeros(1, size(points, 1));
  replicator = ones(size(centroids, 1), 1);

  for i = 1:size(points, 1)
    replicPoint = replicator * points(i, :);
    distances   = distance(replicPoint, centroids);
    [minDist, clustIx] = min(distances);
    clustering(i) = clustIx;
  end
end

function dist = distance(a, b)
  c = a - b;
  dist = sum(c .* c, 2);
end


function [clustering, centroids] = kmeans(points, initCentroids)
  centroids = initCentroids;
  clusterCount = size(initCentroids, 1);

  do
    predCentroids = centroids;

    clustering = cluster(points, centroids);
    centroids  = center(points, clustering, clusterCount);
  until (centroids == predCentroids)

endfunction


function [clustering, centroids, error] = kmeansBest(points, clusterCount, iterationCount)
  clustering = zeros(size(points, 1), 1);
  centroids  = zeros(clusterCount, size(points, 2));
  error = inf;

  for i = 1:iterationCount
    indices = unidrnd(size(points, 1), 1, clusterCount);
    [nclust, ncent] = kmeans(points, points(indices, :));
    nerror = wcss(points, nclust, ncent);

    if (nerror < error)
      error = nerror;
      clustering = nclust;
      centroids  = ncent;
    end
  end  
end


function plotClusters(points, clustering, centroids)
  clf();

  for i = 1:size(points, 1)
    hold on;
    plot(points(i, 1), points(i, 2), sprintf("*%i", clustering(i)));
  end
  
  for i = 1:size(centroids, 1)
    hold on;
    plot(centroids(i, 1), centroids(i, 2), 
	 sprintf("*%i", i), 
	 "markersize", 50);
  end

  hold off;
end


function  test(points, clusterCount)
  [clustering, centroids, error] = kmeansBest(points, 
					      clusterCount, 10);
  plotClusters(points, clustering, centroids);
  printf("error = %i \n", error);
  fread();
end


# within-cluster sum of squares
function error = wcss(points, clustering, centroids)
  error = 0;

  for i = 1:size(points, 1)
    error += distance(points(i, :), centroids(clustering(i), :));
  end
end

sofm.m

function animateSofm(points, clusterCount, 
                     epochs, radius, lrate, time)
#  fm =  init_sofm(points, clusterCount);
  fm = unidrnd(1000, clusterCount, size(points, 2)) ./ 20;
 
  for i = 1:epochs
    fm = train(fm, points, radius, lrate);
    hold on;
    plotSofm(points, fm);
    sleep(time);
  end
end

function category = categorize(sofm, point)
  replicator = ones(size(sofm, 1), 1);
  distances   = distance(replicator * point, sofm);
  [_ category] = min(distances);
end

function fmr = correct(fm, centroid, point, radius, lrate)
  for j = 1:size(fm, 1)
    #if distance(centroid, fm(j, :)) < radius
    dist = sqrt(sum((centroid - fm(j, :)) .^ 2));
    near = exp(- (1 / radius) * dist);
    fm(j, :) += lrate * (point - fm(j, :)) * near;
   # end
  end
  fmr = fm;
end

function dist = distance(a, b)
  c = a - b;
  dist = sum(c .* c, 2);
end

function fm = init_sofm(points, clusterCount)
  idxs = unidrnd(size(points, 1), 1, clusterCount);
  fm = points(idxs, :);
end

function plotSofm(points, centroids)
  clf();
  plot(points(:, 1), points(:, 2), "*");

  hold on;
  plot(centroids(:, 1), centroids(:, 2), "*2",
       "markersize", 50);
end

function fm = sofm(points, radius, lrate, 
                   clusterCount, epochs)
  fm = init_sofm(points, clusterCount);

  for i = 1:epochs
    fm = train(fm, points, radius, lrate);
  end
end

function fm = animateSofm(points, clusterCount, 
                          epochs, radius, lrate)
  fm = init_sofm(points, clusterCount)
  for i = 1:epochs
    fm = train(fm, points, radius, lrate);
    plotSofm(points, fm)
  end
end

function fm = train(init_fm, points, radius, lrate)
  fm = init_fm;

  for i = 1:size(points, 1)
    point = points(i, :);
    category = categorize(fm, point);
    fm = correct(fm, fm(category, :), point, radius, lrate);
  end
end