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 am attempting to Compile this function which is used to calculate a elastic collision in two dimensions. The two first coordinates of each list are the position and the last are the velocity of each colliding particle. It should take an argument like:

list={{1,2,3,4},{5,6,7,8}}

cF = Compile[{{list, _Real, 2}}, 
  If[
    (list[[2]] - list[[1]])[[{1, 2}]] != {0, 0}, 
    Map[
     Flatten[
      {#[[1,1 ;; 2]] 
       + (#[[1,3 ;; 4]].{(#[[1, 2]] - #[[2, 2]]), (#[[2, 1]] - #[[1,1]])}
          / Sqrt[(#[[1, 2]] - #[[2, 2]])^2 + (#[[2,1]] - #[[1, 1]])^2]
         ) {(#[[1, 2]] - #[[2,2]]), (#[[2, 1]] - #[[1, 1]])}
          /Sqrt[(#[[1, 2]] - #[[2, 2]])^2 + (#[[2, 1]] - #[[1,1]])^2] 
       + (#[[2,3 ;; 4]].(#[[2]] - #[[1]])[[{1, 2}]]
          / Sqrt[(#[[2]] - #[[1]])[[1]]^2 + (#[[2]] - #[[1]])[[2]]^2]
         ) (#[[2]] - #[[1]])[[{1, 2}]]
          /Sqrt[(#[[2]] - #[[1]])[[1]]^2 + (#[[2]] - #[[1]])[[2]]^2] 
     , (#[[1, 3 ;; 4]].{(#[[1, 2]] - #[[2, 2]]), (#[[2, 1]] - #[[1, 1]])}
         / Sqrt[(#[[1, 2]] - #[[2, 2]])^2 + (#[[2, 1]] - #[[1,1]])^2]
       ) {(#[[1, 2]] - #[[2, 2]]), (#[[2, 1]] - #[[1, 1]])}
         / Sqrt[(#[[1, 2]] - #[[2, 2]])^2 + (#[[2, 1]] - #[[1,1]])^2] 
       + (#[[2, 3 ;; 4]].(#[[2]] - #[[1]])[[{1, 2}]]
          * Sqrt[(#[[2]] - #[[1]])[[1]]^2 + (#[[2]] - #[[1]])[[2]]^2]
         ) (#[[2]] - #[[1]])[[{1, 2}]]    
          / Sqrt[(#[[2]] - #[[1]])[[1]]^2 + (#[[2]] - #[[1]])[[2]]^2]}
     ] &, (* End Flatten *)

     {{list[[1]], list[[2]]}, {list[[2]], list[[1]]}}
    ], (* End Map; End True branch of If *)
    Map[Flatten[{#[[{1, 2}]] + #[[{3, 4}]], #[[{3, 4}]]}] &, list]
  ] (* End If *)
]; (* End Compile*)

However, it throws an error:

Compile::cplist:"Compile`GetElement[System`Private`CompileSymbol[0],System`Private`CompileSymbol[1]][[1,Compile`$5]] should be a tensor of type Integer, Real, or Complex; evaluation will use the uncompiled function."
share|improve this question
3  
Define "doesn't work," as there could be any number of things happening here. Does the code work prior to trying to compile it? Also, as written, it is insanely complex and practically impossible to determine what it is supposed to do. Have you tried using With inside of Flatten so that you can have named variables instead of #? Similarly, the various forms or Sqrt[(#[[2]] - #[[1]])[[1]]^2 + (#[[2]] - #[[1]])[[2]]^2] look like a dot product, have you tried using Norm, instead? –  rcollyer Mar 13 '12 at 15:59

3 Answers 3

Not a proper answer, but I just want to comment that the procedure carried out by @rcollyer can be automated to a large extent. Here is a code for a simplistic common subexpression eliminator:

ClearAll[csub];
csub[expr_Hold, rules_List, limitCount_] :=
  With[{newrule = 
     Replace[
       If[# =!= {} && #[[-1, -1]] > 1, #[[-1, 1]], {}] &@
         SortBy[Tally[
           Cases[expr,          
             x_ /; Depth[Unevaluated[x]] > 2 && 
               LeafCount[Unevaluated[x]] > limitCount :> Hold[x],
             Infinity]],
           Last],
       Hold[x_] :> (HoldPattern[x] -> Unique[])]
    },
    csub[expr /. newrule, Append[rules, newrule], limitCount] /; 
        newrule =!= {}];

csub[Hold[expr_], rules_, _] :=
 Append[
   Thread[
     (rules /. (Verbatim[HoldPattern][x_] -> var_) :> Hold[var := x]), 
     Hold
   ],
   Unevaluated[expr]
  ] /. Hold[defs_, exp_] :> Hold[LetL[defs, Hold[exp]]]

where the LetL macro implements the sequential With construct, defined as

ClearAll[Let];
SetAttributes[Let, HoldAll];
Let /: (lhs_ := Let[vars_, expr_ /; cond_]) := 
   Let[vars, lhs := expr /; cond]
Let[{}, expr_] := expr;
Let[{head_}, expr_] := With[{head}, expr]
Let[{head_, tail__}, expr_] := With[{head}, Let[{tail}, expr]]

I discussed is more here. Using

held = Hold[If[your-expression]]

setting the threshold to say 10 (I played with this a bit), and executing

csub[held,{},10]

we get this:

Hold[
  letL[{
    $3 := #1[[2]] - #1[[1]], $4 := #1[[1, 2]] - #1[[2, 2]], 
    $5 := #1[[2, 1]] - #1[[1, 1]], $6 :=  Sqrt[$3[[1]]^2 + $3[[2]]^2], 
    $7 := (#1[[1, 3 ;; 4]].{$4, $5} {$4, $5})/
       (Sqrt[$4^2 + $5^2] Sqrt[$4^2 + $5^2]), 
    $8 := #1[[2, 3 ;; 4]].$3[[{1, 2}]]
    }, 
    Hold[
       If[(list[[2]] - list[[1]])[[{1, 2}]] != {0, 0}, 
          (
           Flatten[{
             #1[[1, 1 ;; 2]] + $7 + ($8 $3[[{1, 2}]])/($6 $6), 
             $7 + (($8 $6) $3[[{1, 2}]])/$6
           }] &) /@ {{list[[1]], list[[2]]}, {list[[2]],list[[1]]}}, 
         (Flatten[{#1[[{1, 2}]] + #1[[{3, 4}]], #1[[{3, 4}]]}] &) /@ list]
       ]
    ]
  ]

If you now apply ReleaseHold to the above, you recover your code, wrapped in Hold. Note that while I tried automatic minimization of Leafcount for the above, it turns out to give smaller thresholds than what seems to be the easiest for us to read. The real metric here should include some penalty for long chains of dependencies, since they are hard for us humans to digest.

share|improve this answer
    
very nice. this is also relevant for this sort of thing: library.wolfram.com/infocenter/MathSource/3947 –  acl Mar 13 '12 at 17:43
    
or this Experimental`OptimizeExpression[x^(12 Sin[x])*ArcTan[x^(12 Sin[x])]] (not sure if it can be made to work on unnamed functions) –  acl Mar 13 '12 at 17:49
    
@acl One problem with this approach is that LetL can not be used directly inside Compile. It is however possible to modify it so that it only expands to nested With but does not actually execute them. –  Leonid Shifrin Mar 13 '12 at 17:50
    
@acl I tried the Experimental`OptimizeExpression, but it did not seem to help here. And thanks for the reference, this seems quite interesting. –  Leonid Shifrin Mar 13 '12 at 17:51
    
@acl for your example, my optimizer works with a threshold say 5: csub[Hold[x^(12 Sin[x])*ArcTan[x^(12 Sin[x])]], {}, 5]. I could probably come up with some sensible automatic metric, basically a "tilted" version of LeafCount to penalize long dependency chains, but have no time for it now. –  Leonid Shifrin Mar 13 '12 at 18:13

This is not an answer to your question, but it does address some of the issues with your code. In particular, it is just plain unreadable, and unreadable code cannot be maintained in any meaningful manner. If you came back to this even after a week of not using it, you would not understand how it works. Towards that end, I've simplified it quite a bit, just by replacing the common sub-expressions with named variables using With. (A good reference for how the different scoping constructs work, such as With, is this question.)

cf = Compile[{{list, _Real, 2}},
 If[(list[[2]] - list[[1]])[[{1, 2}]] != {0, 0},
  Map[
   With[
    {z = #[[1, 3 ;; 4]], w = #[[2, 3 ;; 4]], 
     x = #[[1, ;; 2]], y = #[[2, ;; 2]]},
    With[{unit = (x - y)/Sqrt[(x - y).(x - y)], 
          revunit = {{0, 1}, {-1, 0}}.(x - y)/Sqrt[(x - y).(x - y)]},
      Flatten[
       {x + (z.revunit) revunit + (w.unit) unit, 
       (z.revunit) revunit + (w.unit) unit}
      ]
   ]] &,
   {{list[[1]], list[[2]]}, {list[[2]], list[[1]]}}
  ], 
  Map[Flatten[{#[[{1, 2}]] + #[[{3, 4}]], #[[{3, 4}]]}] &, list]
 ] (* End If *)
] (* End Compile *)

Note my use of nested With statements. Within With, you cannot refer to another variable within the initialization list, and I do not recall if this is true with the other scoping constructs, also. So, I resorted to nesting them so that I could use simple variable names later.

Incidentally, this compiles without any warnings generated by

SetSystemOptions["CompileOptions" -> "CompileReportExternal" -> True] 

To the best of my knowledge, I converted the code correctly, but please check it.

Edit: as an added note, in the above code you will notice that I did not use Norm despite my comment, above. For whatever reason, Norm cannot be compiled, so I resorted to the equivalent and compilable form

Sqrt[(x - y).(x - y)]
share|improve this answer
    
Thank you. I didn't know that "With" can be used in compilation –  Luis Valero Mar 13 '12 at 17:24
    
A big +1 for introducing the concepts of readability and maintainability of code, which so often seem to be overlooked in efforts to enhance performance. –  whuber Mar 13 '12 at 19:33

You have something like Map[Flatten[{[[1, 1 ;; 2]]+stuff]]. This is a syntax error. As mathematica puts it, Syntax::tsntxi: "[[1,1;;2]]" is incomplete; more input is needed.

What are you trying to do, and how did you produce this code?

EDIT: Given the comment to this answer, I have replaced [[1, 1 ;; 2]] by #[[1, 1 ;; 2]]. cF[list] now gives {{8., 10., 7., 8.}, {8., 10., 3., 4.}}. Is this not what you expect? What did you expect and why?

Also, note that this does not actually get compiled. It calls the main evaluator. To see this, evaluate

SetSystemOptions["CompileOptions" -> "CompileReportExternal" -> True]

before evaluating the Compile.

EDIT2:

If I understood correctly, you wish to compile this. This may be achieved as follows:

cF2 = Compile[{{list, _Real, 2}}, 
   If[(list[[2]] - list[[1]])[[{1, 2}]] != {0, 0}, 
    Map[Flatten[{#[[1, 
           1 ;; 2]] + (Take[#[[1, All]], {3, 
               4}].{(#[[1, 2]] - #[[2, 2]]), (#[[2, 1]] - #[[1, 1]])}/
            Sqrt[(#[[1, 2]] - #[[2, 2]])^2 + (#[[2, 1]] - #[[1, 
                   1]])^2]) {(#[[1, 2]] - #[[2, 2]]), (#[[2, 
                1]] - #[[1, 1]])}/
           Sqrt[(#[[1, 2]] - #[[2, 2]])^2 + (#[[2, 1]] - #[[1, 
                  1]])^2] + (Take[#[[2, All]], {3, 
               4}].(#[[2]] - #[[1]])[[{1, 2}]]/
            Sqrt[(#[[2]] - #[[1]])[[1]]^2 + (#[[2]] - \
#[[1]])[[2]]^2]) (#[[2]] - #[[1]])[[{1, 2}]]/
           Sqrt[(#[[2]] - #[[1]])[[1]]^2 + (#[[2]] - #[[1]])[[2]]^2], \
(Take[#[[1, All]], {3, 
               4}].{(#[[1, 2]] - #[[2, 2]]), (#[[2, 1]] - #[[1, 1]])}/
            Sqrt[(#[[1, 2]] - #[[2, 2]])^2 + (#[[2, 1]] - #[[1, 
                   1]])^2]) {(#[[1, 2]] - #[[2, 2]]), (#[[2, 
                1]] - #[[1, 1]])}/
           Sqrt[(#[[1, 2]] - #[[2, 2]])^2 + (#[[2, 1]] - #[[1, 
                  1]])^2] + (Take[#[[2, All]], {3, 
               4}].(#[[2]] - #[[1]])[[{1, 2}]]*
            Sqrt[(#[[2]] - #[[1]])[[1]]^2 + (#[[2]] - \
#[[1]])[[2]]^2]) (#[[2]] - #[[1]])[[{1, 2}]]/
           Sqrt[(#[[2]] - #[[1]])[[1]]^2 + (#[[2]] - \
#[[1]])[[2]]^2]}] &,(*End Flatten*){{list[[1]], 
       list[[2]]}, {list[[2]], list[[1]]}}],(*End Map;
    End True branch of If*)
    Map[Flatten[{#[[{1, 2}]] + #[[{3, 4}]], #[[{3, 4}]]}] &, 
     list]] (*End If*),
   CompilationTarget -> "C"];

(I have replaced things like #[[1,3;;4]] by Take[#[[1]],{3,4}]). However, this is not much faster than the first version (uncompiled). Maybe if you explain how you got this we can help more.

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.