% Model solution to task 2 of week 11
clear all; close all;

[mass,grad,nodes,element,interior,centers,ncoord,nelem] = FEMdata(5);

n = 2^15;
one = ones(length(interior),1);
gen = load('offtheshelf2048.txt');
slist = 2.^(1:11);
results = [];
for s = slist
    disp(['s = ',num2str(s)]);
    decay = 2.0;
    indices = (1:s)';
    deterministic = indices.^(-decay) .* sin(pi*indices*centers(:,1)') .* sin(pi* indices * centers(:,2)');
    a = @(y) 2 + deterministic'*y;
    rhs = mass(interior,:)*nodes(:,1);
    z = gen(1:s);
    tmp = 0;
    parfor ii = 1:n
        qmcnode = mod(z*ii/n,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

% Least squares fit
x = slist(1:end-1)';
errors = abs(results(end)-results(1:end-1));
lsq = [ones(10,1),log(x)]\log(errors);
lsq(1) = exp(lsq(1));

% Visualize the results
figure;
loglog(x,errors,'r.','MarkerSize',18), hold on
loglog(x,lsq(1)*x.^lsq(2),'b','LineWidth',2), hold off
set(gca,'FontSize',14,'FontWeight','bold')
legend('truncation errors',['slope: ',num2str(lsq(2))]);
title('Truncation error')
xlabel('s');
ylabel('error')