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

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

Идея этого небольшого эксперимента появилась после поста в блоге 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, если хочется. Если не хочется, то не нужно.

Примем два простых допущения:

  1. длительность цикла у всех светлячков одинаковая;

  2. если светлячок видит чужую вспышку, он немного подстраивает свою фазу - вперёд или назад, смотря что ближе.

Изначально все находятся не в фазе. Дальше локальные взаимодействия постепенно начинают выравнивать систему за каждую вспышку.

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.

Ссылки