wolfram mathematica - prueba - Implementando un Quadtree en Mathematica

wolfram mathematica free online

Question

He implementado un quadtree en mathematica. Soy nuevo en la codificación de un lenguaje de programación funcional como Mathematica, y me preguntaba si podría mejorar esto o hacerlo más compacto mediante un mejor uso de los patrones.

(Entiendo que tal vez pueda optimizar el árbol podando los nodos no utilizados, y podría haber mejores estructuras de datos como los árboles kd para la descomposición espacial).

Además, todavía no me siento cómodo con la idea de copiar todo el árbol / expresión cada vez que se agrega un nuevo punto. Pero mi entendimiento es que operar en la expresión como un todo y no modificar las partes es la forma de programación funcional. Apreciaría cualquier aclaración sobre este aspecto.

MV

El código

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

Uso

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]

Salida

Answer 1

Aquí hay una versión más compacta. Utiliza la misma estructura de datos que la versión original. Las funciones splitBox e insideBox son esencialmente las mismas también (escritas de una manera ligeramente diferente).

En lugar de agregar puntos uno por uno, el cuadro inicial contiene todos los puntos al principio, por lo que no es necesario qtInsert rutinas qtInsert . En cada paso de recursión, los cuadros que contienen más de un punto se dividen y los puntos se distribuyen en los cuadros secundarios. Esto significa que todos los nodos con más de un punto son hojas, por lo que tampoco es necesario verificarlo.

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

Ejemplo (usando la versión original de qtDraw ):

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

Resultado:

Answer 2

Creo que su código no tiene tanta memoria como se podría esperar. Rompe y reforma las listas, pero tiende a mantener intactos a la mayoría de los sublistas.

Como otros comentaron, podría ser mejor seguir usando envoltorios de retención y / o atributos de HoldXXX, para emular la llamada por referencia.

Para un enfoque de núcleo duro para algunas implementaciones de estructuras de datos relacionadas, consulte

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

El código relevante está en el cuaderno Hemmecke-final.nb (llamado así porque implementa el algoritmo de base de Groebner tórico debido a R. Hemmecke y coautores).

Intenté volver a implementar el uso de los atributos de Hold ..., pero no soy muy bueno en eso y lo dejé cuando el código me devolvió una puñalada (fallida, pero maté a mi sesión de Mathematica). Entonces, en lugar de eso, tengo una implementación que utiliza un tipo de datos de Mathematica "sin procesar" no documentados, uno que es inerte y, por lo tanto, susceptible de comportamiento de llamada por referencia.

La estructura en cuestión se denomina "bolsa expr" porque la estructura de datos genérica de Mathematica es la "expr". Es como una lista pero (1) puede crecer en un extremo (aunque no encogerse) y (2) al igual que otros tipos de expresiones en bruto (por ejemplo, gráficos en la versión 8) tiene componentes a los que se puede acceder y / o modificar a través de las funciones proporcionadas (Una API, por así decirlo). Sus "elementos" subyacentes son inertes en el sentido de que pueden hacer referencia a CUALQUIER expr (incluida la bolsa en sí) y pueden manipularse de la manera que se indica a continuación.

El primer artículo anterior proporciona la tecnología subyacente para la implementación de Sow / Reap. Es el segundo que será de interés en el siguiente código. Al final, incluiré algunas observaciones junto con la explicación de la estructura de los datos, ya que no hay documentación formal para esto.

Mantuve el código más o menos en el mismo estilo que el original, y en particular sigue siendo una versión en línea (es decir, no es necesario que todos los elementos entren al principio, sino que se pueden agregar individualmente). Cambió algunos nombres. Hecho la estructura básica similar a

nodo (cuadro delimitador, valor, cero o cuatro subnodos)

Si hay subnodos, el campo de valor está vacío. Los campos de cuadro y valor están representados por la expresión habitual de la Lista de Mathematica, aunque podría tener sentido usar cabezas dedicadas y tener un aspecto más parecido a un estilo C struct. Hice algo así al nombrar las diversas funciones de acceso / configuración de campos.

Una advertencia es que este tipo de datos sin procesar consume sustancialmente más sobrecarga de memoria que, por ejemplo, una lista. Así que mi variante a continuación utilizará más memoria que el código publicado originalmente. No asintóticamente más, solo por un factor constante. También requiere un factor constante en la sobrecarga más que, digamos, una estructura C comparable en términos de acceso o configuración del valor del elemento. Entonces, no es una bala mágica, solo un tipo de datos con comportamiento que no debería dar sorpresas asintóticas.

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

Aquí hay un ejemplo. Observaré que la escala es razonable. Tal vez O (n log (n)) o algo así. Definitivamente mejor que 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}

Notas de la bolsa expr general. Estos son viejos, así que no pretendo que todo esto todavía funcione como se indica.

Estas funciones viven en contexto interno.

Bolsa Crea una bolsa expr, opcionalmente con elementos predefinidos.

BagPart Obtiene partes de una bolsa expr, similar a Parte para exprs ordinarios. También se puede utilizar en un lhs, por ejemplo, para restablecer un valor.

StuffBag Anexa elementos al final de una bolsa.

También tenemos un BagLength. Útil para iterar sobre una bolsa.

Estas funciones son extremadamente útiles por dos razones.

Primero, esta es una buena manera de hacer una tabla extensible en Mathematica.

En segundo lugar, el contenido de las bolsas se evalúa, pero luego se coloca en un expr en bruto, por lo tanto, se blindan. Por lo tanto, uno puede usar estos como "punteros" (en el sentido C) en lugar de como objetos, y esto no requiere retención, etc. Aquí hay algunos ejemplos:

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

Si en cambio usamos bolsas obtenemos una estructura auto-referencial.

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

Esto es difícil de emular de cualquier otra manera en Mathematica. Uno necesitaría usar tablas dispersas (hashing) de alguna manera no muy transparente.

Aquí hay un ejemplo relacionado, no completamente depurado. Básicamente implementamos una lista enlazada por la cual uno puede modificar colas destructivamente, reemplazar sublistas, etc.

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

Así que tt es una lista enlazada, el primer elemento es vv, el siguiente es una lista enlazada, etc. Me abstuve de usar la terminología Lisp (car / cdr y similares) porque no puedo recordar si las operaciones de la lista de Lisp son destructivas. Pero tienes la idea general.

En líneas similares, he usado bolsas expr para implementar árboles binarios. Esto es útil porque podemos hacer cambios destructivos en tiempo constante (asumiendo que ya tenemos un "control" en el punto de inserción / eliminación), y además la naturaleza "cruda" de expr bags significa que evitamos completamente la semántica de evaluación infinita de Mathematica.

Otra aplicación, tal vez.

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

Prueba con

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

Answer 3

Puede que esto no sea lo que está intentando hacer, pero Nearest [] puede crear un NearestFunction [] que es una estructura integrada de quadtree.