%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solution to CSE Your Homework Assignment Project 
% Multidimensional Integration: Partition and Conquer
%
% Compute the area of a 2-dimensional region using 
% Monte-Carlo integration with uniformly distributed random 
% numbers.
% In this example, we estimate the area of a quarter circle 
% of radius 0.5, bounded by the positive x and y axes and 
% the curve sqrt(radius^2 - x^2).
% We compute two estimates of the area:
% Method 1:
%   Test whether random points in the bounding square of the
%   quarter circle fall inside the region or outside.  The
%   estimate is
%      the area of the square 
%      times the number of samples inside the area 
%      divided by the total number of samples.
% Method 2:
%   The estimate is
%    the average value of sqrt(radius^2 - x^2) for x in [0,radius]
%    times the length of the interval [0,radius]
%
% challenge1.m Dianne P. O'Leary   09/2004                
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nn = 100000;  

close all
radius = 0.5;
trueanswer = radius^2*pi/4;

disp(sprintf('Estimate the area of the quarter circle of radius %f', ...
     radius))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We generate and save nn random pairs, after 
% reinitializing the random number generator.
% This saves us some time and makes sure that each of the
% experiments uses the same sequence.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rand('state',0);
Xu = radius* rand(nn,2);
disp(' ')
disp(sprintf('True area: %f5.8',trueanswer));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code produces the figure in the paper.
% Plot the region to be integrated and the first 20 pairs 
% of numbers to illustrate how the methods work.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)
hold on
plot(Xu(1:20,1),Xu(1:20,2),'g*')
fvals = sqrt(radius^2 - Xu(1:20,1).^2);
plot(Xu(1:20,1),fvals,'bo')
xx = 0:.01:radius;
xx = [xx,radius];
axis([0 radius 0 radius])
axis square
yy = sqrt(radius^2 - xx .* xx);
plot(xx,yy,'r')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Method 1: find the fraction of points that lie inside 
%           the region, and approximate the integral by 
%           multiplying the fraction times the area of the 
%           square [0, radius] x [0, radius].
% Compute these Monte Carlo approximations
% for 10, 100, 1000, 10000, 100000 pts and compare
% with the true answer.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(sprintf('Method 1: inside or outside'));
inside = (Xu(:,1).^2 + Xu(:,2).^2 < radius^2);
disp('         n   estimate        error   error*sqrt(n)')
for i=1:5,
  n = 10^i;
  area = radius^2*sum(inside(1:n))/n;
  error = area - trueanswer;
  disp( ...
   sprintf('%10.0f   %f   %e   %e', ...
        n, area, error, error*sqrt(n)) )
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Method 2.  Approximate the integral by the average 
% function value sqrt(r^2 - x^2) times the length of the 
% interval [0,r].
% Compute these Monte Carlo approximations
% for 10, 100, 1000, 10000, 100000 pseudorandom points.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(' ')
disp(sprintf('True area: %f5.8',trueanswer));
disp(sprintf('Method 2: average value'));
fvals = sqrt(radius^2 - Xu(:,1).^2);
disp('         n   estimate        error   error*sqrt(n)')
for i=1:5,
  n = 10^i;
  answermc(i) = radius * sum(fvals(1:n))/n;
  error = answermc(i) - trueanswer;
  disp( ...
   sprintf('%10.0f   %f   %e   %e', ...
        n, answermc(i), error, error*sqrt(n)) )
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compare with Matlab's quad
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[quadest,nf] = quad('sqrt(.5^2 -x.^2)',0,.5,1.e-4);
disp(' ')
disp(sprintf('Function quad gives an error of '))
disp(sprintf('%f with %d function evaluations',  ...
     quadest-trueanswer, nf))


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code is for the figure caption.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n=20;
est1 = radius^2*sum(inside(1:n))/n;
est2 = radius * sum(fvals(1:n))/n;
disp(' ')
disp(sprintf('The area estimates from the figure are'))
disp(sprintf('%f from Method 1 and %f from Method 2.',est1,est2))
print -depsc mcareaest.eps
