Внедрение Quadtree в Mathematica

Ответ 1

Вот более компактная версия. Он использует ту же структуру данных, что и исходная версия. Функции splitBox и insideBox по существу одинаковы (просто написано несколько иначе).

Вместо добавления точек поочередно исходное поле содержит все точки в начале, поэтому нет необходимости в подпрограммах qtInsert. На каждом этапе рекурсии ящики, содержащие более одной точки, разделяются, а точки распределяются по вспомогательным полям. Это означает, что все узлы с более чем одной точкой являются листами, поэтому нет необходимости также проверять это.

qtMakeNode[bb_, pts_] := {{}, {}, {}, {}, qtbb @@ bb, pts}

splitBox[bx_] := splitBox[{min_, max_}] := {min + #, max + #}/2 & /@  
  Tuples[Transpose[{min, max}]]


insideBox[pt_, bb_] := bb[[1, 1]] <= pt[[1]] <= bb[[2, 1]] && 
  bb[[1, 2]] <= pt[[2]] <= bb[[2, 2]]

distribute[qtree_] := Which[
  Length[qtree[[6]]] == 1, 
    (* no points in node -> return node unchanged *)
  qtree,

  Length[qtree[[6]]] == 1, 
    (* one point in node -> replace head of point with qtpt and return node *)
  ReplacePart[qtree, 6 -> qtpt @@ qtree[[6, 1]]],

  Length[qtree[[6]]] > 1, 
    (* multiple points in node -> create sub-nodes and distribute points *)
    (* apply distribute to sub-nodes *) 
  Module[{spl = splitBox[qtree[[5]]], div, newtreelist},
   div = Cases[qtree[[6]], a_ /; insideBox[a, #], 1] & /@ spl;
   ReplacePart[qtree, 
    Join[Table[i -> distribute[qtMakeNode[spl[[i]], div[[i]]]], {i, 4}], 
      {6 -> {}}]]]]

Пример (с использованием исходной версии qtDraw):

len = 50;
pts = RandomReal[{0, 2}, {len, 2}];
qt = makeTree[qtMakeNode[{{0.0, 0.0}, {2.0, 2.0}}, pts]];
qtDraw[qt]

Результат:

Ответ 2

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

Как отмечали другие, возможно, было бы еще лучше использовать обертку Hold и/или атрибуты HoldXXX, чтобы эмулировать call-by-reference.

Для жесткого базового подхода к некоторым связанным реализациям структуры данных см.

http://library.wolfram.com/infocenter/MathSource/7619/

Соответствующий код находится в записной книжке Hemmecke-final.nb(назван так потому, что он реализует торический алгоритм Грёбнера из-за Р. Хеммеке и соавторов).

Я взял удар, переопределяя атрибуты Hold..., но я не очень хорош в этом и отказался от этого, когда код набросился на меня (пропустил, но убил сессию Mathematica). Таким образом, вместо этого у меня есть реализация, которая использует недокументированный "сырой" тип данных Mathematica, который является инертным и, следовательно, поддается управлению позором.

Рассматриваемая структура называется "expr bag", потому что общая структура данных Mathematica является "expr". Это похоже на список, но (1) он может расти с одного конца (хотя и не сжиматься) и (2), как и другие типы сырых выражений (например, графики в версии 8), он имеет компоненты, к которым можно получить доступ и/или изменить с помощью предоставленных функций (API, так сказать). Его лежащие в основе "элементы" являются инертными в том смысле, что они могут ссылаться на ANY expr (включая сам пакет) и могут управляться способами, которые я укажу ниже.

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

Я сохранил код более или менее в том же стиле, что и оригинал, и, в частности, он остается онлайновой версией (то есть, элементы не все должны входить с самого начала, но могут добавляться отдельно). Изменено несколько имен. Создал базовую структуру, сродни

node (ограничивающий прямоугольник, значение, ноль или четыре поднода)

Если есть подносы, то поле значения пусто. Поля box и value представлены обычным выражением Mathematica List, хотя может иметь смысл использовать выделенные главы и иметь более сродни стилю C struct. Я сделал что-то подобное, назвав различные функции доступа/настройки поля.

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

AppendTo[$ContextPath, "Internal`"];

makeQuadTreeNode[bounds_] := Bag[{bounds, {}, {}}]

(*is pt inside box?*)

insideBox[pt_, box_] := 
 And @@ Thread[box[[1]] <= (List @@ pt) <= box[[2]]]

(*split bounding box into 4 children*)

splitBox[{{xmin_, ymin_}, {xmax_, ymax_}}] := 
 Map[makeQuadTreeNode, {{{xmin, (ymin + ymax)/2}, {(xmin + xmax)/2, 
     ymax}}, {{xmin, 
     ymin}, {(xmin + xmax)/2, (ymin + ymax)/2}}, {{(xmin + xmax)/2, 
     ymin}, {xmax, (ymin + ymax)/2}}, {{(xmin + xmax)/
      2, (ymin + ymax)/2}, {xmax, ymax}}}]

bounds[qt_] := BagPart[qt, 1]
value[qt_] := BagPart[qt, 2]
children[qt_] := BagPart[qt, 3]

isLeaf[qt_] := value[qt] =!= {}
isSplit[qt_] := children[qt] =!= {}
emptyNode[qt_] := ! isLeaf[qt] && ! isSplit[qt]

(*qtInsert #1-return input if pt is out of bounds*)

qtInsert[qtree_, pt_] /; ! insideBox[pt, bounds[qtree]] := qtree

(*qtInsert #2-empty node (no value,no children)*)

qtInsert[qtree_, pt_] /; emptyNode[qtree] := value[qtree] = pt

(*qtInsert #2-currently a leaf (has a value and no children)*)

qtInsert[qtree_, pt_] /; isLeaf[qtree] := Module[
  {kids = splitBox[bounds[qtree]], currval = value[qtree]},
  value[qtree] = {};
  children[qtree] = kids;
  Map[(qtInsert[#, currval]; qtInsert[#, pt]) &, kids];
  ]

(*qtInsert #4-not a leaf and has children*)

qtInsert[qtree_, pt_] := Map[qtInsert[#, pt] &, children[qtree]];

getBoxes[ee_Bag] := 
 Join[{bounds[ee]}, Flatten[Map[getBoxes, children[ee]], 1]]
getPoints[ee_Bag] := 
 Join[{value[ee]}, Flatten[Map[getPoints, children[ee]], 1]]

qtDraw[qt_] := Module[
  {pts, bboxes},
  pts = getPoints[qt] /. {} :> Sequence[];
  bboxes = getBoxes[qt];
  Graphics[{EdgeForm[Black], Hue[0.2], Map[Disk[#, 0.01] &, pts], 
    Hue[0.7], EdgeForm[Red], 
    FaceForm[], (Rectangle @@ #) & /@ bboxes}, Frame -> True]]

Вот пример. Отмечу, что масштабирование разумно. Возможно O (n log (n)) или так. Определенно лучше, чем O (n ^ 2).

len = 4000;
pts = RandomReal[{0, 2}, {len, 2}];
qt = makeQuadTreeNode[{{0.0, 0.0}, {2.0, 2.0}}];
Timing[Do[qtInsert[qt, pts[[i]]], {i, 1, len}]]

{1.6, Null}

Общие примечания для заметок. Они старые, поэтому я не утверждаю, что все это работает, как указано.

Эти функции живут во внутреннем контексте.

сумка   Создает пакет expr, необязательно с предустановленными элементами.

BagPart   Получает части сумки expr, похожие на Часть для обычных exprs. Также может использоваться на lhs, например. до reset значение.

StuffBag   Добавляет элементы в конец сумки.

У нас также есть BagLength. Полезно для итерации по сумке.

Эти функции чрезвычайно полезны по двум причинам.

Во-первых, это хороший способ сделать расширяемую таблицу в Mathematica.

Во-вторых, содержимое мешков оценивается, но затем помещается в необработанный expr, следовательно, экранированы. Таким образом, можно использовать их как "указатели" (в смысле С), а не как объекты, и это не требует Hold и т.д. Вот несколько примеров:

a = {1,2,a} (* gives infinite recursion *)

Если мы вместо этого используем сумки, мы получаем самореферентную структуру.

In[1]:= AppendTo[$ContextPath, "Internal`"];

In[2]:= a = Bag[{1,2,a}]
Out[2]= Bag[<3>]

In[3]:= expr1 = BagPart[a, All]
Out[3]= {1, 2, Bag[<3>]}

In[4]:= expr2 = BagPart[BagPart[a, 3], All]
Out[4]= {1, 2, Bag[<3>]}

In[5]:= expr1 === expr2
Out[5]= True

Трудно подражать каким-либо другим способом в Mathematica. Нужно будет использовать разреженные таблицы (хеширование) в некоторых не очень прозрачный способ.

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

tail[ll_] := BagPart[ll,2]
settail[ll_, ll2_] := BagPart[ll,2] = ll2
contents[ll_] := BagPart[ll,1]
setcontents[ll_, elem_] := BagPart[ll,1] = elem

createlinkedlist[elems__] := Module[
    {result, elist={elems}, prev, el},
    result = Bag[{elist[[1]],Bag[]}];
    prev = result;
    Do [el = Bag[{elist[[j]],Bag[]}];
        settail[prev, el];
        prev = el,
        {j,2,Length[elist]}];
    result
    ]

In[18]:= tt = createlinkedlist[vv,ww,xx]
Out[18]= Bag[<2>]

In[20]:= BagPart[tt,All]
Out[20]= {vv, Bag[<2>]}

Итак, tt - связанный список, первый элемент - vv, следующий - сам связанный список и т.д. Я воздержался от использования Lisp терминологии (car/cdr и т.п.), потому что я не могу вспомните, являются ли операции списка Lisp разрушительными. Но вы получите общую идею.

В аналогичных строках я использовал пакеты expr для реализации двоичных деревья. Это полезно, потому что мы можем совершать деструктивные изменения в постоянное время (при условии, что у нас уже есть "дескриптор" в точке вставки/удаления) и, кроме того, "сырой" характер выражения мешки означает, что мы полностью избегаем бесконечной оценки Mathematica семантика.

Другое приложение, возможно.

Pointer = Internal`Bag
Contents[aa_Pointer, j_Integer] /;0<j<=Internal`BagLength[aa] :=
    Internal`BagPart[aa,j]
SetContents[aa_Pointer, j_Integer, e_] /; 0<j<=Internal`BagLength[aa] :=
    Internal`BagPart[aa,j] = e
SetContents[aa_Pointer, j_Integer, e_] /; j>BagLength[aa] :=
    (Do[Internal`StuffBag[aa,Null], {k,Internal`BagLength[aa]+1,j-1}];
    Internal`StuffBag[aa,e])

Попробуйте

a = Bag[{1,2,a,6,t,y,99,Bag[{a,q,3,r,a,5,t}]}]
expr1 = BagPart[a, All]
expr2 = BagPart[BagPart[a, 3], All]

Contents[a, 4]
SetContents[a, 7, Contents[a,7]+5]
SetContents[a,11,33]

Даниэль Лихтблау Wolfram Research

Ответ 3

Возможно, это не то, что вы пытаетесь сделать, но Nearest [] может создать NearestFunction [], которая является встроенной структурой quadtree.