% J. Ramakrishnan, 2014/07/14
% Solves qp problem:
% min. 0.5*x'*H*x + f*x
% s.t. Aineq*x <= bineq
%      Aeq*x    = beq
%      lb <= x <= ub
function [x,fval,modelstatus,output,lambda]=gamsqp(H,f,Aineq,bineq,Aeq,beq,lb,ub)

% Detect input problem structure
if nargin == 0
    error('Not enough input arguments.');
elseif nargin == 1 % for structure input argument
   if isa(H,'struct')
       if (length(fieldnames(H)) ~= sum(isfield(H,{'H',...
               'f','Aineq','bineq','Aeq','beq','lb','ub'})))
           error(['Structure contains field names different',...
               ' from H, f, Aineq, bineq, Aeq, beq, lb, and ub.']);
       end
       if all(isfield(H,{'H','f'}))
           if isfield(H,'Aineq')
               Aineq = H.Aineq;
           else
               Aineq = [];
           end
           if isfield(H,'bineq')
               bineq = H.bineq;
           else
               bineq = [];
           end
           if isfield(H,'Aeq')
               Aeq = H.Aeq;
           else
               Aeq = [];
           end
           if isfield(H,'beq')
               beq = H.beq;
           else
               beq = [];
           end
           if isfield(H,'lb')
               lb = H.lb;
           else
               lb = [];
           end
           if isfield(H,'ub')
               ub = H.ub;
           else
               ub = [];
           end  
           f = H.f;
           H = H.H;
       else
           error('Structure input must contain input arguments H and f.')
       end
   else % Single input and non-structure
        error('Single input must be a structure.');
   end
else
    if nargin < 8
        ub = [];
        if nargin < 7
            lb = [];
            if nargin < 6
                beq = [];
                if nargin < 5
                    Aeq = [];
                    if nargin < 4
                        bineq = [];
                        if nargin < 3
                            Aineq = [];
                        end
                    end
                end
            end
        end
    end
end

% Check for non-numeric inputs
if (~isnumeric(H) || ~isnumeric(f) || ~isnumeric(Aineq) || ...
    ~isnumeric(bineq) || ~isnumeric(Aeq) || ~isnumeric(beq) || ...
    ~isnumeric(lb) || ~isnumeric(ub))
    error('Must enter numeric input arguments.');
end

% Check that all data is real
if ~(isreal(H) && isreal(Aineq) && isreal(Aeq) && isreal(f) && ...
     isreal(bineq) && isreal(beq) && isreal(lb) && isreal(ub))
    error(message('Must enter real input arguments.'))
end

% check appropriate dimensions of H and f 
[m_H,n_H] = size(H);
% ismatrix checks that H is a 2D matrix
if (~ismatrix(H)) || (m_H ~= n_H)
    error('H must be a square matrix.');
end

n_f = length(f);
if (~isvector(f)) || (n_f ~= n_H)
    error('f must be a vector that has same dimensions as H.');
end

% set empty parameters appropriately 
if isempty(Aineq)
    Aineq = zeros(0,n_H);
end
if isempty(bineq)
    bineq = zeros(0,1);
end
if isempty(Aeq)
    Aeq = zeros(0,n_H); 
end
if isempty(beq)
    beq = zeros(0,1);
end
if isempty(lb)
    lb = -inf(n_H,1);
elseif ~isvector(lb)
    error('lb must be a vector.');
elseif length(lb) ~= n_H
    error('lb must be of length as number of variables.'); 
end
if isempty(ub)
    ub = inf(n_H,1);
elseif ~isvector(ub)
    error('ub must be a vector.');
elseif length(ub) ~= n_H
    error('ub must be of length as number of variables.'); 
end
                    
% Set the constraints up: defaults and check size
[m_ineq,n_ineq] = size(Aineq);
[m_eq,n_eq] = size(Aeq);

if (~ismatrix(Aineq)) || (n_ineq ~= n_H)
    error(['Aineq must be a 2D matrix with number of columns',...
        ' that has same dimensions as H.']);
end

if (~ismatrix(Aeq)) || (n_eq ~= n_H)
    error(['Aeq must be a 2D matrix with number of columns',...
        ' that has same dimensions as H.']);
end

m_b = length(bineq);
if (~isvector(bineq)) || (m_b ~= m_ineq)
    error(['bineq must be a vector that has length of the same',...
        ' number of rows as Aineq.']);
end

m_beq = length(beq);
if (~isvector(beq)) || (m_beq ~= m_eq)
    error(['beq must be a vector that has length of the same',...
        ' number of rows as Aeq.']);
end

% write data using iwgdx indexed call
iwgdx('qp_data.gdx','H',H,'f',f,'Aineq:2',Aineq,'bineq:1',bineq,...
    'Aeq:2',Aeq,'beq:1',beq,'lb:1',lb,'ub:1',ub);

% call gams to solve qp model
[modelstatus,output]=system('gams qp_model lo=3');

% read solution results
irgdx('sol.gdx','x_sol','fval','modelstatus','lambda_x',...
    'lambda_ineqlin','lambda_eqlin');

lambda.ineqlin = lambda_ineqlin;
lambda.eqlin = lambda_eqlin;
lambda.lower = max(0,lambda_x);
lambda.upper = max(0,-lambda_x);

% rename read output
x = x_sol;
clear x_sol;
