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.

In this great answer a compiled version of the Nelder-Mead algorithm is presented.

Since it works on arbitrary dimensions (i.e. arbitrary number of arguments), it has to use apply on the objective function. The problem is that Apply is not directly supported inside Compile. To overcome this limitation the following code is used:

(* Produces compiled code for the Nelder-Mead algorithm with the objective function inlined *)
ClearAll[apply];
SetAttributes[apply, HoldAll];
apply[f : (_Function | _CompiledFunction), vars : {__Symbol}] :=
    With[{applied := f @@ vars},
      Function[arglist, Block[vars, vars = arglist; applied]]
     ];

This seems to work (inside the package), but I have no idea, how it works.

Could somebody explain the techniques behind this code snippet? Specifically:

  • What does SetDelayed do inside With?
  • Why does vars appear twice in the arguments to Block?

  • Why doesn't this minimal example work? (I believe it mimics what is done in the mentioned answer...)

    Clear[a, b, c, x, y, z, objectiveFunction, cfunc];
    objectiveFunction = Compile[{a, b, c, x, y, z} ,
       (a - x)^2 + 50 (b - y)^2 + (c - z)^2];
    
    cfunc = With[{f = apply[objectiveFunction, {a, b, c, x, y, z}]},
      Compile[ {{pts, _Real, 1}}, f@pts]
      ]
    << CompiledFunctionTools`
    CompilePrint[cfunc]
    

EDIT:

Here is the working snippet. Note that "InlineCompiledFunctions" must be set to True to avoid the call to MainEvaluate.

ClearAll[a, b, c, x, y, z, apply, objectiveFunction, cfunc];
(* Produces inlinable code for use inside Compile (where Apply is not \
supported directly) *)
SetAttributes[apply, HoldRest];
apply[f : (_Function | _CompiledFunction), vars : {__Symbol}] :=
       Function[arglist, Block[vars, vars = arglist; f @@ vars]];

objectiveFunction = Compile[{a, b, c, x, y, z} ,
   (a - x)^2 + 50 (b - y)^2 + (c - z)^2];
apply[objectiveFunction, {a, b, c, x, y, z}]
cfunc = With[{f = apply[objectiveFunction, {a, b, c, x, y, z}]},
  Compile[ {{pts, _Real, 1}}, f@pts,
    CompilationOptions -> {"InlineCompiledFunctions" -> True}]
  ]
<< CompiledFunctionTools`
CompilePrint[cfunc]
share|improve this question
    
This makes me wonder why is the a system option called ApplyCompileLength and why is it set to Infinity –  Szabolcs May 4 '12 at 9:12
2  
@Szabolcs This is an entirely different story. Apply as discussed here evaluates at (pre) Compile-time and is a tool of code generation. The option you mention governs the auto-compilation of Apply, when it would be explicitly present in the Compile body when auto-compiling. Since Apply is only supported for 3 heads (List, Plus and Times) by Compile, it does not generally make sense to set the mentioned option to a finite number. In case you have one of those heads, you can of course reset it. –  Leonid Shifrin May 4 '12 at 9:18
    
@Leonid (It was really a question, not a comment...) Is it documented somewhere that it's supported for these three heads? Do you know why that option is set to Infinity? –  Szabolcs May 4 '12 at 9:27
    
@Szabolcs It is documented, but I can't tell you where off hand. For the option setting, Infinity is equivalent to saying "don't attempt to auto-compile at all". The fact that only 3 heads are supported means that Compile does not generally support Apply, which is probably the reason for this setting. For those heads, it can be reset manually, as I mentioned. –  Leonid Shifrin May 4 '12 at 9:33

1 Answer 1

up vote 10 down vote accepted

First question : With accepts a syntax like

  With[{var:=value}, expression]

in which case, value is injected into expression unevaluated. As far as I know, this syntax is not documented. You can achieve a similar effect with the replacement rules, by using

Unevaluated[expression]/.HoldPattern[var]:>value

There are some subtle differences between the semantics of With and repalcement rules though, mostly related to the treatment of nested scoping constructs and variable name conflicts in them.

Second question: vars appear twice because they must first be Block-ed, and then there is a massive assignment to them performed in the body of the Block. This is probably the most economical way of blocking a number of variables and assigning values to them simultaneously - otherwise a more complex code-generation will be needed. You can see another example of that in this answer (and if you look at the revision history for that answer, you can find an alternative, harder way to do this, in one of the previous revisions).

Third question: this does not work because the apply function was made HoldAll (which isn't quite necessary), and the pattern-matching does not work. There were some past discussions on this topic on SO, but can't find them right now. But I discussed this topic at length also in my book. The idea is that at the pattern-matching time, all seen by apply is a variable objectiveFunction, and because it does not evaluate it, the pattern _CompiledFunction is not matched. The solution is to make apply HoldRest, and then it works:

ClearAll[apply];
SetAttributes[apply, HoldRest];
apply[...]:=...
share|improve this answer
    
Aha, so this works by wrapping the original function in another function which takes just a single argument (the list) and then mass assigning back to individual arguments? When is the Apply used in With evaluated? –  Ajasja May 4 '12 at 8:55
    
Or to put it another way, why is With even necessary, since apply[f : (_Function | _CompiledFunction), vars : {__Symbol}] := Function[arglist, Block[vars, vars = arglist; f @@ vars]]; seems to work as well. –  Ajasja May 4 '12 at 8:57
    
@Ajasja Only inside Block, and that's the entire point, because by that time, Block guarantees that the variables have no global values. –  Leonid Shifrin May 4 '12 at 8:58
3  
@OleksandrR. No problem :) You use some really nice stuff, which is a lot of fun to read and discuss. By the way, I take this opportunity to say that I am (among many others) quite impressed by that answer of yours which originated this question. I think we need more answers like that one, to show that one can often use Mathematica as an implementation language for some missing functionality, and that the absence or limitations of some built-in function don't mean that "Mathematica can not do this", or "can not do this fast enough". I think you did a great job proving such statements wrong. –  Leonid Shifrin May 4 '12 at 10:44
1  
That SetDelayedish syntax in With is really useful. It's a shame it's not documented. +1! –  Pillsy May 4 '12 at 18:46

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.