Getting coordinates of a level (contour) curve
Hi, I have a fixed bivariate gaussian density function, and I want to compute the coordinates (x(l), y(l)) of a given level curve, where l is the parameter of the curve. I can easily plot the level x curve with matplotlib, using function 'contour', but I have no idea about how to get its coordinates (something like an (N,2) array specifying the coordinates of N points along the curve). With fsolve I can find one of such points, but its is not enough :) Thanks in advance for your help, Best Mico
On Mon, Aug 11, 2008 at 7:14 PM, Mico Filós <elmico.filos@gmail.com> wrote:
Hi,
I have a fixed bivariate gaussian density function, and I want to compute the coordinates (x(l), y(l)) of a given level curve, where l is the parameter of the curve. I can easily plot the level x curve with matplotlib, using function 'contour', but I have no idea about how to get its coordinates (something like an (N,2) array specifying the coordinates of N points along the curve). With fsolve I can find one of such points, but its is not enough :)
For bivariate normal distributions, these equal-density contours are ellipses which you can write down the parametric form of (x,y) from the mean and covariance matrix of your bivariate normal distribution.
Thanks in advance for your help,
Best
Mico
2008/8/11 Bing <bing.jian@gmail.com>:
I have a fixed bivariate gaussian density function, and I want to compute the coordinates (x(l), y(l)) of a given level curve, where l is the parameter of the curve. I can easily plot the level x curve with matplotlib, using function 'contour', but I have no idea about how to get its coordinates (something like an (N,2) array specifying the coordinates of N points along the curve). With fsolve I can find one of such points, but its is not enough :)
For bivariate normal distributions, these equal-density contours are ellipses which you can write down the parametric form of (x,y) from the mean and covariance matrix of your bivariate normal distribution.
For those who are interested in plotting these: http://mentat.za.net/refer/gaussian_intersection.png (See also attached code) Stéfan
Thanks for your quick replies.
For bivariate normal distributions, these equal-density contours are ellipses which you can write down the parametric form of (x,y) from the mean and covariance matrix of your bivariate normal distribution. Yes, you are right. But what if I have a mixture of gaussians, or any other 2D probability density function?
Yes, you are right. But what if I have a mixture of gaussians, or any other 2D probability density function?
Indeed. Isn't the question about how to extract the data points for the curve from the 'contour' object in matplotlib, in the general case? Unfortunately I don't have the answer to that, but maybe introspection of the object would lead to an answer. From the API doc I see a mysterious attribute called 'level'.
On Tue, Aug 12, 2008 at 9:36 AM, Rob Clewley <rob.clewley@gmail.com> wrote:
Yes, you are right. But what if I have a mixture of gaussians, or any other 2D probability density function?
Indeed. Isn't the question about how to extract the data points for the curve from the 'contour' object in matplotlib, in the general case? Unfortunately I don't have the answer to that, but maybe introspection of the object would lead to an answer. From the API doc I see a mysterious attribute called 'level'.
The mpl contour function returns a matplotlib.contour.ContourSet instance which has an attribute "level" array of levels that the contours are drawn on In [57]: CS = plt.contour(X, Y, Z) In [58]: CS.levels Out[58]: array([-1. , -0.5, 0. , 0.5, 1. , 1.5]) It also has an equal length list of line collections (matplotlib.collections.LineCollection) which you can use to extract the x, y vertices of the contour lines at a given level. For a single level, the line collection may contain one or more independent lines in the collections. Here is some example code to get you started: In [59]: level0 = CS.levels[0] In [60]: print level0 -1.0 In [61]: c0 = CS.collections[0] In [62]: paths = c0.get_paths() In [63]: len(paths) Out[63]: 1 In [64]: path0 = paths[0] In [65]: xy = path0.vertices In [66]: xy.shape Out[66]: (237, 2) In [67]: xy[:10,] Out[67]: array([[-0.15 , -0.95150169], [-0.15877627, -0.95 ], [-0.175 , -0.94720234], [-0.2 , -0.94221229], [-0.225 , -0.93652781], [-0.25 , -0.93013814], [-0.26810676, -0.925 ], [-0.275 , -0.9230207 ], [-0.3 , -0.91514134], [-0.325 , -0.90651218]]) In [68]:
Hi all, It's straightforward to *estimate* the level curves of a function that has been evaluated a regularly-spaced grid. (e.g. the "marching cubes" algorithm and it's 2D antecedent, "marching squares".) I suspect that this is what matplotlib is doing. I can send a reasonably-fast C implementation of this for the 2D case if anyone wants. (GPL or BSD, take your pick.) Following a level curve from an arbitrary function is a bit harder. If you have the function's gradient, you could in theory just go around in a direction orthogonal to the gradient, but in the real world that wouldn't work with numerical error and finite step sizes. You could probably take steps orthogonal to the gradient, then correct back to the desired level value by stepping along the gradient, and then repeat until you get back near to where you started. This sounds like far more trouble than it's worth, but if the function is very expensive to evaluate, it might be cheaper and more accurate than evaluating the function on a lattice and then estimating the level curves from that... Zach On Aug 12, 2008, at 10:36 AM, Rob Clewley wrote:
Yes, you are right. But what if I have a mixture of gaussians, or any other 2D probability density function?
Indeed. Isn't the question about how to extract the data points for the curve from the 'contour' object in matplotlib, in the general case? Unfortunately I don't have the answer to that, but maybe introspection of the object would lead to an answer. From the API doc I see a mysterious attribute called 'level'. _______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
Zach,
Following a level curve from an arbitrary function is a bit harder. If you have the function's gradient, you could in theory just go around in a direction orthogonal to the gradient, but in the real world that wouldn't work with numerical error and finite step sizes. You could probably take steps orthogonal to the gradient, then correct back to the desired level value by stepping along the gradient, and then repeat until you get back near to where you started.
You can do that in many circumstances using Newton's method or variants of it. It's called continuation, and it's what packages like AUTO, PyCont, MatCont, and Content do as their bread and butter. There are several steps to getting it set up properly so that it converges. -Rob
Following a level curve from an arbitrary function is a bit harder. If you have the function's gradient, you could in theory just go around in a direction orthogonal to the gradient, but in the real world that wouldn't work with numerical error and finite step sizes. You could probably take steps orthogonal to the gradient, then correct back to the desired level value by stepping along the gradient, and then repeat until you get back near to where you started.
You can do that in many circumstances using Newton's method or variants of it. It's called continuation, and it's what packages like AUTO, PyCont, MatCont, and Content do as their bread and butter. There are several steps to getting it set up properly so that it converges.
Aah, cool... good to know! (I should have mentioned above that I was speculating off the top of my head. Makes sense that there would be formal numerical methodologies for that.) A lot of the PyCont, etc., descriptions seem to be wrapped up in pretty specialized terminology -- what would an example of tracing a level curve look like, given f(x), fprime(x), and some x0? Zach
A lot of the PyCont, etc., descriptions seem to be wrapped up in pretty specialized terminology -- what would an example of tracing a level curve look like, given f(x), fprime(x), and some x0?
I can't really explain it easily off the top of my head. It's a lot like you described, and there are several ways to do it. A popular method is pseudo-arc length continuation, and the idea for it is graphically shown on the wiki page for "numerical continuation". It takes into account how the curve can bend at "fold points," where a naive method based on parameterization of the curve along one axis would lead to a singularity (f' = 0 leading to a 1/0 in the algorithm). You might also read about predictor-corrector methods. However, I might be able to help you set up an example of applying PyCont to finding a level curve - it would be valuable tutorial for the PyCont documentation (I think all our present examples are based on bifurcations in dynamical systems). There are some good book references on the wiki page, maybe the most accessible are [B1], [B5] and [B6] but I don't know a few of them to comment on them all. From the title, [B12] looks like it might be too. -Rob
On Aug 12, 2008, at 2:57 PM, Rob Clewley wrote:
A lot of the PyCont, etc., descriptions seem to be wrapped up in pretty specialized terminology -- what would an example of tracing a level curve look like, given f(x), fprime(x), and some x0?
I can't really explain it easily off the top of my head. It's a lot like you described, and there are several ways to do it. A popular method is pseudo-arc length continuation, and the idea for it is graphically shown on the wiki page for "numerical continuation". It takes into account how the curve can bend at "fold points," where a naive method based on parameterization of the curve along one axis would lead to a singularity (f' = 0 leading to a 1/0 in the algorithm). You might also read about predictor-corrector methods. However, I might be able to help you set up an example of applying PyCont to finding a level curve - it would be valuable tutorial for the PyCont documentation (I think all our present examples are based on bifurcations in dynamical systems). There are some good book references on the wiki page, maybe the most accessible are [B1], [B5] and [B6] but I don't know a few of them to comment on them all. From the title, [B12] looks like it might be too.
Thanks for the background! A basic example on the PyCont page about tracing a level curve given python functions f(x, y), fprime(x, y) (*) and an initial (x0,y0) coordinate would be very useful indeed. I definitely had not realized that PyDSTool/PyCont could be used for these general purposes -- which purposes will be very useful for me and I presume others as well. Zach (*) Or perhaps f_fprime(x, y) which returns the value and gradient? I'm not sure how PyCont is set up...
Attached is a commented example in PyCont for a 2D zero level set that defines an ellipse. It's very easy! I'll add this to PyDSTool's examples. Thanks for giving me the impetus to do this. Let me know what you think. -Rob -- Robert H. Clewley, Ph.D. Assistant Professor Department of Mathematics and Statistics Georgia State University 720 COE, 30 Pryor St Atlanta, GA 30303, USA tel: 404-413-6420 fax: 404-413-6403 http://www2.gsu.edu/~matrhc http://brainsbehavior.gsu.edu/
Aah, cool! Thanks for the example... Zach On Aug 13, 2008, at 1:59 PM, Rob Clewley wrote:
Attached is a commented example in PyCont for a 2D zero level set that defines an ellipse. It's very easy! I'll add this to PyDSTool's examples. Thanks for giving me the impetus to do this. Let me know what you think.
-Rob
-- Robert H. Clewley, Ph.D. Assistant Professor Department of Mathematics and Statistics Georgia State University 720 COE, 30 Pryor St Atlanta, GA 30303, USA
tel: 404-413-6420 fax: 404-413-6403 http://www2.gsu.edu/~matrhc http://brainsbehavior.gsu.edu/ <PyCont_LevelCurve.py>_______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
Hi Rob, I just tried your code after installing PyDSTool from svn. But I get the following error : Velocity around curve is always 1, e.g. look at 100th point norm(Point(sol[100].labels['EP']['data'].V)) =--------------------------------------------------------------------------- NotImplementedError Traceback (most recent call last) /home/cohen/PyCont_LevelCurve.py in <module>() 59 print "\nVelocity around curve is always 1, e.g. look at 100th point" 60 print "norm(Point(sol[100].labels['EP']['data'].V)) =", \ ---> 61 norm(Point(sol[100].labels['EP']['data'].V)) 62 63 print "... at which we have travelled distance ds =", \ /usr/lib/python2.5/site-packages/matplotlib/mlab.pyc in norm(x, y) 1783 Deprecated - see numpy.linalg.norm 1784 """ -> 1785 raise NotImplementedError('Deprecated - see numpy.linalg.norm') 1786 1787 NotImplementedError: Deprecated - see numpy.linalg.norm WARNING: Failure executing file: <PyCont_LevelCurve.py> My version of numpy is '1.2.0.dev5694'. htanks, Johann Rob Clewley wrote:
Attached is a commented example in PyCont for a 2D zero level set that defines an ellipse. It's very easy! I'll add this to PyDSTool's examples. Thanks for giving me the impetus to do this. Let me know what you think.
-Rob
------------------------------------------------------------------------
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
-> 1785 raise NotImplementedError('Deprecated - see numpy.linalg.norm')
Hmm, that's a new one. Well, I guess I'll need to start making sure that the namespace is cleaned up when PyDSTool does its imports (these have been a mixture from numpy, scipy and matplotlib). For now it would appear that you'd get the script to work by importing norm from numpy.linalg explicitly in this script.
thanks, Rob, I confirm that importing explicitely solves this issue. Johann Rob Clewley wrote:
-> 1785 raise NotImplementedError('Deprecated - see numpy.linalg.norm')
Hmm, that's a new one. Well, I guess I'll need to start making sure that the namespace is cleaned up when PyDSTool does its imports (these have been a mixture from numpy, scipy and matplotlib). For now it would appear that you'd get the script to work by importing norm from numpy.linalg explicitly in this script. _______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
participants (7)
-
Bing -
Johann Cohen-Tanugi -
John Hunter -
Mico Filós -
Rob Clewley -
Stéfan van der Walt -
Zachary Pincus