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

    Многие слышали о великом и ужасном быстром преобразовании Фурье (БПФ / 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!

    Similar posts

    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More
    Ads

    Comments 76

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

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

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

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

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

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

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

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

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

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


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

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


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

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


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


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

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

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

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

                                                  Only users with full accounts can post comments. Log in, please.