Когда смотришь на рой светлячков ночью, возникает довольно естественный вопрос: как они вообще понимают, когда нужно вспыхнуть? Почему через некоторое время хаотичные вспышки вдруг начинают напоминать согласованное поведение?
На самом деле здесь быстро приходим к классической задаче о взаимодействии множества тел. У каждого светлячка есть собственный внутренний цикл, но при этом он ещё и реагирует на соседей. Если подобрать простую модель, то всё это начинает выглядеть почти как клеточный автомат – только «живой».
Идея этого небольшого эксперимента появилась после поста в блоге Cian Luke Martin — Fireflies, Magnets and Emergence. Мне захотелось повторить этот подход и посмотреть, получатся ли сверкающие жучки из этого и emergent behaviour.
Жучок
Если упростить задачу до минимума, каждому светлячку нужны всего две вещи:
-
внутреннее состояние, или «часы»;
-
вспышка, которую видят остальные.
Когда внутренний цикл доходит до конца, светлячок вспыхивает и тем самым немного подталкивает часы соседей к синхронизации:

TL;DR
Если хочется сразу увидеть конечную идею, то в самом сжатом виде она выглядит так:
Module[{f = RandomReal[{0,1.01}, {25,25}]}, Refresh[
f = Map[Clip[# - 0.01, {0.0,1.0}, {1.0,1.0}]&, f, {2}];
f = Clip[
f (2.0 Round[f] - 1.0)
Clip[
ListConvolve[
{
{1.0,1.0,1.0},
{1.0,0.0,1.0},
{1.0,1.0,1.0}
},
Floor[f],
2,
0
],
{0, 0.001}
],
{0., 1.0}
];
ArrayPlot[f, Frame->True, ColorFunction->"PlumColors"]
, 1/30.0]]
Его можно запустить без полпинка в WLJS Notebook, и с полпинка в Mathematica, если хочется. Если не хочется, то не нужно.
Примем два простых допущения:
-
длительность цикла у всех светлячков одинаковая;
-
если светлячок видит чужую вспышку, он немного подстраивает свою фазу – вперёд или назад, смотря что ближе.
Изначально все находятся не в фазе. Дальше локальные взаимодействия постепенно начинают выравнивать систему за каждую вспышку.
N светлячков
Сначала просто раскидаем точки по области и назначим каждой случайное состояние:
generateFireFlies[n_:200, region_:Rectangle[{-10,-10}, {10,10}]] := fireFlies = With[{
pos = RandomPoint[region, n]
},
Transpose[{pos, Table[RandomReal[{0,1.0}], {Length[pos]}]}]
];
Каждый светлячок здесь – это список из двух элементов:
-
position -
state
Например:
generateFireFlies[1] // First
вывод:
{{-5.896003641067651, 1.7548145516864455}, 0.1972596724287563}
Запекаем связи между соседями
Теперь нужно понять, кто кого вообще «видит». Для этого заранее вычислим связи между точками в заданном радиусе, полагая, что жуки устали и не двигаются:
bakeConnections[r_:3.7] := (
connectionMatrix = Table[
If[i == j, Infinity, Power[Norm[i[[1]] - j[[1]]], 2]]
, {i, fireFlies}, {j, fireFlies}];
connectionMatrix = MapIndexed[
Function[{value, index},
If[value < r, index[[1]], Nothing]
],
#
] & /@ connectionMatrix
);
После этого у каждого светлячка есть список соседей, на которых он будет реагировать
generateFireFlies[10, Disk[{0,0},10]];
bakeConnections[30.0];
Graphics[{
Blue, Table[
Line[{fireFlies[[i,1]], fireFlies[[#,1]]}]&/@connectionMatrix[[i]],
{i,Length[fireFlies]}
],
Red, Disk@@#&/@fireFlies
}, ImageSize->250]

Дальше всё поведение можно разложить на три функции:
-
затухание;
-
вспышка;
-
подстройка под соседей.
decay[{pos_, state_}] := {pos, Clip[state - 0.01, {0,1}]};
flash[{pos_, state_}] := {pos, If[state == 0, 1, state]};
adjust[{pos_, 0}, index_] := {pos, 0};
adjust[{pos_, state_}, index_] := With[{i = index[[1]]},
{
pos,
If[
Or @@ Thread[fireFlies[[connectionMatrix[[i]], 2]] == 1],
Clip[state Sign[2 Round[state] - 1] 0.0001, {0,1}],
state
]
}
];
Самая интересная часть — это adjust.
Если кто-то из видимых соседей находится в состоянии вспышки (1.0), то мы слегка подталкиваем текущий внутренний таймер. Причём не всегда в одну сторону: если текущее состояние ближе к началу цикла, то логично откатывать его назад, а если ближе к концу - толкаем вперёд. Именно за это отвечает выражение Sign[2 Round[state] - 1].
Теперь собираем один шаг симуляции:
update = Function[Null,
fireFlies = decay /@ fireFlies;
fireFlies = flash /@ fireFlies;
fireFlies = MapIndexed[adjust, fireFlies];
];
2 жучка
Как и в случае с любой подобной моделью, сначала полезно проверить поведение на совсем маленьком числе жуков
generateFireFlies[2, Circle[{0,0},0.5]];
bakeConnections[2.0];
radii = fireFlies[[All,2]];
history = {Table[{i, 0}, {i,200}], Table[{i, 0}, {i,200}]};
{Graphics[{
Circle[fireFlies[[1]][[1]], radii[[1]]//Offload],
Circle[fireFlies[[2]][[1]], radii[[2]]//Offload],
EventHandler[AnimationFrameListener[radii//Offload], Function[Null,
update[];
update[];
radii = 0.5 fireFlies[[All,2]];
history[[1,1,2]] = fireFlies[[1,2]];
history[[2,1,2]] = fireFlies[[2,2]];
history = {
Transpose[{history[[1,All,1]], RotateLeft[history[[1,All,2]], 1]}],
Transpose[{history[[2,All,1]], RotateLeft[history[[2,All,2]], 1]}]
};
]]
}, TransitionType->None, "Controls"->False,
ImageSize->250,
PlotRange->{{-1,1}, {-1,1}}], Graphics[{
ColorData[97][1], Line[history[[1]] // Offload],
ColorData[97][2], Line[history[[2]] // Offload]
}, Axes->True, Frame->True, "Controls"->False,
PlotRange->{{0,200}, {0,1}},
TransitionType->None,
ImageSize->250
]}//Row
Здесь мы также трекаем их состояния и строим график того, как они меняются со временем

Как видно для полной синхронизации нужно десяток циклов
Увеличиваем масштаб
С двумя объектами всё хорошо видно, но хочется посмотреть, что будет, если их уже не два, а, скажем, двести.
Для такой сцены круги удобнее заменить на цветовые вспышки:
цветовая функция
cf = Blend[{Darker[Red], Yellow}, #] &;
А дальше визуализировать поле светлячков как набор цветных дисков на чёрном фоне:
generateFireFlies[200, Disk[{0,0},10]];
bakeConnections[3.7];
colors = cf /@ fireFlies[[All,2]];
Graphics[{
Table[
With[{i = i, xy = fireFlies[[i,1]]},
{
RGBColor[colors[[i]]],
Disk[xy, 0.3]
} // Offload
],
{i, Length[fireFlies]}
],
EventHandler[
AnimationFrameListener[colors // Offload],
Function[Null,
update[];
colors = cf /@ fireFlies[[All,2]];
]
]
}, Background->Black]
Для большего «реализма» поверх можно продублировать слой и слегка размыть его - так появляется приятный светящийся ореол:

Клеточный автомат
На этом месте можно заметить, что сама логика обновления очень похожа на клеточный автомат в роде этого:
ArrayPlot[#, ImageSize -> 40, Mesh -> True] & /@
CellularAutomaton["GameOfLife", {{{0,1,0},{0,0,1},{1,1,1}}, 0}, 8]

Если раскидать жуков по регулярной сетке, то задачу можно разложить на параллельные операции над массивами. Нужно думать как “шейдер”!
Допустим, у нас есть поле:
field = RandomReal[{0,1.01}, {25,25}];
Чтобы определить клетки, находящиеся в состоянии вспышки, достаточно округлить значения вниз:
Floor[field]
или наглядно
ArrayPlot[
field Floor[field],
ColorFunction->Function[Piecewise[{{(RGBColor[1, 1, 0]),#1>0.7`},{GrayLevel[#1],True}}]],
PlotLabel->"Cells in a flashing state",
Frame->True
]

Дальше поиск вспыхнувших соседей превращается в обычную двумерную свёртку с клиппингом:
c = Clip[
ListConvolve[
{
{1.0,1.0,1.0},
{1.0,0.0,1.0},
{1.0,1.0,1.0}
},
Floor[field],
2,
0
],
{0, 0.01}
];
построим:
ArrayPlot[
field 300 c,
ColorFunction->Function[Piecewise[{{(RGBColor[1, 0, 0]),#1>0.7`},{GrayLevel[#1],True}}]],
PlotLabel->"Cells affected by flashes",
Frame->True
]

Красота в том, что здесь мы уже не итерируем по жучкам процедурно. Вместо этого используем простой и очень дешёвый (ну не прям дешевый, но JIT-friendly) для массива оператор. Это масштабируется гораздо лучше.
Остаётся последний шаг: для каждой клетки определить знак коррекции. Его можно вычислить опять операциями над массивами:
c = (2.0 Round[field] - 1.0) * c;
Если значение ближе к одному краю цикла - поправка получается отрицательной, если к другому - положительной. Покажем это цветами:
ArrayPlot[
field 3000 c,
ColorFunction->Function[Piecewise[{{(RGBColor[1, 0, 0]),#1>0.7`},{(RGBColor[0, 0, 1]),#1<0.3`},{GrayLevel[#1],True}}]],
PlotLabel->"Cells affected by flashes (forward and backward)",
Frame->True
]

А затухание и вспышка по достижении нуля реализуются совсем просто:
f = Map[Clip[# - 0.01, {0.0,1.0}, {1.0,1.0}]&, f, {2}];
Финальная версия
Теперь можно собрать всё вместе:
Module[{f = RandomReal[{0,1.01}, {25,25}]}, Refresh[
f = Map[Clip[# - 0.01, {0.0,1.0}, {1.0,1.0}]&, f, {2}];
f = Clip[
f (2.0 Round[f] - 1.0)
Clip[
ListConvolve[
{
{1.0,1.0,1.0},
{1.0,0.0,1.0},
{1.0,1.0,1.0}
},
Floor[f],
2,
0
],
{0, 0.001}
],
{0., 1.0}
];
ArrayPlot[f, Frame->True, ColorFunction->"PlumColors"]
, 1/30.0]]

Через несколько минут итераций начинают проявляться устойчивые волновые структуры

Из-за того, что для соседей здесь используется прямоугольное ядро свёртки, и сами паттерны получаются слегка квадратными. В реальности, конечно, светлячки не живут на идеальной сетке, да и циклы у них могут отличаться. Но как виртуальный эксперимент на тему коллективного поведения - получается вполне себе наглядно.
Оптимизированная версия
Refresh вместе с ArrayPlot отлично подходят для прототипирования, но на каждом кадре создают лишние накладные расходы. Если хочется выжать из этого побольше производительности, лучше перейти к Image и напрямую кормить GPU текстурой через NumericArray.
field = RandomReal[{0,1.0}, {50,50}];
render = Function[Null,
Do[
field = Map[Clip[# - 0.01, {0.0,1.0}, {1.0,1.0}]&, field, {2}];
field = Clip[
field (2.0 Round[field] - 1.0)
Clip[
ListConvolve[
{
{1.0,1.0,1.0},
{1.0,0.0,1.0},
{1.0,1.0,1.0}
},
Floor[field],
2,
0
],
{0, 0.001}
],
{0., 1.0}
];
, {2}];
imageBuffer = NumericArray[255.0 field, "Byte", "ClipAndRound"];
];
render[];
Image[
imageBuffer // Offload,
"Byte",
Epilog -> EventHandler[
AnimationFrameListener[imageBuffer // Offload],
render
],
Magnification -> 20,
Antialiasing -> False
]

На больших числах итераций там начинают возникать уже совсем красивые структуры - почти как экран радара…

Рендерим видео
Раз уж мы уже работаем с изображением, собрать из этого видео совсем несложно:
val = 1.0;
PrintTemporary[ProgressIndicator[val // Offload, {1, 60 10}]];
movie = Table[
val ;
render[];
ImageResize[
Colorize[Image[imageBuffer]],
Scaled[6],
Method->"NearestNeighbor"
],
{60 10}
];
И потом экспортируем в видео:
FrameListVideo[movie, FrameRate->60]
Мне нравится эта модель тем, что она стоит где-то посередине между “агентной” (моделируем состояние отдельных жуков) симуляцией и клеточным автоматом. Или можно назвать это неклассическим клеточным автоматом, где состояния клеток меняются кусочно-непрерывно, но сама сетка всё еще дискретная. Кстати в тему непрерывных клеточных автоматов - посмотрите на Smooth Life.
Ссылки
-
Оригинальный пост на буржуйском: Fireflies or Nature’s Cellular Automaton
-
Текст, который вдохновил этот пост: Cian Luke Martin — Fireflies, Magnets and Emergence
-
WLJS Notebook - среда разработки: wljs.io
-
Wolfram Engine - интерпретатор: wolfram.com/engine


