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.03a';

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;
%     #BEGIN CHANGES by BOSSE
% if ~isa(x,'double') | min(size(x)) > 1,
%  error('Second input argument must be a vector');
% end;
if ~isa(X,'double') | length(size(X)) > 2,
    error('Second input argument must be a vector or a matrix');
else    
    numOfSols=size(X);
    numOfSols=numOfSols(2);
end;
%     #END CHANGES by BOSSE


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(['   Solutions to check: ' num2str(numOfSols)]);
    disp(['   Tolerance = ' num2str(tol)]);    
end;


X=X';
for i = 1:numOfSols,
    if numOfSols >1,
        disp(['-------------------------------------------']);
        disp(['  [solution ',int2str(i),']']);
        
        x=X(i,:)';
        disp(x');
    else
        x=X';
    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;
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;