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.

Compile fun continues today, and I'm trying to get the function below to compile. First the function:

NewC[ymax_, xmax_] := Block[{x1, a2, a0, b0, a1, b1, c1, b2, c2, a3, b3, c3, C0, xb = 0, xe = 100, ve},
  a3 = (xmax + ymax - 100)/(2*(xe - xmax)*(xe - xmax));
  x1 = (200 - xmax - 2*ymax + 2*a3*xe*(xe - xmax))/(2*a3*(xe - xmax) + 1);
  a2 = -(2*a3*(xe - x1) + 1)/(2*(x1 - xmax));
  b2 = -2*a2*xmax;
  b3 = -1 - 2*a3*xe;
  a0 = -(xb + ymax)/((xmax - xb)*(xmax - xb));
  b0 = -2*a0*xmax;
  C0 = ymax + a0*xmax*xmax;
  a1 = (100 - xe - ymax)/((xe - xmax)*(xe - xmax));
  b1 = -2*a1*xmax;
  c1 = ymax + a1*xmax*xmax;
  c2 = ymax + a2*xmax*xmax;
  c3 = a3*xe*xe + 100;
  ve2[ym_, tnlr_] := If[ym <= ((100 - xmax) + (100 - xe))/2,
  If[tnlr <= xmax, a0*tnlr*tnlr + b0*tnlr + xb, a1*tnlr*tnlr + b1*tnlr + c1],
  If[tnlr <= x1, a2*tnlr*tnlr + b2*tnlr + c2, a3*tnlr*tnlr + b3*tnlr + c3]];
  ve = Chop[Map[ve2[ymax, #[[1]]] &, $range] // N]
 ]

So pretty straight forward, except for the embedded function. So I externalize the important parts and global variable, but have no idea how to handle the embedded function.

Here's what I tried:

NewC = With[{$range = $range, xb = 0, xe = 100},
   Compile[{{ymax, _Real, 1}, {xmax, _Real, 1}},
    Block[{x1, a2, a0, b0, a1, b1, c1, b2, c2, a3, b3, c3, C0, ve,ve2},
   a3 = (xmax + ymax - 100)/(2*(xe - xmax)*(xe - xmax)); 
   x1 = (200 - xmax - 2*ymax + 
     2*a3*xe*(xe - xmax))/(2*a3*(xe - xmax) + 1);
   a2 = -(2*a3*(xe - x1) + 1)/(2*(x1 - xmax));
   b2 = -2*a2*xmax;
   b3 = -1 - 2*a3*xe;
   a0 = -(xb + ymax)/((xmax - xb)*(xmax - xb));
   b0 = -2*a0*xmax;
   C0 = ymax + a0*xmax*xmax;
   a1 = (100 - xe - ymax)/((xe - xmax)*(xe - xmax));
   b1 = -2*a1*xmax;
   c1 = ymax + a1*xmax*xmax;
   c2 = ymax + a2*xmax*xmax;
   c3 = a3*xe*xe + 100;
   ve2[ym_, tnlr_] := If[ym <= ((100 - xmax) + (100 - xe))/2,
     If[tnlr <= xmax, a0*tnlr*tnlr + b0*tnlr + xb, a1*tnlr*tnlr + b1*tnlr + c1],
     If[tnlr <= x1, a2*tnlr*tnlr + b2*tnlr + c2,a3*tnlr*tnlr + b3*tnlr + c3]];
   ve = Chop[Map[ve2[ymax, #[[1]]] &, $range] // N]
   ],
  CompilationTarget -> "C"]];

Certainly not the winning ticket. So are there principles to follow when this type of function must be modified? Thanks!

share|improve this question
    
What exactly is $range and does it ever change? Should it be an input to the function or a global variable? –  Szabolcs Dec 18 '13 at 23:09
    
$range is a list of values from 0-1 and is a global variable used in other functions as well. Thanks! –  R Hall Dec 18 '13 at 23:17
    
You are Maping #[[1]] onto $range so I assumed that it was a matrix, not a vector (or it was a vector of vectors). Also, I just noticed Curve2 vs ve2. Is this a typo and are these two the same? I assumed that they are in the answer. –  Szabolcs Dec 18 '13 at 23:19

1 Answer 1

up vote 2 down vote accepted

This is not a general solution, but it should apply to your example. For a general solution the question is going to be: what does the function depend on?

First let's get a few things out of the way: $range appears to be a global variable, possibly a matrix where you only use the first column. Instead let's make range represent only the first column (a simple vector of numbers) and let's make it an input to the function.

The key operation is the mapping. It map a pre-constructed function. We can exploit the fact that Function objects can usually be compiled without problem, and write the function in that form. Also, we must avoid assigning the function to a variable name and write it directly in the Map.

Here's the result:

cf = Compile[
  (* the inputs; range is also an input, I assume it's a vector of reals *)
  {{ymax, _Real}, {xmax, _Real}, {range, _Real, 1}},

  (* no changes here: *)  
  Block[{x1, a2, a0, b0, a1, b1, c1, b2, c2, a3, b3, c3, C0, xb = 0, xe = 100, ve}, 
   a3 = (xmax + ymax - 100)/(2*(xe - xmax)*(xe - xmax));
   x1 = (200 - xmax - 2*ymax + 2*a3*xe*(xe - xmax))/(2*a3*(xe - xmax) + 1);
   a2 = -(2*a3*(xe - x1) + 1)/(2*(x1 - xmax));
   b2 = -2*a2*xmax;
   b3 = -1 - 2*a3*xe;
   a0 = -(xb + ymax)/((xmax - xb)*(xmax - xb));
   b0 = -2*a0*xmax;
   C0 = ymax + a0*xmax*xmax;
   a1 = (100 - xe - ymax)/((xe - xmax)*(xe - xmax));
   b1 = -2*a1*xmax;
   c1 = ymax + a1*xmax*xmax;
   c2 = ymax + a2*xmax*xmax;
   c3 = a3*xe*xe + 100;

   (* embed the function in the Map directly, as a Function[...] expression *)
   Chop@Map[
    Function[tnlr,
     If[ymax <= ((100 - xmax) + (100 - xe))/2, 
       If[tnlr <= xmax, a0*tnlr*tnlr + b0*tnlr + xb, 
          a1*tnlr*tnlr + b1*tnlr + c1], 
       If[tnlr <= x1, a2*tnlr*tnlr + b2*tnlr + c2, 
          a3*tnlr*tnlr + b3*tnlr + c3] ]],

    range]
  ] (* end Block *)
 ] (* end Compile *)

This compiles without issues.

share|improve this answer
    
If you needed to map this function multiple times, and thus couldn't avoid naming it, compiling it would be a more difficult problem ... –  Szabolcs Dec 18 '13 at 23:21
    
Looks like the last line of the function is missing, can I just add it in as is? –  R Hall Dec 18 '13 at 23:23
    
@RHall The last two lines are hidden in the code block, did you scroll down? It looks pretty confusing on a Mac where the scrollbars are not shown ... EDIT: I removed some whitespace to make it shorter now. –  Szabolcs Dec 18 '13 at 23:23
    
@RHall Chop is okay because it is compilable, though you could just apply Chop outside of the compiled function. But N should not be used. It is not compilable and it is not needed here. range is going to be converted to machine precision numbers (the same thing N does) anyway before it's passed to the compiled function internally and the result is also going to be a packed machine precision array. –  Szabolcs Dec 18 '13 at 23:38
    
Correction to "N is not compilable": the compiler recognizes that N is not needed and will automatically remove it if you do add it. –  Szabolcs Dec 18 '13 at 23:39

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.