2011-07-14 76 views
25

Mam zaimplementowane w Mathematica quadtree. Jestem nowy w kodowaniu w funkcjonalnym języku programowania, takim jak Mathematica, i zastanawiałem się, czy mogę to poprawić lub uczynić go bardziej kompaktowym dzięki lepszemu wykorzystaniu wzorców.Implementacja Quadtree w Mathematica

(rozumiem, że mogę być może zoptymalizować drzewa poprzez przycinanie nieużywane węzły, i nie może być lepsze struktury danych, takie jak drzewa kd dla rozkładu przestrzennego.)

Także, ja nadal nie jestem wygodne z ideą kopiowania całe drzewo/wyrażenie za każdym razem, gdy dodawany jest nowy punkt. Ale rozumiem, że operowanie na ekspresji jako całości i niezmodyfikowanie części to funkcjonalny sposób programowania. Byłbym wdzięczny za wszelkie wyjaśnienia dotyczące tego aspektu.

MV

Kodeks

ClearAll[qtMakeNode, qtInsert, insideBox, qtDraw, splitBox, isLeaf, qtbb, qtpt]; 

(* create a quadtree node *) 
qtMakeNode[{{xmin_,ymin_}, {xmax_, ymax_}}] := 
{{}, {}, {}, {}, qtbb[{xmin, ymin}, {xmax, ymax}], {}} 

(* is pt inside box? *) 
insideBox[pt_, bb_] := If[(pt[[1]] <= bb[[2, 1]]) && (pt[[1]] >= bb[[1, 1]]) && 
    (pt[[2]] <= bb[[2, 2]]) && (pt[[2]] >= bb[[1, 2]]), 
    True, False] 

(* split bounding box into 4 children *) 
splitBox[{{xmin_,ymin_}, {xmax_, ymax_}}] := { 
{{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}} 
} 

(* is node a leaf? *) 
isLeaf[qt_] := If[ And @@((# == {})& /@ Join[qt[[1;;4]], {List @@ qt[[6]]}]),True, False] 

(*--- insert methods ---*) 

(* qtInsert #1 - return input if pt is out of bounds *) 
qtInsert[qtree_, pt_] /; !insideBox[pt, List @@ qtree[[5]]]:= qtree 

(* qtInsert #2 - if leaf, just add pt to node *) 
qtInsert[qtree_, pt_] /; isLeaf[qtree] := 
{qtree[[1]],qtree[[2]],qtree[[3]],qtree[[4]],qtree[[5]], qtpt @@ pt} 

(* qtInsert #3 - recursively insert pt *) 
qtInsert[qtree_, pt_] := 
    Module[{cNodes, currPt}, 
    cNodes = qtree[[1;;4]]; 
    (* child nodes not created? *) 
    If[And @@ ((# == {})& /@ cNodes), 
    (* compute child node bounds *) 
    (* create child nodes with above bounds*) 
    cNodes = qtMakeNode[#]& /@ splitBox[List @@ qtree[[5]]]; 
    ]; 
    (* move curr node pt (if not empty) into child *) 
    currPt = List @@ qtree[[6]]; 
    If[currPt != {}, 
    cNodes = qtInsert[#, currPt]& /@ cNodes; 
    ]; 
(* insert new pt into child *) 
cNodes = qtInsert[#, pt]& /@ cNodes; 
(* return new quadtree *) 
{cNodes[[1]],cNodes[[2]], cNodes[[3]], cNodes[[4]], qtree[[5]], {}} 
] 

(* draw quadtree *) 
qtDraw[qt_] := Module[{pts, bboxes}, 
    pts = Cases[qt, _qtpt, Infinity] /. qtpt :> List; 
    bboxes = Cases[qt, _qtbb, Infinity] /. qtbb :> List; 
    Graphics[{ 
    EdgeForm[Black],Hue[0.2], Map[Disk[#, 0.01]&, pts], 
    Hue[0.7],EdgeForm[Red], FaceForm[],(Rectangle @@ #) & /@ bboxes 
    }, 
    Frame->True 
] 
] 

Wykorzystanie

Clear[qt]; 
len = 50; 
pts = RandomReal[{0, 2}, {len, 2}]; 
qt = qtMakeNode[{{0.0, 0.0}, {2.0, 2.0}}]; 
Do[qt = qtInsert[qt, pts[[i]]], {i, 1, len}] 
qtDraw[qt] 

Wyjście

enter image description here

+1

nie mam czasu na zabawę z kodem, ale w odpowiedzi na”... Jestem nadal nie jest wygodnie z ideą kopiowania całego drzewa/wyrażenia za każdym razem, gdy dodawany jest nowy punkt ...'--- czy szukasz sposobu na dodanie pamięci do funkcji (funkcji)? Jeśli tak, szukaj "funkcji, które pamiętają swoje wartości". Jeśli nie tego, czego szukasz, zignoruj, a ja zagram później. Wygląda fajnie. :) – telefunkenvf14

+0

Nie, chodzi mi o to, że kiedy dodaję punkt, zamiast wstawiania węzłów do wychodzącego drzewa, w efekcie generujemy całe nowe drzewo, odrzucając stare. Zastanawiam się tylko nad implikacjami tego zarządzania pamięcią, kiedy takie wyrażenia są naprawdę ogromne. Trudno zejść z myślenia w C++! ;-) Dzięki. –

+1

Jeśli dobrze cię rozumiem, chcesz przejść przez referencję, ponieważ nie lubisz linii 'Return' w' qtInsert', aby utworzyć nową listę zagnieżdżoną z (potencjalnie) ogromnego starego drzewa i kilku nowych węzłów .W Mathematica jest sposób na zrobienie czegoś takiego z 'Attributes [qtInsert] = HoldAll;' z wadą, że wszystkie argumenty 'qtInsert' muszą być zmiennymi, a nie literałami. – bbtrb

Odpowiedz

11

Tutaj jest bardziej kompaktowa wersja. Używa tej samej struktury danych co wersja oryginalna. Funkcje splitBox i są w zasadzie takie same (jak napisane w nieco inny sposób).

Zamiast dodawać punkty jeden po drugim, początkowe pole zawiera wszystkie punkty na początku, więc nie ma potrzeby stosowania procedur qtInsert. W każdym kroku rekursji pola zawierające więcej niż jeden punkt są dzielone, a punkty są rozdzielane na podskładniki. Oznacza to, że wszystkie węzły z więcej niż jednym punktem są liśćmi, więc nie ma również potrzeby sprawdzania tego.

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 -> {}}]]]] 

Przykład (stosując oryginalną wersję qtDraw):

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

Wynik:

Quadtree example

+0

Dzięki za kod - wiele ciekawych zastosowań. Sądzę jednak, że ważna jest umiejętność dodawania punktów jeden po drugim. Weźmy na przykład symulację 2D, w której cząstki pojawiają się jeden po drugim. –

+0

@Heike Albo twój obraz jest pusty, albo moja maszyna ma z tym pewne trudności. –

+0

@belisarius: dziwne, na moim ekranie jest widoczne. Czy widzisz to, gdy podążasz za tym linkiem? http://i.stack.imgur.com/kNyrg.jpg – Heike

3

To nie może być to, co próbujesz zrobić, ale Najbliższe [ ] może utworzyć funkcję NearestFunction [], która jest wbudowana w strukturę quadtree.

35

Myślę, że twój kod nie jest tak głodny pamięci, jak można się spodziewać. To łamie i reformuje listy, ale zwykle sprawia, że ​​większość podlist jest nietknięta.

Jak zauważyli inni, możliwe jest, że lepiej będzie jeszcze używać atrybutów Zawijanie i/lub Atrybuty HoldXXX, aby emulować wywołanie referencyjne.

Dla hard core podejścia do niektórych innych implementacjach struktury danych, patrz

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

Odpowiedni kod znajduje się na notebooka Hemmecke-final.nb (tak nazwany, ponieważ realizuje toryczne Baza Gröbnera algorytm powodu R. Hemmecke i współautorzy).

Zaatakowałem ponownie przy użyciu atrybutów Hold ... ale nie jestem w tym zbyt dobry i poddałem się, gdy kod zabrał mi ukłucie (brakowało, ale zabił moją sesję Mathematica). Zamiast tego mam implementację, która wykorzystuje nieudokumentowany "surowy" typ danych Mathematica, taki, który jest obojętny, a tym samym możliwy do zachowania.

Struktura, o której mowa, nazywa się "workiem expr", ponieważ ogólna struktura danych Mathematica to "expr". Jest jak lista, ale (1) może rosnąć na jednym końcu (choć nie maleć) i (2) jak inne surowe typy ekspresji (np. Wykresy w wersji 8) ma komponenty, do których można uzyskać dostęp i/lub zmienić za pomocą dostarczonych funkcji (API, że tak powiem). Jego "elementy" leżące u podstaw są obojętne w tym sensie, że mogą odwoływać się do KAŻDEGO wyrażenia (w tym do samej torby) i można nimi manipulować w sposób, który wskażę poniżej.

Pierwsza pozycja powyżej przedstawia podstawową technologię wdrażania programu Sow/Reap. Jest to drugie, które będzie interesujące w poniższym kodzie. Na koniec dołączę kilka uwag dotyczących wyjaśnienia struktury danych, ponieważ nie ma na to formalnej dokumentacji.

Zachowałem kod mniej więcej w tym samym stylu, co oryginał, w szczególności pozostaje on-line (tzn. Elementy nie muszą koniecznie wchodzić na początku, ale mogą być dodawane indywidualnie). Zmieniono kilka nazwisk. Wykonane podstawową strukturę zbliżoną do

węzła

(pole, wartość zero lub cztery podwęzłów ograniczająca) jeżeli istnieją podwęzły to pole wartość jest pusta. Pola i pola wartości są reprezentowane przez zwykłe wyrażenie listy Mathematica, chociaż może być sens używanie dedykowanych głowic i bardziej zbliżone do stylu C-struct. Zrobiłem coś takiego w nazewnictwie różnych funkcji dostępu do pól/ustawień.

Jednym z zastrzeżeń jest to, że ten surowy typ danych zużywa znacznie więcej pamięci narzut niż np. lista. Tak więc mój wariant poniżej użyje więcej pamięci niż pierwotnie opublikowany kod. Nie asymptotycznie więcej, tylko przez stały czynnik. Ponadto wymaga on stałego współczynnika narzutowego więcej niż, powiedzmy, porównywalnej struktury C pod względem dostępu lub ustawiania wartości elementu. Więc to nie jest magiczna kula, tylko typ danych z zachowaniem, które nie powinno dawać asymptotycznych niespodzianek.


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]] 

Oto przykład. Zaznaczę, że skalowanie jest rozsądne. Może O (n log (n)) lub tak. Zdecydowanie lepsze niż 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} 

Ogólne wskazówki torba expr. Są stare, więc nie twierdzę, że to wszystko działa tak, jak wskazano.

Te funkcje występują w kontekście wewnętrznym.

Torba Tworzy worek ekspresowy, opcjonalnie z gotowymi elementami.

BagPart Otrzymuje części worka ekspresowego, podobne do części zwykłej . Można go również użyć w Lhs, np. aby zresetować wartość.

Torebka na rzeczy Dołącza elementy do końca torby.

Posiadamy również BagLength. Przydatny do iteracji nad torbą.

Funkcje te są niezwykle użyteczne z dwóch powodów.

Po pierwsze, jest to dobry sposób na tworzenie rozszerzalnego stołu w Mathematica.

Po drugie, zawartość worków jest oceniana, a następnie umieszczana w postaci surowego wyrażenia, w związku z tym są one ekranowane. W ten sposób można z nich korzystać jak „wskaźników” (w sensie C), a nie jak przedmioty, a to nie wymaga zatrzymaniem itd Oto kilka przykładów:

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

Gdybyśmy zamiast używać torby otrzymujemy własny - struktura referencyjna.

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 

Jest to trudne do naśladowania w jakikolwiek inny sposób w Mathematica. Należy użyć rzadkich tabel (skrótów) w jakiś niezauważalny sposób.

Oto pokrewny przykład, nie w pełni debugowany. My w zasadzie wdrożenia połączonej listy przy czym można destrukcyjnie modyfikowania ogony, wymienić list podrzędnych itd

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>]} 

Więc tt jest powiązana lista, pierwszy element jest vv, następny jest sama jest powiązana lista, etc Powstrzymałem się od używania terminologii Lisp (samochód/cdr i tym podobne), ponieważ nie jestem w stanie stwierdzić, czy operacje na liscie List są destrukcyjne. Ale masz ogólny pomysł.

W podobny sposób użyłem worków expr do zaimplementowania drzew binarnych . Jest to przydatne, ponieważ możemy przeprowadzać niszczące zmiany w stałym czasie (zakładając, że mamy już "uchwyt" w punkcie wstawiania/usuwania), a ponadto "surowy" charakter toreb expr oznacza, że ​​całkowicie unikamy nieskończonej oceny Mathematica semantyka.

Być może inna aplikacja.

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]) 

Spróbuj z

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] 

Daniel Lichtblau Wolfram Research

+0

Doskonała informacja! Wydaje się być bardzo szybkim sposobem tworzenia list. – acl

+0

+1, przeczytałem to po raz pierwszy. Dobra informacja, dzięki. – rcollyer