ASKSAGE: Sage Q&A Forum - Individual question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Sun, 10 Sep 2017 06:13:33 -0500Empty graph using implicit_plot3d with contour-optionhttp://ask.sagemath.org/question/38741/empty-graph-using-implicit_plot3d-with-contour-option/Hello everybody,
I wanted to make an isosurface plot from a 3d matrix, which contains random values from 0 to 1 at an equally spaced 3d grid.
Therefore I generated random numbers, and reshaped them in a 3d matrix. I interpolated the 3d matrix linearly with the RegularGridInterpolator from scipy. To make a 3D plot of it I am using the implicit_plot3d function of sage with a given contour value.
I get no errors, but in the end the graph is empty, which should not be in my opinion.
Here is my code:
from scipy.interpolate import griddata
import numpy as np
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm
import scipy.interpolate
numbers=np.random.random_sample((1000,)) #generate random numbers
#print numbers
data = np.reshape(numbers, (10, 10, 10)) #reshape the numbers on a 3D matrix
#print data
xi = yi = zi = np.linspace(1, 10, 10)
interp = scipy.interpolate.RegularGridInterpolator((xi,yi,zi), data) #interpolate the 3d matrix with a function
#print xi
var('x,y,z')
#test = implicit_plot3d(interp,(x,1,10),(y,1,10),(z,1,10),contour=0.5)
test = implicit_plot3d(interp==0.5,(x,1,10), (y,1,10),(z,1,10),plot_points=60, color='seagreen')
test.show() #plot the the function at a certain value
#a=(5.,5.,2.)
#interp(a)
Any ideas on that?
Thank you very much!Tue, 05 Sep 2017 04:31:34 -0500http://ask.sagemath.org/question/38741/empty-graph-using-implicit_plot3d-with-contour-option/Answer by dan_fulea for <p>Hello everybody,</p>
<p>I wanted to make an isosurface plot from a 3d matrix, which contains random values from 0 to 1 at an equally spaced 3d grid.</p>
<p>Therefore I generated random numbers, and reshaped them in a 3d matrix. I interpolated the 3d matrix linearly with the RegularGridInterpolator from scipy. To make a 3D plot of it I am using the implicit_plot3d function of sage with a given contour value.</p>
<p>I get no errors, but in the end the graph is empty, which should not be in my opinion.</p>
<p>Here is my code:</p>
<pre><code>from scipy.interpolate import griddata
import numpy as np
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm
import scipy.interpolate
numbers=np.random.random_sample((1000,)) #generate random numbers
#print numbers
data = np.reshape(numbers, (10, 10, 10)) #reshape the numbers on a 3D matrix
#print data
xi = yi = zi = np.linspace(1, 10, 10)
interp = scipy.interpolate.RegularGridInterpolator((xi,yi,zi), data) #interpolate the 3d matrix with a function
#print xi
var('x,y,z')
#test = implicit_plot3d(interp,(x,1,10),(y,1,10),(z,1,10),contour=0.5)
test = implicit_plot3d(interp==0.5,(x,1,10), (y,1,10),(z,1,10),plot_points=60, color='seagreen')
test.show() #plot the the function at a certain value
#a=(5.,5.,2.)
#interp(a)
</code></pre>
<p>Any ideas on that?</p>
<p>Thank you very much!</p>
http://ask.sagemath.org/question/38741/empty-graph-using-implicit_plot3d-with-contour-option/?answer=38749#post-id-38749The following works for me:
import numpy
import scipy.interpolate
n = 10
numbers = numpy.random.random_sample( ( n**3, ) )
data = numpy.reshape( numbers, (n,n,n) )
xi = yi = zi = np.linspace(1, n, n)
interp = scipy.interpolate.RegularGridInterpolator((xi,yi,zi), data)
var('x,y,z')
test = implicit_plot3d( lambda x,y,z: interp( (x,y,z) )
, (x,1,n)
, (y,1,n)
, (z,1,n)
, plot_points=60, color='seagreen', contour=0.5)
test.show()
My favorite values for `n` were in the tests $3$ or $5$, to make the machine less angry. And for the numbers, it was `numbers = [ 1 for _ in range(n**3) ]` .
The difference is the use of `lambda x,y,z: interp( (x,y,z) )` instead of `interp`.Thu, 07 Sep 2017 07:36:03 -0500http://ask.sagemath.org/question/38741/empty-graph-using-implicit_plot3d-with-contour-option/?answer=38749#post-id-38749Comment by mforets for <p>The following works for me:</p>
<pre><code>import numpy
import scipy.interpolate
n = 10
numbers = numpy.random.random_sample( ( n**3, ) )
data = numpy.reshape( numbers, (n,n,n) )
xi = yi = zi = np.linspace(1, n, n)
interp = scipy.interpolate.RegularGridInterpolator((xi,yi,zi), data)
var('x,y,z')
test = implicit_plot3d( lambda x,y,z: interp( (x,y,z) )
, (x,1,n)
, (y,1,n)
, (z,1,n)
, plot_points=60, color='seagreen', contour=0.5)
test.show()
</code></pre>
<p>My favorite values for <code>n</code> were in the tests $3$ or $5$, to make the machine less angry. And for the numbers, it was <code>numbers = [ 1 for _ in range(n**3) ]</code> .</p>
<p>The difference is the use of <code>lambda x,y,z: interp( (x,y,z) )</code> instead of <code>interp</code>.</p>
http://ask.sagemath.org/question/38741/empty-graph-using-implicit_plot3d-with-contour-option/?comment=38780#post-id-38780this is probably slow not of `implicit_plot3d` itself, but because it goes through the Python function `lambda x,y,z: interp( (x,y,z) )` (which is not statically typed). the plot data could be efficiently generated by a single call to `interp`, which accepts array arguments, eg. try `interp(np.array([[1.0, 2.0, 3.0], [3.0, 5.0, 6.0]]))`. but then i don't know how to interface with Sage's `implicit_plot3d`.Sun, 10 Sep 2017 06:13:33 -0500http://ask.sagemath.org/question/38741/empty-graph-using-implicit_plot3d-with-contour-option/?comment=38780#post-id-38780Comment by os for <p>The following works for me:</p>
<pre><code>import numpy
import scipy.interpolate
n = 10
numbers = numpy.random.random_sample( ( n**3, ) )
data = numpy.reshape( numbers, (n,n,n) )
xi = yi = zi = np.linspace(1, n, n)
interp = scipy.interpolate.RegularGridInterpolator((xi,yi,zi), data)
var('x,y,z')
test = implicit_plot3d( lambda x,y,z: interp( (x,y,z) )
, (x,1,n)
, (y,1,n)
, (z,1,n)
, plot_points=60, color='seagreen', contour=0.5)
test.show()
</code></pre>
<p>My favorite values for <code>n</code> were in the tests $3$ or $5$, to make the machine less angry. And for the numbers, it was <code>numbers = [ 1 for _ in range(n**3) ]</code> .</p>
<p>The difference is the use of <code>lambda x,y,z: interp( (x,y,z) )</code> instead of <code>interp</code>.</p>
http://ask.sagemath.org/question/38741/empty-graph-using-implicit_plot3d-with-contour-option/?comment=38772#post-id-38772Thank you for your ideas!
I do not think it is because of the interpolation; commenting out the lines
test =...
test.show()
and starting sage takes about 2 seconds for the whole process (n=2).
The trick with the brackets [()] doesn't seem to have any impact, for n=2 with and without brackets sage needs about 2 minutes on my computer.
Later I will try the matplotlib contour3D function (not sure if they are suited to this problem), maybe they are faster... (http://matplotlib.org/mpl_toolkits/mplot3d/api.html?highlight=contour3d#mpl_toolkits.mplot3d.axes3d.Axes3D.contour3D)Sat, 09 Sep 2017 08:45:33 -0500http://ask.sagemath.org/question/38741/empty-graph-using-implicit_plot3d-with-contour-option/?comment=38772#post-id-38772Comment by dan_fulea for <p>The following works for me:</p>
<pre><code>import numpy
import scipy.interpolate
n = 10
numbers = numpy.random.random_sample( ( n**3, ) )
data = numpy.reshape( numbers, (n,n,n) )
xi = yi = zi = np.linspace(1, n, n)
interp = scipy.interpolate.RegularGridInterpolator((xi,yi,zi), data)
var('x,y,z')
test = implicit_plot3d( lambda x,y,z: interp( (x,y,z) )
, (x,1,n)
, (y,1,n)
, (z,1,n)
, plot_points=60, color='seagreen', contour=0.5)
test.show()
</code></pre>
<p>My favorite values for <code>n</code> were in the tests $3$ or $5$, to make the machine less angry. And for the numbers, it was <code>numbers = [ 1 for _ in range(n**3) ]</code> .</p>
<p>The difference is the use of <code>lambda x,y,z: interp( (x,y,z) )</code> instead of <code>interp</code>.</p>
http://ask.sagemath.org/question/38741/empty-graph-using-implicit_plot3d-with-contour-option/?comment=38770#post-id-38770Yes, thanks, that `np.linspace(1, n, n)` should be `numpy.linspace(1, n, n)` . (Too much copy+paste.)
For the performance... maybe the reason is the "complicated" type of the result of interpolation. With the above notations, we have for instance (this time):
sage: v = interp( (1,1,1) )
sage: v
array(0.4639585570135113)
and it is not obvious to make a number out of `v`. The trick is to call "v of nothing", `v[ () ]` to get:
sage: v[()]
0.4639585570135113
and the implicit plot also has to find this or his way. So i tried
test = implicit_plot3d( lambda x,y,z: interp( (x,y,z) )[()]
, (x,1,n), (y,1,n), (z,1,n)
, plot_points=60, color='seagreen', contour=0.5 )
Or maybe the interpolation kills the performance.Sat, 09 Sep 2017 07:26:40 -0500http://ask.sagemath.org/question/38741/empty-graph-using-implicit_plot3d-with-contour-option/?comment=38770#post-id-38770Comment by os for <p>The following works for me:</p>
<pre><code>import numpy
import scipy.interpolate
n = 10
numbers = numpy.random.random_sample( ( n**3, ) )
data = numpy.reshape( numbers, (n,n,n) )
xi = yi = zi = np.linspace(1, n, n)
interp = scipy.interpolate.RegularGridInterpolator((xi,yi,zi), data)
var('x,y,z')
test = implicit_plot3d( lambda x,y,z: interp( (x,y,z) )
, (x,1,n)
, (y,1,n)
, (z,1,n)
, plot_points=60, color='seagreen', contour=0.5)
test.show()
</code></pre>
<p>My favorite values for <code>n</code> were in the tests $3$ or $5$, to make the machine less angry. And for the numbers, it was <code>numbers = [ 1 for _ in range(n**3) ]</code> .</p>
<p>The difference is the use of <code>lambda x,y,z: interp( (x,y,z) )</code> instead of <code>interp</code>.</p>
http://ask.sagemath.org/question/38741/empty-graph-using-implicit_plot3d-with-contour-option/?comment=38768#post-id-38768Thank you very much for your answer, I wasn't familiar with lambda functions. However, I am still a little bit confused, why it takes so long to make the implicit_plot3d.
In mathematica the in my opinion equivalent function would be ContourPlot3D (http://reference.wolfram.com/language/ref/ContourPlot3D.html). I already used it to plot interpolated functions to 50x50x50 matrices and it seems to be considerably faster, or in other words with implicit_plot3D it doesn't seem to be possible to do that for n=50.
Any ideas on that?
Ps. I think it should be numpy.linspace instead of np.linspace (line 7)Sat, 09 Sep 2017 02:34:13 -0500http://ask.sagemath.org/question/38741/empty-graph-using-implicit_plot3d-with-contour-option/?comment=38768#post-id-38768