25.06.2018 2D-Linq и оптимизация цифровых фильтров
 
Этот пост, по-хорошему, заслуживает место во flame.comp, т.к. побудительным мотивом изначально являлся вопрос alex_public, заданный в недрах флейма в том форуме:

Но могу для разнообразия подкинуть ещё один примерчик для понимания. Простейшая, но очень нужная во многих областях задачка: имеем двухмерный массив чисел и надо получить из него новый массив, в котором каждая точка является средним от своих соседей (4-ёх или 8-ми, не суть). Что об этом скажет linq? )))

(https://rsdn.org/forum/flame.comp/7137817.1
Автор: alex_public
Дата: 07.05 20:18
)

Но, поскольку речь идёт об пределах возмжностей кода на C#, я всё же запощу это сюда — в интересах тех, кого тоже волнует вопрос о производительности "математического" кода на дотнете, а перечитывать топики глубиной за сотню постов терпения нет.

Итак, что об этом скажет Linq? Вопрос, понятное дело, с подковыркой — типа, да Linq же — это убогая надстройка над IEnumerable, которую потом и соплями прикрутили к RDBMS. А тут у нас — принципиально другой материал, неперечислимый и, в отличие от, скажем, графов, к перечислениям несводимый.

Лично моими устами об этом Linq скажет следующее:
int[,] FourNeighborAverage(int[,] data) => from d in data select (d[-1, 0] + d[1, 0] + d[0, 1] + d[0, -1]) / 4;


Как это работает?
Ну, для начала у нас есть public static класс Array2d, где определён вот такой вот extension метод:
public static T[,] Select<T>(this T[,] source, Kernel<T> kernel)

То есть выражение d => (d[-1, 0] + d[1, 0] + d[0, 1] + d[0, -1]) / 4 трактуется как Kernel<int>:
public delegate T Kernel<T>(IRelativeAccess2d<T> point);

В свою очередь, интерфейс IRelativeAccess2d<T> определён весьма незатейливо:
public interface IRelativeAccess2d<T>
{
   T this[int dx, int dy] { get; } 
}

Теперь всё становится понятным: в наш Kernel передаётся интерфейс для доступа к "соседним" относительно текущего элементам массива, что позволяет нам описывать выражения типа "среднее от своих четырёх соседей".

Реализовать наш метод Select не очень трудно. Наивная реализация выглядит так:
public static T[,] Select<T>(this T[,] source, Kernel<T> kernel)
{
  var h = source.GetLength(0);
  var w = source.GetLength(1);
  var result = new T[h, w];
  for (int i = 0; i < h; i++)
    for (int j = 0; j < w; j++)
      result[i, j] = kernel(new RelativeArrayAccess2d<T>(source, i, j));
  return result;
}

структура RelativeArrayAccess2d<T> слишком тривиальна, чтобы её приводить — она хранит ссылку на массив source и координаты текущего элемента; при обращении по this[dx, dy] мы прибавляем относительные координаты к текущим и возвращаем элемент.
Этого кода уже достаточно для того, чтобы реализовать копирование массива в стиле linq:
from d in data select d[0, 0];

К сожалению, первый же запуск с ядром усреднения вылетит с ArrayIndexOutOfBoundsException. Логично — мы пробуем обращаться за пределы исходного массива. Суровые программисты на плюсах на такие мелочи внимания не обращают — ну, или требуют от программиста ручного указания размеров "каёмки", в которой фильтр не работает.
Правильное решение этой задачи — довольно сложно: нам придётся учиться подпиливать фильтр на ходу таким образом, чтобы он корректно работал на границах массива (у элемента рядом со "стороной" прямоугольника не 4, а 3 соседа, а у углового элемента — только 2. Строго говоря, именно их и нужно усреднять).
Но пока что мы обойдёмся решением в стиле С++ — ограничением области вывода фильтра.
Сделать это не очень сложно — мы просто один раз вызовем kernel со специальной реализацией IRelativeAccess, которая запоминает, куда идёт обращение:
private sealed class KernelMeasure<T> : IRelativeAccess2d<T>
{
  ...
  T IRelativeAccess2d<T>.this[int x, int y]
  {
    get
    {
      xmin = Math.Min(xmin, x);
      xmax = Math.Max(xmax, x);
      ymin = Math.Min(ymin, y);
      ymax = Math.Max(ymax, y);
      return default(T);
    }
  }
}

Вообще говоря, это будет работать не для всех типов Kernel — например, мы могли использовать Math.Random() для выбора смещения. Но для практически интересных случаев, когда смещения всегда постоянны, этот подход сработает.
Итого, у нас получается более-менее работоспособный код, который позволяет удобно записывать формулы фильтрации:
public static T[,] Select<T>(this T[,] source, Kernel<T> kernel)
{
  var h = source.GetLength(0);
  var w = source.GetLength(1);
  var km = new KernelMeasure<T>();
  kernel(km); // провели измерения
  var t = 0 - km.xmin;
  var b = h - km.xmax;
  var l = 0 - km.ymin;
  var r = w - km.ymax;
  var result = new T[h, w];

  for (int i = t; i < b; i++)
    for (int j = l; j < r; j++)
      result[i, j] = kernel(new RelativeArrayAccess2d<T>(source, i, j));
  return result;
}

На этом можно было бы и остановиться, если бы не ехидный комментарий в том же флейме на 30 уровней выше
Автор: alex_public
Дата: 28.03 17:37
:

Ну как же, если мы говорим про работу с коллекциями (а не про СУБД), то обычный императивный код (типа for() ... ) никуда не делся. И при этом он по прежнему является как наиболее обобщённым, так и наиболее производительным.

CHALLENGE ACCEPTED!
По поводу "обобщённости" спорить смысла нет — очевидно, что приведённый from d in data select... можно без изменений заставить работать не только с 2d-массивом, но и с произвольной структурой, которая представляет из себя двумерную сетку, индексированную целочисленными координатами. Например, мы можем работать с "изображением" из double с размерами миллион на миллион, то есть имея 10^12 точек или около 7 терабайт. Перебор их циклом for упрётся даже не во время, а в размер памяти и диска локальной машины. Так что мы пилим всё на тайлы, раскидываем по нескольким тысячам серверов, отправляем код для исполнения, дожидаемся, мёрджим. Ну, то есть получаем опять хэндл чудовищного "холста" для дальнейших операций — например, edge-detect, а потом ещё свёртка, а потом понижение размерности, и уже в конце получаем что-нибудь приемлемое для памяти локальной машины. И всё — в компактном и удобном для чтения (и написания) виде.

Вопрос у нас к тому пункту, где про "наиболее производительным".
Давайте напишем тот идеальный императивный код, который бы написал на C# программист, не доверяющий всем этми новомодным линкам, с их "динамикой и рефлексией":
static int[,] PerfectFourNeighborAverage(this int[,] source)
{
  var h = source.GetLength(0);
  var w = source.GetLength(1);
  var result = new int[h, w];

  for (int i = 1; i < h-1; i++)
    for (int j = 1; j < w-1; j++)
      result[i, j] = (source[i, j-1] + source[i, j+1] + source[i-1, j] + source[i+1,j])/4;
  return result;
}

Возьмём случайный массив размером 1000*1000, и попробуем сравнить производительность нашего "обобщённого" кода со "специфичным":
BenchmarkDotNet=v0.10.14, OS=Windows 10.0.15063.1088 (1703/CreatorsUpdate/Redstone2)
Intel Core i7-6600U CPU 2.60GHz (Skylake), 1 CPU, 4 logical and 2 physical cores
Frequency=2742190 Hz, Resolution=364.6720 ns, Timer=TSC
  [Host]     : .NET Framework 4.6.1 (CLR 4.0.30319.42000), 32bit LegacyJIT-v4.7.2650.0
  DefaultJob : .NET Framework 4.6.1 (CLR 4.0.30319.42000), 32bit LegacyJIT-v4.7.2650.0

     Method |      Mean |     Error |    StdDev | Scaled | ScaledSD |
----------- |----------:|----------:|----------:|-------:|---------:|
    Perfect | 15.008 ms | 0.1706 ms | 0.1512 ms |   1.00 |     0.00 |
       Linq | 33.423 ms | 0.6599 ms | 0.8345 ms |   2.23 |     0.06 |

Ну, на первый взгляд, linq справился не так уж и плохо — всего лишь в 2.23 раза медленнее.
Тем не менее, эта незначительная цена за удобную абстракцию легко может отпугнуть фанатов тактовыжимания.

Поэтому в следующей серии я попробую показать, что можно с этим сделать, ничего не меняя в коде бенчмарка. То есть те самые преимущества, которые бесплатны для прикладного программиста — от него потребуется разве что перекомпилировать свой проект. Никаких хинтов в коде — ни императивного развёртывания циклов, ни ручной векторизации, ни даже умелой расстановки "декларативных" прагм.
Только хардкор, только голый from d in data select (d[-1, 0] + d[1, 0] + d[0, 1] + d[0, -1]) / 4.