Kohonen model for "hearing"

A sound source is located in a region G bounded by  y>0, y < x^2, |x| < 1, as shown below:

Plot[1 - x^2,  {x, -1, 1}]

at locations (a,0), and (-a,0), two microphone are located.   The microphone signals are fed into two logarithmic amplifier, whose output signals (two scalars) {v1,v2} serve as input for the net work.  The signal v1 and v2 is related to the location {x,y}   by the function v defined below.  First we at the distance of microphones from the center a as 0.25.

a  =  0.25

v[{x_, y_}]  :=  { - Log[ (x - a)^2  +  y^2],  - 
          Log[(x + a)^2  +  y^2] }

I.e., the vector v = {v1, v2} is proportional to the distances from the source to the microphones.

We define a function to get sample in such a way so that a point (x,y)is distributed in the domain G uniformly. This can be done by first generating points uniformly in the rectangle -1 < x < 1, and   0 < y < 0, and simply omitting the points outside the domain.

usample[]  :=  
      Module[ {x  =  0,  y  =  2.0}, 
While[
          y  >  1 -  x^2, 
x  =  
            2.0*Random[]  -  1.0; 
y  =  
            Random[];
]; 
{x, 
          y}
]

Let's test output of the usample.

ListPlot[Table[usample[],  {10000}]]

The sample that is fed into neural network is

ListPlot[Table[v[usample[]], {10000}]]

Let the neural network lattice be L by L.  We set the size here:

L  =  15

We initialize the synaptic weight w[[i,j,k]] at random.  Here i and j varies over 1 to L, and k = 1 and 2.  That is, neuron {i,j} has two inputs from w[[i,j,1]] and w[[i,j,2]].

w  =  
        Table[ usample[],  {i,  1,  L},  {j,  1,  L}];

Some explanation of the notation.  w refers to the hole matrx with index i,j, k.   w[[1,1,1]] is an element for the site {1,1} first compoment.  w[[1]] refers to the first row of a table each of which is a vector of two component.  w[[i,j]] is the synaptic weight associated with [i,j]-th nueron.

We define a function that gives a map from v to r = {i,j}. i.e., given the input v ={v1,v2}, which neuron fire.  This is obtained by the condition that | v - w | is a minimum.  The function distancesquare computes the Euclidean distance squared.  The function fire takes a input vector v and the synaptic matrix w and computes the smallest distance small and the index {i,j} corresponding to the smallest distance.

distancesquare[u_,  
        v_]  :=  ( u  -  v)  .  (u  -  v)

fire[v_List, 
        w_List] := Module[{i, j, small,  ind, 
          measure}, 
measure  =  
          Map[ distancesquare[v,  #] &,  
            w,  {2}]; 
    small  =  
          Infinity; 
    ind  =  {1, 
            1}; 
For[i  =  1,  
          i  <=  Length[w],  ++i, 
For[j  =  1,  
            j  <=  Length[w],  ++j, 
If[ 
              small >  
                measure[[i, j]],   {small,  
                  ind }  =  {measure[[i, j]],  {i, 
                    j}}]
]
]; 
{small,  ind}
]

We know define the step that updates the synaptic strength w according to Kohonen's proposal.  That is
w(r)  <- w(r) + epsilon h(r,r') (v - w(r)) .  
h(r,r') is the Gaussian function, for which we first define

h[r_,  s_,  sigma_]  :=  
      Exp[ - distancesquare[r, s]/(2 sigma^2)]

adaptation[w_,  v_,  s_,  epsilon_,  sigma_]  :=  
      Module[ {i, j, k, L  =  Length[w]}, 
 
        Table[ w[[i, j, k]]   +  
            epsilon * 
              h[{i, j}, s,  sigma] * ( 
                v[[k]]  -  w[[i, j, k]]),  {i,  1,  
            L},  {j,  1,  L},  {k,  1, 
            2}]
]

We now can put all the functions defined ealier together to get the main function.  In the main function, we generate a point unformly in the domain G.  The samples are fed to the network for the adaption.  After sufficient samples, we plot the w to inspect the results.   In the program below, we have commented out the transformation from (x,y) to (v1,v2).

main[nmax_,  ϵ0_,  σ0_]  :=  
      Module[{vdata, ϵ, σ, small,  ind, 
          n}, 
Do[ 
 (* vdata  =  
                  v[usample[]]; *) 
vdata  =  
              usample[]; 
ϵ  =  ϵ0 * 
                Exp[ - n/
                    1000.0]; 
σ  =  σ0 * (0.01)^(n/10000.0); 
{small,  ind}  =  
              fire[vdata,  w]; 
w  =  
              adaptation[w, vdata, 
                ind, ϵ,  σ], 
{n,  
              nmax}];
  
]

draw[w_]  :=  
      Show[Graphics[Map[Line,  Join[w,  Transpose[w]]]],  
        Axes  ->  Automatic]

Intial random w

draw[w]

main[5,  0.5, 15.0]; draw[w]

main[100,  0.5, 15.0]; draw[w]

main[1000,  0.5, 15.0]; draw[w]

main[10000,  0.1, 2.0]; draw[w]