%% Aufgabe 2: Mittlere Fläche eines zufälliges Dreiecks im Einheitsquadrat
clear all; clc;

% Funktion fuer QoI: Flächeninhalt des Dreiecks
f = @(x1,x2,x3) 0.5 * norm(cross([x2-x1, 0], [x3-x1,0]));

% mit exakten Erwartungswert
exa = 11/144;

% % ... oder Quadrat des Flächeninhalts
%f = @(x1,x2,x3) (0.5 * norm(cross([x2-x1, 0], [x3-x1,0])))^2;
% % ... entsprechender exakter Erwartungswert
%exa = 3/288;


%% (a) Kollokation bzw. Gauss Quadratur

% Anzahl der stochastischen Variablen
d = 6;

% ... und Orthogonalpolynome bzgl. deren Verteilung
for i=1:d, s{i} = 'legendre'; end;
   
A_Kollo = zeros(1,5);
for p = 0:4,
    % Collocation points and weights
    [X,W] = FullGridSC(d, p*ones(1,d), s);
    nX = size(X,1);
    
    A_Kollo(1+p) = 0;
    for m = 1:nX,
        % Eckpunkte des Dreiecks
        x1 = X(m, 1:2); x2 = X(m, 3:4); x3 = X(m, 5:6);

        % Flächenberechnung mittels Kreuzprodukt
        A_Kollo(1+p) =  A_Kollo(1+p) + W(m) * f(x1, x2, x3);
    end
end

% Plot
figure();
loglog((1:5).^6, abs(A_Kollo-exa), '-ok');


%% (b) Monte Carlo
%

A_MC = zeros(5,5);
% Äußere Schleife über Anzahl der MC-Simulationen
for j=1:5,
    for n = 1:5,
        M = n^6;

        % Vektor mit Flächen der zufälligen Dreiecke
        A = zeros(M,1);

        % Schleife zur Simulation und Flächenberechnung der Dreiecke
        for m = 1:M,
            % Simulation der zufälligen Eckpunkte
            x1 = rand(1,2); 
            x2 = rand(1,2); 
            x3 = rand(1,2);

            % Flächenberechnung mittels Kreuzprodukt
            A(m) =  f(x1, x2, x3);
        end

        % Monte Carlo Schätzer
        A_MC(j,n) = mean(A);
    end
end

hold on;
loglog((1:5).^6, abs(A_MC-exa)', '-+')
hold off;
legend('SC','MC 1','MC 2','MC 3','MC 4','MC 5')
xlabel('Number of Evaluations')
ylabel('Error in mean')
title('Convergence of mean estimate of MC and SC')