(long) Re: Speed Challenge
- To: mathgroup at smc.vnet.net
- Subject: [mg27202] (long) Re: [mg27159] Speed Challenge
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Fri, 9 Feb 2001 03:10:51 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
"Gareth J. Russell" wrote:
>
> For the efficiency and speed demons among you: a challenge.
>
> Take a list of integers that represents, say, the number balls in each of a
> set of boxes. I want the fastest code that removes ONE ball at random. The
> ball is chosen at random from all the balls (rather than a box being
> randomly chosen). Here is my suggestion, and its speed over 10,000
> iterations on my system:
>
> removeone = Compile[{{l, _Integer, 1}}, Module[{chosenball, chosenbox},
> chosenball = Random[Integer, {1, Apply[Plus, l]}];
> chosenbox =
> Position[
> Map[If[# ? chosenball, 1, 0] &, Drop[FoldList[Plus, 0, l],
> 1]],
> 1][[1, 1]];
> MapAt[# - 1 &, l, chosenbox]
> ]
> ];
>
> Do[removeone[{5, 4, 7, 6, 5, 4, 7, 8, 9, 4, 5, 6}], {10000}]; // Timing
>
> {1.06667 Second, Null}
>
> I think there must be faster way to do the middle expression (chosenbox =
> ...). Can you beat it!
>
> Gareth
>
> P.S. I do have a practical use for this...
>
> ==================================================
> Dr. Gareth J. Russell
> University of Tennessee
> Columbia University (as visitor)
>
> MAILING ADDRESS (new as of 15TH JANUARY, 2001)
>
> 110 Geosciences
> Lamont-Doherty Earth Observatory
> Columbia University
> 61 Route 9W
> Palisades, NY 10964, U.S.A.
>
> PHONE: ++1 845 365 8683
> FAX: ++1 845 365 8156 (not conveniently located, avoid if possible)
> E-MAIL: russell at cerc.columbia.edu (NO CHANGE)
> WWW: http://web.utk.edu/~grussell (NO CHANGE)
> ==================================================
I fixed some typos. Also note that I use "box" and "bin" interchangeably
below.
removeone = Compile[{{l,_Integer,1}}, Module[
{chosenball,chosenbox},
chosenball = Random[Integer, {1,Apply[Plus,l]}];
chosenbox = Position[
Map[If[# >= chosenball, 1, 0] &,
Drop[FoldList[Plus,0,l],1]],1][[1,1]];
MapAt[# - 1 &, l, chosenbox]
]];
mylist = {5, 4, 7, 6, 5, 4, 7, 8, 9, 4, 5, 6};
In[82]:= Timing[Do[removeone[mylist], {10000}]]
Out[82]= {1.58 Second, Null}
You can save some of time (factor > 2) by not summing the list each call
and by not doing that FoldList and Position. For this you would prepend
to your list the sum of the box sizes, and decrement it in the result
along with the size of the relevant bin.
removeonetoo = Compile[{{ll,_Integer,1}},
Module[{chosenball,chosenbox=2,sum,res=ll},
chosenball = Random[Integer, {1,ll[[1]]}];
sum = ll[[2]];
While[chosenball>sum,
chosenbox++;
sum += ll[[chosenbox]];
];
res[[1]]--;
res[[chosenbox]]--;
res
]
];
In[84]:= InputForm[newlist = Prepend[mylist, Apply[Plus,mylist]]]
Out[84]//InputForm= {70, 5, 4, 7, 6, 5, 4, 7, 8, 9, 4, 5, 6}
In[85]:= removeonetoo[newlist]
Out[85]= {69, 5, 4, 7, 6, 5, 4, 7, 8, 9, 3, 5, 6}
In[86]:= Timing[Do[removeonetoo[newlist], {10000}]]
Out[86]= {0.59 Second, Null}
There are two other ways to do this sort of thing, and which of the
possibilities you will want to use depends on some input sizes. I'll
show the simpler one first. It is essentially the fast shuffle I posted
in a thread several weeks ago.
To start, you "inflate" your bins by repeating the bin number once for
every ball it contains. Prepend to this a "pointer" telling us where to
start. When called, swap the ball in the start position with a ball
randomly chosen between there and the end, and increment the start
value. We do this one a bit differently to avoid recreating at every
call a large list. Instead we'll alter in place by making the list a
held argument. Unfortunately we now lose the ability to use Compile.
removeonetwostep[ll_] := Module[
{start=ll[[1]]+1,ball,tmp},
(*ll[[1]]++;*)
ball = Random[Integer,{start,Length[ll]}];
ll[[{start,ball}]] = ll[[{ball,start}]];
ll
]
In[113]:= InputForm[newlist2 =
Prepend[Flatten[MapIndexed[Table[#2,{#1}]&, mylist]], 1]]
Out[113]//InputForm=
{1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5,
5,
5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9,
9,
9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 11, 11, 11, 11, 11, 12, 12, 12, 12,
12, 12}
In[114]:= InputForm[removeonetwostep[newlist2]]
Out[114]//InputForm=
{1, 9, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5,
5,
5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9,
9,
9, 9, 9, 9, 1, 9, 10, 10, 10, 10, 11, 11, 11, 11, 11, 12, 12, 12, 12,
12, 12}
For real runs we would uncomment the increment line, and the first
position above would be '2' rather than '1'.
If I use this on newlist2 10000 times it is a bit slow. If instead we do
mylistbig = Table[Random[Integer,{1,100}], {500}];
newlistbig = Prepend[mylistbig, Apply[Plus,mylistbig]];
newlistbig2 =
Prepend[Flatten[MapIndexed[Table[#2,{#1}]&, mylistbig]], 1];
then this is a bit faster than removeonetoo. If you have alot of bins
but few balls per bin, this may be much faster. In that case you might
want to use Compile and abandon the call-by-reference (the HoldAll
attribute, that is), instead explicitly creating a return value
initialized to the input and altering that instead of the input list.
Here are some trouble spots. If you have a large total number of balls
(say 10^9), then removeonetwostep is out of the question because of the
inflation step. Moreover removeone and removeonetoo have a problem when
the number of bins gets large. Each time called they must do a linear
search to find the bin from which they removed a ball. One way to
improve on this is to use a tree structure. As I want to maintain use of
Compile this gets a bit messy. So I'll walk through the data structure
and code.
Our "tree" will be implemented as a rank-2 tensor of nodes. Each will be
a 4-tuple of the form
{left_pointer,right_pointer,number_to_left,number_in_node}
The "pointers" are simply position numbers in the list of tuples. The
"number_to_left" refers only to the total in the left-side subtree.
We augment with a 4-tuple at the end whose elements are
{total_in_set, center_node, chosen_ball, 0}
We first assume we have this structure (I'll discuss its construction
later). Then we get (gotta love the name) removeonetootree as follows.
(i) Pick a ball between 1 and total.
(ii) Decrement total and store the chosen ball in the last tuple in the
tree.
(iii) Start at the center node. If the number to the left is less than
our ball, and the number to left + number in node is greater than or
equal to our ball, then it came from the center node. Else if number to
left is larger, walk to the left subtree. Else walk to the right
subtree.
(iii-a) To walk to left, we must modify "number to the left" It will now
mean "total number in bins less than current bin" rather than just
number in the left subtree. To get this we subtract the number we had to
the left, go left, and add the number to the left of that subtree.
(iii-b) To walk right, add the number in the node we are leaving (as it
resides in a bin less than the one we are entering), and also add the
number in the left subtree of the bin we enter.
(iv) When we arrive at the appropriate bin, remove one element from its
ball count. Now return the altered tree.
Here is the code.
removeonetootree = Compile[{{ll,_Integer,2}}, Module[
{ball,bin,sum,res=ll,node,posn,numtoleft,tlen=Length[ll]},
ball = Random[Integer, {1,res[[tlen,1]]}];
res[[tlen,1]]--;
res[[tlen,3]] = ball;
posn = res[[tlen,2]];
node = res[[posn]];
numtoleft = node[[3]];
While[!(numtoleft<ball && numtoleft+node[[4]]>=ball),
If [numtoleft>=ball, (* walk to left *)
res[[posn,3]]--;
numtoleft -= node[[3]];
posn = node[[1]];
node = res[[posn]];
numtoleft += node[[3]],
(* else walk to right *)
numtoleft += node[[4]];
posn = node[[2]];
node = res[[posn]];
numtoleft += node[[3]];
];
];
res[[posn,4]]--;
res
]]
We'll check that it really "compiled" (first, because my early efforts
did not do so successfully, and second because reading the input line
has a nice ring).
In[173]:= removeonetootree[[4]] // InputForm
Out[173]//InputForm=
{{1, 5}, {12, 0, 1}, {57, 0, 0}, {4, 1, 1}, {64, 1, 0, 0, 0, 1, 0, 2},
{14, 0, 2, 0}, {4, 1, 1}, {14, 0, 1, 1}, {23, 2, 1, 0, 3}, {4, 1, 2},
{64, 1, 0, 0, 0, 2, 0, 1}, {9, 1, 2}, {4, 1, 4}, {4, 1, 5},
{64, 1, 0, 0, 0, 5, 0, 6}, {4, 1, 5}, {32, 5, 7}, {24, 6, 7, 6},
{65, 1, 0, 0, 0, 4, 0, 6}, {4, 3, 1}, {65, 1, 0, 0, 0, 1, 0, 3}, {4, 2,
1},
{64, 1, 0, 0, 0, 1, 0, 6}, {64, 1, 0, 6, 1, 2}, {4, 3, 1},
{64, 2, 0, 1, 0, 4}, {48, 4, 3, 0}, {4, 4, 1}, {64, 2, 0, 1, 0, 7},
{24, 4, 7, 1}, {50, 3, 1, 1}, {52, 0, 1, 2}, {55, 2, 0}, {41, 0, 42},
{50, 3, 4, 2}, {41, 2, 26}, {4, 3, 1}, {64, 1, 0, 6, 0, 1, 0, 7}, {9,
7, 1},
{4, 3, 5}, {4, 3, 8}, {64, 1, 0, 6, 0, 8, 0, 9}, {4, 1, 8}, {32, 8,
10},
{24, 9, 10, 9}, {65, 1, 0, 6, 0, 5, 0, 9}, {4, 3, 7}, {64, 2, 0, 7, 0,
9},
{32, 9, 7}, {24, 4, 7, 9}, {9, 9, 4}, {4, 1, 9}, {64, 2, 0, 9, 0, 7},
{9, 7, 6}, {64, 1, 0, 6, 1, 3}, {12, 3, 2}, {4, 3, 7}, {64, 2, 0, 7, 0,
9},
{24, 4, 9, 7}, {9, 7, 4}, {42, 14}, {4, 4, 1}, {64, 2, 0, 1, 0, 7},
{24, 4, 7, 1}, {9, 1, 4}, {4, 2, 1}, {64, 2, 0, 1, 0, 7}, {9, 7, 6},
{64, 1, 0, 6, 1, 3}, {12, 3, 2}, {4, 3, 7}, {64, 2, 0, 7, 0, 1},
{24, 4, 1, 7}, {9, 7, 4}, {42, -48}, {4, 4, 7}, {64, 1, 0, 6, 0, 7, 0,
1},
{9, 1, 7}, {4, 4, 9}, {4, 4, 5}, {64, 1, 0, 6, 0, 5, 0, 10}, {4, 1, 5},
{32, 5, 8}, {24, 10, 8, 10}, {65, 1, 0, 6, 0, 9, 0, 10}, {2}}
All tuples of numbers, so it compiled clean.
So how do we make the trees? I did this in steps. First is a procedure
that will return a (non-tensor) triple of the form
{total, position_of_center_4_tuple, {4_tuples}}
This is convenient because it can be called recursively on left and
right sublists of the input, with those put together along with a center
node. To avoid degenerate cases we handle lists of length 1 and 2
separately.
toTree[ll_List] /; Length[ll]==1 := {ll[[1,2]], ll[[1,1]],
{{0,0,0,ll[[1,2]]}}}
toTree[ll_List] /; Length[ll]==2 := {ll[[1,2]]+ll[[2,2]], ll[[2,1]],
{{0,0,0,ll[[1,2]]}, {ll[[1,1]],0,ll[[1,2]],ll[[2,2]]}}}
toTree[ll_List] := Module[
{mid=Ceiling[(Length[ll]+1)/2], left, right, node},
left = toTree[Take[ll,mid-1]];
right = toTree[Drop[ll,mid]];
node = {left[[2]], right[[2]], left[[1]], ll[[mid,2]]};
{left[[1]]+right[[1]]+ll[[mid,2]], ll[[mid,1]],
Join[left[[3]], {node}, right[[3]]]}
]
To put this into the tensor we need we augment with the last 4-tuple as
described above.
treetensor[ll_] :=
With[{tree=toTree[Transpose[{Range[Length[ll]],ll}]]},
Developer`ToPackedArray[Join[tree[[3]], {{tree[[1]],tree[[2]],0,0}}]]]
In[181]:= InputForm[forest = treetensor[mylist]]
Out[181]//InputForm=
{{0, 0, 0, 5}, {1, 3, 5, 4}, {0, 0, 0, 7}, {2, 6, 16, 6}, {0, 0, 0, 5},
{5, 0, 5, 4}, {4, 10, 31, 7}, {0, 0, 0, 8}, {8, 0, 8, 9}, {9, 12, 17,
4},
{0, 0, 0, 5}, {11, 0, 5, 6}, {70, 7, 0, 0}}
In[182]:= removeonetootree[forest]//InputForm
Out[182]//InputForm=
{{0, 0, 0, 5}, {1, 3, 5, 4}, {0, 0, 0, 7}, {2, 6, 16, 6}, {0, 0, 0, 5},
{5, 0, 5, 4}, {4, 10, 31, 7}, {0, 0, 0, 8}, {8, 0, 8, 9}, {9, 12, 17,
4},
{0, 0, 0, 5}, {11, 0, 5, 5}, {69, 7, 66, 0}}
In[183]:= Timing[Do[removeonetootree[forest], {10000}];]
Out[183]= {0.98 Second, Null}
We see that this is competitive with the other methods despite the
relatively intricate tree manipulation. Now let's try some examples with
10^5 bins.
mylistverybig = Table[Random[Integer,{1,1000}], {10^5}];
newlistverybig = Prepend[mylistverybig, Apply[Plus,mylistverybig]];
amazonrainforest = treetensor[mylistverybig];
Note that this offline computation of the tree tensor is itself slow,
and of course is a full iteration over the list. So this strategy is not
useful unless there will be alot of iterations. Also the tree building
code relies on Mathematica recursion which can be slow. So for alot of
bins it would be best to rewrite that using iterative code.
In[219]:= Timing[Do[removeonetoo[newlistverybig], {100}];]
Out[219]= {12.45 Second, Null}
In[221]:= Timing[Do[removeonetootree[amazonrainforest], {100}];]
Out[221]= {3.05 Second, Null}
----------------------------------
Getting the absolute best speed can be a bit tricky and I'd not be
surprised if there are yet faster ways to do this sort of thing; I am
hoping I did not miss one that is entirely obvious. That said, I do
claim to have given efficient programming in Mathematica a fair amount
of thought over the past few years. Some useful ideas related to those
discussed in this post may be found at the URLs below. The first covers
both deck shuffling and a similar tree example. The second deals more
generally with effective data structure use in Mathematica.
http://library.wolfram.com/conferences/conference98/abstracts/symbolic_FAQ.html
http://library.wolfram.com/conferences/devconf99/
For this second one, scroll to "Programming" section at the bottom of
the page, and you can click on the Mathematica notebook option under
heading
"Enhance Your Code Performance: Effective Data Structures Techniques in
Mathematica"
An HTML version may be found in the same location, or more directly at:
http://library.wolfram.com/conferences/devconf99/lichtblau/
----------------------------------
Daniel Lichtblau
Wolfram Research
- Follow-Ups:
- Re: (long) Re: Speed Challenge
- From: Daniel Lichtblau <danl@wolfram.com>
- Re: (long) Re: Speed Challenge