function result = testpoly(P,x,tol)
% The instruction
%
%   TESTPOLY(P,X)
%
% evaluates the vector X of polynomial indeterminates on the cell array P
% of polynomial constraint coefficients as used in GLOPTIPOLY
%
% Some tolerance TOL can be specified as an optional input argument
% Warning messages are displayed when either
% . a non-negative inequality is less than -TOL, or
% . a non-positive inequality is greater than +TOL, or
% . the absolute value of an equality is greater than +TOL
% when evaluated at vector X
% By default TOL is set to 1e-6
%
% With the syntax
%
%  RESULT = TESTPOLY(P,X)
%
% the screen output is turned off and a cell array is returned such that
% - RESULT{I}.t is the I-th constraint type (equal to P{I}.t)
% - RESULT{I}.v is the value taken by the I-th constraint at X
  
% Written by D. Henrion, December 10, 2001 (henrion@laas.fr)
% Last modified by D. Henrion, November 26, 2002
version = '1.03';

silent = (nargout > 0); warning = 0; constr = 0;
if silent, result = {}; end;

if ~isa(P,'cell'),
 P = {P};
end;

if nargin < 2,
 error('Invalid number of input arguments');
end;

if ~isa(x,'double') | min(size(x)) > 1,
 error('Second input argument must be a vector');
end;

if nargin < 3,
 tol = 1e-6;
elseif isempty(tol),
 tol = 1e-6;
end;
if ~isa(tol,'double'),
 error('Third input argument must be a non-negative scalar');
elseif length(tol) > 1,
 error('Third input argument must be a non-negative scalar');
elseif tol < 0,
 error('Third input argument must be a non-negative scalar');
end;

N = length(P); 
n = 0; % number of variables
if ~silent,
 disp(['TestPoly ' version]);
 disp(['Tolerance = ' num2str(tol)]);
end;

for k = 1:N,
 Q = P{k};
 if ~isa(Q,'struct'),
   if ~isa(Q,'double') | issparse(Q),
     error('Invalid first input argument');
   end;
   Q.c = Q;
   if k == 1,
    Q.t = 'min';
   else
    Q.t = '>=';
   end;
 end;
 if isfield(Q,'s'),
   newn = length(Q.s);
 else
   newn = length(size(Q.c));
 end;
 if newn > n,
  n = newn;
 end;
 switch Q.t,
  case {'min', 'Min', 'MIN', 'max', 'Max', 'MAX'},
   val = evaluate(Q,x);
   if ~silent,
    disp(['Criterion = ' num2str(val)]);
   end;
  case {'>=', '=>'},
   constr = constr + 1;
   val = evaluate(Q,x);
   if (val < -tol) & ~silent,
    warning = 1;
    disp(['Warning: expression #' int2str(k) ' = ' num2str(val) ...
	  ' is negative at given tolerance']);
   end;
  case {'<=', '=<'},
   constr = constr + 1;
   val = evaluate(Q,x);
   if (val > tol) & ~silent,
    warning = 1;
    disp(['Warning: expression #' int2str(k) ' = ' num2str(val) ...
	  ' is positive at given tolerance']);
   end;
  case '==',
   constr = constr + 1;
   val = evaluate(Q,x);
   if (abs(val) > tol) & ~silent,
    warning = 1;
    disp(['Warning: expression #' int2str(k) ' = ' num2str(val) ...
	  ' is not zero at given tolerance']);
   end;
  otherwise
   error(['Invalid expression #' int2str(k)]);
 end;
 if silent,
  result{k}.t = Q.t;
  result{k}.v = val;
 end;
end;

if ~silent & ~warning,
 if constr,
  disp('Constraints are satisfied at given tolerance');
 else
  disp('No constraints');
 end;
end;

if n < length(x),
 disp('Warning: length of input vector mismatch with number of variables');
 disp('Bottom entries of input vector have been discarded');
end;
  
function P = evaluate(P,x)

if ~isfield(P,'c'),
 error('Incorrect input argument');
end;

if isfield(P,'s'),

 % sparse matrix
 % add non-zero components

 if ~issparse(P.c)
  error('Coefficient matrix must be sparse');
 end;
 nd = length(P.s);
 if nd > length(x),
  error('Invalid size of input vector');
 end;
 d = P.s;
 val = 0;
 nzndx = find(P.c);
 for j = 1:length(nzndx),
  k = [1 cumprod(d(1:nd-1))];
  ndx = nzndx(j) - 1;
  pow = zeros(1,nd);
  for i = nd:-1:1,
   pow(i) = floor(ndx/k(i));
   ndx = rem(ndx,k(i));
  end;
  coef = 1;
  for i = 1:nd,
   coef = coef*x(i)^pow(i);
  end;
  val = val + P.c(nzndx(j))*coef;
 end;
 P = val;
 
else

 % non-sparse matrix
 % Horner's scheme
 P = P.c;
 nd = ndims(P);
 if (nd == 2) & (size(P,2) == 1), nd = 1; end;
 if nd > length(x),
  error('Invalid size of input vector');
 end;
 d = size(P);
 for n = 1:nd, ind{n} = 1:d(n); end;
 for n = nd:-1:1,
  indP = {ind{1:n-1} d(n)};
  newP = P(indP{:});
  for k = d(n)-1:-1:1,
   indP = {ind{1:n-1} k};
   newP = P(indP{:}) + x(n)*newP;
  end;
  P = newP;
 end;

end;


