% Written by Vesa Kaarnioja, 2023
% Demonstration on how to estimate the R.M.S. error for a 100-dimensional integral in the affine and uniform setting with 4 random shifts and a _single_ fixed value for the number of QMC nodes n

clear all; close all;

% Load the FEM data
[mass,grad,nodes,element,interior,centers,ncoord,nelem] = FEMdata(3);

% Set up the diffusion coefficient
s = 100; % stochastic dimension
decay = 2.0; % decay of the input random field
indices = (1:s)';
deterministic = indices.^(-decay) .* sin(pi*indices*centers(:,1)') .* sin(pi* indices * centers(:,2)');
a = @(y) 2 + deterministic'*y;

% Set up the loading
f = @(x) x(:,1);
rhs = mass(interior,:)*f(nodes);

% Import the generating vector
gen = load('offtheshelf2048.txt');
z = gen(1:s); % truncate the generating vector to s terms
n = 2^15; % number of QMC nodes
R = 4; % number of random shifts
one = ones(length(interior),1);

results = [];

% For each random shift, compute the QMC estimate using n cubature nodes
for r = 1:R
    shift = rand(s,1);
    tmp = 0;
    parfor ii = 1:n
        qmcnode = mod(ii*z/n+shift,1)-1/2;
        A = UpdateStiffness(grad,a(qmcnode));
        sol = A(interior,interior)\rhs;
        tmp = tmp + one'*mass(interior,interior)*sol;
    end
    results = [results;tmp/n];
end

% Take the average of the R QMC approximations as the final estimate
qmcavg = mean(results); 
rmserror = norm(qmcavg-results)/sqrt(R*(R-1)); % R.M.S. error estimate
rmserror
