%% Agenda
% Multi threading
% Parallel pools
% cluster profiles (intro)
% Prallel Server cluster profiles
% parfor construct
% spmd construct
% distributed construct
% GPU computing
% OPTIONAL case study exercises
%
%%
%% Multithreading
%
% Matlab can execute a large number of operators and functions in parallel through multi treading
%% Does anyoone know what multi threading is and/or how Matlab does it?
%
% Examples of Matlab functions that are suitable for multi-threading.
% element wise operators (e.g. every dot function)
% standard operators (e.g. * , + ,  -. / )
% This is all done completely transparant to user.
% By default Matlab will use as many threads as are available on the desktop/laptop/node

% to check the number of computational thread
num_threads=feature('NumThreads')

%% Is the answer what you expected? Anyone knows why it returned this value?
%

% to set the number of computational threads
num_threads=feature('NumThreads', 2)

%  In this case we set the numner of threads to 2.
%% Why is it important when working on Grace to control the number of threads?
%
%   Average Desktop has 4 to 8 cores. Grace has 48 core
%
%% Example: simple matrix multiplication
%
% Let's do a simmple matrix-matrix multiplication. This one is very suitable for multi-threading?
%% Add code to vary the number of threads and see how long it takes to perform the matrix-matrix mult? For fun, let's also change the size of the matrix;
%

% create random matrix
A=rand(15000);
feature('NumThreads',8)
for count = 1:3
    tic;
    C=A*A;
    mytime = toc;
    fprintf('Computation time is %8.7f seconds.\n',mytime);
end

%
%% PARALLEL POOLS
%
%% The Main concept in parallel processing is a Matlab worker. A worker is basically an indepedent Matlab process. The idea is to utilize a number of workers to split a program in smaller parts and all the workers do their part. Hopefully speeding up the computation.
%% Starting Matlab workers


% start a parallel pool
N=2;
mypool = parpool(N);

%
%% How to check if they were actually started?
%
% For documentation about parpool and all the options:

%check documentation of parpool
doc parpool

%
%
%

% let's check the parpool object
mypool

%
%%
%% Retrieving the parallel pool object
% Matlab provides an easy way to retrieve way to retrieve the pool object. Handy when the pool object was cleared. There is no real need to assign the parpool object to a variable

mypool = gcp('nocreate') % gcp --> Get Current Pool

%
%% Deleting a pool
% When you are done with the workers, don't forget to clean them

delete(mypool) % or mypool.delete()

%
%% Open a pool again; this time without specifyfing the number of workers, what will happen. How many workers will be started?

parpool()  
delete(gcp('nocreate'));

%
% We  used 'nocreate' twice above, so what does the 'nocreate' mean, what does it do?
% We just deleted the pool. What will happen if we call gcp:
% with 'nocreate'
% without 'nocreate'

% first try with 'nocreate'
mypool = gcp('nocreate')

% now try without 'nocreate'
clearvars
mypool1 = gcp()

%
%
%
%% CLUSTER PROFILES
%
% In the previous section we inttroduced parallel pools. Now we will get a bit deeper into this. So far we created a pool and all the workers were created on the same node. Matlab actaully provides a lot of features to decide where and how workers should run. Some things a cluster profile does:
% How and where to start workers
% keep track of all "jobs" (a parpool is also considered a job)
% So far only used 'local' profile
% all workers wil/must run on the same node
% Part of the Matlab Paralel Toolbox
% limit up to to 512 in theory ( in practice limit is number of cores on node/computer)
%

% check documentation again, for cluster object
doc parpool

%


%to see all available cluster profiles
parallel.clusterProfiles

%to see the default cluster profile
parallel.defaultClusterProfile()

%
%% Now, let's take a look how to find all the cluster profiles in the GUI
%
%
% The next functions show how to get cluster objects, and how to change them

%to see all available cluster profiles
cluster_object = parcluster('local')

%% What do we see here?


% let's make some changes to this cluster profile
cluster_object.NumThreads=4;
cluster_object.NumWorkers=2;
cluster_object


%
%


% to make the changes persistent 
cluster_object.saveProfile()

% to save it as under different profile name
cluster_object.saveAsProfile('mylocal');

% check all the cluster profiles 
parallel.clusterProfiles()

%to change the default cluster profile
parallel.defaultClusterProfile('mylocal');
parallel.defaultClusterProfile()

%let's justset default back to local
parallel.defaultClusterProfile('local');

%
%
%% Let's check the cluster profiles in the GUI again
%
%% Let's also take a look at the parallel preferences
%
%
%%  Parallel Server: Cluster Profiles
%
% workers can be distributed over multiple nodes
% number of workers limited only by MDCS license
% integrates with batch schedulers (SLURM on Grace)
%%
%% Importing cluster profile
% Matlab comes standard with the 'local' profile. If you want to use addtional profiles you will have to create them or import them. HPRC already created a cluster profile for you, so you only need to import it

% to import a cluster profile
% parallel.importProfile("location of profile");

% tamu convenience function
profile = tamuprofile.importProfile();

%% Let's check the newly created cluster profile  in the GUI again
%

% check the profile
profile

% Let's get a cluster object for this pofile, similar to what we saw in the GUI?

%

% now let's check all the 'AdditionalProperties'  we saw in the GUI
% we can guess from the output of cluster object

%
%%
%% TAMUClusterProperties
%
% The TAMUClusterProperties class is a Matlab class developed by HPRC, to easily set requirement and properties (e.g. how many workers to use). The parallel cluster object will use this information when starting a parpool or batch job.

% see the help page
help TAMUClusterProperties

%

% initialize a TAMUClusterProperties object
tp = TAMUClusterProperties()

%
% Set some of the properties. NOTE: some of the properties that can be set are needed to run on our cluster (e.g. wall time and memory)

% set requirements; 
tp.workers(8) % how many workers to start
tp.workers_per_node(2) % how many workers to start per node
tp.walltime('03:00')  % how long to use the workers

%
% Retrieve cluster object. Will use that later to start a parpool or offload work using batch

% to get a fully populated cluster object
% created by HPRC
profile = tamuprofile.parcluster(tp)

%
%%
%% Parpool revisited
% let's open a parpool again. This time using our cluster profile. We will see it's avery similar to starting a pool using the 'local' cluster profile

% open a parallel pool using the cluster profile
delete(gcp)
% is there anything that we need to do before creating pool
pool = parpool(profile) % no need to specify #workers

%% Let's check where the parpool is running
%
%%
%% HPRC convenience functions
% HPRC developed convenience functions that use the TAMYClusterProperties object directly, so there is no need to get a cluster objet
%

% delete the current pool
delete(gcp);

%start the pool again using the convience function
tamuprofile.parpool(tp)

%
%%
%% Cool things you can do with cluster profiles
% We can do some cool things when we start workers on multiple nodes
% every worker can use multi threading for aditional speedup
%  each worker can utilize different GPU (assuming each node has GPU)
% every worker can use all the memory on node, increasing total memory availble
%%
%%
%
%% PARFOR
%
% In the previous section we started a parallel with Matlab workers. now we will actually start utilizing them. The first concept we will look at is parfor
% parfor is probably the most common and straightforward parallel construct in Matlab
%% Example 1

% let's turn this regular for loop into a parfor loop
clearvars;
tic
n = 200; A = 300; a = zeros(n);
for count = 1:n
    a(count) = max(abs(eig(rand(A))));
end
toc

%
%
%
%
%% You can adjust the number of workers in parfor
% (why would you set it)?

clearvars;
n = 200; A = 300; a = zeros(n);
tic;
% use the format below to specify the number of workers to use in parfor
parfor (count = 1:n, 2)
    a(count) = max(abs(eig(rand(A))));
end
toc

% EXERCISE: Let's change the above code to make it a perfomance study by 
% varying the number of workers 1 ,2 ,3 ,4, 8 (depends how many workers you have)

%
%
%% Limitations with parfor
% Loop index must be consecutive increasing integers ( e.g. i= N : 0 is not allowed)
% there is always some overhead in using parfor (both starting and iterating)
% Not every for loop can be changed into a parfor loop.

%classifying variables, what happens if we make this a parfor
N=1000;
a=rand(N);
for i = 2 : N
    a(i) = a(i) + a(i-1);
end
a

%
%
%% DEMO: processing number of input files concurently
% Suppose we have a large number of input files (e.g. input1.mat, input2.mat ,....) that we want to process, normally will run them one by one, as follows:
% Read input file
% Run matlab on this input file
% Write results back to disk
% To process all, we just put a for loop around it and run it (or just run multiple matlab instances). How can we process the input files in parallel?

N=4;

% create dummy input files for our demo
for i=1:N
    myinvar = ones(100,1) * i;
    save(['input' num2str(i) '.mat'],'myinvar');
end

clearvars;
N=4;
% let's process the input matrices one by one
for count=1:N
    load(['input' num2str(count) '.mat']);
    % normally this would be where all processing is done
    myoutvar = myinvar*10;
    % save output to file
  
    save(['output',num2str(count),'.mat'],'myoutvar');
end

% create a parallel version of this below

% The above demo is just one example how to concurrently apply a function on multiple inpus. If the inputs were variables/values instead of files the approach would be similar (even a bit more straight forward).
%
%% SPMD
% The next parallel construct we will discuss is the SPMD construct. In an spmd block all the workers will execute the block concurrently. SPMD blocks are actually extensively used by other parallel constructs.
%% EXAMPLE 1

clearvars
N=2
spmd
    worker_id=labindex; %special variable labindex
    num_workers=numlabs; %special variable numlabs
    A=zeros(N);
    A=A+labindex
end

%
%% Data persistence
% One of the features of spmd blocks is data persistence. Even after the spmd block is finished, ALL variables created on the workers are still available

% worker_id, num_workers, and A are still accessibl on the client
whos

worker_id

% to access individual elments use the {} operator
A{2} % will return the value of array A on worker 2

% can even change the value
A{2}=magic(N);

A{2}



%
% variables created inside a spmd block will stay on the workers after the spmd block finished.

spmd
    B=A*numlabs
end
B

%%
%
%% Communication between workers
%
% Matlab provided a number of functions for Matlab workers to send data to eachother
% labSend(var,id)         \to send variable var to worker with id id
% var = labReceive(id) \to receive variable from worker with id id and assign to var
% vf =  labSendReceive(it,if,vt) \to \;combine lab send and receive
%


% Example worker whith higest id sends id to worker 0
spmd
    if (labindex == numlabs)
        sent=rand(4)*labindex;
        labSend(sent,1);
    elseif (labindex == 1)
        received = labReceive(numlabs);
        fprintf('worker 1 received %f from worker %i',received,numlabs);
    end
end
% what will variable received show
received

%what happens if we try to access received{2}
%received{2}

%
%%
%
%
%
%
%% Distributed Arrays
%% Distributed array are conceptually the same as regular arrays/matrices. However, Under the hood the data is actually distributed among all the workers.
% most regular matrix operators will work on distributed array (overloaded functions)
% matrix elements can be of any type
%% How to create distributed arrays?
% create distributed array from regular array
% use the Matlab constructor functions (similar to ones(), zeros(),... )
% create directly on individual workers (ADVANCED, will do that in next section)

clearvars;
N=4;
M=8;
% create distributed array from regular array
A=ones(N,M) .* (1:M);
Ad=distributed(A);

% create distributed array directly on workers
Bd=distributed.rand(N);

whos

% to check if a variable is distributed
a_distributed = isdistributed(A)
b_distributed = isdistributed(Ad)

classAd = class(Ad)

Ad

%
% To create a regular array from a distributed array:
%

Ac = gather(Ad);
classAc = class(Ac)
whos

%

% to see all available operators and methods
methods('distributed')

%
%% How is the array actually distributed among the workers?
% As mentioned, distributed arrays, although conceptually regular arrays, are actually distributed on the workers. Inside spmd blocks we can access them directly.

spmd
    % to see how an array is distributed
    local = getLocalPart(Ad)
    
    % to see global indices along the dimensions
    global_rows = globalIndices(Ad,1)
    global_columns = globalIndices(Ad,2)
end

%
%%
%
%
%% GPU computing
% GPU computing is getting more and more popular and utilizing GPUs can speedup certain codes dramatically.
% What exactly is a GPU
% Accelerator card
% Thousends of simple computing cores
% Dedicated high speed memory
% How to check if a gpu is available:

% check how many GPU's are available on computer/node
number_gpus = gpuDeviceCount()
% list all the gpus
if (number_gpus > 0)
    gpuDevice()
end

%
%% Matlab provides

methods('gpuArray')

%
%% How to use the GPU for computation.
%
% Matlab will automatically perform a computations on the the GPU is the opexecute a function on the GPU if the operands involved in the computation are GPU variables (i.e. reside on the GPU)
%
% Since a GPU has its own memnory and Matlab only utilizes the GPU when the operands are located on the GPU, we need a way to transfer variables to the GPU. Using Matlab, you can do the following
% copy variable to the GPU
% create variable directly on GPU
% copy variable back to client

% copy variable from client workspace to the GPU
N=100;
A=rand(N);
Ag=gpuArray(A);

% create variable directly on the GPU
Bg = gpuArray.ones(N);

% copying variable back to client workspace
Ac=gather(Ag);

whos

%
%% Let's run the matix mult example from the multithreading section and compare run time

% the cpu version
% 
A=rand(16000);

for count = 1:3
    tic;
    C=A*A;
    mytime = toc;
    fprintf('Computation time is %8.7f seconds.\n',mytime);
end

%

% the gpu version
% what do we need to change:
A=rand(16000);

for count = 1:3
    tic;
    C=A*A;
    mytime = toc;
    fprintf('Computation time is %8.7f seconds.\n',mytime);
end

%
%
%% CASE STUDY: COMPUTE PI USING MONTE CARLO METHOD
%
%% Assume a circle and an enclosing square box . If you pick N random points, approximately N*\pi/4 of these points will be inside the circle. Assuming the number of random points is large enough, the following formala can be used to approximate \pi
%% \pi =4*\frac{\sum_1^N \left({x\left(i\right)}^{2\;} +{y\left(i\right)}^2 \le 1\right)}{N}
%
%%
%% MONTE CARLO PI: PARFOR
%
%% A straight forward approach to code this in Matlab is:
% generate N randon cooordinates x  (all between 0\le x\le 1), do the same for y
% for every random point (x,y) check if it's  inside the circle using the formula \sqrt{x^2 +y^2 \;}\le 1
% \pi can be approximated using the formula \pi =4*\frac{\textrm{#points}\;\textrm{inside}}{N}
% Let's proivde step 1:

N=1000000000
% step 1: generate large number of random points with X and Y coordinate
xcoord=rand(N,1);
ycoord=rand(N,1);

% step 2: iterate over all the random points and check if it's inside the loop
tic;
total_inside=0;
for count = 1 : N
 
end
% step 3: compute pi
pi = 4*total_inside/N
mytime=toc;
fprintf('Computation time is %8.7f seconds\n',mytime);

%
%% MONTE CARLO PI: SPMD
%
% Every worker will create a subset of X and Y coordinates
% Ever worker will compute number of points inside
% add the partial sums ( either client aggreates the composite or workers send their results to each other, i.e. list communication)
% Client computes pi

%N=1000000000;
tic;
spmd
    % worker creates suset of X and Y coordinates
    NS = N/numlabs;
   
    tc=0;
    % compute inside points
    for i=1:NS
        if (xcoord(i)^2 + ycoord(i)^2 <= 1)
            tc=tc+1;
        end
    end
end
% 3 client combines partial results
global_inside=0;
for i=1:length(tc)
    global_inside = global_inside+tc{i};
end
%client computes pi
piresult = 4*global_inside/N;
toc
fprintf('The computed value of pi is %8.7f.\n',piresult);

%
%
%% MONTE CARLO PI: MULTITHREADING
%
% Now let;sclook at multi threading. We saw before many matlab operators and build-in functions can run multi=threaded. So far we have used a for loop, which is not suitable for multi-threading. Let's redo the monte carlo method again, instead of using a for use only matlab operators.

X=rand(N,1);
Y=rand(N,1);
tic;

mytime=toc;
fprintf('Computation time is %8.7f seconds\n',mytime);

%
% Now let's run the code and vary the number of threads.
%
%%
%% MONTE CARLO PI: DISTRIBUTED ARRAYS
%
% in the previous section we computed pi, both using a for loop and using matlab operators. Now we will do it using distributed arrays. We will do it both using matlab oprators and directly one the workers

% directly
clearvars
N=8000000;
% step 1: generate large number of random points with X and Y coordinate

% step 2: Use the element wise operators to compute sum of points in circle
pi = 4*inside/N;

tic;

mytime=toc;
fprintf('Computation time is %8.7f seconds\n',mytime);

%
%%
%% MONTE CARLO PI: GPU
%
% We will compute pi using the monte carlo method. This time we will compute it on the gpu and compare results with the cpu versions.


%PI
N=8000000;
% step 1: generate large number of random points with X and Y coordinate


% step 2: Use the element wise operators to compute
tic;


mytime=toc;
fprintf('Computation time is %8.7f seconds\n',mytime);

% WHAT HAPPENS TO TIME WHEN OUTPUT IS NOT SURPRESSED

% use gputimeit to get precise time measurement
%f = @() pi_gpu(N);
%gtime = gputimeit(f)

%
%%
%% THAT'S ALL FOLKS, BACK TO THE SLIDES
%
%



%



%



%



%
%% THIS SECTION CONTAINS ALL THE FUNCTIONS CALLED IN MATLAB SNIPPETS ABOVE

function [pi_val ,mytime] = pi_fun(N , numthreads)
feature('NumThreads',numthreads);
X=rand(N,1);
Y=rand(N,1);
tic;
pi_val = 4 * sum(X.*X + Y.*Y <= 1)/N;
mytime=toc;
fprintf('Computation time is %8.7f seconds\n',mytime);
pause(20)
end

\
function res= pi_gpu(N)
% step 1: generate large number of random points with X and Y coordinate
X=gpuArray.rand(N,1);
Y=gpuArray.rand(N,1);

% step 2: Use the element wise operators to compute
res = 4 * sum(X.*X + Y.*Y <= 1)/N;
end

