Skip to main content

list manipulation - 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:



  1. Start with an empty list acc.

  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.




Answer



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:


ClearAll[toLinkedList, ll];
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] :=
    merge[toLinkedList@a, toLinkedList@b, ll[], f];


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


Comments

Post a Comment

Popular posts from this blog

plotting - Filling between two spheres in SphericalPlot3D

Manipulate[ SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, Mesh -> None, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], {n, 0, 1}] I cant' seem to be able to make a filling between two spheres. I've already tried the obvious Filling -> {1 -> {2}} but Mathematica doesn't seem to like that option. Is there any easy way around this or ... Answer There is no built-in filling in SphericalPlot3D . One option is to use ParametricPlot3D to draw the surfaces between the two shells: Manipulate[ Show[SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], ParametricPlot3D[{ r {Sin[t] Cos[1.5 Pi], Sin[t] Sin[1.5 Pi], Cos[t]}, r {Sin[t] Cos[0 Pi], Sin[t] Sin[0 Pi], Cos[t]}}, {r, 1, 2 - n}, {t, 0, Pi}, PlotStyle -> Yellow, Mesh -> {2, 15}]], {n, 0, 1}]
Post a Comment
Read more

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]
Post a Comment
Read more

plotting - Adding a thick curve to a regionplot

Suppose we have the following simple RegionPlot: f[x_] := 1 - x^2 g[x_] := 1 - 0.5 x^2 RegionPlot[{y < f[x], f[x] < y < g[x], y > g[x]}, {x, 0, 2}, {y, 0, 2}] Now I'm trying to change the curve defined by $y=g[x]$ into a thick black curve, while leaving all other boundaries in the plot unchanged. I've tried adding the region $y=g[x]$ and playing with the plotstyle, which didn't work, and I've tried BoundaryStyle, which changed all the boundaries in the plot. Now I'm kinda out of ideas... Any help would be appreciated! Answer With f[x_] := 1 - x^2 g[x_] := 1 - 0.5 x^2 You can use Epilog to add the thick line: RegionPlot[{y < f[x], f[x] < y < g[x], y > g[x]}, {x, 0, 2}, {y, 0, 2}, PlotPoints -> 50, Epilog -> (Plot[g[x], {x, 0, 2}, PlotStyle -> {Black, Thick}][[1]]), PlotStyle -> {Directive[Yellow, Opacity[0.4]], Directive[Pink, Opacity[0.4]],
Post a Comment
Read more