Best way to parallelize montecarlo simulations of Misanthrope process Matlab

Question

I am studying the diffusion of the misanthrope process with periodic boundary conditions. In this process, I consider a 1-D discrete lattice with M sites and N repulsive-interacting particles; each site can be occupied by an arbitrary number of particles at any time.

I already have working scripts to compute this on my personal workstation (4 parallel workers no GPU) way too slow for my purposes. I do have access to the following shared hardware: 10 core intel core i9-7900X processors [20 threads with hyper-threading] and 128GB of RAM with Nvidia Titan V GPUs with 12GB of GPU RAM.

I am not an expert on GPUs but I would like to take advantage of the Nvidia Titan available somehow. My code is written in Matlab and the shared machine has the Matlab 2018a. The simulations record the time evolution of the ensemble of particles and loops over the size of the ensemble as well as over an interaction strength parameter.

This is nested as follows:

loop over interaction strength parameter (for loop beta)
    loop over N                          (parallelize with parfor)
        loop over trajectories           (for loop n_traj)
            loop over time               (while loop t_f)

The two largest loops are the trajectories and time ones. Would making any of my variables (maybe n_traj) a GPUarray on Matlab speed up my code?

This is the script (for a quick test use: n_traj=100, t_f=30)

Repulsion=[6 5 2 1 0.5];
%--- lattice size
L=2; M=10; x=1:L*M;
%--- declares hopping rates  
k0=0.1; X=2; k_r=exp(X)*k0;   k_l=exp(-X)*k0; 
%--- non-interacting single particle drift and diffusion
D=(k_l+k_r)/2; v=(k_r-k_l);
%--- pre-allocates interacting particles drift and diffusion  
driftMEx=zeros(L*M,length(Repulsion)); 
diffuMEx=zeros(L*M,length(Repulsion));
%--- gets ensemble size and time duration
prompt=sprintf('Give the value No of trajectories:\n n_traj = ');
n_traj = input(prompt);
prompt=sprintf('\nGive the value time lenght:\n t_f = ');
t_f = input(prompt);

j=1;
for beta=Repulsion
    parfor i=1:L*M 
            [driftMEx(i,j),diffuMEx(i,j)]=simulation(i,M,X,k0,beta,n_traj,t_f);
    end 
    j=j+1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                        Functions used in script                         %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%----       Montecarlo Simulations for Misanthrope Process      ----%%%%
function [drift,diffusion] = simulation(N,M,X,k_o,beta,n_traj,t_f)
%%%%%%%%%%%%%%%%%%%%%%% parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k_r=exp(X)*k_o; k_l=exp(-X)*k_o; 
s_min=1; % lowest state in boundary conditions
%%%%%%%%%%%%%%%%%%%%%%%  Pre-allocation  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x_final=zeros(n_traj,1);
P=zeros(2*N+1,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for jj=1:n_traj
    tc=0; % current time starts at 0
    state0=1:N;
    state0=mod(state0,M);
    state0(state0==0)=M;
    state=histcounts(state0,0.5+(0:M));
    bc_counter=0;
    hop=1;
    while tc<t_f  % loop over time
        % findMisanthrope finds index of where motors are this state
        N_pos=findMisanthrope(state,N);  
        T=zeros(N,M); % reset transition matrix
        % temporarily assign state of hopping fowards and backwards
        N_plus=N_pos+1;
        N_minus=N_pos-1;
        nn=histcounts(N_pos,0.5+(0:M));
        % finds possible moves
        for qq=1:N     
            % findMisanthrope finds index of where motors are this state
            N_pos=findMisanthrope(state,N);              
            if N_pos(qq)==1
                N_minus(qq)=M;
            end            
            if N_pos(qq)==M
                N_plus(qq)=1;
            end
            T(qq, N_minus(qq))=k_l; % rate of hopping left
            T(qq,N_plus(qq))=k_r;   % rate of hopping right
            
            alpha=0;
            bb=exp(-beta.*nn);
            aa=exp(alpha.*(nn-1));

            temp_qq=N_pos(qq);
            T(qq,setdiff(1:end,temp_qq))=T(qq,setdiff(1:end,temp_qq)).*bb(setdiff(1:end,temp_qq)).*aa(temp_qq);           
        end

        for mm=2:2:2*N
            P(mm)=T(mm/2,N_minus(mm/2))/sum(T(:)) + P(mm-1);
            P(mm+1)=T(mm/2,N_plus(mm/2))/sum(T(:)) + P(mm);
        end
        tau=(-log(rand))/(sum(T(:)));  % find time until next hop
        tc=tc+tau; % update current time
        R=rand; % pick another random number from 0 to 1 (r_2)
        iii=find(P<=R,1,'last'); % find which entry in vector correspons to r_2
        N_hopp=round(iii/2);  % which motor changes position
            state(N_pos(N_hopp))=state(N_pos(N_hopp))-1; % annihilate motor to hop
            N_plus=N_pos(N_hopp)+1;
            N_minus=N_pos(N_hopp)-1;
            
            if mod(iii/2,1)~=0   % if odd then motor hops left
                if N_pos(N_hopp)==s_min
                    N_minus=M;
                    bc_counter=bc_counter-1; % take note of how many times a motor crosses the boudary
                end
                state(N_minus)=1+state(N_minus); % create motor where it hops to
            else   % even then hops right
                if N_plus>M
                    N_plus=s_min;
                    bc_counter=bc_counter+1;
                end               
                state(N_plus)=1+state(N_plus);       
            end
        hop=hop+1; % advance hop index
    end
    N_f=findMisanthrope(state,N); % final position of motors
    delta=N_f-state0; % find change in state for drift calculation
    x_final(jj)=(sum(delta)+M*bc_counter); % add boundary condition info on. To find total distance travelled
end
% calculates diffusion
ave=mean(x_final);
diff_n(1:n_traj)=(x_final(1:n_traj)-ave).^2;
diffusion=sum(diff_n)/(2*n_traj*t_f);
% calculates drift
drift=sum(x_final(1:n_traj)/t_f)/n_traj;
end


%%%%----     Find indices of particles for Misanthrope Process     ----%%%%
function [rr]=findMisanthrope(tt,m)

    if 1*isempty(tt(tt>=2))==0    
        ttt=tt(tt~=0);
        rrr=find(tt);
        rr=zeros(1,m);
        rrt=cell(1,length(rrr));    
        for i=1:length(rrr)
            rrt{i}=rrr(i)*ones(1,ttt(i));
        end
        rr(1:m)=cell2mat(rrt);    
    else
        rr=find(tt);    
    end

end

Answer 1

After trying some changes to my code, I came up with this new version that runs on GPU. The speed is marginally better about 20% faster. Considering the changes were minimal I guess it is a good trade-off.

%--- pre-allocates interacting particles drift and diffusion  
driftMEx=gpuArray(zeros(L*M,length(Repulsion))); 
diffuMEx=gpuArray(zeros(L*M,length(Repulsion)));

1st change was to allocate my main variables as gpuArrays

j=1;
for beta=Repulsion
    tic;
    parfor i=1:L*M 
            [driftMEx(i,j),diffuMEx(i,j)]=simulation(i,M,1,k_r,k_l,beta,n_traj,t_f);
    end 
    j=j+1;
    toc;
end

the main loop remains the same way. The biggest change was made to function

[drift,diffusion] = simulation(N,M,s_min,k_r,k_l,beta,n_traj,t_f)

where all instances where I used findMisanthrope(state,N); were changed to repelem(1:numel(state), state); which is compatible with gpu

%%%%----       Montecarlo Simulations for Misanthrope Process      ----%%%%
function [drift,diffusion] = simulation(N,M,s_min,k_r,k_l,beta,n_traj,t_f)
%%%%%%%%%%%%%%%%%%%%%%%  Pre-allocation  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x_final=gpuArray(zeros(n_traj,1));
P=gpuArray(zeros(2*N+1,1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for jj=1:n_traj
    tc=0; % current time starts at 0
    state0=gpuArray(1:N);
    state0=mod(state0,M);
    state0(state0==0)=M;
    state=gpuArray(histcounts(state0,0.5+(0:M)));
    bc_counter=0;
    hop=1;
    while tc<t_f  % loop over time
        % findMisanthrope finds index of where motors are this state
%         N_pos=findMisanthrope(state,N);
        N_pos=repelem(1:numel(state), state);
        T=gpuArray(zeros(N,M)); % reset transition matrix
        % temporarily assign state of hopping fowards and backwards
        N_plus=N_pos+1;
        N_minus=N_pos-1;
        nn=histcounts(N_pos,0.5+(0:M));
        % finds possible moves
        for qq=1:N     
            % findMisanthrope finds index of where motors are this state
%             N_pos=findMisanthrope(state,N); 
            N_pos=repelem(1:numel(state), state);
            if N_pos(qq)==1
                N_minus(qq)=M;
            end            
            if N_pos(qq)==M
                N_plus(qq)=1;
            end
            T(qq, N_minus(qq))=k_l; % rate of hopping left
            T(qq,N_plus(qq))=k_r;   % rate of hopping right
            
            alpha=0;
            bb=exp(-beta.*nn);
            aa=exp(alpha.*(nn-1));

            temp_qq=N_pos(qq);
            T(qq,setdiff(1:end,temp_qq))=T(qq,setdiff(1:end,temp_qq)).*bb(setdiff(1:end,temp_qq)).*aa(temp_qq);           
        end

        for mm=2:2:2*N
            P(mm)=T(mm/2,N_minus(mm/2))/sum(T(:)) + P(mm-1);
            P(mm+1)=T(mm/2,N_plus(mm/2))/sum(T(:)) + P(mm);
        end
        tau=(-log(rand))/(sum(T(:)));  % find time until next hop
        tc=tc+tau; % update current time
        R=rand; % pick another random number from 0 to 1 (r_2)
        iii=find(P<=R,1,'last'); % find which entry in vector correspons to r_2
        N_hopp=round(iii/2);  % which motor changes position
            state(N_pos(N_hopp))=state(N_pos(N_hopp))-1; % annihilate motor to hop
            N_plus=N_pos(N_hopp)+1;
            N_minus=N_pos(N_hopp)-1;
            
            if mod(iii/2,1)~=0   % if odd then motor hops left
                if N_pos(N_hopp)==s_min
                    N_minus=M;
                    bc_counter=bc_counter-1; % take note of how many times a motor crosses the boudary
                end
                state(N_minus)=1+state(N_minus); % create motor where it hops to
            else   % even then hops right
                if N_plus>M
                    N_plus=s_min;
                    bc_counter=bc_counter+1;
                end               
                state(N_plus)=1+state(N_plus);       
            end
        hop=hop+1; % advance hop index
    end
%     N_f=findMisanthrope(state,N); % final position of motors
    N_f=repelem(1:numel(state), state);
    delta=N_f-state0; % find change in state for drift calculation
    x_final(jj)=(sum(delta)+M*bc_counter); % add boundary condition info on. To find total distance travelled
end
% calculates diffusion
ave=mean(x_final);
diff_n(1:n_traj)=(x_final(1:n_traj)-ave).^2;
diffusion=sum(diff_n)/(2*n_traj*t_f);
% calculates drift
drift=sum(x_final(1:n_traj)/t_f)/n_traj;
end