Writing data from parallel processing to cells and merging

I have written a ray-tracing code that tracks intersections of rays with the edges of objects, including Snell's law (i.e. ray deflection at borders). At a specific transmitter position, I always fire a "shotgun" of rays in small angular steps into the direction of the object and record all that come close enough to the receiver. The shotgun is necessary because I do not know in advance which rays hit the receiver (due to deflection). Now I want to use parallel processing to speed this process up. My idea was therefore to split up the "shotgun" into as many parts as I have logical cores and use a parfor loop over theses parts. However, I do not know how to consistently track the variables (such as position) within the parfor iterations, such that I can combine them at the end. A rough outline of the code for clarification (I'm just gonna include the x-positions as a tracked variable here for readability):

%at each transmitter position, we track the sequence of x-positions 
%(i.e. intersections with objects) for several rays which hit the receiver.
%importantly, transmitter 1 and transmitter 2 might have a differing 
%number of rays that contribute, therefore the cells makes it easier 
%to track these properties (since I do not need to make sure that the
%entries of different cells are the same size).
%each element of the cell is going to be a 2D array of size n_i * m_i,
%with i the transmitter position index, n the number of valid rays and m
%the number of intersection points.

%define number of vectors and logical cores
n_vec = 100;
n_processes = 4;
%rows: x,y, columns number of positions
ant_pos = zeros(2,5);

%we initialize the cells
%the x positions where each ray has some interaction
ray_points_x = cell(size(ant_pos,2),1);
%the total lengths of the individual rays
lengths = cell(size(ant_pos,2),1);

%here I define how I will split the arrays tracking the properties of the
%vectors into roughly equal parts
split_idx = fix(linspace(0, n_vec, n_processes+1))

%we loop over the transmitter position
for j = 1:size(ant_pos,2)
    %here we define an array for the sequence of x-positions of all rays
    %the idea is that we delete columns of e.g. rays that do not hit
    %anything, and we extend the first dimension of this array by 1 in each ray
    %tracing iteration
    all_x = zeros(1,n_vec);
    
    %my idea was then to begin the parfor loop here
    parfor k = 1:n_processes
        %I then split up all_x into four (roughly) equal parts
        all_x_cur = all_x(split_idx(k)+1:split_idx(k+1));
        for jj = 1:n_iteration 
            %then I perform a bunch of operations
            %if rays are found that hit the receiver (found_idx), BEFORE I added
            %parfor, I used to save them this way
            if cond
                dummy = ray_points_x{j};
                dummy = [dummy, all_x(:,found_idx)]
                %note that I do not show a part of the code that extends the
                %positions found in previous iterations, i.e. let's say I
                %find the sequence [0.0 2.5 5.0] in the third iteration,
                %and in the fourth I find [0.0 1.4 4.2 5.0], then the
                %former would be extended to [0.0 2.5 5.0 NaN] so the sizes
                %are consistent
                ray_points_x{j} = dummy;               
            end
            %then I would continue with jj until either we are at jj = 10
            %or another condition is met, then:
            if jj == 10
                %here I calculate the "results" of the ray tracing for
                %transmitter index j, similar to this
                lengths{j} = sum(diff(dummy));
                %note that "dummy" here is still ray_points_x{j}, this is
                %just for readability
                break;
                %then we go to the next transmitter position
            end
            %some more operations that I can skip if jj = 10
        end        
    end    
end

Now, the big issue are the two cell arrays ray_points_x and lengths, since the different parfor workers cannot access the same array. I thought about solving the issue by defining a temporary cell and then operate on that instead:

parfor k = 1:4
    lengths_cur = cell{1,1};
    %...
    for jj = 1:n_iteration
        %...
        if jj == 10
            lengths_cur{1} = sum(diff(dummy));
        end
    end
end
%combine? (and would also need to extend to longest cell array)

However, I do not know how I could then combine the four separate "lengths_cur" into one. I guess the same thing could also be done for ray_points_x, but the same issue of combining them persists.

Therefore my question is: how can I properly use parallel processing here? How can I combine these cells at the end of parfor, is that even the most efficient way to do what I want? Is parfor ideal here in the first place?

[Mike Croucher]

It would be much easier to help out here if I could run the code on my machine. Can you post all of it here or on GitHub perhaps?

[Edric Ellis]

As Mike said, it would be helpful to see some working serial code. A general comment though - it is not usually beneficial to try and run a parfor loop over your number of processors - what if that changes. Better to parfor over a natural loop variable in your problem.

Getting data out of parfor is one stumbling block. The other is whether the loop iterations are (provably) independent. Data comes out of parfor either as "sliced" or "reduction" variables. "sliced" variables are those where you assign one (or several) elements in an output array (or cell array...) per loop iteration. "reduction" variables are those where you do something like x = x + thisResult; in the loop.

[Dominik Rhiem]

@Mike Croucher @Edric Ellis I am currently trying out to loop over the transmitter positions (i.e. the loop over j) instead and see if that already solves the issue; I will provide you with a code as soon as I see whether that works or not.

[Dominik Rhiem]

@Mike Croucher @Edric Ellis Hi again, so... it did work now. As I said, I simply replaced the for loop over the transmitter position with parfor (with a few other modifications), and then it worked. Sorry for wasting your time.

Answer 1

[Dominik Rhiem]

I replaced the for loop over the transmitter position (i.e. index j) with a parfor loop instead of trying to insert it as described in the original post. Then it worked.