prodcasctest(TOL)

function prodcasctest(TOL)
%Performs numerous tests of ProdCascade math operations, comparing the results
%to those from equivalent numerical forms.
%
%  prodcasctest(TOL)
%
%TOL is a tolerance value on the percent error (default=inf). 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

%%some useful functions
 ncols=@(M) size(M,2);     
 err=@(x,y) DiscrepancyMeasure(x,y,TOL);

%%Some operands

a=rand(4); a=a*a'+eye(4); %an invertible, positive definite matrix
b=rand(4,5); 
c=rand(5,4)*i;
   
%%Some ProdCascades formed from the above operands.
A=ProdCascade({a,b,c});
S=ProdCascade({a.',a});


%%Their conventional forms for comparison purposes
Afull=a*b*c;
Sfull=a.'*a;



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


%%Test of full()
Error(1)=     err( full(A), Afull  );  
Error(end+1)= err( full(S), Sfull  );  
 

%%Test of sparse()
Error(end+1)= err( sparse(A), sparse(Afull) );
Error(end+1)= err( issparse(sparse(A)), 1 );


   %%No sense in proceding if the tests so far didn't pass - the error
   %%calculations rely on the functionality of full@ProdCascade(A)
    if max(Error)>TOL,
        Error,
        error 'Something wrong with sparse() and full() methods'; 
    end

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


%%Test of size()
Error(end+1)= err( size(A) ,  size(Afull) );


%%Test of inv

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

%%Test of det

Error(end+1)= err( det(S)  , det(Sfull) );
   
%%Test of mtimes

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

Error(end+1)= err( S*A , Sfull*Afull  ); %prod of 2 ProdCascades

    X=rand(size(Afull,2),size(Afull,1));

Error(end+1)= err( A*X , Afull*X  ); %prod with ordinary matrix
Error(end+1)= err( X*A , X*Afull  );




%%Test of horzcat

Error(end+1)= err( [A,A] , [Afull*Afull]  ); 

%%Test of subsref,subsasgn, end

                C=A;
               
Error(end+1)= err( C{end} , c  ); 
Error(end+1)= err( C(end-1:end) , b*c  ); 

              d=rand(size(c));
              C{end}=d;
              
Error(end+1)= err( C, a*b*d  ); 
           



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


MAX_ERROR=max(Error);

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





function errval=DiscrepancyMeasure(X,Y,TOL)

  errval=Discrepancy(X,Y)/Discrepancy(0,Y)*100; %normalize

  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)]);