Игра Жизнь и преобразование Фурье

    Многие слышали о великом и ужасном быстром преобразовании Фурье (БПФ / FFT — fast fourier transform) — но как его можно применять для решения практических задач за исключением JPEG/MPEG сжатия и разложения звука по частотам (эквалайзеры и проч.) — зачастую остается неясным вопросом.

    Недавно я наткнулся на интересную реализацию игры «Жизнь» Конвея, использующую быстрое преобразование Фурье — и надеюсь, оно поможет вам понять применимость этого алгоритма в весьма неожиданных местах.

    Правила

    Вспомним правила классической «Жизни» — на поле с квадратными клетками, живая клетка погибает если у неё больше 3 или меньше 2 соседей, и если у пустой клетки ровно 3 соседей — она рождается. Соответственно, для эффективной реализации алгоритма нужно быстро считать количество соседей вокруг клетки.

    Алгоритмов для этого существует целая куча (в том числе и моя JS реализация). Но есть у задачи и математическое решение, которое может давать хорошую скорость для плотно заполненных полей, и быстро уходит в отрыв с ростом сложности правил и площади/объема суммирования — например в Smoothlife-подобных (1,2), и 3D вариантах.

    Реализация на БПФ

    Идея алгоритма следующая:
    1. Формируем матрицу суммирования (filter), где 1 стоят в ячейках, сумму которых нам нужно получить (8 единиц, остальные нули). Выполняем над матрицей прямое преобразование Фурье (это нужно сделать только 1 раз).
    2. Выполняем прямое преобразование Фурье над матрицей с содержимым игрового поля.
    3. Перемножаем каждый элемент результата на соответствующий элемент матрицы «суммирования» из пункта 1.
    4. Выполняем обратное преобразование Фурье — и получаем матрицу с нужной нам суммой количества соседей для каждой клетки.
    Весь этот процесс называется сверткой / сonvolution.

    Реализация на C++
    //Author: Mark VandeWettering http://brainwagon.org/2012/10/11/crazy-programming-experiment-of-the-evening/
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <complex.h>
    #include <unistd.h>
    #include <fftw3.h>
     
    #define SIZE    (512)
    #define SHIFT   (18)
     
    fftw_complex *filter ;
    fftw_complex *state ;
    fftw_complex *tmp ;
    fftw_complex *sum ;
     
    int
    main(int argc, char *argv[])
    {
        fftw_plan fwd, rev, flt ;
        fftw_complex *ip, *jp ;
        int x, y, g ;
     
        srand48(getpid()) ;
     
        filter = (fftw_complex *) fftw_malloc(SIZE * SIZE * sizeof(fftw_complex)) ;
        state = (fftw_complex *) fftw_malloc(SIZE * SIZE * sizeof(fftw_complex)) ;
        tmp = (fftw_complex *) fftw_malloc(SIZE * SIZE * sizeof(fftw_complex)) ;
        sum = (fftw_complex *) fftw_malloc(SIZE * SIZE * sizeof(fftw_complex)) ;
     
        flt = fftw_plan_dft_2d(SIZE, SIZE, 
                    filter, filter, FFTW_FORWARD, FFTW_ESTIMATE) ;
        fwd = fftw_plan_dft_2d(SIZE, SIZE, 
                    state, tmp, FFTW_FORWARD, FFTW_ESTIMATE) ;
        rev = fftw_plan_dft_2d(SIZE, SIZE, 
                    tmp, sum, FFTW_BACKWARD, FFTW_ESTIMATE) ;
     
        /* initialize the state */
        for (y=0, ip=state; y<SIZE; y++) {
            for (x=0; x<SIZE; x++, ip++) {
                *ip = (fftw_complex) (lrand48() % 2) ;
            }
        }
     
        /* initialize and compute the filter */
     
        for (y=0, ip=filter; y<SIZE; y++, ip++) {
            for (x=0; x<SIZE; x++) {
                *ip = 0. ;
            }
        }
     
     
    #define IDX(x, y) (((x + SIZE) % SIZE) + ((y+SIZE) % SIZE) * SIZE)
        filter[IDX(-1, -1)] = 1. ;
        filter[IDX( 0, -1)] = 1. ;
        filter[IDX( 1, -1)] = 1. ;
        filter[IDX(-1,  0)] = 1. ;
        filter[IDX( 1,  0)] = 1. ;
        filter[IDX(-1,  1)] = 1. ;
        filter[IDX( 0,  1)] = 1. ;
        filter[IDX( 1,  1)] = 1. ;
     
        fftw_execute(flt) ;
         
        for (g = 0; g < 1000; g++) {
            fprintf(stderr, "generation %03d\r", g) ;
             
            fftw_execute(fwd) ;
     
            /* convolve */
            for (y=0, ip=tmp, jp=filter; y<SIZE; y++) {
                for (x=0; x<SIZE; x++, ip++, jp++) {
                    *ip *= *jp ;
                }
            }
     
            /* go back to the sums */
            fftw_execute(rev) ;
     
            for (y=0, ip=state, jp=sum; y<SIZE; y++) {
                for (x=0; x<SIZE; x++, ip++, jp++) {
                    int s = (int) round(creal(*ip)) ;
                    int t = ((int) round(creal(*jp))) >> SHIFT ;
                    if (s) 
                        *ip = (t == 2 || t == 3) ;
                    else
                        *ip = (t == 3) ;
                }
            }
     
            /* that's it!  dump the frame! */
     
            char fname[80] ;
            sprintf(fname, "frame.%04d.pgm", g) ;
            FILE *fp = fopen(fname, "wb") ;
            fprintf(fp, "P5\n%d %d\n%d\n", SIZE, SIZE, 255) ;
     
            for (y=0, ip=state; y<SIZE; y++) {
                for (x=0; x<SIZE; x++, ip++) {
                    int s = ((int) creal(*ip)) ;
                    fputc(255*s, fp) ;
                }
            }
     
            fclose(fp) ;
        }
        fprintf(stderr, "\n") ;
     
        return 0 ;
    }

    Для сборки нужна библиотека FFTW. Ключи для сборки в gcc:
    gcc life.cpp -lfftw3 -lm -lstdc++
    в Visual Studio нужны изменения в работе с комплексными числами.

    Изображения во время выполнения алгоритма:
    Исходная матрица фильтра (начало координат — левый верхний угол, поле «закольцовано»):

    Действительная часть фильтра после прямого БПФ:

    Исходное поле — 1 глайдер:

    БПФ исходного поля, действительная и мнимая части:

    После умножения на матрицу фильтра:

    После обратного БПФ — получаем суммы:


    Результат вполне ожидаемый:


    «Закольцовывание» поля получается автоматически из-за БПФ.

    Плюшки БПФ

    • Вы можете суммировать любое количество элементов с любыми коэффициентами — а время работы остается фиксированным, N2logN. Т.е. если для классической жизни — обычные алгоритмы на заполненных полях все еще достаточно быстрые, то с увеличением площади/объема суммирования — они становятся все медленнее, а скорость работы БПФ остается фиксированной.
    • БПФ — уже написан, отлажен и оптимизирован идеально — с использованием всех возможностей процессора: sse, sse2, avx, altivec и neon.
    • Вы легко можете использовать все процессоры и видеокарты взяв готовые многопроцессорные и CUDA/OpenCL реализации БПФ. Опять же, об оптимизации БПФ вам заботится не нужно. FFTW кстати поддерживает многопроцессорное выполнение (--enable-openmp во время сборки), чем больше поле — тем меньше потери на синхронизацию.
    • Тот же подход применим и к 3D пространству.

    Сравниваем скорость работы с «наивной» реализацией

    FFTW собран в режиме одинарной точности (--enable-float, дает прирост скорости примерно в 1.5 раза). Процессор — Core2Duo E8600. gcc 4.7.2 -O3
    Исходный текст наивной реализации
    for (int i=8; i<1000; i++)
        {
            delta[i][0] = (lrand48() % 21)-10;
            delta[i][1] = (lrand48() % 21)-10;
        }
    
        #define NUMDELTA 18
        for(int frame=0;frame<10;frame++)
        {
            double start = ticks();
            for (y=0; y<SIZE; y++)
                for (x=0; x<SIZE; x++) 
                {
                    int sum=0;
                    for(int i=0;i<NUMDELTA;i++)
                        sum+=array[readpage][(SIZE+y+delta[i][0]) % SIZE][(SIZE+x+delta[i][1]) % SIZE];
                    if(array[readpage][y][x])
                        array[1-readpage][y][x] = (sum>3 || sum<2)?0:1;else
                        array[1-readpage][y][x] = (sum==3)?1:0;
                }
            readpage = 1-readpage;
            printf("NaiveDelta: %10.10f\n", ticks()-start);
        }
    
    Конфигурация теста При каком количестве суммируемых клеток скорость FFT и наивной реализации равна
    1 ядро, 512x512 29
    2 ядра, 512x512 18
    1 ядро, 4096x4096 88
    2 ядра, 4096x4096 65
    Как видим, хоть асимптотика и против FFT — но если суммировать нужно сотни и тысячи клеток — то FFT выиграет.

    Update:
    Как оказалось, FFTW при указании флага FFTW_ESTIMATE использует далеко не оптимальные планы вычисления FFT. При указании FFTW_MEASURE — скорость сильно выросла, и ситуация выглядит уже радостнее (при суммировании менее 18 клеток — наивная реализация резко становится в 3 раза быстрее, потому тут все упирается в 18):
    Конфигурация теста При каком количестве суммируемых клеток скорость FFT и наивной реализации равна
    1 ядро, 512x512 18
    2 ядра, 512x512 18
    1 ядро, 4096x4096 28
    2 ядра, 4096x4096 19


    Чтобы ни говорили — математика в программировании иногда может пригодиться в самых неожиданных местах.
    И да пребудет с вами сила FFT!
    Поделиться публикацией

    Комментарии 76

      +21
      Было бы очень ценно показать «на пальцах».
      Взять маленькую игровую матрицу (4x4 достаточно), показать её ПФ, показать результат умножения и результат обратного преобразования.

      Каждый бы смог повторить в матлабе.
        0
        Хорошая идея, матрицы сейчас сделаю. Вот только скрипты для матлаба предлагаю написать специалистам по матлабу.
        Я же из C++ получу промежуточные изображения.
          +3
          Матлаб фтопку, даёшь православный сионизмъ!
          (кроме шуток, примеры на матлабе уже задрали)
            +4
            Скрипт не важен. Нужны тестовые данные. Чтобы человек, далёкий от c++, хотящий повторить это в своём любимом maple/python/подставить-нужное, мог на каждом шаге видеть, правильно он двигается или сделал ошибку.
              +1
              Скрипт на матлабе не содержал бы ничего, кроме идеи. Вам, видимо, нравится продираться сквозь malloc'и и макросы
            +2
            Промежуточные изображения добавил в пост.
              +2
              Спасибо.

              У меня было предположение, что можно двигать эволюцию простым умножением без обратного преобразования на каждом шаге, и только после N шагов снять результат.
              С примером я понял, что после каждого шага нужно применять пороговые фильтры, которые обнулят клетки с большими суммами и с маленькими, что невозможно сделать, оставаясь в спектральном представлении.
              • НЛО прилетело и опубликовало эту надпись здесь
            +4
            powdertoy.co.uk/ — в этой «игре» симуляция «жизни» аналогичным образом сделана.
              +2
              Насколько я знаю, в The Powder Toy используется обычный подсчёт соседей. Может быть Golly использует БПФ как один из вариантов, хотя там вроде Hashlife главный алгоритм.
            • НЛО прилетело и опубликовало эту надпись здесь
                +4
                все и так знают.

                В этом месте по-подробнее. Все это кто? Все 292 тыс. пользователей хабра? :-)

                Те, кто хорошо учился в университете на мат./физ. факультете — да, крайне вероятно знают.
                Я например не знал, пока не разобрался с этим алгоритмом.
                • НЛО прилетело и опубликовало эту надпись здесь
                    +3
                    У вас логическая ошибка
                    Те кто не прогуливал && вообще имел в расписании курс математического анализа.
                    • НЛО прилетело и опубликовало эту надпись здесь
                        0
                        Не вполне согласен. Если «почувствовать» как работает БПФ, его можно применять, не вникая до самого дна внутренностей.

                        В этом и состоит разница между программированием как ремеслом, и программированием как искусством :-D
                        Зачастую нужно именно ремесло.
                        • НЛО прилетело и опубликовало эту надпись здесь
                            +4
                            Поддерживаю вашего оппонента — после прочтения вашей заметки стало таки интересно, но вовсе не стало понятно. Надо ширше наверное. Спасибо
                              0
                              Думаю, подробный рассказ о том, как работает и реализовывается БПФ и какие у него свойства — может стать хорошей темой для отдельной статьи, но не в моем исполнении. Я не чувствую, что смог бы рассказать об этом понятно.
                                +2
                                Попробуйте осилить кто-нибудь, пожалуйста. Мне лично очень интересно, но данных по теме не хватает однозначно.

                                P.S. гуманитарий же )
                            • НЛО прилетело и опубликовало эту надпись здесь
                                0
                                Что такое «полный смысл»? ДПФ двухточечной последовательности состоит из суммы и разности ее элементов (с точностью до нормирующего множителя); ДПФ одноточечной последовательности совпадает с ней самой. Никто смысла не теряет.
                                • НЛО прилетело и опубликовало эту надпись здесь
                                    0
                                    Вы знаете, терминологию не просто так люди придумали, а чтобы было более-менее понятно то, что говорят другие. Я вот этой фигней занимаюсь, сколько себя помню, но, уж извините, у Вас не понял ни одного слова. Что такое «спектральное поле»? Что такое «0 точка»? Дальше тоже непонятно. Не бывает никакой действующей частоты. Бывает действующее значение.

                                    По существу: при любом количестве точек, начиная с одной, ДПФ сохраняет свой смысл, включая физическую интерпретацию амплитудного и фазового спектров.
                                    • НЛО прилетело и опубликовало эту надпись здесь
                                        0
                                        Мне кажется, что физикам и электронщикам как раз легче понимать смысл спектральных методов, чем математикам и программистам. Это же просто набор полосовых фильтров, не более того. Правда фильтры так хитро подобраны, что их частотные характеристики в сумме дают тождественную единицу, поэтому по их откликам можно восстановить исходный сигнал.

                                        А понимать отдельно смысл ДПФ/БПФ по-моему и не нужно, потому что это просто схемы вычисления настоящего (интегрального) спектра. Причем БПФ вообще является на самом деле очень хитрой эксплуатацией идеи факторизации матрицы, реализующей ДПФ.

                                        Про всю эту Фурьевскую магию написано очень много и есть книжки, где все очень хорошо разжевано (та же книжка Сато). Мне кажется, что вряд ли стоит все это повторно разжевывать на Хабре. А вот про приложения я собираюсь пару статей написать — про тот же Шазам и про передискретизацию картинок. А томография — очень специфичная тема и там само БПФ не так уж и важно на самом деле, ведь первые томографы вообще не обрабатывали получаемые картинки, а использовали аппаратно реализованное т.н. обратное проецирование и спектральные методы просто позволили улучшить качество картинки.
                                        • НЛО прилетело и опубликовало эту надпись здесь
                                          • НЛО прилетело и опубликовало эту надпись здесь
                              0
                              Это совсем не первый курс математического анализа.
                              Первый курс — это пределы и производные.
                              Преобразование Фурье — в лучшем случае четвёртый семестр.
                                0
                                Третий семестр это, второй курс соответственно. По крайней мере, у меня было так.
                        +1
                        может давать хорошую скорость для плотно заполненных полей, и быстро уходит в отрыв с ростом сложности правил и площади/объема суммирования


                        Но для классических правил оно все же будет заметно медленнее работать, чем прямое вычисление. В основном за счет большой константы перед N2 log N и за счет худшей локальности при обращении к памяти.
                          +1
                          В среднем да, если конечно не брать в расчет варианты вроде «код на Python использует библиотеку на GPU/многопроцессорной системе».
                          • НЛО прилетело и опубликовало эту надпись здесь
                              0
                              Нет, имено N2 log2 N. Посмотрите где-нибудь в книжке. Подсказка: в изображении размером N x N пикселов N2, так что за N log2 N их все даже просуммировать не получится.
                                0
                                Не рассмотерл у Вас слов «для одномерного». Но для одномерного оценка тоже другая — N log2 N.
                                • НЛО прилетело и опубликовало эту надпись здесь
                                    +1
                                    Совсем я запутался, извините. Мне почему-то сначала показалась, что у Вас цифирька 2 вверху стоит, тогда как должна внизу. А оказывается, она у Вас посередине :) Сорри.
                                    • НЛО прилетело и опубликовало эту надпись здесь
                              0
                              Вот тут я понимаю Жизнь, даже музыка соответствующая www.youtube.com/watch?feature=player_detailpage&v=C2vgICfQawE
                                0
                                Хотелось бы ещё сравнения скорости с наивным алгоритмом, хотя бы простейшего — чтобы вообще представлять насколько это реально быстрее/медленнее на разных количествах суммируемых элементов. А вообще, конечно, БПФ отличная вещь, столько применений что многие даже в голову не приходят, как например этот, или там поиск подстрок какой-нибудь.
                                  0
                                  Да, думаю было бы полезно, вечером проведу замеры.
                                  • НЛО прилетело и опубликовало эту надпись здесь
                                      0
                                      Это всё конечно хорошо, только речь о другом. В статье сравнивается на словах алгоритм с БПФ и наивный алгоритм для игры «Жизнь», без ДПФ. Вот о их скорости и вопрос.

                                      Кстати, вы как-то странно противопоставляете ДПФ и БПФ, хотя БПФ — один из алгоритмов для вычисления ДПФ.
                                      • НЛО прилетело и опубликовало эту надпись здесь
                                          +1
                                          У БПФ в данном случае асимптотика O(N^2 log N), а у наивного алгоритма — O(N^2), поэтому асимптотически соотношение как раз не в пользу первого, но тут N не такие большие и было бы интересно сравнить реальное время. Вы, видимо, не совсем поняли что с чем сравнивается. Кстати, основание логарифма тут вообще не важно, ни в статье ни в комментариях здесь мы всё равно мультипликативные константы не пишем (а логарифмы по разным основаниям как раз на такие константы отличаются).
                                          • НЛО прилетело и опубликовало эту надпись здесь
                                              0
                                              «здесь мы всё равно мультипликативные константы не пишем (а логарифмы по разным основаниям как раз на такие константы отличаются»
                                              — вот тут я не понял…


                                              Мультипликативная константа — такая, на которую умножается (ну или делится :) ) Так как loga n = log2n / log2a, то в асимптотиках и не пишут обычно основание логарифма, также как и константы (я имею в виду, что не пишут например что количество операций в алгоритме равно 7*n, а пишут что O(n)).

                                              «поэтому асимптотически соотношение как раз не в пользу первого, но тут N не такие большие и было бы интересно сравнить реальное время. Вы, видимо, не совсем поняли что с чем сравнивается.»
                                              вы хотели что сравнить?
                                              1.однобитный целочисленный алгоритм суммирования окружения клетки и БПФ на двойной точности?
                                              2.Асимптота логарифма затормаживается, а N в квадрате — растёт как взрыв. Вы наверное перепутали. С ростом N — БПФ выигрывает всё сильнее…


                                              Запишите асимптотики для БПФ в данном применении и для этого однобитового (хотя это не важно конечно) алгоритма, и увидите (n^2 log n против n^2).
                                              • НЛО прилетело и опубликовало эту надпись здесь
                                      0
                                      Замеры скорости проведены, результаты в конце статьи.
                                        0
                                        Как по мне, так хорошие результаты. Думал, БПФ хуже себя покажет.
                                      +6
                                      Многие слышали о великом и ужасном быстром преобразовании Фурье (БПФ / FFT — fast fourier transform) — но как его можно применять для решения практических задач за исключением JPEG/MPEG сжатия — зачастую остается неясным вопросом.


                                      ну вы сказали! классический алгоритм перевода из time domain в frequency domain, который используется ВЕЗДЕ! Даже если нам нужно разобрать тональный сигнал телефона (а там сумма двух волн)… или там нарисовать частоты эквалайзера… если честно, у меня даже мысль останавливается, потому что он реально везде, а «Жизнь» — это всего лишь очень экзотическое применение.
                                        +1
                                        Согласен, добавил это :-)
                                        +1
                                        А как обстоят дела с точностью БПФ и ошибками округления — не было проблем с этим?
                                          0
                                          Суммы округляются до ближайшего целого при проверках. Из-за того, что результаты БПФ напрямую в следующих итерациях не используются — ошибки не накапливаются.

                                          Проблем не было, что не удивительно — приведенный пример вообще нагло считает все в double. Но думаю, было бы все хорошо и с single, и с фиксированной точкой.
                                            +1
                                            Единственный минус БПФ — потребление памяти при таком подходе очень большое. Примерно в 64-128 раз больше по сравнении с битовым (или двухбитовым) представлением поля. Это может стать существенным минусом при больших размерах.
                                            • НЛО прилетело и опубликовало эту надпись здесь
                                                0
                                                Вы что именно имеете в виду? В данном случае не столь важно, как в процессоре это считается, а сколько занимает массив из этих вещественных чисел. А он как раз занимает в 32 или 64 раза больше, чем одно- или двухбитовый.
                                                • НЛО прилетело и опубликовало эту надпись здесь
                                                    +1
                                                    Так сопроцессор в классическом виде — медленно но верно уходит в прошлое.
                                                    На x87 инструкциях быстрый FFT не написать.

                                                    Современный математический код использует SSEx инструкции — и разрядность вещественных чисел не больше 64.
                                                    80 бит вещественные числа — остались только для совместимости, и очень медленны по сравнению с SSEx.
                                                    • НЛО прилетело и опубликовало эту надпись здесь
                                            • НЛО прилетело и опубликовало эту надпись здесь
                                          • НЛО прилетело и опубликовало эту надпись здесь
                                              +2
                                              ИМХО, библиотека FFTW заоптимизирована вдоль и поперек — одна из лучших по производительности. Ничего существенно быстрее написать не получится
                                                0
                                                Только если заранее для всех размеров входных массивов создать и сохранить профили. Если же считать без профиля, то она почти не обгоняет стандартную реализацию FFT (которая без использования преобразования Уолша).
                                                • НЛО прилетело и опубликовало эту надпись здесь
                                                    +1
                                                    FFTW оптимизирует схему вычислений в алгоритме FFT с учетом конфигурации кэшей конкретного процессора. Если простыми словами, то он решает, до какого размера еще стоит разбивать последовательность на две (из четных и нечетных элементов) и потом сливать их спектры, а начиная с какого размера можно использовать прямое вычисление (на самом деле там все чуть сложнее). Эта оптимизация проводится для каждой длины исходной последовательности путем выполнения неких бенчмарков на конкретной машине и результаты сохраняются в файл профиля. Потом этот профиль можно использовать при вычислении БПФ.
                                                    • НЛО прилетело и опубликовало эту надпись здесь
                                                        0
                                                        Это скорее всего размер, при котором данные перестают целиком помещаться в кэш процессора.
                                                        • НЛО прилетело и опубликовало эту надпись здесь
                                              +1
                                              Сравниваем скорость работы с «наивной» реализацией
                                              А если сравнить с HashLife?
                                                0
                                                HashLife — он бесспорно адски быстрый.
                                              • НЛО прилетело и опубликовало эту надпись здесь
                                                  +2
                                                  Спасибо за статью. Как-то раньше не приходил в голову очевидный факт, что клеточные автоматы — это просто свёртка.
                                                    +3
                                                    Это комбинация свёртки и дискретной (пороговой) функции.
                                                    Второй шаг не позволяет применить классический (непрерывный) мат. аппарат для анализа или прогнозирования.
                                                    Стивен Вольфрам доказал неприводимость клеточного автомата к более простой модели.
                                                      0
                                                      Прошу прощения за некропост
                                                      дискретной (пороговой) функции… не позволяет применить классический (непрерывный) мат. аппарат для анализа или прогнозирования
                                                      А почему нельзя применить какой-нибудь сигмоид?
                                                        +1
                                                        Сигмоид даст лишь приближенное решение — которое обязательно разойдется с исходным через некоторое число итераций.

                                                        Кроме того, в непрерывном случае также существуют такие стабильные «решения» как поле, заполненное одинаковыми клетками со значением примерно 1/4. Ценность этих решений применительно к исходной игре «Жизнь» — сомнительна.

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

                                                  Самое читаемое