%
% function mlqmc_test(mlmc_fn,M, Qdims, N,Nr1,L, Nr,Eps,Lmin,Lmax, fp)
%
% multilevel quasi-Monte Carlo test routine
%
% [Pf Pc] = mlmc_fn(l,dim,N,UQ)     low-level routine
%
% inputs:  l    = level
%          Qdim = dimension of QMC points
%          N    = number of paths/points
%          UQ   = QMC points in unit hypercube (Qdim,N)
%
% output: Pf, Pc
%
% M      = refinement cost factor (2^gamma in general MLMC Thm)
%
% Qdims  = QMC dimensions on different levels,
%          an array of size at least 1+max(L,Lmax)
%
% N      = array of QMC set sizes for convergence tests 
% Nr1    = number of randomisations for convergence tests
% L      = number of levels for convergence tests
%
% Nr     = number of randomisations for MLQMC
% Eps    = desired accuracy array for MLQMC calcs
% Lmin   = minimum number of levels
% Lmax   = maximum number of levels

%
% fp     = file handle for printing output data

function mlqmc_test(mlqmc_fn,M, Qdims, Ns,Nr1,L, Nr,Eps,Lmin,Lmax, fp)

%
% COMMENT: this initial version will use rank-1 lattices
%

%
% first, convergence tests
%

PRINTF2(fp,'\n');
PRINTF2(fp,'**********************************************************\n');
PRINTF2(fp,'*** Convergence tests, kurtosis, telescoping sum check ***\n');
PRINTF2(fp,'**********************************************************\n');

for N = Ns

PRINTF2(fp,'\n*******************************************************\n');
PRINTF2(fp,'*** QMC set size = %5d                            ***\n',N);
PRINTF2(fp,'*******************************************************\n');

PRINTF2(fp,'\n l   ave(Pf-Pc)    ave(Pf)   var(Pf-Pc)    var(Pf)');
PRINTF2(fp,'    kurtosis     check \n-------------------------');
PRINTF2(fp,'--------------------------------------------------\n');

rng('default');    % reset random number generator used for randomisation

mlqmc_pt_gen('init','rank-1');    % select QMC point generator

% mlqmc_pt_gen('init','sobol');     % select QMC point generator

del1 = [];
del2 = [];
var1 = [];
var2 = [];
kur1 = [];
chk1 = [];
cost = [];

for l = 0:L
%  disp(sprintf('l = %d',l))

  Qdim = Qdims(l+1);
  U    = mlqmc_pt_gen(l,Qdim,N,Nr1);

%
% compute all Pf, Pc sample values
%
  tic;
  [dP Pf] = feval(mlqmc_fn,l,Qdim,N*Nr1,U);
  cost = [ cost toc ];

% QMC average over each family

  dP = sum(reshape(dP,N,Nr1),1)/N;
  Pf = sum(reshape(Pf,N,Nr1),1)/N;

% now do usual averages over multiple families

  sums = [sum(dP) sum(dP.^2) sum(dP.^3) sum(dP.^4) ...
          sum(Pf) sum(Pf.^2) ] / Nr1;

  if (l==0)
    kurt = 0.0;
  else
    kurt = (     sums(4)             ...
             - 4*sums(3)*sums(1)     ...
             + 6*sums(2)*sums(1)^2   ...
             - 3*sums(1)*sums(1)^3 ) ...
           / (sums(2)-sums(1)^2)^2;
  end

  del1 = [del1 sums(1)];
  del2 = [del2 sums(5)];
  var1 = [var1 (sums(2)-sums(1)^2)*N*Nr1/(Nr1-1) ];
  var2 = [var2 (sums(6)-sums(5)^2)*N*Nr1/(Nr1-1) ];
  var1 = max(var1, 1e-10);  % fix for cases with var=0
  var2 = max(var2, 1e-10);  % fix for cases with var=0
  kur1 = [kur1 kurt ];

  if l==0
    check = 0;
  else
    check = abs(       del1(l+1)  +      del2(l)  -      del2(l+1)) ...
      /    ( 3.0*(sqrt(var1(l+1)) + sqrt(var2(l)) + sqrt(var2(l+1)) )/sqrt(N*Nr1));
  end
  chk1 = [chk1 check];

  PRINTF2(fp,'%2d   %8.4e  %8.4e  %8.4e  %8.4e  %8.4e  %8.4e \n', ...
          l,del1(l+1),del2(l+1),var1(l+1),var2(l+1),kur1(l+1),chk1(l+1));
end

%
% print out a warning if kurtosis or consistency check looks bad
%

if ( kur1(end) > 100.0 )
  PRINTF2(fp,'\n WARNING: kurtosis on finest level = %f \n',kur1(end));
  PRINTF2(fp,' indicates MLMC correction dominated by a few rare paths; \n');
  PRINTF2(fp,' for information on the connection to variance of sample variances,\n');
  PRINTF2(fp,' see http://mathworld.wolfram.com/SampleVarianceDistribution.html\n\n');
end

if ( max(chk1) > 1.0 )
  PRINTF2(fp,'\n WARNING: maximum consistency error = %f \n',max(chk1));
  PRINTF2(fp,' indicates identity E[Pf-Pc] = E[Pf] - E[Pc] not satisfied \n\n');
end

%
% use linear regression to estimate alpha, beta and gamma
%

L1 = ceil(0.4*L)+1;
L2 = L+1;
pa    = polyfit(L1:L2,log2(abs(del1(L1:L2))),1);  alpha = -pa(1);
pb    = polyfit(L1:L2,log2(abs(var1(L1:L2))),1);  beta  = -pb(1);
gamma = log2(cost(end)/cost(end-1));

PRINTF2(fp,'\n*******************************************************\n');
PRINTF2(fp,'*** Linear regression estimates of MLQMC parameters ***\n');
PRINTF2(fp,'*******************************************************\n');
PRINTF2(fp,'\n alpha = %f  (exponent for MLQMC weak convergence)\n',alpha);
PRINTF2(fp,' beta  = %f  (exponent for MLQMC variance) \n',beta);
PRINTF2(fp,' gamma = %f  (exponent for MLQMC cost) \n',gamma);

end

%
% second, mlqmc complexity tests
%

PRINTF2(fp,'\n');
PRINTF2(fp,'****************************** \n');
PRINTF2(fp,'*** MLQMC complexity tests *** \n');
PRINTF2(fp,'****************************** \n\n');
PRINTF2(fp,'  eps       value  mlqmc_cost   std_cost  savings     N_l \n');
PRINTF2(fp,'--------------------------------------------------------- \n');
 
rng('default');    % reset random number generator

gamma = log2(M);

for i = 1:length(Eps)
  eps = Eps(i);
  [P, Nl] = mlqmc(Lmin,Lmax,Qdims,Nr, eps,mlqmc_fn,alpha,beta,gamma);
  l = length(Nl)-1;
  mlqmc_cost = Nr*( (1+1/M)*sum(Nl.*M.^(0:l)) - (1/M)*Nl(1) );

  [P2, Nl2] = mlqmc(l,-l,Qdims,Nr, eps,mlqmc_fn,alpha,beta,gamma);
  std_cost = Nr*Nl2*M^l;

  PRINTF2(fp,'%.4f  %.4e  %.3e  %.3e  %7.2f ', ...
	  eps, P, mlqmc_cost, std_cost, std_cost/mlqmc_cost);
  PRINTF2(fp,'%6d',Nl2);
  PRINTF2(fp,'%6d',Nl);
  PRINTF2(fp,'\n');
end

PRINTF2(fp,'\n');

end

%
% function to print to both a file and stdout 
%

function PRINTF2(fp,varargin)
  fprintf(fp,varargin{:});
  fprintf( 1,varargin{:});
end
