Message Boards Message Boards

0
|
7987 Views
|
10 Replies
|
4 Total Likes
View groups...
Share
Share this post:

[?] Postprocess a ContourPlot to find minimum and maximum values?

Posted 7 years ago

I have a complicated function that is time consuming to evaluate, and has lots of local minima and maxima

cp = ContourPlot[x y Sin[30 x] Sin[20 y], {x, 0, 1}, {y, 0, 1}]

Is there a simple way to find the coordinates and function value of the maximum and minimum of the function, as found by Mathematica in the process of generating the contour plot?

I looked at the FullForm of cp, but it is a long and cryptic GraphicsComplex.

POSTED BY: James D Hanson
10 Replies
Posted 7 years ago

You can use Reap and Sow to collect function values sampled by ContourPlot using it's option EvaluationMonitor. Then you can either select the minimum value or interpolate the points and apply NMinimize.

POSTED BY: Alexey Popkov

If you want to express the contour lines as graphs of interpolating functions, you have to split the contours that are not monotonic in the x values. The following code seems to work in this specific case:

cp = ContourPlot[x^2 + y^2, {x, -1, 1}, {y, -1, 1}]
lns = Cases[Normal[cp[[1]]], _Line, All];
Flatten[lns //. 
    Line[{pts1___List, {x1_, y1_}, {x2_, y2_}, {x3_, y3_}, pts2___} /;
        And[Not[Less[x1, x2, x3]], 
        Not[Less[x3, x2, x1]]]] :> {Line[{pts1, {x1, y1}, {x2, y2}}], 
      Line[{{x2, y2}, {x3, y3}, pts2}]}] /. 
  Line[{{x_, _}, {x_, _}}] :> Nothing /. Line -> Interpolation

It is not clear to me how this could help finding maxima and minima of the 2 variable function.

POSTED BY: Gianluca Gorni

James,

If you really want to parse the ContourPlot it is doable. The example in the thread that was deleted will be of little help because it only had one contour and you have many. @Gianluca Gorni proposed to parse the graphics primitives and extract your data. If you REALLY want to go there, here is what I know and maybe Gianluca has more insight.

If we call the contourplot cp.

ContourPlot produces an object in the form of Graphics[GraphicsComplex[data]]. First step is to take it out of the Graphics with cp[[1]]. GraphicsComplex uses point units so it is of little use as is. The Normal function converts the units to "real" units so when you parse them they will be recognizable.

We now have Normal[cp[[1]]]. The first part of this is polygon colored regions, The second part of this expression are the lines we care about (contour lines).

We now have Normal[cp[[1]]][[2]].

These contour lines have N parts. Each part corresponds with a different contour value so the Nth contour is in
Normal[cp[[1]]][[2,N]]. The first part of this are the lines, the second part is the numerical value of the contour line.

So you need to parse Normal[cp[[1]]][[2,N,1]] to get the list of contour lines and Normal[cp[[1]]][[2,N,2]] gives you the contour value.

This should be enough information to get you started in processing the raw line data. From this you can find the maximum contour and the minimum contour.

In the example you posted I found N = 10 contours in 0.1 increments, several of which were empty. I still maintain that you may be better off generating your data, interpolating it and plotting it later but this approach will work if you really want to extract contour lines.

I hope this helps.

POSTED BY: Neil Singer
Posted 7 years ago

This was the clue I needed. Thank you.

In my actual problem, I had a legended ContourPlot, and the maximum and minimum values found by Mathematica are used for the ColorFunction scaling.

Here is an example. The command below generates a simple ContourPlot with minimum 0.11111 and maximum 0.23456, with a Legend.

cp = ContourPlot[
  a + (b - a) x /. {a -> .11111, b -> .23456}, {x, 0, 1}, {y, 0, 1}, 
  PlotLegends -> Automatic]

Looking at a Shallow version of the FullForm shows that the [[2,1]] Part of cp is the BarLegend.

Shallow[FullForm[cp], 5]

The first argument to BarLegend is a list, the second element of which is a list of the minimum and maximum values for use in the ColorFunction scaling.

cp[[2, 1, 1, 2]]

and the result is

{0.11111, 0.23456}
POSTED BY: James D Hanson

I see @Gianluca Gorni answered your question about postprocessing a contourplot today in another thread -- in case you want to go that route. Was my post helpful? what did you decide to do?

see: http://community.wolfram.com/groups/-/m/t/1128342

Regards

POSTED BY: Neil Singer
Posted 7 years ago

Thank you for your suggestion to first generate the values in a table, and then generate the plot from that. That approach is not quite what I was looking for, because then I lose the benefits of the adaptive evaluation algorithm the Mathematica uses when generating the contour plot.

Your comment from two days is most intriguing. I have been busy, but when I first got your comment, I looked at the link you provided, http://community.wolfram.com/groups/-/m/t/1128342 and I recall that there were two functions applied to the ContourPlot function, I think the first one was Normalize. It looked like the direction I wanted to go in - post-processing the ContourPlot.

Today, when I navigated to that link, I got the message

We're sorry, the thread you were looking for is no longer available. It could have been deleted by the original author or removed by our moderation team.

Seems strange to me.

Do you have any further information about the post-processing approach?

POSTED BY: James D Hanson

James,

You should generate your data first, and then generate your contour plot.

cpdata = Table[
  x y Sin[30 x] Sin[20 y], {y, 0, 1, .01}, {x, 0, 1, .01}]
ListContourPlot[cpdata, DataRange -> {{0, 1}, {0, 1}}]

you get the same contour plot. Note that Contourplot does not find maxima -- it just plots the digitized function. You can get approximate values (the same ones you would get by processing the contourplot data) by processing the cpdata.

You could use Max and Position to find the max value of the digitized data. You can also use Interpolation to create a fast 2D function for Maximize and other data exploration functions which would give you your values. You could even use the interpolation maxima as a starting point for a FindMaximum on your complex function (with really good bounds) to get more precision, if needed.

For example:

In[23]:= Max[cpdata]
Out[23]= 0.853278
In[24]:= Position[cpdata,Max[cpdata]]
Out[24]= {{88,101}}
In[25]:= N[x y Sin[30 x] Sin[20 y]/. x->100/100/.y->87/100]
Out[25]= 0.853278

Or regenerate the data in Interpolation format: (Note -- I made x go a bit larger than 1 because the max we found was at 1)

cpdata2 = Flatten[Table[{{x,y},x y Sin[30 x] Sin[20 y]},{y,0,1,.01},{x,0,1.1,.01}],1];
cdi = Interpolation[cpdata2]

In[30]:= cdi[1, .87]

Out[30]= 0.853278

In[31]:= FindMaximum[cdi[x,y],{{x,1,.99,1.01},{y,.87,.86,.88}}]
Out[31]= {0.861212,{x->0.995971,y->0.866834}}

Hope this helps

Regards,

Neil

POSTED BY: Neil Singer

yes, global optimization is hard

POSTED BY: Frank Kampas
In[1]:= Minimize[{x y Sin[30 x] Sin[20 y], 0 <= x <= 1, 
  0 <= y <= 1}, {x, y}]

Out[1]= {Root[{Sin[30 #1] + 30 Cos[30 #1] #1 &, 
    0.99595288353691310152}] Sin[20] Sin[
   30 Root[{Sin[30 #1] + 30 Cos[30 #1] #1 &, 
      0.99595288353691310152}]], {x -> 
   Root[{Sin[30 #1] + 30 Cos[30 #1] #1 &, 0.99595288353691310152}], 
  y -> 1}}

In[2]:= N[%]

Out[2]= {-0.908742, {x -> 0.995953, y -> 1.}}

Use Maximize to find the maximum value

POSTED BY: Frank Kampas
Posted 7 years ago

The problem with Maximize and NMaximize is that they 1) revaluate the function (that has already been evaluated in ContourPlot), and 2) they are not terribly robust when there are multiple maxima.

Note that the function x y Sin[30 x] Sin[20 y]` that I used above was a proxy for the function I really want to use, which is significantly more complicated, and takes significantly more time to evaluate.

POSTED BY: James D Hanson
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract