# Conditional styling of contour

My last question on varying plot styles was perhaps too minimal of an example (and a duplicate to boot!). Here's a more realistic example that still gives me trouble.

I'd like to plot an isocline where one function (dg[x1,x2]) equals zero and style it differently depending on the sign of another function (dg2[x1,x2]). For example,

dg[x1_, x2_] := (2 (x1 + E^(4 (x1 - x2)^2) x1 - 2 x1^3 - 2 x2 + 2 x1^2 x2 + 2 E^(2 (x1 - x2)^2) (x1 - x2) (-1 + x2^2)))/((-1 + E^(4 (x1 - x2)^2)) (-1 + x1^2));
dg2[x1_, x2_] := 2/((-1 + E^(4 (x1 - x2)^2)) (-1 + x1^2)^2) (-1 + 7 x1^2 - 10 x1^4 + 8 x1^6 + E^(4 (x1 - x2)^2) (1 - 7 x1^2 + 2 x1^4) - 8 x1 x2 + 24 x1^3 x2 - 16 x1^5 x2 + 8 x2^2 - 16 x1^2 x2^2 + 8 x1^4 x2^2 - 8 E^(2 (x1 - x2)^2) (x1 - x2) (-1 + x2) (1 + x2) (x1^3 + x2 - x1^2 x2));


Here's what the isocline looks like and the regions that should be styled differently:

iso = ContourPlot[dg[x1, x2] == 0, {x1, -1, 1}, {x2, -1, 1},
Exclusions -> {x1 == x2}, PlotPoints -> 50];
cp = ContourPlot[dg2[x1, x2], {x1, -1, 1}, {x2, -1, 1},
PlotPoints -> 50, Contours -> {0}, ContourShading -> {White, LightRed}];
Show[cp, iso] My approach so far is to extract the contours, then ListLinePlot and color them using ColorFunction:

ExtractPlotPoints[plot_Graphics] := Cases[Normal@plot, Line[x_] :> x, \[Infinity]];

isocol = ListLinePlot[ExtractPlotPoints[iso],
ColorFunction -> (If[dg2[#1, #2] > 0, Red, Black] &),
ColorFunctionScaling -> False, Frame -> True, Axes -> False,
AspectRatio -> 1, PlotRange -> {{-1, 1}, {-1, 1}}];
Show[cp, isocol] This works great, but I'd like other style options besides varying the color (say, thick vs thin). This is where I'm stuck. I tried MeshFunctions but got nothing:

ListLinePlot[ExtractPlotPoints[iso],
MeshFunctions -> {dg2[#1, #2] &}, Mesh -> {{0}},
MeshShading -> {Directive[Black, Thick], Directive[Red, Thin]},
MeshStyle -> None, Frame -> True, Axes -> False, AspectRatio -> 1] dg2 is sketchy when x1==x2 due to the denominator equalling zero, but it has a well-defined limit.

Limit[dg2[x1, x2], {x2 -> x1}]
(* (-2 + 4 x1^2)/(-1 + x1^2) *)


Defining a better dg2 when x1==x2 helps, but isn't quite right:

dg2[x1_, x1_] = (-2 + 4 x1^2)/(-1 + x1^2);
isonew = ListLinePlot[ExtractPlotPoints[iso],
MeshFunctions -> {dg2[#1, #2] &}, Mesh -> {{0}},
MeshShading -> {Directive[Black, Thick], Directive[Red, Thin]},
MeshStyle -> None, Frame -> True, Axes -> False, AspectRatio -> 1];
Show[cp, isonew] Any thoughts? The functions dg and dg2 may differ from this example.

(The eventual application is isoclines of evolutionary dynamics as in Fig. 6a of this paper.)

cp2 = Quiet@ ContourPlot[{ConditionalExpression[dg[x1, x2], dg2[x1, x2] < 0] ==  0,
ConditionalExpression[dg[x1, x2], dg2[x1, x2] > 0] ==  0},
{x1, -1, 1}, {x2, -1, 1}, PlotPoints -> 50,  Exclusions -> {x1 == x2},
ContourStyle -> {Directive[Black, Thick, Dashed], Directive[Red, Thin]}];

Show[cp, cp2] Alternatively,

dg2b[x1_, x2_] :=  If[x1 == x2, (-2 + 4 x1^2)/(-1 + x1^2),
2/((-1 + E^(4 (x1 - x2)^2)) (-1 + x1^2)^2) (-1 + 7 x1^2 -
10 x1^4 + 8 x1^6 + E^(4 (x1 - x2)^2) (1 - 7 x1^2 + 2 x1^4) -
8 x1 x2 + 24 x1^3 x2 - 16 x1^5 x2 + 8 x2^2 - 16 x1^2 x2^2 +
8 x1^4 x2^2 - 8 E^(2 (x1 - x2)^2) (x1 - x2) (-1 + x2) (1 + x2) (x1^3 + x2 - x1^2 x2))];

isonewb = ListLinePlot[ExtractPlotPoints[iso],
MeshFunctions -> {dg2b[#1, #2] &}, Mesh -> {{0}},
MeshShading -> {Directive[Black, Thick], Directive[Red, Thin]},
MeshStyle -> None, Frame -> True, Axes -> False, AspectRatio -> 1];
Show[cp, isonewb] • Thanks, that's elegant. Unfortunately for me, it's significantly slower (4.33s vs 1.2s). Of course I can spare 3.13s in this example, but sometimes dg and dg2 are a lot more computationally expensive, so I'd like to save that 3.6X speedup if possible. But maybe it's an inevitable trade-off. – Chris K Feb 21 '18 at 4:28
• Huh, when I run your edited version, I get the same as my last figure where the parts don't line up (M11.2 on MacOS10.13.3). ¯\_(ツ)_/¯ – Chris K Feb 22 '18 at 3:32
• @ChrisK, version 9.0 Windows 10 here. – kglr Feb 22 '18 at 3:54

My colleague suggested the following variation, which is both fast and accurate (avoids potential problems I had with MeshFunctions). Idea: make ContourPlot, extract contours, make a separate ListLinePlot for each desired style using Opacity in ColorFunction to blank out the undesired parts, then combine.

iso1 = ListLinePlot[ExtractPlotPoints[iso], ColorFunctionScaling -> False,
ColorFunction -> (If[dg2[#1, #2] > 0, {Red}, {Opacity}] &),
PlotStyle -> Thin];
iso2 = ListLinePlot[ExtractPlotPoints[iso], ColorFunctionScaling -> False,
ColorFunction -> (If[dg2[#1, #2] < 0, {Black}, {Opacity}] &),
PlotStyle -> Thick];
Show[cp, iso1, iso2, Frame -> True, Axes -> False, AspectRatio -> 1,
PlotRange -> {{-1, 1}, {-1, 1}}] 