Packs search

Table of Contents

Packs (or clusters)

Let's consider the problem of recognizing data clusters. Imagine a space with a metric that can correspond to the similarity of two objects: 

|   |   | | ||| ||  |    | | ||| |  |   |
|   |   | | ||| ||  |    | | ||| |  |   |
        :........:       :.......:
        :        :       :       :
        :        :       :       :
      left     right   left    right

For example, it could be a clusterization of colors, texts, etc. - everything having metric. And we are interested in finding of (left, right) pairs (2 points) of input points where a pack starts, ends.

The idea of the algorithm

The idea of the algorithm is based on a calculation of distance between adjacent points, then to process these distances (called "differences"): 

                1             2                                      3
                                                                +---------+  Result: segments
              +----+   +-------------+                          |rise/fall|     of indexes
1,2,1,5... -->|diff|-->|   Riemann   |-->1.33,1.33,1.33,5.2..-->|  front  |-->(0,2),(3,..)..
              +----+   |approximation|                          | detect  |
                       +-------------+                          +---------+

In the algorithm below units are implemented as functions:

  1. calc_diffs()
  2. approximate_by_rectangles()
  3. find_clusters_centers()

Important nuance is that approximation_by_rectangles takes an argument tolerance which means what values to treat as "still-in-the-same-rectangle" (the function itself takes into accound the min/max values of points group).

find_clusters_centers also have 2 additional arguments: ratio and nesting. To be honest, nesting is redundant and such post-processing can be done outside the function (it's idea just to increment indexes by this value nesting). But ratio is important: it tells when change of "frequency" leads to pack start/end bounds detection. I.e, if the distance is more than this, then treat it as a bound of a pack.

Octave implementation

This is GNU Octave code of the algorithm, and at the end you can see example output of it (found indexes) - also they are shown as thick horizontal segments on the 0x axis (ie, they are (1, 9), (12, 21) if to numerate from 0): 

function algorithm_demo(odir)
  file = fullfile(odir, 'Packs_Search.svg');
  ps = [1,7,9,10,11,13,15,16,18,19,25,30,34,36,37,38,39,40,42,43,45,46,49,52,58];
  # ps = [1,2,10,20,30,40,50,60,70,80,90,100,101]
  # ps = [1,2,1,1,1,10,20,30,40,50,60,70,80,90,100,101]
  display('INPUT');
  disp(ps);
  # clusters - matrix where a raw is indexes of input values where we find packs:
  clusters = calc_clusters(ps);
  disp('CLUST');
  disp(clusters);
  % Create smaller figure
  close all;
  reset(groot);
  fig = figure('visible', 'off', 'position', [0, 0, 350, 350]); % width, height in pixels
  hold on;
  # Draw thick cluster segments
  for k = 1:size(clusters,1)
    x1i = clusters(k,1); # left index of a pack (in input)
    x2i = clusters(k,2); # right index of a pack (in input)
    x1 = ps(x1i);        # get a value of the left bound to draw it on 0x axis (as input values are drawn)
    x2 = ps(x2i);        # get a value of the left bound to draw it on 0x axis (as input values are drawn)
    line([x1, x2],[0, 0], ...
         'Color', [1, 0, 1], ...  # magenta (or fuchsia in Web/CSS terms)
         'LineWidth',8);
  endfor
  # Draw cal sticks:
  for x = ps
    line([x, x],[0, 1], ...
         'Color', [1, 0, 1], ...
         'LineWidth',1.5);
  endfor
  title('Packs detection');
  # Export as SVG
  print(fig, file, '-dsvg', '-painters');
  close(fig);
endfunction

function [diffs] = calc_diffs(ps, absolute=false)
  diffs = diff(ps);
  if (absolute)
    diffs = abs(diffs);
  endif
endfunction

function [ratio] = calc_ratio(x, y)
  abs_x = abs(x);
  abs_y = abs(y);
  res_sign = sign(x) * sign(y);
  if (abs_x == 0 && abs_y == 0)
    ratio = NaN;
  elseif (abs_x > abs_y)
    if (abs_y == 0)
      if (res_sign < 0)
        ratio = -Inf;
      else
        ratio = Inf;
      endif
    else
      ratio = x / y;
    endif
  elseif (abs_x < abs_y)
    if (abs_x == 0)
      if (res_sign < 0)
        ratio = -Inf;
      else
        ratio = Inf;
      endif
    else
      ratio = y/x;
    endif
  else
    if (abs_x == 0)
      ratio = NaN;
    else
      ratio = res_sign;
    endif
  endif
endfunction

function [needed] = new_rect_is_needed(slice, p, tolerance)
  if (isempty(slice))
    error("slice is empty");
  else
    slice_max = max(slice);
    slice_min = min(slice);
    needed = calc_ratio(p, slice_max) > tolerance || ...
             calc_ratio(p, slice_min) > tolerance;
  endif
endfunction

function [res_ps, res_ranges] = approximate_by_rectangles(ps, tolerance=1.5)
  res_ps = [];
  res_ranges = [];
  slice = [];
  for p = ps
    if (isempty(slice))
      slice = [slice, p];
    else
      if (new_rect_is_needed(slice, p, tolerance))
        slice_avg = mean(slice);
        res_ps = [res_ps, slice_avg * ones(1, length(slice))];
        slice = [p];
      else
        slice = [slice, p];
      endif
    endif
  endfor
  if (!isempty(slice))
    slice_avg = mean(slice);
    res_ps = [res_ps, slice_avg * ones(1, length(slice))];
  endif
  pr = NaN;
  for i = 1:length(res_ps)
    r = res_ps(i);
    if (i == 1)
      res_ranges = [res_ranges, i];
    elseif (i != 1 && pr != r)
      res_ranges = [res_ranges, i];
    endif
    pr = r;
  endfor
endfunction

function [res] = find_clusters_centers(ranges, isovalues, ratio=2, nesting=0)
  prev_isovalue = NaN;
  res = [];
  max_irange = length(ranges) - 1;
  for irange = 0:max_irange
    isovalue = isovalues(ranges(irange + 1) + 1);
    if (!isnan(prev_isovalue))
      if (isovalue > prev_isovalue && ...
          calc_ratio(isovalue, prev_isovalue) >= ratio)
        % left___/‾‾
        right = ranges(irange + 1);
        left = ranges(max(irange - 1, 0) + 1);
      elseif (isovalue < prev_isovalue && ...
              calc_ratio(prev_isovalue, isovalue) >= ratio)
        % ‾‾\___right
        left = ranges(irange + 1);
        right = ranges(min(irange + 1, max_irange) + 1);
      else
        prev_isovalue = isovalue;
        continue;
      endif
      if (right == left)
        right = right + 1;
      endif
      left = left + nesting;
      right = right + nesting;
      if (isempty(res) || res(end,1) != left || res(end,2) != right)
        res(end + 1, :) = [left, right];
      endif
    endif
    prev_isovalue = isovalue;
  endfor
endfunction

function [clusters] = calc_clusters(ps)
  diffs = calc_diffs(ps);
  [isovalues, ranges] = approximate_by_rectangles(diffs, 2);
  clusters = find_clusters_centers(ranges, isovalues, 0);
endfunction

algorithm_demo(odir);

INPUT
 Columns 1 through 16:

    1    7    9   10   11   13   15   16   18   19   25   30   34   36   37   38

 Columns 17 through 25:

   39   40   42   43   45   46   49   52   58
CLUST
    2   10
   13   22

Gnuplot Produced by GNUPLOT 6.0 patchlevel 4 0 0.2 0.4 0.6 0.8 1 0 10 20 30 40 50 60 gnuplot_plot_1a gnuplot_plot_2a gnuplot_plot_3a gnuplot_plot_4a gnuplot_plot_5a gnuplot_plot_6a gnuplot_plot_7a gnuplot_plot_8a gnuplot_plot_9a gnuplot_plot_10a gnuplot_plot_11a gnuplot_plot_12a gnuplot_plot_13a gnuplot_plot_14a gnuplot_plot_15a gnuplot_plot_16a gnuplot_plot_17a gnuplot_plot_18a gnuplot_plot_19a gnuplot_plot_20a gnuplot_plot_21a gnuplot_plot_22a gnuplot_plot_23a gnuplot_plot_24a gnuplot_plot_25a gnuplot_plot_26a gnuplot_plot_27a Packs detection

Python implementation

The same algorithm implemented in Python: 

'''Finds clusters (concentrations) of points - their centers, ranges
'''
def calc_avg(ps:list[int]) -> float:
    return sum(ps)/len(ps)


def calc_diffs(ps:list[int], absolute:bool=False) -> list[int]:
    diffs_num = max(len(ps) - 1, 0)
    if diffs_num == 0:
        return []
    else:
        diffs:list[int] = [0] * diffs_num
        i = 0
        while i != len(ps) - 1:
            diff = ps[i + 1] - ps[i]
            if absolute: diff = abs(diff)
            diffs[i] = diff
            i += 1
        print('DIFFR', diffs) #, calc_avg(diffs), calc_rmsd(diffs, normalize_with_iqr=False))
        return diffs


class RangeBounds(list): ...
class Isolines(list): ...


def sign(x):
    return -1 if x < 0 else 1


def calc_ratio(x:int|float, y:int|float) -> float:
    abs_x = abs(x)
    abs_y = abs(y)
    res_sign = sign(x) * sign(y)
    if abs_x == 0 and abs_y == 0:
        return float('nan')
    elif abs_x > abs_y:
        if abs_y == 0:
            if res_sign < 0:
                return float('-inf')
            else:
                return float('inf')
        else:
            return x / y
    elif abs_x < abs_y:
        if abs_x == 0:
            if res_sign < 0:
                return float('-inf')
            else:
                return float('inf')
        else:
            return y / x
    else: #if abs_x == abs_y:
        if abs_x == 0:
            return float('nan')
        else:
            return res_sign


def approximate_by_rectangles(ps:list[int], tolerance:float=1.5) -> tuple[RangeBounds,Isolines]:
    def new_rect_is_needed(slice, p):
        # new rect if current point p is too high/low comparing to min/max of current slice
        if slice:
            slice_max = max(slice)
            slice_min = min(slice)
            return calc_ratio(p, slice_max) > tolerance or calc_ratio(p, slice_min) > tolerance
        raise ValueError('slice is empty')
    res_ps:Isolines = Isolines()
    res_ranges:RangeBounds = RangeBounds()
    slice = []
    for p in ps:
        if not slice:
            slice.append(p)
        else:
            if new_rect_is_needed(slice, p):
                slice_avg = calc_avg(slice)
                res_ps.extend(len(slice) * [slice_avg])
                slice = [p]
            else:
                slice.append(p)
    if slice:
        slice_avg = calc_avg(slice)
        res_ps.extend(len(slice) * [slice_avg])
    pr = None
    for i, r in enumerate(res_ps):
        if i == 0:
            res_ranges.append(i)
        elif i != 0 and pr != r:
            res_ranges.append(i)
        pr = r
    print('APPRX', res_ps, res_ranges)
    return (res_ranges, res_ps)


def find_clusters_centers(
        *,
        ranges:RangeBounds,
        isovalues:Isolines,
        ratio:float=2,
        nesting:int=0) -> list[tuple[int,int]]:
    '''Finds high concentrations of points based on the approximation of the distance
    between neighboring points by rectangles. The input from approximation is ranges and isovalues
    '''
    prev_isovalue = None
    res = []
    max_irange = len(ranges) - 1
    for irange in range(0, max_irange + 1):
        isovalue = isovalues[ranges[irange]]
        if prev_isovalue is not None:
            if isovalue > prev_isovalue and calc_ratio(isovalue, prev_isovalue) >= ratio:
                # left___/‾‾
                right = ranges[irange]
                left = ranges[max(irange - 1, 0)]
            elif isovalue < prev_isovalue and calc_ratio(prev_isovalue, isovalue) >= ratio:
                # ‾‾\___right
                left = ranges[irange]
                right = ranges[min(irange + 1, max_irange)]
            else:
                prev_isovalue = isovalue
                continue
            if right == left:
                right += 1
            left += nesting
            right += nesting
            if not res or res[-1] != (left, right):
                # ‾‾\___, ___/‾‾ such cases would be reported twice w/o this "if" checking that
                # if was not reported in `res` in the previous time:
                res.append((left, right))
        prev_isovalue = isovalue
    return res


def calc_clusters(ps:list[int]):
    diffs = calc_diffs(ps)
    ranges, isovalues = approximate_by_rectangles(diffs, tolerance=2)
    clusters = find_clusters_centers(ranges=ranges, isovalues=isovalues, nesting=0)
    return clusters


ps = [1,7,9,10,11,13,15,16,18,19,25,30,34,36,37,38,39,40,42,43,45,46,49,52,58]
# ps = [1,2,10,20,30,40,50,60,70,80,90,100,101]
# ps = [1,2,1,1,1,10,20,30,40,50,60,70,80,90,100,101]
print('INPUT', ps)
clusters = calc_clusters(ps)
print('CLUST', clusters)

INPUT [1, 7, 9, 10, 11, 13, 15, 16, 18, 19, 25, 30, 34, 36, 37, 38, 39, 40, 42, 43, 45, 46, 49, 52, 58]
DIFFR [6, 2, 1, 1, 2, 2, 1, 2, 1, 6, 5, 4, 2, 1, 1, 1, 1, 2, 1, 2, 1, 3, 3, 6]
APPRX [6.0, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 5.0, 5.0, 5.0, 1.3333333333333333, 1.3333333333333333, 1.3333333333333333, 1.3333333333333333, 1.3333333333333333, 1.3333333333333333, 1.3333333333333333, 1.3333333333333333, 1.3333333333333333, 4.0, 4.0, 4.0] [0, 1, 9, 12, 21]
CLUST [(1, 9), (12, 21)]

Sure, it can be achieved in other ways too (I think about FFT, for example). This algorithm can be used to form clusters of objects that look similar, that are distributed by some parameter.