Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

At some point I noticed that I was using Riffle and Partition together a lot. I would do things like

Partition[Riffle[{1,2,3},{4,5,6}],2]

or

Partition[Riffle[{1,2,3}, 6 , {2, -1, 2}], 2]

The question is: are better alternatives, in terms of memory and speed?

share|improve this question
6  
{{1, 2, 3}, {4, 5, 6}}\[Transpose]? 3x speedup. Thread[{{1, 2, 3}, 6}] 6x – Kuba Apr 23 '14 at 18:55
    
@Kuba yup ;). See also my chat with Silvia. Anyway I was focussing on the LibraryLink solution for the second "idiom" first. Also see the functions betterish at the bottom of my answer. – Jacob Akkerboom Apr 23 '14 at 18:57
up vote 6 down vote accepted

Here is a solution for both "idioms", using LibraryLink and some kind of custom generics. I'm not sure if I like the abstractions, but here is the answer anyway.

<< SymbolicC`
<< Developer`
<< CCompilerDriver`
<< CCodeGenerator`

The definitions below are made to help abstract this function, so that it works for integers, real numbers and complex numbers.

(*lfn is short for library function name*)
(*LL is short for LibraryLink*)

typeTableHeaders = {"type name", "c type", "LL generic name", 
   "abbreviation in lfn"};
typeTable =
  {
   {Integer, "mint", "Integer", "I"},
   {Real, "mreal", "Real", "R"},
   {Complex, "mcomplex", "Complex", "C"}
   };
cTypes = typeTable[[All, 2]];

abstractFunctionName = "parif";

libraryFunctionNames = 
  abstractFunctionName <> "T" <> # <> "_T" & /@ typeTable[[All, 4]];

getters = StringJoin["MArgument_get", #] & /@ typeTable[[All, 3]];

dataGetters = "MTensor_get" <> # <> "Data" & /@ typeTable[[All, 3]];

tensorTypes = "MType_" <> # & /@ typeTable[[All, 3]];

libraryFunctionNamesSingle = # <> "Single" & /@ libraryFunctionNames;
libraryFunctionNamesMulti = # <> "Multi" & /@ libraryFunctionNames;

relations =
  Hold@
   {
    getter -> getters,
    dataGetter -> dataGetters,
    type -> cTypes ,
    tensorType -> tensorTypes
    };

preRules = Thread /@ First@ MapAt[HoldPattern, relations, {1, All, 1}];

The next definitions distinguish between the "single" case where we want to riffle with a constant, and the "multi" case where we want to riffle two lists. I have not taken the effort of making this code more expressive. That is, there is some duplicate code.

preRulesSingle =
  Append[
   preRules,
   Thread[libraryFunctionName -> libraryFunctionNamesSingle]
   ];

singleDIERule =
  setupSecondInputAccess :> 
   Sequence[
    CDeclare[type, "inputSingleton"], 
    CAssign["inputSingleton", CCall[getter, CArray["Args", 1]]]
    ];

rulesSingle =
  Thread[
   Join[
    preRulesSingle
    ,
    {
     HoldPattern@inputElement -> "inputSingleton",
     HoldPattern@determineNextInputElement :> Sequence[],
     disownIfNeeded :> Sequence[],
     singleDIERule
     }
    ]
   ];

preRulesMulti =
  Append[
   preRules,
   Thread[libraryFunctionName -> libraryFunctionNamesMulti]
   ];


multiDIERule =
  setupSecondInputAccess :> 
   Sequence[
    CDeclare["MTensor", "input2"],
    CAssign["input2", 
     CCall["MArgument_getMTensor", CArray["Args", 1]]],
    CDeclare[CPointerType[type], "input2DataPtr"],
    CAssign[
     "input2DataPtr",
     CCall[CPointerMember["libData", dataGetter], {"input2"}]
     ]
    ];

rulesMulti =
  Thread[
   Join[
    preRulesMulti
    ,
    {
     HoldPattern@inputElement -> CDereference["input2DataPtr"],
     HoldPattern@determineNextInputElement :> 
      COperator[Increment, "input2DataPtr"],
     disownIfNeeded :> 
      CCall[CPointerMember["libData", "MTensor_disown"], "input2"],
     multiDIERule
     }
    ]
   ];

Now we construct some SymbolicC with some tokens inside it

parifSCHeld =
 Hold@
  CFunction[
   "int"
   ,
   libraryFunctionName
   ,
   {{"WolframLibraryData", "libData"}, {"mint", 
     "Argc"}, {CPointerType["MArgument"], "Args"}, {"MArgument", 
     "Res"}}
   ,

   CBlock[
    {
     CDeclare["int", CAssign["err", "LIBRARY_NO_ERROR"]],
     CDeclare["MTensor", "input"],
     CDeclare["MTensor", "result"],
     CDeclare["int", "inputLength"],
     CDeclare["mint", CArray["resultDimensions", 2]],
     CAssign[
      "input", 
      CCall["MArgument_getMTensor", CArray["Args", 0]]
      ],
     CAssign[
      "inputLength",
      CDereference[
       CCall[
        CPointerMember["libData", "MTensor_getDimensions"],
        {"input"}
        ]
       ]
      ],
     CAssign[
      CArray["resultDimensions", 0],
      "inputLength"
      ],
     CAssign[
      CArray["resultDimensions", 1],
      2
      ]
     ,
     CAssign["err", 
      CCall[
       CPointerMember["libData", "MTensor_new"],
       {tensorType, 2, "resultDimensions", CAddress["result"]}
       ]
      ]
     ,
     CDeclare[CPointerType[type], "inputDataPtr"],
     CAssign[
      "inputDataPtr",
      CCall[
       CPointerMember["libData", dataGetter],
       {"input"}
       ]
      ],
     CDeclare[CPointerType[type], "resultDataPtr"],
     CAssign[
      "resultDataPtr",
      CCall[
       CPointerMember["libData", dataGetter],
       {"result"}
       ]
      ]
     ,
     CDeclare["long", CAssign["iter", 0]]
     ,
     setupSecondInputAccess
     ,
     CWhile[
      COperator[Less, {"iter", "inputLength"}],
      CBlock[
       {
        CAssign[
         CDereference["resultDataPtr"],
         CDereference["inputDataPtr"]
         ],

        COperator[Increment, "resultDataPtr"],
        COperator[Increment, "inputDataPtr"],
        CAssign[
         CDereference["resultDataPtr"],
         inputElement
         ],
        determineNextInputElement,
        COperator[Increment, "resultDataPtr"],
        COperator[Increment, "iter"]
        }
       ]
      ]
     ,
     CCall["MArgument_setMTensor", {"Res", "result"}],
     CCall[CPointerMember["libData", "MTensor_disown"], "input"],
     disownIfNeeded,
     CReturn["err"]
     }
    ]
   ];

Now, we turn it into a string, that would normally correspond to the contents of a .c file.

symbCWithReplacementsSingle = parifSCHeld //. rulesSingle;
symbCWithReplacementsMulti = parifSCHeld //. rulesMulti;
cStrings = 
  "DLLEXPORT" <> " " <> ToCCodeString@First@# & /@ 
   Join[symbCWithReplacementsSingle, symbCWithReplacementsMulti];
cFunctionsString = StringRiffle[cStrings, "\n\n"];

boilerPlate = "
  #include \"WolframLibrary.h\"

  /* Return the version of Library Link */
  DLLEXPORT mint WolframLibrary_getVersion( ) {
  \treturn WolframLibraryVersion;
  }

  /* Initialize Library */
  DLLEXPORT int WolframLibrary_initialize( WolframLibraryData \
libData) {
  \treturn LIBRARY_NO_ERROR;
  }

  /* Uninitialize Library */
  DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData \
libData) {
  \treturn;
  }

  ";

totalCString = boilerPlate <> cFunctionsString;

The code below creates the library.

libName = "parif";
CreateLibrary[totalCString, libName]

This loads the functions (again duplicate code)

parifFunctionsSingle = 
  LibraryFunctionLoad[
     libName, #2, {{#, 1, "Shared"}, {#}}, {#, 2}] & @@@ 
   Transpose@{typeTable[[All, 1]], libraryFunctionNamesSingle};

Set @@ {ToExpression["parifSingle" <> #, InputForm, 
      Unevaluated], #2} & @@@ 
  Transpose@{typeTable[[All, 4]], parifFunctionsSingle};

parifFunctionsMulti = 
  LibraryFunctionLoad[
     libName, #2, {{#, 1, "Shared"}, {#, 1, "Shared"}}, {#, 2}] & @@@ 
   Transpose@{typeTable[[All, 1]], libraryFunctionNamesMulti};

Set @@ {ToExpression["parifMulti" <> #, InputForm, 
      Unevaluated], #2} & @@@ 
  Transpose@{typeTable[[All, 4]], parifFunctionsMulti};

Timing comparisons

betterish[a_List,const_]:=Transpose[{a,ConstantArray[const,Length@a]}];
arFlatWiz[a_,b_]:= ArrayFlatten[{{{a}\[Transpose],b}}];

nn=10^7;
randomIntegers=RandomInteger[100,nn];

(resB=betterish[randomIntegers,2])//RepeatedTiming//First
(resAFW=arFlatWiz[randomIntegers, 2])//RepeatedTiming//First
(resSI=parifSingleI[randomIntegers,2])//RepeatedTiming//First

resB === resAFW===resSI

gives

0.201
0.213
0.0769
True

And

randomReals1=RandomReal[1,nn];
randomReals2=RandomReal[1,nn];

(resT=Transpose@{randomReals1, randomReals2})//RepeatedTiming//First
(resMRLL=parifMultiR[randomReals1,randomReals2])//RepeatedTiming//First

resMRLL===resT

gives

0.18
0.0838
True

So indeed we are faster than Transpose here.

share|improve this answer
    
I don't like running the code I don't understand, unless I want to get it :) – Kuba Apr 23 '14 at 19:11
    
@Kuba but don't you know how many little computer gnome calculator slaves are killed every time you run Riffle :P? You run that without knowing how it is implemented, don't you ;)? – Jacob Akkerboom Apr 23 '14 at 19:14
    
This kind of logic will lead us to Assembler, won't it? :) – Kuba Apr 23 '14 at 19:16
    
@Kuba ain't nobody got time for that, I have to hunt for my own food! – Jacob Akkerboom Apr 23 '14 at 19:16
    
Exactly! :P (4togo) – Kuba Apr 23 '14 at 19:18

Kind of stolen from Kuba's earlier comment:

Another way to rewrite the original expression if speed is not a huge concern is via Thread Transpose as in:

Transpose[{{1, 2, 3}, {4, 5, 6}}]

which is especially pleasing to the eye if your initial lists are named variables (Transpose@{a, b}). On my machine this variant seems to be 10 times slower than the original Riffle and Partition.

The advantage and reason for posting (besides being easier to read) is that this construct can be used to merge more than two lists which is especially handy for e.g. ListContourPlot

ListContourPlot@Transpose[{xValues, yValues, fValues}]
share|improve this answer

Since you only give two examples of your use of Riffle and Partition, and since I don't make common use of this combination myself, I can only address those specific examples.

The first can be done nearly an order of magnitude faster on packed arrays with a naive application of Transpose (as already mentioned by Kuba):

a = Range[3*^6];
b = Mod[a, 10];
Partition[Riffle[a, b], 2]  // timeAvg
{a, b}\[Transpose]          // timeAvg
0.03928

0.005616

(Search the site for timeAvg code; I have posted it many times.)

The second can be done several times faster with ArrayFlatten:

a = Range[3*^6]; b = 6;

Partition[Riffle[a, b, {2, -1, 2}], 2] // timeAvg
ArrayFlatten[{{{a}\[Transpose], b}}]   // timeAvg
0.04992

0.01684

For other methods and timings see this closely related Q&A: Prepend 0 to sublists

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.