function [xout, lambdaout, howout, success] = bpmpd_qp(H,f,A,b,VLB,VUB,x0,N,verbose)
% BPMPD_QP  Quadratic program solver.
%   Wrapper function for BPMPD_MEX using the calling syntax of qp()
%   from version 1 of the Optimization Toolbox, with the exception
%   that verbose == 0 means no output. 
%
%   From the Optimization Toolbox v.1 docs ...
%     X=QP(H,f,A,b) solves the quadratic programming problem:
%  
%              min 0.5*x'Hx + f'x   subject to:  Ax <= b 
%               x    
%  
%         X=QP(H,f,A,b,VLB,VUB) defines a set of lower and upper
%         bounds on the design variables, X, so that the solution  
%         is always in the range VLB <= X <= VUB.
%  
%         X=QP(H,f,A,b,VLB,VUB,X0) sets the initial starting point to X0.
%  
%         X=QP(H,f,A,b,VLB,VUB,X0,N) indicates that the first N constraints 
%         defined by A and b are equality constraints.
%  
%         X=QP(H,f,A,b,VLB,VUB,X0,N,DISPLAY) controls the level of warning
%         messages displayed.  Warning messages can be turned off with
%         DISPLAY = -1.
%  
%         [x,LAMBDA]=QP(H,f,A,b) returns the set of Lagrangian multipliers,
%         LAMBDA, at the solution.
%  
%         [X,LAMBDA,HOW] = QP(H,f,A,b) also returns a string HOW that 
%         indicates error conditions at the final iteration.
%  
%         QP produces warning messages when the solution is either unbounded
%         or infeasible. 

%   MATPOWER
%   $Id: bpmpd_qp.m,v 1.8 2009/12/04 21:42:55 ray Exp $
%   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
%   and Ray Zimmerman, PSERC Cornell
%   Copyright (c) 1996-2009 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/matpower/ for more info.

n = length(f);
m = length(b);

if nargin < 9
    verbose = 0;
    if nargin < 8
        N = 0;
        if nargin < 7
            x0 = zeros(n,1);    % Until bpmpd features warm start, this is a dummy arg
            if nargin < 6
                VUB = [];
                if nargin < 5
                    VLB = [];
                end
            end
        end
    end
end

if ~isempty(H)
    if ~issparse(H)
        H = sparse(H);
    end
end
if issparse(f)
    f = full(f);                % avoid a crash
end
if ~issparse(A)
    A = sparse(A);
end

e = -ones(m,1);
if N>0
    e(1:N,:) = zeros(N,1);
end

if ~isempty(VLB)
    llist = find(VLB > -1e15);  % interpret limits <= -1e15 as unbounded
    if isempty(llist)
        llist = [];
        lval  = [];
    else
        lval = VLB(llist);
    end
else
    llist = [];
    lval = [];
end

if ~isempty(VUB)
    ulist = find(VUB < 1e15);   % interpret limits >= 1e15 as unbounded
    if isempty(ulist)
        ulist = [];
        uval  = [];
    else
        uval = VUB(ulist);
    end
else
    ulist = [];
    uval = [];
end

if verbose > 1
    prnlev = 1;
else
    prnlev = 0;
end
if strcmp(computer, 'PCWIN')
    if prnlev
        fprintf('Windows version of BPMPD_MEX cannot print to screen.\n');
    end
    prnlev = 0;   % The DLL incarnation of bp was born mute and deaf,
end               % probably because of acute shock after realizing its fate.
                  % Can't be allowed to try to speak or its universe crumbles.

myopt = bpopt;
%myopt(14)= 1e-1;   % TPIV1  first relative pivot tolerance (desired)
%myopt(20)= 1e-10;  % TOPT1  stop if feasible and rel. dual gap less than this
%myopt(22)= 1e-6;   % TFEAS1 relative primal feasibility tolerance
%myopt(23)= 1e-6;   % TFEAS2 relative dual feasibility tolerance
%myopt(29)= 1e-10;  % TRESX  acceptable primal residual
%myopt(30)= 1e-10;  % TRESY  acceptable dual residual
%myopt(38)= 0;      % SMETHOD1 prescaling method

[xout,y,s,w,howout] = bp(H, A, b, f, e, llist, lval, ...
                        ulist, uval, myopt, prnlev);

ilt = find(w<=0);
igt = find(w>0);
mulow = zeros(n,1);
muupp = zeros(n,1);
muupp(ilt) = -w(ilt);
mulow(igt) = w(igt);

lambdaout = -y;
if ~isempty(VLB)
    lambdaout = [lambdaout; mulow];
end
if ~isempty(VUB)
    lambdaout = [lambdaout; muupp];
end

% zero out lambdas smaller than a certain tolerance
ii = find(abs(lambdaout)<1e-9);
lambdaout(ii) = zeros(size(ii));

% The next is necessary for proper operation of constr.m
if strcmp(howout, 'infeasible primal')
    lambdaout = zeros(size(lambdaout));
end
success = strcmp(howout, 'optimal solution');
if success
    % double-check feasibility
    bpmpd_bug_fatal = 0;
    err_tol = 5e-4;
    nb = length(b);
    if ~isempty(VLB)
        lb_violation = VLB - xout;
        if verbose > 1
            fprintf('max variable lower bound violatation: %g\n', max(lb_violation));
        end
    else
        lb_violation = zeros(size(xout));
    end
    if ~isempty(VUB)
        ub_violation = xout - VUB;
        if verbose > 1
            fprintf('max variable upper bound violation: %g\n', max(ub_violation));
        end
    else
        ub_violation = zeros(size(xout));
    end
    if N > 0
        eq_violation = abs( A(1:N,:) * xout - b(1:N) );
        if verbose > 1
            fprintf('max equality constraint violation: %g\n', max(eq_violation));
        end
    else
        eq_violation = zeros(N,1);
    end
    if N < nb
        ineq_violation = A(N+1:nb,:) * xout - b(N+1:nb);
        if verbose > 1
            fprintf('max inequality constraint violation: %g\n', max(ineq_violation));
        end
    else
        ineq_violation = zeros(nb-N,1);
    end

    if any( [ lb_violation;
              ub_violation;
              eq_violation;
              ineq_violation ] > err_tol)
        err_cnt = 0;
        if any( lb_violation > err_tol )
            err_cnt = err_cnt + 1;
            errs{err_cnt} = ...
                sprintf('variable lower bound violated by %g', ...
                    max(lb_violation));
        end
        if any( ub_violation > err_tol )
            err_cnt = err_cnt + 1;
            errs{err_cnt} = ... 
                sprintf('variable upper bound violated by %g', ...
                    max(ub_violation));
        end
        if any( eq_violation > err_tol )
            err_cnt = err_cnt + 1;
            errs{err_cnt} = ... 
                sprintf('equality constraint violated by %g', ...
                    max(eq_violation));
        end
        if any( ineq_violation > err_tol )
            err_cnt = err_cnt + 1;
            errs{err_cnt} = ... 
                sprintf('inequality constraint violated by %g', ...
                    max(ineq_violation));
        end
        fprintf('\nWARNING: This version of BPMPD_MEX has a bug which caused it to return\n');
        fprintf(  '         an incorrect (infeasible) solution for this particular problem.\n');
        for err_idx = 1:err_cnt
            fprintf('         %s\n', errs{err_idx});
        end
        if bpmpd_bug_fatal
            error('%s\n%s', ...
                'To avoid this bug in BPMPD_MEX you will need', ...
                'to use a different QP solver for this problem.');
        end
        success = 0;
        howout = 'BPMPD bug: returned infeasible solution';
    end
end
