%% Agenda
% Parallel pools
% parfor construct
% spmd construct
% distributed construct
% codistributed construct
% GPU computing
% local cluster profiles
% MDCS cluster profiles
% Multi threading ( time permitted)
%%
%
%
%% MATLAB Parallel Toolbox
%% Introducing the MATLAB Parallel Toolbox
%
% Parallel Computing Toolbox™ lets you solve computationally and data-intensive problems using multicore processors, GPUs, and computer cluster.  The toolbox lets you use the full processing power of multicore desktops by executing applications on workers (MATLAB computational engines) that can run locally and on compute clusters.
%% Starting Matlab workers

clearvars
N=2;
tic;
% start a parallel pool
mypool = parpool(N);
toc
whos

%
% (See detailed information slides explaining parpool)
% (use top command to show the workers)
%
% For documentation about parpool and all the options:

doc parpool
% help parpool

%
%
% Use the pool object to get information


mypool
mypool.NumWorkers()

%
%% 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

clearvars
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) and close using gcp

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()

% (show parallel preferences)
%
%%
%% PARFOR
% In the previous section we started additional 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

% (See detailed information slides explaining parfor)
%
%
%
%% 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)

%
%
%%
%% 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 compute \pi
%
%% \pi =4*\frac{\sum_1^N \left({x\left(i\right)}^{2\;} +{y\left(i\right)}^2 \le 1\right)}{N}
% 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}

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

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

%
%% 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

% (See detailed information slides explaining spmd)
%
%% 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}

%
%%
%% EXERCISE: list communication
% Let's write a script where worker with id id waits to receive data from worker id+1 and sends it to worker id-1. Worker 0 will print results.

% Example worker whith higest id sends id to worker 0
spmd
    if (labindex == numlabs)
        
    elseif (labindex > 1)
        
    else
       
        fprintf('worker 1 received %f',received);
    end
end

%
%
%% CASE STUDY: Monte Carlo pi using 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=8000000;
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);

%
%
%% 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

% (See detailed information slides explaining distributed )
%


% 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

%
%% CASE STUDY:  MONTE CARLO PI USING 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);

%
%
%% Codistributed arrays
% Codistributed arrays and distributed arrays are basically the same array, viewed from different perspectives.  Let's use an example to illustrate:

clearvars;
A=ones(8) .* (1:8)
Ad=distributed(A);
class_on_client=class(Ad)
spmd
    class_on_worker=class(Ad)
    %Ad
end

%
%
% Simply put, it's called a distributed array from a global view (from the client) and codistributed from a local view (on the workers). To see if an array is codistributed.

iscodistributed(Ad)
spmd
    iscodistributed(Ad)
end

%
% To see all function that work on codistributed variables

methods(codistributed)

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

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

%
%% CASE STUDY:  MONTE CARLO PI USING DISTRIBUTED ARRAYS
%
% Let's do it again, this time will do all the work directly on the workers/ Let's use Matlab operators again

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

tic;
spmd
    % compute the number of points inside
    NS=length(xcoord);
    tc= sum((xcoord.^2 +ycoord.^2) <= 1);
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);

%
% Another way to do it is using a for loop, similar as we did in the SPMD version above
%
%
%% Creating distributed arrays on workers
% In the previous section we created distribtued arrays on the client (using the "distributed" function or using the constructors. Another option is to create the distributed array directly on the workers. Some reasons are:
% speed : all workers are busy creating array
% size : arrray might be to large to fit into memory (use workers on multiple nodes)
% In the next sections we will create distributed arrays directly on the workers
%% Distributed version of replicated array

clearvars
spmd
    a=ones(4);
    a=a*labindex
    b=codistributed(a);
end
a
b
whos 

%
%%
%% Creating larger dsitributed arrays from local parts
% Next we will create a new larger array from the local parts

clearvars
spmd(4)
    a=ones(4);
    a=a*labindex
    % have to specify how the larger array will be created/partitioned
    % param 1: adimension to partition: 1=row, 2 is colum
    % param 2: #rows(columns) for every worker
    % param 3: dimension of distributed array
    codist1=codistributor1d(1,[4, 4, 4, 4],[16,4]);
    % actually build the distributed array
    c=codistributed.build(a,codist1)
end
c
whos 

% Question: What do we need to change if we want to partition along columns?

%
%%
%% Exercise: Create distributed array where local parts have different dimension
% let's do the same as above, but now the dimensions are not the same. In this exercise, every local part on a worker will asume dimension [labindex,4]. E.g. on worker 1 dimension will be 1x4, on worker 2 will be 2x4, etc

clearvars
spmd(4)
    %a has different dimensions on every worker. 
    a=ones(labindex, 4);
    a=a*labindex;
    % need to set the codistributer
    %codist1=codistributor1d(1, [X, X, X, X],[X,X]);
    % actually build the distributed array
    c=codistributed.build(a,codist1);
end
c
whos 

% The above code does make one assumption? What is it?

%
%
%
%% 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

%
%
%% CASE STUDY:  MONTE CARLO PI USING 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)

%%
%
%% Offloading work using the batch command
% Matlab provides an easy way  offload work to one or more worker(s) in  non-blocking fashion. You can continue your own interactive session while the Matlab worker is executing your offloaded code.


clearvars;
% run a script on one of the workers
% control will be given back immediately
myjob = batch('myscript')

% do some other stuff
disp("hello");
dummy_var=rand(100);
dummy_var = dummy_var * dummy_var;

% let's check the status of the jobs
job_status = myjob.State()

% just wait until the script is finished
myjob.wait()

%this point will only be reached
fprintf('the job is finished');

%%
%% Get information from a Job


startTime = myjob.StartDateTime()
endTime = myjob.FinishDateTime()

% to get all the screen output.
myjob.diary();
whos
% to load all the variables used in the script into current workspace
myjob.load();
whos


%  (show the Job Monitor)
%
%%
%% Run a function on the worker
% besides offloading scripts, you can alsos execute function calls on the workers

N=1000000;
% offload function to workers
j = batch(@pi_fun,2,{N,1})

% can do other things in meantime, but for now , lets wait
j.wait();

% to get the results from the function
output = j.fetchOutputs()
% get run time (second element in cell array)
runtime=output{2}

%
%%
%% Offloading a parallel job
%

N=2000000;
% offload function to workers start a parallel pool
j = batch(@pi_fun,2,{N,1},'Pool',4);

%also try to run two parpool jobs
%j2 = batch(@pi_fun,2,{N,1},'Pool',4);

j.wait();
output = j.fetchOutputs();
runtime = output{2}

%
%% cancel a running job

j = batch(@mydummy,0,{});

% we realize this code has errors or just takes too long
j.cancel();

%look at the state
j.State()

%
%% Don't forget to delete a job once it's not needed anymore


% delete a job
j.delete()

% check if job is still valid
j.isvalid()


%
%% Cluster profiles (local)
% 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()

%
% The next functions show how to get cluster objects, and how to change them

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

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

%


% 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');


% (Show Cluster Profiles)
%
%% Let's revists jobs again
% the cluster profiles keeps track of all jobs (batch and pools) executed/created on it. to retrieve the job information:

% if you know the Job ID
myjob=cluster_object.findJob('ID',1)

% to get all jobs associated with this cluster profile
all_jobs=cluster_object.Jobs()

%
%% MDCS Cluster Profiles: HPRC clusters
% workers can be distributed over multiple nodes
% number of workers limited only by MDCS license
% integrates with batch schedulers (LSF on ada, SLURM on terra)
% (See detailed information slides explaining cluster profile )
%% 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();

% for versions less than R2019b
%tamu_import_TAMU_clusterprofile()

% (Show Cluster Profiles again + show code tamuprofile.importProfile)
%
%% 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
%pool = parpool(profile) % no need to specify #workers

% do little test, print out the names where the workers are running
spmd
    host = getenv('HOSTNAME')
end
for i=1:length(host)
    fprintf('worker %d, on host %s\n',i,host{i});
end

%
%
%%
%% batch revisited
% just like parpool, we can easily offload code to the workers.


N=1000000;
% offload function to workers
% two small differences with before
% 1) specify the cluster object to use (profile)
% 2) set the number of workers to total number of workers -1
j = batch(profile,@pi_fun,2,{N,1},'Pool',profile.NumWorkers-1);

% You can close the current session, the offloaded work will still continue to execute
%
%%
%% HPRC convenience functions
% HPRC developed two convenience functions that use the TAMYClusterProperties object directly, so there is no need to get a cluster object
% tamuprofile.batch(TAMUClusterProperties, script)
% tamuprofile.parpool(TAMUClusterProperties)
%
%%
%% 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
%%
%%
%
%%
%% MATLABSUBMIT
% HPRC created a tool to easily submit matlab jobs from ada and terra login nodes, without the need to write your own batch file/
%
% (demo matlabsubmit)
%
%%
%% Multithreading
% Matlab can execute a large number of operators and functions in parallel through multi treading. For example:
% element wise operators (e.g. every dot function)
% standard operators (e.g. * , + ,  -. / )
% This is all done completely transparant to user.
%% Example: simple matrix multiplication

% create random matrix
A=rand(8000);

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

%   Average Desktop has 4 to 8 cores HPRC cluster terra has 28 cores, HPRC cluster ada has 20 cores
%
%% How to check/set the number of cores
% How many threads did Matlab actually use to execute the above code

% get the number of computational threads
nt = feature('NumThreads')
% set the number of computational threads
nt_new = feature('NumThreads',4)

% alternatively you can use maxNumCompThreads function
nt_alt = maxNumCompThreads()
nt_alt = maxNumCompThreads(8)

%% How to see it's running in parallel?
% on Linux, use the "top" command (might need to open another shell)
% on Windows use the Task manager (will not work on OAL machines)
% Obviously, it should run faster compared to serial
% (demo utilization)
%
%
%
%%
% Let's rewrite the above code using only Matlab operators (i.e. remove the for block). If we can do that we can make use of the multi threading capabilities of Matlab. Vary the number of threads and see how much speedup you can get.

clearvars;
feature('NumThreads',1)
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);

%
%%
%% QUESTIONS?
%
%
%
%% 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 out= mydummy()
pause(210)
out=1;
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

function mysave(fname,myoutvar)
save(fname,'myoutvar','-mat');
end