%
% function U = mlqmc_pt_gen(l,Qdim,N,Nr,rans)
%
% QMC point generator
%
% inputs: l    = level
%         Qdim = dimension of QMC points
%         N    = number of QMC points
%         Nr   = number of randomisations
%
% output: U = (Qdim,N*Nr) array of QMC points
%
% NOTE: on each level, the randomisation is initialised
%       the first time that level is called -- the details
%       vary depending on the generator.

function U = mlqmc_pt_gen(l,Qdim,N,Nr)

persistent generator Qdims Nrs rank_1_z rank_1_rans sobol_gens

%
% select and initialise generator
%

if ischar(l)
  if strcmp(l,'init')
    Qdims = [];

    switch Qdim
      case 'rank-1'
        generator = 'rank-1';
        load z_d4096_m15_gf1
        rank_1_z = z;
        return;

      case 'sobol'
        generator = 'sobol';
        return;

      case 'reset'
        return

      otherwise
        error ('unrecognised initialisation type')
    end
  end

  error ('unrecognised flag')
end

%
% initialise randomisation on each level, if needed
%

if ( (l+1>length(Qdims)) | (Qdims(l+1)==0) )
  Qdims(l+1) = Qdim;  
  Nrs(l+1)   = Nr;

  switch generator
    case 'rank-1'
      rank_1_rans{l+1} = rand(Qdim,Nr);  

    case 'sobol'
      P = sobolset(Qdim);
      for n = 1:Nr
	sobol_gens{l+1}{n} = scramble(P,'MatousekAffineOwen');
      end

    otherwise
      error('internal error: illegal generator type')
  end

%
% otherwise, check no change in Qdim or Nr for this level
%

else
  if Qdims(l+1) ~= Qdim
    error('change in Qdim value')
  end

  if Nrs(l+1) ~= Nr
    error('internal error: change in Nr value')
  end
end

%
% finally, generate points
%

switch generator

%
%------ rank-1 lattice rule -----------------------
%
  case 'rank-1'
%
% version 1: interleaved family sets
%
%     U = rank_1_z(1:Qdim) * reshape(repmat(1:N,Nr,1),1,Nr*N) / N;
%     U = rem(U+repmat(rank_1_rans{l+1},1,N),1);
%
% version 2: grouped into family sets
%
    U    = repmat(rank_1_z(1:Qdim)*(1:N)/N, 1,Nr);
    rans = reshape(repmat(rank_1_rans{l+1},N,1),Qdim,N*Nr);
    U    = rem(U+rans,1);
%
%------ Sobol points with digital scrambling ------
%
  case 'sobol'
    U = zeros(Qdim,N*Nr);
    for n = 1:Nr
      U(:,(n-1)*N+1:n*N) = net(sobol_gens{l+1}{n},N)';
    end

  otherwise
    error('internal error: illegal generator type')
end

end

