Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm trying to use Compile[] to improve the performance of an algorithm. The algorithm works on a set of matrices which are not of a similar shape. Unfortunately to use the Compile function, lists passed to it must be of similar shape.
To circumvent this inconvenience, I propose to Flatten the set of matrices before passing the resultant flat vector to the compiled function and then reconstructing the vector into the set of matrices before applying the algorithm. Sounds simple but actually implementing it in Mathematica is proving difficult (for me) - help appreciated.

The general gist of what I want to do is as follows:

(* Generate a set of irregular shaped matricies *)
mSet = {Table[RandomReal[], {i, 3}, {j, 5}], Table[RandomReal[], {i, 5}, {j, 4}], Table[RandomReal[], {i, 4}, {j, 6}]};

(* Create a mapping that describes the structure of the set of matrices *)
mStruct = FunctionToGenerateVectorThatDescribesTheStructure[mSet]

(* myCompiledFunc uses mStruct to reassemble the flattened mSet before doing it's thing *)
myCompiledFunc = Compile[{{setOfMats, _Real, 1}, {matStructure, _Integer, 1}},
    (* Reconstruct the setOfMats *)
    ReconstructedSetOfMats = someReconstructionFunction[setOfMats, matStructure];
    (* do stuff with the reconstructed matricies *)
];


(* Use the compiled function *)
(* Flatten the set of matricies so we can pass it and the structure to a compiled function *)
myCompiledFunc[Flatten[mSet], mStruct];    
share|improve this question
2  
Perhaps you can use the answers in Unflattening a list? –  kguler Apr 16 at 13:10
1  
If you attempt to reconstruct the non-tensor object inside Compile it will generate a callback to the main evaluator. –  Simon Woods Apr 16 at 14:10
    
Would you please give an example of the type of operation you wish to perform on the reconstructed matrices? –  Mr.Wizard Apr 16 at 20:58
    
Hi Mr.Wizard. I purposely left out the detail of what was inside the compiled function as it would no doubt generate a whole separate topic of debate - it basically has a number of nested For loops with a number of sequential matrix operations per iteration of the loop. Each subsequent iteration depends on the result of the previous iteration and each iteration requires modification of the shape of the resultant matrix from the calculations. I felt that adding this detail would have been a distraction from my main question. Happy to add it or start another question if you would like. –  CrustyNoodle Apr 16 at 23:13
    
As far as I know you can't dynamically reshape arrays within a single CompiledFunction block without losing (some of?) the advantage of compiling. I am wondering if instead it is possible to create a series of separate compiled functions, perhaps dynamically generated (meta-programming), then pass the data from one to another. I am surely not an expert on compilation so I'm not sure if there is any merit to this idea (e.g. overhead may be significant) but is one possible approach. –  Mr.Wizard Apr 18 at 4:27

3 Answers 3

This is in case of set of 2D arrays. I hope I've not missed the point.

mSet = {RandomReal[1, {3, 5}], RandomReal[1, {5, 4}], RandomReal[1, {4, 6}]};


reco[flatten_, dims_] := Composition[
  MapThread[Partition, {#, dims[[ ;; , 2]]}] &,
  Take[flatten, #] & /@ # &,
  Transpose,
  {Most[Join[{0}, #]] + 1, #} &,
  Accumulate
  ][
    Times @@@ dims]

 reco[Flatten[mSet], Dimensions /@ mSet] == mSet
True
share|improve this answer
    
Nice style with Composition. I'll have to replace the horrible @-chains I currently write sometimes with this idiom. Actually, also makes me think of x // Reverse@Composition[h, g, f] as a nice idiom closer to F#'s forward-pipe operator, which I like. Anyway, +1 –  billisphere Apr 16 at 18:16
    
indeed, Composition is great for code readability :) I'm glad you like it. Tell me if it works at the end, I have not compiled it, I don't compile almost at all so I left here only an approach. :) –  Kuba Apr 16 at 18:20

I'll assume your list of 2D examples, this can easily be extended to arbitrary dimensions, and for that matter to a list with elements of differing depths. On a quick test using

Table[RandomInteger[100, {RandomInteger[{50, 100}], RandomInteger[{50, 100}]}], {2000}];

to generate 2000 randomly sized 2D arrays, over 30X faster than reco:

ranger[list_, lens_] := With[{x = Accumulate@lens},
  Inner[list[[# ;; #2]] &, Most@Prepend[x, 0] + 1, x, List]]

Use example:

test = {RandomReal[10, {2, 3}], RandomReal[20, {5, 5}]};

dims = {{3, 3, 5, 5, 5, 5, 5}, {2, 5}};

target = Flatten[test];

Fold[ranger, target, dims] == test

(* True *)

Note the "dimensions" argument is specified from the "bottom" up, per construction element. This allows, among other things, the targets themselves to be ragged.

Btw- ranger is simply a ragged partitioner I built long ago: given a flat list and a list of lengths, it returns the original list partitioned by the lengths. IIRC, faster than the (undocumented) built-in.

share|improve this answer
1  
Regarding ranger see this Q&A and this comment and link. As I like to say to Leonid in these cases: "Great minds think alike." :-) –  Mr.Wizard Apr 16 at 20:54
    
@Mr.Wizard: Ah, that's neat - favorited. There's actually a way using Extract that can be marginally faster, but uses a weird (bug?) "feature" of extract for pulling spans. Thanks for links! –  rasher Apr 16 at 21:01
    
Please tell me about this Extract behavior. I hope you will make use of my dynamicPartition function; I put some thought into its construction, e.g. partitioning non-List expressions, and I find it quite useful myself. And please let me know if you can improve it! –  Mr.Wizard Apr 16 at 21:04
2  
@Mr.Wizard: Try, e.g., target = Join[Range[100], {{1, 2, 3, Range[20]}}] Rest@Extract[target, {{{}}, {2 ;; 10 ;; 2}, {-1, -1, 1 ;; -1 ;; 2}}] - note the use of Rest and prepended {{}}. This suppresses the error you get normally trying to use spans, and surprisingly allows you to use pretty much and valid part/span spec to get pieces of lists/sublists/etc. returned as a list of results. I've found it faster than say mapping over a list of spans, but don't use it in important code since using "weird" constructs gives me the Willies... –  rasher Apr 16 at 21:13
1  
Interesting; in version 7 I just get Extract::argtu: Extract called with 1 argument; 2 or 3 arguments are expected. >> and the expression unevaluated. I suspect they are extending Extract but for some reason the functionality wasn't finalized. MapAt was silently extended to work with Span as well. You should post a self-Q&A about this behavior in my opinion. I bet it will get a lot of votes. –  Mr.Wizard Apr 16 at 21:17

Here is a compilable function that demonstrates how to use Map to produce a tensor of intermediate results from the set of input matrices. The intermediate results are then processed together to produce a final result:

myCompiledFunc1 =
 Compile[
  {{flat, _Real, 1},
   {dims, _Integer, 2}},
  Block[{first, last, matrix, tensor},

   last = 0;
   tensor = (

       (* Get the next matrix *)
       first = last + 1;
       last = last + #[[1]] #[[2]];
       matrix = Partition[flat[[first ;; last]], #[[2]]];

       (* 
       Get some appropriate intermediate result of fixed dimensions \
from this matrix, e.g. 
       with Dimensions \[Rule] {2} *)
       {Max[Mean /@ matrix],
        Min[Mean /@ Transpose@matrix]}

       ) & /@ dims;

   (* Process the contents of the intermediate tensor *)
   {Min@
     tensor[[All, 1]],
    Max@tensor[[All, 2]]}
   ]
  ]

And here is another compilable function, computing the same thing, that demonstrates how to use Fold to deal with the matrices in pairs, in case there might not be a way to compute a nice tensor of intermediate results to get the final result:

myCompiledFunc2 =
 Compile[
  {{flat, _Real, 1},
   {dims, _Integer, 2}},
  Block[{first, last, prevMatrix, currMatrix},

   last = 0;
   prevMatrix = (
     first = last + 1;
     last = last + dims[[1, 1]] dims[[1, 2]];
     Partition[flat[[first ;; last]], dims[[1, 2]]]
     );

   Block[{minMaxRowMean, maxMinColMean, currMaxRowMean, 
     currMinColMean},
    Fold[
     Function[{acc, dim},

      currMatrix = (
        first = last + 1;
        last = last + dim[[1]] dim[[2]];
        Partition[flat[[first ;; last]], dim[[2]]]
        );

      (* This contrived processing doesn't do anything with \
prevMatrix, but if you wanted to, you could. 
      This is just a "frame" to modify and fill in as you see fit. *)
\
      {minMaxRowMean, maxMinColMean} = acc;
      currMaxRowMean = Max[Mean /@ currMatrix];
      currMinColMean = Min[Mean /@ Transpose@currMatrix];

      prevMatrix = currMatrix;
      {Min[minMaxRowMean, currMaxRowMean],
       Max[maxMinColMean, currMinColMean]}
      ],

     {Max[Mean /@ prevMatrix], Min[Mean /@ Transpose@prevMatrix]},

     Rest@dims
     ]
    ]
   ]
  ]

E.g., an accumulator could be used for any intermediate value you might like (as long as it has a consistent Compile-type throughout evaluation), and the approach could also be modified to deal with triples of matrices rather than just pairs.

The main "trick" in using a compiled function to deal with all your matrices at once is that they'd have to be stored in an array, and arrays in Compile must be tensors (basically not "ragged" in any way). This is essentially the Compile-typing problem you originally encountered that I believe prompted you to go this route, flattening and passing the dimensions of each matrix.

Another "trick" that is not so easy to come by without this site is a list of compilable top-level functions. Fortunately, we have some of this information here:

List of compilable functions

I think the functions I've provided here are good examples of being compilable, but they may not be very efficient. I was mainly going for the former quality, since I don't know the details of the computation you'll be doing.

For convenient reference from this answer:

http://reference.wolfram.com/mathematica/Compile/tutorial/Overview.html

http://reference.wolfram.com/mathematica/Compile/tutorial/Operation.html

The former is where to go for all things Compile, and the latter (a child section of the former) describes the type system and other details of operation for Compile in Mathematica.

Finally, note that you actually could manage to deal with all matrices at once if you were to handle them "virtually", i.e. dynamically indexing into the flattened list based on which matrix and which indexed element you wanted. That probably wouldn't be very efficient, but I haven't tried it.

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.