% Calling program for main.m - to solve   % revising nov20/2014
%         approx closest Euclidean distance matrix problem.
%  Generates random matrices A,H symmetric and >= 0 elementwise.
%   Solves  min ||H.*(D-A)|| s.t. D is a distance matrix
%
% USER:  Can change the following parameters below:
%   n  -  dimension of matrices H and A
%   spH  -  percentage nonzeros in H and A (excluding diagonal)
%   noprobs  - number of problems to solve
%   toler    - tolerance value for relative duality gap, i.e. stopping
%                   criteria for interior-point method
%
%  OUTPUT: Output is saved, after each problem, in data.mat
%            X      -  primal optimum matrix
%            Lam    -  dual optimum matrix
%             D     -  optimal distance matrix
%          itersum  -  iteration count for all problems - saved
%          optvals   -  optimal values for all problems  - saved
%%%%%%%%%%%%%%%%%%%%%%%
clear
%warning off    % in matlab 5 only
n=48;   % dimensions of data distance matrix A 
spH=.3;   % sparsity in graph  (done using   H.*A)
noprobs=1;   % number of problems to solve
toler=1e-11;   % e-11 is a very 'small' tolerance
kind=1;     %   1 indicates random A generated
            %   2 indicates that A is a distance matrix so optval=0
            %      in this case H = ones.
            %   2 indicates that A has banded structure
%%%%%%%%%%%%%%%%%%%%%%%
%Initialization:
%rand('seed',sum(100*clock));  
%randn('seed',sum(100*clock));
%% comment out next two lines to repeat the same 'random' problem
rng('shuffle')
%rng('default');     % for same problem
saverng=rng;
rng(saverng);
itersum=[];
optvals=[];
rankXs=[];
%%%%%%%%%%%%%%%%%%%%%%%
%Generate and solve noprobs problems
for ii=1:noprobs,
   % for fixed seed, use next 3 lines
   % iseed=ii;             
   % randn('seed',iseed);
   % rand('seed',iseed);
  disp(['generating test problem number and dim : ' ...
   ' prob no. ',num2str(ii),'  of  ',num2str(noprobs), ' and dim=',num2str(n)]);
  tic;
  if kind==1,   % random problem
    A=sprandsym(n,spH);
    A=abs(A);
    A=sparse(A-diag(diag(A)));
    A=.5*sparse(A+A');
    if spH ==1,  % dense case
         H=randn(n);
       else
         H=sprandsym(n,spH);
    end
    H=abs(H);   % nonnegative weights
    H=sparse(H+H')/2;
  %%%%%%%%%%%%beginning of special structure problems
  elseif kind==2,   % A is a distance matrix
    H=ones(n);
    kk=floor(n/2);
    x=rand(n,kk);
    A=K(x*x');
    A=.5*(A+A');
%   A=ones(n)-eye(n);  % to generate the simplex 
  elseif kind==3  % band structure for A
    A=zeros(n) ;
    for iii=1:n-1, 
     A(iii,iii+1)= 1+floor(5*rand(1)) ;
    end ;
    A=A+A' ;
    for iii=1:n-1, 
     H(iii,iii+1)= 10 ;
    end ;
    H(1,n)= 10 ;
    H=H+H' ;
  end
  %%%%%%%%%%%%end of special structure problems
  
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
  save probdata      % save for debugging
  fprintf('problem saved in   probdata.mat \n');
  fprintf('cpu time for problem generation is:  %g  \n',toc);
profile on    % in matlab 5 only
     [optval,Lam,oX,X,oD,D,V,iter]=main(A,H,toler);  % solve the problem
                       % using p-d i-p method
profile report
  rankX=rank(full(X),.0001);
% probi=num2str(ii);
  dimi=num2str(n);
  spHi=num2str(spH);
  toleri=num2str(toler);
  embdim=num2str(rankX);
  title(['dim ' dimi ', H density ' spHi ', dual. gap toler. ' toleri ', rankX ' embdim ]);
% probitext=['Problem ' probi ' with dimension ' dimi ' and density ' ...
%         spHi ' and tolerance ' toleri];
  hold off
  figure(1)
  itersum=[itersum iter];
  optvals=[optvals optval];
  rankXs=[rankXs rankX];
%'at end of problem in gendist'
%keyboard
% save data
% save tempdata n spH noprobs toler itersum optvals rankXs kind
% emerg
% clear
% load tempdata
end
