How can I plot this one-dimensional Self Organizing (Kohonen) Map?

Question

Mathematica used to have a package called NeuralNetworks that is not currently working. With that package, you could create a Kohonen or Self Organizing Map (https://en.wikipedia.org/wiki/Self-organizing_map) which I want to use as a nonlinear alternative to PCA.

The points of the cloud are generated like this:

angles = Table[Pi/200. i, {i, 0, 99}];

x = Map[{Cos[#], Sin[#]} RandomReal[{0.9, 1.1}] &, angles];

p1 = ListPlot[x, PlotRange -> {{0, 1}, {0, 1}}, 
  PlotStyle -> {PointSize[0.01], RGBColor[0, 0, 0]}, AspectRatio -> 1]

But I need to create the nodes (red dots) like in this plot:

This could be donde with the next code (not available anymore):

unsup = IntializeUnsupervisedNet[x, 6];
{unsup, fitrecord} = 
   UnsupervisedNetFit[x, unsup, 100, ReportFrequency -> 1]; // Quiet

And plotting:

p3=ListPlot[unsup[[1]],PlotRange {{0,1},{0,1}},
PlotStyle {PointSize[0.02],RGBColor[1,0,0]}];
Show[{p3,p1},AspectRatio 1,PlotRange All]

Documentation for the functions is available at http://media.wolfram.com/documents/NeuralNetworksDocumentation.pdf, chapter 10 Unsupervised Networks.

There is also an example in the Wikipedia page:

Answer 1

I implemented the algorithm from the Wikipedia page. I hope it can be of help. I decided to write it down as a sort of guide, but read on Wikipedia for the details and for ways to further tune it.

A SOM has two types of data, weights and units. The weights are data points, we use the following:

img = Import["https://i.sstatic.net/eoAEt.png"] // Binarize;
weights = Join[
   N@Rescale@PixelValuePositions[img, 1],
   RandomReal[1, {500, 2}]
   ];
ListPlot[weights]

The units are objects that have a weight and a position in a grid. The weight is in the same space as the weights previously defined, the grid describes the spatial relationships between the units. We can initialize units by picking weights randomly from the general area where the data points are:

units = N@MapThread[unit, {Subdivide[39], RandomReal[{0.2, 0.8}, {40, 2}]}];

We proceed now to update the units by feeding them weights from the dataset one by one. One update consists of:

1. Identifying the unit whose weight is closest to the input weight. This the so-called best matching unit, the BMU.
2. Identifying units that are within a pre-defined range of the BMU on the grid (not in the weight space), the so-called activated units.
3. Applying the update formula.

We start by defining the following helper functions:

weightDistance[unit[_, w_], weight_] := EuclideanDistance[w, weight]
findBMU[units_, weight_] := First@MinimalBy[units, weightDistance[#, weight] &]
gridDistance[unit[pos1_, _], unit[pos2_, _]] := EuclideanDistance[pos1, pos2]
findActivatedUnits[units_, bmu_, range_] := Transpose@Nearest[
   units -> {"Index", "Element"}, bmu, {All, range},
   DistanceFunction -> gridDistance
   ]
update[unit[pos_, w_], weight_, eta_] := unit[pos, w + eta (weight - w)]
getWeights[units_] := Last /@ units

The following is the iteration with eta the initial learning rate, nIterations the number of iterations, neighborhoodRadius the range that decides how close units must be to the BMU to be affected, and nUnits the number of units:

eta = 1;
neighborhoodRadius = 0.1;
nIterations = 20000;
nUnits = 50;

units = N@MapThread[unit, {Subdivide[nUnits - 1], RandomReal[{0.45, 0.55}, {nUnits, 2}]}];
Do[
  weight = RandomChoice[weights];
  bmu = findBMU[units, weight];
  {indices, activatedUnits} = findActivatedUnits[units, bmu, neighborhoodRadius];
  units = MapAt[update[#, weight, eta - s eta/nIterations] &, units, List /@ indices];
  ,
  {s, nIterations}
  ];

This produces the following output, where the dark yellow/orange line represents the SOM units:

ListPlot[{
  weights,
  getWeights[units]
  }, Joined -> {False, True}]

This looks pretty good, but note that I had to play with the parameters to get it to work. I needed a lot of iterations, a sufficient number of units, and the neighborhood radius needed to be small enough to make the updates localized. If the parameters are off it will probably still return a solution that makes sense in that it has spread out to cover the points pretty well, but it won't have found the optimal solution, which is the S shape.