Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

Hello,

I want to visualize a cost function of a localization problem. However if i use fimplicit3 the following warning appears: "Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments" and the plot is empty. Here is my code:

syms x y z
Jd=(-(((x - 10)^2 + (y - 8)^2 + (z - 1/10)^2)^(1/2)/343 - ((x - 9)^2 + (z - 1/10)^2 + y^2)^(1/2)/343 + 166153461756515/288230376151711744)*((1227133513142857*((y - 3)^2 + (z - 1)^2 + x^2)^(1/2))/841813590016 + (8589934592000001*((x - 9)^2 + (z - 1/10)^2 + y^2)^(1/2))/5892695130112 + (500000*((x - 5)^2 + (y - 10)^2 + (z - 3/2)^2)^(1/2))/343 - (1500000*((x - 10)^2 + (y - 8)^2 + (z - 1/10)^2)^(1/2))/343 + 47361719490181303919776257030055/4951760157141521099596496896) - (((y - 3)^2 + (z - 1)^2 + x^2)^(1/2)/343 - ((x - 9)^2 + (z - 1/10)^2 + y^2)^(1/2)/343 + 3755404798969945/288230376151711744)*((500000*((x - 9)^2 + (z - 1/10)^2 + y^2)^(1/2))/343 - (1500000*((y - 3)^2 + (z - 1)^2 + x^2)^(1/2))/343 + (8589934592000001*((x - 5)^2 + (y - 10)^2 + (z - 3/2)^2)^(1/2))/5892695130112 + (1227133513142857*((x - 10)^2 + (y - 8)^2 + (z - 1/10)^2)^(1/2))/841813590016 - 75964017393466287210169242730447/4951760157141521099596496896) - (((x - 5)^2 + (y - 10)^2 + (z - 3/2)^2)^(1/2)/343 - ((x - 9)^2 + (z - 1/10)^2 + y^2)^(1/2)/343 + 564171350756517/72057594037927936)*((8589934592000001*((y - 3)^2 + (z - 1)^2 + x^2)^(1/2))/5892695130112 + (1227133513142857*((x - 9)^2 + (z - 1/10)^2 + y^2)^(1/2))/841813590016 - (1500000*((x - 5)^2 + (y - 10)^2 + (z - 3/2)^2)^(1/2))/343 + (500000*((x - 10)^2 + (y - 8)^2 + (z - 1/10)^2)^(1/2))/343 - 4.9414e+03));
fs=@(x,y,z)double(Jd);
fimplicit3(fs,[-30 30 -30 30 -30 30]);

Due to the fact that the cost function is calculated out of a localization problem it's so huge.

Answer 1

[Walter Roberson]

syms x y z
Jd = (-(((x - 10)^2 + (y - 8)^2 + (z - 1/10)^2)^(1/2)/343 - ((x - 9)^2 + (z - 1/10)^2 + y^2)^(1/2)/343 + 166153461756515/288230376151711744)*((1227133513142857*((y - 3)^2 + (z - 1)^2 + x^2)^(1/2))/841813590016 + (8589934592000001*((x - 9)^2 + (z - 1/10)^2 + y^2)^(1/2))/5892695130112 + (500000*((x - 5)^2 + (y - 10)^2 + (z - 3/2)^2)^(1/2))/343 - (1500000*((x - 10)^2 + (y - 8)^2 + (z - 1/10)^2)^(1/2))/343 + 47361719490181303919776257030055/4951760157141521099596496896) - (((y - 3)^2 + (z - 1)^2 + x^2)^(1/2)/343 - ((x - 9)^2 + (z - 1/10)^2 + y^2)^(1/2)/343 + 3755404798969945/288230376151711744)*((500000*((x - 9)^2 + (z - 1/10)^2 + y^2)^(1/2))/343 - (1500000*((y - 3)^2 + (z - 1)^2 + x^2)^(1/2))/343 + (8589934592000001*((x - 5)^2 + (y - 10)^2 + (z - 3/2)^2)^(1/2))/5892695130112 + (1227133513142857*((x - 10)^2 + (y - 8)^2 + (z - 1/10)^2)^(1/2))/841813590016 - 75964017393466287210169242730447/4951760157141521099596496896) - (((x - 5)^2 + (y - 10)^2 + (z - 3/2)^2)^(1/2)/343 - ((x - 9)^2 + (z - 1/10)^2 + y^2)^(1/2)/343 + 564171350756517/72057594037927936)*((8589934592000001*((y - 3)^2 + (z - 1)^2 + x^2)^(1/2))/5892695130112 + (1227133513142857*((x - 9)^2 + (z - 1/10)^2 + y^2)^(1/2))/841813590016 - (1500000*((x - 5)^2 + (y - 10)^2 + (z - 3/2)^2)^(1/2))/343 + (500000*((x - 10)^2 + (y - 8)^2 + (z - 1/10)^2)^(1/2))/343 - 4.9414e+03));

N = 20;
[X, Y, Z] = ndgrid(linspace(-30, 30, N));
result = double(subs(Jd, {x, y, z}, {X, Y, Z}));
isosurface(X, Y, Z, result, 0)
xlabel('x'); ylabel('y'); zlabel('z');

Unfortunately at the moment, my system is busy so I cannot tell how well the graphic worked. I suspect you will need to increase N to get better resolution.

[Max Lacher]

First of all thank you for your answer. I tried the code but it seems that the isovalue is not correct. If isovalue < 1 the plot is empty. If isovalue > 1 an error occurs, but at least a two dimensional plot is created. The error that occurs is: Input grid is not a valid MESHGRID. I don't really understand what the isovalue is good for in the isosurface function so maybe it's just because I use false values.

[Walter Roberson]

The above did not work (empty plot), so I investigated further.

The minimum value of Jd is approximately -4E-10, over a relatively small region.

It appears to me that the implicit surface is an ellipsoid, that fits approximately in the box

x = 2.89203876360063195 to 2.89204737484551222
y = 5.02378524365726964 to 5.02379993198380337
z = 1.857973242462718   to 1.85810414034411076]

isosurface would generally be able to plot that, but your grid would have to be very fine indeed to have caught this.

And unfortunately, vpasolve(Jd) will not find solutions unless you start very close to that range -- on the order of +/- 1/1000

Fortunately, if you take

fminsearch( matlabFunction(Jd.^2,'vars',{[x,y,z]}), rand(1,3) )

then you get into the right area quickly.

[Max Lacher]

Sorry but I don't understand the proposed solution. fminsearch finds the minimum of the cost function. But how can I plot it now?

Thanks again and I'm sorry for the inconvenience, I don't have that much experience in Matlab.

[Walter Roberson]

When you apply fminsearch to the square of the cost function, then provided the cost function values cannot be imaginary, what you are searching for is zeros of the cost function. This fminsearch is relatively fast, and gives you a point on (near) the surface.

A big problem with the fimplicit3 and the ndgrid approaches is that they depend upon finding very close starting points, since the formula is too complex for areas outside to be of much use in finding the surface.

Once you know a point on (near) the surface, and knowing that you are dealing with an ellipsoid, you should be able to probe around to find the extent of the ellispsoid, such as

temp = 10.^(-10:0);
delta = [-fliplr(temp), 0, temp];

now use delta as offsets to each coordinate in turn, and look for positive (outside), negative(inside), positive (outside) pattern in each dimension. Then pick a grid from that, ndgrid, evaluate, isosurface on value 0

[Max Lacher]

Thank you very much that helps me a lot.

Answer 2

[madhan ravi]

 fs = matlabFunction(Jd)

[Max Lacher]

Ok thanks. Now the warning does not appear anymore, but the figure is still empty.

[Max Lacher]

Or is there another function I can use instead of fimplicit3?

[Walter Roberson]

fimplicit3 takes a long time to compute the faces

[Max Lacher]

I have already waited for a long time and I don't think that the computer is still calculating.