Thread Subject: Mapping a random walk to another spatial grid

Hello all
I calculate a number of random walks, generating x and y coordinates. I wish to map these to a spatial grid (X,Y) of a different resolution. What I want to do is count the number of 'stops' (x,y pairs) that 'land' in each cell of the spatial grid (X,Y). The code below might make it a bit easier to see what I mean. There is a severe bottleneck in the last part of the code where i calculate euclidean distances. If I increase the size and number of the walks (which I really need to do, by a great deal), then this part of the code slows everything down to a crawl. Can anyone think of a way to speed that part up, or replace it?

--------------------------------------------------------------------------

n = 100; % number of steps in each walk
N = 100; % number of separate walks
minstep = 10; % minimum length of each step
mu = 2; % exponent of a power law distribution
x = zeros(N,n); % x coordinates of the walks
y = zeros(N,n); % y coordinates of the walks

angle = 2*pi*rand(N,n); % random trajectory angles
length = minstep*rand(N,n).^(1/(1-mu)); % random trajectory lengths
x = length.*cos(angle); % steps
y = length.*sin(angle); % steps
X = cumsum(x,2); % x coordinates of walk
Y = cumsum(y,2); % y coordinates of walk

% Plot walk coordinates
shg
figure(1),h = plot(X,Y,'b.');
axis square
box on
hold off;
zzz=1;

% Cartesian grid
ls = -5000; % [m]. extent to left of origin
rs = 5000; % [m]. extent to right of origin
ts = 5000; % [m]. extent above origin
bs = -5000; % [m]. extent below origin
xnp = 100; % number of x coordinates
ynp = 100; % number of y coordinates
xk = linspace(ls,rs,xnp); % vector of x coordinates
yk = linspace(ts,bs,ynp); % vector of y coordinates
[Xk,Yk] = meshgrid(xk,yk); % X and Y coordinates

% Storage array for counts
grid = zeros(size(Xk)); % [#]. Will contain a count of the
                                             % number of walk coordinates that 'land'
                                             % in each cell defined by the cartesian grid

% Remove random points that lie outside cartesian grid boundaries
[i] = find(X>rs);X(i) = [];Y(i) = [];
[j] = find(X<ls);X(j) = [];Y(j) = [];
[k] = find(Y>ts);X(k) = [];Y(k) = [];
[l] = find(Y<bs);X(l) = [];Y(l) = [];
nrp = numel(X);

% Assign each random point to its nearest cartesian coordinate
for i = 1:nrp
    euc = sqrt((X(i)-Xk).^2+(Y(i)-Yk).^2);
    grid(euc==min(min(euc))) = grid(euc==min(min(euc)))+1;
end
---------------------------------------------------------------

I have tried the following also, but the arrays get too large when I increase n and N to the values I want (e.g., n = 1000, N = 1000):

% Assign each random point to its nearest cartesian coordinate
Xrep = repmat(reshape(X,[1 1 nrp]),[ynp xnp 1]);
Yrep = repmat(reshape(Y,[1 1 nrp]),[ynp xnp 1]);
Xkrep = repmat(Xk,[1 1 nrp]);
Ykrep = repmat(Yk,[1 1 nrp]);
euc = sqrt((Xrep-Xkrep).^2+(Yrep-Ykrep).^2);

for i = 1:nrp
    euk = euc(:,:,i);
    grid(euk==min(min(euk))) = ...
        grid(euk==min(min(euk)))+1;
end

I really don't understand that last for loop. You have X and Y which
is a list of the X and Y coordinates. These are floating point
numbers so you just have to figure out what index they correspond to
and increment it. Something like this
for k =1:nrp
  indexX = X(k);
  indexY = Y(k);
  grid(indexY, indexX) = grid(indexY, indexX) + 1
end

You'd need to modify that if your x and y were in the range (-5000 to
+5000) instead of (1 to nrp) but that's a simple mapping, for example
just add 5000 and divide by 100 or something like that to get the
index in the range of 1-100. This would be many, many fewer
computations than what you are doing.

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <13e0efdd-f5d1-4dbb-9a14-0bffd471b4a3@m25g2000yqc.googlegroups.com>...
> I really don't understand that last for loop. You have X and Y which
> is a list of the X and Y coordinates. These are floating point
> numbers so you just have to figure out what index they correspond to
> and increment it. Something like this
> for k =1:nrp
> indexX = X(k);
> indexY = Y(k);
> grid(indexY, indexX) = grid(indexY, indexX) + 1
> end
>
> You'd need to modify that if your x and y were in the range (-5000 to
> +5000) instead of (1 to nrp) but that's a simple mapping, for example
> just add 5000 and divide by 100 or something like that to get the
> index in the range of 1-100. This would be many, many fewer
> computations than what you are doing.

Hey ImageAnalyst
The reality is that X and Y are not created using meshgrid - I just put that in to try and make my example a bit clearer. X and Y are non-uniformly spaced coordinates - covering a wide range of resolutions and a very large extent.........so it does not seem possible to create a finer grid that encompasses all X,Y pairs (to assist in indexing). I was hoping someone would know a nice trick or command that could help...........

Pete:
Yes I realized that. And I thought my last paragraph was the nice
"trick" that would speed it up for you.

Why calculate the min of a huge bunch of numbers, and then compare it
to that same bunch of numbers just to find an index that you already
have? And then you do this unnecessary stuff not once, but twice.
And it's not only twice, it's twice for every one of your ten thousand
numbers - what a huge waste of time and effort for something that can
be done so simply and efficiently.
-ImageAnalyst

I think I see - however, X and Y form circles of points......each column of X and Y (from left to right) roughly forms a circle of points (using positive and negative coordinates)....as the columns of coordinates go from left to right, the radii of the circles increases.........and the radii do not increase by a constant amount from column to column. Probably i am missing the point, but i can't see a way to simply add on and divide etc. to change the random walk coordinates (x,y) so that they lie between 1:nrp and map to such a spatial grid in a meaningful way.

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <ede0ad28-c333-4083-82c4-420434440c4c@t42g2000vba.googlegroups.com>...
> Pete:
> Yes I realized that. And I thought my last paragraph was the nice
> "trick" that would speed it up for you.
>
> Why calculate the min of a huge bunch of numbers, and then compare it
> to that same bunch of numbers just to find an index that you already
> have? And then you do this unnecessary stuff not once, but twice.
> And it's not only twice, it's twice for every one of your ten thousand
> numbers - what a huge waste of time and effort for something that can
> be done so simply and efficiently.
> -ImageAnalyst

I think I see - however, X and Y form circles of points......each column of X and Y (from left to right) roughly forms a circle of points (using positive and negative coordinates)....as the columns of coordinates go from left to right, the radii of the circles increases.........and the radii do not increase by a constant amount from column to column. Probably i am missing the point, but i can't see a way to simply add on and divide etc. to change the random walk coordinates (x,y) so that they lie between 1:nrp and map to such a spatial grid in a meaningful way.

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <ede0ad28-c333-4083-82c4-420434440c4c@t42g2000vba.googlegroups.com>...
> Pete:
> Yes I realized that. And I thought my last paragraph was the nice
> "trick" that would speed it up for you.
>
> Why calculate the min of a huge bunch of numbers, and then compare it
> to that same bunch of numbers just to find an index that you already
> have? And then you do this unnecessary stuff not once, but twice.
> And it's not only twice, it's twice for every one of your ten thousand
> numbers - what a huge waste of time and effort for something that can
> be done so simply and efficiently.
> -ImageAnalyst

What's wrong with the way I suggested? You basically send out a bunch
of random paths. Then you want to calculate a histogram of the number
that stop at various radii, right? Why then have a 2D histogram? Why
not just histogram along distance? If, instead, you want a histogram
(count) of how many stop at each 2D coordinate (because radius alone
is not enough, and you need direction also), then you can use a 2D
array and use the method I suggested to increment the number that stop
at each location.

Hi ImageAnalyst:
What I need to do is calculate how many points stop at each 2D coordinate.
Imagine the radii of the circles that X and Y produce increase like so:
 r = [0.1 2 13 200 7000 36000 50000.......etc], so we have a bunch of X and Y coordinates covering a whole host of magnitudes, that are non-uniformly spaced, and positive and negative.
How can I obtain the index of random walk coordinates in order to map them to an X and a Y in such a radically non-uniform grid by just adding on a number and dividing by another? If that's possible, I seem to be missing a piece.........
Cheers
Pete



ImageAnalyst <imageanalyst@mailinator.com> wrote in message <78bb22bc-1649-4e5a-a945-c185a5d4d75d@o28g2000yqh.googlegroups.com>...
> What's wrong with the way I suggested? You basically send out a bunch
> of random paths. Then you want to calculate a histogram of the number
> that stop at various radii, right? Why then have a 2D histogram? Why
> not just histogram along distance? If, instead, you want a histogram
> (count) of how many stop at each 2D coordinate (because radius alone
> is not enough, and you need direction also), then you can use a 2D
> array and use the method I suggested to increment the number that stop
> at each location.