How to vectorise the code below

I have asked about fsurf two or three times on this platform, and received two replies, neither of which extracted me from my difficulties. I have, however, made slow progress. I think I have found one bug in both execution and documentation of fsurf, and whether or not this is really a bug is still under consideration by Mathworks. I have finally succeeded in writing some code that actually outputs an image that seems to be related to what I am aiming at. The code is given the "green light" by the editor, but, when run produces many warnings. These all relate to the complaint that my code does not vectorize. Since the functions concerned are functions of two variables, I am not sure how fsurf would like me to change the code. Here is my code (called SE3):

close all; clear all; %#ok
R = 5; % centre of square is at r=R, z=0, where r is radial coord in (x,y)
%vertices of square in (x,z) plane.
A = [1,1,-1,-1,1;1,-1,-1,1,1];
figure; hold;
for x = 1:4
    funx = @(ang,s) cos(ang)*(R + [cos(ang),sin(ang)]*(s*A(:,x)...
         + (1-s)*A(:,x+1)));
    funy = @(ang,s) sin(ang)*(R + [cos(ang),sin(ang)]*(s*A(:,x)...
         + (1-s)*A(:,x+1)));
    funz = @(ang,s)[-sin(ang),cos(ang)]*(s*A(:,x)+(1-s)*A(:,x+1));
    fsurf(funx,funy,funz,[0,2*pi,0,1]);
end
hold off;

Here is the warning (only the first few lines of many warnings)

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.

> In matlab.graphics.function.ParameterizedFunctionSurface>checkVectorization

In matlab.graphics.function.ParameterizedFunctionSurface/set.XFunction

In matlab.graphics.function.ParameterizedFunctionSurface

In fsurf>singleFsurf (line 263)

In fsurf>@(f1,f2,f3)singleFsurf(cax,{f1,f2,f3},extraOpts,args) (line 236)

In fsurf>vectorizeFsurf (line 236)

In fsurf (line 200)

In SE3 (line 12)

[Rik]

You supressed an m-lint warning, but you shouldn't have: you are clearing a lot of things with clear all, not just variables.

Outside of debugging you should use functions to keep workspaces clean, and while debugging you should use clear variables instead.

As for the other warnings: your functions don't support vectors as inputs, which would probably speed up the process, but not fundamentally change anything in your output. I have no real suggestions for you how you could edit your functions to make them support vector (or array) inputs.

[Stephen23]

"Since the functions concerned are functions of two variables,..."

Vectorization has nothing to do with how many input arguments a function has.

"...I am not sure how fsurf would like me to change the code."

Vectorization is explained in the MATLAB documentation:

https://www.mathworks.com/help/matlab/matlab_prog/vectorization.html

But unfortunately matrix multiplication cannot be easily vectorized, so I don't see any trivial way to vectorize that code. The warning refers to performance, so if that is not your main concern than you can ignore it or suppress it:

https://www.mathworks.com/help/matlab/matlab_prog/suppress-warnings.html

[David Epstein]

@Rik Wisselink: At an earlier stage, I was getting unexplained errors from lint, which I tracked down to the use of symbolic variables. (This was the result of trying to use code from the Mathworks documentatation that I believe is incorrect.) It was indeed necessary to clear the symbolic expressions before running my code. And the inappropriate warnings from m-lint were an irritation in the presence of actual errors I was trying to find.

[David Epstein]

@Stephen Cobeldick: Thanks for the reference to Mathworks documentation on vectorization. Here is an example, pasted from the page you cited:

for n = 1:10000
   V(n) = 1/12*pi*(D(n)^2)*H(n));
end
With MATLAB, you can perform the calculation for each element of a vector with similar syntax as the scalar case:
% Vectorized Calculation
V = 1/12*pi*(D.^2).*H;

Note that this is really an example of a single variable vectorization, namely the variable n. It's unpleasant of fsurf to complain about lack of vectorization in my code, while omitting any mention of vectorization in the documentation.

It seems to me that, in my code, one would have to use meshgrid on my variables ang and s, and this would entail using something like linspace first. Unfortunately, the documentation of fsurf gives no example. I don't think matrix multiplication presents a problem as the matrices involved are 1x2 or 2x1, and so one can use multiplication and addition of scalars. But I will probably have to use an ordinary function, rather than an anonymous one. I'll see if that provokes warnings or errors from fsurf.

[Stephen23]

"Note that this is really an example of a single variable vectorization, namely the variable n."

Not at all, n is just an index. That example has two variables, D and H.

"It's unpleasant of fsurf to complain about lack of vectorization in my code, while omitting any mention of vectorization in the documentation"

Actually the fsurf documentation clearly states: "The function must accept two matrix input arguments and return a matrix output argument of the same size." You ignored this requirement (your functions assume scalar inputs). That requirement could be met in various ways: loops, vectorized code, arrayfun, interpolation, etc. etc.: it is not possible for the fsurf help to know what every user is doing and what is the best way to code their algorithm.

"...one would have to use meshgrid on my variables ang and s, and this would entail using something like linspace first."

Using that approach you might as well just call surf.

[Rik]

I suspect that fsurf uses something like ndgrid internally if it can, so it can have a large speed increase compared to when it has to use loops internally to generate the output.

So, yes, what fsurf is asking for is a function that supports array inputs. As long as your function does not, it doesn't matter if it is an anonymous function, or a regular function.

[David Epstein]

@Stephen Cobeldick: Yes, I did ignore the instruction of have array inputs and outputs. I did this unconsciously, presumably because I couldn't make sense of the instruction. After this discussion, I do now understand what is required.

I had been ruling out the use of surf, because I was assuming it could only be used to display the graph of a real valued function of two variables. I looked again, after seeing your comment, and saw that one can use the syntax surf(X,Y,Z,C). This looks like a much better, and much simpler, approach than the one I was originally adopting, and on which you commented. And I can use matrices and matrix multiplication freely when defining X, Y and Z, which is a great relief.

@Rik Wisselink: Thanks for the further illumination. It's still a puzzle to me how I would use regular functions. (Anonymous functions are no good, because any natural code would require more than one executable.) In my code there are implicitly 12 distinct functions. A factor 4 for x=1:4 and a factor 3 for funx, funy and funz. As I understand Matlab (not very deeply) I cannot define regular functions in a loop. So this would be very ugly, very repetitive, and every improvement in the style of one function definition might be incorporated 12 times---which is very error prone.

I think I can now have a much more successful go at making this image. It has been a somewhat demoralizing experience trying and failing to get this to work, but I think that I finally know what I'm supposed to do.

[Stephen23]

@David Epstein: personally I have never found a use for fsurf and fplot and all those other magical plotting functions. I like to know what my data is doing, which means using meshgrid or ndgrid and generating the values myself. Don't feel demoralized, just chalk it up to experience.

" It's still a puzzle to me how I would use regular functions.... In my code there are implicitly 12 distinct functions. A factor 4 for x=1:4 ..."

One way would be to note that x (and possibly others) are also variables: you don't just have a function of two variables, but a function of three (or four or five) variables, so you can write a function that with those input arguments:

function out = funX(ang,S,x,...)
out = ...
end

and then within the code you can specify those values explicitly or use define a parametrized function, e.g. within your loop:

fx = @(a,s)funX(a,s,x);
...
fsurf(fx, ...)

https://www.mathworks.com/help/matlab/math/parameterizing-functions.html

Thus it would be easy to bring it down to two functions by specifying the X/Y functions and inputting a function handle to distinguish:

function out = funXY(ang,S,x,fun)
out = fun(ang) .* ...;
end

So then you can specify the first factor when the function is parameterized:

fx = @(a,s)funXY(a,s,x,@cos);
fy = @(a,s)funXY(a,s,x,@sin);
...
fsurf(fx, fy, ...)

Answer 1

[Stephen23]

Here are some vectorized versions of your functions:

R = 5;
A = [1,1,-1,-1,1;1,-1,-1,1,1];
B = cat(3,A(1,:),A(2,:));
for x = 1:4
    % Vectorized:
    fz = @(ang,s) sum(cat(3,-sin(ang),+cos(ang)).*(s.*B(1,x,:)+(1-s).*B(1,x+1,:)),3);
    fa = @(ang,s) sum(cat(3,+cos(ang),+sin(ang)).*(s.*B(1,x,:)+(1-s).*B(1,x+1,:)),3);
    fx = @(ang,s) cos(ang).*(R + fa(ang,s));
    fy = @(ang,s) sin(ang).*(R + fa(ang,s));
    % Original:
    funx = @(ang,s) cos(ang)*(R + [cos(ang),sin(ang)]*(s*A(:,x) + (1-s)*A(:,x+1)));
    funy = @(ang,s) sin(ang)*(R + [cos(ang),sin(ang)]*(s*A(:,x) + (1-s)*A(:,x+1)));
    funz = @(ang,s)[-sin(ang),cos(ang)]*(s*A(:,x)+(1-s)*A(:,x+1));
    % Compare:
    new = fx(G,S); % can be called by FSURF.
    old = arrayfun(funx,G,S);
    all(abs(new(:)-old(:))<1e-10)
end