%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MCXLAB - Monte Carlo eXtreme for MATLAB/Octave by Qianqina Fang
%
% Speed benchmarks for MCXLAB (equivalent to the examples in mcx/examples/benchmark)
%
% This file is part of Monte Carlo eXtreme (MCX) URL:https://mcx.space
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% only clear cfg to avoid accidentally clearing other useful data
clear cfg;

gpuinfo = mcxlabcl('gpuinfo');

cfg.gpuid = 1;  % gpuid can be an integer 1-N to specify the n-th GPU
% or it can be a string containing only '0's and '1's.
% exp: if gpuid='1010', it means the 1st and 3rd GPU are both used.

% cfg.gpuid='111';          % on wazu, GPU#1,2,3 are 980Ti, 590 Core1 and Core 2
% cfg.workload=[90,10,10];  % workload distribution between the 3 GPUs

gpuid = cfg.gpuid;

if (ischar(gpuid) && regexp(gpuid, '^[01]+$'))
    gpuid = char(gpuid) - '0';
    lastdev = find(gpuid, 1, 'last');
    if (lastdev > length(gpuinfo))
        error('you specified a non-existant GPU');
    end
    gpuid = gpuid(1:lastdev);
elseif (isnumeric(gpuid) && gpuid > 0)
    gpuid = [zeros(1, gpuid - 1) 1];
else
    error('the specified GPU id is invalid');
end

dates = datestr(now, 'yyyy-mm-dd HH:MM:SS');
mcxbenchmark = struct('date', dates, 'gpu', gpuinfo(gpuid == 1));

try
    hbar = waitbar(0, 'Running benchmarks #1');
catch
    hbar = [];
end
count = 0;

cfg.nphoton = 1e8;
cfg.vol = uint8(ones(60, 60, 60));
cfg.srcpos = [29 29 0];
cfg.srcdir = [0 0 1];
cfg.issrcfrom0 = 1;
cfg.autopilot = 1;
cfg.prop = [0 0 1 1; 0.005 1 0.01 1.37; 0.002 5.0 0.9, 1.0];
cfg.tstart = 0;
cfg.tend = 5e-9;
cfg.tstep = 5e-9;
cfg.seed = 29012392;
cfg.detpos = [29 19 0 1; 29 39 0 1; 19 29 0 1; 39 29 0 1];
cfg.isreflect = 0;
% calculate the flux distribution with the given config

speed = zeros(3);

for i = 1:size(speed, 1)
    [flux, detps] = mcxlabcl(cfg);
    if (abs(flux.stat.energyabs / flux.stat.energytot - 0.1769) > 0.005 || ...
        abs(flux.stat.energytot - flux.stat.nphoton) > 10)
        flux.stat;
        error('output absorption fraction is incorrect');
    end
    count = count + 1;
    if (~isempty(hbar))
        waitbar(count / 9, hbar, 'Running benchmarks #1');
    end
    speed(i, 1) = flux.stat.nphoton / flux.stat.runtime;
end
mcxbenchmark.benchmark1.stat = flux.stat;

if (~isempty(hbar))
    waitbar(count / 9, hbar, 'Running benchmarks #2');
end

cfg.shapes = '{"Shapes":[{"Sphere":{"Tag":2, "O":[30,30,30],"R":15}}]}';
cfg.isreflect = 1;

for i = 1:size(speed, 1)
    [flux, detps] = mcxlabcl(cfg);
    if (abs(flux.stat.energyabs / flux.stat.energytot - 0.2701) > 0.005 || ...
        abs(flux.stat.energytot - flux.stat.nphoton) > 10)
        flux.stat;
        error('output absorption fraction is incorrect');
    end
    count = count + 1;
    if (~isempty(hbar))
        waitbar(count / 9, hbar, 'Running benchmarks #2');
    end
    speed(i, 2) = flux.stat.nphoton / flux.stat.runtime;
end
mcxbenchmark.benchmark2.stat = flux.stat;

if (~isempty(hbar))
    waitbar(count / 9, hbar, 'Running benchmarks #3');
end

cfg = rmfield(cfg, 'shapes');
cfg.srctype = 'planar';
cfg.srcparam1 = [40 0 0 0];
cfg.srcparam2 = [0 40 0 0];
cfg.srcpos = [10 10 -10];

for i = 1:size(speed, 1)
    [flux, detps] = mcxlabcl(cfg);
    if (abs(flux.stat.energyabs / flux.stat.energytot - 0.2551) > 0.005)
        flux.stat;
        error('output absorption fraction is incorrect');
    end
    count = count + 1;
    if (~isempty(hbar))
        waitbar(count / 9, hbar, 'Running benchmarks #3');
    end
    speed(i, 3) = flux.stat.nphoton / flux.stat.runtime;
end
mcxbenchmark.benchmark3.stat = flux.stat;

speed = mean(speed);

for i = 1:length(speed)
    mcxbenchmark.(sprintf('benchmark%d', i)).speed = speed(i);
end

mcxbenchmark.speedsum = sum(speed);

if (~isempty(hbar))
    delete(hbar);
end
