A method to vectorise a call to prod() for lots of arrays of varying length?

Question

So my problem is, I'd like to do this without the for loop. Geting the prod() of multiple vectors but of different lengths.

I am dealing with rays intersecting voxels. I typically have 1e6 rays and 1e5 voxels, but this can vary.

intxRays is a list of rays that have intersected voxels.

gainList is a one dimensional vector, each element has a value that corresponds to a specific ray voxel intersection calculated previously (actually with the help of you lovely lot here).

rayIntxStart and rayIntxEnd are vectors of indices for, where in the gainlist array, each ray's corresponding values start and end (they're all in order).

Here is the code and some examples and expected outputs.

gainSum = zeros(1, 5);

% only interested in the intx uniques
intxSegCtr = 1;

% loop through all of the unique segments
for rayCtr = 1:max(intxRays)

    if rayCtr == intxRays(intxSegCtr)

        startInd = rayIntxStart(intxSegCtr);
        endInd = rayIntxEnd(intxSegCtr);

        % find which rows correspoond to those segements
        gainVals = gainList(startInd:endInd);
        gainProd = prod(gainVals);

        % get the product of the gains for those voxels
        gainSumIdx = intxRays(intxSegCtr);
        gainSum(gainSumIdx) = gainProd;

        % increment counter
        intxSegCtr = intxSegCtr + 1;

    end
end

Example data for five rays and nine voxels. Assume the voxel gain array looked like this (for simplicity) for nine voxels (used in previous step).

voxelGains = 10:10:90;

Now say rays 1 and 3 don't hit anything, ray 2 hits voxels 1 and 2, ray 4 hits voxels 2:7 and ray 5 hits voxels 6:9

intxRays = [2, 4, 5];

gainList = [10, 20, 20, 30, 40, 50, 60, 70, 60 70, 80, 90];

rayIntxStart = [1, 3, 9];

rayIntxEnd = [2, 8, 12];

For these numbers the above code would give as a result:

gainSum = [0, 200, 0, 5.0400e+09, 3.024e+07];

I hope this all makes sense.

When I developed it I was using far smaller ray and voxel numbers and it worked fine. As I'm moving up though, the major bottleneck in my code is this loop. Actually just the gainVals and gainProd assignment is like 80% and 15% of my runtime on their own.

This is the only method I can find that works, padding and the like won't work due to the sizes involved.

Is there a way to get the value I want, without this loop?

Many thanks!

Answer 1

ok this is a very small performance boost, but it might help. for testing the matrix way without the loop a bigger data sample is needed.

These are 3 soultions, your original, an optimized and the optimized way as a oneliner. could you please try if this is already doing something for you?

clear all
% orignial loop through all Rays
intxRays = [2, 4, 5];
gainList = [10, 20, 20, 30, 40, 50, 60, 70, 60 70, 80, 90];
rayIntxStart = [1, 3, 9];
rayIntxEnd = [2, 8, 12];
gainSum = zeros(1, 5);
tic
% only interested in the intx uniques
intxSegCtr = 1;

% loop through all of the unique segments
for rayCtr = 1:max(intxRays)

    if rayCtr == intxRays(intxSegCtr)

        startInd = rayIntxStart(intxSegCtr);
        endInd = rayIntxEnd(intxSegCtr);

        % find which rows correspoond to those segements
        gainVals = gainList(startInd:endInd);
        gainProd = prod(gainVals);

        % get the product of the gains for those voxels
        gainSumIdx = intxRays(intxSegCtr);
        gainSum(gainSumIdx) = gainProd;

        % increment counter
        intxSegCtr = intxSegCtr + 1;

    end
end
toc

clear all
%loop insted of every single one to max just through the intxRays
intxRays = [2, 4, 5];
gainList = [10, 20, 20, 30, 40, 50, 60, 70, 60 70, 80, 90];
rayIntxStart = [1, 3, 9];
rayIntxEnd = [2, 8, 12];
gainSum = zeros(1, 5);
tic
for rayCtr=1:length(intxRays)
    %no if as you just go through them
    %intxRays(rayCtr) is the corresponding element

     startInd = rayIntxStart(rayCtr);
     endInd = rayIntxEnd(rayCtr);
     % find which rows correspoond to those segements
     gainVals = gainList(startInd:endInd);
     gainProd = prod(gainVals);     

    % get the product of the gains for those voxels and set them to the ray
    gainSum(intxRays(rayCtr)) = gainProd;
end
%disp(gainSum);
toc

clear all
%same as above, but down to 1 line so no additional values are generated
intxRays = [2, 4, 5];
gainList = [10, 20, 20, 30, 40, 50, 60, 70, 60 70, 80, 90];
rayIntxStart = [1, 3, 9];
rayIntxEnd = [2, 8, 12];
gainSum = zeros(1, 5);
tic
for rayCtr=1:length(intxRays)
gainSum(intxRays(rayCtr))=prod(gainList(rayIntxStart(rayCtr):rayIntxEnd(rayCtr)));
end
toc