Создание образа Мона Лизы в Игре «Жизнь»

Клеточные автоматы представляют большой интерес и являются предметом исследования во многих областях, включая математику, физику, биологию, программирование и прочие. В статье мы разберем базовую реализацию и оптимизацию алгоритма для поиска состояния Жизни, из которого в течение нескольких поколений будет генерироваться образ Мона Лизы.

На загрузку этого видео может уйти несколько секунд. Рекомендую прищуриться – так удастся лучше разглядеть результат:)

Как видите, итог эксперимента получился не совсем запланированный, но я посчитал, что об этом все же стоит написать. Мой ум долго будоражила идея относительно Игры «Жизнь»: «Интересно, можно ли задействовать какой-нибудь стохастический алгоритм, дающий начальное состояние, которое бы через множество циклов формировало разборчивый текст».

Недавно мне попалась одноименная статья Кевина Галлигана, которая натолкнула меня на мысль, что можно добиться аналогичных результатов, но спомощью другого подхода. Что если вместо решателя задач выполнимости булевых формул применить определенный эвристический алгоритм, который сможет «запрограммировать» огромный мир Игры «Жизнь» так, чтобы спустя несколько поколений формировалось изображение?

Этого реально добиться другими путями. Например, поместить натюрморты (устойчивые фигуры) в определенных пикселях, как описано в вопросе на Codegolf.

Мой же замысел был отобразить картину Мона Лизы для одного кадра/поколения Жизни без натюрморта.

Алгоритм

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

Вот весь алгоритм:

best_score := infinity
    target := mona lisa with dimensions m x n
    canvas := random matrix of m x n
    best_result := canvas
    do
        modified_canvas := Copy of canvas with a single random cell inverted
        nth_modified_canvas := Run N generations of Game of Life modified_canvas
        Compute a score of how close nth_modified_canvas is with target
        if score < best_score then
        	best_score := score
            best_result := modified_canvas
        canvas := best_result
    while(max_iterations limit passed or best_score < threshold)

Изначально я создал прототип для выполнения на одном ядре.

def modify(canvas, shape):
    x,y = shape
    px = int(np.random.uniform(x+1))-1
    py = int(np.random.uniform(y+1))-1
    canvas[px][py] = not canvas[px][py]
    return canvas

def rmse(predictions,targets):
    return np.sqrt(np.mean((predictions-targets)**2))

while best_score>limit:
    canvases = np.tile(np.copy(best_seed), (batch_size, 1, 1))
    rms_errors = []
    for canvas in range(len(canvases)):
        canvases[canvas] = modify(states[state], (m,n))
        rmse_val = rmse(target, nth_generation(np.copy(canvases[canvas])))
        rms_errors.append(rmse_val)
    lowest = min(rms_errors)
    if lowest < best_score:
        best_score = lowest
        best_result = canvases[rms_errors.index(lowest)]

Принцип работы поиска с восхождением к вершине заключается в определении ближайшего к текущему соседнего состояния с наименьшей ошибкой по сравнению с target_state, представляющим Мона Лизу. Ближайший сосед на каждом шаге находится через создание копии лучшего текущего решения и инвертирования случайной клетки. Столь небольшое изменение не вызывает риска выхода из локального минимума. В качестве метрики для сравнения лучшего состояния с целевым используется среднеквадратичная ошибка (RMSE).

Спустя несколько дней CPU-вычислений я смог получить изображение, напоминающее Мона Лизу, через воспроизведение 4 поколений Жизни.

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

Препроцессинг

В качестве целевого изображения Мона Лизы, с которым алгоритм сравнивает случайное состояние, использовалась версия среднего разрешения, взятая из Википедии и преобразованная в черно-белую форму с помощью инструкции Python Imaging Library Image.open('target.png').convert('L').

При сравнении с логическими переменными лучше использовать бинарную матрицу, а не весь диапазон оттенков серого. Для этого я просто округлил значения оттенков серого до 0 и 1. Все же это было ошибкой, потому что в итоге утратилось много деталей.

Можно было вообще не округлять и сравнивать с версией в оттенках серого, но есть способ получше.

Состояния Сада Эдема

Не каждая случайная матрица 0 и 1 представляет допустимое состояние Игры «Жизнь». Состояния, которые никогда не могут быть n-поколением (n>0) любого клеточного автомата называются Сад Эдема.

Практически нереально, чтобы наша округленная черно-белая вариация Мона Лизы оказалась возможным поколением Жизни. Поэтому остается искать решение, которое окажется просто максимально приближенным к цели.

Это часть 4-го поколения только что подготовленного состояния

Проанализировав текстуру, способ развития паттернов Жизни и просто поэкспериментировав с изображениями, я выяснил, что если делать сравнение с прошедшей дизеринг 1-битной версией цели, то результат должен получаться лучше.

1-битная версия Мона Лизы

В таком изображении приблизительно равно распределены клетки 0 и 1, что в некоторой степени соответствует случайно инициализированному состоянию Игры, получаемому спустя несколько поколений. Это свойство также сохраняется при масштабировании изображения, которое мы вскоре реализуем.

Получить такой вариант изображения можно с помощью функции дизеринга Флойда-Стейнберга библиотеки PIL, использовав инструкцию Image.open('target.png').convert('1').

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

Векторизация с помощью JAX

Версия для одного ядра CPU без использования векторизации оказывается чрезвычайно медленной. Я пробовал выполнять ее и на своем Core i7 восьмого поколения, и на машинах Google Colab, но для получения похожего на цель результата приходится всякий раз ждать часы или даже дни, в зависимости от разрешения.

К счастью, эта задача отлично поддается распараллеливанию.

JAX – это библиотека Python, которая позволяет использовать numpy-версию и компилировать ее в высоко векторизованный код, который можно выполнять на графических и тензорных процессорах. В нашем случае нужно переработать этот алгоритм конкретно для графического (GPU).

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

Мы экструдируем target (Мона Лиза) и canvas (начальное случайное состояние) в 3-е измерение, в котором они формируют многослойные тензорные структуры (loafs) с длиной равной batch_size.




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

Дополнительно для каждой итерации этого цикла нужно производить случайный тензор, называемый мутатором (той же формы, что target), со следующим свойством: в каждом слое должна присутствовать только одна 1, расположенная в случайном месте.

Вот пример мутатора формы 5, 3, 2, где 5 — это размер batch-size:

array([[[1, 0],
        [0, 0],
        [0, 0]],

       [[1, 0],
        [0, 0],
        [0, 0]],

       [[0, 0],
        [0, 1],
        [0, 0]],

       [[0, 1],
        [0, 0],
        [0, 0]],

       [[1, 0],
        [0, 0],
        [0, 0]]])

Идея в том, чтобы в каждом цикле с помощью мутатора вычислять ближайший набор соседних с best_canvas состояний: canvas = (best_canvas + mutator)%2.

Мы вычисляем N поколений Игры в каждом слое этого модифицируемого холста, после чего определяем показатель RMSE (среднее вычисляется только для слоя) холста N-поколения в отношении к Мона Лизе и находим слой с наименьшей ошибкой. Далее этот слой экструдируется и устанавливается как best_canvas, после чего данный цикл повторяется в течение конечного числа итераций.

Код

Весь код этого проекта доступен на GitHub. В текущем же разделе я поясню, что в нем означает каждый блок. Если вы хотите сразу увидеть результаты, просто промотайте к концу статьи. Основа проекта, а именно функция Игры «Жизнь» взята из этого поста.

Выражаю благодарность Бояну Николичу. Я следовал его условию импорта jax.numpy в качестве N, а jax.lax в качестве L.

%matplotlib inline 
import jax
N=jax.numpy
L=jax.lax
from jax.experimental import loops
from jax import ops
import matplotlib.pyplot as plt
import numpy as onp
import time
from PIL import Image 
from google.colab import files

Далее выполняем wget Мона Лизы:

!wget -O target.png https://upload.wikimedia.org/wikipedia/commons/thumb/e/ec/Mona_Lisa%2C_by_Leonardo_da_Vinci%2C_from_C2RMF_retouched.jpg/483px-Mona_Lisa%2C_by_Leonardo_da_Vinci%2C_from_C2RMF_retouched.jpg?download

Это скромная версия шириной в 483px.

batch_size = 100
image_file = Image.open("target.png")
image_file = image_file.convert('1')
lisa = N.array(image_file, dtype=N.int32)
width,height = lisa.shape
lisa_loaf = onp.repeat(lisa[onp.newaxis, :, :,], batch_size, axis = 0)

Здесь выполняется дизеринг изображения Мона Лизы и его экструдирование до глубины равной batch_size:

key = jax.random.PRNGKey(42)
canvas_loaf = jax.random.randint(key, (batch_size, width, height), 0, 2, dtype= N.int32) 
#для тестов инициализируем случайное изображение Лизы

Здесь создается начальное значение JAX PRNG (вскоре я поясню, о чем речь), а также начальный случайный canvas_loaf с целыми числами 0 и 1.

@jax.jit
def rgen(a):
    # Это сокращение превосходит численностью соседей живых клеток, так как включает в себя #саму центральную клетку. Чтобы это исправить, нужно вычесть массив
    nghbrs=L.reduce_window(a, 0, L.add, (3,3), (1,1), "SAME")-a
    birth=N.logical_and(a==0, nghbrs==3)
    underpop=N.logical_and(a==1, nghbrs<2)
    overpop=N.logical_and(a==1, nghbrs>3)
    death=N.logical_or(underpop, overpop)

    na=L.select(birth,
                N.ones(a.shape, N.int32),
                a)

    na=L.select(death,
                N.zeros(a.shape, N.int32),
                na)
    return na

vectorized_rgen = jax.vmap(rgen)

@jax.jit
def nv_rgen(state):
  for _ in range(n_generations):
      state = vectorized_rgen(state)
  return state

Пояснение функции rgen, которая воспроизводит одно поколение Жизни, можете найти в статье упомянутого мной Бояна Николича.

jax.vmap позволяет создать функцию, которая отображает входную функцию на оси аргументов (выполняет векторизацию). Это дает возможность воспроизвести поколение Жизни для каждого слоя холста.

nv_rgen воспроизводит на холсте N-поколений Жизни.

Помимо этого, декоратор @jax.jit обуславливает JIT-компиляцию данной функции (Just-in-time, то есть компиляцию в процессе выполнения программы). Не уверен, произошли ли в данном случае какие-либо улучшения, так как nv_rgen просто складывается из других функций, также компилируемых при выполнении.

def mutate_nj(b, w, h, subkey):
  a = jax.random.normal(subkey, (b, w, h))
  return (a == a.max(axis=(1,2))[:,None,None]).astype(int)

mutate = jax.jit(mutate_nj, static_argnums=(0,1,2))

mutate_nj (nj означает без JIT) генерирует тензор мутатора с помощью jax.random.normal, устанавливая для каждого слоя максимум одну 1 и остальные 0. Аргумент subkey я объясню чуть позже.

Далее мы JIT-компилируем эту функцию как mutate, а также отмечаем аргументы b, w, h как статические, указывая компилятору, что в процессе выполнения они изменяться не будут.

def rmse_nj(original, canvas, b, w, h):
  return N.sqrt(N.mean(L.reshape((original-canvas)**2,(b,w*h)) , axis=1))

rmse = jax.jit(rmse_nj, static_argnums=(2,3,4))

rmse говорит сама за себя. Единственное существенное отличие от CPU-версии в том, что среднее вычисляется по 1-й оси (оси длины экструдированной структуры).

def hill_climb(original, canvas, prng_key, iterations):
  with loops.Scope() as s:
    s.best_score = N.inf
    s.best_canvas = canvas
    s.canvas = canvas
    s.prng_key = prng_key
    for run in s.range(iterations):
      s.prng_key, subkey = jax.random.split(s.prng_key)
      s.canvas+=modify(batch_size, width, height, subkey)
      s.canvas%=2
      rmse_vals = rmse(original, nv_rgen( s.canvas ), batch_size, width, height)
      curr_min = N.min(rmse_vals)
      for _ in s.cond_range(curr_min < s.best_score):
        s.best_score = curr_min
        s.best_canvas = N.repeat((s.canvas[N.argmin(rmse_vals)])[N.newaxis, :, :,], batch_size, axis = 0)
      s.canvas = s.best_canvas
    return s.canvas

hill_climb является основной функцией программы. Это одна большая конструкция цикла JAX. Можно было использовать здесь стандартные циклы Python, но нужно по максимуму задействовать возможности JAX.

Циклы в JAX (модуль jax.experimental.loops) – это функции, представляющие синтаксический сахар, например lax.fori_loop и lax.cond. Циклы lax (фактически, циклы XLA), которые содержат больше нескольких инструкций и вложения, становятся очень сложными. Тем не мене в JAX они приближены к стандартным циклам Python. Единственная загвоздка в том, что состояние цикла, то есть все, что изменяется по ходу итераций, должно храниться как член его области. В нашем случае сюда относятся best_score, best_canvas, временный холст, где мы воспроизводим Жизнь, а также ключ PRNG.

JAX PRNG

В NumPy для всех функций, предполагающих наличие случайных чисел, используется управляемый генератор псевдослучайных чисел (PRNG). То есть NumPy полностью отвечает за определение для него начального числа (seed) и последующую обработку состояния. Как я понимаю, при параллельном выполнении, а также в ситуациях, требующих большого количества случайных чисел, этот метод имеет свои недостатки. С его помощью сложно обеспечить достаточную энтропию для производства нужного количества случайных значений.

В отличие от NumPy, JAX не управляет случайной генерацией чисел. Каждой функции jax.random в качестве первого аргумента необходимо предоставлять текущее состояние PRNG, и при каждом выполнении одной из этих функций состояние PRNG должно обновляться через jax.random.split.

Если состояние не обновлять, то это быстро приведет к повторяющейся генерации одинаковых наборов случайных чисел. Я не до конца понимал этот момент, когда первый раз писал цикл, и в результате алгоритм перестал находить новые вариации состояния холста. Причина в том, что генерируется одинаковый тензор-мутатор.

Генерацию каждым параллельным компонентом алгоритма различных случайных чисел можно также обеспечить с помощью разделения состояния PRNG. Подробнее о структуре PRNG в JAX можно почитать здесь.

cond_range

Почему в JAX условные конструкции также являются циклами? Честно, я не совсем это понимаю. У cond_range должна быть возможность выводить стандартное логическое значение вместо итератора длиной 0/1, но по какой-то причине спроектирована конструкция именно так.

Если мы находим более соответствующий слой холста, то экструдируем его и устанавливаем как best_canvas, а его показатель ошибки как best_score.

Спустя конечное число итераций мы получаем состояние Жизни, которое через N поколений формирует образ Мона Лизы.

Итоги

Выполнение ~1000 итераций для изображения шириной 483px на GPU Google Colab занимает всего ~40 секунд. Если сравнивать с версией программы для CPU, которая даже изображение меньшего разрешения обрабатывала несколько часов, можно считать, что цели мне достичь удалось.

Состояние Жизни, дающее наибольшее сходство с целевым изображением, было достигнуто через ~23 000 итераций, что заняло 10 минут. После 23 000 повторений прирост качества прекращается, и даже при достижении 100 000 итераций видимых улучшений не наблюдается.
Также отмечу, что изображения, получаемые при меньшем числе поколений, как и ожидалось, проявляют большее сходство с целью.

Мона Лиза, 10 поколений:

Тестовый образец шахматной доски, 7 поколений:

Тестовый образец текста, 5 поколений:

Давид (работа Микеланджело), 3 поколения:

Луна, 7 поколений (https://unsplash.com/photos/pd4lo70LdbI):

Нил Армстронг, 7 поколений:

Заключение

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

Джон Хортон Конвей (26 декабря 1937 — 11 апреля 2020), RIP: