Consensus in CxSOM

The letter c in cxsom stands for “consensus”, since cxsom enable to define “consensus driven self organizing maps”. Consensus is, at time \(t\), a relaxation process that runs until all the map “agree” on values which depend ones on the others.

A basic relaxation process

Consensus between two maps

Let us define two maps “A” and “B”. For the sake of visualization, 1D maps are considered. They are self-organizing maps, with the specificity that the input of one map is the BMU of the other. So the weights of the map represent prototypes of the BMU values seen in the other.

Let us ommit the learning process, and start with arbitrary weights in the two maps. If we set the two BMUs (in map “A” and map “B”) as 0, we get the following unstabe update:

The idea underlying cxsom BMU updates is, as opposed to classical SOMs, to move the \(\mathrm{BMU}^t\) position toward the best matching position, with a small steps (we use 0.005, a very small one, for the sake of illustration), rather than directly setting it to the best matching position. Such a process often reaches stability. The whole relaxation aims at computing the \(\mathrm{BMU}^t\), the relaxation cycle occur inside a timestep \(t\).

The `cxsom` implementation

The relaxation is done internally, when computing time steps. The following code implementrs the rules corresponding to related “A” and “B” maps.

// This is the relax.cpp file, describing the relaxation example updating rules.

#include <cxsom-rules.hpp>
using namespace cxsom::rules;
context* cxsom::rules::ctx = nullptr;

#define CACHE          2 // nb of instances kept in processor memory
#define WTRACE         1 // Size of the weight buffer
#define ITRACE      1000 // Size of the BMU buffer
#define FORGET         0 // Size of the 0-sized buffers.
#define OPENED      true // We keep the files opened at processor side.

std::string current(const std::string& map) {return std::string("/") + map + "/";}
std::string other  (const std::string& map) {if(map == "A") return current("B"); return current("A");}

int main(int argc, char* argv[]) {
  context c(argc, argv);

  kwd::parameters p_main, p_relax;
  p_main  | kwd::use("walltime", -1), kwd::use("epsilon", 1e-3);
  p_relax | p_main, kwd::use("delta", .005), kwd::use("deadline", 500);
  {
    timeline t("weights");
    for(std::string map : {"A", "B"}) {// weights/A/W and weights/B/W definitions.
      name_space ns(map);
      kwd::type("W", "Map1D<Pos1D>=500", CACHE, WTRACE, OPENED);
    }
  }
  {
    timeline t("inits");
    for(std::string map : {"A", "B"}) { // inits/A/BMU and inits/B/BMU definitions.
      name_space ns(map);
      kwd::type("BMU", "Pos1D", CACHE, ITRACE, OPENED);
      
      "BMU" << fx::random()                                                                                | kwd::use("walltime", ITRACE-1);
    }
  }  
  {
    timeline t("relax");
    kwd::type("Cvg", "Scalar", CACHE, ITRACE, OPENED); // For counting relaxation steps.
    for(std::string map : {"A", "B"}) { // inits/A/BMU and inits/Y/BMU definitions.
      name_space ns(map);
      kwd::type("A",   "Map1D<Scalar>=500", CACHE, FORGET, OPENED);
      kwd::type("BMU", "Pos1D",             CACHE, ITRACE, OPENED);
      
      "A"   << fx::match_gaussian(other(map) + "BMU", kwd::at(kwd::var("weights", current(map) + "W"), 0)) | p_relax, kwd::use("sigma", .16);     
      "BMU" <= fx::copy(kwd::var("inits", "BMU"))                                                          | p_relax;
      "BMU" << fx::toward_argmax("A", "BMU")                                                               | p_relax; 
    }
    "Cvg" << fx::converge({"/A/BMU", "/B/BMU"})                                                            | p_relax;                                                                   
  }
  return 0;
}

It describes the following updates:

The relaxation first steps of updating rules

Run `cxsom` to make statistics about that relaxation

On a terminal, let us restart a processor (see the Basic SOM section)… or use the one we have already started previously if there is any… let us consider this in the following. We have to clear the previous data (can be dangerous) and to clear the updating rules that may be still loaded on our processor.

mylogin@mymachine:~$ make cxsom-is-processor-running
mylogin@mymachine:~$ make cxsom-clear-processor
mylogin@mymachine:~$ make cxsom-clear-rootdir

Let us send the rules we have defined

mylogin@mymachine:~$ g++ -o relax -std=c++17 relax.cpp $(pkg-config --cflags --libs cxsom-rules)
mylogin@mymachine:~$ ./relax send localhost 10000

Then run that script to get the statistics

mylogin@mymachine:~$ python3 relax-stats.py root_dir/ localhost 10000

# This is the relax-stats.py script

import matplotlib.pyplot as plt
import numpy as np
import pycxsom as cx
import sys


if len(sys.argv) < 2:
    print('Usage: {} <root_dir> [<hostname> <port>]'.format(sys.argv[0]))
    sys.exit(0)

root_dir = sys.argv[1]

if len(sys.argv) > 2:
    hostname = sys.argv[2]
    port     = int(sys.argv[3])
else:
    hostname = None

MAP_SIZE = 500
pos      = np.linspace(0, 1, MAP_SIZE)
WA       = .1 + .5*pos + .15*np.sin(13*pos)
WB       = 1 - WA

with cx.variable.Realize(cx.variable.path_from(root_dir, 'weights', 'A/W')) as W:
    r = W.time_range()
    if r is None: # W[0] not initialized yet
        W += WA
        
with cx.variable.Realize(cx.variable.path_from(root_dir, 'weights', 'B/W')) as W:
    r = W.time_range()
    if r is None: # W[0] not initialized yet
        W += WB

BMU_A_path = cx.variable.path_from(root_dir, 'relax', 'A/BMU')
BMU_B_path = cx.variable.path_from(root_dir, 'relax', 'B/BMU')
CVG_path   = cx.variable.path_from(root_dir, 'relax', 'Cvg'  )
    
if hostname is not None:
    print('Sending "ping" to {}:{}'.format(hostname, port))
    if cx.client.ping(hostname, port):
        print('Something went wrong.')

    ### We wait for the computation to be achieved (checking A/BMU)
    with cx.variable.Realize(BMU_A_path) as BMU:
        last = 0
        r = BMU.time_range()
        print('BMU range {}'.format(r))
        if r is not None:
            last = r[1]
        while last != 999:
            BMU.sync_init()
            BMU.wait_next(sleep_duration=2.0)
            r = BMU.time_range()
            print('BMU range {}'.format(r))
            if r is not None:
                last = r[1]
            
    ### We wait for the computation to be achieved (checking B/BMU)
    with cx.variable.Realize(BMU_B_path) as BMU:
        r = BMU.time_range()
        while r[1] != 999:
            BMU.sync_init()
            BMU.wait_next(sleep_duration=1.0)


### We are ready to plot BMUs
BMU_A = np.array([bmu for t, bmu in cx.variable.data_range_full(BMU_A_path)])
BMU_B = np.array([bmu for t, bmu in cx.variable.data_range_full(BMU_B_path)])
CVG   = np.array([nb  for t, nb  in cx.variable.data_range_full(CVG_path  )])
converged = CVG < 500 # Deadline is 500 for relaxation updates/
print(CVG)
plt.figure(figsize=(5,5))
ax = plt.gca()
ax.set_aspect(1)
ax.set_xlim(-.1, 1.1)
ax.set_ylim(-.1, 1.1)
ax.set_title('BMU A vs BMU B once relaxed')
ax.scatter(BMU_A, BMU_B, alpha=.2)
# We could have displayed BMU_A[converged], BMU_B[converged] but here, there is
# no convergence....
plt.savefig('relaxation-result.png')
print()
print('File "relaxation-result.png" generated')
print()

and you get that plot, once the computation is done… Remind that relaxation has a deadline, so these points are not all fixed points (do we have cycles ?).

A similar script can be used to get this:

Hervé Frezza-Buet,

Consensus in CxSOM

A basic relaxation process

Consensus between two maps

The cxsom implementation

Run cxsom to make statistics about that relaxation

The `cxsom` implementation

Run `cxsom` to make statistics about that relaxation