%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solution to CS&E Your Homework Assignment Project 10
% Multidimensional Integration: Partition and Conquer
%
% Compute (poor) estimates of the partition function
% for the Harmonic oscillator potential for
% myd = 1,2,3,8, and 16 particles.
%
% problem5.m Dianne P. O'Leary   09/2004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Global variables:
%  myk = Planck's constant
%  myd = number of particles
%  mydelta = L/(myd+1)
%  zsamples = random or quasirandom samples
%  mynsamples = number of samples
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global myk myd mydelta mybnd
global zsamples mynsamples

mybnd = 3;
myk = 1.38;   % angstroms^2  g /sec^2 K
L = 1;
sz = [100 1000 10000 100000];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Repeat the computation for random and quasi-random
% sampling.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for j=1:2,
   if (j==1)
     disp('Generating quasirandom numbers.')
     zsamples = quasirand(1,sz(end),16);
   else
     disp('Generating random numbers.')
     zsamples = rand(sz(end),16);
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Repeat the computation for various numbers  n2  of random
% samples.
% Use Matlab's quad to evaluate the integral with respect
% to a, providing function values to it using the 
% Monte Carlo estimates produced in quadmcqr.
%
% est(i,j) is the estimate of the 2^{i-1} particle function
% using sz(j) points in each function evaluation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   kk = 0;
   for n2=sz
      disp(sprintf('Computing estimates using %d points.',n2))
      kk = kk + 1;
      for dexp=0:4
        myd = 2^dexp;
        mydelta = L/(myd+1);
        d = myd;
        mynsamples = n2;
        est(dexp+1,kk) = quad(@quadmcqr,-mybnd,mybnd);
      end % for dexp
   end % for n2

   if (j==1)
      estq = est;
   else
      estr = est;
   end
   clear est

end % for j

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display the results.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(' ')
disp(sprintf('Estimates of Partition Function Using Quasirandom Samples'));
disp(' ')
disp('         n        d=1         d=2         d=4         d=8         d=16')
for i=1:4,
  disp( ...
   sprintf('%10.0f  %10.2e  %10.2e  %10.2e  %10.2e  %10.2e', ...
        sz(i), estr(:,i)'))
end

disp(' ')
disp(sprintf('Estimates of Partition Function Using Quasirandom Samples'));
disp(' ')
disp('         n        d=1         d=2         d=4         d=8         d=16')
for i=1:4,
  disp( ...
   sprintf('%10.0f  %10.2e  %10.2e  %10.2e  %10.2e  %10.2e', ...
        sz(i), estq(:,i)'))
end