ode parallel processing help

Hi everyone! I was wondering if anyone could help me out. I am trying to map parameter space for a set of coupled oscillators. I can map it, but for one run it sometimes takes the whole day. I have reading about parfor that can be used for parallel processing of iterations to speed up the processing time.

My code is set up such that it stores a phase difference in a multidimensional array for each combination of parameter values, and passes through many parameters (eventually there will be nine independent parameters) and a large range of values for each parameter. Unfortunately I am not understanding the logic of the parfor statement...it gives me an error when I try to store things in my array. I understand that you cannot nest parfors; however, I guess without nested for statements, I can't really understand how the indexing works.

Also I have Matlab 2007b (7b), when I think that parfor changed in some way. However, even if you dont know matlab 2007 that is fine. I have access to 2012 at work.

Here is example of my SLOW code:

A1_vector=(0:0.1:10); number_of_A1s=length(A1_vector);

B1_vector=(-10:0.1:10); number_of_B1s=length(B1_vector);

my_data1=zeros(number_of_A1s, number_of_B1s);

    for i=1:1:number_of_A1s 
        for j=1:1:number_of_B1s 
                tspan=0:pi/96:30;
                x0=2*pi*(rand(1,number_of_oscillators));
                A=A1_vector(i);
                B=B1_vector(j);
                [t,s]=ode15s(@(t,a)coupledvdp2c(t,a,A,B),tspan,x0); 
                if std(s(:,1)-s(:,2))>0.1;
                    my_data1(i,j)=nan;
                else
                    my_data1(i,j)=mod(mean(s(:,2)-s(:,1)),2*pi); 
                end    
        end
    end
    pcolor(my_data1)

Thanks! Amelia

PS can you put a for statement in a parfor statment?

[JJ]

PPS if you could comment also on how to replace more than two nested for statements with parfor that would also be great!--Amelia

[José-Luis]

Could you include the error you get?

[JJ]

what do you mean?

[José-Luis]

"Unfortunately I am not understanding the logic of the parfor statement...it gives me an error when I try to store things in my array."

[JJ]

I am not even sure how to script the parfor to get a similar piece of code as above...I tried to follow http://blogs.mathworks.com/seth/2010/10/17/parallel-computing-with-simulink-running-thousands-of-simulations/ but how to use the same index for multiple parameters...I don't really understand!

When I tried changing the first 'for' into a 'parfor', I got the error:

??? Error: File: solvecoupledvdp2c.m Line: 35 Column: 21 The variable my_data1 in a parfor cannot be classified. See Parallel for Loops in MATLAB, "Overview"

Which I read, and I semi-understand that my_data, which stores data from the parfor index and the for index is not a valid variable inside this for statement.

Answer 1

[José-Luis]

PARFOR loops can not be nested. So rewrite your code to avoid the double for loop. For that you can create a matrix of indices:

idx = NaN * ones(number_of_A1s * number_of_A2s,2);
    counter = 0;
    for k = 1:1:number_of_A1s 
        for l=1:1:number_of_B1s 
             counter = counter + 1;
             idx(counter,1) = A1_vector(k);
             idx(counter,2) = A2_vector(k); 
        end
    end

There is a better way to do the above, but i'm soon on company time... And then you can do your simulations in a single parfor:

parfor k = 1:number_of_A1s * number_of_A2s 
      value_ofA1 = idx(i,1);
      value_ofA2 = idx(i,2); 
      %do your stuff here...     
end

Cheers!

[JJ]

This is helpful!But isn't this parfor loop just giving me parameter combination of when A1 and A2 both are at index i, and not all the other combinations of indices? Maybe I don't understand the counter? And also the two for loops appear to take a long time anyway!

Should I create an index matrix like:

|function practice

B1_vector=(-50:0.1:50);             
number_of_B1s=length(B1_vector);
B2_vector=(-50:0.1:50);              
number_of_B2s=length(B2_vector);     

my_data=zeros(number_of_B1s, number_of_B2s);

idx=zeros(number_of_B1s, number_of_B2s, 2)
    for i=1:1:number_of_B1s
        for j=1:1:number_of_B2s
            idx(i,:,1)=B1_vector(i);
            idx(:,j,2)=B2_vector(j);
            idx(i,j,:)=B;
        end
    end

parfor (k=1:1:number_of_B1s*number_of_B2s)
       A=1;
       B1=B(k,1);
       B2=B(k,2);
       tspan=0:pi/96:30;
       x0=[1 1 2 2];
       [t,s]=ode45(@(t,a)coupledvdp2c(t,a,A,B1,B2),tspan,x0);
       my_data=zeros(number_of_B1s,number_of_B2s);
       if std(s(:,1)-s(:,2))>0.1;
           my_data(k)=nan; 
       else
           my_data(k)=mean(s(:,2)-s(:,1));
       end
       pcolor(my_data)
       toc|
end
Can I nest this parfor loops within these for loops?

--Amelia

[José-Luis]

You are completely right. Sorry, there was a typo in my previous post, by

parfor k = 1:number_of_A1s * number_of_A2s 
      value_ofA1 = idx(i,1);
      value_ofA2 = idx(i,2); 
      %do your stuff here...     
end

I meant

parfor k = 1:number_of_A1s * number_of_A2s 
      value_ofA1 = idx(k,1);
      value_ofA2 = idx(k,2); 
      %do your stuff here...     
end

Then it should work

Cheers!

[JJ]

But it does not seem to have sped anything up maybe because it has these two for loops at the beginning...

It also says in the warning labels that "idx is indexed but not sliced, which might cause unnecessary communication overhead."

My other thought was to try to create a number_of_B1s by number_of_B2s by 2 matrix. and index just the x and y part, so that each index refers to a 1x2 vector, which can be plugged into the set of differential equations. I guess I will experiment with this...

[José-Luis]

It is very strange that the first loops ar so slow, especially if preallocated. Anyway, if you don't want the two for loops in the beginning:

A1_vector=(0:0.1:10)'; nA1s=length(A1_vector);

A2_vector=(-10:0.1:10)'; nA2s=length(A2_vector);

idxVec = NaN * ones(nA1s * nA2s,2);
idxVec(:,1) = cell2mat(arrayfun(@(x) repmat(x,nA2s,1) , A1_vector , 'uniformoutput',0));
idxVec(:,2) = repmat(A2_vector,nA1s,1);

But note that you don't really have to create that matrix of indices. You just need to create A1_vector and B1_vector. Once you are inside your parfor loop, then you can use ind2sub:

parfor k = 1:number_of_A1s * number_of_A2s
      [v b] = ind2sub([number_of_A1s,number_of_A2s],k);

        value_ofA1 = A1_vector(v);
        value_ofA2 = A2_vector(b); 
        %do your stuff here...     
  end

You can test it out for yourself to see which one is faster. I would place my money on creating the vector of indices beforehand...

Cheers!