Randomized Algorithms Control Toolbox Code

function Result = utest3(hObject, varargin)
%UTEST3 3-steps wrapper for calling uncertain property testing
%
%  ObjectSamples = utest3(hObject, nSamples)
%  ObjectSamples = utest3(hObject, nSamples, VecObj)
%  CritSamples = utest3(hObject, hCrit, nSamples)
%  CritSamples = utest3(hObject, hCrit, nSamples, VecObj)
%  CritSamples = utest3(hObject, hCrit, nSamples, VecObj, VecCrit)
%  OverallResult = utest3(hObject, hCrit, hSecCrit, nSamples)
%  OverallResult = utest3(hObject, hCrit, hSecCrit, nSamples, VecObj, VecCrit)
%
%  Gathers object samples by hObject function, 
%  for each sample calculates value of criterion (by hCrit function)
%  finally checks all the criterion values with over-all function (by hSecCrit)
%  (the last one used for getting worst-case criterion (performance), for example)
%
%  instead of function handle hCrit there can be used built-in function names
%    'ident' - for case when sample of object returned is already value of criterion
%    'hurwitz' - check if sample returned by hObject is a Hurwitz polynomial
%    'schur' - check if sample returned by hObject is a Schur polynomial
%    'posdef' - check if sample returned is a positive definite matrix
%    'possemidef' - check if sample returned is a positive semi-definite matrix
%  
%  instead of function handle hSecCrit can be used built-in function names
%    'ident' - for case when sample of object returned is already value of criterion
%    'maxval' - for maximum value of criterion values (worst-case probability analysis)
%    'meanval' - for mean criterion values (probability calculation)
%    'nonzero' - for number of non-zero criterion values
% 
%  nSamples is the number of samples taken
%
%  VecXXX are flags letting know if corresponding funcion can accept cell array of samples
%
%   See also UBASE, ROOTS.

% $Id$

% parameter check
if ~any([nargin == 2, nargin == 3, nargin == 4, nargin == 5, nargin == 6])
  error(['RACT:utest3:InputSize', ['Wrong number of parameters: ' num2str(nargin)]]);
elseif ~isa(hObject, 'function_handle')
  error(['RACT:utest3:InputType', 'First parameter must be function handle!']);
end

% fill undefined parameters
switch nargin 
  case 2 % (hObject, nSamples)
  nSamples = varargin{1}; % second parameter
  VecObj = -1; % default value - autodetect via try...catch
  hCrit = 'ident';  % default value
  VecCrit = -1; % default autodetect value
  hSecCrit = 'ident'; % default value 
  case 3 % (hObject, hCrit, nSamples) or (hObject, nSamples, VecObj)
    if or(isa(varargin{1}, 'function_handle'), isa(varargin{1}, 'char')) % (hObject, hCrit, nSamples)
      nSamples = varargin{2}; % third parameter
      VecObj = -1; % default value - autodetect via try...catch
      hCrit = varargin{1};  % second parameter
      VecCrit = 1; % default value
      hSecCrit = 'ident'; % default value 
    elseif isa(varargin{1}, 'numeric') % (hObject, nSamples, VecObj)
      nSamples = varargin{1}; % second parameter
      VecObj = varargin{2};  % third parameter
      hCrit = 'ident';  % default value
      VecCrit = -1; % default autodetect value
      hSecCrit = 'ident'; % default value 
    end
  case 4 % (hObject, hCrit, hSecCrit, nSamples) or (hObject, hCrit, nSamples, VecObj)
    if or(isa(varargin{2}, 'function_handle'), isa(varargin{2}, 'char')) % (hObject, hCrit, hSecCrit, nSamples)
      nSamples = varargin{3}; % fourth parameter
      VecObj = -1; % default autodetect value
      hCrit = varargin{1};  % second parameter
      VecCrit = -1; % default autodetect value
      hSecCrit = varargin{2}; % third parameter 
    elseif isa(varargin{2}, 'numeric') % (hObject, hCrit, nSamples, VecObj)
      nSamples = varargin{2}; % third parameter
      VecObj = varargin{3}; % fourth parameter
      hCrit = varargin{1};  % second parameter
      VecCrit = -1; % default autodetect value
      hSecCrit = 'ident'; % default value 
    end
  case 5 % (hObject, hCrit, nSamples, VecObj, VecCrit) 
    nSamples = varargin{2}; % third parameter
    VecObj = varargin{3}; % fourth parameter
    hCrit = varargin{1};  % second parameter
    VecCrit = varargin{4}; % fifth parameter
    hSecCrit = 'ident'; % default value 
  otherwise % (hObject, hCrit, hSecCrit, nSamples, VecObj, VecCrit) - full set
    nSamples = varargin{3}; % forth parameter
    VecObj = varargin{4}; % fifth parameter
    hCrit = varargin{1};  % second parameter
    VecCrit = varargin{5}; % sixth parameter
    hSecCrit = varargin{2}; % third parameter 
end

if ~isa(hCrit, 'function_handle')
  if isa(hCrit, 'char') % check if second parameter a name of built-in criterion
    switch hCrit
      case 'ident'
       hCrit = @ident;
       VecCrit = 1; 
      case 'hurwitz'
        hCrit = @ishurwitz;
        VecCrit = 1; 
      case 'schur'
        hCrit = @isschur;
        VecCrit = 1; 
      case 'posdef'
        hCrit = @isposdef;
        VecCrit = 1; 
      case 'possemidef'
        hCrit = @ispossemidef;
        VecCrit = 1; 
      otherwise
        error(['RACT:utest3:InputType', ['Unknown built-in criterion name: ' hCrit]]);
    end
  else
    error(['RACT:utest3:InputType', 'Second parameter must be function handle or built-in criterion name!']);
  end
end

if ~isa(hSecCrit, 'function_handle')
  if isa(hSecCrit, 'char') % check if second parameter a name of built-in criterion
    switch hSecCrit
      case 'ident'
        hSecCrit = @ident;
      case 'maxval'
        hSecCrit = @maxval; % use internal function
      case 'meanval'
        hSecCrit = @meanval; % use internal function
      case 'nonzero'
        hSecCrit = @nonzero;
      otherwise
        error(['RACT:utest3:InputType', ['Unknown built-in criterion name: ' hSecCrit]]);
    end
  else
    error(['RACT:utest3:InputType', 'Third parameter must be function handle or built-in criterion name!']);
  end
end

if ~(nSamples > 0)
  error(['RACT:utest3:InputType', ['Invalid number of samples: ' num2str(nSamples)]]);
end

% get samples of user-defined objects as array of cells
uObjects = uevaluate(hObject, nSamples, VecObj, 'cell');

% vectorization for criterion
if VecCrit == 1 
  CritVals = feval(hCrit, uObjects);
elseif VecCrit == 0
  CritVals = cell(size(uObjects));
  for k = 1:nSamples
    CritVals{k} = feval(hCrit, uObjects{k});
  end
elseif VecCrit == -1 % autodetect
  try
    CritVals = feval(hCrit, uObjects);
  catch
    CritVals = cell(size(uObjects));
    for k = 1:nSamples
      CritVals{k} = feval(hCrit, uObjects{k});
    end
  end
end

if isa(hSecCrit, 'char')
  if strcmp(hSecCrit, 'ident')
    Result = CritVals;
    return;
  end
end

% SecCrit should allow input, packed in cells
Result = feval(hSecCrit, CritVals);

return;
%% END OF THE FUNCTION

% Some built-in functions:
% 1 - 'ident'
% just return the argument
function CritVals = ident(Var)
CritVals = Var;

% 2 - 'hurwitz'
% check if polynomial(s) is Hurwitz
function CritVals = ishurwitz(Var)
CritVals = cell(size(Var));
N0 = max(size(Var));
N = length(Var{1}); % N = n + 1
N2 = ceil(N/2); % (n+1)/2 for odd degree, n/2 + 1 for even
N_add = N-2*N2; % for even/odd degree diff: = -1 if N is odd
for n=1:N0
  p = Var{n}; % should be row to the Routh algorihtm
	% make first coefficient positive
	if p(1) < 0 
    p = -p;
	end
	% neccessary condition check
	if ~all(p > 0)
    CritVals{n} = 0;
    continue;
	end
	% full table 
	T = zeros(N, N2);
	T(1, 1:N2) = p(1:2:end);
	T(2, 1:N2+N_add) = p(2:2:end+N_add); % N-2*N2 
	for k=3:N
    if T(k-1, 1) <= 0 % a zero or negative value is found in first line
      CritVals{n} = 0; 
      continue;
    else
     r = T(k-2, 1)/T(k-1, 1);
    end
    for m = 1:N2-1
      T(k, m) = T(k-2, m+1) - r*T(k-1, m+1);
    end
	end
  CritVals{n} = (T(N, 1) > 0); % other elems are already checked
%   if max(real(roots(Var{k}))) < 0 - old straightforward realization
%     CritVals{k} = 1;
%   else
%     CritVals{k} = 0;
%   end
end

% 3 - 'schur'
function CritVals = isschur(Var)
CritVals = cell(size(Var));
N0 = max(size(Var));
N = size(Var{1}, 2);
for k=1:N0
  CritVals{k} = 1;
  p = Var{k}; % should be row to the Routh algorihtm
  T1 = p;
  T2 = conj(fliplr(p)); % z^n conj(P(1/z))
  for n=N:-1:2
    if abs(T1(end)) >= abs(T1(1));
      CritVals{k} = 0;
      break;
    else
      r = T1(end)/conj(T1(1));
      T1 = T1 - r*T2;
      T1 = T1(1:end-1); % free element should be zero 
      T2 = conj(fliplr(T1)); % z^n conj(P(1/z))
    end
  end
%  if max(abs(roots(Var{k}))) < 1 - old straightforward realization
%    CritVals{k} = 1;
%  else
%    CritVals{k} = 0;
%  end
end

% 4 - 'posdef'
function CritVals = isposdef(cellA)
% check if matrix A is positive definite
Eps1 = 1e-8;
Eps2 = 1e-7;
if size(cellA{1},1) ~= size(cellA{1},2)
  error(['RACT:utest3:posdef:InputType', 'The argument should be a square matrix cell array.']);
end
CritVals = cell(size(cellA));
N = max(size(cellA));
for k=1:N
  A = cellA{k};
  if norm(A - A', 'fro') > Eps1
    warning('RACT:utest3:posdef:InputType', 'The matrix A should be symmetric (or ... for complex?AT? ).');
  end
  T = schur(A); % A = U*T*U.' decomposition
  t = diag(T);
  if all(t >= Eps2)
    CritVals{k} = 1;
  else
    CritVals{k} = 0;
  end
end

% 5 - 'possemidef'
function CritVals = ispossemidef(cellA)
% check if matrix A is positive semi-definite
Eps1 = 1e-8;
if size(cellA{1},1) ~= size(cellA{1},2)
  error(['RACT:utest3:possemidef:InputSize', 'The argument should be a square matrix cell array.']);
end
CritVals = cell(size(cellA));
N = max(size(cellA));
for k=1:N
  A = cellA{k};
  if norm(A - A', 'fro') > Eps1
    warning('RACT:utest3:possemidef:InputSize', 'The matrix A should be symmetric (or ... for complex?AT? ).');
  end
  T = schur(A); % A = U*T*U.' decomposition
  t = diag(T);
  if all(t >= 0)
    CritVals{k} = 1;
  else
    CritVals{k} = 0;
  end
end

% auxilary - if numeric content of cell array is a-priori known, 
% convert to numeric matrix. No error check.
function NumArray = cell2num(CellArray)
NumArray = zeros(size(CellArray));
for k=1:max(size(NumArray))
  NumArray(k) = CellArray{k};
end

% 6 - 'maxval'
function MaxVal = maxval(Var)
MaxVal = max(cell2num(Var));

% 7 - 'meanval'
function Mean = meanval(Var)
Mean = mean(cell2num(Var));

% 8 - 'nonzero'
function NZ = nonzero(Var)
NZ = sum(cell2num(Var) ~= 0);