function HW05

rng('default')
rng(0)

DRIVER_monte_carlo
DRIVER_RSVD

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function estimates the value of pi by performing a Monte Carlo
% simulation using "n" points.
%
% It draw n points at random in the box [-1,1] x [-1,1]. Let n_in denote 
% the number of points inside and outside the unit circle. As n % grows, 
% we should get
%
%    n_in / n -> pi/4
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_monte_carlo

%%% Set parameters
n_vec  = 100*(2.^(0:14));
pi_est = zeros(size(n_vec));

%%% Loop over all values of n to be tested.
for i = 1:length(n_vec)
   
  n = n_vec(i);

  %%% Draw points in the box [-1,1] x [-1,1]:
  xx = 1 - 2*rand(2,n);

  %%% Compute the distance of the points from the origin
  rsquare = sum(xx.*xx);

  %%% Find the points inside the circle.
  J_in = find(rsquare <= 1);
  %plot(xx(1,:),xx(2,:),'k.',...
  %     xx(1,J_in),xx(2,J_in),'r.')
  %axis equal

  %%% Compute the estimate for pi:
  pi_est(i) = 4*length(J_in)/n;
  fprintf(1,'n = %6d      Estimate of pi  = %14.10f\n',n,pi_est(i))
  
end

%%% Tabulate the results
fprintf(1,'      n         pi_est         error        error*sqrt(n)\n')
for i = 1:length(n_vec)
  fprintf(1,'%7d   %14.10f   %14.10f   %14.10f\n',...
          n_vec(i),pi_est(i),abs(pi_est(i) - pi),abs(pi_est(i) - pi)*sqrt(n_vec(i)))
end
keyboard

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_RSVD

%%% Set parameters
m = 300;
n = 200;
k = 40;
p = 5;
q = 1;

%%% Compute a test matrix
A = LOCAL_fast_decay(m,n,0.8*n);
%A = LOCAL_slow_decay(m,n);
ss = svd(A);

%%% Execute RSVD
[U,D,V] = LOCAL_RSVD(A,k,p,q);

fprintf(1,'Error in RSVD = %12.5e\n',norm(A - U*D*V'))
fprintf(1,'Minimal error = %12.5e\n',ss(k+1))

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = LOCAL_fast_decay(m,n,k)

% Determine the "inner dimension"
r = min([m,n,2000]);

[U,~] = qr(randn(m,r),0);
[V,~] = qr(randn(n,r),0);
alpha = (1e-15)^(1/(k-1));
ss    = alpha.^(0:(r-1));
A     = (U.*(ones(m,1)*ss))*V';

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = LOCAL_slow_decay(m,n)

% Determine the "inner dimension"
r = min([m,n]);

[U,~] = qr(randn(m,r),0);
[V,~] = qr(randn(n,r),0);
ss    = 1./(1:r);
A     = (U.*(ones(m,1)*ss))*V';

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [U,D,V] = LOCAL_RSVD(A,k,p,q)

n = size(A,2);
G = randn(n,k+p);
disp(G(1:3,1:3))
Y = A*G;
[Q,~,~] = qr(Y,0);
for i = 1:q
  [Q,~,~] = qr(A*(A'*Q),0);
end
B = Q'*A;
[Uhat,D,V] = svd(B,'econ');
U = Q*Uhat(:,1:k);
D = D(1:k,1:k);
V = V(:,1:k);

return