% Yi-Kai Liu
% 5/10/01

% N * N real symmetric matrices, SPARSE: 
% elements are +1 or -1 with small probability p, 
% and 0 otherwise

% Implemented using regular matrix representation

% OK: Mehta normalizes by sqrt(2*N), but I use 2*sqrt(N) ??? 
% Apparently this is due to some notation trouble. In Porter's 
% book, p.197, the semicircle is normalized by 2*sqrt(N). 

% DONE: Exporting data
% FIXME: Examine debugging output

% FIXME: Criteria for using sparse storage
% NOTE: The sparse matrix functions appear to save space 
% at the expense of running time. At these sizes (up to 
% matrices of size 800*800), this hurts the performance; 
% see "readme-sparse-mtx.txt". Also note that when using 
% sparse matrices, it is better to allocate a fresh matrix 
% instead of trying to write over an already-used matrix. 
% Whenever possible, avoid iterating over the entries in 
% a sparse matrix (the data structure is not designed for 
% this kind of access). 

% TO DO: Review this code!
% TO DO: Other sparse matrices, e.g., band-diagonal

distribution_name = 'Sparse';

N = input('N * N matrices, choose N: ');
num_matrices = input('Number of matrices: ');

p = input('p = probability of +1 or -1, choose p: ');
distribution_name = [distribution_name, ' (p=', num2str(p), ')'];

save_to_disk = input('Save output to disk? (1=yes, 0=no) ');
if save_to_disk == 1
   base_filename = input('Filename (without suffix): ', 's');
end;

d = zeros(N, num_matrices);     % Eigenvalues
sp = zeros(N-1, num_matrices);  % Spacings


for m = 1 : num_matrices
  
  A = zeros(N);   % Get a new, all-zero matrix before each iteration

  % Generate random matrix
  % for i = 1 : N
  %   for j = i : N
  %     % This entry has probability p of being nonzero, 
  %     % in which case it is equally likely to be +1 or -1
  %     if rand < p
  %       if rand > 0.5
  %         A(i,j) = 1.0;
  %         A(j,i) = 1.0; 
  %       else
  %         A(i,j) = -1.0;
  %         A(j,i) = -1.0;
  %       end
  %     end
  %   end
  % end
  
  % Generate random matrix--vectorized
  for i = 1 : N
    for j = i : N
      % This entry has probability p of being nonzero, 
      % in which case it is equally likely to be +1 or -1
      if rand < p
        if rand > 0.5
          A(i,j) = 1.0;
        else
          A(i,j) = -1.0;
        end
      end
    end
  end
  A = A + (triu(A,1))';

  % Normalized eigenvalues in interval [-1,1]
  % d(:,m) = eig(A)./(1.0*sqrt(N));
  % CHANGED TO: 
  if p==0.05
     fudge_factor = 0.5;
  elseif p==0.15
     fudge_factor = 0.8;
  elseif p==0.25
     fudge_factor = 1.0;
  elseif p==0.35
     fudge_factor = 1.2;
  else
     fudge_factor = 1.0;
  end;
  d(:,m) = eig(A)./(fudge_factor*sqrt(N));

  sd = sort(d(:,m));   % Sorted eigenvalues

  % Compute spacings
  % for k = 1 : N-1
  %   sp(k,m) = sd(k+1) - sd(k);
  % end
  
  % Compute spacings--vectorized
  sp(:,m) = sd(2:N) - sd(1:N-1);

  sp(:,m) = sp(:,m) ./ mean(sp(:,m));   % Normalized spacings
                                        % with mean=1
  
end


% Save variables to disk
if save_to_disk == 1
   save([base_filename, '-data'],... 
        'distribution_name', 'save_to_disk', 'base_filename',... 
        'N', 'num_matrices', 'd', 'sp');
end;

% Plot distributions
plotdistrib;