Matlab Examples using MDCS

A few examples for Matlab applications using MDCS (prepared using Matlab version R2011b) are illustrated below.

Example application

Consider the Matlab .m-file myExample_2DRandWalk.m (listed below), which among other things illustrates the use of sliced variables and independent stremas of random numbers for use with parfor-loops.

This example program generates a number of N independent 2D random walks (a single step has steplength 1 and a random direction). Each random walk performs tMax steps. At each step t, the radius of gyration (Rgyr) of walk i is stored in the array Rgyr_t in the entry Rgyr_t(i,t). While the whole data is availabe for further postprocessing, only the average radius of gyration Rgyr_av and the respective standard error Rgyr_sErr for the time steps 1...tMax are computed immediately and stored in an output file on HERO.

 
%% FILE:  myExample_2DRandWalk.m
%  BRIEF: illustrate sliced variables and independent streams 
%         of random numbers for use with parfor-loops
%         
%  DEPENDENCIES:
%    singleRandWalk.m   - implements single random walk
%    averageRgyr.m      - computes average radius of gyration
%                         for time steps 1...tMax
%
%  AUTHOR: Oliver Melchert
%  DATE:   2013-06-05
%

N       = 10000;         % number of independent walks
tMax    = 100;           % number of steps in individual walk
Rgyr_t  = zeros(N,tMax); % matrix to hold results: row=radius 
                         % of gyration as fct of time; 
                         % col=independent random walk instances

parfor n=1:N

        % create random number stream seeded by the
        % current value of n; you can obtain a list
        % of all possible random number streams by
        % typing RandStream.list in the command window
        myStream = RandStream('mt19937ar','Seed',n);

        % obtain radius of gyration as fct of time for 
        % different independent random walks (indepence 
        % of RWs is ensured by connsidering different 
        % random number streams for each RW instance)
        Rgyr_t(n,:) = singleRandWalk(myStream,tMax);

end

% compute average Rgyr and its standard error for all steps
[Rgyr_av,Rgyr_sErr] = averageRgyr(Rgyr_t);

  

As liste above, the .m-file depends on the following files:

• singleRandWalk.m, implementing a single random walk, reading:

 
function [Rgyr_t]=singleRandWalk(randStream,tMax)
% Usage:  [Rgyr_t]=singleRandWalk(randStream,tMax)
% Input:
%   randStream - random number stream
%   tMax       - number of steps in random walk
% Output:
%   Rgyr_r     - array holding the radius of gyration
%                for all considered time steps

  x=0.;y=0.; % initial walker position
  Rgyr_t = zeros(tMax,1);
  for t = 1:tMax
    % implement random step
    phi=2.*pi*rand(randStream);
    x = x+cos(phi);
    y = y+sin(phi);
    % record radius of gyration for current time
    Rgyr_t(t)=sqrt(x*x+y*y);
  end
end
  

• averageRgyr.m, which computes the average radius of gyration of the random walks for time steps 1...tMax, reading:

 
function [avList,stdErrList]=averageRgyr(rawDat)
% Usage: [av]=averageRgyr(rawDat)
% Input:
%  rawData   - array of size [N,tMax] where N is the
%              number of independent random walks and
%              tMax is the number of steps taken by an
%              individual walk
% Returns:
%  av         - aveage radius of gyration for the steps

 [Lx,Ly]=size(rawDat);
  avList = zeros(Ly,1);
  stdErrList = zeros(Ly,1);
  for i = 1:Ly
    [av,var,stdErr] = basicStats(rawDat(:,i));
    avList(i) = av;
    stdErrList(i) = stdErr; 
  end
end

function [av,var,stdErr]=basicStats(x)
% usage: [av,var,stdErr]=basicStats(x)
% Input: 
%  x   - list of numbers
% Returns:
%  av  - average 
%  var - variance
%  stdErr - standard error
  av=sum(x)/length(x);
  var=sum((x-av).^2)/(length(x)-1);
  stdErr=sqrt(var/length(x));
end
  

For test purposes one might execute the myExample_2DRandWalk.m directly from within a Matlab session on a local Desktop PC. So as to sumbit the respective job to the local HPC system one might assemble the following job submission script, called mySubmitScript.m:

 
sched = findResource('scheduler','Configuration','HERO');

jobRW =...
   batch(...
        sched,...    
        'myExample_2DRandWalk',...
        'matlabpool',2,...
        'FileDependencies',{...
                           'singleRandWalk.m',...
                           'averageRgyr.m'...
                           }...
        );