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 considered the following function:

sin[x_] := Module[{},
  Print["x=", x];
  Sin[x]
]

in Mathematica. Next, I tried to plot it using:

Plot[sin[t], {t, 0, 2 Pi}]

Surprisingly, the first three lines of output are:

x=0.000128356
x=t
x=1.28228*10^-7

Can someone explain this behavior? In this case it doesn't cause a problem, but in my "real" case it does.

Summary

acl's answer below offers, at its very beginning a solution to the specific problem. In very-short, the reason that this x=t appears is hidden somewhere in the way Mathematica evaluates the functions. The answers below provide interesting insight into the way it works.

The interested reader should read all the answers and details below, they are invaluable, although might be behind the reach of some of the readers (like, partially, in my case).

share|improve this question

migrated from stackoverflow.com Jul 28 '13 at 6:31

This question came from our site for professional and enthusiast programmers.

2 Answers 2

up vote 19 down vote accepted

If the problem is that a symbolic argument is passed, you can avoid it thus:

ClearAll[sin];
sin[x_?NumericQ] := Module[{},
  Print[x];
  Sin[x]
  ]

which simply defines sin so that it only matches for numeric arguments.

To see what it does, try sin[3.] and sin[x] and notice that the second evaluates to itself, as the definition above does not match.

You can also see what values of x are being evaluated by ClearAll[sin]; sin[x_] := Module[{}, Sow[x]; Sin[x] ] and then Plot[sin[x], {x, 0, 10}]; // Reap. x now appears.

However,

lst = {};
ClearAll[sin2];
sin2[x_] := Module[{},
  AppendTo[lst, x];
  Sin[x]
  ]

and then Plot[sin2[x], {x, 0, 10}]; and then we have no symbols in lst at the end.

EDIT: The explanation for this discrepancy between Sow/Reap and using a list is explained by Leonid in the comments. To test his proposal, I tried using a Bag instead of a list (this is undocumented, see Daniel Lichtblau's description) as follows:

AppendTo[$ContextPath, "Internal`"];
lst = Bag[];
sin3[x_] := Module[{},
  StuffBag[lst, x];
  Sin[x]
]

followed by Plot[sin3[x], {x, 0, 10}];. We now inspect the contents of the bag by BagPart[lst, All] and observe that there is indeed a symbol x in there.

Presumably it has to do with the way scoping constructs interact with the evaluations performed by AppendTo and StuffBag.

EDIT 2 (by Leonid Shifrin)

We can also demonstrate the same using more usual tools. In particular, instead of a Bag that has its own API, we can use any HoldAll wrapper (just not a list), and then the code for the function itself we need not change at all:

In[51]:= 
ClearAll[h];
SetAttributes[h,HoldAll];
lst=h[];
ClearAll[sin2];
sin2[x_]:=Module[{},AppendTo[lst,x];Sin[x]]

In[58]:= 
Plot[sin2[x],{x,0,10}];
lst//Short

Out[59]//Short= h[0.000204286,x,<<1131>>,9.99657]

This clarifies what happens. The x inside List is substituted by the numerical value as a result of evaluation in AppendTo, roughly as follows:

In[60]:= 
Clear[x];
lst = {0.000204,x};
Block[{x = 2.04*10^(-7)},
  AppendTo[lst,x]];
lst

Out[63]= {0.000204,2.04*10^-7,2.04*10^-7}

while the HoldAll attribute of h prevents the evaluation from happening (this will be more clear yet if we write AppendTo as lst = Append[lst,x]. It is the evaluation of the r.h.s. (Append), where lst is evaluated and x is substituted by its bound value). For h, x inside it does not evaluate, and is therefore kept symbolic. Similar thing happens with Reap-Sow, although the mechanism Reap-Sow is using to store the results is obviously different (but, whatever it is, it bypasses the main evaluation loop, and that is what matters).

EDIT 3 (acl):

There was a question in the comments as to why the numbers returned by Sow/Reap are not in ascending order. The reason is that Plot apparently uses an adaptive algorithm, in the same spirit as one does in adaptive integration (see en.wikipedia.org/wiki/Adaptive_quadrature, for instance).

Do Plot[sin[x], {x, 0, 10}]; // Reap // Last // Last // ListPlot to see it spend more effort at the turning points:

enter image description here

If you add the option MaxRecursion -> 0 to the Plot command, the algorithm does not subdivide steps that it deems inaccurate and the values are in order:

enter image description here

Maybe it is clearer to do it interactively. Let us play with MaxRecursion and PlotPoints:

ClearAll[sin];
sin[x_?NumericQ] := (Sow[x]; Sin[x])

Manipulate[
 pts = ((plt = Plot[
          sin[x],
          {x, 0, 10},
          PlotStyle \[Rule] {Red, Thin},
          PlotPoints \[Rule] n,
          MaxRecursion \[Rule] m
          ];) // Reap // Last // Last);

 Show[
  {
   ListPlot[
    Transpose@{pts, Sin[pts]},
    PlotMarkers \[Rule] {Automatic, 3}
    ],
   plt
   }
  ],
 {m, Range[0, 5]},
 {{n, 10}, Range[1, 50]}
 ]

m is the value of MaxRecursion, n that of PlotPoints. The plot shows the resulting plot of Sin and, overlaid, the points that have been evaluated to produce it. Play with the numbers and it should become clear what is happening: PlotPoints tells Plot how many points to evaluate initially, MaxRecursion tells Plot how many times it may subdivide the regions thus defined if necessary (see here for a discussion of what "necessary" means).

enter image description here

share|improve this answer
1  
Sounds like there is more added to the puzzle than solved :D The Evaluated option of Plot is of relevance as well (it can be True or False or Automatic). I can't explain your results either. –  Szabolcs Aug 26 '11 at 12:54
    
@Szabolcs yes I tried that too (I believe I first heard about this option from you). Didn't help to explain much, so I thought I'd avoid further muddling the waters. –  acl Aug 26 '11 at 13:15
    
can you reproduce the OPs results? –  rcollyer Aug 26 '11 at 14:11
2  
@Sjoerd, rcollyer Here is what I think is happening. The crucial thing is that at some point, the binding of symbol x to a numerical value is happening. But first, symbolic argument x is tried (presumably as a part of an attempt to symbolically preprocess the function). At that point, symbol x does get into the list. To see that, use sin2[x_] := Module[{}, AppendTo[lst, x]; If[MatchQ[lst, {___, _Symbol, ___}], Print[lst]]; Sin[x]] (for example). Then, I suspect, Plot us using dynamic scoping (like Block) to localize the variable x. Therefore, once numeric computation ... –  Leonid Shifrin Aug 26 '11 at 15:49
3  
@Sjoerd continuing ... starts, x gets dynamically bound to a given numerical value (first, the first one, etc). The next AppendTo, being within that dynamic scope, evaluates lst = {firstnumber, x} to lst = {firstnumber, secondnumber}, where secondnumber is what x is bound to now. After that, we get purely numerical list of course. The reason that this does not happen with Reap- Sow is that they use a different mechanism to store results, which is immune to the dynamic scope imposed by Plot. At least, this is what I think is happening. –  Leonid Shifrin Aug 26 '11 at 15:53

It is interesting to compare the Plot algorithms of Mathematica 5.2 and Mathematica 6+. Based on acl's code:

In Mathematica 5.2 we get:

Plot[Sow[x]; Sin[x], {x, 0, 10}, 
    DisplayFunction -> (Null &)] // Reap // Last // Last // ListPlot

Mathematica 5.2 output

In Mathematica 7.0.1:

Plot[Sow[x]; Sin[x], {x, 0, 10}] // Reap // Last // Last // ListPlot

Mathematica 7.0.1 output

One can see that Mathematica 5.2 computes only 105 points while version 7.0.1 computes 567 points (including one symbolic evaluation). At the same time, the plots produced by two versions are visually indistinguishable. Only very careful comparison reveals tiny differences. Here is ListPlot of both sets of points: generated by version 5.2 (dashed black line) and version 7.0.1 (blue line) with resolution 600 dpi (click to enlarge!):

comparison



Edit

In Mathematica 5.2 the default value for PlotPoints is 25 and for PlotDivision is 30 as documented in Section 1.9.2 of "The Mathematica Book" for version 5.2. I do not know where the default value for PlotPoints in version 7 is documented but we can find it by setting MaxRecursion to zero:

In[1]:= Cases[Plot[x, {x, 0, 10}, MaxRecursion -> 0], 
                            Line[x_] :> Length[x], Infinity]

Out[1]= {50}

From the other side, using Reap and Sow we get different value:

In[2]:= Select[
  Plot[Sow[x], {x, 0, 10}, MaxRecursion -> 0]; // Reap // 
    Last // First, NumericQ] // Length

Out[2]= 51

For Mathematica 5.2 both methods give the same result (25). It seems that in Mathematica 7.0.1 the first point is calculated just for the check that the objective function gives numerical value for numerical argument but this point is not included in the final plot:

In[4]:= Complement[
 Plot[Sow[x], {x, 0, 10}, MaxRecursion -> 0, PlotPoints -> 25] // 
    Reap // Last // First, 
 Cases[Plot[x, {x, 0, 10}, MaxRecursion -> 0, PlotPoints -> 25], 
   Line[x_] :> x, Infinity][[1, All, 1]]]

Out[4]= {0.000417083, x}

Edit 2

In Mathematica 7 increasing MaxRecursion just adds new layers of points to the ListPlot:

v7points[r_] := 
  Module[{i = 1}, 
   Last@Last@
     Reap@Plot[Sow[{i++, x}]; Sin[x], {x, 0, 10}, PlotPoints -> 25, 
       MaxRecursion -> r]];
v7plot = ListPlot[Join[{v7points[0]}, 
   Complement[v7points[# + 1], v7points[#]] & /@ Range[0, 10]], 
  PlotMarkers -> (Graphics`PlotMarkers[] /. {m_, s_} :> {m, s/2}), 
  PlotStyle -> ColorData[60, "ColorList"]]]]

screenshot

In Mathematica 5.2 we have PlotDivision instead of MaxRecursion:

v5Points[d_] := 
 krn5Eval[Module[{i = 1}, 
   Last@Last@
     Reap@Plot[Sow[{i++, x}]; Sin[x], {x, 0, 10}, PlotPoints -> 25, 
       PlotDivision -> d, DisplayFunction -> (Null &)]]]
v5plot = ListPlot[v5Points[#] & /@ Range[1, 9], 
  PlotMarkers -> (Graphics`PlotMarkers[] /. {m_, s_} :> {m, s/2}), 
  PlotStyle -> ColorData[60, "ColorList"]]

(here krn5Eval[expr] is a MathLink function which evaluates expr in the kernel of Mathematica 5 from inside of Mathematica 7)

screenshot

From the point of view of the number of evaluations in the case of plotting of Sin[x] PlotPoints -> 2, PlotDivision -> 30 is roughly equivalent to PlotPoints -> 2, MaxRecursion -> 5. So we can compare:

v7Ps = Last@
   Last@Reap@
     Plot[Sow[{x, Sin[x]}]; Sin[x], {x, 0, 10}, PlotPoints -> 2, 
      MaxRecursion -> 5];
v5Ps = Last@
   Last@krn5Eval[
     Reap@Plot[Sow[{x, Sin[x]}]; Sin[x], {x, 0, 10}, PlotPoints -> 2, 
       DisplayFunction -> (Null &), PlotDivision -> 30]];
ListLinePlot[Sort /@ {v7Ps, v5Ps}, 
 PlotStyle -> {Blue, {Black, Dashed}}]

plot (click to enlarge!)


Edit 3

Both algorithms have the same drawback: adaptive refinement is made by minimizing the angles at the nodes connecting the segments of the approximating polyline, but if a corner was once close enough to 180 degrees (180 degrees minus this angle is less than the corresponding parameter: MaxBend or ControlValue), the corresponding node is dropped and no longer checked. Here is an illustration of what happens in version 5.2:

screenshot

The same case for Mathematica 6+ was already investigated by Yaroslav Bulatov: "Strange Sin[x] graph in Mathematica."

In such cases increasing of MaxRecursion in MMa 6+ does nothing since the node is already dropped from the list of nodes to check. In Mathematica 5 the problem is more subtle: changing any of the control parameters (PlotPoints, MaxBend, PlotDivision) shifts all sample points and as a result the problematic node disappear but now it may emerge in another place. And increasing PlotDivision will not reduce the probability to face this problem again if you have already faced it. The only reliable solution is to considerably increase PlotPoints.

Edit 4

The MaxBend option of Plot in Mathematica 5 has a completely equivalent undocumented sub-sub-option ControlValue in Mathematica 6+ with only difference: the latter should be specified in radians while the former in degrees. At the same time, Mathematica 6+ still has the old MaxBend option moved inside of Method option. I have found it accidentally by evaluating

Plot[x,{x,0,1},MaxBend->1,PlotDivision->1];
MaxBend::deprec: MaxBend->1 is deprecated and will not be supported in future versions of Mathematica. Use Method->{MaxBend->1} instead.
PlotDivision::deprec: PlotDivision->1 is deprecated and will not be supported in future versions of Mathematica. Use Method->{PlotDivision->1} instead. >>

I have tested it and found that the following two ways to specify MaxBend are completely equivalent in Mathematica 7.0.1 and 8.0.4:

In[1]:= With[{maxBend = 5}, 
 First@Plot[Sin[x], {x, -42 Pi, 42 Pi}, 
    PlotRange -> {{-110, -90}, All}, Method -> {MaxBend -> maxBend}] ===
   First@Plot[Sin[x], {x, -42 Pi, 42 Pi}, 
    PlotRange -> {{-110, -90}, All}, 
    Method -> {"Refinement" -> {"ControlValue" -> maxBend*\[Degree]}}]]

Out[1]= True

When these options are specified together the ControlValue is used.

Note that MaxBend in version 5 by has default value 10. (degrees) while in version 6+ it has default value 5*Degree (radians). The less this value - the more precise plot will be generated, so in really it is not correct to compare these algorithms with no attention to this option.

One important feature (bug?) of Plot in Mathematica 6+ is that it does not stop adding new levels of recursion when the "ControlValue" condition is already satisfied:

l[mr_] := 
 Length@Reap[Plot[Sow[x], {x, -Pi, Pi}, MaxRecursion -> mr]][[2, 1]]
ListPlot[l /@ Range[1, 15], Axes -> False, Frame -> True, 
 FrameLabel -> {"MaxRecursion", "Number of evaluation points "}, 
 PlotLabel -> StandardForm@HoldForm[Plot[x, {x, -Pi, Pi}]]]

screenshot

At the same time, in Mathematica 5 Plot stops recursion when bend angles become less than MaxBend:

screenshot

share|improve this answer
    
very cool, +1. could you repeat this with both versions and with explicit setting for PlotPoints and MaxRecursion? –  acl Aug 30 '11 at 12:52
    
It looks like v5.2 steps through the range until it finds interesting (big 2nd derivative) parts - in this case the turning points - it then investigates that part in more detail until moving on. Versions >= 6 do an even sampling (given by PlotPoints) before sectioning the graphics into interesting areas to look closer at. I wonder if the v>=6 algorithm has more points because computers are assumed to be better (what's the default for MaxRecursion in v5.2?) or if it can afford more points because it's inherently faster? Which algorithm gives the better value for cycles? –  Simon Aug 30 '11 at 12:53
1  
@acl Please see updated answer. –  Alexey Popkov Sep 1 '11 at 11:51
    
Upvoted your comment, since I can't upvote your answer twice! –  Simon Sep 1 '11 at 12:35
    
@Alexey thanks. v8 shows apparently the same behaviour as v7. –  acl Sep 1 '11 at 12:36

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.