Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Ways to skip evaluation of function for specific coordinates?

Subject: Ways to skip evaluation of function for specific coordinates?

From: AJP

Date: 5 Jan, 2011 15:49:05

Message: 1 of 12

I am evaluating a velocity potential function on a polar grid (theta, r) produced from a cart2pol(X,Y) where [X,Y]=meshgrid(x,y) of two vectors x and y=linspace(min,max,n). The grids used in computation are therefore n x n.

Since the velocity potential function has no useful value for (theta, r) coordinates where r<R (R = the radius of a cylindrical structure) and the grid is very fine, I want to skip evaluation of the function at these "useless" coordinates.

To try and achieve this I set r(r<R) = NaN in the hope MATLAB would skip computation at these locations. However, the computations are taking the same time as before, so MATLAB is evidently still trying to evaluate the NaNs.

I then tried r(r<R)=[] which was just as slow and caused problems when plotting the results.

Can suggest a way to either define to grid in a better way so I don't create grid points with r<R or some way to skip evaluation of the function for these grid points?

I still want to retain the ability to plot the results from the function, evaluated using the polar coordinates, on the X,Y cartesian grid.

Subject: Ways to skip evaluation of function for specific coordinates?

From: Walter Roberson

Date: 5 Jan, 2011 16:13:41

Message: 2 of 12

On 05/01/11 9:49 AM, AJP wrote:
> I am evaluating a velocity potential function on a polar grid (theta, r)
> produced from a cart2pol(X,Y) where [X,Y]=meshgrid(x,y) of two vectors x
> and y=linspace(min,max,n). The grids used in computation are therefore n
> x n.
>
> Since the velocity potential function has no useful value for (theta, r)
> coordinates where r<R (R = the radius of a cylindrical structure) and
> the grid is very fine, I want to skip evaluation of the function at
> these "useless" coordinates.

Use logical indexing and only apply the function to the "useful" positions.

You might perhaps find one of the recent additions to the FAQ to be
useful as a guide to constructing the appropriate code:

http://matlab.wikia.com/wiki/FAQ#How_do_I_create_a_circle.3F

Subject: Ways to skip evaluation of function for specific coordinates?

From: Matt J

Date: 5 Jan, 2011 16:27:05

Message: 3 of 12

"AJP" wrote in message <ig23th$n2b$1@fred.mathworks.com>...
>
> I then tried r(r<R)=[] which was just as slow and caused problems when plotting the results.
=====

And you also did theta(r<R)=[] presumably?

If eliminating points did not help, you either have too few points with r<R for speed to be affected significantly, or else there are bottlenecks in your function that make its execution time non-linear in length(r).

We would have to see the function you are passing this to to get a better idea.

Subject: Ways to skip evaluation of function for specific coordinates?

From: AJP

Date: 5 Jan, 2011 17:11:21

Message: 4 of 12

"Matt J" wrote in message <ig264p$j99$1@fred.mathworks.com>...
> And you also did theta(r<R)=[] presumably?
>
> If eliminating points did not help, you either have too few points with r<R for speed to be affected significantly, or else there are bottlenecks in your function that make its execution time non-linear in length(r).
>
> We would have to see the function you are passing this to to get a better idea.


Yes I also put in theta(r<R)=[]. The grid is large and the function is time-consuming to evaluate, but is written reasonably efficiently. I am sure it is the number of grid points, not the function, that is causing the slow perfomance.

The problem with setting r(r<R)=[] and similarly for theta this is that it turns my n x n (theta, r) grid into two row vectors that are 1 x (n^2 - numel(r(r<R))).

When I go to plot these there's trouble as the plot commands (countourf, surf etc.) are expecting arrays.

My friend tells me that I can use griddata or triscatteredinterp to recreate the arrays, using a condition to NaN results that are on r<R. However, since I want to plot on a cartesian grid using my original X,Y arrays, and gridata or triscattered interp need to interpolate a polar funtion, I am not sure how this works.

Subject: Ways to skip evaluation of function for specific coordinates?

From: Matt J

Date: 5 Jan, 2011 17:43:05

Message: 5 of 12

"AJP" wrote in message <ig28np$9r7$1@fred.mathworks.com>...
>
> Yes I also put in theta(r<R)=[]. The grid is large and the function is time-consuming to evaluate, but is written reasonably efficiently. I am sure it is the number of grid points, not the function, that is causing the slow perfomance.
================

If so, then why even bother trying to eliminate the points r<R? There are clearly not enough of them to make a significant impact on the computation, or at least you said there was no impact in your original post.


> The problem with setting r(r<R)=[] and similarly for theta this is that it turns my n x n (theta, r) grid into two row vectors that are 1 x (n^2 - numel(r(r<R))).
=====

That's a minor problem. Similar to what Walter was saying, you can extract the values you need for processing and then afterward re-embed the processed values in the original array for plotting. It just involves logical indexing manipulations like the following:

idx=(r>=R); %points to process

[r_new,theta_new]=ProcessingFunction( r(idx), theta(idx));

r(idx)=r_new;
theta(idx)=theta_new;

Subject: Ways to skip evaluation of function for specific coordinates?

From: Matt J

Date: 5 Jan, 2011 18:00:23

Message: 6 of 12

"Matt J" wrote in message <ig2aj9$b4k$1@fred.mathworks.com>...
>
> That's a minor problem. Similar to what Walter was saying, you can extract the values you need for processing and then afterward re-embed the processed values in the original array for plotting. It just involves logical indexing manipulations like the following:
>
> idx=(r>=R); %points to process
>
> [r_new,theta_new]=ProcessingFunction( r(idx), theta(idx));
>
> r(idx)=r_new;
> theta(idx)=theta_new;
=========

This might be a better description:

 idx=(r>=R); %points to process

 V=nan(size(r)); %stores the velocity potentials

V(idx)=VelocityPotential( r(idx), theta(idx)); %evaluate

Subject: Ways to skip evaluation of function for specific coordinates?

From: AJP

Date: 6 Jan, 2011 08:53:04

Message: 7 of 12

"AJP" wrote in message <ig28np$9r7$1@fred.mathworks.com>...
>
> My friend tells me that I can use griddata or triscatteredinterp to recreate the arrays, using a condition to NaN results that are on r<R. However, since I want to plot on a cartesian grid using my original X,Y arrays, and gridata or triscattered interp need to interpolate a polar funtion, I am not sure how this works.

==================================================

Right, I've cracked it using a combination of the suggestions above and those of a friend on a separate messageboard. Here's the process for reference:

First I create 2 linspace x and y vectors, each length n, then meshgrid them to get n x n matrices X,Y. I then cart2pol(X,Y) to get matrices (theta,r). In order to get varying grid resolutions close to and far from the cylinder I could also define separate linspaces and join them together or use logspaces.

Then I get ind = r<R and set theta(ind) = [] and r(ind)=[] giving me two "truncated" polar coordinate row vectors theta and r that are 1 x (n^2 - numel(ind)). Since numel(ind) is fairly large due to small grid spacings (i.e. large n) and R being large relative to the grid extent, this step eliminates a large number of points, all of which are, as noted above, "useless" for the function evaluation.

I also create two correpsonding X1,Y1 coordinate row vectors where X1=X(ind) and Y1=Y(ind). This allows me to keep equivalent truncated polar and cartesian coordinate systems ready for later.

I then evaluate the velocity potential function VP = fVP(theta,r) on the truncated polar coordinate row vectors, to produce a corresponding row vector of VP values.

Since I want to plot the VP values on the original meshgrid matrices X,Y the row vector must be made back into a matrix and the "hole" made by the truncation operations above filled with NaN.

To do this I use F = TriScatteredInterp(X1',Y1',VP') which generates an interpolant function, F, of the column vector data VP', corresponding to coordinate column vectors X1' and Y1', assuming VP=f(X1,Y1).

F can then be evaluated on X,Y to produce Z = F(X,Y), where Z is an n x n matrix of VP values interpolated from the scattered set (X1,Y1,VP).

F will produce "dummy" interpolated values the fill the "hole" left by the coordinate truncation. Setting Z(ind)=NaN removes these "dummy" interpolated results from the matrix.

I then simply contourf(X,Y,Z) to finish.

This process still isn't ideal as for maximum efficiency it would be best to increase the radial stepsize with increasing r, not the cartesian stepsize with increasing x or y. However, since I want to end up with a rectangular cartesian plot, it's probably easiest to create the X,Y grid first then the theta,r grid from this using cart2pol as above.

I had never used TriScatteredInterp before, but it works very well and is extremely fast in this application. Thus main objective, to speed up the evaluation of VP by skipping useless gridpoints, has been achieved.

Thanks to all for your inspiring suggestions!

Subject: Ways to skip evaluation of function for specific coordinates?

From: Matt J

Date: 6 Jan, 2011 11:55:21

Message: 8 of 12

"AJP" wrote in message <ig3vtg$gfj$1@fred.mathworks.com>...
> "AJP" wrote in message <ig28np$9r7$1@fred.mathworks.com>...
> >
> > My friend tells me that I can use griddata or triscatteredinterp to recreate the arrays, using a condition to NaN results that are on r<R. However, since I want to plot on a cartesian grid using my original X,Y arrays, and gridata or triscattered interp need to interpolate a polar funtion, I am not sure how this works.
>
> ==================================================
>
> Right, I've cracked it using a combination of the suggestions above and those of a friend on a separate messageboard. Here's the process for reference:
>
> First I create 2 linspace x and y vectors, each length n, then meshgrid them to get n x n matrices X,Y. I then cart2pol(X,Y) to get matrices (theta,r). In order to get varying grid resolutions close to and far from the cylinder I could also define separate linspaces and join them together or use logspaces.
>
> Then I get ind = r<R and set theta(ind) = [] and r(ind)=[] giving me two "truncated" polar coordinate row vectors theta and r that are 1 x (n^2 - numel(ind)). Since numel(ind) is fairly large due to small grid spacings (i.e. large n) and R being large relative to the grid extent, this step eliminates a large number of points, all of which are, as noted above, "useless" for the function evaluation.
>
> I also create two correpsonding X1,Y1 coordinate row vectors where X1=X(ind) and Y1=Y(ind). This allows me to keep equivalent truncated polar and cartesian coordinate systems ready for later.
>
> I then evaluate the velocity potential function VP = fVP(theta,r) on the truncated polar coordinate row vectors, to produce a corresponding row vector of VP values.
>
> Since I want to plot the VP values on the original meshgrid matrices X,Y the row vector must be made back into a matrix and the "hole" made by the truncation operations above filled with NaN.
>
> To do this I use F = TriScatteredInterp(X1',Y1',VP') which generates an interpolant function, F, of the column vector data VP', corresponding to coordinate column vectors X1' and Y1', assuming VP=f(X1,Y1).
>
> F can then be evaluated on X,Y to produce Z = F(X,Y), where Z is an n x n matrix of VP values interpolated from the scattered set (X1,Y1,VP).
>
> F will produce "dummy" interpolated values the fill the "hole" left by the coordinate truncation. Setting Z(ind)=NaN removes these "dummy" interpolated results from the matrix.
===========

Well, if you're happy, I guess that's the important thing, and I suppose it's also possible that I still don't understand something about your situation.

However, it seems really absurd to use TriScatteredInterp purely for the purpose of creating "dummy" values that you then remove later with Z(ind)=NaN. As we've been arguing, you could have just done

Z=nan(size(X));
Z(~ind)=VP(:);

and skipped the considerable computational effort used in TriScatteredInterp.

Subject: Ways to skip evaluation of function for specific coordinates?

From: AJP

Date: 6 Jan, 2011 12:09:04

Message: 9 of 12

"Matt J" wrote in message <ig4aj9$edn$1@fred.mathworks.com>...
>
> However, it seems really absurd to use TriScatteredInterp purely for the purpose of creating "dummy" values that you then remove later with Z(ind)=NaN. As we've been arguing, you could have just done
>
> Z=nan(size(X));
> Z(~ind)=VP(:);
>
> and skipped the considerable computational effort used in TriScatteredInterp.

=================================================

VP was calculated from the truncated polar row vectors, and is hence a row vector, whereas Z and ind are n x n matrices. Therefore wouldn't the operation Z(~ind)=VP(:) fail as the dimensions of these arrays do not match?

Subject: Ways to skip evaluation of function for specific coordinates?

From: AJP

Date: 6 Jan, 2011 12:22:05

Message: 10 of 12

"AJP" wrote in message <ig4bd0$6up$1@fred.mathworks.com>...
>
> VP was calculated from the truncated polar row vectors, and is hence a row vector, whereas Z and ind are n x n matrices. Therefore wouldn't the operation Z(~ind)=VP(:) fail as the dimensions of these arrays do not match?

=============================================

Yes I have tried your suggestion and indeed, while Z=nan(size(X)) produces an n x n matrix, the operation Z(~ind)=VP(:) changes it to a 1 x (n^2 - numel(ind)) row vector.

This is because VP(:) is a row vector also.

The triscatteredinterp calculations are required to create a plot on a uniform grid of a surface sampled at scattered points. This is what I have because of the "hole" in the middle of the grid created due to eliminating these points before calculating VP.

You are right in that triscatteredinterp unfortunately produces "dummy" values within the region r<R, and these must be set to NaN. But the computation of the dummy values is by interpolation of the surrounding data points, which is much faster than evaluating VP within R, as was my original problem.

Subject: Ways to skip evaluation of function for specific coordinates?

From: AJP

Date: 6 Jan, 2011 12:59:08

Message: 11 of 12

"AJP" wrote in message <ig4c5d$pks$1@fred.mathworks.com>...
> "AJP" wrote in message <ig4bd0$6up$1@fred.mathworks.com>...

> =============================================
>
> Yes I have tried your suggestion and indeed, while Z=nan(size(X)) produces an n x n matrix, the operation Z(~ind)=VP(:) changes it to a 1 x (n^2 - numel(ind)) row vector.
>

MattJ,

I have ssen the error of my ways.

You are indeed right that you can recover the original n x n form of the truncated row vectors using logial indexing, as you suggested above, but that I had implemented incorrectly so believed it not to work.

Thank you, now I can avoid the use of triscatteredinterp.

Thanks for your perservering comments, otherwise I would not have reached this efficient solution!

Adam

Subject: Ways to skip evaluation of function for specific coordinates?

From: Matt J

Date: 6 Jan, 2011 14:53:05

Message: 12 of 12

"AJP" wrote in message <ig4eas$fpp$1@fred.mathworks.com>...
>
> You are indeed right that you can recover the original n x n form of the truncated row vectors using logial indexing, as you suggested above, but that I had implemented incorrectly so believed it not to work.
>
> Thank you, now I can avoid the use of triscatteredinterp.
=====

And in fact, you should be able to just do

Z(~ind)=VP;

without the '(:)'. When you do a logically indexed assignment, the shape of the right hand side doesn't matter. It only matters that it have the correct number of elements, e.g.,

>> A=[1 3 5; 2 4 6]

A =

     1 3 5
     2 4 6

 
>> A(A>2)=[20 30 40 50]

A =

     1 20 40
     2 30 50


>> A(A>2)=[200 300 ;400 500]

A =

     1 200 300
     2 400 500

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us