rstest(TOL)

function rstest(TOL)
%Performs numerous tests of RobustSparse math operations, 
%
%  rstest(TOL)
%
%TOL is a tolerance value on the percent error. Execution will pause in debug
%mode for inspection if any one of the tests exhibits an error greater than
%TOL.

if nargin<1
 TOL=inf; %default tolerance value on discrepancies
end

CHECKTYPES=false;

%%function for measuring error

 err=@(x,y) DiscrepancyMeasure(x,y,TOL,CHECKTYPES);

 
Asp=sparse(rand(3,2));  
Bsp=sparse(rand(3,2));
Ssp=sparse(rand(3));
Vsp=sparse(rand(3,1));


%%Representations of the above as RobustSparse objects 
A=rsparse(Asp);   % =Asp
B=rsparse(Bsp);   % =Bsp
S=rsparse(Ssp);   % =Ssp
V=rsparse(Vsp);   % =Vsp


%%%%%%%%%%%%%%%%%%%%%%%%%%%TESTS%%%%%%%%%%%%%%%%%%%%%%%%%%


 
%%Test of transpose, ctranspose
Error(1)= err( A.' ,  Asp.'   );
Error(end+1)= err( A'  ,  Asp'   );

%%Test of uplus, uminus

Error(end+1)= err( +S  , +Ssp );
Error(end+1)= err( -S  , -Ssp );

%%Test of plus, minus

Error(end+1)= err( A+B  , Asp+Bsp );
Error(end+1)= err( A-B ,  Asp-Bsp );


%%Test of inv

Error(end+1)= err( inv(S)  , inv(Ssp) );


%%Test of inv

Error(end+1)= err( inv(S)  , inv(Ssp) );


%%Test of sum

Error(end+1)= err( sum(A,1)  , sum(Asp,1) );
Error(end+1)= err( sum(B,2)  , sum(Bsp,2) );
Error(end+1)= err( sum(S)    , sum(Ssp) );


%%Test of mtimes

Error(end+1)= err( A*3 , Asp*3  );  %scalar with RobustSparse
Error(end+1)= err( 3*A , 3*Asp  );

Error(end+1)= err( S*A , Ssp*Asp  ); %prod of 2 RobustSparses

   x=Asp(end,:).';
   y=Asp(:,end).';

Error(end+1)= err( A*[x,x] , Asp*[x,x]  ); %pre-mult with columnized data
Error(end+1)= err( [y;y]*A , [y;y]*Asp  ); %post-mult with columnized data




%%Test of mldivide
Error(end+1)= err( 3\B , 3\Bsp  );  %scalar with RobustSparse

Error(end+1)= err( B\A , Bsp\Asp  ); %prod of 2 RobustSparses



%%Test of mrdivide
  
    Bt=B.'; Btsp=Bsp.';

Error(end+1)= err( Bt/3 , Btsp/3  );%scalar with RobustSparse

Error(end+1)= err( A.'/Bt , Asp.'/Btsp ); %prod of 2 RobustSparses



%%Test of times

Error(end+1)= err( A.*3 , Asp.*3  );  %scalar with RobustSparse
Error(end+1)= err( 3.*A , 3.*Asp  );

Error(end+1)= err( A.*B , Asp.*Bsp  ); %prod of 2 RobustSparses


%%Test of rdivide

Error(end+1)= err( B./3 , Bsp./3  );  %scalar with RobustSparse
Error(end+1)= err( 3./B , 3./Bsp  );

Error(end+1)= err( A./B , Asp./Bsp  ); %2 RobustSparses
Error(end+1)= err( B./A , Bsp./Asp  ); %2 RobustSparses

%%Test of ldivide

Error(end+1)= err( B.\3 , Bsp.\3  );  %scalar with RobustSparse
Error(end+1)= err( 3.\B , 3.\Bsp  );

Error(end+1)= err( B.\A , Bsp.\Asp  ); %2 RobustSparses
Error(end+1)= err( A.\B , Asp.\Bsp  ); %2 RobustSparses

%%Test of power, mpower
Error(end+1)= err( A.^2 ,  Asp.^2   );
Error(end+1)= err( S^2  ,  Ssp^2   );


%%Test of relops

Error(end+1)= err( A>A ,  Asp>Asp   );
Error(end+1)= err( A>B  , Asp>Bsp   );

Error(end+1)= err( A>=A ,  Asp>=Asp   );
Error(end+1)= err( A>=B  , Asp>=Bsp   );

Error(end+1)= err( A<A ,  Asp<Asp   );
Error(end+1)= err( A<B  , Asp<Bsp   );

Error(end+1)= err( A<=A ,  Asp<=Asp   );
Error(end+1)= err( A<=B  , Asp<=Bsp   );


Error(end+1)= err( A==A ,  Asp==Asp   );
Error(end+1)= err( A==B  , Asp==Bsp   );

Error(end+1)= err( A~=A ,  Asp~=Asp   );
Error(end+1)= err( A~=B  , Asp~=Bsp   );

%%Test of logical ops

Error(end+1)= err( A&A ,  Asp&Asp   );
Error(end+1)= err( A|B  , Asp|Bsp   );

Error(end+1)= err( ~A ,   ~Asp   );
Error(end+1)= err( ~B  ,  ~Bsp   );

%%Test of spfun

 f=@(x) cos(x).^2;
 
Error(end+1)= err( spfun(f,A) ,    spfun(f,Asp)  );


%Test of logical indexing
Error(end+1)= err( A(A<.5) ,  Asp(Asp<.5)   );

%Test of linear indexing
Error(end+1)= err( A(:) ,  Asp(:)   );
Error(end+1)= err( A(1:4) ,  Asp(1:4)   );
Error(end+1)= err( V(1:2) ,  Vsp(1:2)   );

%Test of subscript indexing
Error(end+1)= err( A(2,1) ,  Asp(2,1)   );
Error(end+1)= err( A(2,:) ,  Asp(2,:)   );
Error(end+1)= err( A(:,1) ,  Asp(:,1)   );

%%%%%%%%%%%%%%%%%%%%%%%END OF TESTS%%%%%%%%%%%%%%%%%%%%%%%%%%


MAX_ERROR=max(Error);

disp(['Maximum observed error was   ' num2str(MAX_ERROR) ' percent.'])





function errval=DiscrepancyMeasure(X,Y,TOL,CHECKTYPES)

  errval=Discrepancy(X,Y)/Discrepancy(0,Y)*100; %normalize
  
  isrobust=@(c) ~isempty(strfind(class(c),'Robust'));
  
  if CHECKTYPES
    if ~isrobust(X),
       warning(['X is not Robust class, but rather class ' class(X)]) 
    end
  end

  if errval>TOL, 
      disp ' '; disp 'Discrepancy detected'
      errval, 
      x=full(X); y=full(Y);
      keyboard;
  end 
  
  
function errval=Discrepancy(X,Y) 
%Primary error measurement function

  fin=@(a) a(isfinite(a));
  nonfin=@(a) a(~isfinite(a)); 

  x=full(X); y=full(Y);

  errval= norm( fin(x-y) , inf)+...   
       ~isequalwithequalnans(nonfin(x),nonfin(y))*...
       ~isempty([nonfin(x);nonfin(y)]);