Large Scale Parameter Sweep Studies Using Distributed Matlab
Vikas Argod 225C Computer Building vikasargod@psu.edu Ph: (814) Anshuman Gupta Research Engineer / Scientist 222B Computer Building axg218@psu.edu Ph: (814) Research Computing and Cyberinfrastructure (RCC) Information Technology Services (ITS) The Pennsylvania State University, University Park PA 16802 http://gears.aset.psu.edu/

Overview

Typically Ordinary Differential Equations (ODEs) consists of several constants. One of the interesting problems is to find out the nature of variation in values of species for different values of constants. This involves, solving the set of ODEs with different combination of values of the constants till steady state and observing the result for every combination. In this exercise, one of the efficient approaches of solving is discussed. The implementation is done in Matlab. The discussion extends to large scale problems of similar type using distributed matlab. Use of distributed matlab reduced computation time significantly saving 90% of running time using 16 processors compared to single processor runs. The whole exercise has been carried out in LION-XO, a cluster run by Research computing and Cyberinfrastructure Group of ITS at Penn State. Keyword: Set of ODEs, Sensitivity analysis, Distributed matlab
Last updated on 23rd Nov 2008. Feel free to drop an email if you have any questions regarding the content in this article.
Examples are always the best way of explanation. Hence, this article is completely based on one simple example of a set of ODEs. The example considered is explained below followed with solving methodology using matlab. This is followed by large scale implementation using distributed matlab exploiting the beautifully parallel nature of the study. The study is then extended for a bigger set of ODEs representing a reaction network. All the m files used is given in the index.
Now, let us consider a set of ODEs.
dx k1 x k 2 y dt dy k3 y k 4 xz dt dz k5 z k6 xy k7 dt
Initial conditions: x(0) =-8; y(0) = 8; z(0)=27. This consists of 7 constants. Let us consider that these constants can take one of the value from the set S= {10, 25}. Approach and implementation: First step is writing a M file of the ODE set: (F.m) function xprime = F(t,x,k) xprime =zeros(3,1); xprime(1) = - k(1)*x(1) + k(2)*x(2); xprime(2) = -k(3)*x(2)-k(4)*x(1)*x(3); xprime(3) = -k(5)*x(3) +k(6)* x(1)*x(2) - k(7); (The constants ki s are represented in an one-dimensional array.) Overview of the solution methodology is as follows.
Get the combination of constants

Solve set of ODEs

Write the output in desired format
Combination of constants is getting all possible sets of [k1 , k2 , k3 ], where ki S = {10,25}. Most obvious way of getting various combinations is by using nested loops. In the inner most loop, a set of combination gets generated. It will look as follows: for i = 1:2 if (i == 1) k(1) = 10; else k(1) = 25; end for j = 1:2 Consists of nested 5 for loops
. if (p == 1) k(1) = 10; else k(1) = 25; end sol = ode45(@F,T,x0,[ ],k) ; s = deval(sol,15) ; output(q) = [k s]; q = q+1; end end.series of ends to close all the loops. Full listing of the code is given in the appendix-A. The highlighted section is the portion where ODEs are solved and solution is extracted. More explanation is given in the sections that follow.

As you might have observed, this approach gets quite complicated when number of constants increases and also number of values in the set S increases. More importantly, this code can not be run in a parallel processor system which needs independent tasks running in parallel (as solution to the ODEs is done and then loop advances). So, we present a more simple approach which is beautifully parallel. 3
1. Forming various combinations of constants : There are 7 constants and each can take 2 different values. Hence, there are 27 different constant sets. Each number from 1 to 27 corresponds to one unique combination of 7 constants. So there exists a relation between the sets {i, 1 i 27} and combination of ks i.e [k1 , k2 , k3 ]. In order to achieve this, each number is converted to base 2 with 7 digits. p = dec2base (i, 2, 7); where i is the index of loop that runs from i = 1 to i = 27. We associate every digit of this base 2 number with one possible value of constants. Hence, p represents one unique combination of constants. Let us consider for i = 6 for which p = 0000110 Consider that, 1 represents 10, 0 represents 25, from the set S. So, p represents following combination of constants: k = [25]
2. Solving ODE set : Once we have one combination of k, we use ODE45 function of matlab to solve the ODE set. T = [0 15]; % Time interval to solve x0 = [-27]; %Initial value of species sol = ode45(@F,T,x0,[ ],k); This returns a structure that can be used with DEVAL to evaluate the solution in the defined T. s = deval(sol,15); 3. Writing the solution to a text file fid = fopen('results_file.txt','a+t'); output = [i k s']; fprintf(fid,'%4.5f \t',output); fprintf(fid, '\n'); fclose(fid);
Use of distributed matlab for large scale problem:
In the above example, the set S consists of only two values. Assume a case of 5 values in the set, we will have 57 combinations. In the same manner, if number of equations in ODE set increases, computation time for the solution increases. In this case, its a good idea to use Distributed matlab. In distributed matlab, independent tasks are created and submitted for analysis. According to the number of processors (or workers as termed in Distributed matlab), these tasks are executed. Output is collected by a job manager. Once all the tasks are completed, its written to the output file as directed by the user. For more detailed information of the functioning of job manager and workers, please visit http://www.mathworks.com/access/helpdesk/help/toolbox/mdce/ As we are familiar with the above set of ODEs, let us see how we implement this analysis in Distributed matlab. A wrapper script is to be written where in we set environment variables for job manager. A detailed description is given in http://gears.aset.psu.edu/hpc/software/math/matlab/dist.shtml. As our tasks require access to few M files, specifying Path dependencies is done as follows

j = createJob(jm,'FileDependencies',{'~/F.m','~/many_ode.m'});
many_ode.m is the M-file has a function which is called while creating tasks. Tasks are created in a for loop, the index of which is the input argument for the function many_ode. The input argument is converted to base-2 numbering system with 7 digits. This is then parsed and respective combination is generated. These actions are similar to serial matlab code as explained in the previous section. Tasks are created as follows; for i = 1:2^7 createTask(j, @many_ode, 1, {i}); end This is followed by usual Distributed matlab syntax for submitting the tasks, waiting for the results and getting all output arguments, followed by destroying the jobs.
results = getAllOutputArguments(j); fid = fopen('results_file.txt','w+t'); for i = 1:2^7 fprintf(fid,'%4.5f \t',results{i}); fprintf(fid, '\n'); end fclose(fid);
Full listing of the codes is given in the appendix - C Real strength of distributed matlab is seen when we go to a higher scale. Here we present two cases. In the first case we consider a case of solving set of ODEs representing a reacting network It consists of 16 ODEs in 18 variables with 28 coefficients. Sensitivity analysis of 10 coefficients out of these 28 coefficients is carried out, with S = {10, 100}. We have 210 combinations with which ODEs has to be solved. The implementation is in the same manner as explained in the previous section. Effect of batch size is also studied by varying the number of tasks created (in other words, changing the number of iterations solved per task). The time reduction is shown below. Speed up, efficiency and split up of computation time and overhead (consumed by the job manger). No. of Processors 16 Time taken in sec (with batch

size = 32 combinations)

Time taken in sec (with batch
size = 1024/no.of processors)

Time taken in sec (with

batch size = 16)

batch size = 8)

1024 tasks)
5263.243 2876.142 1526.093 812.324 472.573
5459.996 3120.564 1525.884 869.477 456.092
5428.283 3068.050 1499.885 796.308 420.791
5605.596 2802.247 1478.454 754.799 395.185
5567.671 2901.518 1552.727 889.694 552.397
Figure 1: Graph showing reduction in the time for increasing number of processors for different batch size of the input.
Figure 2: Speedup profile for the distributed matlab in LION-XO.
Figure 3: Efficiency of distributed matlab for the example considered and its implementation in LION-XO.
In an attempt of scaling up the exercise, we consider 3 different possibilities for the same coefficients. That means, we have 310 combinations to solve. The time taken for different approaches is shown below: No. of Processors 16 Time taken in min

(with batch size = 9)

Time taken in min

(with batch size = 243)

1540.216 756.422 402.389

1546.149 755.714 414.546

(It was not tried with 1 or 2 processors as it takes not less than 102 hrs with one processor). So, 102 hrs of computation is reduced to 6.7 hrs.
Appendix - A Achieving combinations of constants with nested loops:
function main v = [10,25]; k = [7]; q = 1; x0 = [-27]; for i = 1:2 if(i == 1) k(1) = 10; else k(1) = 25; end for j = 1:2 if(j == 1) k(2) = 10; else k(2) = 25; end for k = 1:2 if(k == 1) k(3) = 10; else k(3) = 25; end for l = 1:2 if(l == 1) k(4) = 10; else k(4) = 25; end for m = 1:2 if(m == 1) k(5) = 10; else k(5) = 25; end for n = 1:2 if(n == 1) k(6) = 10; else k(6) = 25; end for p = 1:2 if(p == 1) k(7) = 10; else k(7) = 25; end sol = ode45(@F,[0 15],x0,[ ],k); s = deval(sol,15); output(q) = [k s']; q = q+1; end end end end end end end
Appendix B Serial matlab code:
function [ ] = many_ode() a n c l u = = = = = 2; 10; [10 25]; 1 2^7;
tic; %Time interval T = [0 15]; %Initial value of species x0 = [8 -8 27]; for i=l:u p = dec2base(i,a,n); for j = 1:n for m=1:a if p(j) == num2str(a-m) ck(j) = c(m); break; else continue; end end end fid = fopen('results_v4.txt','w+t'); %Create a structure of solution sol = ode45(@F,T,x0,[ ],ck); %Extract the solution at specific time output = deval(sol,15); s = [i ck output']; fprintf(fid,'%4.5f \t',s); fprintf(fid, '\n'); end fclose(fid); toc
The input file for running the above serial code in batch (input.m) many_ode The PBS script :
#PBS -l nodes=1:ppn=1 #PBS -l walltime=08:30:00 #PBS -j oe cd $PBS_O_WORKDIR matlab < input
m file for the function F: function xprime = F(t,x,k) xprime =zeros(3,1); xprime(1) = - k(1)*x(1) + k(2)*x(2); xprime(2) = -k(3)*x(2)-k(4)*x(1)*x(3); xprime(3) = -k(5)*x(3) +k(6)* x(1)*x(2) - k(7);
Appendix C Distributed matlab code with wrapper script:
nprocs = [ getenv('DMATLAB_NPROCS') ] np = sscanf( nprocs, '%d' ) mgr_name = [ getenv('JOBMANAGER') ] mgr_host = [ getenv('JOBMANAGERHOST') ] jm = findResource('jobmanager','Name',mgr_name,'LookupURL',mgr_host); j = createJob(jm,'FileDependencies',{'/home1/vxa137/sensitivity/F.m','/home 1/vxa137/sensitivity/many_ode.m'}); i=0; tic; for i = 1:1024 createTask(j, @many_ode, 1, {i}); end submit(j); get(jm); waitForState(j) results = getAllOutputArguments(j); fid = fopen('results_file.txt','w+t'); for i = 1:2^7 fprintf(fid,'%4.5f \t',results{i}); fprintf(fid, '\n'); end fclose(fid); toc destroy(j);