% Input:
% C: Matrix (C*x <= d defines the constraints)
% d: vector
% algorithm: String 'LargeScale/MediumScale/Simplex' for linprog
% 'active-set','interior-point','sqp','trust-region-reflective' for
% fmincon, 'CDD-Criss-Cross','GLPK','SeDuMi','CDD-Dual-Simplex' for
% mpt
% optim: String 'linprog,'fmincon','cvx','mpt'
% startv: vector (start vector or [] for optimization)
% TolFun: double (tolerance or negative value for default)
% MaxIter: integer (max number of iterations or negative value for default)

% Output:
% r: function handle (Support function of a polytope)

% Description:
% Generates the support function representation for the polytope Cx = d.
% The accuracy and speed of later evaluations depend heavily on the chosen
% algorithms.

function [r] = calcPolytopeSP(C,d,algorithm,optim,startv,TolFun,MaxIter)

if strcmp(optim,'linprog')
    option = optimset(@linprog);

    if strcmp(algorithm,'LargeScale')
        option.LargeScale = 'on';
        option.Simplex = 'off';
    elseif strcmp(algorithm,'MediumScale')
        option.LargeScale = 'off';
        option.Simplex = 'off';
    elseif strcmp(algorithm,'Simplex')
        %use simplex 
        option.LargeScale = 'off';
        option.Simplex = 'on';
    end;

    if TolFun > 0
        option.TolFun = TolFun;
    end;

    if MaxIter > 0
        option.MaxIter = MaxIter;
    end;

    option.Display = 'off';
    r = @(l) callPolytopeSPlinprog(l,@linprog,startv,option,C,d);
    
elseif strcmp(optim,'fmincon')
    option = optimset('Algorithm',algorithm);

    if TolFun > 0
        option.TolFun = TolFun;
    end;

    if MaxIter > 0
        option.MaxIter = MaxIter;
    end;

    option.Display = 'off';
    
    r = @(l) callPolytopeSPfmincon(l,@fmincon,startv,option,C,d);
    
elseif strcmp(optim,'mpt')
    
    opt = @mpt_solve;
    
    %S.f = -l';
    S.A = C;
    S.b = d;
    S.n = size(C,2);
    S.m = size(C,1);
    
    if strcmp(algorithm,'GLPK')
        %option = 4;
         S.solver= 'GLPK';
    elseif strcmp(algorithm,'CDD Criss-Cross')
        %option = 3;
        opt = @mpt_call_cdd;
        S.solver= 'criss-cross';
        S.problem_type ='LP';
        S.Ae = [];
        S.be = [];
        S.lb = [];
        S.ub = [];
        S.me = 0;
        S.test = false;
%         S.solver= 'CDD';
    elseif strcmp(algorithm,'linprog')
        %option = 1;
        S.solver= 'LINPROG';
    elseif strcmp(algorithm,'CDD Dual Simplex')
        %option = 5;
        opt = @mpt_call_cdd;
        S.solver= 'dual simplex';
        S.problem_type ='LP';
        S.Ae = [];
        S.be = [];
        S.lb = [];
        S.ub = [];
        S.me = 0;
        S.test = false;
    elseif strcmp(algorithm,'SeDuMi')
        %option = 6;
        S.solver= 'SEDUMI';
    else
        %option = 1; % linprog
        S.solver= 'LINPROG';
    end;
    
    %r = @(l) callPolytopeSPmpt(l,@mpt_solve,startv,option,C,d);
    
   
    r = @(l) callPolytopeSPmpt(S,l,opt);
 
elseif strcmp(optim,'cvx')
    r = @(l) callPolytopeSPcvx(l,C,d);
end;



% linprog
function [x,exitflag,supportvector] = callPolytopeSPlinprog(l,sp,startv,options,C,d)

[supportvector,ft,exitflag] = sp(-l,C,d,[],[],[],[],startv,options);
x = -ft;
if exitflag <= 0
    fprintf('Warning: the support function might be unevaluationable due to restrictions of the optimizing algorithm ! exitflag: %i\n',exitflag);
    if exitflag ~= -7
        computationerror();
    end;
end;

% fmincon
function [x,exitflag,supportvector] = callPolytopeSPfmincon(l,sp,startv,options,C,d)
f = @(x) -l' * x;
[supportvector,ft,exitflag] = sp(f,startv,C,d,[],[],[],[],[],options);
x = -ft;
if exitflag <= 0
    fprintf('Warning: the support function might be unevaluationable due to restrictions of the optimizing algorithm ! exitflag: %i\n',exitflag);
    if exitflags ~= -3
        computationerror();
    end;
end;

% MPT
function [x,exitflag,supportvector] = callPolytopeSPmpt(S,l,sp)%l,sp,~,options,C,d)

%[supportvector,ft,~,exitflag,how] = sp(-l',C,d,[],[],[],options,[],[]);
S.f = -l';
R = sp(S);

x = -R.obj;
exitflag = R.exitflag;
supportvector = R.xopt;
switch exitflag 
    case 0
        fprintf('Warning: the support function is %s ! exitflag: %i\n',R.how,exitflag);
    case -1
        fprintf('Warning: the support function is %s ! exitflag: %i\n',R.how,exitflag);
        computationerror();
end;

% CVX
function [res,exitflag,supportvector] = callPolytopeSPcvx(l,C,d)

exitflag = 1;

cvx_begin quiet
    variable x(size(C,2));
    maximize( l'*x );
    subject to
        C*x <= d;
cvx_end

supportvector = x;
res = l'*x;

function [] = computationerror()

msgbox('A constructed support function seems to be unevaluationable. Please try another algorithm for the support function representation');
error('Problem evaluating the support function with the given algorithm for support function representation.');