% function [P, Nl] = mlqmc(Lmin,Lmax,Qdims,Nr,eps,mlqmc_l, alpha,beta,gamma)
%
% multi-level Monte Carlo estimation
%
% P     = value
% Nl    = QMC set size on each level
%
% Lmin  = minimum level of refinement       >= 2
% Lmax  = maximum level of refinement       >= Lmin
% Qdims = QMC dimensions on each level
% Nr    = number of QMC randomisations      > 0
% eps   = desired accuracy (rms error)      > 0 
%
% alpha -> weak error is  O(2^{-alpha*l})
% beta  -> variance is    O(2^{-beta*l})
% gamma -> sample cost is O(2^{gamma*l})    > 0
%
% if alpha is not positive it will be estimated
%
% mlqmc_l = function for level l estimator 
%
% [Pf Pc] = mlqmc_l(l,Qdim,N,U)     low-level routine
%
% inputs:  l    = level
%          Qdim = QMC dimension
%          N    = number of paths
%          U    = QMC points (Qdim,N)
%
% output: Pf, Pc = fine, coarse path values 

function [P, Nl] = mlqmc(Lmin,Lmax,Qdims,Nr,eps,mlqmc_l, alpha_0,beta_0,gamma)

  mlqmc_pt_gen('init','reset');  % initialise QMC point generator

  theta = 0.25; % parameter for MSE split between bias^2 and variance

%
% special section for std QMC for comparisons
%

  if (Lmax<0)
    l  = Lmin;
    Qdim = Qdims(l+1);
    Nl = 0;
    Vl = eps^2;

    while Vl > (1-theta)*eps^2
      Nl = max(2*Nl,1);

      U       = mlqmc_pt_gen(l,Qdim,Nl,Nr);
      [dP, P] = feval(mlqmc_l,l,Qdim,Nl*Nr,U);
      P       = sum(reshape(P,Nl,Nr),1)/Nl;

      Vl = max(0, sum(P.^2)/Nr - (sum(P)/Nr)^2) / (Nr-1);
    end

    return
  end

%
% check input parameters
%
  if (Lmin<2)
    error('error: needs Lmin >= 2');
  end

  if (Lmax<Lmin)
    error('error: needs Lmax >= Lmin');
  end

  if (Nr<=0 || eps<=0 || gamma <= 0)
    error('error: needs Nr>0, eps>0, gamma>0 \n');
  end

%
% initialisation
%
  alpha = max(0, alpha_0);
  beta  = max(0, beta_0);

  Nl(1) = 0;  % initial set size on level 0

  L = 0;  % finest level
  l = 0;  % level to work on with larger QMC set size

  while l>=0
    Qdim = Qdims(l+1);
    N    = max(2*Nl(l+1),1);

    U    = mlqmc_pt_gen(l,Qdim,N,Nr);
    dP   = feval(mlqmc_l,l,Qdim,N*Nr,U);
    dP   = sum(reshape(dP,N,Nr),1)/N;

    Nl(l+1) = N;
    Ml(l+1) =        sum(dP)   /Nr ;
    ml(l+1) = abs(   sum(dP)   /Nr);
    Vl(l+1) = max(0, sum(dP.^2)/Nr - ml(l+1)^2) / (Nr-1);

%
% fix to cope with possible zero values for ml
% (can happen in some applications when there are few samples)
%
    for l = 3:L+1
      ml(l) = max(ml(l), 0.5*ml(l-1)/2^alpha);
      Vl(l) = max(Vl(l), 0.5*Vl(l-1)/2^beta);
    end

%
% use linear regression to estimate alpha, beta if not given
%
    if L>1
      if alpha_0 <= 0
        A     = repmat((1:L)',1,2).^repmat(1:-1:0,L,1);
        x     = A \ log2(ml(2:end))';
        alpha = max(0.5,-x(1));
      end

      if beta_0 <= 0
        A    = repmat((1:L)',1,2).^repmat(1:-1:0,L,1);
        x    = A \ log2(Vl(2:end))';
        beta = max(0.5,-x(1));
      end
    end
%
% double set size if too few QMC points compared to level above
%
    if (L>0 & Nl(1)/Nl(2)>8 )
      l = 1;

    elseif (L>1 & max(Nl(2:L)./Nl(3:L+1))>4 )
      [foo l] = max(Nl(2:L)./Nl(3:L+1));
      l = l+1;
%
% double set size if too few QMC points compared to level below
%
    elseif (L>1 & max(Nl(3:L+1)./Nl(2:L))>1 )
      [foo l] = max(Nl(3:L+1)./Nl(2:L));
%
% double set size if variance condition not met
%
    elseif sum(Vl(1:L+1)) > (1-theta)*eps^2
      [foo l] = max( Vl(1:L+1) ./ (Nl(1:L+1).*2.^(gamma*(0:L))) );
      l = l-1;
%
% increase maximum level if bias condition not met
%
    else
      if L>1
        range = 1;
      else
        range = 0;
      end
        
      rem = max(ml(L+1-range:L+1).*2.^(alpha*(-range:0))) / (2^alpha-1);

%      rem = ml(L+1) / (2^alpha - 1);

      if rem > sqrt(theta)*eps
        if (L==Lmax)
          fprintf(1,'*** failed to achieve weak convergence *** \n');
          l = -1;
        else
          L = L+1;
          l = L;
          Nl(L+1) = 0;
        end
      else
        l = -1;
      end
    end
  end

%
% finally, evaluate multilevel estimator
%
  P = sum(Ml(1:L+1));
end
