- PVSM.RU - https://www.pvsm.ru -

Julia и клеточные автоматы

Julia и клеточные автоматы - 1

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

Для практики и лучшего понимания следует попытаться самому реализовать приведенные ниже алгоритмы. Ну а если лень, или просто из интереса, можно поиграться с различными имплементациями:

  • Golly [1] — небольшая программка использующая оптимальнейшие методы для моделирования клеточных автоматов с различными правилами и сетками. В комплекте есть база шаблонов и предустановок.
  • Lenia [2] — набор инструментов для изучения непрерывного варианта "жизни".
  • Онлайн симулятор [3] для муравья Лэнгтона. Особо примечательна возможность выбора триангулярной сетки, в которой стандартные правила ведут себя как генератор пиксельных космических кораблей

Julia и клеточные автоматы - 2

Одномерный случай

Одномерный клеточный автомат представляет с собой отрезок разбитый на равные части — ячейки. Каждая ячейка может принимать одно из двух состояний, условно 0 или 1, черный или белый, живой или мертвый и т. д. Состояние каждой ячейки меняется по определенному правилу, учитывающему состояние соседей. Если принимать во внимание только ближайших соседей, то мы получим комбинацию из восьми состояний искомой ячейки и пары соседей. Каждой из этих комбинаций можно поставить в соответствие одно новое состояние для искомой ячейки, из чего следует общее количество правил равное 256 (количество 8-битных слов).

Julia и клеточные автоматы - 3

Чтобы надежно закрепить образ в своей нейросети, реализуем функцию визуализирующую одномерные клеточные автоматы.

Для начала следует установить и подгрузить все необходимые пакеты:

Код

using Pkg
pkgs = ["Images", "ColorSchemes", "FFTW"]

for p in pkgs
    Pkg.add(p)
end

using Images, ColorSchemes, FFTW, LinearAlgebra: kron
using Random: bitrand
cd("C:\Users\User\Desktop\Mycop")

Правила для одномерных клеточных автоматов выражают числами от 0 до 255 (десятичные представления восьмибитных слов). Значит наша функция должна уметь представлять номер правила в двоичном виде, и составлять для него паттерн. Затем, при наивной реализации, происходит проход по всем ячейкам с заменой состояний:

Код

# размеры, правило, множитель для масштаба картинки
function cellauto( n::Int64, m::Int64, rule::Int64, s::Int64 = 1 )

    ptrn = digits(Bool, rule, base = 2, pad = 8)
    bt = [ bitstring(i)[end-2:end] for i = 0:7 ]
    d = Dict( bt[i] => ptrn[i] for i = 1:8 ) # словарь

    M = falses(n,m)
    M[1,m ÷ 2] = true
    #M[1,:] .= bitrand(m)
    for i = 1:n-1, j = 2:m-1
        key = "$(M[i, j-1]*1)$(M[i, j]*1)$(M[i, j+1]*1)"
        M[i+1,j] = d[key]
    end
    kron(M, ones(Int,s,s) ) # произведение Кронекера
end

M0 = cellauto(100, 200, 30, 4)
Gray.( M0 )

Julia и клеточные автоматы - 4

Это правило 30, одно из любимых Стивена Вольфрама. Он, кстати, назначил денежные призы за решения задач [4] связанных с этим правилом. Действительно, в этом что-то есть: периодичность граничащая с хаосом, плюс спонтанное зарождение упорядоченных структур… Но лично меня зацепили треугольники Серпинского, которые часто встречаются здесь и там [5]. А для одномерных КА они довольно частое явление:

Arr = cellauto.(40, 40, [0:23;], 2);
Imgs = [ Gray.(a) for a in Arr ]
reshape(Imgs, 4,6)

Julia и клеточные автоматы - 5

В кругах и в цвете

Приверженцам объектно ориентированного программирования можно предложить вариант [6], где клеточный автомат выполнен в виде структуры. В этом случае будет удобно реализовать варьирование цвета и даже круглые поля:

Julia и клеточные автоматы - 6

Julia и клеточные автоматы - 7

Визуализация выполнена с помощью пакета Luxor, который мы разбирали ранее [7].

Игра в жизнь

"Жизнь" Конвея — самый популярный из двумерных (да и трехмерных) клеточных автоматов. Эта модель отдаленно напоминает копошение бактерий в чашке Петри, и обладая Тьюринг полнотой, может послужить основой для чего угодно: моделирования мышления [8], компьютеров (архитектуры тетриса [9]) и самой себя [10].

Julia и клеточные автоматы - 8
играбельно [11]

Множество энтузиастов находят все новые структуры [12]: осцилляторы, глайдеры, корабли (на гусеничном ходу! [13]) и фабрики. То есть, вся прелесть здесь заключается в том, что на основе простых правил можно создавать сложные структуры, откуда возникает соблазн попытаться распространить данную модель на фундаментальные законы Естества. В частности С. Вольфрам балует нас подобными высказываниями [14]. Вопрос, на сколько это правдоподобно и, главное, полезно, всегда вызывает споры, так что, не будем на нем заостряться. Лучше реализуем "Жизнь" использующую преобразования Фурье.

Жизнь без Фурье и не жизнь вовсе

Тем, кто занимался обработкой изображений [15], работа с БПФ [16] в таком аспекте должна быть знакома. Остальным можно посоветовать попробовать заняться обработкой изображений, так как это очень интересно, и снабдит вас полезными абстракциями для понимания линейной алгебры.

За основу использовался алгоритм из хабровской статьи [17]

Код

function makefilter(N::Int64)
    filter = zeros(Complex, N, N);
    IDX(x, y) = ( (x + N) % N ) + ( (y+N) % N ) * N + 1
        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. ;

    return fft( real.(filter) )
end

function fftlife(N = 16, steps = 100, dx = 0, glider = true)

    if glider
        state = falses(N, N)
        state[4,5] = state[5,6] = state[6,6] = true
        state[6,5] = state[6,4] = true
    else
        state = bitrand(N, N)
    end

    filter = makefilter(N)

    for i in 1:steps
        tmp = fft( real.(state) ) # forward
        tmp .*= filter

        summ = ifft(tmp) # reverse

        for i in eachindex(state)
            t = round( Int, real(summ[i]) ) >> dx
            state[i] = ( state[i] ? t == 2 || t == 3 : t == 3 )
        end
        save("KonLife_$(N)x$(N)_$i.png", kron( Gray.(state), ones(8,8) ) )
    end

end

fftlife(16, 60)

Julia и клеточные автоматы - 9

Приятные бонусы этого подхода — это сравнительно быстрая отработка БПФ в силу его отполированности, и закольцованность граничных условий (действо происходит как бы на поверхности тора). Собственно, теперь можно обобщить "жизнь" непрерывным случаем.

Гладкая жизнь

Вот так в общем виде выглядит модель:

$ frac{df(x,t)}{dt}=Sleft(N(x,t),M(x,t) right ) - f(x,t) \ \ S(n,m)=sigma_nleft(sigma_m(b_1, b_2), sigma_m(d_1,d_2)right)\ M(x,t)=frac{1}{pi h^2} intlimits_{0le |y| le h} f(x-y, t) dy \ N(x,t)=frac{1}{8pi h^2} intlimits_{0le |y| le 3h} f(x-y, t) dy $

Состояние теперь описывается не 1 и 0, а сигмоидой, соседи учитываются в некотором радиусе, а еще понадобится решать дифуры.

Код

Начальные условия

function clamp(x)
    y = copy(x)
    y[x.>1] .= 1
    y[x.<0] .= 0
    y
end

function func_linear(X, a, b)
    Y = [ (x-a + 0.5b)/b for x in X ]

    Y[X.<a-0.5b] .= 0
    Y[X.>a+0.5b] .= 1
    return Y
end

function splat!(aa, ny, nx, ra)
    x = round(Int, rand()*nx ) + 1
    y = round(Int, rand()*ny ) + 1
    c = rand() > 0.5

    for dx = -ra:ra, dy = -ra:ra
        ix = x+dx
        iy = y+dy
        if ix>=1 && ix<=nx && iy>=1 && iy<=ny
            aa[iy,ix] = c
        end
    end
end

function initaa(ny, nx, ra)
    aa = zeros(ny, nx)
    for t in 0:((nx/ra)*(ny/ra))
        splat!(aa, ny, nx, ra);
    end
    aa
end

Сигмоиды

func_smooth(x::Float64, a, b)  = 1 / ( 1 + exp(-4(x-a)/b) ) 

sigmoid_a(x, a, ea) = func_smooth(x, a, ea)
sigmoid_b(x, b, eb) = 1 - sigmoid_a(x, b, eb)
sigmoid_ab(x, a, b, ea, eb) = sigmoid_a(x, a, ea) * sigmoid_b(x, b, eb)
sigmoid_mix(x, y, m, em) = x - x * func_smooth(m, 0.5, em) + y * func_smooth(m, 0.5, em)

function snm(N, M, en, em, b1, b2, d1, d2)

    [ sigmoid_mix( sigmoid_ab(N[i,j], b1, b2, en, en), 
                   sigmoid_ab(N[i,j], d1, d2, en, en), M[i,j], em ) 
                    for i = 1:size(N, 1), j = 1:size(N, 2) ]
end

Главная функция

function smoothlife(NX = 128, NY = 128, tfin = 10, scheme = 1)

    function derivative(aa)
        aaf = fft(aa)
        nf = aaf .* krf
        mf = aaf .* kdf
        n = real.(ifft(nf)) / kflr # irfft
        m = real.(ifft(mf)) / kfld
        2snm(n, m, alphan, alpham, b1, b2, d1, d2) .- 1
    end

    ra = 10
    ri = ra/3
    b = 1
    b1 = 0.257
    b2 = 0.336
    d1 = 0.365
    d2 = 0.551
    alphan = 0.028
    alpham = 0.147

    kd = zeros(NY,NX)
    kr = zeros(NY,NX)
    aa = zeros(NY,NX)

    x = [ j - 1 - NX/2 for i=1:NY, j=1:NX ]
    y = [ i - 1 - NY/2 for i=1:NY, j=1:NX ]

    r = sqrt.(x.^2 + y.^2)
    kd = 1 .- func_linear(r, ri, b) # ones(NY, NX)
    kr = func_linear(r, ri, b) .* ( 1 .- func_linear(r, ra, b) )
    #aa = snm (ix/NX, iy/NY, alphan, alpham, b1, b2, d1, d2);
    kflr = sum(kr)
    kfld = sum(kd)
    krf = fft(fftshift(kr))
    kdf = fft(fftshift(kd))

    #for td in 2 .^(10:-1:3)
    for td = 64

        aa = initaa(NY,NX,ra)
        dt = 1/td;
        l = 0
        nx = 0
        for t = 0:dt:tfin
# 1=euler  2=improved euler  3=adams bashforth 3  4=runge kutta 4
            if scheme==1

                aa += dt*derivative(aa)

            elseif scheme==2

                da = derivative(aa);
                aa1 = clamp(aa + dt*da)
                for h = 0:20
                    alt = aa1
                    aa1 = clamp(aa + dt*(da + derivative(aa1))/2)
                    if maximum(abs.(alt-aa1))<1e-8 
                        break
                    end
                end
                aa = copy(aa1)

            elseif scheme==3

                n0 = 1+mod(l,3)
                n1 = 1+mod(l-1,3)
                n2 = 1+mod(l-2,3)

                f = zeros(NY, NX, 3) # ?
                f[:,:,n0] = derivative(aa)
                if l==0
                    aa += dt*f[:,:,n0]
                elseif l==1
                    aa += dt*(3*f[:,:,n0] - f[:,:,n1])/2
                elseif l>=2
                    aa += dt*(23*f[:,:,n0] - 16*f[:,:,n1] + 5*f[:,:,n2])/12
                end

            elseif scheme==4

                k1 = derivative(aa)
                k2 = derivative(clamp(aa + dt/2*k1))
                k3 = derivative(clamp(aa + dt/2*k2))
                k4 = derivative(clamp(aa + dt*k3))
                aa += dt*(k1 + 2*k2 + 2*k3 + k4)/6

            end

            aa = clamp(aa)

            if t >= nx
                save("$(scheme)\$(td)_$t.png", Gray.(kron(aa, ones(2, 2) ) ) )
                nx += 1;
            end

            l += 1;
        end

    end     # for td
end

@time smoothlife(256, 256, 20, 3)

Julia и клеточные автоматы - 11

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

Julia и клеточные автоматы - 12

Тьюрмиты

Еще один интересный вид машин Тьюринга это тьюрмиты [21] — клеточные автоматы, состояние которых определяется неким агентом, который будет перемещаться по полю в зависимости от этих самых состояний. Наиболее известный представитель этого вида — муравей Лэнгтона [22].

Julia и клеточные автоматы - 13

На клеточное поле, каждая ячейка которого может быть либо черной либо белой, сажается мураш, который будет действовать по следующему алгоритму:

  • На чёрном квадрате — повернуть на $90^circ$ влево, изменить цвет квадрата на белый, сделать шаг вперед на следующую клетку.
  • На белом квадрате — повернуть на $90^circ$ вправо, изменить цвет квадрата на чёрный, сделать шаг вперед на следующую клетку.

На rosetta code есть довольно хитрая реализация [23], где направление муравья задается комплексным числом:

function ant(width, height)
    y, x = fld(height, 2), fld(width, 2)
    M = trues(height, width)

    dir = im
    for i in 0:1000000
        x in 1:width && y in 1:height || break
        dir *= M[y, x] ? im : -im
        M[y, x] = !M[y, x]
        x, y = reim(x + im * y + dir)
        i%100==0 && save("LR//zLR_$i.png", Gray.( kron(M,ones(4,4) ) ) )
    end

    Gray.(M)
end

ant(100, 100)

Julia и клеточные автоматы - 16

Поначалу, поведение агента кажется хаотичным. Но потом появляется магистраль — система выходит на стационар. Прекраснейший пример рождения порядка из кажимого хаоса!

Но самый смак начинается если расширить состояние ячейки за пределы одного бита:

Код

function ant(width, height, comnds;
        steps = 10, cxema = reverse(ColorSchemes.hot), 
        savevery = 100, pixfactor = 20)

    ma3x2pix() = [ clrs[k%n+1] for k in M ]
      bigpix() = kron( ma3x2pix(), ones(Int64,m,m) )
    save2pix(i) = save("$(comnds)_$i.png", bigpix() )

    m = pixfactor
    n = length(comnds)
    colorsinscheme = length(cxema)
    M = zeros(Int64, height, width)
    y, x = fld(height, 2), fld(width, 2)

    st = colorsinscheme ÷ n
    clrs = [ cxema.colors[i] for i in 1:st:colorsinscheme ]

    dir = im
    for i in 0:steps
        x in 1:width && y in 1:height || (save2pix(i); break)
        j = M[y, x] % n + 1
        dir *= comnds[j] == 'L' ? -im : im
        M[y, x] += 1
        x, y = reim(x + im * y + dir)

        i % savevery==0 && save2pix(i)
    end
    ma3x2pix()
end

@time ant(16, 16, "LLRR", steps = 100, savevery = 1, pixfactor = 20)

Julia и клеточные автоматы - 17

Послойными наращиваниями получаем симметричную структуру, похожую на мозг [24]! Для некоторых цветовых схем получается довольно футуристично

Julia и клеточные автоматы - 18

Julia и клеточные автоматы - 19

Правда потом он вырождается в жопу кардиоиду [25]. Кстати, ничего еще не напоминает? А если добавить фрактальных наростов?..

Код

# http://rosettacode.org/wiki/Mandelbrot_set
function hsv2rgb(h, s, v)
    c = v * s
    x = c * (1 - abs(((h/60) % 2) - 1) )
    m = v - c

    r,g,b =
        if h < 60
            (c, x, 0)
        elseif h < 120
            (x, c, 0)
        elseif h < 180
            (0, c, x)
        elseif h < 240
            (0, x, c)
        elseif h < 300
            (x, 0, c)
        else
            (c, 0, x)
        end

    (r + m), (b + m), (g + m)
end

function mandelbrot()

    w, h = 1000, 1000

    zoom  = 1.0
    moveX = 0
    moveY = 0

    img = Array{RGB{Float64}, 2}(undef,h, w)
    maxIter = 30

    for x in 1:w
        for y in 1:h
            i = maxIter
            c = Complex(
                (2*x - w) / (w * zoom) + moveX,
                (2*y - h) / (h * zoom) + moveY
            )
            z = c
            while abs(z) < 2 && (i -= 1) > 0
                z = z^2 + c
            end
            r,g,b = hsv2rgb(i / maxIter * 360, 1, i / maxIter)
            img[y,x] = RGB{Float64}(r, g, b)
        end
    end

    save("mandelbrot_set2.png", img)
end

mandelbrot()

Julia и клеточные автоматы - 20

и еще ряд интересных правил:

Julia и клеточные автоматы - 21

Julia и клеточные автоматы - 22

Магистрали довольно частое явление

Julia и клеточные автоматы - 23

Julia и клеточные автоматы - 24

Julia и клеточные автоматы - 25

Julia и клеточные автоматы - 26

Julia и клеточные автоматы - 27

Мозг [24] в ящике — довольно стабильная структура. Здесь результат тридцати миллиардов шагов:

Julia и клеточные автоматы - 28

Julia и клеточные автоматы - 29

Julia и клеточные автоматы - 30

Julia и клеточные автоматы - 31

Julia и клеточные автоматы - 32

Короче, с такого рода экспериментами можно засесть надолго. И это мы только комбинировали два направления (лево-право), а ведь можно добавить еще вперед-назад, к тому же, никто не говорил, что муравей обязан продвигаться на одну клетку! Держу пари, что такого рода расширение стандартных правил породит еще много интересных закономерностей!

Эпилог

К слову о жизни. Немало интереса вызывают самовоспроизводящиеся структуры, а ля циклы Лэнгтона [26]

Julia и клеточные автоматы - 33

Можно почитать статью [27], где данную модель расширяют наличием пола и передачей генов. (Тут уж передаю эстафету кому-нибудь другому. Очень хотелось бы посмотреть на реализацию подобной штуки)

Тема клеточных автоматов с ростом вычислительных мощностей и улучшением алгоритмов становится все популярней. У Вольфрама в блоге небольшой отчет [28] по этой теме, да и каждый может сам убедиться [29] — статей много, причем самых разнообразных: там вам и дружба с нейронными сетями, и генераторы случайных чисел, и квантовые точки, и сжатие и много всякого другого...

Ну и напоследок самокопирующаяся структура созданная Тэдом Коддом на спор за пинту пива [30]

Julia и клеточные автоматы - 34

Автор: Yermack

Источник [31]


Сайт-источник PVSM.RU: https://www.pvsm.ru

Путь до страницы источника: https://www.pvsm.ru/programmirovanie/348573

Ссылки в тексте:

[1] Golly: http://golly.sourceforge.net

[2] Lenia: https://github.com/Chakazul/Lenia

[3] Онлайн симулятор: http://lpcs.math.msu.su/~zolin/ant/

[4] решения задач: https://habr.com/ru/company/wolfram/blog/470425/

[5] здесь и там: https://habr.com/ru/post/471908/#comment_20786518

[6] вариант: https://github.com/YermolenkoIgor/Julia_tutorial_rus/blob/master/Toyz/automata.ipynb

[7] разбирали ранее: https://habr.com/ru/post/459842/

[8] моделирования мышления: https://habr.com/ru/post/308370/

[9] архитектуры тетриса: https://habr.com/ru/post/338584/

[10] самой себя: https://www.youtube.com/watch?time_continue=3&v=xP5-iIeKXE8&feature=emb_logo

[11] играбельно: https://copy.sh/life/?gist=f3413564b1fa9c69f2bad4b0400b8090&step=512

[12] все новые структуры: https://www.youtube.com/watch?feature=player_detailpage&v=C2vgICfQawE

[13] на гусеничном ходу!: https://www.conwaylife.com/wiki/Caterpillar

[14] подобными высказываниями: https://habr.com/ru/company/wolfram/blog/303836/

[15] обработкой изображений: https://habr.com/ru/company/yandex/blog/254249/

[16] БПФ: https://ru.wikipedia.org/wiki/%D0%91%D1%8B%D1%81%D1%82%D1%80%D0%BE%D0%B5_%D0%BF%D1%80%D0%B5%D0%BE%D0%B1%D1%80%D0%B0%D0%B7%D0%BE%D0%B2%D0%B0%D0%BD%D0%B8%D0%B5_%D0%A4%D1%83%D1%80%D1%8C%D0%B5

[17] хабровской статьи: https://habr.com/ru/post/180135/

[18] Репозиторий проекта: https://github.com/duckythescientist/SmoothLife

[19] Доступное объяснение: https://0fps.net/2012/11/19/conways-game-of-life-for-curved-surfaces-part-1/

[20] Дисер с доскональными исследованиями: https://www.csun.edu/~kme52026/thesis.html

[21] тьюрмиты: https://en.wikipedia.org/wiki/Turmite

[22] муравей Лэнгтона: https://ru.wikipedia.org/wiki/%D0%9C%D1%83%D1%80%D0%B0%D0%B2%D0%B5%D0%B9_%D0%9B%D1%8D%D0%BD%D0%B3%D1%82%D0%BE%D0%BD%D0%B0

[23] довольно хитрая реализация: http://rosettacode.org/wiki/Langton%27s_ant

[24] мозг: http://www.braintools.ru

[25] кардиоиду: https://ru.wikipedia.org/wiki/%D0%9A%D0%B0%D1%80%D0%B4%D0%B8%D0%BE%D0%B8%D0%B4%D0%B0

[26] циклы Лэнгтона: https://en.wikipedia.org/wiki/Langton%27s_loops

[27] статью: https://www.researchgate.net/publication/4249781_Sexyloop_Self-Reproduction_Evolution_and_Sex_in_Cellular_Automata

[28] небольшой отчет: https://writings.stephenwolfram.com/2012/05/its-been-10-years-whats-happened-with-a-new-kind-of-science/

[29] сам убедиться: https://arxiv.org/search/?query="cellular+automata"&searchtype=title&abstracts=show&order=-announced_date_first&size=50

[30] на спор за пинту пива: https://timhutton.github.io/2010/03/10/30984.html

[31] Источник: https://habr.com/ru/post/490454/?utm_source=habrahabr&utm_medium=rss&utm_campaign=490454