32

Mathematicaで四分木を実装しました。Mathematica のような関数型プログラミング言語でコーディングするのは初めてで、これを改善したり、パターンをより適切に使用してよりコンパクトにしたりできないかと考えていました。

(未使用のノードを剪定することでツリーを最適化できる可能性があることを理解しています。また、kd ツリーのような空間分解のためのより良いデータ構造があるかもしれません。)

また、新しいポイントが追加されるたびにツリー/式全体をコピーするという考えにはまだ満足していません。しかし、私の理解では、部分を変更せずに式全体を操作することが関数型プログラミングの方法です。この点について説明をいただければ幸いです。

MV

コード

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

使用法

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]

出力

ここに画像の説明を入力

4

3 に答える 3

43

あなたのコードは、あなたが期待するほどメモリを消費していないと思います。リストを壊したり再構築したりしますが、ほとんどのサブリストをそのまま維持する傾向があります。

他の人が指摘したように、call-by-reference をエミュレートするために、Hold ラッパーや HoldXXX 属性を使用する方が良いかもしれません。

関連するデータ構造の実装に対するハードコアなアプローチについては、以下を参照してください。

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

関連するコードは、ノートブック Hemmecke-final.nb にあります (R. Hemmecke と共著者によるトーリック Groebner 基底アルゴリズムを実装しているため、この名前が付けられています)。

Hold... 属性を使用して再実装することに挑戦しましたが、それがあまり得意ではなく、コードが突き返したときにあきらめました (見逃したが、Mathematica セッションを終了させました)。その代わりに、文書化されていない「生の」Mathematica データ型を使用する実装を用意しました。このデータ型は不活性であり、したがって参照渡しの動作に適しています。

一般的な Mathematica データ構造が「expr」であるため、問題の構造は「expr bag」と呼ばれます。これは List に似ていますが、(1) 一端で拡大できます (縮小はしません)。(2) 他の生の式タイプ (バージョン 8 のグラフなど) と同様に、提供された関数を介してアクセスおよび/または変更できるコンポーネントがあります。 (いわばAPI)。その基礎となる「要素」は、任意の expr (バッグ自体を含む) を参照でき、以下に示す方法で操作できるという意味で不活性です。

上記の最初の項目は、Sow/Reap を実装するための基盤となるテクノロジを提供します。以下のコードで注目されるのは 2 番目のコードです。最後に、データ構造を説明する行に沿っていくつかのコメントを含めます。これについては正式なドキュメントがないためです。

私はコードを元のコードとほぼ同じスタイルに保ちました。特に、オンライン バージョンのままです (つまり、最初からすべての要素を入れる必要はなく、個別に追加することができます)。いくつかの名前を変更しました。基本構造を似たものにしました

ノード (バウンディング ボックス、値、ゼロまたは 4 つのサブノード)

サブノードがある場合、値フィールドは空です。ボックスと値フィールドは、通常の Mathematica List 式で表されますが、専用のヘッドを使用して C 構造体スタイルに似たものにするのが理にかなっているかもしれません。さまざまなフィールドのアクセス/設定関数に名前を付ける際に、そのようなことをしました。

1 つの注意点は、この raw データ型は、List などよりもかなり多くのメモリ オーバーヘッドを消費することです。したがって、以下の私のバリアントは、最初に投稿されたコードよりも多くのメモリを使用します。漸近的にそれ以上ではなく、一定の要因によるものです。また、要素値へのアクセスまたは設定に関して、たとえば、同等の 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 バッグのメモ。これらは古いので、これがすべてまだ示されているように機能するとは主張しません。

これらの関数は Internal` コンテキストに存在します。

Bag 必要に応じて事前設定された要素を使用して expr バッグを作成します。

BagPart 通常の expr の Part と同様に、expr バッグのパーツを取得します。値をリセットするなど、lhs で使用することもできます。

StuffBag 要素をバッグの末尾に追加します。

BagLength もあります。バッグを反復処理するのに便利です。

これらの関数は、2 つの理由で非常に便利です。

まず、これは Mathematica で拡張可能なテーブルを作成する良い方法です。

次に、バッグの内容が評価されますが、生の expr に配置されるため、シールドされます。したがって、これらをオブジェクトとしてではなく「ポインタ」(C の意味で) として使用でき、これには 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 のリスト操作が破壊的であるかどうかを思い出すことができないため、Lisp 用語 (car/cdr など) の使用を控えました。しかし、あなたは一般的な考えを得る。

同様に、バイナリ ツリーを実装するために expr バッグを使用しました。これは、一定の時間内に破壊的な変更を行うことができるため (挿入/削除のポイントに「ハンドル」が既にあると仮定して) 便利です。さらに、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]

Daniel Lichtblau Wolfram Research

于 2011-07-22T20:28:45.240 に答える
13

これは、よりコンパクトなバージョンです。元のバージョンと同じデータ構造を使用します。関数splitBoxinsideBoxも本質的に同じです (わずかに異なる方法で記述されているだけです)。

ポイントを 1 つずつ追加する代わりに、最初のボックスには最初にすべてのポイントが含まれているため、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]

結果:

四分木の例

于 2011-07-14T17:40:25.527 に答える
3

これはあなたがしようとしていることではないかもしれませんが、Nearest[] は組み込みの四分木構造である NearestFunction[] を作成できます。

于 2011-07-15T00:03:12.327 に答える