5
$\begingroup$

Q: Is there a general way to remove particular functions from a previously stored call to a plot function?

Here is a specific example:

plot1 = Plot[{x, x^2}, {x, 0, 2}];
plot2 = plot1;
plot3 = plot1;
plot2[[1, 1, 1, 3]] = {};
plot3[[1, 1, 1, 4]] = {};
{plot1, plot2, plot3}

Plot dissection

But I'd like to do it on other plot functions, like Plot3D.

Here's an attempt at that:

plot1 = Plot3D[{x, x^2}, {x, 0, 2}, {y, 0, 2}];
plot2 = plot1;
plot3 = plot1;
plot2[[1, 2, 1, 1]] = {};
plot3[[1, 2, 1, 2]] = {};
{plot1, plot2, plot3}

Plot3D dissection

You can see that I can get rid of the most of the plots of x and x^2 individually again, but now the mesh remains. Perhaps I could poke around inside the Plot3D object some more and get rid of those, but you may also notice that the relevant indices for Plot and Plot3D are different. Rather than doing microsurgery on my plots every time I want to do this, I was hoping there was a general method to do this cleanly.

Ideally, the solution should also work on something like:

mat={{x^2+y^2-1,0},{0,-x^2-y^2+1}};
plot1=Plot3D[Sort[Eigenvalues[mat]],{x,-2,2},{y,-2,2}]

(but with a way bigger matrix!)

Here's why I care: I have a reasonably nasty matrix (216x216) which is a function of two variables. Only numerical methods have any hope of getting its eigenvalues, and it's sufficiently nasty that Eigenvalues does not return the eigenvalues in any particular order. That means that I have to call Sort on the results of Eigenvalues before plotting in order to get anything that doesn't look like it's been run through the wood-chipper.

But that means that if I just want to plot one band, I first have to get the numerical eigenvalues corresponding to that x and y for ALL bands.

Well, I really want to plot all bands, just not all on the same graph. Now, if I plot each band separately, it's going to find all those eigenvalues every single time. That's a huge waste of time, since in principle, Mathematica has all the info it needs on the first call.

So I would like to make one call to Plot3D on the full spectrum of eigenvalues which I would store like in the examples above, and then make copies of that which have all but one or a few bands removed from the plot.

Thanks for getting into the nitty gritty with me!

$\endgroup$
  • 1
    $\begingroup$ unlike polygons mesh lines are note separated into groups. So it is difficult to eliminate the right set of lines. To get just the surfaces (without the lines) you can do plot2 =Normal@ plot1;plot3 =Normal@ plot1;plot2[[1,1, 2]] = plot2[[1,1, 1,2]] = {};plot3[[1,1,2]] = plot3[[1,1,1,1]] = {}; $\endgroup$ – kglr Jul 18 at 5:03
  • $\begingroup$ why not move Sort@Eigenvalues@mat outside Plot3D, that is, use, ev = Sort[Eigenvalues[mat]];plots = Plot3D[#, {x, -2, 2}, {y, -2, 2}] & /@ ev? or plots = Plot3D[#, {x, -2, 2}, {y, -2, 2}, PlotStyle -> #2] & @@@ Transpose[{ev, {Red, Blue}}]? $\endgroup$ – kglr Jul 18 at 5:07
  • 2
    $\begingroup$ if it is ok to have just the surfaces you can also do : {plot2,plot3} =Graphics3D/@Cases[Normal[plot1], {__,_GraphicsGroup}, All] $\endgroup$ – kglr Jul 18 at 5:20
  • 1
    $\begingroup$ @kglr, as OP explained, it is only possible to compute ev numerically, so simply moving Sort@Eigenvalues@mat outside of Plot3D is hardly possible. One thing to try is to precompute Eigenvalues in Table-like way (on a mesh of x and y with some small enough steps) and then to ListPlot3D these Sorted bands. Or even make then some Interpolation. $\endgroup$ – Alx Jul 18 at 8:34
  • 1
    $\begingroup$ @kglr, thanks for the extra info on the guts of Plot3D. I do in fact move the Sort@Eigenvalues@mat outside of the plot in my own code, but as @Alx points out, that serves only for notational convenience. After asking this question yesterday, I went the route of ListPlot3D. It looks at least as good as Plot3D, and it is MUCH faster. $\endgroup$ – thejacksonjive Jul 18 at 20:55
6
$\begingroup$
plot1 = Plot3D[{x, x^2}, {x, 0, 2}, {y, 0, 2}, ImageSize -> Medium];

We get the surfaces by extracting GraphicsGroups:

surfaces = Graphics3D[#, Options@plot1]&/@Cases[Normal[plot1], {__, _GraphicsGroup}, All];

Getting correct set of lines for each surface is challenging. The following seems to work (not sure how robust it is in general):

lines = Cases[Normal[plot1], _Line, All];
partition = Partition[lines[[1 + Length @ surfaces ;;]], 15];
meshlines = Graphics3D /@ MapIndexed[Flatten[Prepend[#, lines[[#2]]]] &, 
    Transpose[Partition[partition, Length @ surfaces]]];

Row[Prepend[MapThread[Show, {surfaces, meshlines}], plot1]]

enter image description here

Note: the magic number 15 above is the default number of mesh lines in x and y directions. The lines are organized as follows: if there a k surfaces, the first k lines are the boundary lines for the k surfaces. The next 15 k lines are mesh lines in x direction and the last 15 k lines are mesh lines in y direction. In each group, the first 15 lines belong to the first surface, the second 15 lines belong to the second surface etc.

For

plot1 = Plot3D[{x, x^2, 2+Sqrt[x]}, {x, 0, 2}, {y, 0, 2}, ImageSize -> 300];

we get

enter image description here

For

mat = {{x^2 + y^2 - 1, 0}, {0, -x^2 - y^2 + 1}};
plot1 = Plot3D[Evaluate[Eigenvalues[mat]],  {x, -2, 2}, {y, -2, 2}, 
  PlotRange -> All, ImageSize -> Medium];

we get

enter image description here

If we use Sort[Eigenvalues[mat]] in the previous example we get a single surface and the approach breaks down. Using Evaluate[{RankedMin[#,1],RankedMin[#,2]}&[Eigenvalues[mat]]] instead of Sort[Eigenvalues[mat]]

plot1 = Plot3D[Evaluate[{RankedMin[#,1],RankedMin[#,2]}&[Eigenvalues[mat]]],  
  {x, -2, 2}, {y, -2, 2}, PlotRange -> All, ImageSize -> Medium];

we do get multiple surfaces but we need to change the magic number 15 to 31 to get

enter image description here

So for a more robust approach:

ClearAll[separatePlots]
separatePlots = Module[{opts = Options[#], 
     gc = Cases[#, GraphicsComplex[c_, ___] :> c, All][[1]], 
     gcsurfaces = Cases[#, {___, _GraphicsGroup}, All], 
     gclines = Cases[#, _Line, All], pts, surfaceMeshAssoc},
    pts = DeleteDuplicates[Flatten[Cases[#, Polygon[x_, ___] :> x, All][[1]]]] & /@ 
      gcsurfaces;
    surfaceMeshAssoc = GroupBy[Function[x, 
     {First @ MaximalBy[Range[Length@gcsurfaces], 
       Length[Intersection[pts[[#]], x[[1]]]] &], x}] /@ gclines, First -> Last];
    Prepend[Graphics3D[GraphicsComplex[gc,
      {surfaceMeshAssoc[#], gcsurfaces[[#]]}], opts] & /@ Range[Length@ gc surfaces], #]]&;

Examples:

Row @ separatePlots @ Plot3D[{x, x^2, 2 + Sqrt[x]},  {x, 0, 2}, {y, 0, 2}, 
  PlotRange -> All, ImageSize -> 300]

enter image description here

Row @ separatePlots @ Plot3D[Evaluate[{RankedMin[#, 1], RankedMin[#, 2]} &[
     Eigenvalues[mat]]],  {x, -2, 2}, {y, -2, 2},  
    PlotRange -> All, ImageSize -> Medium, Exclusions -> None]

enter image description here

mat2 = {{-1 + x^2 + y^2, 0, x y }, {0, 1 - x^2 - y^2 , x + y}, {x y, 
    x + y, 1 - x^2 - y }} ;
Row @ separatePlots @ Plot3D[
  Evaluate[{RankedMin[#, 1], RankedMin[#, 2], RankedMin[#, 3]} &[Eigenvalues[mat2]]],  
  {x, -2, 2}, {y, -2, 2},   PlotRange -> All, ImageSize -> 300]

enter image description here

Reap/Sow

An alternative approach is to use Reap/Sow to get the evaluations inside Plot3D and use ListPlot3D on the reaped evaluations:

Examples:

{plot, evaluations} =  Reap[ Plot3D[(Sow[Thread[{x, y, #}]]; #) &@
           {x, x^2}, {x, 0, 2}, {y, 0, 2}, 
         PlotRange -> All, ImageSize -> Medium]];

plots = MapThread[ListPlot3D[#, ImageSize -> Medium, PlotStyle -> #2] &,
      {Transpose[evaluations[[1]]], {Red, Blue}}];
Row[Prepend[plots, plot]]

enter image description here

mat = {{x^2 + y^2 - 1, 0}, {0, -x^2 - y^2 + 1}};
{plot, evaluations} =  Reap[
   Plot3D[(Sow[Thread[{x, y, #}]]; #) & @ Sort[Eigenvalues[mat]], 
    {x, -2, 2}, {y, -2, 2}, PlotRange -> All, ImageSize -> Medium]];

plots = MapThread[ListPlot3D[#, ImageSize -> Medium, PlotStyle -> #2] &,
      {Transpose[evaluations[[1]]], {Red, Blue}}];
Row[Prepend[plots, plot]]

enter image description here

For

mat2 = {{-1 + x^2 + y^2, 0, x y }, {0, 1 - x^2 - y^2 ,x+y},{x y, x+y,1 - x^2 - y }} ;

(changing {Red, Blue} to {Red, Blue, Green}) we get

enter image description here

$\endgroup$
  • 1
    $\begingroup$ +1 for Rip / Sow $\endgroup$ – Alx Jul 18 at 9:58
  • $\begingroup$ Thank you @Alx. $\endgroup$ – kglr Jul 18 at 9:58
1
$\begingroup$

This is not exactly what you asked for, but it might come in handy regardless. Check out the resource function CheckboxLegended, which makes it easy to create plots with controls for turning data on and off:

https://resources.wolframcloud.com/FunctionRepository/resources/CheckboxLegended

The central idea behind CheckboxLegended is to style your plot with a Dynamic opacity so that you can toggle the opacity between 0 and 1. The following code snippet illustrates the idea:

opacity = 1;
Plot[Style[Sin[x], Opacity[Dynamic[opacity]]], {x, 0, 5}]
Slider[Dynamic[opacity], {0, 1}]
$\endgroup$
  • $\begingroup$ Thanks! I look forward to trying out this new toy! $\endgroup$ – thejacksonjive Jul 18 at 20:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.