# MatCarlo

## Matlab and Monte Carlo simulations

This example shows one possible way to run monte carlo simulations on the barley cluster using matlab. The example was lifted from http://en.literateprograms.org/Monte_Carlo_simulation_%28Matlab%29. The technique outlined here is actually much more general than this, as the essence of it is just to run many jobs with different parameters.

The technique shown below makes use of a python "driver script" to generate all of the unique files needed to run all of the jobs. If you follow the steps you can see where the python script creates a directory for a job, populates it with a matlab source file and a qsub script, and then executes qsub. Doing it this way has the advantage that the source code and the output for a given job are all kept as a unit. This makes it very easy to analyze later and to reproduce any given result.

Also listed below is an example of saving the matlab variables to a .mat file (so that at some future time you can run matlab, load up the .mat file and manipulate the data), as well as a picture to visualize the data.

### setup

mkdir /farmshare/user_data/bishopj/mattest
cd /farmshare/user_data/bishopj/mattest

save this python script as driver.py and place it in the directory you created above.

import os

n = 0
for mymu in [ 0.05, 0.06, 0.07, 0.08, 0.09 ]:
for myT in range(10,50,10):
n += 1

matlabstartup = '''
S0      = 1000
mu      = %f
sigma   = .12
T       = %d
nb_traj = 100
step    = 1/255

[S, t]  = simplest_montecarlo0( sigma, T, nb_traj, S0, mu, step);

listOfVariables = {'S', 't'};
save('run%d/mydata.mat', listOfVariables{:});

mcplot = figure('Color',[0.9412 0.9412 0.9412 ]);
plot(t,S, 'linewidth',2);
axis([t(1) t(end) min(S(:)) max(S(:))]);
xlabel('time');ylabel('value');
title(sprintf('%%d trajectories', size(S,2)));
print(mcplot,'-dpng','run%d/figure.png')
''' % (mymu, myT, n, n)

qsubscript = '''#!/bin/bash

#\$ -N run%d
#\$ -o run%d/job.out
#\$ -e run%d/job.error
#\$ -cwd
#\$ -S /bin/bash
##\$ -l testq=1

''' % (n,n,n,n)

os.mkdir('run%d' % n)

runfile = open('run%d/run.m' % n, 'w')
runfile.write(matlabstartup)
runfile.close()

qsubfile = open('run%d/run.submit' % n, 'w')
qsubfile.write(qsubscript)
qsubfile.close()

os.system('qsub run%d/run.submit' % n)

save this file as simplest_montecarlo0.m in the same directory as driver.py:

function [S, t] = simplest_montecarlo0( sigma, T, nb_traj, S0, mu, step)
% SIMPLEST_MONTECARLO - the simplest Monte Carlo simulation with Matlab
% use:
%  S = simplest_montecarlo( sigma, T, nb_traj, S0, mu, step)
% example:
%  [S, t] = simplest_montecarlo( .12, 5, 50, 100, .05, 1/255);

nT = ceil(T/step);
W  = sigma * sqrt(step) * cumsum(randn(nT, nb_traj));
c  = repmat((mu - sigma^2/2) *step * (1:nT)',1,nb_traj);
S  = [repmat(S0,1,nb_traj); S0 * exp( c + W)];
if nargout > 1
t = [0;step * (1:nT)'];
end

alternatively, you can copy the two files from /farmshare/software/examples/montecarlomatlab into your directory.

### running

Lets run it

bishopj@corn:/farmshare/user_data/bishopj/mattest\$ python driver.py
Your job 669323 ("run1") has been submitted
Your job 669324 ("run2") has been submitted
Your job 669325 ("run3") has been submitted
Your job 669326 ("run4") has been submitted
Your job 669327 ("run5") has been submitted
Your job 669328 ("run6") has been submitted
Your job 669329 ("run7") has been submitted
Your job 669330 ("run8") has been submitted
Your job 669331 ("run9") has been submitted
Your job 669332 ("run10") has been submitted
Your job 669333 ("run11") has been submitted
Your job 669334 ("run12") has been submitted
Your job 669335 ("run13") has been submitted
Your job 669336 ("run14") has been submitted
Your job 669337 ("run15") has been submitted
Your job 669338 ("run16") has been submitted
Your job 669339 ("run17") has been submitted
Your job 669340 ("run18") has been submitted
Your job 669341 ("run19") has been submitted
Your job 669342 ("run20") has been submitted

bishopj@corn08:/farmshare/user_data/bishopj/mattest\$ qstat
job-ID  prior   name       user         state submit/start at     queue                          slots ja-task-ID
-----------------------------------------------------------------------------------------------------------------
669323 0.25000 run1       bishopj      r     02/19/2013 14:08:10 test.q@barley-testq.stanford.e     1
669324 0.25000 run2       bishopj      r     02/19/2013 14:08:10 test.q@barley-testq.stanford.e     1
669325 0.25000 run3       bishopj      r     02/19/2013 14:08:10 test.q@barley-testq.stanford.e     1
669326 0.25000 run4       bishopj      r     02/19/2013 14:08:10 test.q@barley-testq.stanford.e     1
669327 0.25000 run5       bishopj      qw    02/19/2013 14:08:07                                    1
669328 0.25000 run6       bishopj      qw    02/19/2013 14:08:08                                    1
669329 0.25000 run7       bishopj      qw    02/19/2013 14:08:10                                    1
669330 0.25000 run8       bishopj      qw    02/19/2013 14:08:11                                    1
669331 0.25000 run9       bishopj      qw    02/19/2013 14:08:12                                    1
669332 0.25000 run10      bishopj      qw    02/19/2013 14:08:13                                    1
669333 0.25000 run11      bishopj      qw    02/19/2013 14:08:14                                    1
669334 0.25000 run12      bishopj      qw    02/19/2013 14:08:16                                    1
669335 0.25000 run13      bishopj      qw    02/19/2013 14:08:17                                    1
669336 0.25000 run14      bishopj      qw    02/19/2013 14:08:18                                    1
669337 0.25000 run15      bishopj      qw    02/19/2013 14:08:19                                    1
669338 0.25000 run16      bishopj      qw    02/19/2013 14:08:20                                    1
669339 0.25000 run17      bishopj      qw    02/19/2013 14:08:22                                    1
669340 0.25000 run18      bishopj      qw    02/19/2013 14:08:23                                    1
669341 0.25000 run19      bishopj      qw    02/19/2013 14:08:24                                    1
669342 0.00000 run20      bishopj      qw    02/19/2013 14:08:25                                    1

Since job run1 is started, lets see if things are running as expected:

bishopj@corn08:/farmshare/user_data/bishopj/mattest/run1\$ cd run1
bishopj@corn08:/farmshare/user_data/bishopj/mattest/run1\$ cat job.out
Warning: No display specified.  You will not be able to display graphics on the screen.
Warning: No window system found.  Java option 'MWT' ignored.

< M A T L A B (R) >
R2012b (8.0.0.783) 64-bit (glnxa64)
August 22, 2012

To get started, type one of these: helpwin, helpdesk, or demo.
For product information, visit www.mathworks.com.

>> >>
S0 =

1000

>>
mu =

0.0500

>>
sigma =

0.1200

>>
T =

10

>>
nb_traj =

100

>>
step =

0.0039

>> >> >> >> >> >> >>

Good. Wait for some time for the jobs to finish...

Then we can examine the job directory in detail:

bishopj@corn08:/farmshare/user_data/bishopj/mattest/run18\$ ls -l
total 3736
-rw-r--r-- 1 bishopj root   22803 Feb 19 14:13 figure.png
-rw-r--r-- 1 bishopj root      52 Feb 19 14:13 job.error
-rw-r--r-- 1 bishopj root     762 Feb 19 14:13 job.out
-rw-r--r-- 1 bishopj root 3781664 Feb 19 14:13 mydata.mat
-rw-r--r-- 1 bishopj root     475 Feb 19 14:08 run.m
-rw-r--r-- 1 bishopj root     195 Feb 19 14:08 run.submit

the run.submit file contains the commands to start matlab:

bishopj@corn08:/farmshare/user_data/bishopj/mattest/run18\$ cat run.submit
#!/bin/bash

#\$ -N run18
#\$ -o run18/job.out
#\$ -e run18/job.error
#\$ -cwd
#\$ -S /bin/bash
##\$ -l testq=1

the run.m file contains the matlab source program. Note that the save function is used to save the state of the S and t arrays as well as a plot of the data is saved to figure.png

bishopj@corn08:/farmshare/user_data/bishopj/mattest/run18\$ cat run.m

S0      = 1000
mu      = 0.090000
sigma   = .12
T       =   20
nb_traj = 100
step    = 1/255

[S, t]  = simplest_montecarlo0( sigma, T, nb_traj, S0, mu, step);

listOfVariables = {'S', 't'};
save('run18/mydata.mat', listOfVariables{:});

mcplot = figure('Color',[0.9412 0.9412 0.9412 ]);
plot(t,S, 'linewidth',2);
axis([t(1) t(end) min(S(:)) max(S(:))]);
xlabel('time');ylabel('value');
title(sprintf('%d trajectories', size(S,2)));
print(mcplot,'-dpng','run18/figure.png')

We can look at the png via:

bishopj@corn08:/farmshare/user_data/bishopj/mattest/run18\$ display figure.png

We can also peek inside the mydata file and access the S and t data saved there.

bishopj@corn08:/farmshare/user_data/bishopj/mattest/run18\$ matlab -nodesktop

Warning: No display specified.  You will not be
able to display graphics on the screen.
Warning: No window system found.  Java option 'MWT' ignored.

< M A T L A B (R) >
R2012b (8.0.0.783) 64-bit (glnxa64)
August 22, 2012

No window system found.  Java option 'MWT' ignored.

To get started, type one of these: helpwin, helpdesk, or demo.
For product information, visit www.mathworks.com.

>> whos -file mydata
Name         Size               Bytes  Class     Attributes

S         5101x100            4080800  double
t         5101x1                40808  double

>> t

t =

0
0.0039
0.0078
0.0118
0.0157
0.0196
0.0235
0.0275
0.0314
0.0353
0.0392
0.0431
0.0471
0.0510
0.0549
0.0588
0.0627
0.0667
0.0706
0.0745
0.0784
0.0824
0.0863
0.0902
0.0941
0.0980
0.1020
0.1059
0.1098
0.1137
0.1176
0.1216
0.1255
0.1294
0.1333
0.1373
0.1412
0.1451
0.1490
0.1529
0.1569
0.1608
0.1647
0.1686
0.1725
0.1765
0.1804
0.1843
0.1882
0.1922
0.1961
0.2000
0.2039
0.2078
0.2118
0.2157
0.2196
0.2235
0.2275
0.2314
0.2353
0.2392
0.2431
0.2471
0.2510
0.2549
0.2588
0.2627
0.2667
0.2706
0.2745
0.2784
.
. [data trimmed]
.
20.0000

Don't forget to look for errors

bishopj@corn08:/farmshare/user_data/bishopj/mattest/run18\$ cat job.error
No window system found.  Java option 'MWT' ignored.