Skip to main content

list manipulation - Riffle and Partition


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?



Answer



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.


Comments

Post a Comment

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...
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