I have a small function that gets a vector of N linear indices locations in a matrix and returns N new locations in the same matrix. For each of the input locations, a random angle tet and a distance d from a power-law distribution are randomized. The new ith location is the cell on the matrix that is in distance d(i) in direction tet(i) from the ith input location. The matrix represents a torus map, so if the new location is out of the matrix boundaries, it flows back from the other side.
This function is a part of a much larger program that simulates ecological community. It is called hundreds of millions of times and it is the slowest part (by far) in the program, taking about 10 times more than the next slowest operation.
I have tried the following:
- Convert it to a MEX file
- Casting all variables to the minimum size (i.e. uintX instead of double)
- Using the built-in
GeneralizedParetoDistributionto make a PDF object for the randomization
None of them reduced execution time, all of them did the opposite. Also, I cannot use parallel computing as I already use it at 'higher' level of the program.
I'll be glad for any advice/idea that will make this run faster.
Thanks to several advise from @Cris Luengo I was able to cut some runtime:
- I use another random generator called 'simdTwister' which is apparently significantly faster than the default 'twister'.
- I have changed the modulo conversion so I won't have to index and find the zero values after it.
- I use a predefined list of values from power-law distribution to avoid the power computation on every run.
All other ideas within the comments and the answer was slower than the original code.
This may look like an answer to the question, but it isn't. The reason I'm still puzzled about it is that when I implement the third part (using a list of distances) in my outer function and use it in within (part) of the full program, it suddenly takes longer with this pre-computed list then without it, although I compute the list only once for all runs. I have checked this in different ways and have no idea why it is so. Please let me know if there is something wrong (or inefficient) with the way I implement part 3 in the code below.
Also, I would like to know if there is a way to randomize numbers from a 2D Pareto distribution, so I can avoid randomizing angles and converting them.
I have updated the code below to reflect this changes.
Here is the function code:
function [destinations] = randisp_so(Origins,L,alpha,distances) % select random direction with a distance from power-law distribution and % convert back to vectorized form of the matrix % % ORIGINS is a list of origin cells % L is the size of one side of the map % ALPHA is the scale parameter of the distribution % The destinations vector is the same size as ORIGINS OriginSize = numel(Origins); maxd = floor(sqrt(2)*L/2); % maximum distance for dispersal % converting vetorized represatation to matrix: c = ceil(double(Origins)./L); % ind2sub(q) - columns r = double(Origins)-(c-1).*L; % ind2sub(q) - rows if nargin==4 % choose random length from a given distribution (as array): distance = maxd.*distances(randi(numel(distances),OriginSize,1)); else % List of distances to disperse from power-law distribution: e = (-1./(alpha-1)); RND = (1-rand(OriginSize,1)); distance = maxd.*RND.^e; end % choose a random angle from uniform distribution [0 2*pi] angles = rand(OriginSize,1)*2*pi; new_r = round(r + distance.*cos(angles)); % convert length to cell row new_c = round(c + distance.*sin(angles)); % convert length to cell column % transform rows and columns to a 'close-map' new_r = mod(new_r-1,L)+1; new_c = mod(new_c-1,L)+1; % convert back to vectotized form: destinations = new_r+(new_c-1).*L; %sub2ind(r,c); end
And a little testbench for running this function:
% randisp_so testbench n = 10; L = 2^n+1; % the size of one side of the map MapSize = L^2; % total number of cells in the map N = MapSize*30; alpha = 0.7; rng('shuffle','simdTwister') % the fastest random generator q = uint32(randi(MapSize,N,1)); % List of distances to disperse from power-law distribution: distances = linspace(1,0,MapSize).^(-1./(alpha-1)).'; timeit(@() randisp_so(q,L,alpha)) % 2.5153 sec timeit(@() randisp_so(q,L,alpha,distances)) % 1.9008 sec
And here is a comparison of the profile of a typical run: