% Perform a divand benchmark
% A 2d-variational analysis is performed over a square domain
% as described in http://www.geosci-model-dev.net/7/225/2014/gmd-7-225-2014.pdf
% with different number of grid points (100 x 100, 200 x 200, ... 600 x 600).
% For every domain size, the benchmark is run 10 times and the median time is 
% computed.
%
% Input (optional):
%   name: if name is present then the results will be saved in the file called 
%     name
% 
% Output:
%   median_time: the median of the run-time in seconds
%   ng: number of grid points along one dimension
%   time: 2-d array with all results
%
% To run this benchmark you need to install divand, for example
%
% pkg install -forge divand
% pkg load divand

% Alexander Barth
% GPLv2 or later

function [median_time,ng,time] = test_2dvar_benchmark(name)

  fprintf('Running divand benchmark in 2 dimensions\n');

  % domain sizes
  ng = [100:100:600];

  for j=1:10
    for i=1:length(ng)
      [time(i,j),RMS(i,j)] = benchmark2d(ng(i));
      if (RMS(i,j) > 0.2)
        error('unexpected large RMS. Results might be wrong');
      end

      fprintf('size %5d time %10.4f \n',ng(i),time(i,j));
    end
  end

  fprintf('\nMedian results\n');

  median_time = median(time,2);

  for i=1:length(ng)
    fprintf('size %5d time %10.4f \n',ng(i),median_time(i));
  end

  if nargin == 1
    fname = ['test_2dvar_benchmark_' name '.mat'];
    fprintf('save result in file %s\n',fname)
    save(fname,'time','RMS','ng')
  end

end


function [time,RMS] = benchmark2d(ng)

  mg = ng/5;
  len = 10/ng;

  lambda = 20;

  f = @(x,y) cos(2*pi*ng*x/20) .* cos(2*pi*ng*y/20);


  % grid of background field
  [xi,yi] = ndgrid(linspace(0,1,ng));
  vi = f(xi,yi);

  mask = ones(size(xi));
  pm = ones(size(xi)) / (xi(2,1)-xi(1,1));
  pn = ones(size(xi)) / (yi(1,2)-yi(1,1));

  % grid of observations
  [x,y] = ndgrid(linspace(1e-6,1-1e-6,mg));
  x = x(:);
  y = y(:);
  v = f(x,y);

  xx = repmat(x,[1 numel(x)]);
  yy = repmat(y,[1 numel(y)]);

  tic()
  [va,s] = divand(mask,{pm,pn},{xi,yi},{x,y},v,len,lambda,'diagnostics',1,'primal',1);
  time = toc();
  RMS = rms(va,vi);

end


function r = rms(x,y,norm)

  d = x-y;

  if nargin == 3
    d = d .* sqrt(norm);
  end

  m = ~isnan(d);
  r = mean(d(m).^2);

  if nargin == 3
    r = r/mean(norm(m));
  end

  r = sqrt(r);
end