%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solution to CSE Your Homework Assignment Project 10
% Multidimensional Integration: Partition and Conquer
%
% Compute the area of a 2-dimensional region using 
% Monte-Carlo integration with importance sampling,
% comparing with the results of Problem 1.
% 
% See also
% Beichl and Sullivan, Computing in Science and
% Engineering, March/April 1999, p.72
%
% problem2.m Dianne P. O'Leary   09/2004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nn = 100000;
ni = 10;

radius = 0.5;
disp(sprintf('Compute the area of the quarter circle of radius %f', ...
     radius))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We generate and save nn random numbers.
% This saves us some time and makes sure that each of the
% experiments uses the same sequence.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rand('state',0);
Xu = rand(nn,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Method 3:
% Compute the Monte Carlo approximations
% (using importance sampling)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 1: Determine the distribution p(x) by sampling
%         the function we are integrating at equally-spaced points.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x = linspace(0,radius,ni+1);
h = radius/ni;   % length of interval between consecutive x's
x = x(1:ni);
f = sqrt(radius^2-x.^2);
p = f / sum(f);  % probability of choosing each interval

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 2: We choose points by importance sampling:
%         The number of samples in each interval is 
%         proportional to p.
%         We iterate for 10, 100, 1000, 10000 pts.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:5,
  n = 10^i/2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine how many points npts will be chosen in each 
% interval.  We keep this close to n*p(i), with a 
% correction to convert to integers and still have a total
% of n points. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  npts = round(n*p);
  diff = n-sum(npts);
  if (diff>0)
    npts(1:diff) = npts(1:diff)+1;
  else
    npts(1:-diff) = npts(1:-diff)-1;
  end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sum the values of sampling uniformly in each interval.
% We use npts(j) points in interval j, and their coordinates
% are stored in xpts.  The estimate is stored in answeris,
% after multiplying by the distance between limits of 
% integration and dividing by the number of samples and 
% the number of intervals.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  first = 0;
  answeris(i) = 0;
  for j=1:ni,
    xpts = (j-1)*h+Xu(first+1:first+npts(j))*h; 
    answeris(i) = answeris(i) ...
           + sum(sqrt(radius^2-xpts.^2))/p(j);
    first = first + npts(j);
  end
  answeris(i) = radius*answeris(i)/(n*ni);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 3: Display the results.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('            Monte Carlo                       Importance Sampling')
disp('      n estimate     error   error*sqrt(n) estimate     error   error*sqrt(n)')

errormc = abs(answermc - trueanswer);
erroris = abs(answeris - trueanswer);

for i=1:5,
  disp(sprintf('%7d %f %e %f     %f %e %f', ...
        10^i, answermc(i), errormc(i), errormc(i)*sqrt(10^i), ...
              answeris(i), erroris(i), erroris(i)*sqrt(10^i)))
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following code completes the figure in the paper.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

hold on
n = 20;
npts = round(n*p);
diff = n-sum(npts);
if (diff>0)
 npts(1:diff) = npts(1:diff)+1;
else
 npts(1:-diff) = npts(1:-diff)-1;
end

first = 0;
estis = 0;
xx = [];
for interval=1:ni,
  xpts = (interval-1)*h+Xu(first+1:first+npts(interval))*h;
  estis = estis + sum(sqrt(radius^2-xpts.^2))/p(interval);
  first = first + npts(interval);
  xx = [xx,xpts];
end
  estis = radius*estis/(n*ni);

figure(1)
hold on
plot(xx(1:n),zeros(n),'rx')
print -depsc mcareaest.eps
pt = Xu(n+1:2*n)*h;        % Choose points uniformly
disp(sprintf('The area estimate from the figure are'))
disp(sprintf('%f from Method 3.',estis))
print -depsc mcareaest.eps