Предисловие: Я пишу на Python более 6 лет и могу назвать себя профессионалом в этом языке. Недавно я даже написал о нем книгу. Однако последние 8 месяцев я переключился на D и уже 4 месяца активно участвую в разработке этого языка по части расширения стандартной библиотеки Phobos. Так же я участвовал в код-ревью модуля std.ndslice о котором и пойдет речь.

std.ndslice так же как и Numpy предназначен для работы с многомерными массивами. Однако в отличие от Numpy ndslice имет крайне низкий оверхэд так как базируется на ranges (диапазонах), которые используются в штатной библиотеке повсеместно. Ranges позволяют избежать лишние процедуры копирования, а так же позволяют красиво организовать ленивые вычисления.

В этой статье мне хотелось бы рассказать о том какие преимущества std.ndslice дает по сравнению с Numpy.

Почему программисты Python должны посмотреть в сторону D?


Потому что D позволяет писать практически такой же простой и понятный код как Python, при этом этот код будет работать значительно быстрее чем код на Python.

Следующий пример создает диапазон чисел от 0 до 999 используя функцию iota (аналог в Python называется xrange) и возвращает массив размерностью 5x5x40.

import std.range : iota;
import std.experimental.ndslice;

void main() {
    auto slice = sliced(iota(1000), 5, 5, 40);
}

Хотя D статически типизируемый язык, и явное указание типа повышает читаемость кода, чтобы программистам на Python было проще мы воспользуемся автоматической типизацией с использованием ключевого слова auto.

Функций sliced возвращает многомерный срез. На входе она может принимать как простые массивы так и ranges. В итоге у нас получился куб размером 5x5x40 состоящий из чисел от 0 до 999.

Пару слов о том что такое Ranges. На русский язык их правильнее переводить как диапазоны. Ranges позволяют описать правила перебора данных любого типа данных, бу��ь то класс или структура. Для этого достаточно чтобы вы реализовали функцию: front, возвращающую первый элемент диапазона, popFront, переходящий к следующему элементу и empty возвращающую булевое значение показывающую, что перебираемая последовательность пуста. Ranges позволяют выполнять ленивый перебор обращаясь к данным по мере их необходимости. Более подробно концепция Ranges освещена тут.

Обратите внимание, что никаких пустых аллокаций памяти! Это происходит потому что iota так же позволяет генерировать ленивые диапазоны, а sliced в ленивом режиме так же принимает данные из iota и обрабатывает их по мере поступления.

Как вы видите std.ndslice работает несколько иначе чем Numpy. Numpy создает собственный тип для данных, тогда как std.ndslice представляет способ манипуляции с уже существующими данными. Это позволяет вам во всей вашей программе использовать одни и те же данные не тратя ресурсы на бесполезный мемори аллокейшен! Не сложно догадаться, что это крайне серьезно сказывается на производительности ваших решений.

Давайте посмотрим на следующий пример. В нем мы будем получать данные из stdin, фильтровать только уникальные строки, сортировать их, и затем возвращать в stdout.

import std.stdio;
import std.array;
import std.algorithm;

void main() {
    stdin
        // get stdin as a range
        .byLine(KeepTerminator.yes)
        .uniq
        // stdin is immutable, so we need a copy
        .map!(a => a.idup)
        .array
        .sort
        // stdout.lockingTextWriter() is an output range, meaning values can be
        // inserted into to it, which in this case will be sent to stdout
        .copy(stdout.lockingTextWriter());
}

Желающим более подробно разобраться с ленивой генерацией диапазонов рекомендую почитать эту статью.

Так как slice имеет три измерения это диапазон который возвращает диапазон диапазонов (ranges of ranges). Это хорошо видно на следующем примере:

import std.range : iota;
import std.stdio : writeln;
import std.experimental.ndslice;

void main() {
    auto slice = sliced(iota(1000), 5, 5, 40);

    foreach (item; slice) {
        writeln(item);
    }
}

Результат его работы будет следующим (троеточия для краткости):

[[0, 1, ... 38, 39], [40, 41, ... 78, 79], [80, 81, ... 118, 119], [120, 121, ... 158, 159], [160, 161, ... 198, 199]]

...

[[800, 801, ... 838, 839], [840, 841, ... 878, 879], [880, 881, ... 918, 919], [920, 921, ... 958, 959], [960, 961, ... 998, 999]]

Цикл foreach работает почти так же как for в Python. Однако в D вы можете его использовать так в стиле Cи, так и на манер циклов в Python, но без мороки с enumerate или xrange.

Используя UFCS (Uniform Function Call Syntax) вы можете переписать код на следующий манер:

import std.range : iota;
import std.experimental.ndslice;

void main() {
    auto slice = 1000.iota.sliced(5, 5, 40);
}

UFCS позволяет записывать вызов методов по цепочке и писать:

a.func(b)

вместо: 

func(a, b)

Давайте теперь сгенирируем пустой проект при помощи менеджера пакетов dub. Командой dub init и в файле \source\app.d напишем:

import std.experimental.ndslice;

void main() {
}

В настоящий момент std.experimental.ndslice; находится в секции std.experimental. Это не значит, что он сырой. Это значит, что ему нужно настояться.

Теперь соберем ��роект командой:

dub

Модуль D ndslice весьма похож на Numpy:

a = numpy.arange(1000).reshape((5, 5, 40))
b = a[2:-1, 1, 10:20]

Равнозначно:

auto a = 1000.iota.sliced(5, 5, 40);
auto b = a[2 .. $, 1, 10 .. 20];

Теперь давайте создадим двухмерный массив и получим середину каждой его колонки.

Python: 

import numpy

data = numpy.arange(100000).reshape((100, 1000))
means = numpy.mean(data, axis=0)

D:

import std.range;
import std.algorithm.iteration;
import std.experimental.ndslice;
import std.array : array;

void main() {
    auto means = 100_000.iota
        .sliced(100, 1000)
        .transposed
        .map!(r => sum(r) / r.length)
        .array;
}

Для того, чтобы данный код не работал в ленивом режиме мне пришлось вызывать метод array. Однако в реальном приложении результат не будет рассчитан, пока он используется в другой части программы.

В настоящий момент в Phobos нет встроенных модулей работы со статистикой. Поэтому в примере используется простая лямбда для нахождения среднего значения. Функция map! имеет в конце восклицательный знак. Это значит, что это шаблонная функция. Она позволяет на этапе компиляции генерировать код основанный на выражении указанном в ее теле. Вот хорошая статья о том как работают сами шаблоны в D.

Хотя код на D и получился чуть более многословным, чем на Python, но благодаря map! код будет работать с любыми входными данными являющимися диапазоном (range). В то время как код на Python будет работать только со специальными массивами из Numpy.

Тут нужно сказать, что после проведения этого теста оказалось, что Python проиграл D в разы и после долгих дисскусий на Hacker News, я понял что допустил ошибку и сравнение оказалось не совсем корректным. iota динамически создает данные которые принимает функция sliced. И соответственно мы не трогаем память до момента последней ее релокации. Так же D возвращает массив с типом данных long в то время как Numpy из double. В итоге переписал код и довел значение массива до 1000 000 вместо 10 000. Вот что получилось:

import std.range : iota;
import std.array : array;
import std.algorithm;
import std.datetime;
import std.conv : to;
import std.stdio;
import std.experimental.ndslice;

enum test_count = 10_000;

double[] means;
int[] data;

void f0() {
    means = data
        .sliced(100, 10000)
        .transposed
        .map!(r => sum(r, 0L) / cast(double) r.length)
        .array;
}

void main() {
    data = 1_000_000.iota.array;

    auto r = benchmark!(f0)(test_count);
    auto f0Result = to!Duration(r[0] / test_count);
    f0Result.writeln;
}

Тесты проводил на 2015 MacBook Pro with a 2.9 GHz Intel Core Broadwell i5. В Python для замера скорости использовал функцию %timeit в D std.datetime.benchmark. Компилировал все при помощи LDC v0.17 со следующими ключами: ldmd2 -release -inline -boundscheck=off -O. Или если вы используете dub то аналогом этих ключей будут опции dub --build=release-nobounds --compiler=ldmd2.

Вот результаты первого теста:

Python: 145 µs
LDC:      5 µs

D is 29x faster

Вот теста исправленной версии:

Python: 1.43 msec
LDC:  628    μs

D is 2.27x faster

Согласитесь весьма не плохая разница учитывая то что Numpy написан на С, а в D все ругают сборщик мусора.

Как D позволяет избежать проблем Numpy?


Да Numpy быстр, но быстр он лишь в сравнении с простыми масисвами Python. Но проблема в том, что он с этими простыми массивами совместим лишь частично.

Библиотека Numpy находится где то сбоку от остального Python. Она живет своей жизнью. В ней используются свои функции, она работает со своими типами. К примеру если нам потребуется использовать массив созданный в Python в NumPy нам нужно будет использовать np.asarray который скопирует его в новую переменную. Беглый поиск по github показал что тысячи проектов используют этот костыль. В то время как данные могли бы быть просто переданы из одной функции в другую без этих пустых копирований.

import numpy as np

a = [[0.2,0.5,0.3], [0.2,0.5,0.3]]
p = np.asarray(a)
y = np.asarray([0,1])

Эту проблему пытаются решить переписав часть стандартной библиотеки Python на использование типов Numpy. Однако это все равно не полноценное решение, которое приводит к замечательным приколам, когда написав:

sum(a)

вместо:

a.sum()

мы получаем 10x кратное падение в скорости.

У D таких проблем просто нет by design. Это компилируемый, статически типизируемый язык. Во время генерации кода известны все типы переменных. В самом std.ndslice вы получаете полный доступ ко всем функциям библиотеки Phobos к примеру к таким замечательным вещам как std.algorithm и std.range. Ах да, при этом вы сможете использовать шаблоны D позволяющие генерировать код прямо на этапе компиляции.

Вот пример:

import std.range : iota;
import std.algorithm.iteration : sum;
import std.experimental.ndslice;

void main() {
    auto slice = 1000.iota.sliced(5, 5, 40);
    auto result = slice
        // sum expects an input range of numerical values, so to get one
        // we call std.experimental.ndslice.byElement to get the unwound
        // range
        .byElement
        .sum;
}

Вы берете и просто используете функцию sum и она просто берет и работает, ровно как и любая другая функция из базовой библиотеки.

В самом Python для того чтобы получить список опреленной длинны инициализированный определенным значением нам нужно написать:

a = [0] * 1000

В Numpy уже будет совсем по-другому:

a = numpy.zeros((1000))

И если вы не вспользуетесь Numpy то ваш код будет в 4 раза медленне не считая лишней операции копирования съедающей память. В D нам на помощь приходят range, который позволяют сделать эту же операцию быстро и без пустых операций копирования:

auto a = repeat(0, 1000).array;

И если нужно мы можем тут же вызвать ndslice:

auto a = repeat(0, 1000).array.sliced(5, 5, 40);

Главное преимущество Numpy в настоящее время это его распространенность. Сейчас он используется реально очень широко от банковских систем до машинного обучения. По нему очень много книг, примеров и статей. Однако математические возможности в D очевидно будут уже очень скоро расширены. Так автор ndslice заявил, что сейчас ведет работы над BLAS (Basic Linear Algebra Subprograms) для Phobos, который так же будет полностью интегрирован с ndslice и со стандартной библиотекой.

Мощная математическая подсистема позволит очень эффективно решать ряд задач где необходима работы с большими данными. К примеру системы компьютерного видения. Прототип одной из таких систем уже разрабатывается и называется DCV.

Обработка изображений на D


Следующий пример покажет как работает медианный фильтр позволяющий ликвидировать шумы на изображении. Функция movingWindowByChannel может быть использована так же в других фильтрах где требуется использование скользящее окно. movingWindowByChannel позволяет перемещаться по изображению используя скользящее окно. Каждое такое окно передается в фильтр, который рассчитывает значения пикселей на основании выбранной зоны.

Эта функция не обрабатывает зоны с частичным перекрытием. Однако с ее помощью можно вычеслить значения в них тоже. Для этого нужно создать увеличенное изображение с краями отражающими границы оригинального изображения и затем обработать его.

/**
Params:
    filter = unary function. Dimension window 2D is the argument.
    image = image dimensions `(h, w, c)`,
        where с is the number of channels in the image
    nr = number of rows in the window
    nс = number of columns in the window

Returns:
    image dimensions `(h - nr + 1, w - nc + 1, c)`,
        where с is the number of channels in the image.
        Dense data layout is guaranteed.
*/
Slice!(3, C*) movingWindowByChannel(alias filter, C)
(Slice!(3, C*) image, size_t nr, size_t nc)
{
    // local imports in D work much like Python's local imports,
    // meaning if your code never runs this function, these will
    // never be imported because this function wasn't compiled
    import std.algorithm.iteration: map;
    import std.array: array;

    // 0. 3D
    // The last dimension represents the color channel.
    auto wnds = image
        // 1. 2D composed of 1D
        // Packs the last dimension.
        .pack!1
        // 2. 2D composed of 2D composed of 1D
        // Splits image into overlapping windows.
        .windows(nr, nc)
        // 3. 5D
        // Unpacks the windows.
        .unpack
        // 4. 5D
        // Brings the color channel dimension to the third position.
        .transposed!(0, 1, 4)
        // 5. 3D Composed of 2D
        // Packs the last two dimensions.
        .pack!2;

    return wnds
        // 6. Range composed of 2D
        // Gathers all windows in the range.
        .byElement
        // 7. Range composed of pixels
        // 2D to pixel lazy conversion.
        .map!filter
        // 8. `C[]`
        // The only memory allocation in this function.
        .array
        // 9. 3D
        // Returns slice with corresponding shape.
        .sliced(wnds.shape);
}

Следующая функция пример того как можно рассчитать медиану у объекта. Функция сильно упрощена в целях повышения читаемости.

/**
Params:
    r = input range
    buf = buffer with length no less than the number of elements in `r`
Returns:
    median value over the range `r`
*/
T median(Range, T)(Range r, T[] buf)
{
    import std.algorithm.sorting: sort;

    size_t n;

    foreach (e; r) {
        buf[n++] = e;
    }

    buf[0 .. n].sort();
    immutable m = n >> 1;
    return n & 1 ? buf[m] : cast(T)((buf[m - 1] + buf[m]) / 2);
}

Ну а теперь собственно сам Main:

void main(string[] args)
{
    import std.conv: to;
    import std.getopt: getopt, defaultGetoptPrinter;
    import std.path: stripExtension;

    // In D, getopt is part of the standard library
    uint nr, nc, def = 3;
    auto helpInformation = args.getopt(
        "nr", "number of rows in window, default value is " ~ def.to!string, &nr,
        "nc", "number of columns in window, default value is equal to nr", &nc);

    if (helpInformation.helpWanted)
    {
        defaultGetoptPrinter(
            "Usage: median-filter [<options...>] [<file_names...>]\noptions:",
            helpInformation.options);
        return;
    }

    if (!nr) nr = def;
    if (!nc) nc = nr;

    auto buf = new ubyte[nr * nc];

    foreach (name; args[1 .. $])
    {
        import imageformats; // can be found at code.dlang.org

        IFImage image = read_image(name);

        auto ret = image.pixels
            .sliced(cast(size_t)image.h, cast(size_t)image.w, cast(size_t)image.c)
            .movingWindowByChannel
                !(window => median(window.byElement, buf))
                 (nr, nc);

        write_image(
            name.stripExtension ~ "_filtered.png",
            ret.length!1,
            ret.length!0,
            (&ret[0, 0, 0])[0 .. ret.elementsCount]);
    }
}

Если не все примеры показались вам понятными рекомендую вам прочитать бесплатную версию замечательной книги Programming in D.

P.S. Если знаете как данную публикацию перевести в статус "пе��еводы", то напишите в приват.