# Implementing a function which generalizes the merging step in merge sort

One of the key steps in merge sort is the merging step. Given two sorted lists

sorted1={2,6,10,13,16,17,19};
sorted2={1,3,4,5,7,8,9,11,12,14,15,18,20};

of integers, we want to produce a new list as follows:

2. Compare the first elements of sorted1 and sorted2. Append the smaller one to acc.
3. Remove the element used in step 2 from either sorted1 or sorted2.
4. If neither sorted1 nor sorted2 is empty, go to step 2. Otherwise append the remaining list to acc and output the value of acc.

Applying this process to sorted1 and sorted2, we get

acc={1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20}

Added in response to Rojo's question: We can carry out this procedure even if the two lists are not pre-sorted. So list1 and list2 below are not assumed to be sorted.

If there were a built-in function MergeList which carries out this process, it would probably take three arguments list1, list2, and f. Here f is a Boolean function of two arguments used to decide which element to pick. In the case of merge sort, f = LessEqual. I feel that MergeList is a fundamental list operation, so

Question 1: Is there such a built-in function or one very close to that?

If I were to write such a function in Scheme, I would use a recursive definition equivalent to the following:

MergeList[list1_,{},f_,acc_:{}]:=Join[acc,list1];
MergeList[{},list2_,f_,acc_:{}]:=Join[acc,list2];
MergeList[list1_,list2_,f_,acc_:{}]:=
If[
f@@First/@{list1,list2},
MergeList[Rest[list1],list2,f,Append[acc,First[list1]]],
MergeList[list1,Rest[list2],f,Append[acc,First[list2]]]
]

Sample output with unsorted lists:

In[2]:= MergeList[{2,5,1},{3,6,4},LessEqual]
Out[2]= {2,3,5,1,6,4}

My impression is that recursive solutions tend to be inefficient in Mathematica, so

Question 2: What would be a better way to implement MergeList?

If you have tips about converting loops into their functional equivalents, feel free to mention them as well.

-
So you want a function that returns the same result as Sort[Join[#1, #2], #3]&? –  Rojo Jun 17 '12 at 7:43
@Rojo For MergeList, the two lists need not be sorted. For example, MergeList[{2, 5, 1}, {3, 4}, LessEqual] evaluates to {2, 3, 4, 5, 1}. –  Michael Wijaya Jun 17 '12 at 7:49
Oh, I see. Then I'll write something up. –  Rojo Jun 17 '12 at 7:56
@Rojo Thanks for pointing out the ambiguity. I edited the question. –  Michael Wijaya Jun 17 '12 at 8:00
If you care for efficiency a bit as you seem to, and you will use it for numeric quantities, this is one of those cases where it's probably very worth it to go with Compile and a procedural solution –  Rojo Jun 17 '12 at 8:05

## Preamble

Since I agree that it would be nice to have a generic function of this type, I will provide a general implementation. First, I will give a generic one based on linked lists, then I will add a JIT-compiled one for special numeric types, and lastly, I will bring it all together.

## Top-level implementation based on linked lists

Here is a reasonably efficient implementation based on linked lists:

SetAttributes[ll, HoldAllComplete];
toLinkedList[s_List] := Fold[ll[#2, #1] &, ll[], Reverse[s]];

and the main function:

ClearAll[merge];
merge[a_ll, ll[], s_, _] := List @@ Flatten[ll[s, a], Infinity, ll];
merge[ll[], b_ll, s_, _] := List @@ Flatten[ll[s, b], Infinity, ll];
merge[ll[a1_, atail_], b : ll[b1_, _], s_, f_: LessEqual] /;f[a1, b1] :=
merge[atail, b, ll[s, a1], f];
merge[a : ll[a1_, _], ll[b1_, brest_], s_, f_: LessEqual] :=
merge[a, brest, ll[s, b1], f];
merge[a_List, b_List, f_: LessEqual] :=

For example:

merge[{2,5,1},{3,6,4},LessEqual]
{2,3,5,1,6,4}
merge[{2,5,1},{3,6,4},Greater]
{3,6,4,2,5,1}

And also for large lists:

large1 = RandomInteger[100, 10000];
large2 = RandomInteger[100, 10000];

Block[{$IterationLimit = Infinity}, merge[large1,large2,LessEqual]]//Short//AbsoluteTiming {0.0751953,{70,54,78,84,11,21,41,49,78,93,90,70,19, <<19975>>,42,2,10,40,53,12,2,47,89,40,2,80}} For a complete implementation of merge sort algorithm based on linked lists, see this post (the difference there is that I used repeated rule application instead of recursion. Originally, the goal of that example was to show that ReplaceRepeated is not necessarily slow if the patterns are constructed efficiently). ## Full implementation including JIT-compilation I'd like to show here how one could implement a fairly complete function which would automatically dispatch to an efficient JIT-compiled code when the arguments are appropriate. Compilation will work not just for numeric lists, but for lists of tensors in general, as long as they are of the same shape. ### JIT - compilation First comes the JIT-compiled version, done along the lines discussed in this answer, section "Making JIT-compiled functions" ClearAll[mergeJIT]; mergeJIT[pred_, listType_, target : ("MVM" | "C") : "MVM"] := mergeJIT[pred, Verbatim[listType], target] = Block[{fst, sec}, With[{decl = {Prepend[listType, fst], Prepend[listType, sec]}}, Compile @@ Hold[decl, Module[{result = Table[0, {Length[fst] + Length[sec]}], i = 0, fctr = 1, sctr = 1}, While[fctr <= Length[fst] && sctr <= Length[sec], If[pred[fst[[fctr]], sec[[sctr]]], result[[++i]] = fst[[fctr++]], (* else *) result[[++i]] = sec[[sctr++]] ] ]; If[fctr > Length[fst], result[[i + 1 ;; -1]] = sec[[sctr ;; -1]], (* else *) result[[i + 1 ;; -1]] = fst[[fctr ;; -1]] ]; result ], CompilationTarget -> target ]]]; You can use this in isolation: mergeJIT[LessEqual,{_Integer,1},"MVM"][{2,5,1},{3,6,4}] {2,3,5,1,6,4} but it is much better to use as a part of the generic function, which would figure out the types for you automatically. ### Generic function implementation This is a function to find the type of our lists: Clear[getType,$useCompile];
getType[arg_List] /; $useCompile && ArrayQ[arg, _, IntegerQ] := {_Integer, Length@Dimensions@arg}; getType[arg_List] /;$useCompile && ArrayQ[arg, _, NumericQ] &&
Re[arg] == arg :=
{_Real, Length@Dimensions@arg};
getType[_] := General;

This is a function to dispatch to a right type:

Clear[mergeDispatch];
SetAttributes[mergeDispatch, Orderless];
mergeDispatch[{Verbatim[_Integer], n_}, {Verbatim[_Real], n_}, pred_] :=
mergeDispatch[{Verbatim[_Real], n}, {Verbatim[_Real], n}, pred];

mergeDispatch[f : {Verbatim[_Real], n_}, {Verbatim[_Real], n_}, pred_] :=
mergeJIT[pred, f, $target]; mergeDispatch[f : {Verbatim[_Integer], n_}, {Verbatim[_Integer], n_}, pred_] := mergeJIT[pred, f,$target];

mergeDispatch[_, _, pred_] :=
Function[{fst, sec},
Block[{$IterationLimit = Infinity}, merge[fst, sec, pred]]]; and this is a function to bring it all together: ClearAll[mergeList]; Options[mergeList] = { CompileToC -> False, Compiled -> True }; mergeList[f_, s_, pred_, opts : OptionsPattern[]] := Block[{$target = If[TrueQ[OptionValue[CompileToC]], "C", "MVM"],
$useCompile = TrueQ[OptionValue[Compiled]] }, mergeDispatch[getType@f, getType@s, pred][f, s] ]; Finally, a helper function to clear the cache of mergeJIT, if that would be desired: ClearAll[clearMergeJITCache]; clearMergeJITCache[] := DownValues[mergeJIT] = {Last@DownValues[mergeJIT]}; ### Benchmarks and tests First, create test data: clearMergeJITCache[]; huge1 = RandomInteger[1000,1000000]; huge2 = RandomInteger[1000,1000000]; A first call to the function with C compilation target is expensive: mergeList[huge1,huge2,Less,CompileToC -> True]//Short//AbsoluteTiming {3.8652344,{267,461,66,607,797,116,197,474,852,805,135, <<1999978>>,266,667,799,280,261,930,241,83,594,904,894}} But then, for the same types of lists, it will pay off for huge lists: mergeList[huge1,huge2,Less,CompileToC -> True]//Short//AbsoluteTiming {0.0468750,{267,461,66,607,797,116,197,474,852,805,135, <<1999978>>,266,667,799,280,261,930,241,83,594,904,894}} On the other hand, the call with MVM target is fast out of the box, but not as fast as the one with the C target after the "warm-up": mergeList[huge1,huge2,Less]//Short//AbsoluteTiming {0.2138672,{267,461,66,607,797,116,197,474,852,805,135, <<1999978>>,266,667,799,280,261,930,241,83,594,904,894}} The call to generic one is general but comparatively very slow: mergeList[huge1,huge2,Less,Compiled->False]//Short//AbsoluteTiming {5.015,{267,461,66,607,797,116,197,474,852,805,135, <<1999978>>,266,667,799,280,261,930,241,83,594,904,894}} - The difference in performance is astounding: 0.12s vs 1.23s for lists of lengths 10000 and 1.13s vs 146.60s for lists of length 100000. I find it reassuring that even without making significant changes to the straightforward recursion, I can get major performance boost by switching to linked lists. Perhaps reading The Little Schemer did not cause too much harm. – Michael Wijaya Jun 17 '12 at 10:22 @MichaelWijaya Check out my new addition :). And yes, linked lists are great and under-appreciated in Mathematica data srtucture. I intend to make a blog post on this topic once we have a blog. – Leonid Shifrin Jun 17 '12 at 10:24 You evaluated mergeList[huge1,huge2,Less,CompileToC -> True]//Short//AbsoluteTiming twice at the beginning of the Benchmarks section, but the timing results are significantly different. Am I missing something here? – Michael Wijaya Jun 17 '12 at 11:09 I think I know what is going on now. There is some caching going on in the previous section, thus the comment on clearMergeJITCache. – Michael Wijaya Jun 17 '12 at 11:10 I was showing that, when you use the JIT-version with the compilation to C, the first call for arguments of a given type takes a while since we have to inslude the time to compile to C. All subsequent calls are instant, because the compiled function is cached. – Leonid Shifrin Jun 17 '12 at 11:11 Here's another approach. mergeLists[lista_, listb_, crit_: LessEqual] := Module[{merge}, merge[list1_, list2_] /; crit[First[list1], First[list2]] := With[{part = TakeWhile[list1, crit[#, First[list2]] &]}, Sow[part]; If[Length[part] == Length[list1], Sow[list2], merge[list1[[Length[part] + 1 ;;]], list2]]]; merge[list2_, list1_] /; crit[First[list1], First[list2]] := merge[list1, list2]; merge[list1_, list2_] := With[ {part = TakeWhile[list1, Not[crit[First[list2], #]] &]}, Sow[part]; If[Length[part] == Length[list1], Sow[list2], merge[list1[[Length[part] + 1 ;;]], list2]]]; Flatten[Reap[merge[lista, listb];][[2]]]] It does give slightly different results from Leonid's code though. For example for list1 = {1, 4, 3}; list2 = {2, 3, 4}; I get with my code mergeLists[{1, 4, 3}, {2, 3, 4}, LessEqual] (* out: {1, 2, 3, 4, 4, 3} *) whereas with Leonid's code I get Block[{$IterationLimit = Infinity}, merge[{1, 4, 3}, {2, 3, 4}, LessEqual]]

(* out: {1, 2, 3, 4, 3, 4} *)

If I take Less instead of LessEqual I get the same result for both codes, so I expect that it has to do with a different treatment of border cases where the two sublists start with the same element.

Taking this issue aside, my code does seem to be faster than Leonid's solution. Consider for example (I'm choosing large1 and large2 such that their intersection is empty to avoid the issue above)

{large1, large2} = Partition[RandomSample[Range[20000]], 10000];

then with Leonid's code I get

Block[{$IterationLimit = Infinity}, merge[large1, large2, LessEqual]] // Short // AbsoluteTiming (* {0.070483,{9941,7246,4261,11184,10148,1867,12324, <<19986>>,6927,17973,10762,9165,19379,11449,7735}} *) and with my code mergeLists[large1, large2, LessEqual] // Short // AbsoluteTiming (* {0.039470,{9941,7246,4261,11184,10148,1867,12324, <<19986>>,6927,17973,10762,9165,19379,11449,7735}} *) - Actually, your code is remarkably fast, I can only beat it when I enable the compilation to C, and it is more compact, since what I do explicitly with type-identification and Dispatch, is done implicitly by LengthWhile here. I was initially suspecting that your version might be slower because typically rule-based varsions not using linked lists are slow due to excessive copying of elements. But in your case, for some reason which I don't yet understand this only becomes noticable for truly huge lists, and is not dominating the running time in any case. Will have to think more about it. – Leonid Shifrin Jun 17 '12 at 11:24 Ok, I got it. Statistically, the purely random data is such that very soon you get the number close to a maximal, so your function does not perform many recursive steps. This will be a killer though: ls1 = Flatten[RandomSample /@ Partition[Range[100000], 5]]; ls2 = Flatten[RandomSample /@ Partition[Range[100000], 5]];Block[{$RecursionLimit = Infinity}, mergeLists[ls1, ls2, Less]]; // AbsoluteTiming (Note that you have to lift the $RecursionLimit in general - which was another thing I initially did not understand - how could you avoid that with your function not being tail-recursive). – Leonid Shifrin Jun 17 '12 at 11:35 Heike, I did not have a motivation to kill your code, really. I just could not understand why it worked so fast on the rather large data I tried. Also, make sure you don't have anything worth saving in the session when trying my test code from the previous comment - it eats up my 8Gb of RAM very quickly and then quits the kernel. – Leonid Shifrin Jun 17 '12 at 11:36 @LeonidShifrin I didn't take it that way. I was a bit surprised by the speed myself so I suspected that there was catch somewhere anyway. – Heike Jun 17 '12 at 12:08 You can easily cure your solution by passing just the starting positions from where to read the lists, and using Do loops instead of LengthWhile to read elements - so that you will never copy lists but always read from the same lista and listb. One other thing though would be to make the function tail-recursive, but I think this is also doable, given that it anyway works by side-effects. – Leonid Shifrin Jun 17 '12 at 12:15 Forget Leonid and Heike's recursive stuff (okay, actually I upvoted both as they are both good responses). But here is a simple, direct version. Note that it will not sort, so if the inputs are unsorted the result will be as well. mergeSortedLists[lista_, listb_, crit_: LessEqual] := Module[ {result, len1 = Length[lista], len2 = Length[listb], i = 1, j = 1}, result = Reap[While[i <= len1 && j <= len2, If[TrueQ[crit[lista[[i]], listb[[j]]]], Sow[lista[[i]]]; i++, Sow[listb[[j]]]; j++]; ]][[2, 1]]; If[i <= len1, result = Join[result, lista[[i ;; -1]]], If[j <= len2, result = Join[result, listb[[j ;; -1]]]]]; result] Here is Heike's example modified slightly: {large1, large2} = Partition[Sort[RandomSample[Range[100000]]], 50000]; In[2236]:= Timing[ml = mergeSortedLists[large1, large2, LessEqual];] ml === Range[100000] Out[2236]= {0.21, Null} Out[2237]= True - +1, a fine solution. I retained recursion because this was closer to the original one in the question, and bacause I find it elegant. But note that on tests ls1 = Flatten[RandomSample /@ Partition[Range[100000], 5]]; ls2 = Flatten[RandomSample /@ Partition[Range[100000], 5]];, my compiled version is about 20 times faster on my machine (mergeList[ls1, ls2, LessEqual]), the one with C target is about 50 times faster (after warm-up), mergeList[ls1, ls2, LessEqual, CompileToC -> True], and even recursive one is about 20 percent faster: mergeList[ls1, ls2, LessEqual, Compiled -> False]. – Leonid Shifrin Jun 18 '12 at 14:33 But if I would go for Reap and Sow, then (for what it's worth) I'd do that just like you did here, rather than use recursion like Heike did, bacause recursion is IMO ill-suited for non-linked lists, since it leads to excessive copying of lists as a part of parameter-passing, and requires extra effort to make a function tail-recursive (both issues can be well seen in Heike-s answer - which actually serves as a good example to illustrate those - as per my comments there). – Leonid Shifrin Jun 18 '12 at 14:43 @Leonid Shifrin Mine is definitely not "top in class" in terms of speed. One could probably do better by filling in a Table. Also using Compile, if/when appropriate, would make it faster. Whether it would then become competitive I do not know. – Daniel Lichtblau Jun 18 '12 at 14:50 I actually was expecting an even much larger overhead of explicit list indexing in the top-level code, than what your code exhibits. Your code is definitely better suited for compilation than the recursive one, which can not be compiled at all (at least my version), because it uses linked lists. – Leonid Shifrin Jun 18 '12 at 15:01 Not to different to Heike's I think, because I haven't followed it line by line. Please let me know if it's too similar to be a separate answer merge[l1_, l2_, f_] := Block[{mergeAux}, mergeAux[list1_, list2_] := mergeAux[list2, Function[fr, Drop[list1, Length@Sow[TakeWhile[list1, f[#, fr] &]]]][ First@list2]]; mergeAux[{}, l_] := Sow[l]; mergeAux[l_, {}] := Sow[l]; Flatten[Reap[mergeAux[l1, l2]][[2, 1]], 1] ] EDIT Same idea but with a custom takeWhile that allows for setting from what position to start counting. Given that, lengthWhile[l_, cond_, from_] := If[# === Null, Length[l] - from + 1, #] &[ Do[If[! cond@l[[c]], Return[c - from, Do]], {c, from, Length[l]}]] takeWhile[l_, cond_, from_] := l[[from ;; from - 1 + lengthWhile[l, cond, from]]] ClearAll[merge]; merge[l1_, l2_, f_] := Block[{mergeAux}, mergeAux[list1_, list2_, in1_, in2_] := mergeAux[list2, list1, in2, in1 + Function[fr, Length@Sow@takeWhile[list1, f[#, fr] &, in1]][list2[[in2]]]]; mergeAux[_, l_, Length[l1] + 1, i_] := Sow[l[[i ;;]] ]; mergeAux[l_, _, i_, Length[l2] + 1] := Sow[l[[i ;;]]]; Flatten[Reap[mergeAux[l1, l2, 1, 1]][[2, 1]], 1]] - It is certainly different from Heike-s answer, because it is tail-recursive, so that one has only to lift$IterationLimit rather than $RecursionLimit. Your solution still suffers from quadratic complexity when you have a large number of flips, but at least it finishes in finite time: ls1 = Flatten[RandomSample /@ Partition[Range[100000], 5]]; ls2 = Flatten[RandomSample /@ Partition[Range[100000], 5]];Block[{$IterationLimit = Infinity}, mergeRojo[ls1, ls2, LessEqual]]; // AbsoluteTiming (takes about 9 secs on my machine), while her code leads to a disaster in such cases. –  Leonid Shifrin Jun 18 '12 at 15:17
In fact, in her case, this is so for a subtle reason: because her function is not tail-recursive, the allocated memory for partial lists she constructs can not be claimed by the garbage-collector, because it is referenced from the function calls (expressions) on expression (call) stack. This leads to a memory exhaustion faster than to a stack blowing, so kernel quits for the lack of memory rather than crashes (don't know which is worse). Anyways, your code is better in this respect, but still sub-optimal performance for this type of data, so, sorry, no upvote this time. But, please do keep it! –  Leonid Shifrin Jun 18 '12 at 15:21
I am still developing my feel for efficient programs, so it helps to look at other implementations even if they are not the fastest possible. Leonid's commentary also adds value to your answer. +1 –  Michael Wijaya Jun 18 '12 at 18:06
@MichaelWijaya then let's keep this going :). I just edited –  Rojo Jun 18 '12 at 18:55
@LeonidShifrin, feel free to check out the edit and comment on it. It's still not the fastest possible, but have I avoided the quadratic complexity? –  Rojo Jun 18 '12 at 18:56