function dmc = demo_MC()
%
% Integrace 2D funkce metodou Monte Carlo
% Integrujeme na intervalu <0,1> x <0,1>
%

test = 1;

% integrovana funkce a presne reseni (jen pro vykresleni)
if(test==1)
    f = inline('abs(sin(pi^2.*((x-0.5).^2+(y-0.5).^2)))','x','y');
    presny = 0.655172456; % by Matlab. By Maple: 0.655172424
elseif(test==2)
    f = inline('x+y','x','y');
    presny = 1;  
end;

%% pocet kroku
kroku = 10000;


% pole pro ukladani vysledku
vysl = zeros(kroku,1);

tic % zacni merit cas
z = 0;
for i=1:kroku,
  r       = rand(1,2);    % dve nahodna cisla v (0,1)
  value   = f(r(1),r(2)); % vypocteme funkci ve dvou nahodnych bodech
  z       = z + value;    % pricteme ke kumulovane hodnote
  vysl(i) = z / i;        % prumer (kumul. hodnota delena poctem kroku)
end;
integral = z / kroku;
cas = toc; % skonci mereni casu


%% spocitej relativni chybu
err_rel = abs(integral-presny)/max(abs(presny),1e-30);


%% vypis vysledek
disp(['Vypocteny integral je ' num2str(integral) ...
      ',   rel.chyba ' num2str(err_rel)]);
disp(['Vypocet trval ' num2str(cas) ' sekund']);


%% nakresli prubeh a chyby

x    = 1:kroku;
xlog = log10(1:kroku);
chyb = abs(vysl-presny);
teor = x.^(-1/2);

fig=figure;

%  vypoctena hodnota integralu a presne reseni
subplot(1,3,1)
plot( x,vysl,'b', [1 kroku],[presny presny],'k--' )
xlabel('K')
ylabel('Int     ,     exact')

%  abs. chyba a krivka k^(-1/2) 
subplot(1,3,2)
plot( x,chyb,'b', x,teor ,'k--');
xlabel('K')
ylabel('|Err|     ,     K^{-1/2}         (detail)')
%  ukaz detail
yl = get(gca,'yLim');
set(gca,'yLim',[yl(1) 0.2*yl(2)]);

%  to same logaritmicky
subplot(1,3,3)
plot( xlog,log10(chyb),'b', xlog,log10(teor),'k--' );
xlabel('log_{10}(K)')
ylabel('log_{10}(|Err|)     ,     log_{10}(K^{-1/2})')
