Help me make my code faster

I'm making a grid of 3D positions and using colour to express the number of particles/agents/whatever at that position. Currently I'm using 3 for-loops to cycle through each point and find how many particles at that position. Is there a more MATLAB way to do this using indices or quicker functions? Full code attached, important bit below:

% lim = 6
% LIMS = -lim:lim
[xm,ym,zm] = meshgrid(LIMS,LIMS,LIMS);
density = ones(2*lim+1,2*lim+1,2*lim+1); % 2n+1 to capture entirety of LIMS. (-6 to 6 spans 2*6+1 points)

for a = [10 20] % outermost for-loop because I want to capture probability at two separate time steps, namely t=10 and t=20.
    for i = LIMS
        for j = LIMS
            for k = LIMS
                % i+lim+1 translates -6:6 -> 1:13
                density(i+lim+1,j+lim+1,k+lim+1) = sum(ismember(A(:,:,a),[i j k],'rows'));
                % ^^ this ^^ bit counts how many agents, from A, have that exact position.
                % A is 100000 x 3, with 100 sheets (one per timestep).
                % That's 100,000 agents, where the first column is x dimension, second is y, third is z of course.
            end
        end
    end
    figure; slice(xm,ym,zm,density,0,0,0);
    colormap('hot'); grid minor;
end

The full story is that I'm trying to simulate 100,000 agents starting at the origin and propagating outwards in 1,2,3,4, ... dimensions, with an equal probability of moving in any one direction. This'll just help me visualise this distribution in 3D.

Answer 1

[Guillaume]

Unless I'm misreading your code, you're simply computing the histogram of the rows of A. So, assuming that A is made up solely of integers in the set LIM, then:

for a = [10 20]
    density = accumarray(A(:, :, a) + lim + 1, 1, repelem(2*lim+1, 3));  
    %plotting code
end

[Gilad Gur Harush]

That's beautiful thank you! Can you briefly explain the logic there? (so I can apply it to the 2D version)

[Guillaume]

It's simply computing the histogram. A(:, :, a) + lim + 1 simply offsets the values of A to make it a valid index in the destination matrix (as you've done with density), and accumarray simply sums as many 1 (2nd input) as each combination is found.

accumarray(subs, 1, size_of_matrix)

simply sets the location subs to how many times it is repeated.

[Gilad Gur Harush]

Thanks a lot :)

[Gilad Gur Harush]

Can you perhaps help me normalise the bins? Right now it just represents the count in each co-ordinate.

[Guillaume]

Normalize how?

Isn't it just:

%...
density = accumarray(..);
density = density / sum(density, 'all');  %'all' option requires R2018b or later. Earlier versions: sum(density(:))

[Gilad Gur Harush]

pdf normalisation is done by diving by the area not the sum of all values I don't think. When I use your line

density = density / sum(density, 'all');

then the shape of the distribution becomes flatter. It's basically the 3D version of what's happening here: https://stackoverflow.com/questions/5320677/how-to-normalize-a-histogram-in-matlab.

[Gilad Gur Harush]

It seems like because my bin areas are all equal this should give the same answer actually - thanks :)

[Guillaume]

It would be just as easy to normalise by the bin area. Assuming you have a 3D matrix of bin area, it's simply:

density = density .* binarea / sum(density .* binarea, 'all');

[Gilad Gur Harush]

What would this 3D matrix look like?

[Guillaume]

Well, it would be 2*lim+1 x 2*lim+1 x 2*lim+1 matrix giving you the area of each bin. Of course, if all the bins have the same area, the original formula works and you don't need to bother with that.

[Gilad Gur Harush]

Thought so but good to know anyway thank you :)