function [xout, lambdaout, howout, success] = mp_lp(f,A,b,VLB,VUB,x0,N,verbose,alg,opt)
% MP_LP  Linear program solver.
%
%   [xout, lambdaout, howout, success] = ...
%       mp_lp(f, A, b, VLB, VUB, x0, N, verbose, alg, opt)
%
%   A common wrapper for various LP solvers, using the calling syntax of
%   lp() from version 1 of the Optimization Toolbox, with the exception
%   that verbose == 0 means no output. The optional argument alg
%   determines the solver.
%     alg = 100  :  BPMPD_MEX
%     alg = 200  :  PDIPM, pure Matlab implementation of a primal-dual
%                   interior point method
%     alg = 250  :  SCPDIPM, pure Matlab implementation of a
%                   step-controlledprimal-dual interior point method
%     alg = 300  :  Optimization Toolbox, linprog() or lp()
%   If alg is missing or equal to zero, the first available solver is used.
%   An additional optional argument opt can be used to set algorithm
%   specific options.
%
%   From the Optimization Toolbox v.1 docs ...
%     X=LP(f,A,b) solves the linear programming problem:
%
%          min f'x    subject to:   Ax <= b 
%           x
%
%     X=LP(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=LP(f,A,b,VLB,VUB,X0) sets the initial starting point to X0.
%  
%     X=LP(f,A,b,VLB,VUB,X0,N) indicates that the first N constraints defined
%     by A and b are equality constraints.
%  
%     X=LP(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]=LP(f,A,b) returns the set of Lagrangian multipliers,
%     LAMBDA, at the solution. 
%  
%     [X,LAMBDA,HOW] = LP(f,A,b) also returns a string how that indicates 
%     error conditions at the final iteration.
%  
%     LP produces warning messages when the solution is either unbounded
%     or infeasible. 

%   MATPOWER
%   $Id: mp_lp.m,v 1.23 2009/11/04 15:58:42 ray Exp $
%   by 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.

if nargin < 9
    alg = 0;
end
if alg == 0
    if have_fcn('bpmpd')
        alg = 100;      %% BPMPD_MEX
    else
        alg = 200;      %% PDIPM (pure Matlab)
    end
end

if alg == 100
	[xout, lambdaout, howout, success] = bpmpd_qp([],f,A,b,VLB,VUB,x0,N,verbose);
	if ~success && strncmp(howout, 'BPMPD bug', 9)
		fprintf('         Retrying with PDIPM-based QP solver ...\n\n');
		mlver = ver('matlab');
		if str2double(mlver.Version(1)) < 7    %% anonymous functions not available
		  qpf = @pdipm6_qp;
		else
		  qpf = @pdipm_qp;
		end
		[xout, lambdaout, howout, success] = ...
			feval(qpf,[],f,A,b,VLB,VUB,x0,N,verbose,opt);
	end
elseif alg == 200 || alg == 250
    mlver = ver('matlab');
    if str2double(mlver.Version(1)) < 7    %% anonymous functions not available
      qpf = @pdipm6_qp;
    else
      qpf = @pdipm_qp;
    end
    if alg == 200
        opt.step_control = 0;
    else
        opt.step_control = 1;
    end
    [xout, lambdaout, howout, success] = ...
        feval(qpf,[],f,A,b,VLB,VUB,x0,N,verbose,opt);
elseif alg == 300                   %% LP solver from Opt Tbx
    success = 0;
    if have_fcn('linprog')          %% use linprog() from Opt Tbx ver 2.x+
        if N
            nb = length(b);
            Aeq = A(1:N,:);
            beq = b(1:N);
            Aieq = A(N+1:nb,:);
            bieq = b(N+1:nb);
        else
            Aeq = [];
            beq = [];
            Aieq = A;
            bieq = b;
        end
        if verbose > 1
            lpopt = optimset('Display', 'iter');
        elseif verbose == 1
            lpopt = optimset('Display', 'final');
        else
            lpopt = optimset('Display', 'off');
        end
        [xout,fval,exitflag,output,lambda] = linprog(f,Aieq,bieq,Aeq,beq,VLB,VUB,x0,lpopt);
        howout = exitflag;
        if exitflag == 1
            success = 1;
        end
        lambdaout = [   lambda.eqlin;
                        lambda.ineqlin;
                        lambda.lower;
                        lambda.upper    ];
    elseif have_fcn('lp')           %% use lp() from Opt Tbx ver 1.x/2.x
        [xout, lambdaout, howout] = lp(f,A,b,VLB,VUB,x0,N,verbose-1);
        if strcmp(howout, 'ok')
            success = 1;
        end
    else
        error('mp_lp: alg code %d requires the Optimization Toolbox.', alg);
    end
else
    error('mp_lp: %d is not a valid algorithm code', alg);
end
